Mathematics.h

Go to the documentation of this file.
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) 1999-2009 Soeren Sonnenburg
00008  * Written (W) 1999-2008 Gunnar Raetsch
00009  * Written (W) 2007 Konrad Rieck
00010  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
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 /* Size of RNG seed */
00073 #define RNG_SEED_SIZE 256
00074 
00075 /* Maximum stack size */
00076 #define RADIX_STACK_SIZE        512
00077 
00078 /* Stack macros */
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                 // can't be a>=0?(a):(-a), because compiler complains about
00168                 // 'comparison always true' when T is unsigned
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             //fall back to double precision sqrt if sqrtl is not
00259             //available
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             //fall back to double precision pow if powl is not
00272             //available
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         //C := alpha*op( A )*op( B ) + beta*C
00358         //op( X ) = X   or   op( X ) = X',
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         //y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
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             //seed=42
00400             //SG_SPRINT("initializing random number generator with %d (seed size %d)\n", seed, RNG_SEED_SIZE);
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          * Inline function to extract the byte at position p (from left)
00743          * of an 64 bit integer. The function is somewhat identical to 
00744          * accessing an array of characters via [].
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                 /* Push initial array to stack */
00771                 sp = s;
00772                 radix_push(array, size, i);
00773 
00774                 /* Loop until all digits have been sorted */
00775                 while (sp>s) {
00776                     radix_pop(array, size, i);
00777                     an = array + size;
00778 
00779                     /* Make character histogram */
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                                 /* Determine smallest character */
00787                                 if (c < cmin)
00788                                     cmin = c;
00789                                 nc++;
00790                             }
00791                         }
00792 
00793                         /* Sort recursively if stack size too small */
00794                         if (sp + nc > s + RADIX_STACK_SIZE) {
00795                             radix_sort_helper(array, size, i);
00796                             continue;
00797                         }
00798                     }
00799 
00800                     /*
00801                      * Set pile[]; push incompletely sorted bins onto stack.
00802                      * pile[] = pointers to last out-of-place element in bins.
00803                      * Before permuting: pile[c-1] + count[c] = pile[c];
00804                      * during deal: pile[c] counts down to pile[c-1].
00805                      */
00806                     olds = bigs = sp;
00807                     cmax = 2;
00808                     ak = array;
00809                     pile[0xff] = an;
00810                     for (cp = count + cmin; nc > 0; cp++) {
00811                         /* Find next non-empty pile */
00812                         while (*cp == 0)
00813                             cp++;
00814                         /* Pile with several entries */
00815                         if (*cp > 1) {
00816                             /* Determine biggest pile */
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                     /* Play it safe -- biggest bin last. */
00830                     swap(*olds, *bigs);
00831 
00832                     /*
00833                      * Permute misplacements home. Already home: everything
00834                      * before aj, and in pile[c], items from pile[c] on.
00835                      * Inner loop:
00836                      *      r = next element to put in place;
00837                      *      ak = pile[r[i]] = location to put the next element.
00838                      *      aj = bottom of 1st disordered bin.
00839                      * Outer loop:
00840                      *      Once the 1st disordered bin is done, ie. aj >= ak,
00841                      *      aj<-aj + count[c] connects the bins in array linked list;
00842                      *      reset count[c].
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                 //T split=output[random(0,size-1)];
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         /* finds the smallest element in output and puts that element as the 
00988            first element  */
00989         template <class T>
00990             static void min(float64_t* output, T* index, int32_t size);
00991 
00992         /* finds the n smallest elements in output and puts these elements as the 
00993            first n elements  */
00994         template <class T>
00995             static void nmin(
00996                 float64_t* output, T* index, int32_t size, int32_t n);
00997 
00998         /* performs a inplace unique of a vector of type T using quicksort
00999          * returns the new number of elements */
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         /* finds an element in a sorted array via binary search
01012          * returns -1 if not found */
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         /* finds an element in a sorted array via binary search
01041          *     * returns -1 if not found */
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         /* finds an element in a sorted array via binary search 
01052          * if it exists, else the index the largest smaller element
01053          * is returned
01054          * note: a successor is not mandatory */
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 //implementations of template functions
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         //T1 split=output[random(0,size-1)];
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(&lthread, 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     //T1 split=output[random(0,size-1)];
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     //T1 split=output[random(0,size-1)];
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 /* move the smallest entry in the array to the beginning */
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

SHOGUN Machine Learning Toolbox - Documentation