SHOGUN
v3.2.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) 2013 Thoralf Klein 00008 * Written (W) 2013 Soumyajit De 00009 * Written (W) 2012 Fernando Jose Iglesias Garcia 00010 * Written (W) 2011 Siddharth Kherada 00011 * Written (W) 2011 Justin Patera 00012 * Written (W) 2011 Alesis Novik 00013 * Written (W) 2011-2013 Heiko Strathmann 00014 * Written (W) 1999-2009 Soeren Sonnenburg 00015 * Written (W) 1999-2008 Gunnar Raetsch 00016 * Written (W) 2007 Konrad Rieck 00017 * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society 00018 */ 00019 00020 #ifndef __MATHEMATICS_H_ 00021 #define __MATHEMATICS_H_ 00022 00023 #include <shogun/base/SGObject.h> 00024 #include <shogun/lib/common.h> 00025 #include <shogun/io/SGIO.h> 00026 #include <shogun/base/Parallel.h> 00027 #include <shogun/mathematics/Random.h> 00028 00029 #ifndef _USE_MATH_DEFINES 00030 #define _USE_MATH_DEFINES 00031 #endif 00032 00033 #include <math.h> 00034 #include <stdio.h> 00035 #include <float.h> 00036 #include <sys/types.h> 00037 #ifndef _WIN32 00038 #include <unistd.h> 00039 #endif 00040 00041 #ifdef HAVE_PTHREAD 00042 #include <pthread.h> 00043 #endif 00044 00045 #ifdef SUNOS 00046 #include <ieeefp.h> 00047 #endif 00048 00050 #ifdef log2 00051 #define cygwin_log2 log2 00052 #undef log2 00053 #endif 00054 00055 #ifndef M_PI 00056 #define M_PI 3.14159265358979323846 00057 #endif 00058 00059 #ifdef _WIN32 00060 #ifndef isnan 00061 #define isnan _isnan 00062 #endif 00063 00064 #ifndef isfinite 00065 #define isfinite _isfinite 00066 #endif 00067 #endif //_WIN32 00068 00069 /* Size of RNG seed */ 00070 #define RNG_SEED_SIZE 256 00071 00072 /* Maximum stack size */ 00073 #define RADIX_STACK_SIZE 512 00074 00075 /* Stack macros */ 00076 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i 00077 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si 00078 00079 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00080 00081 template <class T> struct radix_stack_t 00082 { 00084 T *sa; 00086 size_t sn; 00088 uint16_t si; 00089 }; 00090 00092 00094 template <class T1, class T2> struct thread_qsort 00095 { 00097 T1* output; 00099 T2* index; 00101 uint32_t size; 00102 00104 int32_t* qsort_threads; 00106 int32_t sort_limit; 00108 int32_t num_threads; 00109 }; 00110 #endif // DOXYGEN_SHOULD_SKIP_THIS 00111 00112 #define COMPLEX128_ERROR_ONEARG(function) \ 00113 static inline complex128_t function(complex128_t a) \ 00114 { \ 00115 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\ 00116 #function);\ 00117 return complex128_t(0.0, 0.0); \ 00118 } 00119 00120 #define COMPLEX128_STDMATH(function) \ 00121 static inline complex128_t function(complex128_t a) \ 00122 { \ 00123 return std::function(a); \ 00124 } 00125 00126 namespace shogun 00127 { 00129 extern CRandom* sg_rand; 00130 class CSGObject; 00133 class CMath : public CSGObject 00134 { 00135 public: 00139 00140 CMath(); 00141 00143 virtual ~CMath(); 00145 00149 00151 // 00152 template <class T> 00153 static inline T min(T a, T b) 00154 { 00155 return (a<=b) ? a : b; 00156 } 00157 00159 template <class T> 00160 static inline T max(T a, T b) 00161 { 00162 return (a>=b) ? a : b; 00163 } 00164 00166 template <class T> 00167 static inline T clamp(T value, T lb, T ub) 00168 { 00169 if (value<=lb) 00170 return lb; 00171 else if (value>=ub) 00172 return ub; 00173 else 00174 return value; 00175 } 00176 00178 template <class T> 00179 static inline T abs(T a) 00180 { 00181 // can't be a>=0?(a):(-a), because compiler complains about 00182 // 'comparison always true' when T is unsigned 00183 if (a==0) 00184 return 0; 00185 else if (a>0) 00186 return a; 00187 else 00188 return -a; 00189 } 00190 00192 static inline float64_t abs(complex128_t a) 00193 { 00194 float64_t a_real=a.real(); 00195 float64_t a_imag=a.imag(); 00196 return (CMath::sqrt(a_real*a_real+a_imag*a_imag)); 00197 } 00199 00202 00203 static inline float64_t round(float64_t d) 00204 { 00205 return ::floor(d+0.5); 00206 } 00207 00208 static inline float64_t floor(float64_t d) 00209 { 00210 return ::floor(d); 00211 } 00212 00213 static inline float64_t ceil(float64_t d) 00214 { 00215 return ::ceil(d); 00216 } 00217 00219 template <class T> 00220 static inline T sign(T a) 00221 { 00222 if (a==0) 00223 return 0; 00224 else return (a<0) ? (-1) : (+1); 00225 } 00226 00228 template <class T> 00229 static inline void swap(T &a,T &b) 00230 { 00231 /* register for fast swaps */ 00232 register T c=a; 00233 a=b; 00234 b=c; 00235 } 00236 00238 template <class T> 00239 static inline T sq(T x) 00240 { 00241 return x*x; 00242 } 00243 00245 static inline float32_t sqrt(float32_t x) 00246 { 00247 return ::sqrtf(x); 00248 } 00249 00251 static inline float64_t sqrt(float64_t x) 00252 { 00253 return ::sqrt(x); 00254 } 00255 00257 static inline floatmax_t sqrt(floatmax_t x) 00258 { 00259 //fall back to double precision sqrt if sqrtl is not 00260 //available 00261 #ifdef HAVE_SQRTL 00262 return ::sqrtl(x); 00263 #else 00264 return ::sqrt(x); 00265 #endif 00266 } 00267 00269 COMPLEX128_STDMATH(sqrt) 00270 00271 00272 static inline float32_t invsqrt(float32_t x) 00273 { 00274 union float_to_int 00275 { 00276 float32_t f; 00277 int32_t i; 00278 }; 00279 00280 float_to_int tmp; 00281 tmp.f=x; 00282 00283 float32_t xhalf = 0.5f * x; 00284 // store floating-point bits in integer tmp.i 00285 // and do initial guess for Newton's method 00286 tmp.i = 0x5f3759d5 - (tmp.i >> 1); 00287 x = tmp.f; // convert new bits into float 00288 x = x*(1.5f - xhalf*x*x); // One round of Newton's method 00289 return x; 00290 } 00291 00293 static inline floatmax_t powl(floatmax_t x, floatmax_t n) 00294 { 00295 //fall back to double precision pow if powl is not 00296 //available 00297 #ifdef HAVE_POWL 00298 return ::powl((long double) x, (long double) n); 00299 #else 00300 return ::pow((double) x, (double) n); 00301 #endif 00302 } 00303 00304 static inline int32_t pow(bool x, int32_t n) 00305 { 00306 return 0; 00307 } 00308 00309 static inline int32_t pow(int32_t x, int32_t n) 00310 { 00311 ASSERT(n>=0) 00312 int32_t result=1; 00313 while (n--) 00314 result*=x; 00315 00316 return result; 00317 } 00318 00319 static inline float64_t pow(float64_t x, int32_t n) 00320 { 00321 if (n>=0) 00322 { 00323 float64_t result=1; 00324 while (n--) 00325 result*=x; 00326 00327 return result; 00328 } 00329 else 00330 return ::pow((double)x, (double)n); 00331 } 00332 00333 static inline float64_t pow(float64_t x, float64_t n) 00334 { 00335 return ::pow((double) x, (double) n); 00336 } 00337 00339 static inline complex128_t pow(complex128_t x, int32_t n) 00340 { 00341 return std::pow(x, n); 00342 } 00343 00344 static inline complex128_t pow(complex128_t x, complex128_t n) 00345 { 00346 return std::pow(x, n); 00347 } 00348 00349 static inline complex128_t pow(complex128_t x, float64_t n) 00350 { 00351 return std::pow(x, n); 00352 } 00353 00354 static inline complex128_t pow(float64_t x, complex128_t n) 00355 { 00356 return std::pow(x, n); 00357 } 00358 00359 static inline float64_t exp(float64_t x) 00360 { 00361 return ::exp((double) x); 00362 } 00363 00365 COMPLEX128_STDMATH(exp) 00366 00367 00368 static inline float64_t tan(float64_t x) 00369 { 00370 return ::tan((double) x); 00371 } 00372 00374 COMPLEX128_STDMATH(tan) 00375 00376 00377 static inline float64_t atan(float64_t x) 00378 { 00379 return ::atan((double) x); 00380 } 00381 00383 COMPLEX128_ERROR_ONEARG(atan) 00384 00385 00386 static inline float64_t atan2(float64_t x, float64_t y) 00387 { 00388 return ::atan2((double) x, (double) y); 00389 } 00390 00392 COMPLEX128_ERROR_ONEARG(atan2) 00393 00394 00395 static inline float64_t tanh(float64_t x) 00396 { 00397 return ::tanh((double) x); 00398 } 00399 00401 COMPLEX128_STDMATH(tanh) 00402 00403 static inline float64_t log10(float64_t v) 00404 { 00405 return ::log(v)/::log(10.0); 00406 } 00407 00409 COMPLEX128_STDMATH(log10) 00410 00411 static inline float64_t log2(float64_t v) 00412 { 00413 #ifdef HAVE_LOG2 00414 return ::log2(v); 00415 #else 00416 return ::log(v)/::log(2.0); 00417 #endif //HAVE_LOG2 00418 } 00419 00420 static inline float64_t log(float64_t v) 00421 { 00422 return ::log(v); 00423 } 00424 00426 COMPLEX128_STDMATH(log) 00427 00428 static inline index_t floor_log(index_t n) 00429 { 00430 index_t i; 00431 for (i = 0; n != 0; i++) 00432 n >>= 1; 00433 00434 return i; 00435 } 00436 00437 static inline float64_t sin(float64_t x) 00438 { 00439 return ::sin(x); 00440 } 00441 00443 COMPLEX128_STDMATH(sin) 00444 00445 static inline float64_t asin(float64_t x) 00446 { 00447 return ::asin(x); 00448 } 00449 00451 COMPLEX128_ERROR_ONEARG(asin) 00452 00453 static inline float64_t sinh(float64_t x) 00454 { 00455 return ::sinh(x); 00456 } 00457 00459 COMPLEX128_STDMATH(sinh) 00460 00461 static inline float64_t cos(float64_t x) 00462 { 00463 return ::cos(x); 00464 } 00465 00467 COMPLEX128_STDMATH(cos) 00468 00469 static inline float64_t acos(float64_t x) 00470 { 00471 return ::acos(x); 00472 } 00473 00475 COMPLEX128_ERROR_ONEARG(acos) 00476 00477 static inline float64_t cosh(float64_t x) 00478 { 00479 return ::cosh(x); 00480 } 00481 00483 COMPLEX128_STDMATH(cosh) 00484 00485 static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed) 00486 { 00487 ASSERT(len>0 && xy) 00488 00489 float64_t area = 0.0; 00490 00491 if (!reversed) 00492 { 00493 for (int i=1; i<len; i++) 00494 area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]); 00495 } 00496 else 00497 { 00498 for (int i=1; i<len; i++) 00499 area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]); 00500 } 00501 00502 return area; 00503 } 00504 00505 static bool strtof(const char* str, float32_t* float_result); 00506 static bool strtod(const char* str, float64_t* double_result); 00507 static bool strtold(const char* str, floatmax_t* long_double_result); 00508 00509 static inline int64_t factorial(int32_t n) 00510 { 00511 int64_t res=1; 00512 for (int i=2; i<=n; i++) 00513 res*=i ; 00514 return res ; 00515 } 00516 00517 static void init_random(uint32_t initseed=0) 00518 { 00519 if (initseed==0) 00520 seed = CRandom::generate_seed(); 00521 else 00522 seed=initseed; 00523 00524 sg_rand->set_seed(seed); 00525 } 00526 00527 static inline uint64_t random() 00528 { 00529 return sg_rand->random_64(); 00530 } 00531 00532 static inline uint64_t random(uint64_t min_value, uint64_t max_value) 00533 { 00534 return sg_rand->random(min_value, max_value); 00535 } 00536 00537 static inline int64_t random(int64_t min_value, int64_t max_value) 00538 { 00539 return sg_rand->random(min_value, max_value); 00540 } 00541 00542 static inline uint32_t random(uint32_t min_value, uint32_t max_value) 00543 { 00544 return sg_rand->random(min_value, max_value); 00545 } 00546 00547 static inline int32_t random(int32_t min_value, int32_t max_value) 00548 { 00549 return sg_rand->random(min_value, max_value); 00550 } 00551 00552 static inline float32_t random(float32_t min_value, float32_t max_value) 00553 { 00554 return sg_rand->random(min_value, max_value); 00555 } 00556 00557 static inline float64_t random(float64_t min_value, float64_t max_value) 00558 { 00559 return sg_rand->random(min_value, max_value); 00560 } 00561 00562 static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value) 00563 { 00564 return sg_rand->random(min_value, max_value); 00565 } 00566 00570 static inline float32_t normal_random(float32_t mean, float32_t std_dev) 00571 { 00572 // sets up variables & makes sure rand_s.range == (0,1) 00573 float32_t ret; 00574 float32_t rand_u; 00575 float32_t rand_v; 00576 float32_t rand_s; 00577 do 00578 { 00579 rand_u = CMath::random(-1.0, 1.0); 00580 rand_v = CMath::random(-1.0, 1.0); 00581 rand_s = rand_u*rand_u + rand_v*rand_v; 00582 } while ((rand_s == 0) || (rand_s >= 1)); 00583 00584 // the meat & potatos, and then the mean & standard deviation shifting... 00585 ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s); 00586 ret = std_dev*ret + mean; 00587 return ret; 00588 } 00589 00593 static inline float64_t normal_random(float64_t mean, float64_t std_dev) 00594 { 00595 return sg_rand->normal_distrib(mean, std_dev); 00596 } 00597 00600 static inline float32_t randn_float() 00601 { 00602 return normal_random(0.0, 1.0); 00603 } 00604 00607 static inline float64_t randn_double() 00608 { 00609 return sg_rand->std_normal_distrib(); 00610 } 00611 00612 template <class T> 00613 static int32_t get_num_nonzero(T* vec, int32_t len) 00614 { 00615 int32_t nnz = 0; 00616 for (index_t i=0; i<len; ++i) 00617 nnz += vec[i] != 0; 00618 00619 return nnz; 00620 } 00621 00622 static int32_t get_num_nonzero(complex128_t* vec, int32_t len) 00623 { 00624 int32_t nnz=0; 00625 for (index_t i=0; i<len; ++i) 00626 nnz+=vec[i]!=0.0; 00627 return nnz; 00628 } 00629 00630 static inline int64_t nchoosek(int32_t n, int32_t k) 00631 { 00632 int64_t res=1; 00633 00634 for (int32_t i=n-k+1; i<=n; i++) 00635 res*=i; 00636 00637 return res/factorial(k); 00638 } 00639 00647 static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100); 00648 00655 template <class T> 00656 static T log_sum_exp(SGVector<T> values) 00657 { 00658 REQUIRE(values.vector, "Values are empty"); 00659 00660 /* find minimum element index */ 00661 index_t min_index=0; 00662 T X0=values[0]; 00663 if (values.vlen>1) 00664 { 00665 for (index_t i=1; i<values.vlen; ++i) 00666 { 00667 if (values[i]<X0) 00668 { 00669 X0=values[i]; 00670 min_index=i; 00671 } 00672 } 00673 } 00674 00675 /* remove element from vector copy and compute log sum exp */ 00676 SGVector<T> values_without_X0(values.vlen-1); 00677 index_t from_idx=0; 00678 index_t to_idx=0; 00679 for (from_idx=0; from_idx<values.vlen; ++from_idx) 00680 { 00681 if (from_idx!=min_index) 00682 { 00683 values_without_X0[to_idx]=exp(values[from_idx]-X0); 00684 to_idx++; 00685 } 00686 } 00687 00688 return X0+log(SGVector<T>::sum(values_without_X0)+1); 00689 } 00690 00697 template <class T> 00698 static T log_mean_exp(SGVector<T> values) 00699 { 00700 return log_sum_exp(values) - log(values.vlen); 00701 } 00702 00706 static void sort(int32_t *a, int32_t cols, int32_t sort_col=0); 00707 static void sort(float64_t *a, int32_t*idx, int32_t N); 00708 00711 template <class T> 00712 static void qsort(T* output, int32_t size) 00713 { 00714 if (size<=1) 00715 return; 00716 00717 if (size==2) 00718 { 00719 if (output[0] > output [1]) 00720 CMath::swap(output[0],output[1]); 00721 return; 00722 } 00723 //T split=output[random(0,size-1)]; 00724 T split=output[size/2]; 00725 00726 int32_t left=0; 00727 int32_t right=size-1; 00728 00729 while (left<=right) 00730 { 00731 while (output[left] < split) 00732 left++; 00733 while (output[right] > split) 00734 right--; 00735 00736 if (left<=right) 00737 { 00738 CMath::swap(output[left],output[right]); 00739 left++; 00740 right--; 00741 } 00742 } 00743 00744 if (right+1> 1) 00745 qsort(output,right+1); 00746 00747 if (size-left> 1) 00748 qsort(&output[left],size-left); 00749 } 00750 00753 template <class T> 00754 static void insertion_sort(T* output, int32_t size) 00755 { 00756 for (int32_t i=0; i<size-1; i++) 00757 { 00758 int32_t j=i-1; 00759 T value=output[i]; 00760 while (j >= 0 && output[j] > value) 00761 { 00762 output[j+1] = output[j]; 00763 j--; 00764 } 00765 output[j+1]=value; 00766 } 00767 } 00768 00770 template <class T> 00771 inline static void radix_sort(T* array, int32_t size) 00772 { 00773 radix_sort_helper(array,size,0); 00774 } 00775 00776 /* 00777 * Inline function to extract the byte at position p (from left) 00778 * of an 64 bit integer. The function is somewhat identical to 00779 * accessing an array of characters via []. 00780 */ 00781 template <class T> 00782 static inline uint8_t byte(T word, uint16_t p) 00783 { 00784 return (word >> (sizeof(T)-p-1) * 8) & 0xff; 00785 } 00786 00788 static inline uint8_t byte(complex128_t word, uint16_t p) 00789 { 00790 SG_SERROR("CMath::byte():: Not supported for complex128_t\n"); 00791 return uint8_t(0); 00792 } 00793 00794 template <class T> 00795 static void radix_sort_helper(T* array, int32_t size, uint16_t i) 00796 { 00797 static size_t count[256], nc, cmin; 00798 T *ak; 00799 uint8_t c=0; 00800 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs; 00801 T *an, *aj, *pile[256]; 00802 size_t *cp, cmax; 00803 00804 /* Push initial array to stack */ 00805 sp = s; 00806 radix_push(array, size, i); 00807 00808 /* Loop until all digits have been sorted */ 00809 while (sp>s) { 00810 radix_pop(array, size, i); 00811 an = array + size; 00812 00813 /* Make character histogram */ 00814 if (nc == 0) { 00815 cmin = 0xff; 00816 for (ak = array; ak < an; ak++) { 00817 c = byte(*ak, i); 00818 count[c]++; 00819 if (count[c] == 1) { 00820 /* Determine smallest character */ 00821 if (c < cmin) 00822 cmin = c; 00823 nc++; 00824 } 00825 } 00826 00827 /* Sort recursively if stack size too small */ 00828 if (sp + nc > s + RADIX_STACK_SIZE) { 00829 radix_sort_helper(array, size, i); 00830 continue; 00831 } 00832 } 00833 00834 /* 00835 * Set pile[]; push incompletely sorted bins onto stack. 00836 * pile[] = pointers to last out-of-place element in bins. 00837 * Before permuting: pile[c-1] + count[c] = pile[c]; 00838 * during deal: pile[c] counts down to pile[c-1]. 00839 */ 00840 olds = bigs = sp; 00841 cmax = 2; 00842 ak = array; 00843 pile[0xff] = an; 00844 for (cp = count + cmin; nc > 0; cp++) { 00845 /* Find next non-empty pile */ 00846 while (*cp == 0) 00847 cp++; 00848 /* Pile with several entries */ 00849 if (*cp > 1) { 00850 /* Determine biggest pile */ 00851 if (*cp > cmax) { 00852 cmax = *cp; 00853 bigs = sp; 00854 } 00855 00856 if (i < sizeof(T)-1) 00857 radix_push(ak, *cp, (uint16_t) (i + 1)); 00858 } 00859 pile[cp - count] = ak += *cp; 00860 nc--; 00861 } 00862 00863 /* Play it safe -- biggest bin last. */ 00864 swap(*olds, *bigs); 00865 00866 /* 00867 * Permute misplacements home. Already home: everything 00868 * before aj, and in pile[c], items from pile[c] on. 00869 * Inner loop: 00870 * r = next element to put in place; 00871 * ak = pile[r[i]] = location to put the next element. 00872 * aj = bottom of 1st disordered bin. 00873 * Outer loop: 00874 * Once the 1st disordered bin is done, ie. aj >= ak, 00875 * aj<-aj + count[c] connects the bins in array linked list; 00876 * reset count[c]. 00877 */ 00878 aj = array; 00879 while (aj<an) 00880 { 00881 T r; 00882 00883 for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);) 00884 swap(*ak, r); 00885 00886 *aj = r; 00887 aj += count[c]; 00888 count[c] = 0; 00889 } 00890 } 00891 } 00892 00894 static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i) 00895 { 00896 SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n"); 00897 } 00898 00908 template <class T> 00909 static void qsort(T** vector, index_t length) 00910 { 00911 if (length<=1) 00912 return; 00913 00914 if (length==2) 00915 { 00916 if (*vector[0]>*vector[1]) 00917 swap(vector[0],vector[1]); 00918 return; 00919 } 00920 T* split=vector[length/2]; 00921 00922 int32_t left=0; 00923 int32_t right=length-1; 00924 00925 while (left<=right) 00926 { 00927 while (*vector[left]<*split) 00928 ++left; 00929 while (*vector[right]>*split) 00930 --right; 00931 00932 if (left<=right) 00933 { 00934 swap(vector[left],vector[right]); 00935 ++left; 00936 --right; 00937 } 00938 } 00939 00940 if (right+1>1) 00941 qsort(vector,right+1); 00942 00943 if (length-left>1) 00944 qsort(&vector[left],length-left); 00945 } 00946 00948 static void qsort(complex128_t** vector, index_t length) 00949 { 00950 SG_SERROR("CMath::qsort():: Not supported for complex128_t\n"); 00951 } 00952 00954 template <class T> static void display_bits(T word, int32_t width=8*sizeof(T)) 00955 { 00956 ASSERT(width>=0) 00957 for (int i=0; i<width; i++) 00958 { 00959 T mask = ((T) 1)<<(sizeof(T)*8-1); 00960 while (mask) 00961 { 00962 if (mask & word) 00963 SG_SPRINT("1") 00964 else 00965 SG_SPRINT("0") 00966 00967 mask>>=1; 00968 } 00969 } 00970 } 00971 00973 static void display_bits(complex128_t word, 00974 int32_t width=8*sizeof(complex128_t)) 00975 { 00976 SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n"); 00977 } 00978 00984 template <class T1,class T2> 00985 static void qsort_index(T1* output, T2* index, uint32_t size); 00986 00988 template <class T> 00989 static void qsort_index(complex128_t* output, T* index, uint32_t size) 00990 { 00991 SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n"); 00992 } 00993 00999 template <class T1,class T2> 01000 static void qsort_backward_index( 01001 T1* output, T2* index, int32_t size); 01002 01004 template <class T> 01005 static void qsort_backword_index( 01006 complex128_t* output, T* index, uint32_t size) 01007 { 01008 SG_SERROR("CMath::qsort_backword_index():: \ 01009 Not supported for complex128_t\n"); 01010 } 01011 01019 template <class T1,class T2> 01020 inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144) 01021 { 01022 int32_t n=0; 01023 thread_qsort<T1,T2> t; 01024 t.output=output; 01025 t.index=index; 01026 t.size=size; 01027 t.qsort_threads=&n; 01028 t.sort_limit=limit; 01029 t.num_threads=n_threads; 01030 parallel_qsort_index<T1,T2>(&t); 01031 } 01032 01034 template <class T> 01035 inline static void parallel_qsort_index(complex128_t* output, T* index, 01036 uint32_t size, int32_t n_threads, int32_t limit=0) 01037 { 01038 SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n"); 01039 } 01040 01041 01042 template <class T1,class T2> 01043 static void* parallel_qsort_index(void* p); 01044 01045 01046 /* finds the smallest element in output and puts that element as the 01047 first element */ 01048 template <class T> 01049 static void min(float64_t* output, T* index, int32_t size); 01050 01052 static void min(float64_t* output, complex128_t* index, int32_t size) 01053 { 01054 SG_SERROR("CMath::min():: Not supported for complex128_t\n"); 01055 } 01056 01057 /* finds the n smallest elements in output and puts these elements as the 01058 first n elements */ 01059 template <class T> 01060 static void nmin( 01061 float64_t* output, T* index, int32_t size, int32_t n); 01062 01064 static void nmin(float64_t* output, complex128_t* index, 01065 int32_t size, int32_t n) 01066 { 01067 SG_SERROR("CMath::nmin():: Not supported for complex128_t\n"); 01068 } 01069 01070 01071 01072 /* finds an element in a sorted array via binary search 01073 * returns -1 if not found */ 01074 template <class T> 01075 static int32_t binary_search_helper(T* output, int32_t size, T elem) 01076 { 01077 int32_t start=0; 01078 int32_t end=size-1; 01079 01080 if (size<1) 01081 return 0; 01082 01083 while (start<end) 01084 { 01085 int32_t middle=(start+end)/2; 01086 01087 if (output[middle]>elem) 01088 end=middle-1; 01089 else if (output[middle]<elem) 01090 start=middle+1; 01091 else 01092 return middle; 01093 } 01094 01095 return start; 01096 } 01097 01099 static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem) 01100 { 01101 SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n"); 01102 return int32_t(0); 01103 } 01104 01105 /* finds an element in a sorted array via binary search 01106 * * returns -1 if not found */ 01107 template <class T> 01108 static inline int32_t binary_search(T* output, int32_t size, T elem) 01109 { 01110 int32_t ind = binary_search_helper(output, size, elem); 01111 if (ind >= 0 && output[ind] == elem) 01112 return ind; 01113 return -1; 01114 } 01115 01117 static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem) 01118 { 01119 SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n"); 01120 return int32_t(-1); 01121 } 01122 01123 01124 /* Finds an element in a sorted array of pointers via binary search 01125 * Every element is dereferenced once before being compared 01126 * 01127 * @param array array of pointers to search in (assumed being sorted) 01128 * @param length length of array 01129 * @param elem pointer to element to search for 01130 * @return index of elem, -1 if not found */ 01131 template<class T> 01132 static inline int32_t binary_search(T** vector, index_t length, 01133 T* elem) 01134 { 01135 int32_t start=0; 01136 int32_t end=length-1; 01137 01138 if (length<1) 01139 return -1; 01140 01141 while (start<end) 01142 { 01143 int32_t middle=(start+end)/2; 01144 01145 if (*vector[middle]>*elem) 01146 end=middle-1; 01147 else if (*vector[middle]<*elem) 01148 start=middle+1; 01149 else 01150 { 01151 start=middle; 01152 break; 01153 } 01154 } 01155 01156 if (start>=0&&*vector[start]==*elem) 01157 return start; 01158 01159 return -1; 01160 } 01161 01163 static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem) 01164 { 01165 SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n"); 01166 return int32_t(-1); 01167 } 01168 01169 template <class T> 01170 static int32_t binary_search_max_lower_equal( 01171 T* output, int32_t size, T elem) 01172 { 01173 int32_t ind = binary_search_helper(output, size, elem); 01174 01175 if (output[ind]<=elem) 01176 return ind; 01177 if (ind>0 && output[ind-1] <= elem) 01178 return ind-1; 01179 return -1; 01180 } 01181 01183 static int32_t binary_search_max_lower_equal(complex128_t* output, 01184 int32_t size, complex128_t elem) 01185 { 01186 SG_SERROR("CMath::binary_search_max_lower_equal():: \ 01187 Not supported for complex128_t\n"); 01188 return int32_t(-1); 01189 } 01190 01193 static float64_t Align( 01194 char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost); 01195 01196 01198 01200 static float64_t real(complex128_t c) 01201 { 01202 return c.real(); 01203 } 01204 01206 static float64_t imag(complex128_t c) 01207 { 01208 return c.imag(); 01209 } 01210 01212 inline static uint32_t get_seed() 01213 { 01214 return CMath::seed; 01215 } 01216 01218 inline static uint32_t get_log_range() 01219 { 01220 return CMath::LOGRANGE; 01221 } 01222 01223 #ifdef USE_LOGCACHE 01224 01225 inline static uint32_t get_log_accuracy() 01226 { 01227 return CMath::LOGACCURACY; 01228 } 01229 #endif 01230 01232 static int is_finite(double f); 01233 01235 static int is_infinity(double f); 01236 01238 static int is_nan(double f); 01239 01250 #ifdef USE_LOGCACHE 01251 static inline float64_t logarithmic_sum(float64_t p, float64_t q) 01252 { 01253 float64_t diff; 01254 01255 if (!CMath::is_finite(p)) 01256 return q; 01257 01258 if (!CMath::is_finite(q)) 01259 { 01260 SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q) 01261 return NOT_A_NUMBER; 01262 } 01263 diff = p - q; 01264 if (diff > 0) 01265 return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)]; 01266 return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)]; 01267 } 01268 01270 static void init_log_table(); 01271 01273 static int32_t determine_logrange(); 01274 01276 static int32_t determine_logaccuracy(int32_t range); 01277 #else 01278 static inline float64_t logarithmic_sum( 01279 float64_t p, float64_t q) 01280 { 01281 float64_t diff; 01282 01283 if (!CMath::is_finite(p)) 01284 return q; 01285 if (!CMath::is_finite(q)) 01286 return p; 01287 diff = p - q; 01288 if (diff > 0) 01289 return diff > LOGRANGE? p : p + log(1 + exp(-diff)); 01290 return -diff > LOGRANGE? q : q + log(1 + exp(diff)); 01291 } 01292 #endif 01293 #ifdef USE_LOGSUMARRAY 01294 01299 static inline float64_t logarithmic_sum_array( 01300 float64_t *p, int32_t len) 01301 { 01302 if (len<=2) 01303 { 01304 if (len==2) 01305 return logarithmic_sum(p[0],p[1]) ; 01306 if (len==1) 01307 return p[0]; 01308 return -INFTY ; 01309 } 01310 else 01311 { 01312 register float64_t *pp=p ; 01313 if (len%2==1) pp++ ; 01314 for (register int32_t j=0; j < len>>1; j++) 01315 pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ; 01316 } 01317 return logarithmic_sum_array(p,len%2 + (len>>1)) ; 01318 } 01319 #endif //USE_LOGSUMARRAY 01320 01321 01323 virtual const char* get_name() const { return "Math"; } 01324 public: 01327 01328 static const float64_t NOT_A_NUMBER; 01330 static const float64_t INFTY; 01331 static const float64_t ALMOST_INFTY; 01332 01334 static const float64_t ALMOST_NEG_INFTY; 01335 01337 static const float64_t PI; 01338 01340 static const float64_t MACHINE_EPSILON; 01341 01342 /* largest and smallest possible float64_t */ 01343 static const float64_t MAX_REAL_NUMBER; 01344 static const float64_t MIN_REAL_NUMBER; 01345 01346 protected: 01348 static int32_t LOGRANGE; 01349 01351 static uint32_t seed; 01352 01353 #ifdef USE_LOGCACHE 01354 01356 static int32_t LOGACCURACY; 01358 01359 static float64_t* logtable; 01360 #endif 01361 }; 01362 01363 //implementations of template functions 01364 template <class T1,class T2> 01365 void* CMath::parallel_qsort_index(void* p) 01366 { 01367 struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p; 01368 T1* output=ps->output; 01369 T2* index=ps->index; 01370 int32_t size=ps->size; 01371 int32_t* qsort_threads=ps->qsort_threads; 01372 int32_t sort_limit=ps->sort_limit; 01373 int32_t num_threads=ps->num_threads; 01374 01375 if (size<2) 01376 { 01377 return NULL; 01378 } 01379 01380 if (size==2) 01381 { 01382 if (output[0] > output [1]) 01383 { 01384 swap(output[0], output[1]); 01385 swap(index[0], index[1]); 01386 } 01387 return NULL; 01388 } 01389 //T1 split=output[random(0,size-1)]; 01390 T1 split=output[size/2]; 01391 01392 int32_t left=0; 01393 int32_t right=size-1; 01394 01395 while (left<=right) 01396 { 01397 while (output[left] < split) 01398 left++; 01399 while (output[right] > split) 01400 right--; 01401 01402 if (left<=right) 01403 { 01404 swap(output[left], output[right]); 01405 swap(index[left], index[right]); 01406 left++; 01407 right--; 01408 } 01409 } 01410 bool lthread_start=false; 01411 bool rthread_start=false; 01412 pthread_t lthread; 01413 pthread_t rthread; 01414 struct thread_qsort<T1,T2> t1; 01415 struct thread_qsort<T1,T2> t2; 01416 01417 if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1)) 01418 qsort_index(output,index,right+1); 01419 else if (right+1> 1) 01420 { 01421 (*qsort_threads)++; 01422 lthread_start=true; 01423 t1.output=output; 01424 t1.index=index; 01425 t1.size=right+1; 01426 t1.qsort_threads=qsort_threads; 01427 t1.sort_limit=sort_limit; 01428 t1.num_threads=num_threads; 01429 if (pthread_create(<hread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0) 01430 { 01431 lthread_start=false; 01432 (*qsort_threads)--; 01433 qsort_index(output,index,right+1); 01434 } 01435 } 01436 01437 01438 if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1)) 01439 qsort_index(&output[left],&index[left], size-left); 01440 else if (size-left> 1) 01441 { 01442 (*qsort_threads)++; 01443 rthread_start=true; 01444 t2.output=&output[left]; 01445 t2.index=&index[left]; 01446 t2.size=size-left; 01447 t2.qsort_threads=qsort_threads; 01448 t2.sort_limit=sort_limit; 01449 t2.num_threads=num_threads; 01450 if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0) 01451 { 01452 rthread_start=false; 01453 (*qsort_threads)--; 01454 qsort_index(&output[left],&index[left], size-left); 01455 } 01456 } 01457 01458 if (lthread_start) 01459 { 01460 pthread_join(lthread, NULL); 01461 (*qsort_threads)--; 01462 } 01463 01464 if (rthread_start) 01465 { 01466 pthread_join(rthread, NULL); 01467 (*qsort_threads)--; 01468 } 01469 01470 return NULL; 01471 } 01472 01473 template <class T1,class T2> 01474 void CMath::qsort_index(T1* output, T2* index, uint32_t size) 01475 { 01476 if (size<=1) 01477 return; 01478 01479 if (size==2) 01480 { 01481 if (output[0] > output [1]) 01482 { 01483 swap(output[0],output[1]); 01484 swap(index[0],index[1]); 01485 } 01486 return; 01487 } 01488 //T1 split=output[random(0,size-1)]; 01489 T1 split=output[size/2]; 01490 01491 int32_t left=0; 01492 int32_t right=size-1; 01493 01494 while (left<=right) 01495 { 01496 while (output[left] < split) 01497 left++; 01498 while (output[right] > split) 01499 right--; 01500 01501 if (left<=right) 01502 { 01503 swap(output[left],output[right]); 01504 swap(index[left],index[right]); 01505 left++; 01506 right--; 01507 } 01508 } 01509 01510 if (right+1> 1) 01511 qsort_index(output,index,right+1); 01512 01513 if (size-left> 1) 01514 qsort_index(&output[left],&index[left], size-left); 01515 } 01516 01517 template <class T1,class T2> 01518 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size) 01519 { 01520 if (size<=1) 01521 return; 01522 01523 if (size==2) 01524 { 01525 if (output[0] < output [1]) 01526 { 01527 swap(output[0],output[1]); 01528 swap(index[0],index[1]); 01529 } 01530 return; 01531 } 01532 01533 //T1 split=output[random(0,size-1)]; 01534 T1 split=output[size/2]; 01535 01536 int32_t left=0; 01537 int32_t right=size-1; 01538 01539 while (left<=right) 01540 { 01541 while (output[left] > split) 01542 left++; 01543 while (output[right] < split) 01544 right--; 01545 01546 if (left<=right) 01547 { 01548 swap(output[left],output[right]); 01549 swap(index[left],index[right]); 01550 left++; 01551 right--; 01552 } 01553 } 01554 01555 if (right+1> 1) 01556 qsort_backward_index(output,index,right+1); 01557 01558 if (size-left> 1) 01559 qsort_backward_index(&output[left],&index[left], size-left); 01560 } 01561 01562 template <class T> 01563 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n) 01564 { 01565 if (6*n*size<13*size*CMath::log(size)) 01566 for (int32_t i=0; i<n; i++) 01567 min(&output[i], &index[i], size-i) ; 01568 else 01569 qsort_index(output, index, size) ; 01570 } 01571 01572 /* move the smallest entry in the array to the beginning */ 01573 template <class T> 01574 void CMath::min(float64_t* output, T* index, int32_t size) 01575 { 01576 if (size<=1) 01577 return; 01578 float64_t min_elem=output[0]; 01579 int32_t min_index=0; 01580 for (int32_t i=1; i<size; i++) 01581 { 01582 if (output[i]<min_elem) 01583 { 01584 min_index=i; 01585 min_elem=output[i]; 01586 } 01587 } 01588 swap(output[0], output[min_index]); 01589 swap(index[0], index[min_index]); 01590 } 01591 01592 #define COMPLEX128_ERROR_ONEARG_T(function) \ 01593 template <> \ 01594 inline complex128_t CMath::function<complex128_t>(complex128_t a) \ 01595 { \ 01596 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\ 01597 #function);\ 01598 return complex128_t(0.0, 0.0); \ 01599 } 01600 01601 #define COMPLEX128_ERROR_TWOARGS_T(function) \ 01602 template <> \ 01603 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \ 01604 { \ 01605 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\ 01606 #function);\ 01607 return complex128_t(0.0, 0.0); \ 01608 } 01609 01610 #define COMPLEX128_ERROR_THREEARGS_T(function) \ 01611 template <> \ 01612 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \ 01613 { \ 01614 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\ 01615 #function);\ 01616 return complex128_t(0.0, 0.0); \ 01617 } 01618 01619 #define COMPLEX128_ERROR_SORT_T(function) \ 01620 template <> \ 01621 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \ 01622 { \ 01623 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\ 01624 #function);\ 01625 } 01626 01628 COMPLEX128_ERROR_TWOARGS_T(min) 01629 01630 01631 COMPLEX128_ERROR_TWOARGS_T(max) 01632 01634 COMPLEX128_ERROR_THREEARGS_T(clamp) 01635 01637 // COMPLEX128_ERROR_ONEARG_T(sign) 01638 01640 COMPLEX128_ERROR_SORT_T(qsort) 01641 01643 COMPLEX128_ERROR_SORT_T(insertion_sort) 01644 01646 COMPLEX128_ERROR_SORT_T(radix_sort) 01647 01648 } 01649 #undef COMPLEX128_ERROR_ONEARG 01650 #undef COMPLEX128_ERROR_ONEARG_T 01651 #undef COMPLEX128_ERROR_TWOARGS_T 01652 #undef COMPLEX128_ERROR_THREEARGS_T 01653 #undef COMPLEX128_STDMATH 01654 #undef COMPLEX128_ERROR_SORT_T 01655 #endif