00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #ifndef __MATHEMATICS_H_
00017 #define __MATHEMATICS_H_
00018
00019 #include <shogun/lib/common.h>
00020 #include <shogun/io/SGIO.h>
00021 #include <shogun/mathematics/lapack.h>
00022 #include <shogun/base/SGObject.h>
00023 #include <shogun/base/Parallel.h>
00024
00025 #include <math.h>
00026 #include <stdio.h>
00027 #include <float.h>
00028 #include <pthread.h>
00029 #include <unistd.h>
00030 #include <sys/types.h>
00031 #include <sys/time.h>
00032 #include <time.h>
00033
00034 #ifdef SUNOS
00035 #include <ieeefp.h>
00036 #endif
00037
00039 #ifdef log2
00040 #define cygwin_log2 log2
00041 #undef log2
00042 #endif
00043
00044
00045
00047 #ifdef _GLIBCXX_CMATH
00048 #if _GLIBCXX_USE_C99_MATH
00049 #if !_GLIBCXX_USE_C99_FP_MACROS_DYNAMIC
00050
00052 using std::signbit;
00053
00054 using std::fpclassify;
00055
00056 using std::isfinite;
00057 using std::isinf;
00058 using std::isnan;
00059 using std::isnormal;
00060
00061 using std::isgreater;
00062 using std::isgreaterequal;
00063 using std::isless;
00064 using std::islessequal;
00065 using std::islessgreater;
00066 using std::isunordered;
00067 #endif
00068 #endif
00069 #endif
00070
00071
00072 #ifdef _WIN32
00073 #ifndef isnan
00074 #define isnan _isnan
00075 #endif
00076
00077 #ifndef isfinite
00078 #define isfinite _isfinite
00079 #endif
00080 #endif //_WIN32
00081
00082 #ifndef NAN
00083 #include <stdlib.h>
00084 #define NAN (strtod("NAN",NULL))
00085 #endif
00086
00087
00088 #define RNG_SEED_SIZE 256
00089
00090
00091 #define RADIX_STACK_SIZE 512
00092
00093
00094 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
00095 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
00096
00097 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00098
00099 template <class T> struct radix_stack_t
00100 {
00102 T *sa;
00104 size_t sn;
00106 uint16_t si;
00107 };
00108
00110
00112 template <class T1, class T2> struct thread_qsort
00113 {
00115 T1* output;
00117 T2* index;
00119 uint32_t size;
00120
00122 int32_t* qsort_threads;
00124 int32_t sort_limit;
00126 int32_t num_threads;
00127 };
00128 #endif // DOXYGEN_SHOULD_SKIP_THIS
00129
00130 namespace shogun
00131 {
00134 class CMath : public CSGObject
00135 {
00136 public:
00140
00141 CMath();
00142
00144 virtual ~CMath();
00146
00150
00152
00153 template <class T>
00154 static inline T min(T a, T b)
00155 {
00156 return (a<=b) ? a : b;
00157 }
00158
00160 template <class T>
00161 static inline T max(T a, T b)
00162 {
00163 return (a>=b) ? a : b;
00164 }
00165
00167 template <class T>
00168 static inline T clamp(T value, T lb, T ub)
00169 {
00170 if (value<=lb)
00171 return lb;
00172 else if (value>=ub)
00173 return ub;
00174 else
00175 return value;
00176 }
00177
00179 template <class T>
00180 static inline T abs(T a)
00181 {
00182
00183
00184 if (a==0)
00185 return 0;
00186 else if (a>0)
00187 return a;
00188 else
00189 return -a;
00190 }
00192
00195
00196 static inline float64_t round(float64_t d)
00197 {
00198 return ::floor(d+0.5);
00199 }
00200
00201 static inline float64_t floor(float64_t d)
00202 {
00203 return ::floor(d);
00204 }
00205
00206 static inline float64_t ceil(float64_t d)
00207 {
00208 return ::ceil(d);
00209 }
00210
00212 template <class T>
00213 static inline T sign(T a)
00214 {
00215 if (a==0)
00216 return 0;
00217 else return (a<0) ? (-1) : (+1);
00218 }
00219
00221 template <class T>
00222 static inline void swap(T &a,T &b)
00223 {
00224 T c=a;
00225 a=b;
00226 b=c;
00227 }
00228
00232 template <class T>
00233 static inline void resize(T* &data, int64_t old_size, int64_t new_size)
00234 {
00235 if (old_size==new_size)
00236 return;
00237 T* new_data = SG_MALLOC(T, new_size);
00238 for (int64_t i=0; i<old_size && i<new_size; i++)
00239 new_data[i]=data[i];
00240 SG_FREE(data);
00241 data=new_data;
00242 }
00243
00245 template <class T>
00246 static inline T twonorm(T* x, int32_t len)
00247 {
00248 float64_t result=0;
00249 for (int32_t i=0; i<len; i++)
00250 result+=x[i]*x[i];
00251
00252 return CMath::sqrt(result);
00253 }
00254
00256 template <class T>
00257 static inline T qsq(T* x, int32_t len, float64_t q)
00258 {
00259 float64_t result=0;
00260 for (int32_t i=0; i<len; i++)
00261 result+=CMath::pow(x[i], q);
00262
00263 return result;
00264 }
00265
00267 template <class T>
00268 static inline T qnorm(T* x, int32_t len, float64_t q)
00269 {
00270 ASSERT(q!=0);
00271 return CMath::pow(qsq(x, len, q), 1/q);
00272 }
00273
00275 template <class T>
00276 static inline T sq(T x)
00277 {
00278 return x*x;
00279 }
00280
00282 static inline float32_t sqrt(float32_t x)
00283 {
00284 return ::sqrtf(x);
00285 }
00286
00288 static inline float64_t sqrt(float64_t x)
00289 {
00290 return ::sqrt(x);
00291 }
00292
00294 static inline floatmax_t sqrt(floatmax_t x)
00295 {
00296
00297
00298 #ifdef HAVE_SQRTL
00299 return ::sqrtl(x);
00300 #else
00301 return ::sqrt(x);
00302 #endif
00303 }
00304
00306 static inline float32_t invsqrt(float32_t x)
00307 {
00308 union float_to_int
00309 {
00310 float32_t f;
00311 int32_t i;
00312 };
00313
00314 float_to_int tmp;
00315 tmp.f=x;
00316
00317 float32_t xhalf = 0.5f * x;
00318
00319
00320 tmp.i = 0x5f3759d5 - (tmp.i >> 1);
00321 x = tmp.f;
00322 x = x*(1.5f - xhalf*x*x);
00323 return x;
00324 }
00325
00327 static inline floatmax_t powl(floatmax_t x, floatmax_t n)
00328 {
00329
00330
00331 #ifdef HAVE_POWL
00332 return ::powl((long double) x, (long double) n);
00333 #else
00334 return ::pow((double) x, (double) n);
00335 #endif
00336 }
00337
00338 static inline int32_t pow(int32_t x, int32_t n)
00339 {
00340 ASSERT(n>=0);
00341 int32_t result=1;
00342 while (n--)
00343 result*=x;
00344
00345 return result;
00346 }
00347
00348 static inline float64_t pow(float64_t x, int32_t n)
00349 {
00350 if (n>=0)
00351 {
00352 float64_t result=1;
00353 while (n--)
00354 result*=x;
00355
00356 return result;
00357 }
00358 else
00359 return ::pow((double)x, (double)n);
00360 }
00361
00362 static inline float64_t pow(float64_t x, float64_t n)
00363 {
00364 return ::pow((double) x, (double) n);
00365 }
00366
00367 static inline float64_t exp(float64_t x)
00368 {
00369 return ::exp((double) x);
00370 }
00371
00373 static inline float64_t lgamma(float64_t x)
00374 {
00375 return ::lgamma((double) x);
00376 }
00377
00379 static inline float64_t tgamma(float64_t x)
00380 {
00381 return ::tgamma((double) x);
00382 }
00383
00385 static inline float64_t atan(float64_t x)
00386 {
00387 return ::atan((double) x);
00388 }
00389
00390 static inline floatmax_t lgammal(floatmax_t x)
00391 {
00392 return ::lgammal((long double) x);
00393 }
00394
00395 static inline float64_t log10(float64_t v)
00396 {
00397 return ::log(v)/::log(10.0);
00398 }
00399
00400 static inline float64_t log2(float64_t v)
00401 {
00402 #ifdef HAVE_LOG2
00403 return ::log2(v);
00404 #else
00405 return ::log(v)/::log(2.0);
00406 #endif //HAVE_LOG2
00407 }
00408
00409 static inline float64_t log(float64_t v)
00410 {
00411 return ::log(v);
00412 }
00413
00414 static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
00415 {
00416 ASSERT(len>0 && xy);
00417
00418 float64_t area = 0.0;
00419
00420 if (!reversed)
00421 {
00422 for (int i=1; i<len; i++)
00423 area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
00424 }
00425 else
00426 {
00427 for (int i=1; i<len; i++)
00428 area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
00429 }
00430
00431 return area;
00432 }
00433
00434 template <class T>
00435 static void transpose_matrix(
00436 T*& matrix, int32_t& num_feat, int32_t& num_vec)
00437 {
00438 T* transposed=SG_MALLOC(T, num_vec*num_feat);
00439 for (int32_t i=0; i<num_vec; i++)
00440 {
00441 for (int32_t j=0; j<num_feat; j++)
00442 transposed[i+j*num_vec]=matrix[i*num_feat+j];
00443 }
00444
00445 SG_FREE(matrix);
00446 matrix=transposed;
00447
00448 CMath::swap(num_feat, num_vec);
00449 }
00450
00451 #ifdef HAVE_LAPACK
00452
00453
00454 static float64_t* pinv(
00455 float64_t* matrix, int32_t rows, int32_t cols,
00456 float64_t* target=NULL);
00457
00458
00459
00460
00461 static inline void dgemm(
00462 double alpha, const double* A, int rows, int cols,
00463 CBLAS_TRANSPOSE transposeA, double *B, int cols_B,
00464 CBLAS_TRANSPOSE transposeB, double beta, double *C)
00465 {
00466 cblas_dgemm(CblasColMajor, transposeA, transposeB, rows, cols, cols_B,
00467 alpha, A, cols, B, cols_B, beta, C, cols);
00468 }
00469
00470
00471 static inline void dgemv(
00472 double alpha, const double *A, int rows, int cols,
00473 const CBLAS_TRANSPOSE transposeA, const double* X, double beta,
00474 double* Y)
00475 {
00476 cblas_dgemv(CblasColMajor, transposeA,
00477 rows, cols, alpha, A, cols,
00478 X, 1, beta, Y, 1);
00479 }
00480 #endif
00481
00482 static inline int64_t factorial(int32_t n)
00483 {
00484 int64_t res=1;
00485 for (int i=2; i<=n; i++)
00486 res*=i ;
00487 return res ;
00488 }
00489
00490 static void init_random(uint32_t initseed=0)
00491 {
00492 if (initseed==0)
00493 {
00494 struct timeval tv;
00495 gettimeofday(&tv, NULL);
00496 seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
00497 }
00498 else
00499 seed=initseed;
00500 #if !defined(CYGWIN) && !defined(__INTERIX)
00501
00502
00503 initstate(seed, CMath::rand_state, RNG_SEED_SIZE);
00504 #endif
00505 }
00506
00507 static inline int64_t random()
00508 {
00509 #if defined(CYGWIN) || defined(__INTERIX)
00510 return rand();
00511 #else
00512 return ::random();
00513 #endif
00514 }
00515
00516 static inline int32_t random(int32_t min_value, int32_t max_value)
00517 {
00518 int32_t ret = min_value + (int32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00519 ASSERT(ret>=min_value && ret<=max_value);
00520 return ret ;
00521 }
00522
00523 static inline float32_t random(float32_t min_value, float32_t max_value)
00524 {
00525 float32_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00526
00527 if (ret<min_value || ret>max_value)
00528 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00529 ASSERT(ret>=min_value && ret<=max_value);
00530 return ret;
00531 }
00532
00533 static inline float64_t random(float64_t min_value, float64_t max_value)
00534 {
00535 float64_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00536
00537 if (ret<min_value || ret>max_value)
00538 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00539 ASSERT(ret>=min_value && ret<=max_value);
00540 return ret;
00541 }
00542
00546 static inline float32_t normal_random(float32_t mean, float32_t std_dev)
00547 {
00548
00549 float32_t ret;
00550 float32_t rand_u;
00551 float32_t rand_v;
00552 float32_t rand_s;
00553 do
00554 {
00555 rand_u = random(-1.0, 1.0);
00556 rand_v = random(-1.0, 1.0);
00557 rand_s = rand_u*rand_u + rand_v*rand_v;
00558 } while ((rand_s == 0) || (rand_s >= 1));
00559
00560
00561 ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
00562 ret = std_dev*ret + mean;
00563 return ret;
00564 }
00565
00569 static inline float64_t normal_random(float64_t mean, float64_t std_dev)
00570 {
00571 float64_t ret;
00572 float64_t rand_u;
00573 float64_t rand_v;
00574 float64_t rand_s;
00575 do
00576 {
00577 rand_u = random(-1.0, 1.0);
00578 rand_v = random(-1.0, 1.0);
00579 rand_s = rand_u*rand_u + rand_v*rand_v;
00580 } while ((rand_s == 0) || (rand_s >= 1));
00581
00582 ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
00583 ret = std_dev*ret + mean;
00584 return ret;
00585 }
00586
00589 static inline float32_t randn_float()
00590 {
00591 return normal_random(0.0, 1.0);
00592 }
00593
00596 static inline float64_t randn_double()
00597 {
00598 return normal_random(0.0, 1.0);
00599 }
00600
00601 template <class T>
00602 static T* clone_vector(const T* vec, int32_t len)
00603 {
00604 T* result = SG_MALLOC(T, len);
00605 for (int32_t i=0; i<len; i++)
00606 result[i]=vec[i];
00607
00608 return result;
00609 }
00610 template <class T>
00611 static void fill_vector(T* vec, int32_t len, T value)
00612 {
00613 for (int32_t i=0; i<len; i++)
00614 vec[i]=value;
00615 }
00616 template <class T>
00617 static void range_fill_vector(T* vec, int32_t len, T start=0)
00618 {
00619 for (int32_t i=0; i<len; i++)
00620 vec[i]=i+start;
00621 }
00622
00623 template <class T>
00624 static void random_vector(T* vec, int32_t len, T min_value, T max_value)
00625 {
00626 for (int32_t i=0; i<len; i++)
00627 vec[i]=CMath::random(min_value, max_value);
00628 }
00629
00630 static inline int32_t* randperm(int32_t n)
00631 {
00632 int32_t* perm = SG_MALLOC(int32_t, n);
00633
00634 if (!perm)
00635 return NULL;
00636 for (int32_t i = 0; i < n; i++)
00637 perm[i] = i;
00638 for (int32_t i = 0; i < n; i++)
00639 swap(perm[random(0, n - 1)], perm[i]);
00640 return perm;
00641 }
00642
00643 static inline int64_t nchoosek(int32_t n, int32_t k)
00644 {
00645 int64_t res=1;
00646
00647 for (int32_t i=n-k+1; i<=n; i++)
00648 res*=i;
00649
00650 return res/factorial(k);
00651 }
00652
00654 template <class T>
00655 static inline void vec1_plus_scalar_times_vec2(T* vec1,
00656 T scalar, const T* vec2, int32_t n)
00657 {
00658 for (int32_t i=0; i<n; i++)
00659 vec1[i]+=scalar*vec2[i];
00660 }
00661
00663 static inline float64_t dot(const bool* v1, const bool* v2, int32_t n)
00664 {
00665 float64_t r=0;
00666 for (int32_t i=0; i<n; i++)
00667 r+=((v1[i]) ? 1 : 0) * ((v2[i]) ? 1 : 0);
00668 return r;
00669 }
00670
00672 static inline floatmax_t dot(const floatmax_t* v1, const floatmax_t* v2, int32_t n)
00673 {
00674 floatmax_t r=0;
00675 for (int32_t i=0; i<n; i++)
00676 r+=v1[i]*v2[i];
00677 return r;
00678 }
00679
00680
00682 static inline float64_t dot(const float64_t* v1, const float64_t* v2, int32_t n)
00683 {
00684 float64_t r=0;
00685 #ifdef HAVE_LAPACK
00686 int32_t skip=1;
00687 r = cblas_ddot(n, v1, skip, v2, skip);
00688 #else
00689 for (int32_t i=0; i<n; i++)
00690 r+=v1[i]*v2[i];
00691 #endif
00692 return r;
00693 }
00694
00696 static inline float32_t dot(
00697 const float32_t* v1, const float32_t* v2, int32_t n)
00698 {
00699 float64_t r=0;
00700 #ifdef HAVE_LAPACK
00701 int32_t skip=1;
00702 r = cblas_sdot(n, v1, skip, v2, skip);
00703 #else
00704 for (int32_t i=0; i<n; i++)
00705 r+=v1[i]*v2[i];
00706 #endif
00707 return r;
00708 }
00709
00711 static inline float64_t dot(
00712 const uint64_t* v1, const uint64_t* v2, int32_t n)
00713 {
00714 float64_t r=0;
00715 for (int32_t i=0; i<n; i++)
00716 r+=((float64_t) v1[i])*v2[i];
00717
00718 return r;
00719 }
00721 static inline float64_t dot(
00722 const int64_t* v1, const int64_t* v2, int32_t n)
00723 {
00724 float64_t r=0;
00725 for (int32_t i=0; i<n; i++)
00726 r+=((float64_t) v1[i])*v2[i];
00727
00728 return r;
00729 }
00730
00732 static inline float64_t dot(
00733 const int32_t* v1, const int32_t* v2, int32_t n)
00734 {
00735 float64_t r=0;
00736 for (int32_t i=0; i<n; i++)
00737 r+=((float64_t) v1[i])*v2[i];
00738
00739 return r;
00740 }
00741
00743 static inline float64_t dot(
00744 const uint32_t* v1, const uint32_t* v2, int32_t n)
00745 {
00746 float64_t r=0;
00747 for (int32_t i=0; i<n; i++)
00748 r+=((float64_t) v1[i])*v2[i];
00749
00750 return r;
00751 }
00752
00754 static inline float64_t dot(
00755 const uint16_t* v1, const uint16_t* v2, int32_t n)
00756 {
00757 float64_t r=0;
00758 for (int32_t i=0; i<n; i++)
00759 r+=((float64_t) v1[i])*v2[i];
00760
00761 return r;
00762 }
00763
00765 static inline float64_t dot(
00766 const int16_t* v1, const int16_t* v2, int32_t n)
00767 {
00768 float64_t r=0;
00769 for (int32_t i=0; i<n; i++)
00770 r+=((float64_t) v1[i])*v2[i];
00771
00772 return r;
00773 }
00774
00776 static inline float64_t dot(
00777 const char* v1, const char* v2, int32_t n)
00778 {
00779 float64_t r=0;
00780 for (int32_t i=0; i<n; i++)
00781 r+=((float64_t) v1[i])*v2[i];
00782
00783 return r;
00784 }
00785
00787 static inline float64_t dot(
00788 const uint8_t* v1, const uint8_t* v2, int32_t n)
00789 {
00790 float64_t r=0;
00791 for (int32_t i=0; i<n; i++)
00792 r+=((float64_t) v1[i])*v2[i];
00793
00794 return r;
00795 }
00796
00798 static inline float64_t dot(
00799 const int8_t* v1, const int8_t* v2, int32_t n)
00800 {
00801 float64_t r=0;
00802 for (int32_t i=0; i<n; i++)
00803 r+=((float64_t) v1[i])*v2[i];
00804
00805 return r;
00806 }
00807
00809 static inline float64_t dot(
00810 const float64_t* v1, const char* v2, int32_t n)
00811 {
00812 float64_t r=0;
00813 for (int32_t i=0; i<n; i++)
00814 r+=((float64_t) v1[i])*v2[i];
00815
00816 return r;
00817 }
00818
00820 template <class T>
00821 static inline void vector_multiply(
00822 T* target, const T* v1, const T* v2,int32_t len)
00823 {
00824 for (int32_t i=0; i<len; i++)
00825 target[i]=v1[i]*v2[i];
00826 }
00827
00828
00830 template <class T>
00831 static inline void add(
00832 T* target, T alpha, const T* v1, T beta, const T* v2,
00833 int32_t len)
00834 {
00835 for (int32_t i=0; i<len; i++)
00836 target[i]=alpha*v1[i]+beta*v2[i];
00837 }
00838
00840 template <class T>
00841 static inline void add_scalar(T alpha, T* vec, int32_t len)
00842 {
00843 for (int32_t i=0; i<len; i++)
00844 vec[i]+=alpha;
00845 }
00846
00848 template <class T>
00849 static inline void scale_vector(T alpha, T* vec, int32_t len)
00850 {
00851 for (int32_t i=0; i<len; i++)
00852 vec[i]*=alpha;
00853 }
00854
00856 template <class T>
00857 static inline T sum(T* vec, int32_t len)
00858 {
00859 T result=0;
00860 for (int32_t i=0; i<len; i++)
00861 result+=vec[i];
00862
00863 return result;
00864 }
00865
00867 template <class T>
00868 static inline T max(T* vec, int32_t len)
00869 {
00870 ASSERT(len>0);
00871 T maxv=vec[0];
00872
00873 for (int32_t i=1; i<len; i++)
00874 maxv=CMath::max(vec[i], maxv);
00875
00876 return maxv;
00877 }
00878
00880 template <class T>
00881 static inline T sum_abs(T* vec, int32_t len)
00882 {
00883 T result=0;
00884 for (int32_t i=0; i<len; i++)
00885 result+=CMath::abs(vec[i]);
00886
00887 return result;
00888 }
00889
00891 template <class T>
00892 static inline bool fequal(T x, T y, float64_t precision=1e-6)
00893 {
00894 return CMath::abs(x-y)<precision;
00895 }
00896
00898 static inline float64_t mean(float64_t* vec, int32_t len)
00899 {
00900 SG_SDEPRECATED;
00901 ASSERT(vec);
00902 ASSERT(len>0);
00903
00904 float64_t mean=0;
00905 for (int32_t i=0; i<len; i++)
00906 mean+=vec[i];
00907 return mean/len;
00908 }
00909
00910 static inline float64_t trace(
00911 float64_t* mat, int32_t cols, int32_t rows)
00912 {
00913 float64_t trace=0;
00914 for (int32_t i=0; i<rows; i++)
00915 trace+=mat[i*cols+i];
00916 return trace;
00917 }
00918
00922 static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
00923 static void sort(float64_t *a, int32_t*idx, int32_t N);
00924
00925
00926
00927
00928
00929
00930
00932 template <class T>
00933 inline static void radix_sort(T* array, int32_t size)
00934 {
00935 radix_sort_helper(array,size,0);
00936 }
00937
00938 template <class T>
00939 static inline uint8_t byte(T word, uint16_t p)
00940 {
00941 return (word >> (sizeof(T)-p-1) * 8) & 0xff;
00942 }
00943
00944 template <class T>
00945 static void radix_sort_helper(T* array, int32_t size, uint16_t i)
00946 {
00947 static size_t count[256], nc, cmin;
00948 T *ak;
00949 uint8_t c=0;
00950 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
00951 T *an, *aj, *pile[256];
00952 size_t *cp, cmax;
00953
00954
00955 sp = s;
00956 radix_push(array, size, i);
00957
00958
00959 while (sp>s) {
00960 radix_pop(array, size, i);
00961 an = array + size;
00962
00963
00964 if (nc == 0) {
00965 cmin = 0xff;
00966 for (ak = array; ak < an; ak++) {
00967 c = byte(*ak, i);
00968 count[c]++;
00969 if (count[c] == 1) {
00970
00971 if (c < cmin)
00972 cmin = c;
00973 nc++;
00974 }
00975 }
00976
00977
00978 if (sp + nc > s + RADIX_STACK_SIZE) {
00979 radix_sort_helper(array, size, i);
00980 continue;
00981 }
00982 }
00983
00984
00985
00986
00987
00988
00989
00990 olds = bigs = sp;
00991 cmax = 2;
00992 ak = array;
00993 pile[0xff] = an;
00994 for (cp = count + cmin; nc > 0; cp++) {
00995
00996 while (*cp == 0)
00997 cp++;
00998
00999 if (*cp > 1) {
01000
01001 if (*cp > cmax) {
01002 cmax = *cp;
01003 bigs = sp;
01004 }
01005
01006 if (i < sizeof(T)-1)
01007 radix_push(ak, *cp, (uint16_t) (i + 1));
01008 }
01009 pile[cp - count] = ak += *cp;
01010 nc--;
01011 }
01012
01013
01014 swap(*olds, *bigs);
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028 aj = array;
01029 while (aj<an)
01030 {
01031 T r;
01032
01033 for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
01034 swap(*ak, r);
01035
01036 *aj = r;
01037 aj += count[c];
01038 count[c] = 0;
01039 }
01040 }
01041 }
01042
01045 template <class T>
01046 static void insertion_sort(T* output, int32_t size)
01047 {
01048 for (int32_t i=0; i<size-1; i++)
01049 {
01050 int32_t j=i-1;
01051 T value=output[i];
01052 while (j >= 0 && output[j] > value)
01053 {
01054 output[j+1] = output[j];
01055 j--;
01056 }
01057 output[j+1]=value;
01058 }
01059 }
01060
01061
01064 template <class T>
01065 static void qsort(T* output, int32_t size)
01066 {
01067 if (size==1)
01068 return;
01069
01070 if (size==2)
01071 {
01072 if (output[0] > output [1])
01073 swap(output[0],output[1]);
01074 return;
01075 }
01076
01077 T split=output[size/2];
01078
01079 int32_t left=0;
01080 int32_t right=size-1;
01081
01082 while (left<=right)
01083 {
01084 while (output[left] < split)
01085 left++;
01086 while (output[right] > split)
01087 right--;
01088
01089 if (left<=right)
01090 {
01091 swap(output[left],output[right]);
01092 left++;
01093 right--;
01094 }
01095 }
01096
01097 if (right+1> 1)
01098 qsort(output,right+1);
01099
01100 if (size-left> 1)
01101 qsort(&output[left],size-left);
01102 }
01103
01112 template <class T>
01113 static void qsort(SGVector<T*> array)
01114 {
01115 if (array.vlen==1)
01116 return;
01117
01118 if (array.vlen==2)
01119 {
01120 if (*array.vector[0]>*array.vector[1])
01121 swap(array.vector[0],array.vector[1]);
01122 return;
01123 }
01124 T* split=array.vector[array.vlen/2];
01125
01126 int32_t left=0;
01127 int32_t right=array.vlen-1;
01128
01129 while (left<=right)
01130 {
01131 while (*array.vector[left]<*split)
01132 ++left;
01133 while (*array.vector[right]>*split)
01134 --right;
01135
01136 if (left<=right)
01137 {
01138 swap(array.vector[left],array.vector[right]);
01139 ++left;
01140 --right;
01141 }
01142 }
01143
01144 if (right+1>1)
01145 qsort(SGVector<T*>(array.vector,right+1));
01146
01147 if (array.vlen-left>1)
01148 qsort(SGVector<T*>(&array.vector[left],array.vlen-left));
01149 }
01150
01152 template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
01153 {
01154 ASSERT(width>=0);
01155 for (int i=0; i<width; i++)
01156 {
01157 T mask = ((T) 1)<<(sizeof(T)*8-1);
01158 while (mask)
01159 {
01160 if (mask & word)
01161 SG_SPRINT("1");
01162 else
01163 SG_SPRINT("0");
01164
01165 mask>>=1;
01166 }
01167 }
01168 }
01169
01171 template <class T> static void display_vector(
01172 const T* vector, int32_t n, const char* name="vector");
01173
01175 template <class T> static void display_matrix(
01176 const T* matrix, int32_t rows, int32_t cols, const char* name="matrix");
01177
01183 template <class T1,class T2>
01184 static void qsort_index(T1* output, T2* index, uint32_t size);
01185
01191 template <class T1,class T2>
01192 static void qsort_backward_index(
01193 T1* output, T2* index, int32_t size);
01194
01202 template <class T1,class T2>
01203 inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
01204 {
01205 int32_t n=0;
01206 thread_qsort<T1,T2> t;
01207 t.output=output;
01208 t.index=index;
01209 t.size=size;
01210 t.qsort_threads=&n;
01211 t.sort_limit=limit;
01212 t.num_threads=n_threads;
01213 parallel_qsort_index<T1,T2>(&t);
01214 }
01215
01216
01217 template <class T1,class T2>
01218 static void* parallel_qsort_index(void* p);
01219
01220
01221
01222
01223 template <class T>
01224 static void min(float64_t* output, T* index, int32_t size);
01225
01226
01227
01228 template <class T>
01229 static void nmin(
01230 float64_t* output, T* index, int32_t size, int32_t n);
01231
01232
01233
01234 template <class T>
01235 static int32_t unique(T* output, int32_t size)
01236 {
01237 qsort(output, size);
01238 int32_t i,j=0 ;
01239 for (i=0; i<size; i++)
01240 if (i==0 || output[i]!=output[i-1])
01241 output[j++]=output[i];
01242 return j ;
01243 }
01244
01252 static double* compute_eigenvectors(double* matrix, int n, int m)
01253 {
01254 #ifdef HAVE_LAPACK
01255 ASSERT(n == m);
01256
01257 char V='V';
01258 char U='U';
01259 int info;
01260 int ord=n;
01261 int lda=n;
01262 double* eigenvalues=SG_CALLOC(float64_t, n+1);
01263
01264
01265 wrap_dsyev(V, U, ord, matrix, lda,
01266 eigenvalues, &info);
01267
01268 if (info!=0)
01269 SG_SERROR("DSYEV failed with code %d\n", info);
01270
01271 return eigenvalues;
01272 #else
01273 SG_SERROR("Function not available - Lapack/Atlas not enabled at compile time!\n");
01274 return NULL;
01275 #endif
01276 }
01277
01278
01279 template <class T>
01280 static T* get_row_sum(T* matrix, int32_t m, int32_t n)
01281 {
01282 T* rowsums=SG_CALLOC(T, n);
01283
01284 for (int32_t i=0; i<n; i++)
01285 {
01286 for (int32_t j=0; j<m; j++)
01287 rowsums[i]+=matrix[j+int64_t(i)*m];
01288 }
01289 return rowsums;
01290 }
01291
01292
01293 template <class T>
01294 static T* get_column_sum(T* matrix, int32_t m, int32_t n)
01295 {
01296 T* colsums=SG_CALLOC(T, m);
01297
01298 for (int32_t i=0; i<n; i++)
01299 {
01300 for (int32_t j=0; j<m; j++)
01301 colsums[j]+=matrix[j+int64_t(i)*m];
01302 }
01303 return colsums;
01304 }
01305
01306
01307 template <class T>
01308 static void center_matrix(T* matrix, int32_t m, int32_t n)
01309 {
01310 float64_t num_data=n;
01311
01312 T* colsums=get_column_sum(matrix, m,n);
01313 T* rowsums=get_row_sum(matrix, m,n);
01314
01315 for (int32_t i=0; i<m; i++)
01316 colsums[i]/=num_data;
01317 for (int32_t j=0; j<n; j++)
01318 rowsums[j]/=num_data;
01319
01320 T s=sum(rowsums, n)/num_data;
01321
01322 for (int32_t i=0; i<n; i++)
01323 {
01324 for (int32_t j=0; j<m; j++)
01325 matrix[int64_t(i)*m+j]+=s-colsums[j]-rowsums[i];
01326 }
01327
01328 SG_FREE(rowsums);
01329 SG_FREE(colsums);
01330 }
01331
01332
01333
01334 template <class T>
01335 static int32_t binary_search_helper(T* output, int32_t size, T elem)
01336 {
01337 int32_t start=0;
01338 int32_t end=size-1;
01339
01340 if (size<1)
01341 return 0;
01342
01343 while (start<end)
01344 {
01345 int32_t middle=(start+end)/2;
01346
01347 if (output[middle]>elem)
01348 end=middle-1;
01349 else if (output[middle]<elem)
01350 start=middle+1;
01351 else
01352 return middle;
01353 }
01354
01355 return start;
01356 }
01357
01358
01359
01360 template <class T>
01361 static inline int32_t binary_search(T* output, int32_t size, T elem)
01362 {
01363 int32_t ind = binary_search_helper(output, size, elem);
01364 if (ind >= 0 && output[ind] == elem)
01365 return ind;
01366 return -1;
01367 }
01368
01369
01370
01371
01372
01373
01374
01375 template<class T>
01376 static inline int32_t binary_search(SGVector<T*> array, T* elem)
01377 {
01378 int32_t start=0;
01379 int32_t end=array.vlen-1;
01380
01381 if (array.vlen<1)
01382 return -1;
01383
01384 while (start<end)
01385 {
01386 int32_t middle=(start+end)/2;
01387
01388 if (*array.vector[middle]>*elem)
01389 end=middle-1;
01390 else if (*array.vector[middle]<*elem)
01391 start=middle+1;
01392 else
01393 {
01394 start=middle;
01395 break;
01396 }
01397 }
01398
01399 if (start>=0&&*array.vector[start]==*elem)
01400 return start;
01401
01402 return -1;
01403 }
01404
01405 template <class T>
01406 static int32_t binary_search_max_lower_equal(
01407 T* output, int32_t size, T elem)
01408 {
01409 int32_t ind = binary_search_helper(output, size, elem);
01410
01411 if (output[ind]<=elem)
01412 return ind;
01413 if (ind>0 && output[ind-1] <= elem)
01414 return ind-1;
01415 return -1;
01416 }
01417
01420 static float64_t Align(
01421 char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
01422
01423
01425
01428 static float64_t mutual_info(float64_t* p1, float64_t* p2, int32_t len);
01429
01432 static float64_t relative_entropy(
01433 float64_t* p, float64_t* q, int32_t len);
01434
01436 static float64_t entropy(float64_t* p, int32_t len);
01437
01439 inline static uint32_t get_seed()
01440 {
01441 return CMath::seed;
01442 }
01443
01445 inline static uint32_t get_log_range()
01446 {
01447 return CMath::LOGRANGE;
01448 }
01449
01450 #ifdef USE_LOGCACHE
01451
01452 inline static uint32_t get_log_accuracy()
01453 {
01454 return CMath::LOGACCURACY;
01455 }
01456 #endif
01457
01459 inline static int is_finite(double f)
01460 {
01461 #if defined(isfinite) && !defined(SUNOS)
01462 return isfinite(f);
01463 #else
01464 return finite(f);
01465 #endif
01466 }
01467
01469 inline static int is_infinity(double f)
01470 {
01471 #ifdef SUNOS
01472 if (fpclass(f) == FP_NINF || fpclass(f) == FP_PINF)
01473 return 1;
01474 else
01475 return 0;
01476 #else
01477 return isinf(f);
01478 #endif
01479 }
01480
01482 inline static int is_nan(double f)
01483 {
01484 #ifdef SUNOS
01485 return isnand(f);
01486 #else
01487 return isnan(f);
01488 #endif
01489 }
01490
01494 static SGVector<float64_t> fishers_exact_test_for_multiple_2x3_tables(SGMatrix<float64_t> tables);
01495
01499 static float64_t fishers_exact_test_for_2x3_table(SGMatrix<float64_t> table);
01500
01511 #ifdef USE_LOGCACHE
01512 static inline float64_t logarithmic_sum(float64_t p, float64_t q)
01513 {
01514 float64_t diff;
01515
01516 if (!CMath::is_finite(p))
01517 return q;
01518
01519 if (!CMath::is_finite(q))
01520 {
01521 SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
01522 return NAN;
01523 }
01524 diff = p - q;
01525 if (diff > 0)
01526 return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
01527 return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
01528 }
01529
01531 static void init_log_table();
01532
01534 static int32_t determine_logrange();
01535
01537 static int32_t determine_logaccuracy(int32_t range);
01538 #else
01539 static inline float64_t logarithmic_sum(
01540 float64_t p, float64_t q)
01541 {
01542 float64_t diff;
01543
01544 if (!CMath::is_finite(p))
01545 return q;
01546 if (!CMath::is_finite(q))
01547 return p;
01548 diff = p - q;
01549 if (diff > 0)
01550 return diff > LOGRANGE? p : p + log(1 + exp(-diff));
01551 return -diff > LOGRANGE? q : q + log(1 + exp(diff));
01552 }
01553 #endif
01554 #ifdef USE_LOGSUMARRAY
01555
01560 static inline float64_t logarithmic_sum_array(
01561 float64_t *p, int32_t len)
01562 {
01563 if (len<=2)
01564 {
01565 if (len==2)
01566 return logarithmic_sum(p[0],p[1]) ;
01567 if (len==1)
01568 return p[0];
01569 return -INFTY ;
01570 }
01571 else
01572 {
01573 register float64_t *pp=p ;
01574 if (len%2==1) pp++ ;
01575 for (register int32_t j=0; j < len>>1; j++)
01576 pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
01577 }
01578 return logarithmic_sum_array(p,len%2 + (len>>1)) ;
01579 }
01580 #endif //USE_LOGSUMARRAY
01581
01582
01584 inline virtual const char* get_name() const { return "Mathematics"; }
01585 public:
01588
01589 static const float64_t INFTY;
01590 static const float64_t ALMOST_INFTY;
01591
01593 static const float64_t ALMOST_NEG_INFTY;
01594
01596 static const float64_t PI;
01597
01599 static const float64_t MACHINE_EPSILON;
01600
01601
01602 static const float64_t MAX_REAL_NUMBER;
01603 static const float64_t MIN_REAL_NUMBER;
01604
01605 protected:
01607 static int32_t LOGRANGE;
01608
01610 static uint32_t seed;
01611 static char* rand_state;
01612
01613 #ifdef USE_LOGCACHE
01614
01616 static int32_t LOGACCURACY;
01618
01619 static float64_t* logtable;
01620 #endif
01621 };
01622
01623
01624 template <class T1,class T2>
01625 void* CMath::parallel_qsort_index(void* p)
01626 {
01627 struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
01628 T1* output=ps->output;
01629 T2* index=ps->index;
01630 uint32_t size=ps->size;
01631 int32_t* qsort_threads=ps->qsort_threads;
01632 int32_t sort_limit=ps->sort_limit;
01633 int32_t num_threads=ps->num_threads;
01634
01635 if (size==2)
01636 {
01637 if (output[0] > output [1])
01638 {
01639 swap(output[0], output[1]);
01640 swap(index[0], index[1]);
01641 }
01642 return NULL;
01643 }
01644
01645 T1 split=output[size/2];
01646
01647 int32_t left=0;
01648 int32_t right=size-1;
01649
01650 while (left<=right)
01651 {
01652 while (output[left] < split)
01653 left++;
01654 while (output[right] > split)
01655 right--;
01656
01657 if (left<=right)
01658 {
01659 swap(output[left], output[right]);
01660 swap(index[left], index[right]);
01661 left++;
01662 right--;
01663 }
01664 }
01665 bool lthread_start=false;
01666 bool rthread_start=false;
01667 pthread_t lthread;
01668 pthread_t rthread;
01669 struct thread_qsort<T1,T2> t1;
01670 struct thread_qsort<T1,T2> t2;
01671
01672 if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
01673 qsort_index(output,index,right+1);
01674 else if (right+1> 1)
01675 {
01676 (*qsort_threads)++;
01677 lthread_start=true;
01678 t1.output=output;
01679 t1.index=index;
01680 t1.size=right+1;
01681 t1.qsort_threads=qsort_threads;
01682 t1.sort_limit=sort_limit;
01683 t1.num_threads=num_threads;
01684 if (pthread_create(<hread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
01685 {
01686 lthread_start=false;
01687 (*qsort_threads)--;
01688 qsort_index(output,index,right+1);
01689 }
01690 }
01691
01692
01693 if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
01694 qsort_index(&output[left],&index[left], size-left);
01695 else if (size-left> 1)
01696 {
01697 (*qsort_threads)++;
01698 rthread_start=true;
01699 t2.output=&output[left];
01700 t2.index=&index[left];
01701 t2.size=size-left;
01702 t2.qsort_threads=qsort_threads;
01703 t2.sort_limit=sort_limit;
01704 t2.num_threads=num_threads;
01705 if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
01706 {
01707 rthread_start=false;
01708 (*qsort_threads)--;
01709 qsort_index(&output[left],&index[left], size-left);
01710 }
01711 }
01712
01713 if (lthread_start)
01714 {
01715 pthread_join(lthread, NULL);
01716 (*qsort_threads)--;
01717 }
01718
01719 if (rthread_start)
01720 {
01721 pthread_join(rthread, NULL);
01722 (*qsort_threads)--;
01723 }
01724
01725 return NULL;
01726 }
01727
01728 template <class T1,class T2>
01729 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
01730 {
01731 if (size==1)
01732 return;
01733
01734 if (size==2)
01735 {
01736 if (output[0] > output [1])
01737 {
01738 swap(output[0],output[1]);
01739 swap(index[0],index[1]);
01740 }
01741 return;
01742 }
01743
01744 T1 split=output[size/2];
01745
01746 int32_t left=0;
01747 int32_t right=size-1;
01748
01749 while (left<=right)
01750 {
01751 while (output[left] < split)
01752 left++;
01753 while (output[right] > split)
01754 right--;
01755
01756 if (left<=right)
01757 {
01758 swap(output[left],output[right]);
01759 swap(index[left],index[right]);
01760 left++;
01761 right--;
01762 }
01763 }
01764
01765 if (right+1> 1)
01766 qsort_index(output,index,right+1);
01767
01768 if (size-left> 1)
01769 qsort_index(&output[left],&index[left], size-left);
01770 }
01771
01772 template <class T1,class T2>
01773 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
01774 {
01775 if (size==2)
01776 {
01777 if (output[0] < output [1])
01778 {
01779 swap(output[0],output[1]);
01780 swap(index[0],index[1]);
01781 }
01782 return;
01783 }
01784
01785
01786 T1 split=output[size/2];
01787
01788 int32_t left=0;
01789 int32_t right=size-1;
01790
01791 while (left<=right)
01792 {
01793 while (output[left] > split)
01794 left++;
01795 while (output[right] < split)
01796 right--;
01797
01798 if (left<=right)
01799 {
01800 swap(output[left],output[right]);
01801 swap(index[left],index[right]);
01802 left++;
01803 right--;
01804 }
01805 }
01806
01807 if (right+1> 1)
01808 qsort_backward_index(output,index,right+1);
01809
01810 if (size-left> 1)
01811 qsort_backward_index(&output[left],&index[left], size-left);
01812 }
01813
01814 template <class T>
01815 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
01816 {
01817 if (6*n*size<13*size*CMath::log(size))
01818 for (int32_t i=0; i<n; i++)
01819 min(&output[i], &index[i], size-i) ;
01820 else
01821 qsort_index(output, index, size) ;
01822 }
01823
01824
01825 template <class T>
01826 void CMath::min(float64_t* output, T* index, int32_t size)
01827 {
01828 if (size<=1)
01829 return;
01830 float64_t min_elem=output[0];
01831 int32_t min_index=0;
01832 for (int32_t i=1; i<size; i++)
01833 {
01834 if (output[i]<min_elem)
01835 {
01836 min_index=i;
01837 min_elem=output[i];
01838 }
01839 }
01840 swap(output[0], output[min_index]);
01841 swap(index[0], index[min_index]);
01842 }
01843 }
01844 #endif