00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MATHEMATICS_H_
00014 #define __MATHEMATICS_H_
00015
00016 #include "lib/common.h"
00017 #include "lib/io.h"
00018 #include "lib/lapack.h"
00019 #include "base/SGObject.h"
00020 #include "base/Parallel.h"
00021
00022 #include <math.h>
00023 #include <stdio.h>
00024 #include <float.h>
00025 #include <pthread.h>
00026 #include <unistd.h>
00027 #include <sys/types.h>
00028 #include <sys/time.h>
00029 #include <time.h>
00030
00031 #ifdef SUNOS
00032 #include <ieeefp.h>
00033 #endif
00034
00036 #ifdef _GLIBCXX_CMATH
00037 #if _GLIBCXX_USE_C99_MATH
00038 #if !_GLIBCXX_USE_C99_FP_MACROS_DYNAMIC
00039
00041 using std::signbit;
00042
00043 using std::fpclassify;
00044
00045 using std::isfinite;
00046 using std::isinf;
00047 using std::isnan;
00048 using std::isnormal;
00049
00050 using std::isgreater;
00051 using std::isgreaterequal;
00052 using std::isless;
00053 using std::islessequal;
00054 using std::islessgreater;
00055 using std::isunordered;
00056 #endif
00057 #endif
00058 #endif
00060
00061 #ifdef _WIN32
00062 #ifndef isnan
00063 #define isnan _isnan
00064 #endif
00065
00066 #ifndef isfinite
00067 #define isfinite _isfinite
00068 #endif
00069 #endif //_WIN32
00070
00071 #ifndef NAN
00072 #include <stdlib.h>
00073 #define NAN (strtod("NAN",NULL))
00074 #endif
00075
00076
00077 #define RNG_SEED_SIZE 256
00078
00079
00080 #define RADIX_STACK_SIZE 512
00081
00082
00083 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
00084 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
00085
00086 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00087
00088 template <class T> struct radix_stack_t
00089 {
00091 T *sa;
00093 size_t sn;
00095 uint16_t si;
00096 };
00097
00099
00101 template <class T1, class T2> struct thread_qsort
00102 {
00104 T1* output;
00106 T2* index;
00108 uint32_t size;
00109
00111 int32_t* qsort_threads;
00113 int32_t sort_limit;
00115 int32_t num_threads;
00116 };
00117 #endif // DOXYGEN_SHOULD_SKIP_THIS
00118
00119 namespace shogun
00120 {
00123 class CMath : public CSGObject
00124 {
00125 public:
00129
00130 CMath();
00131
00133 virtual ~CMath();
00135
00139
00141
00142 template <class T>
00143 static inline T min(T a, T b)
00144 {
00145 return (a<=b) ? a : b;
00146 }
00147
00149 template <class T>
00150 static inline T max(T a, T b)
00151 {
00152 return (a>=b) ? a : b;
00153 }
00154
00156 template <class T>
00157 static inline T clamp(T value, T lb, T ub)
00158 {
00159 if (value<=lb)
00160 return lb;
00161 else if (value>=ub)
00162 return ub;
00163 else
00164 return value;
00165 }
00166
00168 template <class T>
00169 static inline T abs(T a)
00170 {
00171
00172
00173 if (a==0)
00174 return 0;
00175 else if (a>0)
00176 return a;
00177 else
00178 return -a;
00179 }
00181
00184
00185 static inline float64_t round(float64_t d)
00186 {
00187 return ::floor(d+0.5);
00188 }
00189
00190 static inline float64_t floor(float64_t d)
00191 {
00192 return ::floor(d);
00193 }
00194
00195 static inline float64_t ceil(float64_t d)
00196 {
00197 return ::ceil(d);
00198 }
00199
00201 template <class T>
00202 static inline T sign(T a)
00203 {
00204 if (a==0)
00205 return 0;
00206 else return (a<0) ? (-1) : (+1);
00207 }
00208
00210 template <class T>
00211 static inline void swap(T &a,T &b)
00212 {
00213 T c=a;
00214 a=b;
00215 b=c;
00216 }
00217
00221 template <class T>
00222 static inline void resize(T* &data, int64_t old_size, int64_t new_size)
00223 {
00224 if (old_size==new_size)
00225 return;
00226 T* new_data = new T[new_size];
00227 for (int64_t i=0; i<old_size && i<new_size; i++)
00228 new_data[i]=data[i];
00229 delete[] data;
00230 data=new_data;
00231 }
00232
00234 template <class T>
00235 static inline T twonorm(T* x, int32_t len)
00236 {
00237 float64_t result=0;
00238 for (int32_t i=0; i<len; i++)
00239 result+=x[i]*x[i];
00240
00241 return CMath::sqrt(result);
00242 }
00243
00245 template <class T>
00246 static inline T qsq(T* x, int32_t len, float64_t q)
00247 {
00248 float64_t result=0;
00249 for (int32_t i=0; i<len; i++)
00250 result+=CMath::pow(x[i], q);
00251
00252 return result;
00253 }
00254
00256 template <class T>
00257 static inline T qnorm(T* x, int32_t len, float64_t q)
00258 {
00259 ASSERT(q!=0);
00260 return CMath::pow(qsq(x, len, q), 1/q);
00261 }
00262
00264 template <class T>
00265 static inline T sq(T x)
00266 {
00267 return x*x;
00268 }
00269
00271 static inline float32_t sqrt(float32_t x)
00272 {
00273 return ::sqrtf(x);
00274 }
00275
00277 static inline float64_t sqrt(float64_t x)
00278 {
00279 return ::sqrt(x);
00280 }
00281
00283 static inline floatmax_t sqrt(floatmax_t x)
00284 {
00285
00286
00287 #ifdef HAVE_SQRTL
00288 return ::sqrtl(x);
00289 #else
00290 return ::sqrt(x);
00291 #endif
00292 }
00293
00294
00296 static inline floatmax_t powl(floatmax_t x, floatmax_t n)
00297 {
00298
00299
00300 #ifdef HAVE_POWL
00301 return ::powl((long double) x, (long double) n);
00302 #else
00303 return ::pow((double) x, (double) n);
00304 #endif
00305 }
00306
00307 static inline int32_t pow(int32_t x, int32_t n)
00308 {
00309 ASSERT(n>=0);
00310 int32_t result=1;
00311 while (n--)
00312 result*=x;
00313
00314 return result;
00315 }
00316
00317 static inline float64_t pow(float64_t x, int32_t n)
00318 {
00319 ASSERT(n>=0);
00320 float64_t result=1;
00321 while (n--)
00322 result*=x;
00323
00324 return result;
00325 }
00326
00327 static inline float64_t pow(float64_t x, float64_t n)
00328 {
00329 return ::pow((double) x, (double) n);
00330 }
00331
00332 static inline float64_t exp(float64_t x)
00333 {
00334 return ::exp((double) x);
00335 }
00336
00337 static inline float64_t log10(float64_t v)
00338 {
00339 return ::log(v)/::log(10.0);
00340 }
00341
00342 #ifdef HAVE_LOG2
00343 static inline float64_t log2(float64_t v)
00344 {
00345 return ::log2(v);
00346 }
00347 #else
00348 static inline float64_t log2(float64_t v)
00349 {
00350 return ::log(v)/::log(2.0);
00351 }
00352 #endif //HAVE_LOG2
00353
00354 static inline float64_t log(float64_t v)
00355 {
00356 return ::log(v);
00357 }
00358
00359 template <class T>
00360 static void transpose_matrix(
00361 T*& matrix, int32_t& num_feat, int32_t& num_vec)
00362 {
00363 T* transposed=new T[num_vec*num_feat];
00364 for (int32_t i=0; i<num_vec; i++)
00365 {
00366 for (int32_t j=0; j<num_feat; j++)
00367 transposed[i+j*num_vec]=matrix[i*num_feat+j];
00368 }
00369
00370 delete[] matrix;
00371 matrix=transposed;
00372
00373 CMath::swap(num_feat, num_vec);
00374 }
00375
00376 #ifdef HAVE_LAPACK
00379 static float64_t* pinv(
00380 float64_t* matrix, int32_t rows, int32_t cols,
00381 float64_t* target=NULL);
00382
00383
00384
00385
00386 static inline void dgemm(
00387 double alpha, const double* A, int rows, int cols,
00388 CBLAS_TRANSPOSE transposeA, double *B, int cols_B,
00389 CBLAS_TRANSPOSE transposeB, double beta, double *C)
00390 {
00391 cblas_dgemm(CblasColMajor, transposeA, transposeB, rows, cols, cols_B,
00392 alpha, A, cols, B, cols_B, beta, C, cols);
00393 }
00394
00395
00396 static inline void dgemv(
00397 double alpha, const double *A, int rows, int cols,
00398 const CBLAS_TRANSPOSE transposeA, const double* X, double beta,
00399 double* Y)
00400 {
00401 cblas_dgemv(CblasColMajor, transposeA,
00402 rows, cols, alpha, A, cols,
00403 X, 1, beta, Y, 1);
00404 }
00405 #endif
00406
00407 static inline int64_t factorial(int32_t n)
00408 {
00409 int64_t res=1;
00410 for (int i=2; i<=n; i++)
00411 res*=i ;
00412 return res ;
00413 }
00414
00415 static void init_random(uint32_t initseed=0)
00416 {
00417 if (initseed==0)
00418 {
00419 struct timeval tv;
00420 gettimeofday(&tv, NULL);
00421 seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
00422 }
00423 else
00424 seed=initseed;
00425 #if !defined(CYGWIN) && !defined(__INTERIX)
00426
00427
00428 initstate(seed, CMath::rand_state, RNG_SEED_SIZE);
00429 #endif
00430 }
00431
00432 static inline int64_t random()
00433 {
00434 #if defined(CYGWIN) || defined(__INTERIX)
00435 return rand();
00436 #else
00437 return ::random();
00438 #endif
00439 }
00440
00441 static inline int32_t random(int32_t min_value, int32_t max_value)
00442 {
00443 int32_t ret = min_value + (int32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00444 ASSERT(ret>=min_value && ret<=max_value);
00445 return ret ;
00446 }
00447
00448 static inline float32_t random(float32_t min_value, float32_t max_value)
00449 {
00450 float32_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00451
00452 if (ret<min_value || ret>max_value)
00453 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00454 ASSERT(ret>=min_value && ret<=max_value);
00455 return ret;
00456 }
00457
00458 static inline float64_t random(float64_t min_value, float64_t max_value)
00459 {
00460 float64_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00461
00462 if (ret<min_value || ret>max_value)
00463 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00464 ASSERT(ret>=min_value && ret<=max_value);
00465 return ret;
00466 }
00467
00468 template <class T>
00469 static T* clone_vector(const T* vec, int32_t len)
00470 {
00471 T* result = new T[len];
00472 for (int32_t i=0; i<len; i++)
00473 result[i]=vec[i];
00474
00475 return result;
00476 }
00477 template <class T>
00478 static void fill_vector(T* vec, int32_t len, T value)
00479 {
00480 for (int32_t i=0; i<len; i++)
00481 vec[i]=value;
00482 }
00483 template <class T>
00484 static void range_fill_vector(T* vec, int32_t len, T start=0)
00485 {
00486 for (int32_t i=0; i<len; i++)
00487 vec[i]=i+start;
00488 }
00489
00490 template <class T>
00491 static void random_vector(T* vec, int32_t len, T min_value, T max_value)
00492 {
00493 for (int32_t i=0; i<len; i++)
00494 vec[i]=CMath::random(min_value, max_value);
00495 }
00496
00497 static inline int32_t* randperm(int32_t n)
00498 {
00499 int32_t* perm = new int32_t[n];
00500
00501 if (!perm)
00502 return NULL;
00503 for (int32_t i = 0; i < n; i++)
00504 perm[i] = i;
00505 for (int32_t i = 0; i < n; i++)
00506 swap(perm[random(0, n - 1)], perm[i]);
00507 return perm;
00508 }
00509
00510 static inline int64_t nchoosek(int32_t n, int32_t k)
00511 {
00512 int64_t res=1;
00513
00514 for (int32_t i=n-k+1; i<=n; i++)
00515 res*=i;
00516
00517 return res/factorial(k);
00518 }
00519
00521 template <class T>
00522 static inline void vec1_plus_scalar_times_vec2(T* vec1,
00523 T scalar, const T* vec2, int32_t n)
00524 {
00525 for (int32_t i=0; i<n; i++)
00526 vec1[i]+=scalar*vec2[i];
00527 }
00528
00530 static inline float64_t dot(const bool* v1, const bool* v2, int32_t n)
00531 {
00532 float64_t r=0;
00533 for (int32_t i=0; i<n; i++)
00534 r+=((v1[i]) ? 1 : 0) * ((v2[i]) ? 1 : 0);
00535 return r;
00536 }
00537
00539 static inline floatmax_t dot(const floatmax_t* v1, const floatmax_t* v2, int32_t n)
00540 {
00541 floatmax_t r=0;
00542 for (int32_t i=0; i<n; i++)
00543 r+=v1[i]*v2[i];
00544 return r;
00545 }
00546
00548 static inline float64_t dot(const float64_t* v1, const float64_t* v2, int32_t n)
00549 {
00550 float64_t r=0;
00551 #ifdef HAVE_LAPACK
00552 int32_t skip=1;
00553 r = cblas_ddot(n, v1, skip, v2, skip);
00554 #else
00555 for (int32_t i=0; i<n; i++)
00556 r+=v1[i]*v2[i];
00557 #endif
00558 return r;
00559 }
00560
00562 static inline float32_t dot(
00563 const float32_t* v1, const float32_t* v2, int32_t n)
00564 {
00565 float64_t r=0;
00566 #ifdef HAVE_LAPACK
00567 int32_t skip=1;
00568 r = cblas_sdot(n, v1, skip, v2, skip);
00569 #else
00570 for (int32_t i=0; i<n; i++)
00571 r+=v1[i]*v2[i];
00572 #endif
00573 return r;
00574 }
00575
00577 static inline float64_t dot(
00578 const uint64_t* v1, const uint64_t* v2, int32_t n)
00579 {
00580 float64_t r=0;
00581 for (int32_t i=0; i<n; i++)
00582 r+=((float64_t) v1[i])*v2[i];
00583
00584 return r;
00585 }
00587 static inline float64_t dot(
00588 const int64_t* v1, const int64_t* v2, int32_t n)
00589 {
00590 float64_t r=0;
00591 for (int32_t i=0; i<n; i++)
00592 r+=((float64_t) v1[i])*v2[i];
00593
00594 return r;
00595 }
00596
00598 static inline float64_t dot(
00599 const int32_t* v1, const int32_t* v2, int32_t n)
00600 {
00601 float64_t r=0;
00602 for (int32_t i=0; i<n; i++)
00603 r+=((float64_t) v1[i])*v2[i];
00604
00605 return r;
00606 }
00607
00609 static inline float64_t dot(
00610 const uint32_t* v1, const uint32_t* v2, int32_t n)
00611 {
00612 float64_t r=0;
00613 for (int32_t i=0; i<n; i++)
00614 r+=((float64_t) v1[i])*v2[i];
00615
00616 return r;
00617 }
00618
00620 static inline float64_t dot(
00621 const uint16_t* v1, const uint16_t* v2, int32_t n)
00622 {
00623 float64_t r=0;
00624 for (int32_t i=0; i<n; i++)
00625 r+=((float64_t) v1[i])*v2[i];
00626
00627 return r;
00628 }
00629
00631 static inline float64_t dot(
00632 const int16_t* v1, const int16_t* v2, int32_t n)
00633 {
00634 float64_t r=0;
00635 for (int32_t i=0; i<n; i++)
00636 r+=((float64_t) v1[i])*v2[i];
00637
00638 return r;
00639 }
00640
00642 static inline float64_t dot(
00643 const char* v1, const char* v2, int32_t n)
00644 {
00645 float64_t r=0;
00646 for (int32_t i=0; i<n; i++)
00647 r+=((float64_t) v1[i])*v2[i];
00648
00649 return r;
00650 }
00651
00653 static inline float64_t dot(
00654 const uint8_t* v1, const uint8_t* v2, int32_t n)
00655 {
00656 float64_t r=0;
00657 for (int32_t i=0; i<n; i++)
00658 r+=((float64_t) v1[i])*v2[i];
00659
00660 return r;
00661 }
00662
00664 static inline float64_t dot(
00665 const float64_t* v1, const char* v2, int32_t n)
00666 {
00667 float64_t r=0;
00668 for (int32_t i=0; i<n; i++)
00669 r+=((float64_t) v1[i])*v2[i];
00670
00671 return r;
00672 }
00673
00675 template <class T>
00676 static inline void add(
00677 T* target, T alpha, const T* v1, T beta, const T* v2,
00678 int32_t len)
00679 {
00680 for (int32_t i=0; i<len; i++)
00681 target[i]=alpha*v1[i]+beta*v2[i];
00682 }
00683
00685 template <class T>
00686 static inline void add_scalar(T alpha, T* vec, int32_t len)
00687 {
00688 for (int32_t i=0; i<len; i++)
00689 vec[i]+=alpha;
00690 }
00691
00693 template <class T>
00694 static inline void scale_vector(T alpha, T* vec, int32_t len)
00695 {
00696 for (int32_t i=0; i<len; i++)
00697 vec[i]*=alpha;
00698 }
00699
00701 template <class T>
00702 static inline T sum(T* vec, int32_t len)
00703 {
00704 T result=0;
00705 for (int32_t i=0; i<len; i++)
00706 result+=vec[i];
00707
00708 return result;
00709 }
00710
00712 template <class T>
00713 static inline T max(T* vec, int32_t len)
00714 {
00715 ASSERT(len>0);
00716 T maxv=vec[0];
00717
00718 for (int32_t i=1; i<len; i++)
00719 maxv=CMath::max(vec[i], maxv);
00720
00721 return maxv;
00722 }
00723
00725 template <class T>
00726 static inline T sum_abs(T* vec, int32_t len)
00727 {
00728 T result=0;
00729 for (int32_t i=0; i<len; i++)
00730 result+=CMath::abs(vec[i]);
00731
00732 return result;
00733 }
00734
00736 template <class T>
00737 static inline bool fequal(T x, T y, float64_t precision=1e-6)
00738 {
00739 return CMath::abs(x-y)<precision;
00740 }
00741
00742 static inline float64_t mean(float64_t* vec, int32_t len)
00743 {
00744 ASSERT(vec);
00745 ASSERT(len>0);
00746
00747 float64_t mean=0;
00748 for (int32_t i=0; i<len; i++)
00749 mean+=vec[i];
00750 return mean/len;
00751 }
00752
00753 static inline float64_t trace(
00754 float64_t* mat, int32_t cols, int32_t rows)
00755 {
00756 float64_t trace=0;
00757 for (int32_t i=0; i<rows; i++)
00758 trace+=mat[i*cols+i];
00759 return trace;
00760 }
00761
00765 static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
00766 static void sort(float64_t *a, int32_t*idx, int32_t N);
00767
00768
00769
00770
00771
00772
00773
00775 template <class T>
00776 inline static void radix_sort(T* array, int32_t size)
00777 {
00778 radix_sort_helper(array,size,0);
00779 }
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
00787 template <class T>
00788 static void radix_sort_helper(T* array, int32_t size, uint16_t i)
00789 {
00790 static size_t count[256], nc, cmin;
00791 T *ak;
00792 uint8_t c=0;
00793 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
00794 T *an, *aj, *pile[256];
00795 size_t *cp, cmax;
00796
00797
00798 sp = s;
00799 radix_push(array, size, i);
00800
00801
00802 while (sp>s) {
00803 radix_pop(array, size, i);
00804 an = array + size;
00805
00806
00807 if (nc == 0) {
00808 cmin = 0xff;
00809 for (ak = array; ak < an; ak++) {
00810 c = byte(*ak, i);
00811 count[c]++;
00812 if (count[c] == 1) {
00813
00814 if (c < cmin)
00815 cmin = c;
00816 nc++;
00817 }
00818 }
00819
00820
00821 if (sp + nc > s + RADIX_STACK_SIZE) {
00822 radix_sort_helper(array, size, i);
00823 continue;
00824 }
00825 }
00826
00827
00828
00829
00830
00831
00832
00833 olds = bigs = sp;
00834 cmax = 2;
00835 ak = array;
00836 pile[0xff] = an;
00837 for (cp = count + cmin; nc > 0; cp++) {
00838
00839 while (*cp == 0)
00840 cp++;
00841
00842 if (*cp > 1) {
00843
00844 if (*cp > cmax) {
00845 cmax = *cp;
00846 bigs = sp;
00847 }
00848
00849 if (i < sizeof(T)-1)
00850 radix_push(ak, *cp, (uint16_t) (i + 1));
00851 }
00852 pile[cp - count] = ak += *cp;
00853 nc--;
00854 }
00855
00856
00857 swap(*olds, *bigs);
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871 aj = array;
00872 while (aj<an)
00873 {
00874 T r;
00875
00876 for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
00877 swap(*ak, r);
00878
00879 *aj = r;
00880 aj += count[c];
00881 count[c] = 0;
00882 }
00883 }
00884 }
00885
00888 template <class T>
00889 static void insertion_sort(T* output, int32_t size)
00890 {
00891 for (int32_t i=0; i<size-1; i++)
00892 {
00893 int32_t j=i-1;
00894 T value=output[i];
00895 while (j >= 0 && output[j] > value)
00896 {
00897 output[j+1] = output[j];
00898 j--;
00899 }
00900 output[j+1]=value;
00901 }
00902 }
00903
00904
00907 template <class T>
00908 static void qsort(T* output, int32_t size)
00909 {
00910 if (size==2)
00911 {
00912 if (output[0] > output [1])
00913 swap(output[0],output[1]);
00914 return;
00915 }
00916
00917 T split=output[size/2];
00918
00919 int32_t left=0;
00920 int32_t right=size-1;
00921
00922 while (left<=right)
00923 {
00924 while (output[left] < split)
00925 left++;
00926 while (output[right] > split)
00927 right--;
00928
00929 if (left<=right)
00930 {
00931 swap(output[left],output[right]);
00932 left++;
00933 right--;
00934 }
00935 }
00936
00937 if (right+1> 1)
00938 qsort(output,right+1);
00939
00940 if (size-left> 1)
00941 qsort(&output[left],size-left);
00942 }
00943
00945 template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
00946 {
00947 ASSERT(width>=0);
00948 for (int i=0; i<width; i++)
00949 {
00950 T mask = ((T) 1)<<(sizeof(T)*8-1);
00951 while (mask)
00952 {
00953 if (mask & word)
00954 SG_SPRINT("1");
00955 else
00956 SG_SPRINT("0");
00957
00958 mask>>=1;
00959 }
00960 }
00961 }
00962
00964 template <class T> static void display_vector(
00965 const T* vector, int32_t n, const char* name="vector");
00966
00968 template <class T> static void display_matrix(
00969 const T* matrix, int32_t rows, int32_t cols, const char* name="matrix");
00970
00976 template <class T1,class T2>
00977 static void qsort_index(T1* output, T2* index, uint32_t size);
00978
00984 template <class T1,class T2>
00985 static void qsort_backward_index(
00986 T1* output, T2* index, int32_t size);
00987
00995 template <class T1,class T2>
00996 inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
00997 {
00998 int32_t n=0;
00999 thread_qsort<T1,T2> t;
01000 t.output=output;
01001 t.index=index;
01002 t.size=size;
01003 t.qsort_threads=&n;
01004 t.sort_limit=limit;
01005 t.num_threads=n_threads;
01006 parallel_qsort_index<T1,T2>(&t);
01007 }
01008
01009
01010 template <class T1,class T2>
01011 static void* parallel_qsort_index(void* p);
01012
01013
01014
01015
01016 template <class T>
01017 static void min(float64_t* output, T* index, int32_t size);
01018
01019
01020
01021 template <class T>
01022 static void nmin(
01023 float64_t* output, T* index, int32_t size, int32_t n);
01024
01025
01026
01027 template <class T>
01028 static int32_t unique(T* output, int32_t size)
01029 {
01030 qsort(output, size);
01031 int32_t i,j=0 ;
01032 for (i=0; i<size; i++)
01033 if (i==0 || output[i]!=output[i-1])
01034 output[j++]=output[i];
01035 return j ;
01036 }
01037
01038
01039
01040 template <class T>
01041 static int32_t binary_search_helper(T* output, int32_t size, T elem)
01042 {
01043 int32_t start=0;
01044 int32_t end=size-1;
01045
01046 if (size<1)
01047 return 0;
01048
01049 while (start<end)
01050 {
01051 int32_t middle=(start+end)/2;
01052
01053 if (output[middle]>elem)
01054 end=middle-1;
01055 else if (output[middle]<elem)
01056 start=middle+1;
01057 else
01058 return middle;
01059 }
01060
01061 return start;
01062 }
01063
01064
01065
01066 template <class T>
01067 static inline int32_t binary_search(T* output, int32_t size, T elem)
01068 {
01069 int32_t ind = binary_search_helper(output, size, elem);
01070 if (ind >= 0 && output[ind] == elem)
01071 return ind;
01072 return -1;
01073 }
01074
01075
01076
01077
01078
01079 template <class T>
01080 static int32_t binary_search_max_lower_equal(
01081 T* output, int32_t size, T elem)
01082 {
01083 int32_t ind = binary_search_helper(output, size, elem);
01084
01085 if (output[ind]<=elem)
01086 return ind;
01087 if (ind>0 && output[ind-1] <= elem)
01088 return ind-1;
01089 return -1;
01090 }
01091
01094 static float64_t Align(
01095 char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
01096
01101 static int32_t calcroc(
01102 float64_t* fp, float64_t* tp, float64_t* output, int32_t* label,
01103 int32_t& size, int32_t& possize, int32_t& negsize,
01104 float64_t& tresh, FILE* rocfile);
01106
01109 static float64_t mutual_info(float64_t* p1, float64_t* p2, int32_t len);
01110
01113 static float64_t relative_entropy(
01114 float64_t* p, float64_t* q, int32_t len);
01115
01117 static float64_t entropy(float64_t* p, int32_t len);
01118
01120 inline static uint32_t get_seed()
01121 {
01122 return CMath::seed;
01123 }
01124
01126 inline static int is_finite(double f)
01127 {
01128 #if defined(isfinite) && !defined(SUNOS)
01129 return isfinite(f);
01130 #else
01131 return finite(f);
01132 #endif
01133 }
01134
01136 inline static int is_infinity(double f)
01137 {
01138 #ifdef SUNOS
01139 if (fpclass(f) == FP_NINF || fpclass(f) == FP_PINF)
01140 return 1;
01141 else
01142 return 0;
01143 #else
01144 return isinf(f);
01145 #endif
01146 }
01147
01149 inline static int is_nan(double f)
01150 {
01151 #ifdef SUNOS
01152 return isnand(f);
01153 #else
01154 return isnan(f);
01155 #endif
01156 }
01157
01158
01169 #ifdef USE_LOGCACHE
01170 static inline float64_t logarithmic_sum(float64_t p, float64_t q)
01171 {
01172 float64_t diff;
01173
01174 if (!CMath::finite(p))
01175 return q;
01176
01177 if (!CMath::finite(q))
01178 {
01179 SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
01180 return NAN;
01181 }
01182 diff = p - q;
01183 if (diff > 0)
01184 return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
01185 return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
01186 }
01187
01189 static void init_log_table();
01190
01192 static int32_t determine_logrange();
01193
01195 static int32_t determine_logaccuracy(int32_t range);
01196 #else
01197 static inline float64_t logarithmic_sum(
01198 float64_t p, float64_t q)
01199 {
01200 float64_t diff;
01201
01202 if (!CMath::is_finite(p))
01203 return q;
01204 if (!CMath::is_finite(q))
01205 return p;
01206 diff = p - q;
01207 if (diff > 0)
01208 return diff > LOGRANGE? p : p + log(1 + exp(-diff));
01209 return -diff > LOGRANGE? q : q + log(1 + exp(diff));
01210 }
01211 #endif
01212 #ifdef LOG_SUM_ARRAY
01213
01218 static inline float64_t logarithmic_sum_array(
01219 float64_t *p, int32_t len)
01220 {
01221 if (len<=2)
01222 {
01223 if (len==2)
01224 return logarithmic_sum(p[0],p[1]) ;
01225 if (len==1)
01226 return p[0];
01227 return -INFTY ;
01228 }
01229 else
01230 {
01231 register float64_t *pp=p ;
01232 if (len%2==1) pp++ ;
01233 for (register int32_t j=0; j < len>>1; j++)
01234 pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
01235 }
01236 return logarithmic_sum_array(p,len%2+len>>1) ;
01237 }
01238 #endif
01239
01240
01242 inline virtual const char* get_name() const { return "Mathematics"; }
01243 public:
01246
01247 static const float64_t INFTY;
01248 static const float64_t ALMOST_INFTY;
01249
01251 static const float64_t ALMOST_NEG_INFTY;
01252
01254 static int32_t LOGRANGE;
01255
01257 static uint32_t seed;
01258 static char* rand_state;
01259
01260 #ifdef USE_LOGCACHE
01261
01263 static int32_t LOGACCURACY;
01265 protected:
01267 static float64_t* logtable;
01268 #endif
01269 };
01270
01271
01272 template <class T1,class T2>
01273 void* CMath::parallel_qsort_index(void* p)
01274 {
01275 struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
01276 T1* output=ps->output;
01277 T2* index=ps->index;
01278 uint32_t size=ps->size;
01279 int32_t* qsort_threads=ps->qsort_threads;
01280 int32_t sort_limit=ps->sort_limit;
01281 int32_t num_threads=ps->num_threads;
01282
01283 if (size==2)
01284 {
01285 if (output[0] > output [1])
01286 {
01287 swap(output[0], output[1]);
01288 swap(index[0], index[1]);
01289 }
01290 return NULL;
01291 }
01292
01293 T1 split=output[size/2];
01294
01295 int32_t left=0;
01296 int32_t right=size-1;
01297
01298 while (left<=right)
01299 {
01300 while (output[left] < split)
01301 left++;
01302 while (output[right] > split)
01303 right--;
01304
01305 if (left<=right)
01306 {
01307 swap(output[left], output[right]);
01308 swap(index[left], index[right]);
01309 left++;
01310 right--;
01311 }
01312 }
01313 bool lthread_start=false;
01314 bool rthread_start=false;
01315 pthread_t lthread;
01316 pthread_t rthread;
01317 struct thread_qsort<T1,T2> t1;
01318 struct thread_qsort<T1,T2> t2;
01319
01320 if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
01321 qsort_index(output,index,right+1);
01322 else if (right+1> 1)
01323 {
01324 (*qsort_threads)++;
01325 lthread_start=true;
01326 t1.output=output;
01327 t1.index=index;
01328 t1.size=right+1;
01329 t1.qsort_threads=qsort_threads;
01330 t1.sort_limit=sort_limit;
01331 t1.num_threads=num_threads;
01332 if (pthread_create(<hread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
01333 {
01334 lthread_start=false;
01335 (*qsort_threads)--;
01336 qsort_index(output,index,right+1);
01337 }
01338 }
01339
01340
01341 if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
01342 qsort_index(&output[left],&index[left], size-left);
01343 else if (size-left> 1)
01344 {
01345 (*qsort_threads)++;
01346 rthread_start=true;
01347 t2.output=&output[left];
01348 t2.index=&index[left];
01349 t2.size=size-left;
01350 t2.qsort_threads=qsort_threads;
01351 t2.sort_limit=sort_limit;
01352 t2.num_threads=num_threads;
01353 if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
01354 {
01355 rthread_start=false;
01356 (*qsort_threads)--;
01357 qsort_index(&output[left],&index[left], size-left);
01358 }
01359 }
01360
01361 if (lthread_start)
01362 {
01363 pthread_join(lthread, NULL);
01364 (*qsort_threads)--;
01365 }
01366
01367 if (rthread_start)
01368 {
01369 pthread_join(rthread, NULL);
01370 (*qsort_threads)--;
01371 }
01372
01373 return NULL;
01374 }
01375
01376 template <class T1,class T2>
01377 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
01378 {
01379 if (size==2)
01380 {
01381 if (output[0] > output [1])
01382 {
01383 swap(output[0],output[1]);
01384 swap(index[0],index[1]);
01385 }
01386 return;
01387 }
01388
01389 T1 split=output[size/2];
01390
01391 int32_t left=0;
01392 int32_t right=size-1;
01393
01394 while (left<=right)
01395 {
01396 while (output[left] < split)
01397 left++;
01398 while (output[right] > split)
01399 right--;
01400
01401 if (left<=right)
01402 {
01403 swap(output[left],output[right]);
01404 swap(index[left],index[right]);
01405 left++;
01406 right--;
01407 }
01408 }
01409
01410 if (right+1> 1)
01411 qsort_index(output,index,right+1);
01412
01413 if (size-left> 1)
01414 qsort_index(&output[left],&index[left], size-left);
01415 }
01416
01417 template <class T1,class T2>
01418 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
01419 {
01420 if (size==2)
01421 {
01422 if (output[0] < output [1])
01423 {
01424 swap(output[0],output[1]);
01425 swap(index[0],index[1]);
01426 }
01427 return;
01428 }
01429
01430
01431 T1 split=output[size/2];
01432
01433 int32_t left=0;
01434 int32_t right=size-1;
01435
01436 while (left<=right)
01437 {
01438 while (output[left] > split)
01439 left++;
01440 while (output[right] < split)
01441 right--;
01442
01443 if (left<=right)
01444 {
01445 swap(output[left],output[right]);
01446 swap(index[left],index[right]);
01447 left++;
01448 right--;
01449 }
01450 }
01451
01452 if (right+1> 1)
01453 qsort_backward_index(output,index,right+1);
01454
01455 if (size-left> 1)
01456 qsort_backward_index(&output[left],&index[left], size-left);
01457 }
01458
01459 template <class T>
01460 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
01461 {
01462 if (6*n*size<13*size*CMath::log(size))
01463 for (int32_t i=0; i<n; i++)
01464 min(&output[i], &index[i], size-i) ;
01465 else
01466 qsort_index(output, index, size) ;
01467 }
01468
01469
01470 template <class T>
01471 void CMath::min(float64_t* output, T* index, int32_t size)
01472 {
01473 if (size<=1)
01474 return;
01475 float64_t min_elem=output[0];
01476 int32_t min_index=0;
01477 for (int32_t i=1; i<size; i++)
01478 {
01479 if (output[i]<min_elem)
01480 {
01481 min_index=i;
01482 min_elem=output[i];
01483 }
01484 }
01485 swap(output[0], output[min_index]);
01486 swap(index[0], index[min_index]);
01487 }
01488 }
01489 #endif