SHOGUN
v1.1.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) 2011 Siddharth Kherada 00008 * Written (W) 2011 Justin Patera 00009 * Written (W) 2011 Alesis Novik 00010 * Written (W) 1999-2009 Soeren Sonnenburg 00011 * Written (W) 1999-2008 Gunnar Raetsch 00012 * Written (W) 2007 Konrad Rieck 00013 * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society 00014 */ 00015 00016 #ifndef __MATHEMATICS_H_ 00017 #define __MATHEMATICS_H_ 00018 00019 #include <shogun/lib/common.h> 00020 #include <shogun/io/SGIO.h> 00021 #include <shogun/mathematics/lapack.h> 00022 #include <shogun/base/SGObject.h> 00023 #include <shogun/base/Parallel.h> 00024 #include <shogun/lib/DataType.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 00046 00048 #ifdef _GLIBCXX_CMATH 00049 #if _GLIBCXX_USE_C99_MATH 00050 #if !_GLIBCXX_USE_C99_FP_MACROS_DYNAMIC 00051 00053 using std::signbit; 00054 00055 using std::fpclassify; 00056 00057 using std::isfinite; 00058 using std::isinf; 00059 using std::isnan; 00060 using std::isnormal; 00061 00062 using std::isgreater; 00063 using std::isgreaterequal; 00064 using std::isless; 00065 using std::islessequal; 00066 using std::islessgreater; 00067 using std::isunordered; 00068 #endif 00069 #endif 00070 #endif 00071 00072 00073 #ifdef _WIN32 00074 #ifndef isnan 00075 #define isnan _isnan 00076 #endif 00077 00078 #ifndef isfinite 00079 #define isfinite _isfinite 00080 #endif 00081 #endif //_WIN32 00082 00083 #ifndef NAN 00084 #include <stdlib.h> 00085 #define NAN (strtod("NAN",NULL)) 00086 #endif 00087 00088 /* Size of RNG seed */ 00089 #define RNG_SEED_SIZE 256 00090 00091 /* Maximum stack size */ 00092 #define RADIX_STACK_SIZE 512 00093 00094 /* Stack macros */ 00095 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i 00096 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si 00097 00098 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00099 00100 template <class T> struct radix_stack_t 00101 { 00103 T *sa; 00105 size_t sn; 00107 uint16_t si; 00108 }; 00109 00111 00113 template <class T1, class T2> struct thread_qsort 00114 { 00116 T1* output; 00118 T2* index; 00120 uint32_t size; 00121 00123 int32_t* qsort_threads; 00125 int32_t sort_limit; 00127 int32_t num_threads; 00128 }; 00129 #endif // DOXYGEN_SHOULD_SKIP_THIS 00130 00131 namespace shogun 00132 { 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 T c=a; 00226 a=b; 00227 b=c; 00228 } 00229 00233 template <class T> 00234 static inline void resize(T* &data, int64_t old_size, int64_t new_size) 00235 { 00236 if (old_size==new_size) 00237 return; 00238 T* new_data = SG_MALLOC(T, new_size); 00239 for (int64_t i=0; i<old_size && i<new_size; i++) 00240 new_data[i]=data[i]; 00241 SG_FREE(data); 00242 data=new_data; 00243 } 00244 00246 template <class T> 00247 static inline T twonorm(T* x, int32_t len) 00248 { 00249 float64_t result=0; 00250 for (int32_t i=0; i<len; i++) 00251 result+=x[i]*x[i]; 00252 00253 return CMath::sqrt(result); 00254 } 00255 00257 template <class T> 00258 static inline T qsq(T* x, int32_t len, float64_t q) 00259 { 00260 float64_t result=0; 00261 for (int32_t i=0; i<len; i++) 00262 result+=CMath::pow(x[i], q); 00263 00264 return result; 00265 } 00266 00268 template <class T> 00269 static inline T qnorm(T* x, int32_t len, float64_t q) 00270 { 00271 ASSERT(q!=0); 00272 return CMath::pow(qsq(x, len, q), 1/q); 00273 } 00274 00276 template <class T> 00277 static inline T sq(T x) 00278 { 00279 return x*x; 00280 } 00281 00283 static inline float32_t sqrt(float32_t x) 00284 { 00285 return ::sqrtf(x); 00286 } 00287 00289 static inline float64_t sqrt(float64_t x) 00290 { 00291 return ::sqrt(x); 00292 } 00293 00295 static inline floatmax_t sqrt(floatmax_t x) 00296 { 00297 //fall back to double precision sqrt if sqrtl is not 00298 //available 00299 #ifdef HAVE_SQRTL 00300 return ::sqrtl(x); 00301 #else 00302 return ::sqrt(x); 00303 #endif 00304 } 00305 00307 static inline float32_t invsqrt(float32_t x) 00308 { 00309 union float_to_int 00310 { 00311 float32_t f; 00312 int32_t i; 00313 }; 00314 00315 float_to_int tmp; 00316 tmp.f=x; 00317 00318 float32_t xhalf = 0.5f * x; 00319 // store floating-point bits in integer tmp.i 00320 // and do initial guess for Newton's method 00321 tmp.i = 0x5f3759d5 - (tmp.i >> 1); 00322 x = tmp.f; // convert new bits into float 00323 x = x*(1.5f - xhalf*x*x); // One round of Newton's method 00324 return x; 00325 } 00326 00328 static inline floatmax_t powl(floatmax_t x, floatmax_t n) 00329 { 00330 //fall back to double precision pow if powl is not 00331 //available 00332 #ifdef HAVE_POWL 00333 return ::powl((long double) x, (long double) n); 00334 #else 00335 return ::pow((double) x, (double) n); 00336 #endif 00337 } 00338 00339 static inline int32_t pow(int32_t x, int32_t n) 00340 { 00341 ASSERT(n>=0); 00342 int32_t result=1; 00343 while (n--) 00344 result*=x; 00345 00346 return result; 00347 } 00348 00349 static inline float64_t pow(float64_t x, int32_t n) 00350 { 00351 if (n>=0) 00352 { 00353 float64_t result=1; 00354 while (n--) 00355 result*=x; 00356 00357 return result; 00358 } 00359 else 00360 return ::pow((double)x, (double)n); 00361 } 00362 00363 static inline float64_t pow(float64_t x, float64_t n) 00364 { 00365 return ::pow((double) x, (double) n); 00366 } 00367 00368 static inline float64_t exp(float64_t x) 00369 { 00370 return ::exp((double) x); 00371 } 00372 00374 static inline float64_t lgamma(float64_t x) 00375 { 00376 return ::lgamma((double) x); 00377 } 00378 00380 static inline float64_t tgamma(float64_t x) 00381 { 00382 return ::tgamma((double) x); 00383 } 00384 00386 static inline float64_t atan(float64_t x) 00387 { 00388 return ::atan((double) x); 00389 } 00390 00391 static inline floatmax_t lgammal(floatmax_t x) 00392 { 00393 #ifdef HAVE_LGAMMAL 00394 return ::lgammal((long double) x); 00395 #else 00396 return ::lgamma((double) x); 00397 #endif // HAVE_LGAMMAL 00398 } 00399 00400 static inline float64_t log10(float64_t v) 00401 { 00402 return ::log(v)/::log(10.0); 00403 } 00404 00405 static inline float64_t log2(float64_t v) 00406 { 00407 #ifdef HAVE_LOG2 00408 return ::log2(v); 00409 #else 00410 return ::log(v)/::log(2.0); 00411 #endif //HAVE_LOG2 00412 } 00413 00414 static inline float64_t log(float64_t v) 00415 { 00416 return ::log(v); 00417 } 00418 00419 static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed) 00420 { 00421 ASSERT(len>0 && xy); 00422 00423 float64_t area = 0.0; 00424 00425 if (!reversed) 00426 { 00427 for (int i=1; i<len; i++) 00428 area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]); 00429 } 00430 else 00431 { 00432 for (int i=1; i<len; i++) 00433 area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]); 00434 } 00435 00436 return area; 00437 } 00438 00439 template <class T> 00440 static void transpose_matrix( 00441 T*& matrix, int32_t& num_feat, int32_t& num_vec) 00442 { 00443 T* transposed=SG_MALLOC(T, num_vec*num_feat); 00444 for (int32_t i=0; i<num_vec; i++) 00445 { 00446 for (int32_t j=0; j<num_feat; j++) 00447 transposed[i+j*num_vec]=matrix[i*num_feat+j]; 00448 } 00449 00450 SG_FREE(matrix); 00451 matrix=transposed; 00452 00453 CMath::swap(num_feat, num_vec); 00454 } 00455 00456 #ifdef HAVE_LAPACK 00457 00458 00459 static float64_t* pinv( 00460 float64_t* matrix, int32_t rows, int32_t cols, 00461 float64_t* target=NULL); 00462 00463 00464 //C := alpha*op( A )*op( B ) + beta*C 00465 //op( X ) = X or op( X ) = X', 00466 static inline void dgemm( 00467 double alpha, const double* A, int rows, int cols, 00468 CBLAS_TRANSPOSE transposeA, double *B, int cols_B, 00469 CBLAS_TRANSPOSE transposeB, double beta, double *C) 00470 { 00471 cblas_dgemm(CblasColMajor, transposeA, transposeB, rows, cols, cols_B, 00472 alpha, A, cols, B, cols_B, beta, C, cols); 00473 } 00474 00475 //y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, 00476 static inline void dgemv( 00477 double alpha, const double *A, int rows, int cols, 00478 const CBLAS_TRANSPOSE transposeA, const double* X, double beta, 00479 double* Y) 00480 { 00481 cblas_dgemv(CblasColMajor, transposeA, 00482 rows, cols, alpha, A, cols, 00483 X, 1, beta, Y, 1); 00484 } 00485 #endif 00486 00487 static inline int64_t factorial(int32_t n) 00488 { 00489 int64_t res=1; 00490 for (int i=2; i<=n; i++) 00491 res*=i ; 00492 return res ; 00493 } 00494 00495 static void init_random(uint32_t initseed=0) 00496 { 00497 if (initseed==0) 00498 { 00499 struct timeval tv; 00500 gettimeofday(&tv, NULL); 00501 seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec); 00502 } 00503 else 00504 seed=initseed; 00505 #if !defined(CYGWIN) && !defined(__INTERIX) 00506 //seed=42 00507 //SG_SPRINT("initializing random number generator with %d (seed size %d)\n", seed, RNG_SEED_SIZE); 00508 initstate(seed, CMath::rand_state, RNG_SEED_SIZE); 00509 #endif 00510 } 00511 00512 static inline int64_t random() 00513 { 00514 #if defined(CYGWIN) || defined(__INTERIX) 00515 return rand(); 00516 #else 00517 return ::random(); 00518 #endif 00519 } 00520 00521 static inline int32_t random(int32_t min_value, int32_t max_value) 00522 { 00523 int32_t ret = min_value + (int32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0))); 00524 ASSERT(ret>=min_value && ret<=max_value); 00525 return ret ; 00526 } 00527 00528 static inline float32_t random(float32_t min_value, float32_t max_value) 00529 { 00530 float32_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX))); 00531 00532 if (ret<min_value || ret>max_value) 00533 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value); 00534 ASSERT(ret>=min_value && ret<=max_value); 00535 return ret; 00536 } 00537 00538 static inline float64_t random(float64_t min_value, float64_t max_value) 00539 { 00540 float64_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX))); 00541 00542 if (ret<min_value || ret>max_value) 00543 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value); 00544 ASSERT(ret>=min_value && ret<=max_value); 00545 return ret; 00546 } 00547 00551 static inline float32_t normal_random(float32_t mean, float32_t std_dev) 00552 { 00553 // sets up variables & makes sure rand_s.range == (0,1) 00554 float32_t ret; 00555 float32_t rand_u; 00556 float32_t rand_v; 00557 float32_t rand_s; 00558 do 00559 { 00560 rand_u = random(-1.0, 1.0); 00561 rand_v = random(-1.0, 1.0); 00562 rand_s = rand_u*rand_u + rand_v*rand_v; 00563 } while ((rand_s == 0) || (rand_s >= 1)); 00564 00565 // the meat & potatos, and then the mean & standard deviation shifting... 00566 ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s); 00567 ret = std_dev*ret + mean; 00568 return ret; 00569 } 00570 00574 static inline float64_t normal_random(float64_t mean, float64_t std_dev) 00575 { 00576 float64_t ret; 00577 float64_t rand_u; 00578 float64_t rand_v; 00579 float64_t rand_s; 00580 do 00581 { 00582 rand_u = random(-1.0, 1.0); 00583 rand_v = random(-1.0, 1.0); 00584 rand_s = rand_u*rand_u + rand_v*rand_v; 00585 } while ((rand_s == 0) || (rand_s >= 1)); 00586 00587 ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s); 00588 ret = std_dev*ret + mean; 00589 return ret; 00590 } 00591 00594 static inline float32_t randn_float() 00595 { 00596 return normal_random(0.0, 1.0); 00597 } 00598 00601 static inline float64_t randn_double() 00602 { 00603 return normal_random(0.0, 1.0); 00604 } 00605 00606 template <class T> 00607 static T* clone_vector(const T* vec, int32_t len) 00608 { 00609 T* result = SG_MALLOC(T, len); 00610 for (int32_t i=0; i<len; i++) 00611 result[i]=vec[i]; 00612 00613 return result; 00614 } 00615 template <class T> 00616 static void fill_vector(T* vec, int32_t len, T value) 00617 { 00618 for (int32_t i=0; i<len; i++) 00619 vec[i]=value; 00620 } 00621 template <class T> 00622 static void range_fill_vector(T* vec, int32_t len, T start=0) 00623 { 00624 for (int32_t i=0; i<len; i++) 00625 vec[i]=i+start; 00626 } 00627 00628 template <class T> 00629 static void random_vector(T* vec, int32_t len, T min_value, T max_value) 00630 { 00631 for (int32_t i=0; i<len; i++) 00632 vec[i]=CMath::random(min_value, max_value); 00633 } 00634 00635 static inline int32_t* randperm(int32_t n) 00636 { 00637 int32_t* perm = SG_MALLOC(int32_t, n); 00638 00639 if (!perm) 00640 return NULL; 00641 for (int32_t i = 0; i < n; i++) 00642 perm[i] = i; 00643 for (int32_t i = 0; i < n; i++) 00644 swap(perm[random(0, n - 1)], perm[i]); 00645 return perm; 00646 } 00647 00648 static inline int64_t nchoosek(int32_t n, int32_t k) 00649 { 00650 int64_t res=1; 00651 00652 for (int32_t i=n-k+1; i<=n; i++) 00653 res*=i; 00654 00655 return res/factorial(k); 00656 } 00657 00659 template <class T> 00660 static inline void vec1_plus_scalar_times_vec2(T* vec1, 00661 T scalar, const T* vec2, int32_t n) 00662 { 00663 for (int32_t i=0; i<n; i++) 00664 vec1[i]+=scalar*vec2[i]; 00665 } 00666 00668 static inline float64_t dot(const bool* v1, const bool* v2, int32_t n) 00669 { 00670 float64_t r=0; 00671 for (int32_t i=0; i<n; i++) 00672 r+=((v1[i]) ? 1 : 0) * ((v2[i]) ? 1 : 0); 00673 return r; 00674 } 00675 00677 static inline floatmax_t dot(const floatmax_t* v1, const floatmax_t* v2, int32_t n) 00678 { 00679 floatmax_t r=0; 00680 for (int32_t i=0; i<n; i++) 00681 r+=v1[i]*v2[i]; 00682 return r; 00683 } 00684 00685 00687 static inline float64_t dot(const float64_t* v1, const float64_t* v2, int32_t n) 00688 { 00689 float64_t r=0; 00690 #ifdef HAVE_LAPACK 00691 int32_t skip=1; 00692 r = cblas_ddot(n, v1, skip, v2, skip); 00693 #else 00694 for (int32_t i=0; i<n; i++) 00695 r+=v1[i]*v2[i]; 00696 #endif 00697 return r; 00698 } 00699 00701 static inline float32_t dot( 00702 const float32_t* v1, const float32_t* v2, int32_t n) 00703 { 00704 float64_t r=0; 00705 #ifdef HAVE_LAPACK 00706 int32_t skip=1; 00707 r = cblas_sdot(n, v1, skip, v2, skip); 00708 #else 00709 for (int32_t i=0; i<n; i++) 00710 r+=v1[i]*v2[i]; 00711 #endif 00712 return r; 00713 } 00714 00716 static inline float64_t dot( 00717 const uint64_t* v1, const uint64_t* v2, int32_t n) 00718 { 00719 float64_t r=0; 00720 for (int32_t i=0; i<n; i++) 00721 r+=((float64_t) v1[i])*v2[i]; 00722 00723 return r; 00724 } 00726 static inline float64_t dot( 00727 const int64_t* v1, const int64_t* v2, int32_t n) 00728 { 00729 float64_t r=0; 00730 for (int32_t i=0; i<n; i++) 00731 r+=((float64_t) v1[i])*v2[i]; 00732 00733 return r; 00734 } 00735 00737 static inline float64_t dot( 00738 const int32_t* v1, const int32_t* v2, int32_t n) 00739 { 00740 float64_t r=0; 00741 for (int32_t i=0; i<n; i++) 00742 r+=((float64_t) v1[i])*v2[i]; 00743 00744 return r; 00745 } 00746 00748 static inline float64_t dot( 00749 const uint32_t* v1, const uint32_t* v2, int32_t n) 00750 { 00751 float64_t r=0; 00752 for (int32_t i=0; i<n; i++) 00753 r+=((float64_t) v1[i])*v2[i]; 00754 00755 return r; 00756 } 00757 00759 static inline float64_t dot( 00760 const uint16_t* v1, const uint16_t* v2, int32_t n) 00761 { 00762 float64_t r=0; 00763 for (int32_t i=0; i<n; i++) 00764 r+=((float64_t) v1[i])*v2[i]; 00765 00766 return r; 00767 } 00768 00770 static inline float64_t dot( 00771 const int16_t* v1, const int16_t* v2, int32_t n) 00772 { 00773 float64_t r=0; 00774 for (int32_t i=0; i<n; i++) 00775 r+=((float64_t) v1[i])*v2[i]; 00776 00777 return r; 00778 } 00779 00781 static inline float64_t dot( 00782 const char* v1, const char* v2, int32_t n) 00783 { 00784 float64_t r=0; 00785 for (int32_t i=0; i<n; i++) 00786 r+=((float64_t) v1[i])*v2[i]; 00787 00788 return r; 00789 } 00790 00792 static inline float64_t dot( 00793 const uint8_t* v1, const uint8_t* v2, int32_t n) 00794 { 00795 float64_t r=0; 00796 for (int32_t i=0; i<n; i++) 00797 r+=((float64_t) v1[i])*v2[i]; 00798 00799 return r; 00800 } 00801 00803 static inline float64_t dot( 00804 const int8_t* v1, const int8_t* v2, int32_t n) 00805 { 00806 float64_t r=0; 00807 for (int32_t i=0; i<n; i++) 00808 r+=((float64_t) v1[i])*v2[i]; 00809 00810 return r; 00811 } 00812 00814 static inline float64_t dot( 00815 const float64_t* v1, const char* v2, int32_t n) 00816 { 00817 float64_t r=0; 00818 for (int32_t i=0; i<n; i++) 00819 r+=((float64_t) v1[i])*v2[i]; 00820 00821 return r; 00822 } 00823 00825 template <class T> 00826 static inline void vector_multiply( 00827 T* target, const T* v1, const T* v2,int32_t len) 00828 { 00829 for (int32_t i=0; i<len; i++) 00830 target[i]=v1[i]*v2[i]; 00831 } 00832 00833 00835 template <class T> 00836 static inline void add( 00837 T* target, T alpha, const T* v1, T beta, const T* v2, 00838 int32_t len) 00839 { 00840 for (int32_t i=0; i<len; i++) 00841 target[i]=alpha*v1[i]+beta*v2[i]; 00842 } 00843 00845 template <class T> 00846 static inline void add_scalar(T alpha, T* vec, int32_t len) 00847 { 00848 for (int32_t i=0; i<len; i++) 00849 vec[i]+=alpha; 00850 } 00851 00853 template <class T> 00854 static inline void scale_vector(T alpha, T* vec, int32_t len) 00855 { 00856 for (int32_t i=0; i<len; i++) 00857 vec[i]*=alpha; 00858 } 00859 00861 template <class T> 00862 static inline T sum(T* vec, int32_t len) 00863 { 00864 T result=0; 00865 for (int32_t i=0; i<len; i++) 00866 result+=vec[i]; 00867 00868 return result; 00869 } 00870 00872 template <class T> 00873 static inline T max(T* vec, int32_t len) 00874 { 00875 ASSERT(len>0); 00876 T maxv=vec[0]; 00877 00878 for (int32_t i=1; i<len; i++) 00879 maxv=CMath::max(vec[i], maxv); 00880 00881 return maxv; 00882 } 00883 00885 template <class T> 00886 static inline T sum_abs(T* vec, int32_t len) 00887 { 00888 T result=0; 00889 for (int32_t i=0; i<len; i++) 00890 result+=CMath::abs(vec[i]); 00891 00892 return result; 00893 } 00894 00896 template <class T> 00897 static inline bool fequal(T x, T y, float64_t precision=1e-6) 00898 { 00899 return CMath::abs(x-y)<precision; 00900 } 00901 00903 static inline float64_t mean(float64_t* vec, int32_t len) 00904 { 00905 SG_SDEPRECATED; 00906 ASSERT(vec); 00907 ASSERT(len>0); 00908 00909 float64_t mean=0; 00910 for (int32_t i=0; i<len; i++) 00911 mean+=vec[i]; 00912 return mean/len; 00913 } 00914 00915 static inline float64_t trace( 00916 float64_t* mat, int32_t cols, int32_t rows) 00917 { 00918 float64_t trace=0; 00919 for (int32_t i=0; i<rows; i++) 00920 trace+=mat[i*cols+i]; 00921 return trace; 00922 } 00923 00927 static void sort(int32_t *a, int32_t cols, int32_t sort_col=0); 00928 static void sort(float64_t *a, int32_t*idx, int32_t N); 00929 00930 /* 00931 * Inline function to extract the byte at position p (from left) 00932 * of an 64 bit integer. The function is somewhat identical to 00933 * accessing an array of characters via []. 00934 */ 00935 00937 template <class T> 00938 inline static void radix_sort(T* array, int32_t size) 00939 { 00940 radix_sort_helper(array,size,0); 00941 } 00942 00943 template <class T> 00944 static inline uint8_t byte(T word, uint16_t p) 00945 { 00946 return (word >> (sizeof(T)-p-1) * 8) & 0xff; 00947 } 00948 00949 template <class T> 00950 static void radix_sort_helper(T* array, int32_t size, uint16_t i) 00951 { 00952 static size_t count[256], nc, cmin; 00953 T *ak; 00954 uint8_t c=0; 00955 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs; 00956 T *an, *aj, *pile[256]; 00957 size_t *cp, cmax; 00958 00959 /* Push initial array to stack */ 00960 sp = s; 00961 radix_push(array, size, i); 00962 00963 /* Loop until all digits have been sorted */ 00964 while (sp>s) { 00965 radix_pop(array, size, i); 00966 an = array + size; 00967 00968 /* Make character histogram */ 00969 if (nc == 0) { 00970 cmin = 0xff; 00971 for (ak = array; ak < an; ak++) { 00972 c = byte(*ak, i); 00973 count[c]++; 00974 if (count[c] == 1) { 00975 /* Determine smallest character */ 00976 if (c < cmin) 00977 cmin = c; 00978 nc++; 00979 } 00980 } 00981 00982 /* Sort recursively if stack size too small */ 00983 if (sp + nc > s + RADIX_STACK_SIZE) { 00984 radix_sort_helper(array, size, i); 00985 continue; 00986 } 00987 } 00988 00989 /* 00990 * Set pile[]; push incompletely sorted bins onto stack. 00991 * pile[] = pointers to last out-of-place element in bins. 00992 * Before permuting: pile[c-1] + count[c] = pile[c]; 00993 * during deal: pile[c] counts down to pile[c-1]. 00994 */ 00995 olds = bigs = sp; 00996 cmax = 2; 00997 ak = array; 00998 pile[0xff] = an; 00999 for (cp = count + cmin; nc > 0; cp++) { 01000 /* Find next non-empty pile */ 01001 while (*cp == 0) 01002 cp++; 01003 /* Pile with several entries */ 01004 if (*cp > 1) { 01005 /* Determine biggest pile */ 01006 if (*cp > cmax) { 01007 cmax = *cp; 01008 bigs = sp; 01009 } 01010 01011 if (i < sizeof(T)-1) 01012 radix_push(ak, *cp, (uint16_t) (i + 1)); 01013 } 01014 pile[cp - count] = ak += *cp; 01015 nc--; 01016 } 01017 01018 /* Play it safe -- biggest bin last. */ 01019 swap(*olds, *bigs); 01020 01021 /* 01022 * Permute misplacements home. Already home: everything 01023 * before aj, and in pile[c], items from pile[c] on. 01024 * Inner loop: 01025 * r = next element to put in place; 01026 * ak = pile[r[i]] = location to put the next element. 01027 * aj = bottom of 1st disordered bin. 01028 * Outer loop: 01029 * Once the 1st disordered bin is done, ie. aj >= ak, 01030 * aj<-aj + count[c] connects the bins in array linked list; 01031 * reset count[c]. 01032 */ 01033 aj = array; 01034 while (aj<an) 01035 { 01036 T r; 01037 01038 for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);) 01039 swap(*ak, r); 01040 01041 *aj = r; 01042 aj += count[c]; 01043 count[c] = 0; 01044 } 01045 } 01046 } 01047 01050 template <class T> 01051 static void insertion_sort(T* output, int32_t size) 01052 { 01053 for (int32_t i=0; i<size-1; i++) 01054 { 01055 int32_t j=i-1; 01056 T value=output[i]; 01057 while (j >= 0 && output[j] > value) 01058 { 01059 output[j+1] = output[j]; 01060 j--; 01061 } 01062 output[j+1]=value; 01063 } 01064 } 01065 01066 01069 template <class T> 01070 static void qsort(T* output, int32_t size) 01071 { 01072 if (size==1) 01073 return; 01074 01075 if (size==2) 01076 { 01077 if (output[0] > output [1]) 01078 swap(output[0],output[1]); 01079 return; 01080 } 01081 //T split=output[random(0,size-1)]; 01082 T split=output[size/2]; 01083 01084 int32_t left=0; 01085 int32_t right=size-1; 01086 01087 while (left<=right) 01088 { 01089 while (output[left] < split) 01090 left++; 01091 while (output[right] > split) 01092 right--; 01093 01094 if (left<=right) 01095 { 01096 swap(output[left],output[right]); 01097 left++; 01098 right--; 01099 } 01100 } 01101 01102 if (right+1> 1) 01103 qsort(output,right+1); 01104 01105 if (size-left> 1) 01106 qsort(&output[left],size-left); 01107 } 01108 01117 template <class T> 01118 static void qsort(SGVector<T*> array) 01119 { 01120 if (array.vlen==1) 01121 return; 01122 01123 if (array.vlen==2) 01124 { 01125 if (*array.vector[0]>*array.vector[1]) 01126 swap(array.vector[0],array.vector[1]); 01127 return; 01128 } 01129 T* split=array.vector[array.vlen/2]; 01130 01131 int32_t left=0; 01132 int32_t right=array.vlen-1; 01133 01134 while (left<=right) 01135 { 01136 while (*array.vector[left]<*split) 01137 ++left; 01138 while (*array.vector[right]>*split) 01139 --right; 01140 01141 if (left<=right) 01142 { 01143 swap(array.vector[left],array.vector[right]); 01144 ++left; 01145 --right; 01146 } 01147 } 01148 01149 if (right+1>1) 01150 qsort(SGVector<T*>(array.vector,right+1)); 01151 01152 if (array.vlen-left>1) 01153 qsort(SGVector<T*>(&array.vector[left],array.vlen-left)); 01154 } 01155 01157 template <class T> static void display_bits(T word, int32_t width=8*sizeof(T)) 01158 { 01159 ASSERT(width>=0); 01160 for (int i=0; i<width; i++) 01161 { 01162 T mask = ((T) 1)<<(sizeof(T)*8-1); 01163 while (mask) 01164 { 01165 if (mask & word) 01166 SG_SPRINT("1"); 01167 else 01168 SG_SPRINT("0"); 01169 01170 mask>>=1; 01171 } 01172 } 01173 } 01174 01176 template <class T> static void display_vector( 01177 const T* vector, int32_t n, const char* name="vector"); 01178 01180 template <class T> static void display_matrix( 01181 const T* matrix, int32_t rows, int32_t cols, const char* name="matrix"); 01182 01188 template <class T1,class T2> 01189 static void qsort_index(T1* output, T2* index, uint32_t size); 01190 01196 template <class T1,class T2> 01197 static void qsort_backward_index( 01198 T1* output, T2* index, int32_t size); 01199 01207 template <class T1,class T2> 01208 inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144) 01209 { 01210 int32_t n=0; 01211 thread_qsort<T1,T2> t; 01212 t.output=output; 01213 t.index=index; 01214 t.size=size; 01215 t.qsort_threads=&n; 01216 t.sort_limit=limit; 01217 t.num_threads=n_threads; 01218 parallel_qsort_index<T1,T2>(&t); 01219 } 01220 01221 01222 template <class T1,class T2> 01223 static void* parallel_qsort_index(void* p); 01224 01225 01226 /* finds the smallest element in output and puts that element as the 01227 first element */ 01228 template <class T> 01229 static void min(float64_t* output, T* index, int32_t size); 01230 01231 /* finds the n smallest elements in output and puts these elements as the 01232 first n elements */ 01233 template <class T> 01234 static void nmin( 01235 float64_t* output, T* index, int32_t size, int32_t n); 01236 01237 /* performs a inplace unique of a vector of type T using quicksort 01238 * returns the new number of elements */ 01239 template <class T> 01240 static int32_t unique(T* output, int32_t size) 01241 { 01242 qsort(output, size); 01243 int32_t i,j=0 ; 01244 for (i=0; i<size; i++) 01245 if (i==0 || output[i]!=output[i-1]) 01246 output[j++]=output[i]; 01247 return j ; 01248 } 01249 01257 static double* compute_eigenvectors(double* matrix, int n, int m) 01258 { 01259 #ifdef HAVE_LAPACK 01260 ASSERT(n == m); 01261 01262 char V='V'; 01263 char U='U'; 01264 int info; 01265 int ord=n; 01266 int lda=n; 01267 double* eigenvalues=SG_CALLOC(float64_t, n+1); 01268 01269 // lapack sym matrix eigenvalues+vectors 01270 wrap_dsyev(V, U, ord, matrix, lda, 01271 eigenvalues, &info); 01272 01273 if (info!=0) 01274 SG_SERROR("DSYEV failed with code %d\n", info); 01275 01276 return eigenvalues; 01277 #else 01278 SG_SERROR("Function not available - Lapack/Atlas not enabled at compile time!\n"); 01279 return NULL; 01280 #endif 01281 } 01282 01283 /* Sums up all rows of a matrix and returns the resulting rowvector */ 01284 template <class T> 01285 static T* get_row_sum(T* matrix, int32_t m, int32_t n) 01286 { 01287 T* rowsums=SG_CALLOC(T, n); 01288 01289 for (int32_t i=0; i<n; i++) 01290 { 01291 for (int32_t j=0; j<m; j++) 01292 rowsums[i]+=matrix[j+int64_t(i)*m]; 01293 } 01294 return rowsums; 01295 } 01296 01297 /* Sums up all columns of a matrix and returns the resulting columnvector */ 01298 template <class T> 01299 static T* get_column_sum(T* matrix, int32_t m, int32_t n) 01300 { 01301 T* colsums=SG_CALLOC(T, m); 01302 01303 for (int32_t i=0; i<n; i++) 01304 { 01305 for (int32_t j=0; j<m; j++) 01306 colsums[j]+=matrix[j+int64_t(i)*m]; 01307 } 01308 return colsums; 01309 } 01310 01311 /* Centers matrix (e.g. kernel matrix in feature space INPLACE */ 01312 template <class T> 01313 static void center_matrix(T* matrix, int32_t m, int32_t n) 01314 { 01315 float64_t num_data=n; 01316 01317 T* colsums=get_column_sum(matrix, m,n); 01318 T* rowsums=get_row_sum(matrix, m,n); 01319 01320 for (int32_t i=0; i<m; i++) 01321 colsums[i]/=num_data; 01322 for (int32_t j=0; j<n; j++) 01323 rowsums[j]/=num_data; 01324 01325 T s=sum(rowsums, n)/num_data; 01326 01327 for (int32_t i=0; i<n; i++) 01328 { 01329 for (int32_t j=0; j<m; j++) 01330 matrix[int64_t(i)*m+j]+=s-colsums[j]-rowsums[i]; 01331 } 01332 01333 SG_FREE(rowsums); 01334 SG_FREE(colsums); 01335 } 01336 01337 /* finds an element in a sorted array via binary search 01338 * returns -1 if not found */ 01339 template <class T> 01340 static int32_t binary_search_helper(T* output, int32_t size, T elem) 01341 { 01342 int32_t start=0; 01343 int32_t end=size-1; 01344 01345 if (size<1) 01346 return 0; 01347 01348 while (start<end) 01349 { 01350 int32_t middle=(start+end)/2; 01351 01352 if (output[middle]>elem) 01353 end=middle-1; 01354 else if (output[middle]<elem) 01355 start=middle+1; 01356 else 01357 return middle; 01358 } 01359 01360 return start; 01361 } 01362 01363 /* finds an element in a sorted array via binary search 01364 * * returns -1 if not found */ 01365 template <class T> 01366 static inline int32_t binary_search(T* output, int32_t size, T elem) 01367 { 01368 int32_t ind = binary_search_helper(output, size, elem); 01369 if (ind >= 0 && output[ind] == elem) 01370 return ind; 01371 return -1; 01372 } 01373 01374 /* Finds an element in a sorted array of pointers via binary search 01375 * Every element is dereferenced once before being compared 01376 * 01377 * @param array array of pointers to search in (assumed being sorted) 01378 * @param elem pointer to element to search for 01379 * @return index of elem, -1 if not found */ 01380 template<class T> 01381 static inline int32_t binary_search(SGVector<T*> array, T* elem) 01382 { 01383 int32_t start=0; 01384 int32_t end=array.vlen-1; 01385 01386 if (array.vlen<1) 01387 return -1; 01388 01389 while (start<end) 01390 { 01391 int32_t middle=(start+end)/2; 01392 01393 if (*array.vector[middle]>*elem) 01394 end=middle-1; 01395 else if (*array.vector[middle]<*elem) 01396 start=middle+1; 01397 else 01398 { 01399 start=middle; 01400 break; 01401 } 01402 } 01403 01404 if (start>=0&&*array.vector[start]==*elem) 01405 return start; 01406 01407 return -1; 01408 } 01409 01410 template <class T> 01411 static int32_t binary_search_max_lower_equal( 01412 T* output, int32_t size, T elem) 01413 { 01414 int32_t ind = binary_search_helper(output, size, elem); 01415 01416 if (output[ind]<=elem) 01417 return ind; 01418 if (ind>0 && output[ind-1] <= elem) 01419 return ind-1; 01420 return -1; 01421 } 01422 01425 static float64_t Align( 01426 char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost); 01427 01428 01430 01433 static float64_t mutual_info(float64_t* p1, float64_t* p2, int32_t len); 01434 01437 static float64_t relative_entropy( 01438 float64_t* p, float64_t* q, int32_t len); 01439 01441 static float64_t entropy(float64_t* p, int32_t len); 01442 01444 inline static uint32_t get_seed() 01445 { 01446 return CMath::seed; 01447 } 01448 01450 inline static uint32_t get_log_range() 01451 { 01452 return CMath::LOGRANGE; 01453 } 01454 01455 #ifdef USE_LOGCACHE 01456 01457 inline static uint32_t get_log_accuracy() 01458 { 01459 return CMath::LOGACCURACY; 01460 } 01461 #endif 01462 01464 inline static int is_finite(double f) 01465 { 01466 #if defined(isfinite) && !defined(SUNOS) 01467 return isfinite(f); 01468 #else 01469 return finite(f); 01470 #endif 01471 } 01472 01474 inline static int is_infinity(double f) 01475 { 01476 #ifdef SUNOS 01477 if (fpclass(f) == FP_NINF || fpclass(f) == FP_PINF) 01478 return 1; 01479 else 01480 return 0; 01481 #else 01482 return isinf(f); 01483 #endif 01484 } 01485 01487 inline static int is_nan(double f) 01488 { 01489 #ifdef SUNOS 01490 return isnand(f); 01491 #else 01492 return isnan(f); 01493 #endif 01494 } 01495 01499 static SGVector<float64_t> fishers_exact_test_for_multiple_2x3_tables(SGMatrix<float64_t> tables); 01500 01504 static float64_t fishers_exact_test_for_2x3_table(SGMatrix<float64_t> table); 01505 01516 #ifdef USE_LOGCACHE 01517 static inline float64_t logarithmic_sum(float64_t p, float64_t q) 01518 { 01519 float64_t diff; 01520 01521 if (!CMath::is_finite(p)) 01522 return q; 01523 01524 if (!CMath::is_finite(q)) 01525 { 01526 SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q); 01527 return NAN; 01528 } 01529 diff = p - q; 01530 if (diff > 0) 01531 return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)]; 01532 return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)]; 01533 } 01534 01536 static void init_log_table(); 01537 01539 static int32_t determine_logrange(); 01540 01542 static int32_t determine_logaccuracy(int32_t range); 01543 #else 01544 static inline float64_t logarithmic_sum( 01545 float64_t p, float64_t q) 01546 { 01547 float64_t diff; 01548 01549 if (!CMath::is_finite(p)) 01550 return q; 01551 if (!CMath::is_finite(q)) 01552 return p; 01553 diff = p - q; 01554 if (diff > 0) 01555 return diff > LOGRANGE? p : p + log(1 + exp(-diff)); 01556 return -diff > LOGRANGE? q : q + log(1 + exp(diff)); 01557 } 01558 #endif 01559 #ifdef USE_LOGSUMARRAY 01560 01565 static inline float64_t logarithmic_sum_array( 01566 float64_t *p, int32_t len) 01567 { 01568 if (len<=2) 01569 { 01570 if (len==2) 01571 return logarithmic_sum(p[0],p[1]) ; 01572 if (len==1) 01573 return p[0]; 01574 return -INFTY ; 01575 } 01576 else 01577 { 01578 register float64_t *pp=p ; 01579 if (len%2==1) pp++ ; 01580 for (register int32_t j=0; j < len>>1; j++) 01581 pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ; 01582 } 01583 return logarithmic_sum_array(p,len%2 + (len>>1)) ; 01584 } 01585 #endif //USE_LOGSUMARRAY 01586 01587 01589 inline virtual const char* get_name() const { return "Mathematics"; } 01590 public: 01593 01594 static const float64_t INFTY; 01595 static const float64_t ALMOST_INFTY; 01596 01598 static const float64_t ALMOST_NEG_INFTY; 01599 01601 static const float64_t PI; 01602 01604 static const float64_t MACHINE_EPSILON; 01605 01606 /* largest and smallest possible float64_t */ 01607 static const float64_t MAX_REAL_NUMBER; 01608 static const float64_t MIN_REAL_NUMBER; 01609 01610 protected: 01612 static int32_t LOGRANGE; 01613 01615 static uint32_t seed; 01616 static char* rand_state; 01617 01618 #ifdef USE_LOGCACHE 01619 01621 static int32_t LOGACCURACY; 01623 01624 static float64_t* logtable; 01625 #endif 01626 }; 01627 01628 //implementations of template functions 01629 template <class T1,class T2> 01630 void* CMath::parallel_qsort_index(void* p) 01631 { 01632 struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p; 01633 T1* output=ps->output; 01634 T2* index=ps->index; 01635 uint32_t size=ps->size; 01636 int32_t* qsort_threads=ps->qsort_threads; 01637 int32_t sort_limit=ps->sort_limit; 01638 int32_t num_threads=ps->num_threads; 01639 01640 if (size==2) 01641 { 01642 if (output[0] > output [1]) 01643 { 01644 swap(output[0], output[1]); 01645 swap(index[0], index[1]); 01646 } 01647 return NULL; 01648 } 01649 //T1 split=output[random(0,size-1)]; 01650 T1 split=output[size/2]; 01651 01652 int32_t left=0; 01653 int32_t right=size-1; 01654 01655 while (left<=right) 01656 { 01657 while (output[left] < split) 01658 left++; 01659 while (output[right] > split) 01660 right--; 01661 01662 if (left<=right) 01663 { 01664 swap(output[left], output[right]); 01665 swap(index[left], index[right]); 01666 left++; 01667 right--; 01668 } 01669 } 01670 bool lthread_start=false; 01671 bool rthread_start=false; 01672 pthread_t lthread; 01673 pthread_t rthread; 01674 struct thread_qsort<T1,T2> t1; 01675 struct thread_qsort<T1,T2> t2; 01676 01677 if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1)) 01678 qsort_index(output,index,right+1); 01679 else if (right+1> 1) 01680 { 01681 (*qsort_threads)++; 01682 lthread_start=true; 01683 t1.output=output; 01684 t1.index=index; 01685 t1.size=right+1; 01686 t1.qsort_threads=qsort_threads; 01687 t1.sort_limit=sort_limit; 01688 t1.num_threads=num_threads; 01689 if (pthread_create(<hread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0) 01690 { 01691 lthread_start=false; 01692 (*qsort_threads)--; 01693 qsort_index(output,index,right+1); 01694 } 01695 } 01696 01697 01698 if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1)) 01699 qsort_index(&output[left],&index[left], size-left); 01700 else if (size-left> 1) 01701 { 01702 (*qsort_threads)++; 01703 rthread_start=true; 01704 t2.output=&output[left]; 01705 t2.index=&index[left]; 01706 t2.size=size-left; 01707 t2.qsort_threads=qsort_threads; 01708 t2.sort_limit=sort_limit; 01709 t2.num_threads=num_threads; 01710 if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0) 01711 { 01712 rthread_start=false; 01713 (*qsort_threads)--; 01714 qsort_index(&output[left],&index[left], size-left); 01715 } 01716 } 01717 01718 if (lthread_start) 01719 { 01720 pthread_join(lthread, NULL); 01721 (*qsort_threads)--; 01722 } 01723 01724 if (rthread_start) 01725 { 01726 pthread_join(rthread, NULL); 01727 (*qsort_threads)--; 01728 } 01729 01730 return NULL; 01731 } 01732 01733 template <class T1,class T2> 01734 void CMath::qsort_index(T1* output, T2* index, uint32_t size) 01735 { 01736 if (size==1) 01737 return; 01738 01739 if (size==2) 01740 { 01741 if (output[0] > output [1]) 01742 { 01743 swap(output[0],output[1]); 01744 swap(index[0],index[1]); 01745 } 01746 return; 01747 } 01748 //T1 split=output[random(0,size-1)]; 01749 T1 split=output[size/2]; 01750 01751 int32_t left=0; 01752 int32_t right=size-1; 01753 01754 while (left<=right) 01755 { 01756 while (output[left] < split) 01757 left++; 01758 while (output[right] > split) 01759 right--; 01760 01761 if (left<=right) 01762 { 01763 swap(output[left],output[right]); 01764 swap(index[left],index[right]); 01765 left++; 01766 right--; 01767 } 01768 } 01769 01770 if (right+1> 1) 01771 qsort_index(output,index,right+1); 01772 01773 if (size-left> 1) 01774 qsort_index(&output[left],&index[left], size-left); 01775 } 01776 01777 template <class T1,class T2> 01778 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size) 01779 { 01780 if (size==2) 01781 { 01782 if (output[0] < output [1]) 01783 { 01784 swap(output[0],output[1]); 01785 swap(index[0],index[1]); 01786 } 01787 return; 01788 } 01789 01790 //T1 split=output[random(0,size-1)]; 01791 T1 split=output[size/2]; 01792 01793 int32_t left=0; 01794 int32_t right=size-1; 01795 01796 while (left<=right) 01797 { 01798 while (output[left] > split) 01799 left++; 01800 while (output[right] < split) 01801 right--; 01802 01803 if (left<=right) 01804 { 01805 swap(output[left],output[right]); 01806 swap(index[left],index[right]); 01807 left++; 01808 right--; 01809 } 01810 } 01811 01812 if (right+1> 1) 01813 qsort_backward_index(output,index,right+1); 01814 01815 if (size-left> 1) 01816 qsort_backward_index(&output[left],&index[left], size-left); 01817 } 01818 01819 template <class T> 01820 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n) 01821 { 01822 if (6*n*size<13*size*CMath::log(size)) 01823 for (int32_t i=0; i<n; i++) 01824 min(&output[i], &index[i], size-i) ; 01825 else 01826 qsort_index(output, index, size) ; 01827 } 01828 01829 /* move the smallest entry in the array to the beginning */ 01830 template <class T> 01831 void CMath::min(float64_t* output, T* index, int32_t size) 01832 { 01833 if (size<=1) 01834 return; 01835 float64_t min_elem=output[0]; 01836 int32_t min_index=0; 01837 for (int32_t i=1; i<size; i++) 01838 { 01839 if (output[i]<min_elem) 01840 { 01841 min_index=i; 01842 min_elem=output[i]; 01843 } 01844 } 01845 swap(output[0], output[min_index]); 01846 swap(index[0], index[min_index]); 01847 } 01848 } 01849 #endif