00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MATHMATICS_H_
00014 #define __MATHMATICS_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
00032 #ifdef _WIN32
00033 #include <float.h>
00034 #ifndef finite
00035 #define finite _finite
00036 #endif
00037
00038 #ifndef isnan
00039 #define isnan _isnan
00040 #endif
00041 #endif
00042
00043
00044 #ifdef CYGWIN
00045 #ifndef finite
00046 #include <ieeefp.h>
00047 #endif
00048 #endif
00049
00050 #ifndef NAN
00051 #include <stdlib.h>
00052 #define NAN (strtod("NAN",NULL))
00053 #endif
00054
00055 #ifdef SUNOS
00056 extern "C" int finite(double);
00057 #endif
00058
00059
00060 #define RNG_SEED_SIZE 256
00061
00062
00063 #define RADIX_STACK_SIZE 512
00064
00065
00066 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
00067 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
00068
00070 template <class T> struct radix_stack_t
00071 {
00073 T *sa;
00075 size_t sn;
00077 uint16_t si;
00078 };
00079
00081
00083 struct thread_qsort
00084 {
00086 float64_t* output;
00088 int32_t* index;
00090 int32_t size;
00091
00093 int32_t* qsort_threads;
00095 int32_t sort_limit;
00097 int32_t num_threads;
00098 };
00099
00103 class CMath : public CSGObject
00104 {
00105 public:
00109
00110 CMath();
00111
00113 virtual ~CMath();
00115
00119
00121
00122 template <class T>
00123 static inline T min(T a, T b)
00124 {
00125 return ((a)<=(b))?(a):(b);
00126 }
00127
00129 template <class T>
00130 static inline T max(T a, T b)
00131 {
00132 return ((a)>=(b))?(a):(b);
00133 }
00134
00136 template <class T>
00137 static inline T clamp(T value, T lb, T ub)
00138 {
00139 if (value<=lb)
00140 return lb;
00141 else if (value>=ub)
00142 return ub;
00143 else
00144 return value;
00145 }
00146
00148 template <class T>
00149 static inline T abs(T a)
00150 {
00151
00152
00153 if (a==0)
00154 return 0;
00155 else if (a>0)
00156 return a;
00157 else
00158 return -a;
00159 }
00161
00164
00166 static uint32_t crc32(uint8_t *data, int32_t len);
00167
00168 static inline float64_t round(float64_t d)
00169 {
00170 return ::floor(d+0.5);
00171 }
00172
00173 static inline float64_t floor(float64_t d)
00174 {
00175 return ::floor(d);
00176 }
00177
00178 static inline float64_t ceil(float64_t d)
00179 {
00180 return ::ceil(d);
00181 }
00182
00184 template <class T>
00185 static inline T sign(T a)
00186 {
00187 if (a==0)
00188 return 0;
00189 else return (a<0) ? (-1) : (+1);
00190 }
00191
00193 template <class T>
00194 static inline void swap(T & a,T &b)
00195 {
00196 T c=a;
00197 a=b;
00198 b=c;
00199 }
00200
00202 template <class T>
00203 static inline T sq(T x)
00204 {
00205 return x*x;
00206 }
00207
00209 static inline float32_t sqrt(float32_t x)
00210 {
00211 return ::sqrtf(x);
00212 }
00213
00215 static inline float64_t sqrt(float64_t x)
00216 {
00217 return ::sqrt(x);
00218 }
00219
00221 static inline float128_t sqrt(float128_t x)
00222 {
00223
00224
00225 #ifdef HAVE_SQRTL
00226 return ::sqrtl(x);
00227 #else
00228 return ::sqrt(x);
00229 #endif
00230 }
00231
00232
00234 static inline float128_t powl(float128_t x, float128_t n)
00235 {
00236
00237
00238 #ifdef HAVE_POWL
00239 return ::powl((long double) x, (long double) n);
00240 #else
00241 return ::pow((double) x, (double) n);
00242 #endif
00243 }
00244
00245 static inline int32_t pow(int32_t x, int32_t n)
00246 {
00247 ASSERT(n>=0);
00248 int32_t result=1;
00249 while (n--)
00250 result*=x;
00251
00252 return result;
00253 }
00254
00255 static inline float64_t pow(float64_t x, int32_t n)
00256 {
00257 ASSERT(n>=0);
00258 float64_t result=1;
00259 while (n--)
00260 result*=x;
00261
00262 return result;
00263 }
00264
00265 static inline float64_t pow(float64_t x, float64_t n)
00266 {
00267 return ::pow((double) x, (double) n);
00268 }
00269
00270 static inline float64_t log10(float64_t v)
00271 {
00272 return ::log(v)/::log(10.0);
00273 }
00274
00275 #ifdef HAVE_LOG2
00276 static inline float64_t log2(float64_t v)
00277 {
00278 return ::log2(v);
00279 }
00280 #else
00281 static inline float64_t log2(float64_t v)
00282 {
00283 return ::log(v)/::log(2.0);
00284 }
00285 #endif //HAVE_LOG2
00286
00287 static inline float64_t log(float64_t v)
00288 {
00289 return ::log(v);
00290 }
00291
00292 template <class T>
00293 static void transpose_matrix(
00294 T*& matrix, int32_t& num_feat, int32_t& num_vec)
00295 {
00296 T* transposed=new T[num_vec*num_feat];
00297 for (int32_t i=0; i<num_vec; i++)
00298 {
00299 for (int32_t j=0; j<num_feat; j++)
00300 transposed[i+j*num_vec]=matrix[i*num_feat+j];
00301 }
00302
00303 delete[] matrix;
00304 matrix=transposed;
00305
00306 CMath::swap(num_feat, num_vec);
00307 }
00308
00309 #ifdef HAVE_LAPACK
00312 static float64_t* pinv(
00313 float64_t* matrix, int32_t rows, int32_t cols,
00314 float64_t* target=NULL);
00315
00316
00317
00318
00319 static inline void dgemm(
00320 double alpha, const double* A, int rows, int cols,
00321 CBLAS_TRANSPOSE transposeA, double *B, int cols_B,
00322 CBLAS_TRANSPOSE transposeB, double beta, double *C)
00323 {
00324 cblas_dgemm(CblasColMajor, transposeA, transposeB, rows, cols, cols_B,
00325 alpha, A, cols, B, cols_B, beta, C, cols);
00326 }
00327
00328
00329 static inline void dgemv(
00330 double alpha, const double *A, int rows, int cols,
00331 const CBLAS_TRANSPOSE transposeA, const double* X, double beta,
00332 double* Y)
00333 {
00334 cblas_dgemv(CblasColMajor, transposeA,
00335 rows, cols, alpha, A, cols,
00336 X, 1, beta, Y, 1);
00337 }
00338 #endif
00339
00340 static inline int64_t factorial(int32_t n)
00341 {
00342 int64_t res=1;
00343 for (int i=2; i<=n; i++)
00344 res*=i ;
00345 return res ;
00346 }
00347
00348 static void init_random(uint32_t initseed=0)
00349 {
00350 if (initseed==0)
00351 {
00352 struct timeval tv;
00353 gettimeofday(&tv, NULL);
00354 seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
00355 }
00356 else
00357 seed=initseed;
00358 #if !defined(CYGWIN) && !defined(__INTERIX)
00359
00360
00361 initstate(seed, CMath::rand_state, RNG_SEED_SIZE);
00362 #endif
00363 }
00364
00365 static inline int64_t random()
00366 {
00367 #if defined(CYGWIN) || defined(__INTERIX)
00368 return rand();
00369 #else
00370 return ::random();
00371 #endif
00372 }
00373
00374 static inline int32_t random(int32_t min_value, int32_t max_value)
00375 {
00376 int32_t ret = min_value + (int32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00377 ASSERT(ret>=min_value && ret<=max_value);
00378 return ret ;
00379 }
00380
00381 static inline float32_t random(float32_t min_value, float32_t max_value)
00382 {
00383 float32_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00384
00385 if (ret<min_value || ret>max_value)
00386 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00387 ASSERT(ret>=min_value && ret<=max_value);
00388 return ret;
00389 }
00390
00391 static inline float64_t random(float64_t min_value, float64_t max_value)
00392 {
00393 float64_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00394
00395 if (ret<min_value || ret>max_value)
00396 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00397 ASSERT(ret>=min_value && ret<=max_value);
00398 return ret;
00399 }
00400
00401 template <class T>
00402 static void fill_vector(T* vec, int32_t len, T value)
00403 {
00404 for (int32_t i=0; i<len; i++)
00405 vec[i]=value;
00406 }
00407 template <class T>
00408 static void range_fill_vector(T* vec, int32_t len, T start=0)
00409 {
00410 for (int32_t i=0; i<len; i++)
00411 vec[i]=i+start;
00412 }
00413
00414 template <class T>
00415 static void random_vector(T* vec, int32_t len, T min_value, T max_value)
00416 {
00417 for (int32_t i=0; i<len; i++)
00418 vec[i]=CMath::random(min_value, max_value);
00419 }
00420
00421 static inline int32_t* randperm(int32_t n)
00422 {
00423 int32_t* perm = new int32_t[n];
00424
00425 if (!perm)
00426 return NULL;
00427 for (int32_t i = 0; i < n; i++)
00428 perm[i] = i;
00429 for (int32_t i = 0; i < n; i++)
00430 swap(perm[random(0, n - 1)], perm[i]);
00431 return perm;
00432 }
00433
00434 static inline int64_t nchoosek(int32_t n, int32_t k)
00435 {
00436 int64_t res=1;
00437
00438 for (int32_t i=n-k+1; i<=n; i++)
00439 res*=i;
00440
00441 return res/factorial(k);
00442 }
00443
00445 template <class T>
00446 static inline void vec1_plus_scalar_times_vec2(T* vec1,
00447 T scalar, const T* vec2, int32_t n)
00448 {
00449 for (int32_t i=0; i<n; i++)
00450 vec1[i]+=scalar*vec2[i];
00451 }
00452
00454 static inline float64_t dot(const float64_t* v1, const float64_t* v2, int32_t n)
00455 {
00456 float64_t r=0;
00457 #ifdef HAVE_LAPACK
00458 int32_t skip=1;
00459 r = cblas_ddot(n, v1, skip, v2, skip);
00460 #else
00461 for (int32_t i=0; i<n; i++)
00462 r+=v1[i]*v2[i];
00463 #endif
00464 return r;
00465 }
00466
00468 static inline float32_t dot(
00469 const float32_t* v1, const float32_t* v2, int32_t n)
00470 {
00471 float64_t r=0;
00472 #ifdef HAVE_LAPACK
00473 int32_t skip=1;
00474 r = cblas_sdot(n, v1, skip, v2, skip);
00475 #else
00476 for (int32_t i=0; i<n; i++)
00477 r+=v1[i]*v2[i];
00478 #endif
00479 return r;
00480 }
00481
00483 static inline float64_t dot(
00484 const int64_t* v1, const int64_t* v2, int32_t n)
00485 {
00486 float64_t r=0;
00487 for (int32_t i=0; i<n; i++)
00488 r+=((float64_t) v1[i])*v2[i];
00489
00490 return r;
00491 }
00492
00494 static inline float64_t dot(
00495 const uint16_t* v1, const uint16_t* v2, int32_t n)
00496 {
00497 float64_t r=0;
00498 for (int32_t i=0; i<n; i++)
00499 r+=((float64_t) v1[i])*v2[i];
00500
00501 return r;
00502 }
00503
00505 static inline float64_t dot(
00506 const char* v1, const char* v2, int32_t n)
00507 {
00508 float64_t r=0;
00509 for (int32_t i=0; i<n; i++)
00510 r+=((float64_t) v1[i])*v2[i];
00511
00512 return r;
00513 }
00514
00516 static inline float64_t dot(
00517 const uint8_t* v1, const uint8_t* v2, int32_t n)
00518 {
00519 float64_t r=0;
00520 for (int32_t i=0; i<n; i++)
00521 r+=((float64_t) v1[i])*v2[i];
00522
00523 return r;
00524 }
00525
00527 static inline float64_t dot(
00528 const float64_t* v1, const char* v2, int32_t n)
00529 {
00530 float64_t r=0;
00531 for (int32_t i=0; i<n; i++)
00532 r+=((float64_t) v1[i])*v2[i];
00533
00534 return r;
00535 }
00536
00538 template <class T>
00539 static inline void add(
00540 T* target, T alpha, const T* v1, T beta, const T* v2,
00541 int32_t len)
00542 {
00543 for (int32_t i=0; i<len; i++)
00544 target[i]=alpha*v1[i]+beta*v2[i];
00545 }
00546
00548 template <class T>
00549 static inline void add_scalar(T alpha, T* vec, int32_t len)
00550 {
00551 for (int32_t i=0; i<len; i++)
00552 vec[i]+=alpha;
00553 }
00554
00556 template <class T>
00557 static inline void scale_vector(T alpha, T* vec, int32_t len)
00558 {
00559 for (int32_t i=0; i<len; i++)
00560 vec[i]*=alpha;
00561 }
00562
00564 template <class T>
00565 static inline T sum(T* vec, int32_t len)
00566 {
00567 T result=0;
00568 for (int32_t i=0; i<len; i++)
00569 result+=vec[i];
00570
00571 return result;
00572 }
00573
00575 template <class T>
00576 static inline T max(T* vec, int32_t len)
00577 {
00578 ASSERT(len>0);
00579 T maxv=vec[0];
00580
00581 for (int32_t i=1; i<len; i++)
00582 maxv=CMath::max(vec[i], maxv);
00583
00584 return maxv;
00585 }
00586
00588 template <class T>
00589 static inline T sum_abs(T* vec, int32_t len)
00590 {
00591 T result=0;
00592 for (int32_t i=0; i<len; i++)
00593 result+=CMath::abs(vec[i]);
00594
00595 return result;
00596 }
00597
00598 static inline float64_t mean(float64_t* vec, int32_t len)
00599 {
00600 ASSERT(vec);
00601 ASSERT(len>0);
00602
00603 float64_t mean=0;
00604 for (int32_t i=0; i<len; i++)
00605 mean+=vec[i];
00606 return mean/len;
00607 }
00608
00609 static inline float64_t trace(
00610 float64_t* mat, int32_t cols, int32_t rows)
00611 {
00612 float64_t trace=0;
00613 for (int32_t i=0; i<rows; i++)
00614 trace+=mat[i*cols+i];
00615 return trace;
00616 }
00617
00621 static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
00622 static void sort(float64_t *a, int32_t*idx, int32_t N);
00623
00624
00625
00626
00627
00628
00629
00631 template <class T>
00632 inline static void radix_sort(T* array, size_t size)
00633 {
00634 radix_sort_helper(array,size,0);
00635 }
00636
00637 template <class T>
00638 static inline uint8_t byte(T word, uint16_t p)
00639 {
00640 return (word >> (sizeof(T)-p-1) * 8) & 0xff;
00641 }
00642
00643 template <class T>
00644 static void radix_sort_helper(T* array, size_t size, uint16_t i)
00645 {
00646 static size_t count[256], nc, cmin;
00647 T *ak;
00648 uint8_t c=0;
00649 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
00650 T *an, *aj, *pile[256];
00651 size_t *cp, cmax;
00652
00653
00654 sp = s;
00655 radix_push(array, size, i);
00656
00657
00658 while (sp>s) {
00659 radix_pop(array, size, i);
00660 an = array + size;
00661
00662
00663 if (nc == 0) {
00664 cmin = 0xff;
00665 for (ak = array; ak < an; ak++) {
00666 c = byte(*ak, i);
00667 count[c]++;
00668 if (count[c] == 1) {
00669
00670 if (c < cmin)
00671 cmin = c;
00672 nc++;
00673 }
00674 }
00675
00676
00677 if (sp + nc > s + RADIX_STACK_SIZE) {
00678 radix_sort_helper(array, size, i);
00679 continue;
00680 }
00681 }
00682
00683
00684
00685
00686
00687
00688
00689 olds = bigs = sp;
00690 cmax = 2;
00691 ak = array;
00692 pile[0xff] = an;
00693 for (cp = count + cmin; nc > 0; cp++) {
00694
00695 while (*cp == 0)
00696 cp++;
00697
00698 if (*cp > 1) {
00699
00700 if (*cp > cmax) {
00701 cmax = *cp;
00702 bigs = sp;
00703 }
00704
00705 if (i < sizeof(T)-1)
00706 radix_push(ak, *cp, (uint16_t) (i + 1));
00707 }
00708 pile[cp - count] = ak += *cp;
00709 nc--;
00710 }
00711
00712
00713 swap(*olds, *bigs);
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727 aj = array;
00728 while (aj<an)
00729 {
00730 T r;
00731
00732 for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
00733 swap(*ak, r);
00734
00735 *aj = r;
00736 aj += count[c];
00737 count[c] = 0;
00738 }
00739 }
00740 }
00743 template <class T>
00744 static void insertion_sort(T* output, int32_t size)
00745 {
00746 for (int32_t i=0; i<size-1; i++)
00747 {
00748 int32_t j=i-1;
00749 T value=output[i];
00750 while (j >= 0 && output[j] > value)
00751 {
00752 output[j+1] = output[j];
00753 j--;
00754 }
00755 output[j+1]=value;
00756 }
00757 }
00758
00759
00762
00763
00764 template <class T>
00765 static void qsort(T* output, int32_t size)
00766 {
00767 if (size==2)
00768 {
00769 if (output[0] > output [1])
00770 swap(output[0],output[1]);
00771 return;
00772 }
00773 T split=output[random(0,size-1)];
00774
00775 int32_t left=0;
00776 int32_t right=size-1;
00777
00778 while (left<=right)
00779 {
00780 while (output[left] < split)
00781 left++;
00782 while (output[right] > split)
00783 right--;
00784
00785 if (left<=right)
00786 {
00787 swap(output[left],output[right]);
00788 left++;
00789 right--;
00790 }
00791 }
00792
00793 if (right+1> 1)
00794 qsort(output,right+1);
00795
00796 if (size-left> 1)
00797 qsort(&output[left],size-left);
00798 }
00799
00801 template <class T> static void display_vector(
00802 T* vector, int32_t n, const char* name="vector");
00803
00805 template <class T> static void display_matrix(
00806 T* matrix, int32_t rows, int32_t cols, const char* name="matrix");
00807
00813 template <class T1,class T2>
00814 static void qsort_index(T1* output, T2* index, uint32_t size);
00815
00816 template <class T1,class T2>
00817 static void* parallel_qsort_index(void* p);
00818
00824 template <class T1,class T2>
00825 static void qsort_backward_index(
00826 T1* output, T2* index, int32_t size);
00827
00828
00829
00830 template <class T>
00831 static void min(float64_t* output, T* index, int32_t size);
00832
00833
00834
00835 template <class T>
00836 static void nmin(
00837 float64_t* output, T* index, int32_t size, int32_t n);
00838
00839
00840
00841 template <class T>
00842 static int32_t unique(T* output, int32_t size)
00843 {
00844 qsort(output, size);
00845 int32_t i,j=0 ;
00846 for (i=0; i<size; i++)
00847 if (i==0 || output[i]!=output[i-1])
00848 output[j++]=output[i];
00849 return j ;
00850 }
00851
00852
00853
00854 template <class T>
00855 static int32_t binary_search_helper(T* output, int32_t size, T elem)
00856 {
00857 int32_t start=0;
00858 int32_t end=size-1;
00859
00860 if (size<1)
00861 return -1;
00862
00863 while (start<end)
00864 {
00865 int32_t middle=(start+end)/2;
00866
00867 if (output[middle]>elem)
00868 end=middle-1;
00869 else if (output[middle]<elem)
00870 start=middle+1;
00871 else
00872 return middle;
00873 }
00874
00875 if (output[start]==elem)
00876 return start;
00877 else
00878 return -1;
00879 }
00880
00881
00882
00883 template <class T>
00884 static inline int32_t binary_search(T* output, int32_t size, T elem)
00885 {
00886 int32_t ind = binary_search_helper(output, size, elem);
00887 if (ind >= 0 && output[ind] == elem)
00888 return ind;
00889 return -1;
00890 }
00891
00892
00893
00894
00895
00896 template <class T>
00897 static int32_t binary_search_max_lower_equal(
00898 T* output, int32_t size, T elem)
00899 {
00900 int32_t ind = binary_search_helper(output, size, elem);
00901
00902 if (output[ind]<=elem)
00903 return ind;
00904 if (ind>0 && output[ind-1] <= elem)
00905 return ind-1;
00906 return -1;
00907 }
00908
00911 static float64_t Align(
00912 char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
00913
00918 static int32_t calcroc(
00919 float64_t* fp, float64_t* tp, float64_t* output, int32_t* label,
00920 int32_t& size, int32_t& possize, int32_t& negsize,
00921 float64_t& tresh, FILE* rocfile);
00923
00926 static float64_t mutual_info(float64_t* p1, float64_t* p2, int32_t len);
00927
00930 static float64_t relative_entropy(
00931 float64_t* p, float64_t* q, int32_t len);
00932
00934 static float64_t entropy(float64_t* p, int32_t len);
00935
00937 inline static uint32_t get_seed()
00938 {
00939 return CMath::seed;
00940 }
00941
00942
00953 #ifdef USE_LOGCACHE
00954 static inline float64_t logarithmic_sum(float64_t p, float64_t q)
00955 {
00956 float64_t diff;
00957
00958 if (!finite(p))
00959 return q;
00960 if (!finite(q))
00961 {
00962 SG_SWARNING(("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
00963 return NAN;
00964 }
00965 diff = p - q;
00966 if (diff > 0)
00967 return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
00968 return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
00969 }
00970
00972 static void init_log_table();
00973
00975 static int32_t determine_logrange();
00976
00978 static int32_t determine_logaccuracy(int32_t range);
00979 #else
00980 static inline float64_t logarithmic_sum(
00981 float64_t p, float64_t q)
00982 {
00983 float64_t diff;
00984
00985 if (!finite(p))
00986 return q;
00987 if (!finite(q))
00988 return p;
00989 diff = p - q;
00990 if (diff > 0)
00991 return diff > LOGRANGE? p : p + log(1 + exp(-diff));
00992 return -diff > LOGRANGE? q : q + log(1 + exp(diff));
00993 }
00994 #endif
00995 #ifdef LOG_SUM_ARRAY
00996
01001 static inline float64_t logarithmic_sum_array(
01002 float64_t *p, int32_t len)
01003 {
01004 if (len<=2)
01005 {
01006 if (len==2)
01007 return logarithmic_sum(p[0],p[1]) ;
01008 if (len==1)
01009 return p[0];
01010 return -INFTY ;
01011 }
01012 else
01013 {
01014 register float64_t *pp=p ;
01015 if (len%2==1) pp++ ;
01016 for (register int32_t j=0; j < len>>1; j++)
01017 pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
01018 }
01019 return logarithmic_sum_array(p,len%2+len>>1) ;
01020 }
01021 #endif
01022
01023 public:
01026
01027 static const float64_t INFTY;
01028 static const float64_t ALMOST_INFTY;
01029
01031 static const float64_t ALMOST_NEG_INFTY;
01032
01034 static int32_t LOGRANGE;
01035
01037 static uint32_t seed;
01038
01039 #ifdef USE_LOGCACHE
01040
01042 static int32_t LOGACCURACY;
01044 protected:
01046 static float64_t* logtable;
01047 #endif
01048 static char* rand_state;
01049 };
01050
01051 template <class T1,class T2>
01052 void* CMath::parallel_qsort_index(void* p)
01053 {
01054 struct thread_qsort* ps=(thread_qsort*) p;
01055 T1* output=ps->output;
01056 T2* index=ps->index;
01057 int32_t size=ps->size;
01058 int32_t* qsort_threads=ps->qsort_threads;
01059 int32_t sort_limit=ps->sort_limit;
01060 int32_t num_threads=ps->num_threads;
01061
01062 if (size==2)
01063 {
01064 if (output[0] > output [1])
01065 {
01066 swap(output[0], output[1]);
01067 swap(index[0], index[1]);
01068 }
01069 return NULL;
01070 }
01071
01072 float64_t split=output[size/2];
01073
01074 int32_t left=0;
01075 int32_t right=size-1;
01076
01077 while (left<=right)
01078 {
01079 while (output[left] < split)
01080 left++;
01081 while (output[right] > split)
01082 right--;
01083
01084 if (left<=right)
01085 {
01086 swap(output[left], output[right]);
01087 swap(index[left], index[right]);
01088 left++;
01089 right--;
01090 }
01091 }
01092 bool lthread_start=false;
01093 bool rthread_start=false;
01094 pthread_t lthread;
01095 pthread_t rthread;
01096 struct thread_qsort t1;
01097 struct thread_qsort t2;
01098
01099 if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
01100 qsort_index(output,index,right+1);
01101 else if (right+1> 1)
01102 {
01103 *qsort_threads++;
01104 lthread_start=true;
01105 t1.output=output;
01106 t1.index=index;
01107 t1.size=right+1;
01108 if (pthread_create(<hread, NULL, parallel_qsort_index<T1,T2>, (void*) &t1) != 0)
01109 {
01110 lthread_start=false;
01111 *qsort_threads--;
01112 qsort_index(output,index,right+1);
01113 }
01114 }
01115
01116
01117 if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
01118 qsort_index(&output[left],&index[left], size-left);
01119 else if (size-left> 1)
01120 {
01121 *qsort_threads++;
01122 rthread_start=true;
01123 t2.output=&output[left];
01124 t2.index=&index[left];
01125 t2.size=size-left;
01126 if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, (void*)&t2) != 0)
01127 {
01128 rthread_start=false;
01129 *qsort_threads--;
01130 qsort_index(&output[left],&index[left], size-left);
01131 }
01132 }
01133
01134 if (lthread_start)
01135 {
01136 pthread_join(lthread, NULL);
01137 *qsort_threads--;
01138 }
01139
01140 if (rthread_start)
01141 {
01142 pthread_join(rthread, NULL);
01143 *qsort_threads--;
01144 }
01145
01146 return NULL;
01147 }
01148
01149
01150
01151 template <class T1,class T2>
01152 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
01153 {
01154 if (size==2)
01155 {
01156 if (output[0] > output [1])
01157 {
01158 swap(output[0],output[1]);
01159 swap(index[0],index[1]);
01160 }
01161 return;
01162 }
01163 T1 split=output[random(0,size-1)];
01164
01165 int32_t left=0;
01166 int32_t right=size-1;
01167
01168 while (left<=right)
01169 {
01170 while (output[left] < split)
01171 left++;
01172 while (output[right] > split)
01173 right--;
01174
01175 if (left<=right)
01176 {
01177 swap(output[left],output[right]);
01178 swap(index[left],index[right]);
01179 left++;
01180 right--;
01181 }
01182 }
01183
01184 if (right+1> 1)
01185 qsort_index(output,index,right+1);
01186
01187 if (size-left> 1)
01188 qsort_index(&output[left],&index[left], size-left);
01189 }
01190
01191 template <class T1,class T2>
01192 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
01193 {
01194 if (size==2)
01195 {
01196 if (output[0] < output [1])
01197 {
01198 swap(output[0],output[1]);
01199 swap(index[0],index[1]);
01200 }
01201 return;
01202 }
01203
01204 T1 split=output[random(0,size-1)];
01205
01206 int32_t left=0;
01207 int32_t right=size-1;
01208
01209 while (left<=right)
01210 {
01211 while (output[left] > split)
01212 left++;
01213 while (output[right] < split)
01214 right--;
01215
01216 if (left<=right)
01217 {
01218 swap(output[left],output[right]);
01219 swap(index[left],index[right]);
01220 left++;
01221 right--;
01222 }
01223 }
01224
01225 if (right+1> 1)
01226 qsort_backward_index(output,index,right+1);
01227
01228 if (size-left> 1)
01229 qsort_backward_index(&output[left],&index[left], size-left);
01230 }
01231
01232 template <class T>
01233 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
01234 {
01235 if (6*n*size<13*size*CMath::log(size))
01236 for (int32_t i=0; i<n; i++)
01237 min(&output[i], &index[i], size-i) ;
01238 else
01239 qsort_index(output, index, size) ;
01240 }
01241
01242
01243 template <class T>
01244 void CMath::min(float64_t* output, T* index, int32_t size)
01245 {
01246 if (size<=1)
01247 return;
01248 float64_t min_elem=output[0];
01249 int32_t min_index=0;
01250 for (int32_t i=1; i<size; i++)
01251 {
01252 if (output[i]<min_elem)
01253 {
01254 min_index=i;
01255 min_elem=output[i];
01256 }
01257 }
01258 swap(output[0], output[min_index]);
01259 swap(index[0], index[min_index]);
01260 }
01261 #endif