00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MATHEMATICS_H_
00014 #define __MATHEMATICS_H_
00015
00016 #include "lib/common.h"
00017 #include "lib/io.h"
00018 #include "lib/lapack.h"
00019 #include "base/SGObject.h"
00020 #include "base/Parallel.h"
00021
00022 #include <math.h>
00023 #include <stdio.h>
00024 #include <float.h>
00025 #include <pthread.h>
00026 #include <unistd.h>
00027 #include <sys/types.h>
00028 #include <sys/time.h>
00029 #include <time.h>
00030
00031 #ifdef SUNOS
00032 #include <ieeefp.h>
00033 #endif
00034
00036 #ifdef log2
00037 #define cygwin_log2 log2
00038 #undef log2
00039 #endif
00040
00041
00042
00044 #ifdef _GLIBCXX_CMATH
00045 #if _GLIBCXX_USE_C99_MATH
00046 #if !_GLIBCXX_USE_C99_FP_MACROS_DYNAMIC
00047
00049 using std::signbit;
00050
00051 using std::fpclassify;
00052
00053 using std::isfinite;
00054 using std::isinf;
00055 using std::isnan;
00056 using std::isnormal;
00057
00058 using std::isgreater;
00059 using std::isgreaterequal;
00060 using std::isless;
00061 using std::islessequal;
00062 using std::islessgreater;
00063 using std::isunordered;
00064 #endif
00065 #endif
00066 #endif
00067
00068
00069 #ifdef _WIN32
00070 #ifndef isnan
00071 #define isnan _isnan
00072 #endif
00073
00074 #ifndef isfinite
00075 #define isfinite _isfinite
00076 #endif
00077 #endif //_WIN32
00078
00079 #ifndef NAN
00080 #include <stdlib.h>
00081 #define NAN (strtod("NAN",NULL))
00082 #endif
00083
00084
00085 #define RNG_SEED_SIZE 256
00086
00087
00088 #define RADIX_STACK_SIZE 512
00089
00090
00091 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
00092 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
00093
00094 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00095
00096 template <class T> struct radix_stack_t
00097 {
00099 T *sa;
00101 size_t sn;
00103 uint16_t si;
00104 };
00105
00107
00109 template <class T1, class T2> struct thread_qsort
00110 {
00112 T1* output;
00114 T2* index;
00116 uint32_t size;
00117
00119 int32_t* qsort_threads;
00121 int32_t sort_limit;
00123 int32_t num_threads;
00124 };
00125 #endif // DOXYGEN_SHOULD_SKIP_THIS
00126
00127 namespace shogun
00128 {
00131 class CMath : public CSGObject
00132 {
00133 public:
00137
00138 CMath();
00139
00141 virtual ~CMath();
00143
00147
00149
00150 template <class T>
00151 static inline T min(T a, T b)
00152 {
00153 return (a<=b) ? a : b;
00154 }
00155
00157 template <class T>
00158 static inline T max(T a, T b)
00159 {
00160 return (a>=b) ? a : b;
00161 }
00162
00164 template <class T>
00165 static inline T clamp(T value, T lb, T ub)
00166 {
00167 if (value<=lb)
00168 return lb;
00169 else if (value>=ub)
00170 return ub;
00171 else
00172 return value;
00173 }
00174
00176 template <class T>
00177 static inline T abs(T a)
00178 {
00179
00180
00181 if (a==0)
00182 return 0;
00183 else if (a>0)
00184 return a;
00185 else
00186 return -a;
00187 }
00189
00192
00193 static inline float64_t round(float64_t d)
00194 {
00195 return ::floor(d+0.5);
00196 }
00197
00198 static inline float64_t floor(float64_t d)
00199 {
00200 return ::floor(d);
00201 }
00202
00203 static inline float64_t ceil(float64_t d)
00204 {
00205 return ::ceil(d);
00206 }
00207
00209 template <class T>
00210 static inline T sign(T a)
00211 {
00212 if (a==0)
00213 return 0;
00214 else return (a<0) ? (-1) : (+1);
00215 }
00216
00218 template <class T>
00219 static inline void swap(T &a,T &b)
00220 {
00221 T c=a;
00222 a=b;
00223 b=c;
00224 }
00225
00229 template <class T>
00230 static inline void resize(T* &data, int64_t old_size, int64_t new_size)
00231 {
00232 if (old_size==new_size)
00233 return;
00234 T* new_data = new T[new_size];
00235 for (int64_t i=0; i<old_size && i<new_size; i++)
00236 new_data[i]=data[i];
00237 delete[] data;
00238 data=new_data;
00239 }
00240
00242 template <class T>
00243 static inline T twonorm(T* x, int32_t len)
00244 {
00245 float64_t result=0;
00246 for (int32_t i=0; i<len; i++)
00247 result+=x[i]*x[i];
00248
00249 return CMath::sqrt(result);
00250 }
00251
00253 template <class T>
00254 static inline T qsq(T* x, int32_t len, float64_t q)
00255 {
00256 float64_t result=0;
00257 for (int32_t i=0; i<len; i++)
00258 result+=CMath::pow(x[i], q);
00259
00260 return result;
00261 }
00262
00264 template <class T>
00265 static inline T qnorm(T* x, int32_t len, float64_t q)
00266 {
00267 ASSERT(q!=0);
00268 return CMath::pow(qsq(x, len, q), 1/q);
00269 }
00270
00272 template <class T>
00273 static inline T sq(T x)
00274 {
00275 return x*x;
00276 }
00277
00279 static inline float32_t sqrt(float32_t x)
00280 {
00281 return ::sqrtf(x);
00282 }
00283
00285 static inline float64_t sqrt(float64_t x)
00286 {
00287 return ::sqrt(x);
00288 }
00289
00291 static inline floatmax_t sqrt(floatmax_t x)
00292 {
00293
00294
00295 #ifdef HAVE_SQRTL
00296 return ::sqrtl(x);
00297 #else
00298 return ::sqrt(x);
00299 #endif
00300 }
00301
00302
00304 static inline floatmax_t powl(floatmax_t x, floatmax_t n)
00305 {
00306
00307
00308 #ifdef HAVE_POWL
00309 return ::powl((long double) x, (long double) n);
00310 #else
00311 return ::pow((double) x, (double) n);
00312 #endif
00313 }
00314
00315 static inline int32_t pow(int32_t x, int32_t n)
00316 {
00317 ASSERT(n>=0);
00318 int32_t result=1;
00319 while (n--)
00320 result*=x;
00321
00322 return result;
00323 }
00324
00325 static inline float64_t pow(float64_t x, int32_t n)
00326 {
00327 ASSERT(n>=0);
00328 float64_t result=1;
00329 while (n--)
00330 result*=x;
00331
00332 return result;
00333 }
00334
00335 static inline float64_t pow(float64_t x, float64_t n)
00336 {
00337 return ::pow((double) x, (double) n);
00338 }
00339
00340 static inline float64_t exp(float64_t x)
00341 {
00342 return ::exp((double) x);
00343 }
00344
00345 static inline float64_t log10(float64_t v)
00346 {
00347 return ::log(v)/::log(10.0);
00348 }
00349
00350 static inline float64_t log2(float64_t v)
00351 {
00352 #ifdef HAVE_LOG2
00353 return ::log2(v);
00354 #else
00355 return ::log(v)/::log(2.0);
00356 #endif //HAVE_LOG2
00357 }
00358
00359 static inline float64_t log(float64_t v)
00360 {
00361 return ::log(v);
00362 }
00363
00364 template <class T>
00365 static void transpose_matrix(
00366 T*& matrix, int32_t& num_feat, int32_t& num_vec)
00367 {
00368 T* transposed=new T[num_vec*num_feat];
00369 for (int32_t i=0; i<num_vec; i++)
00370 {
00371 for (int32_t j=0; j<num_feat; j++)
00372 transposed[i+j*num_vec]=matrix[i*num_feat+j];
00373 }
00374
00375 delete[] matrix;
00376 matrix=transposed;
00377
00378 CMath::swap(num_feat, num_vec);
00379 }
00380
00381 #ifdef HAVE_LAPACK
00382
00383
00384 static float64_t* pinv(
00385 float64_t* matrix, int32_t rows, int32_t cols,
00386 float64_t* target=NULL);
00387
00388
00389
00390
00391 static inline void dgemm(
00392 double alpha, const double* A, int rows, int cols,
00393 CBLAS_TRANSPOSE transposeA, double *B, int cols_B,
00394 CBLAS_TRANSPOSE transposeB, double beta, double *C)
00395 {
00396 cblas_dgemm(CblasColMajor, transposeA, transposeB, rows, cols, cols_B,
00397 alpha, A, cols, B, cols_B, beta, C, cols);
00398 }
00399
00400
00401 static inline void dgemv(
00402 double alpha, const double *A, int rows, int cols,
00403 const CBLAS_TRANSPOSE transposeA, const double* X, double beta,
00404 double* Y)
00405 {
00406 cblas_dgemv(CblasColMajor, transposeA,
00407 rows, cols, alpha, A, cols,
00408 X, 1, beta, Y, 1);
00409 }
00410 #endif
00411
00412 static inline int64_t factorial(int32_t n)
00413 {
00414 int64_t res=1;
00415 for (int i=2; i<=n; i++)
00416 res*=i ;
00417 return res ;
00418 }
00419
00420 static void init_random(uint32_t initseed=0)
00421 {
00422 if (initseed==0)
00423 {
00424 struct timeval tv;
00425 gettimeofday(&tv, NULL);
00426 seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
00427 }
00428 else
00429 seed=initseed;
00430 #if !defined(CYGWIN) && !defined(__INTERIX)
00431
00432
00433 initstate(seed, CMath::rand_state, RNG_SEED_SIZE);
00434 #endif
00435 }
00436
00437 static inline int64_t random()
00438 {
00439 #if defined(CYGWIN) || defined(__INTERIX)
00440 return rand();
00441 #else
00442 return ::random();
00443 #endif
00444 }
00445
00446 static inline int32_t random(int32_t min_value, int32_t max_value)
00447 {
00448 int32_t ret = min_value + (int32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00449 ASSERT(ret>=min_value && ret<=max_value);
00450 return ret ;
00451 }
00452
00453 static inline float32_t random(float32_t min_value, float32_t max_value)
00454 {
00455 float32_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00456
00457 if (ret<min_value || ret>max_value)
00458 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00459 ASSERT(ret>=min_value && ret<=max_value);
00460 return ret;
00461 }
00462
00463 static inline float64_t random(float64_t min_value, float64_t max_value)
00464 {
00465 float64_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00466
00467 if (ret<min_value || ret>max_value)
00468 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00469 ASSERT(ret>=min_value && ret<=max_value);
00470 return ret;
00471 }
00472
00473 template <class T>
00474 static T* clone_vector(const T* vec, int32_t len)
00475 {
00476 T* result = new T[len];
00477 for (int32_t i=0; i<len; i++)
00478 result[i]=vec[i];
00479
00480 return result;
00481 }
00482 template <class T>
00483 static void fill_vector(T* vec, int32_t len, T value)
00484 {
00485 for (int32_t i=0; i<len; i++)
00486 vec[i]=value;
00487 }
00488 template <class T>
00489 static void range_fill_vector(T* vec, int32_t len, T start=0)
00490 {
00491 for (int32_t i=0; i<len; i++)
00492 vec[i]=i+start;
00493 }
00494
00495 template <class T>
00496 static void random_vector(T* vec, int32_t len, T min_value, T max_value)
00497 {
00498 for (int32_t i=0; i<len; i++)
00499 vec[i]=CMath::random(min_value, max_value);
00500 }
00501
00502 static inline int32_t* randperm(int32_t n)
00503 {
00504 int32_t* perm = new int32_t[n];
00505
00506 if (!perm)
00507 return NULL;
00508 for (int32_t i = 0; i < n; i++)
00509 perm[i] = i;
00510 for (int32_t i = 0; i < n; i++)
00511 swap(perm[random(0, n - 1)], perm[i]);
00512 return perm;
00513 }
00514
00515 static inline int64_t nchoosek(int32_t n, int32_t k)
00516 {
00517 int64_t res=1;
00518
00519 for (int32_t i=n-k+1; i<=n; i++)
00520 res*=i;
00521
00522 return res/factorial(k);
00523 }
00524
00526 template <class T>
00527 static inline void vec1_plus_scalar_times_vec2(T* vec1,
00528 T scalar, const T* vec2, int32_t n)
00529 {
00530 for (int32_t i=0; i<n; i++)
00531 vec1[i]+=scalar*vec2[i];
00532 }
00533
00535 static inline float64_t dot(const bool* v1, const bool* v2, int32_t n)
00536 {
00537 float64_t r=0;
00538 for (int32_t i=0; i<n; i++)
00539 r+=((v1[i]) ? 1 : 0) * ((v2[i]) ? 1 : 0);
00540 return r;
00541 }
00542
00544 static inline floatmax_t dot(const floatmax_t* v1, const floatmax_t* v2, int32_t n)
00545 {
00546 floatmax_t r=0;
00547 for (int32_t i=0; i<n; i++)
00548 r+=v1[i]*v2[i];
00549 return r;
00550 }
00551
00553 static inline float64_t dot(const float64_t* v1, const float64_t* v2, int32_t n)
00554 {
00555 float64_t r=0;
00556 #ifdef HAVE_LAPACK
00557 int32_t skip=1;
00558 r = cblas_ddot(n, v1, skip, v2, skip);
00559 #else
00560 for (int32_t i=0; i<n; i++)
00561 r+=v1[i]*v2[i];
00562 #endif
00563 return r;
00564 }
00565
00567 static inline float32_t dot(
00568 const float32_t* v1, const float32_t* v2, int32_t n)
00569 {
00570 float64_t r=0;
00571 #ifdef HAVE_LAPACK
00572 int32_t skip=1;
00573 r = cblas_sdot(n, v1, skip, v2, skip);
00574 #else
00575 for (int32_t i=0; i<n; i++)
00576 r+=v1[i]*v2[i];
00577 #endif
00578 return r;
00579 }
00580
00582 static inline float64_t dot(
00583 const uint64_t* v1, const uint64_t* v2, int32_t n)
00584 {
00585 float64_t r=0;
00586 for (int32_t i=0; i<n; i++)
00587 r+=((float64_t) v1[i])*v2[i];
00588
00589 return r;
00590 }
00592 static inline float64_t dot(
00593 const int64_t* v1, const int64_t* v2, int32_t n)
00594 {
00595 float64_t r=0;
00596 for (int32_t i=0; i<n; i++)
00597 r+=((float64_t) v1[i])*v2[i];
00598
00599 return r;
00600 }
00601
00603 static inline float64_t dot(
00604 const int32_t* v1, const int32_t* v2, int32_t n)
00605 {
00606 float64_t r=0;
00607 for (int32_t i=0; i<n; i++)
00608 r+=((float64_t) v1[i])*v2[i];
00609
00610 return r;
00611 }
00612
00614 static inline float64_t dot(
00615 const uint32_t* v1, const uint32_t* v2, int32_t n)
00616 {
00617 float64_t r=0;
00618 for (int32_t i=0; i<n; i++)
00619 r+=((float64_t) v1[i])*v2[i];
00620
00621 return r;
00622 }
00623
00625 static inline float64_t dot(
00626 const uint16_t* v1, const uint16_t* v2, int32_t n)
00627 {
00628 float64_t r=0;
00629 for (int32_t i=0; i<n; i++)
00630 r+=((float64_t) v1[i])*v2[i];
00631
00632 return r;
00633 }
00634
00636 static inline float64_t dot(
00637 const int16_t* v1, const int16_t* v2, int32_t n)
00638 {
00639 float64_t r=0;
00640 for (int32_t i=0; i<n; i++)
00641 r+=((float64_t) v1[i])*v2[i];
00642
00643 return r;
00644 }
00645
00647 static inline float64_t dot(
00648 const char* v1, const char* v2, int32_t n)
00649 {
00650 float64_t r=0;
00651 for (int32_t i=0; i<n; i++)
00652 r+=((float64_t) v1[i])*v2[i];
00653
00654 return r;
00655 }
00656
00658 static inline float64_t dot(
00659 const uint8_t* v1, const uint8_t* v2, int32_t n)
00660 {
00661 float64_t r=0;
00662 for (int32_t i=0; i<n; i++)
00663 r+=((float64_t) v1[i])*v2[i];
00664
00665 return r;
00666 }
00667
00669 static inline float64_t dot(
00670 const int8_t* v1, const int8_t* v2, int32_t n)
00671 {
00672 float64_t r=0;
00673 for (int32_t i=0; i<n; i++)
00674 r+=((float64_t) v1[i])*v2[i];
00675
00676 return r;
00677 }
00678
00680 static inline float64_t dot(
00681 const float64_t* v1, const char* v2, int32_t n)
00682 {
00683 float64_t r=0;
00684 for (int32_t i=0; i<n; i++)
00685 r+=((float64_t) v1[i])*v2[i];
00686
00687 return r;
00688 }
00689
00691 template <class T>
00692 static inline void add(
00693 T* target, T alpha, const T* v1, T beta, const T* v2,
00694 int32_t len)
00695 {
00696 for (int32_t i=0; i<len; i++)
00697 target[i]=alpha*v1[i]+beta*v2[i];
00698 }
00699
00701 template <class T>
00702 static inline void add_scalar(T alpha, T* vec, int32_t len)
00703 {
00704 for (int32_t i=0; i<len; i++)
00705 vec[i]+=alpha;
00706 }
00707
00709 template <class T>
00710 static inline void scale_vector(T alpha, T* vec, int32_t len)
00711 {
00712 for (int32_t i=0; i<len; i++)
00713 vec[i]*=alpha;
00714 }
00715
00717 template <class T>
00718 static inline T sum(T* vec, int32_t len)
00719 {
00720 T result=0;
00721 for (int32_t i=0; i<len; i++)
00722 result+=vec[i];
00723
00724 return result;
00725 }
00726
00728 template <class T>
00729 static inline T max(T* vec, int32_t len)
00730 {
00731 ASSERT(len>0);
00732 T maxv=vec[0];
00733
00734 for (int32_t i=1; i<len; i++)
00735 maxv=CMath::max(vec[i], maxv);
00736
00737 return maxv;
00738 }
00739
00741 template <class T>
00742 static inline T sum_abs(T* vec, int32_t len)
00743 {
00744 T result=0;
00745 for (int32_t i=0; i<len; i++)
00746 result+=CMath::abs(vec[i]);
00747
00748 return result;
00749 }
00750
00752 template <class T>
00753 static inline bool fequal(T x, T y, float64_t precision=1e-6)
00754 {
00755 return CMath::abs(x-y)<precision;
00756 }
00757
00758 static inline float64_t mean(float64_t* vec, int32_t len)
00759 {
00760 ASSERT(vec);
00761 ASSERT(len>0);
00762
00763 float64_t mean=0;
00764 for (int32_t i=0; i<len; i++)
00765 mean+=vec[i];
00766 return mean/len;
00767 }
00768
00769 static inline float64_t trace(
00770 float64_t* mat, int32_t cols, int32_t rows)
00771 {
00772 float64_t trace=0;
00773 for (int32_t i=0; i<rows; i++)
00774 trace+=mat[i*cols+i];
00775 return trace;
00776 }
00777
00781 static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
00782 static void sort(float64_t *a, int32_t*idx, int32_t N);
00783
00784
00785
00786
00787
00788
00789
00791 template <class T>
00792 inline static void radix_sort(T* array, int32_t size)
00793 {
00794 radix_sort_helper(array,size,0);
00795 }
00796
00797 template <class T>
00798 static inline uint8_t byte(T word, uint16_t p)
00799 {
00800 return (word >> (sizeof(T)-p-1) * 8) & 0xff;
00801 }
00802
00803 template <class T>
00804 static void radix_sort_helper(T* array, int32_t size, uint16_t i)
00805 {
00806 static size_t count[256], nc, cmin;
00807 T *ak;
00808 uint8_t c=0;
00809 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
00810 T *an, *aj, *pile[256];
00811 size_t *cp, cmax;
00812
00813
00814 sp = s;
00815 radix_push(array, size, i);
00816
00817
00818 while (sp>s) {
00819 radix_pop(array, size, i);
00820 an = array + size;
00821
00822
00823 if (nc == 0) {
00824 cmin = 0xff;
00825 for (ak = array; ak < an; ak++) {
00826 c = byte(*ak, i);
00827 count[c]++;
00828 if (count[c] == 1) {
00829
00830 if (c < cmin)
00831 cmin = c;
00832 nc++;
00833 }
00834 }
00835
00836
00837 if (sp + nc > s + RADIX_STACK_SIZE) {
00838 radix_sort_helper(array, size, i);
00839 continue;
00840 }
00841 }
00842
00843
00844
00845
00846
00847
00848
00849 olds = bigs = sp;
00850 cmax = 2;
00851 ak = array;
00852 pile[0xff] = an;
00853 for (cp = count + cmin; nc > 0; cp++) {
00854
00855 while (*cp == 0)
00856 cp++;
00857
00858 if (*cp > 1) {
00859
00860 if (*cp > cmax) {
00861 cmax = *cp;
00862 bigs = sp;
00863 }
00864
00865 if (i < sizeof(T)-1)
00866 radix_push(ak, *cp, (uint16_t) (i + 1));
00867 }
00868 pile[cp - count] = ak += *cp;
00869 nc--;
00870 }
00871
00872
00873 swap(*olds, *bigs);
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887 aj = array;
00888 while (aj<an)
00889 {
00890 T r;
00891
00892 for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
00893 swap(*ak, r);
00894
00895 *aj = r;
00896 aj += count[c];
00897 count[c] = 0;
00898 }
00899 }
00900 }
00901
00904 template <class T>
00905 static void insertion_sort(T* output, int32_t size)
00906 {
00907 for (int32_t i=0; i<size-1; i++)
00908 {
00909 int32_t j=i-1;
00910 T value=output[i];
00911 while (j >= 0 && output[j] > value)
00912 {
00913 output[j+1] = output[j];
00914 j--;
00915 }
00916 output[j+1]=value;
00917 }
00918 }
00919
00920
00923 template <class T>
00924 static void qsort(T* output, int32_t size)
00925 {
00926 if (size==1)
00927 return;
00928
00929 if (size==2)
00930 {
00931 if (output[0] > output [1])
00932 swap(output[0],output[1]);
00933 return;
00934 }
00935
00936 T split=output[size/2];
00937
00938 int32_t left=0;
00939 int32_t right=size-1;
00940
00941 while (left<=right)
00942 {
00943 while (output[left] < split)
00944 left++;
00945 while (output[right] > split)
00946 right--;
00947
00948 if (left<=right)
00949 {
00950 swap(output[left],output[right]);
00951 left++;
00952 right--;
00953 }
00954 }
00955
00956 if (right+1> 1)
00957 qsort(output,right+1);
00958
00959 if (size-left> 1)
00960 qsort(&output[left],size-left);
00961 }
00962
00964 template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
00965 {
00966 ASSERT(width>=0);
00967 for (int i=0; i<width; i++)
00968 {
00969 T mask = ((T) 1)<<(sizeof(T)*8-1);
00970 while (mask)
00971 {
00972 if (mask & word)
00973 SG_SPRINT("1");
00974 else
00975 SG_SPRINT("0");
00976
00977 mask>>=1;
00978 }
00979 }
00980 }
00981
00983 template <class T> static void display_vector(
00984 const T* vector, int32_t n, const char* name="vector");
00985
00987 template <class T> static void display_matrix(
00988 const T* matrix, int32_t rows, int32_t cols, const char* name="matrix");
00989
00995 template <class T1,class T2>
00996 static void qsort_index(T1* output, T2* index, uint32_t size);
00997
01003 template <class T1,class T2>
01004 static void qsort_backward_index(
01005 T1* output, T2* index, int32_t size);
01006
01014 template <class T1,class T2>
01015 inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
01016 {
01017 int32_t n=0;
01018 thread_qsort<T1,T2> t;
01019 t.output=output;
01020 t.index=index;
01021 t.size=size;
01022 t.qsort_threads=&n;
01023 t.sort_limit=limit;
01024 t.num_threads=n_threads;
01025 parallel_qsort_index<T1,T2>(&t);
01026 }
01027
01028
01029 template <class T1,class T2>
01030 static void* parallel_qsort_index(void* p);
01031
01032
01033
01034
01035 template <class T>
01036 static void min(float64_t* output, T* index, int32_t size);
01037
01038
01039
01040 template <class T>
01041 static void nmin(
01042 float64_t* output, T* index, int32_t size, int32_t n);
01043
01044
01045
01046 template <class T>
01047 static int32_t unique(T* output, int32_t size)
01048 {
01049 qsort(output, size);
01050 int32_t i,j=0 ;
01051 for (i=0; i<size; i++)
01052 if (i==0 || output[i]!=output[i-1])
01053 output[j++]=output[i];
01054 return j ;
01055 }
01056
01057
01058
01059 template <class T>
01060 static int32_t binary_search_helper(T* output, int32_t size, T elem)
01061 {
01062 int32_t start=0;
01063 int32_t end=size-1;
01064
01065 if (size<1)
01066 return 0;
01067
01068 while (start<end)
01069 {
01070 int32_t middle=(start+end)/2;
01071
01072 if (output[middle]>elem)
01073 end=middle-1;
01074 else if (output[middle]<elem)
01075 start=middle+1;
01076 else
01077 return middle;
01078 }
01079
01080 return start;
01081 }
01082
01083
01084
01085 template <class T>
01086 static inline int32_t binary_search(T* output, int32_t size, T elem)
01087 {
01088 int32_t ind = binary_search_helper(output, size, elem);
01089 if (ind >= 0 && output[ind] == elem)
01090 return ind;
01091 return -1;
01092 }
01093
01094
01095
01096
01097
01098 template <class T>
01099 static int32_t binary_search_max_lower_equal(
01100 T* output, int32_t size, T elem)
01101 {
01102 int32_t ind = binary_search_helper(output, size, elem);
01103
01104 if (output[ind]<=elem)
01105 return ind;
01106 if (ind>0 && output[ind-1] <= elem)
01107 return ind-1;
01108 return -1;
01109 }
01110
01113 static float64_t Align(
01114 char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
01115
01120 static int32_t calcroc(
01121 float64_t* fp, float64_t* tp, float64_t* output, int32_t* label,
01122 int32_t& size, int32_t& possize, int32_t& negsize,
01123 float64_t& tresh, FILE* rocfile);
01125
01128 static float64_t mutual_info(float64_t* p1, float64_t* p2, int32_t len);
01129
01132 static float64_t relative_entropy(
01133 float64_t* p, float64_t* q, int32_t len);
01134
01136 static float64_t entropy(float64_t* p, int32_t len);
01137
01139 inline static uint32_t get_seed()
01140 {
01141 return CMath::seed;
01142 }
01143
01145 inline static int is_finite(double f)
01146 {
01147 #if defined(isfinite) && !defined(SUNOS)
01148 return isfinite(f);
01149 #else
01150 return finite(f);
01151 #endif
01152 }
01153
01155 inline static int is_infinity(double f)
01156 {
01157 #ifdef SUNOS
01158 if (fpclass(f) == FP_NINF || fpclass(f) == FP_PINF)
01159 return 1;
01160 else
01161 return 0;
01162 #else
01163 return isinf(f);
01164 #endif
01165 }
01166
01168 inline static int is_nan(double f)
01169 {
01170 #ifdef SUNOS
01171 return isnand(f);
01172 #else
01173 return isnan(f);
01174 #endif
01175 }
01176
01177
01188 #ifdef USE_LOGCACHE
01189 static inline float64_t logarithmic_sum(float64_t p, float64_t q)
01190 {
01191 float64_t diff;
01192
01193 if (!CMath::finite(p))
01194 return q;
01195
01196 if (!CMath::finite(q))
01197 {
01198 SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
01199 return NAN;
01200 }
01201 diff = p - q;
01202 if (diff > 0)
01203 return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
01204 return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
01205 }
01206
01208 static void init_log_table();
01209
01211 static int32_t determine_logrange();
01212
01214 static int32_t determine_logaccuracy(int32_t range);
01215 #else
01216 static inline float64_t logarithmic_sum(
01217 float64_t p, float64_t q)
01218 {
01219 float64_t diff;
01220
01221 if (!CMath::is_finite(p))
01222 return q;
01223 if (!CMath::is_finite(q))
01224 return p;
01225 diff = p - q;
01226 if (diff > 0)
01227 return diff > LOGRANGE? p : p + log(1 + exp(-diff));
01228 return -diff > LOGRANGE? q : q + log(1 + exp(diff));
01229 }
01230 #endif
01231 #ifdef LOG_SUM_ARRAY
01232
01237 static inline float64_t logarithmic_sum_array(
01238 float64_t *p, int32_t len)
01239 {
01240 if (len<=2)
01241 {
01242 if (len==2)
01243 return logarithmic_sum(p[0],p[1]) ;
01244 if (len==1)
01245 return p[0];
01246 return -INFTY ;
01247 }
01248 else
01249 {
01250 register float64_t *pp=p ;
01251 if (len%2==1) pp++ ;
01252 for (register int32_t j=0; j < len>>1; j++)
01253 pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
01254 }
01255 return logarithmic_sum_array(p,len%2+len>>1) ;
01256 }
01257 #endif
01258
01259
01261 inline virtual const char* get_name() const { return "Mathematics"; }
01262 public:
01265
01266 static const float64_t INFTY;
01267 static const float64_t ALMOST_INFTY;
01268
01270 static const float64_t ALMOST_NEG_INFTY;
01271
01273 static int32_t LOGRANGE;
01274
01276 static uint32_t seed;
01277 static char* rand_state;
01278
01279 #ifdef USE_LOGCACHE
01280
01282 static int32_t LOGACCURACY;
01284 protected:
01286 static float64_t* logtable;
01287 #endif
01288 };
01289
01290
01291 template <class T1,class T2>
01292 void* CMath::parallel_qsort_index(void* p)
01293 {
01294 struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
01295 T1* output=ps->output;
01296 T2* index=ps->index;
01297 uint32_t size=ps->size;
01298 int32_t* qsort_threads=ps->qsort_threads;
01299 int32_t sort_limit=ps->sort_limit;
01300 int32_t num_threads=ps->num_threads;
01301
01302 if (size==2)
01303 {
01304 if (output[0] > output [1])
01305 {
01306 swap(output[0], output[1]);
01307 swap(index[0], index[1]);
01308 }
01309 return NULL;
01310 }
01311
01312 T1 split=output[size/2];
01313
01314 int32_t left=0;
01315 int32_t right=size-1;
01316
01317 while (left<=right)
01318 {
01319 while (output[left] < split)
01320 left++;
01321 while (output[right] > split)
01322 right--;
01323
01324 if (left<=right)
01325 {
01326 swap(output[left], output[right]);
01327 swap(index[left], index[right]);
01328 left++;
01329 right--;
01330 }
01331 }
01332 bool lthread_start=false;
01333 bool rthread_start=false;
01334 pthread_t lthread;
01335 pthread_t rthread;
01336 struct thread_qsort<T1,T2> t1;
01337 struct thread_qsort<T1,T2> t2;
01338
01339 if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
01340 qsort_index(output,index,right+1);
01341 else if (right+1> 1)
01342 {
01343 (*qsort_threads)++;
01344 lthread_start=true;
01345 t1.output=output;
01346 t1.index=index;
01347 t1.size=right+1;
01348 t1.qsort_threads=qsort_threads;
01349 t1.sort_limit=sort_limit;
01350 t1.num_threads=num_threads;
01351 if (pthread_create(<hread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
01352 {
01353 lthread_start=false;
01354 (*qsort_threads)--;
01355 qsort_index(output,index,right+1);
01356 }
01357 }
01358
01359
01360 if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
01361 qsort_index(&output[left],&index[left], size-left);
01362 else if (size-left> 1)
01363 {
01364 (*qsort_threads)++;
01365 rthread_start=true;
01366 t2.output=&output[left];
01367 t2.index=&index[left];
01368 t2.size=size-left;
01369 t2.qsort_threads=qsort_threads;
01370 t2.sort_limit=sort_limit;
01371 t2.num_threads=num_threads;
01372 if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
01373 {
01374 rthread_start=false;
01375 (*qsort_threads)--;
01376 qsort_index(&output[left],&index[left], size-left);
01377 }
01378 }
01379
01380 if (lthread_start)
01381 {
01382 pthread_join(lthread, NULL);
01383 (*qsort_threads)--;
01384 }
01385
01386 if (rthread_start)
01387 {
01388 pthread_join(rthread, NULL);
01389 (*qsort_threads)--;
01390 }
01391
01392 return NULL;
01393 }
01394
01395 template <class T1,class T2>
01396 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
01397 {
01398 if (size==1)
01399 return;
01400
01401 if (size==2)
01402 {
01403 if (output[0] > output [1])
01404 {
01405 swap(output[0],output[1]);
01406 swap(index[0],index[1]);
01407 }
01408 return;
01409 }
01410
01411 T1 split=output[size/2];
01412
01413 int32_t left=0;
01414 int32_t right=size-1;
01415
01416 while (left<=right)
01417 {
01418 while (output[left] < split)
01419 left++;
01420 while (output[right] > split)
01421 right--;
01422
01423 if (left<=right)
01424 {
01425 swap(output[left],output[right]);
01426 swap(index[left],index[right]);
01427 left++;
01428 right--;
01429 }
01430 }
01431
01432 if (right+1> 1)
01433 qsort_index(output,index,right+1);
01434
01435 if (size-left> 1)
01436 qsort_index(&output[left],&index[left], size-left);
01437 }
01438
01439 template <class T1,class T2>
01440 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
01441 {
01442 if (size==2)
01443 {
01444 if (output[0] < output [1])
01445 {
01446 swap(output[0],output[1]);
01447 swap(index[0],index[1]);
01448 }
01449 return;
01450 }
01451
01452
01453 T1 split=output[size/2];
01454
01455 int32_t left=0;
01456 int32_t right=size-1;
01457
01458 while (left<=right)
01459 {
01460 while (output[left] > split)
01461 left++;
01462 while (output[right] < split)
01463 right--;
01464
01465 if (left<=right)
01466 {
01467 swap(output[left],output[right]);
01468 swap(index[left],index[right]);
01469 left++;
01470 right--;
01471 }
01472 }
01473
01474 if (right+1> 1)
01475 qsort_backward_index(output,index,right+1);
01476
01477 if (size-left> 1)
01478 qsort_backward_index(&output[left],&index[left], size-left);
01479 }
01480
01481 template <class T>
01482 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
01483 {
01484 if (6*n*size<13*size*CMath::log(size))
01485 for (int32_t i=0; i<n; i++)
01486 min(&output[i], &index[i], size-i) ;
01487 else
01488 qsort_index(output, index, size) ;
01489 }
01490
01491
01492 template <class T>
01493 void CMath::min(float64_t* output, T* index, int32_t size)
01494 {
01495 if (size<=1)
01496 return;
01497 float64_t min_elem=output[0];
01498 int32_t min_index=0;
01499 for (int32_t i=1; i<size; i++)
01500 {
01501 if (output[i]<min_elem)
01502 {
01503 min_index=i;
01504 min_elem=output[i];
01505 }
01506 }
01507 swap(output[0], output[min_index]);
01508 swap(index[0], index[min_index]);
01509 }
01510 }
01511 #endif