00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef __MATHEMATICS_H_
00019 #define __MATHEMATICS_H_
00020
00021 #include <shogun/base/SGObject.h>
00022 #include <shogun/lib/common.h>
00023 #include <shogun/io/SGIO.h>
00024 #include <shogun/base/Parallel.h>
00025
00026 #include <math.h>
00027 #include <stdio.h>
00028 #include <float.h>
00029 #include <pthread.h>
00030 #include <unistd.h>
00031 #include <sys/types.h>
00032 #include <sys/time.h>
00033 #include <time.h>
00034
00035 #ifdef SUNOS
00036 #include <ieeefp.h>
00037 #endif
00038
00040 #ifdef log2
00041 #define cygwin_log2 log2
00042 #undef log2
00043 #endif
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 {
00132 class CSGObject;
00135 class CMath : public CSGObject
00136 {
00137 public:
00141
00142 CMath();
00143
00145 virtual ~CMath();
00147
00151
00153
00154 template <class T>
00155 static inline T min(T a, T b)
00156 {
00157 return (a<=b) ? a : b;
00158 }
00159
00161 template <class T>
00162 static inline T max(T a, T b)
00163 {
00164 return (a>=b) ? a : b;
00165 }
00166
00168 template <class T>
00169 static inline T clamp(T value, T lb, T ub)
00170 {
00171 if (value<=lb)
00172 return lb;
00173 else if (value>=ub)
00174 return ub;
00175 else
00176 return value;
00177 }
00178
00180 template <class T>
00181 static inline T abs(T a)
00182 {
00183
00184
00185 if (a==0)
00186 return 0;
00187 else if (a>0)
00188 return a;
00189 else
00190 return -a;
00191 }
00193
00196
00197 static inline float64_t round(float64_t d)
00198 {
00199 return ::floor(d+0.5);
00200 }
00201
00202 static inline float64_t floor(float64_t d)
00203 {
00204 return ::floor(d);
00205 }
00206
00207 static inline float64_t ceil(float64_t d)
00208 {
00209 return ::ceil(d);
00210 }
00211
00213 template <class T>
00214 static inline T sign(T a)
00215 {
00216 if (a==0)
00217 return 0;
00218 else return (a<0) ? (-1) : (+1);
00219 }
00220
00222 template <class T>
00223 static inline void swap(T &a,T &b)
00224 {
00225
00226 register T c=a;
00227 a=b;
00228 b=c;
00229 }
00230
00232 template <class T>
00233 static inline T sq(T x)
00234 {
00235 return x*x;
00236 }
00237
00239 static inline float32_t sqrt(float32_t x)
00240 {
00241 return ::sqrtf(x);
00242 }
00243
00245 static inline float64_t sqrt(float64_t x)
00246 {
00247 return ::sqrt(x);
00248 }
00249
00251 static inline floatmax_t sqrt(floatmax_t x)
00252 {
00253
00254
00255 #ifdef HAVE_SQRTL
00256 return ::sqrtl(x);
00257 #else
00258 return ::sqrt(x);
00259 #endif
00260 }
00261
00263 static inline float32_t invsqrt(float32_t x)
00264 {
00265 union float_to_int
00266 {
00267 float32_t f;
00268 int32_t i;
00269 };
00270
00271 float_to_int tmp;
00272 tmp.f=x;
00273
00274 float32_t xhalf = 0.5f * x;
00275
00276
00277 tmp.i = 0x5f3759d5 - (tmp.i >> 1);
00278 x = tmp.f;
00279 x = x*(1.5f - xhalf*x*x);
00280 return x;
00281 }
00282
00284 static inline floatmax_t powl(floatmax_t x, floatmax_t n)
00285 {
00286
00287
00288 #ifdef HAVE_POWL
00289 return ::powl((long double) x, (long double) n);
00290 #else
00291 return ::pow((double) x, (double) n);
00292 #endif
00293 }
00294
00295 static inline int32_t pow(int32_t x, int32_t n)
00296 {
00297 ASSERT(n>=0);
00298 int32_t result=1;
00299 while (n--)
00300 result*=x;
00301
00302 return result;
00303 }
00304
00305 static inline float64_t pow(float64_t x, int32_t n)
00306 {
00307 if (n>=0)
00308 {
00309 float64_t result=1;
00310 while (n--)
00311 result*=x;
00312
00313 return result;
00314 }
00315 else
00316 return ::pow((double)x, (double)n);
00317 }
00318
00319 static inline float64_t pow(float64_t x, float64_t n)
00320 {
00321 return ::pow((double) x, (double) n);
00322 }
00323
00324 static inline float64_t exp(float64_t x)
00325 {
00326 return ::exp((double) x);
00327 }
00328
00330 static inline float64_t atan(float64_t x)
00331 {
00332 return ::atan((double) x);
00333 }
00334
00335 static inline float64_t log10(float64_t v)
00336 {
00337 return ::log(v)/::log(10.0);
00338 }
00339
00340 static inline float64_t log2(float64_t v)
00341 {
00342 #ifdef HAVE_LOG2
00343 return ::log2(v);
00344 #else
00345 return ::log(v)/::log(2.0);
00346 #endif //HAVE_LOG2
00347 }
00348
00349 static inline float64_t log(float64_t v)
00350 {
00351 return ::log(v);
00352 }
00353
00354 static inline float64_t sin(float64_t x)
00355 {
00356 return ::sin(x);
00357 }
00358
00359 static inline float64_t cos(float64_t x)
00360 {
00361 return ::cos(x);
00362 }
00363
00364 static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
00365 {
00366 ASSERT(len>0 && xy);
00367
00368 float64_t area = 0.0;
00369
00370 if (!reversed)
00371 {
00372 for (int i=1; i<len; i++)
00373 area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
00374 }
00375 else
00376 {
00377 for (int i=1; i<len; i++)
00378 area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
00379 }
00380
00381 return area;
00382 }
00383
00384 static inline int64_t factorial(int32_t n)
00385 {
00386 int64_t res=1;
00387 for (int i=2; i<=n; i++)
00388 res*=i ;
00389 return res ;
00390 }
00391
00392 static void init_random(uint32_t initseed=0)
00393 {
00394 if (initseed==0)
00395 {
00396 struct timeval tv;
00397 gettimeofday(&tv, NULL);
00398 seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
00399 }
00400 else
00401 seed=initseed;
00402 #if !defined(CYGWIN) && !defined(__INTERIX)
00403
00404
00405 initstate(seed, CMath::rand_state, RNG_SEED_SIZE);
00406 #endif
00407 }
00408
00409 static inline int64_t random()
00410 {
00411 #if defined(CYGWIN) || defined(__INTERIX)
00412 return rand();
00413 #else
00414 return ::random();
00415 #endif
00416 }
00417
00418 static inline uint64_t random(uint64_t min_value, uint64_t max_value)
00419 {
00420 uint64_t ret = min_value + (uint64_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00421 ASSERT(ret>=min_value && ret<=max_value);
00422 return ret ;
00423 }
00424
00425 static inline int64_t random(int64_t min_value, int64_t max_value)
00426 {
00427 int64_t ret = min_value + (int64_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00428 ASSERT(ret>=min_value && ret<=max_value);
00429 return ret ;
00430 }
00431
00432 static inline uint32_t random(uint32_t min_value, uint32_t max_value)
00433 {
00434 uint32_t ret = min_value + (uint32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00435 ASSERT(ret>=min_value && ret<=max_value);
00436 return ret ;
00437 }
00438
00439 static inline int32_t random(int32_t min_value, int32_t max_value)
00440 {
00441 int32_t ret = min_value + (int32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00442 ASSERT(ret>=min_value && ret<=max_value);
00443 return ret ;
00444 }
00445
00446 static inline float32_t random(float32_t min_value, float32_t max_value)
00447 {
00448 float32_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00449
00450 if (ret<min_value || ret>max_value)
00451 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00452 ASSERT(ret>=min_value && ret<=max_value);
00453 return ret;
00454 }
00455
00456 static inline float64_t random(float64_t min_value, float64_t max_value)
00457 {
00458 float64_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00459
00460 if (ret<min_value || ret>max_value)
00461 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00462 ASSERT(ret>=min_value && ret<=max_value);
00463 return ret;
00464 }
00465
00466 static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
00467 {
00468 floatmax_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00469
00470 if (ret<min_value || ret>max_value)
00471 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00472 ASSERT(ret>=min_value && ret<=max_value);
00473 return ret;
00474 }
00475
00479 static inline float32_t normal_random(float32_t mean, float32_t std_dev)
00480 {
00481
00482 float32_t ret;
00483 float32_t rand_u;
00484 float32_t rand_v;
00485 float32_t rand_s;
00486 do
00487 {
00488 rand_u = random(-1.0, 1.0);
00489 rand_v = random(-1.0, 1.0);
00490 rand_s = rand_u*rand_u + rand_v*rand_v;
00491 } while ((rand_s == 0) || (rand_s >= 1));
00492
00493
00494 ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
00495 ret = std_dev*ret + mean;
00496 return ret;
00497 }
00498
00502 static inline float64_t normal_random(float64_t mean, float64_t std_dev)
00503 {
00504 float64_t ret;
00505 float64_t rand_u;
00506 float64_t rand_v;
00507 float64_t rand_s;
00508 do
00509 {
00510 rand_u = random(-1.0, 1.0);
00511 rand_v = random(-1.0, 1.0);
00512 rand_s = rand_u*rand_u + rand_v*rand_v;
00513 } while ((rand_s == 0) || (rand_s >= 1));
00514
00515 ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
00516 ret = std_dev*ret + mean;
00517 return ret;
00518 }
00519
00522 static inline float32_t randn_float()
00523 {
00524 return normal_random(0.0, 1.0);
00525 }
00526
00529 static inline float64_t randn_double()
00530 {
00531 return normal_random(0.0, 1.0);
00532 }
00533
00534 template <class T>
00535 static int32_t get_num_nonzero(T* vec, int32_t len)
00536 {
00537 int32_t nnz = 0;
00538 for (index_t i=0; i<len; ++i)
00539 nnz += vec[i] != 0;
00540
00541 return nnz;
00542 }
00543
00548 static inline SGVector<int32_t> randperm_vec(int32_t n)
00549 {
00550 SGVector<int32_t> result(randperm(n), n);
00551 return result;
00552 }
00553
00560 static inline int32_t* randperm(int32_t n)
00561 {
00562 int32_t* perm = SG_MALLOC(int32_t, n);
00563
00564 if (!perm)
00565 return NULL;
00566 for (int32_t i = 0; i < n; i++)
00567 perm[i] = i;
00568 for (int32_t i = 0; i < n; i++)
00569 swap(perm[random(0, n - 1)], perm[i]);
00570 return perm;
00571 }
00572
00573 static inline int64_t nchoosek(int32_t n, int32_t k)
00574 {
00575 int64_t res=1;
00576
00577 for (int32_t i=n-k+1; i<=n; i++)
00578 res*=i;
00579
00580 return res/factorial(k);
00581 }
00582
00586 static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
00587 static void sort(float64_t *a, int32_t*idx, int32_t N);
00588
00589
00590
00591
00592
00593
00594
00596 template <class T>
00597 inline static void radix_sort(T* array, int32_t size)
00598 {
00599 radix_sort_helper(array,size,0);
00600 }
00601
00602 template <class T>
00603 static inline uint8_t byte(T word, uint16_t p)
00604 {
00605 return (word >> (sizeof(T)-p-1) * 8) & 0xff;
00606 }
00607
00608 template <class T>
00609 static void radix_sort_helper(T* array, int32_t size, uint16_t i)
00610 {
00611 static size_t count[256], nc, cmin;
00612 T *ak;
00613 uint8_t c=0;
00614 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
00615 T *an, *aj, *pile[256];
00616 size_t *cp, cmax;
00617
00618
00619 sp = s;
00620 radix_push(array, size, i);
00621
00622
00623 while (sp>s) {
00624 radix_pop(array, size, i);
00625 an = array + size;
00626
00627
00628 if (nc == 0) {
00629 cmin = 0xff;
00630 for (ak = array; ak < an; ak++) {
00631 c = byte(*ak, i);
00632 count[c]++;
00633 if (count[c] == 1) {
00634
00635 if (c < cmin)
00636 cmin = c;
00637 nc++;
00638 }
00639 }
00640
00641
00642 if (sp + nc > s + RADIX_STACK_SIZE) {
00643 radix_sort_helper(array, size, i);
00644 continue;
00645 }
00646 }
00647
00648
00649
00650
00651
00652
00653
00654 olds = bigs = sp;
00655 cmax = 2;
00656 ak = array;
00657 pile[0xff] = an;
00658 for (cp = count + cmin; nc > 0; cp++) {
00659
00660 while (*cp == 0)
00661 cp++;
00662
00663 if (*cp > 1) {
00664
00665 if (*cp > cmax) {
00666 cmax = *cp;
00667 bigs = sp;
00668 }
00669
00670 if (i < sizeof(T)-1)
00671 radix_push(ak, *cp, (uint16_t) (i + 1));
00672 }
00673 pile[cp - count] = ak += *cp;
00674 nc--;
00675 }
00676
00677
00678 swap(*olds, *bigs);
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692 aj = array;
00693 while (aj<an)
00694 {
00695 T r;
00696
00697 for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
00698 swap(*ak, r);
00699
00700 *aj = r;
00701 aj += count[c];
00702 count[c] = 0;
00703 }
00704 }
00705 }
00706
00709 template <class T>
00710 static void insertion_sort(T* output, int32_t size)
00711 {
00712 for (int32_t i=0; i<size-1; i++)
00713 {
00714 int32_t j=i-1;
00715 T value=output[i];
00716 while (j >= 0 && output[j] > value)
00717 {
00718 output[j+1] = output[j];
00719 j--;
00720 }
00721 output[j+1]=value;
00722 }
00723 }
00724
00732 template <class T>
00733 static index_t find_position_to_insert(SGVector<T> vector, T element)
00734 {
00735 index_t i;
00736 for (i=0; i<vector.vlen; ++i)
00737 {
00738 if (vector[i]>element)
00739 break;
00740 }
00741 return i;
00742 }
00743
00749 template <class T>
00750 static void qsort(SGVector<T> vector)
00751 {
00752 qsort(vector.vector, vector.vlen);
00753 }
00754
00757 template <class T>
00758 static void qsort(T* output, int32_t size)
00759 {
00760 if (size==1)
00761 return;
00762
00763 if (size==2)
00764 {
00765 if (output[0] > output [1])
00766 swap(output[0],output[1]);
00767 return;
00768 }
00769
00770 T split=output[size/2];
00771
00772 int32_t left=0;
00773 int32_t right=size-1;
00774
00775 while (left<=right)
00776 {
00777 while (output[left] < split)
00778 left++;
00779 while (output[right] > split)
00780 right--;
00781
00782 if (left<=right)
00783 {
00784 swap(output[left],output[right]);
00785 left++;
00786 right--;
00787 }
00788 }
00789
00790 if (right+1> 1)
00791 qsort(output,right+1);
00792
00793 if (size-left> 1)
00794 qsort(&output[left],size-left);
00795 }
00796
00806 template <class T>
00807 static void qsort(T** vector, index_t length)
00808 {
00809 if (length==1)
00810 return;
00811
00812 if (length==2)
00813 {
00814 if (*vector[0]>*vector[1])
00815 swap(vector[0],vector[1]);
00816 return;
00817 }
00818 T* split=vector[length/2];
00819
00820 int32_t left=0;
00821 int32_t right=length-1;
00822
00823 while (left<=right)
00824 {
00825 while (*vector[left]<*split)
00826 ++left;
00827 while (*vector[right]>*split)
00828 --right;
00829
00830 if (left<=right)
00831 {
00832 swap(vector[left],vector[right]);
00833 ++left;
00834 --right;
00835 }
00836 }
00837
00838 if (right+1>1)
00839 qsort(vector,right+1);
00840
00841 if (length-left>1)
00842 qsort(&vector[left],length-left);
00843 }
00844
00846 template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
00847 {
00848 ASSERT(width>=0);
00849 for (int i=0; i<width; i++)
00850 {
00851 T mask = ((T) 1)<<(sizeof(T)*8-1);
00852 while (mask)
00853 {
00854 if (mask & word)
00855 SG_SPRINT("1");
00856 else
00857 SG_SPRINT("0");
00858
00859 mask>>=1;
00860 }
00861 }
00862 }
00863
00869 template <class T1,class T2>
00870 static void qsort_index(T1* output, T2* index, uint32_t size);
00871
00877 template <class T1,class T2>
00878 static void qsort_backward_index(
00879 T1* output, T2* index, int32_t size);
00880
00888 template <class T1,class T2>
00889 inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
00890 {
00891 int32_t n=0;
00892 thread_qsort<T1,T2> t;
00893 t.output=output;
00894 t.index=index;
00895 t.size=size;
00896 t.qsort_threads=&n;
00897 t.sort_limit=limit;
00898 t.num_threads=n_threads;
00899 parallel_qsort_index<T1,T2>(&t);
00900 }
00901
00902
00903 template <class T1,class T2>
00904 static void* parallel_qsort_index(void* p);
00905
00906
00907
00908
00909 template <class T>
00910 static void min(float64_t* output, T* index, int32_t size);
00911
00912
00913
00914 template <class T>
00915 static void nmin(
00916 float64_t* output, T* index, int32_t size, int32_t n);
00917
00918
00919
00920
00921
00922
00923 template <class T>
00924 static int32_t binary_search_helper(T* output, int32_t size, T elem)
00925 {
00926 int32_t start=0;
00927 int32_t end=size-1;
00928
00929 if (size<1)
00930 return 0;
00931
00932 while (start<end)
00933 {
00934 int32_t middle=(start+end)/2;
00935
00936 if (output[middle]>elem)
00937 end=middle-1;
00938 else if (output[middle]<elem)
00939 start=middle+1;
00940 else
00941 return middle;
00942 }
00943
00944 return start;
00945 }
00946
00947
00948
00949 template <class T>
00950 static inline int32_t binary_search(T* output, int32_t size, T elem)
00951 {
00952 int32_t ind = binary_search_helper(output, size, elem);
00953 if (ind >= 0 && output[ind] == elem)
00954 return ind;
00955 return -1;
00956 }
00957
00958
00959
00960
00961
00962
00963
00964
00965 template<class T>
00966 static inline int32_t binary_search(T** vector, index_t length,
00967 T* elem)
00968 {
00969 int32_t start=0;
00970 int32_t end=length-1;
00971
00972 if (length<1)
00973 return -1;
00974
00975 while (start<end)
00976 {
00977 int32_t middle=(start+end)/2;
00978
00979 if (*vector[middle]>*elem)
00980 end=middle-1;
00981 else if (*vector[middle]<*elem)
00982 start=middle+1;
00983 else
00984 {
00985 start=middle;
00986 break;
00987 }
00988 }
00989
00990 if (start>=0&&*vector[start]==*elem)
00991 return start;
00992
00993 return -1;
00994 }
00995
00996 template <class T>
00997 static int32_t binary_search_max_lower_equal(
00998 T* output, int32_t size, T elem)
00999 {
01000 int32_t ind = binary_search_helper(output, size, elem);
01001
01002 if (output[ind]<=elem)
01003 return ind;
01004 if (ind>0 && output[ind-1] <= elem)
01005 return ind-1;
01006 return -1;
01007 }
01008
01011 static float64_t Align(
01012 char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
01013
01014
01016
01018 inline static uint32_t get_seed()
01019 {
01020 return CMath::seed;
01021 }
01022
01024 inline static uint32_t get_log_range()
01025 {
01026 return CMath::LOGRANGE;
01027 }
01028
01029 #ifdef USE_LOGCACHE
01030
01031 inline static uint32_t get_log_accuracy()
01032 {
01033 return CMath::LOGACCURACY;
01034 }
01035 #endif
01036
01038 inline static int is_finite(double f)
01039 {
01040 #if defined(isfinite) && !defined(SUNOS)
01041 return isfinite(f);
01042 #else
01043 return finite(f);
01044 #endif
01045 }
01046
01048 inline static int is_infinity(double f)
01049 {
01050 #ifdef SUNOS
01051 if (fpclass(f) == FP_NINF || fpclass(f) == FP_PINF)
01052 return 1;
01053 else
01054 return 0;
01055 #else
01056 return isinf(f);
01057 #endif
01058 }
01059
01061 inline static int is_nan(double f)
01062 {
01063 #ifdef SUNOS
01064 return isnand(f);
01065 #else
01066 return isnan(f);
01067 #endif
01068 }
01069
01080 #ifdef USE_LOGCACHE
01081 static inline float64_t logarithmic_sum(float64_t p, float64_t q)
01082 {
01083 float64_t diff;
01084
01085 if (!CMath::is_finite(p))
01086 return q;
01087
01088 if (!CMath::is_finite(q))
01089 {
01090 SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
01091 return NAN;
01092 }
01093 diff = p - q;
01094 if (diff > 0)
01095 return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
01096 return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
01097 }
01098
01100 static void init_log_table();
01101
01103 static int32_t determine_logrange();
01104
01106 static int32_t determine_logaccuracy(int32_t range);
01107 #else
01108 static inline float64_t logarithmic_sum(
01109 float64_t p, float64_t q)
01110 {
01111 float64_t diff;
01112
01113 if (!CMath::is_finite(p))
01114 return q;
01115 if (!CMath::is_finite(q))
01116 return p;
01117 diff = p - q;
01118 if (diff > 0)
01119 return diff > LOGRANGE? p : p + log(1 + exp(-diff));
01120 return -diff > LOGRANGE? q : q + log(1 + exp(diff));
01121 }
01122 #endif
01123 #ifdef USE_LOGSUMARRAY
01124
01129 static inline float64_t logarithmic_sum_array(
01130 float64_t *p, int32_t len)
01131 {
01132 if (len<=2)
01133 {
01134 if (len==2)
01135 return logarithmic_sum(p[0],p[1]) ;
01136 if (len==1)
01137 return p[0];
01138 return -INFTY ;
01139 }
01140 else
01141 {
01142 register float64_t *pp=p ;
01143 if (len%2==1) pp++ ;
01144 for (register int32_t j=0; j < len>>1; j++)
01145 pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
01146 }
01147 return logarithmic_sum_array(p,len%2 + (len>>1)) ;
01148 }
01149 #endif //USE_LOGSUMARRAY
01150
01151
01153 virtual const char* get_name() const { return "Mathematics"; }
01154 public:
01157
01158 static const float64_t INFTY;
01159 static const float64_t ALMOST_INFTY;
01160
01162 static const float64_t ALMOST_NEG_INFTY;
01163
01165 static const float64_t PI;
01166
01168 static const float64_t MACHINE_EPSILON;
01169
01170
01171 static const float64_t MAX_REAL_NUMBER;
01172 static const float64_t MIN_REAL_NUMBER;
01173
01174 protected:
01176 static int32_t LOGRANGE;
01177
01179 static uint32_t seed;
01180 static char* rand_state;
01181
01182 #ifdef USE_LOGCACHE
01183
01185 static int32_t LOGACCURACY;
01187
01188 static float64_t* logtable;
01189 #endif
01190 };
01191
01192
01193 template <class T1,class T2>
01194 void* CMath::parallel_qsort_index(void* p)
01195 {
01196 struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
01197 T1* output=ps->output;
01198 T2* index=ps->index;
01199 uint32_t size=ps->size;
01200 int32_t* qsort_threads=ps->qsort_threads;
01201 int32_t sort_limit=ps->sort_limit;
01202 int32_t num_threads=ps->num_threads;
01203
01204 if (size==2)
01205 {
01206 if (output[0] > output [1])
01207 {
01208 swap(output[0], output[1]);
01209 swap(index[0], index[1]);
01210 }
01211 return NULL;
01212 }
01213
01214 T1 split=output[size/2];
01215
01216 int32_t left=0;
01217 int32_t right=size-1;
01218
01219 while (left<=right)
01220 {
01221 while (output[left] < split)
01222 left++;
01223 while (output[right] > split)
01224 right--;
01225
01226 if (left<=right)
01227 {
01228 swap(output[left], output[right]);
01229 swap(index[left], index[right]);
01230 left++;
01231 right--;
01232 }
01233 }
01234 bool lthread_start=false;
01235 bool rthread_start=false;
01236 pthread_t lthread;
01237 pthread_t rthread;
01238 struct thread_qsort<T1,T2> t1;
01239 struct thread_qsort<T1,T2> t2;
01240
01241 if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
01242 qsort_index(output,index,right+1);
01243 else if (right+1> 1)
01244 {
01245 (*qsort_threads)++;
01246 lthread_start=true;
01247 t1.output=output;
01248 t1.index=index;
01249 t1.size=right+1;
01250 t1.qsort_threads=qsort_threads;
01251 t1.sort_limit=sort_limit;
01252 t1.num_threads=num_threads;
01253 if (pthread_create(<hread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
01254 {
01255 lthread_start=false;
01256 (*qsort_threads)--;
01257 qsort_index(output,index,right+1);
01258 }
01259 }
01260
01261
01262 if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
01263 qsort_index(&output[left],&index[left], size-left);
01264 else if (size-left> 1)
01265 {
01266 (*qsort_threads)++;
01267 rthread_start=true;
01268 t2.output=&output[left];
01269 t2.index=&index[left];
01270 t2.size=size-left;
01271 t2.qsort_threads=qsort_threads;
01272 t2.sort_limit=sort_limit;
01273 t2.num_threads=num_threads;
01274 if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
01275 {
01276 rthread_start=false;
01277 (*qsort_threads)--;
01278 qsort_index(&output[left],&index[left], size-left);
01279 }
01280 }
01281
01282 if (lthread_start)
01283 {
01284 pthread_join(lthread, NULL);
01285 (*qsort_threads)--;
01286 }
01287
01288 if (rthread_start)
01289 {
01290 pthread_join(rthread, NULL);
01291 (*qsort_threads)--;
01292 }
01293
01294 return NULL;
01295 }
01296
01297 template <class T1,class T2>
01298 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
01299 {
01300 if (size==1)
01301 return;
01302
01303 if (size==2)
01304 {
01305 if (output[0] > output [1])
01306 {
01307 swap(output[0],output[1]);
01308 swap(index[0],index[1]);
01309 }
01310 return;
01311 }
01312
01313 T1 split=output[size/2];
01314
01315 int32_t left=0;
01316 int32_t right=size-1;
01317
01318 while (left<=right)
01319 {
01320 while (output[left] < split)
01321 left++;
01322 while (output[right] > split)
01323 right--;
01324
01325 if (left<=right)
01326 {
01327 swap(output[left],output[right]);
01328 swap(index[left],index[right]);
01329 left++;
01330 right--;
01331 }
01332 }
01333
01334 if (right+1> 1)
01335 qsort_index(output,index,right+1);
01336
01337 if (size-left> 1)
01338 qsort_index(&output[left],&index[left], size-left);
01339 }
01340
01341 template <class T1,class T2>
01342 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
01343 {
01344 if (size==2)
01345 {
01346 if (output[0] < output [1])
01347 {
01348 swap(output[0],output[1]);
01349 swap(index[0],index[1]);
01350 }
01351 return;
01352 }
01353
01354
01355 T1 split=output[size/2];
01356
01357 int32_t left=0;
01358 int32_t right=size-1;
01359
01360 while (left<=right)
01361 {
01362 while (output[left] > split)
01363 left++;
01364 while (output[right] < split)
01365 right--;
01366
01367 if (left<=right)
01368 {
01369 swap(output[left],output[right]);
01370 swap(index[left],index[right]);
01371 left++;
01372 right--;
01373 }
01374 }
01375
01376 if (right+1> 1)
01377 qsort_backward_index(output,index,right+1);
01378
01379 if (size-left> 1)
01380 qsort_backward_index(&output[left],&index[left], size-left);
01381 }
01382
01383 template <class T>
01384 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
01385 {
01386 if (6*n*size<13*size*CMath::log(size))
01387 for (int32_t i=0; i<n; i++)
01388 min(&output[i], &index[i], size-i) ;
01389 else
01390 qsort_index(output, index, size) ;
01391 }
01392
01393
01394 template <class T>
01395 void CMath::min(float64_t* output, T* index, int32_t size)
01396 {
01397 if (size<=1)
01398 return;
01399 float64_t min_elem=output[0];
01400 int32_t min_index=0;
01401 for (int32_t i=1; i<size; i++)
01402 {
01403 if (output[i]<min_elem)
01404 {
01405 min_index=i;
01406 min_elem=output[i];
01407 }
01408 }
01409 swap(output[0], output[min_index]);
01410 swap(index[0], index[min_index]);
01411 }
01412 }
01413 #endif