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