SHOGUN
v2.0.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2012 Fernando José Iglesias García 00008 * Written (W) 2011 Siddharth Kherada 00009 * Written (W) 2011 Justin Patera 00010 * Written (W) 2011 Alesis Novik 00011 * Written (W) 2011-2012 Heiko Strathmann 00012 * Written (W) 1999-2009 Soeren Sonnenburg 00013 * Written (W) 1999-2008 Gunnar Raetsch 00014 * Written (W) 2007 Konrad Rieck 00015 * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society 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 /* Size of RNG seed */ 00088 #define RNG_SEED_SIZE 256 00089 00090 /* Maximum stack size */ 00091 #define RADIX_STACK_SIZE 512 00092 00093 /* Stack macros */ 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 // can't be a>=0?(a):(-a), because compiler complains about 00184 // 'comparison always true' when T is unsigned 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 /* register for fast swaps */ 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 //fall back to double precision sqrt if sqrtl is not 00254 //available 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 // store floating-point bits in integer tmp.i 00276 // and do initial guess for Newton's method 00277 tmp.i = 0x5f3759d5 - (tmp.i >> 1); 00278 x = tmp.f; // convert new bits into float 00279 x = x*(1.5f - xhalf*x*x); // One round of Newton's method 00280 return x; 00281 } 00282 00284 static inline floatmax_t powl(floatmax_t x, floatmax_t n) 00285 { 00286 //fall back to double precision pow if powl is not 00287 //available 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 //seed=42 00404 //SG_SPRINT("initializing random number generator with %d (seed size %d)\n", seed, RNG_SEED_SIZE); 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 // sets up variables & makes sure rand_s.range == (0,1) 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 // the meat & potatos, and then the mean & standard deviation shifting... 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 * Inline function to extract the byte at position p (from left) 00591 * of an 64 bit integer. The function is somewhat identical to 00592 * accessing an array of characters via []. 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 /* Push initial array to stack */ 00619 sp = s; 00620 radix_push(array, size, i); 00621 00622 /* Loop until all digits have been sorted */ 00623 while (sp>s) { 00624 radix_pop(array, size, i); 00625 an = array + size; 00626 00627 /* Make character histogram */ 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 /* Determine smallest character */ 00635 if (c < cmin) 00636 cmin = c; 00637 nc++; 00638 } 00639 } 00640 00641 /* Sort recursively if stack size too small */ 00642 if (sp + nc > s + RADIX_STACK_SIZE) { 00643 radix_sort_helper(array, size, i); 00644 continue; 00645 } 00646 } 00647 00648 /* 00649 * Set pile[]; push incompletely sorted bins onto stack. 00650 * pile[] = pointers to last out-of-place element in bins. 00651 * Before permuting: pile[c-1] + count[c] = pile[c]; 00652 * during deal: pile[c] counts down to pile[c-1]. 00653 */ 00654 olds = bigs = sp; 00655 cmax = 2; 00656 ak = array; 00657 pile[0xff] = an; 00658 for (cp = count + cmin; nc > 0; cp++) { 00659 /* Find next non-empty pile */ 00660 while (*cp == 0) 00661 cp++; 00662 /* Pile with several entries */ 00663 if (*cp > 1) { 00664 /* Determine biggest pile */ 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 /* Play it safe -- biggest bin last. */ 00678 swap(*olds, *bigs); 00679 00680 /* 00681 * Permute misplacements home. Already home: everything 00682 * before aj, and in pile[c], items from pile[c] on. 00683 * Inner loop: 00684 * r = next element to put in place; 00685 * ak = pile[r[i]] = location to put the next element. 00686 * aj = bottom of 1st disordered bin. 00687 * Outer loop: 00688 * Once the 1st disordered bin is done, ie. aj >= ak, 00689 * aj<-aj + count[c] connects the bins in array linked list; 00690 * reset count[c]. 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 //T split=output[random(0,size-1)]; 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 /* finds the smallest element in output and puts that element as the 00908 first element */ 00909 template <class T> 00910 static void min(float64_t* output, T* index, int32_t size); 00911 00912 /* finds the n smallest elements in output and puts these elements as the 00913 first n elements */ 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 /* finds an element in a sorted array via binary search 00922 * returns -1 if not found */ 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 /* finds an element in a sorted array via binary search 00948 * * returns -1 if not found */ 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 /* Finds an element in a sorted array of pointers via binary search 00959 * Every element is dereferenced once before being compared 00960 * 00961 * @param array array of pointers to search in (assumed being sorted) 00962 * @param length length of array 00963 * @param elem pointer to element to search for 00964 * @return index of elem, -1 if not found */ 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 inline 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 /* largest and smallest possible float64_t */ 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 //implementations of template functions 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 //T1 split=output[random(0,size-1)]; 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 //T1 split=output[random(0,size-1)]; 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 //T1 split=output[random(0,size-1)]; 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 /* move the smallest entry in the array to the beginning */ 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