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-2008 Soeren Sonnenburg
00008  * Written (W) 1999-2008 Gunnar Raetsch
00009  * Written (W) 2007 Konrad Rieck
00010  * Copyright (C) 1999-2008 Fraunhofer Institute FIRST and Max-Planck-Society
00011  */
00012 
00013 #ifndef __MATHMATICS_H_
00014 #define __MATHMATICS_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 //define finite for win32
00032 #ifdef _WIN32
00033 #include <float.h>
00034 #ifndef finite
00035 #define finite _finite
00036 #endif
00037 
00038 #ifndef isnan
00039 #define isnan _isnan
00040 #endif
00041 #endif
00042 
00043 //define finite/nan for CYGWIN
00044 #ifdef CYGWIN
00045 #ifndef finite
00046 #include <ieeefp.h>
00047 #endif
00048 #endif
00049 
00050 #ifndef NAN
00051 #include <stdlib.h>
00052 #define NAN (strtod("NAN",NULL))
00053 #endif
00054 
00055 #ifdef SUNOS
00056 extern "C" int  finite(double);
00057 #endif
00058 
00059 /* Size of RNG seed */
00060 #define RNG_SEED_SIZE 256
00061 
00062 /* Maximum stack size */
00063 #define RADIX_STACK_SIZE        512
00064 
00065 /* Stack macros */
00066 #define radix_push(a, n, i)     sp->sa = a, sp->sn = n, (sp++)->si = i
00067 #define radix_pop(a, n, i)      a = (--sp)->sa, n = sp->sn, i = sp->si
00068 
00070 template <class T> struct radix_stack_t
00071 {
00073     T *sa;
00075     size_t sn;
00077     uint16_t si;
00078 };
00079 
00081 
00083 struct thread_qsort
00084 {
00086     float64_t* output;
00088     int32_t* index;
00090     int32_t size;
00091 
00093     int32_t* qsort_threads;
00095     int32_t sort_limit;
00097     int32_t num_threads;
00098 };
00099 
00103 class CMath : public CSGObject
00104 {
00105     public:
00109 
00110         CMath();
00111 
00113         virtual ~CMath();
00115 
00119 
00121         //
00122         template <class T>
00123             static inline T min(T a, T b)
00124             {
00125                 return ((a)<=(b))?(a):(b);
00126             }
00127 
00129         template <class T>
00130             static inline T max(T a, T b) 
00131             {
00132                 return ((a)>=(b))?(a):(b);
00133             }
00134 
00136         template <class T>
00137             static inline T clamp(T value, T lb, T ub) 
00138             {
00139                 if (value<=lb)
00140                     return lb;
00141                 else if (value>=ub)
00142                     return ub;
00143                 else
00144                     return value;
00145             }
00146 
00148         template <class T>
00149             static inline T abs(T a)
00150             {
00151                 // can't be a>=0?(a):(-a), because compiler complains about
00152                 // 'comparison always true' when T is unsigned
00153                 if (a==0)
00154                     return 0;
00155                 else if (a>0)
00156                     return a;
00157                 else
00158                     return -a;
00159             }
00161 
00164 
00166         static uint32_t crc32(uint8_t *data, int32_t len);
00167 
00168         static inline float64_t round(float64_t d)
00169         {
00170             return ::floor(d+0.5);
00171         }
00172 
00173         static inline float64_t floor(float64_t d)
00174         {
00175             return ::floor(d);
00176         }
00177 
00178         static inline float64_t ceil(float64_t d)
00179         {
00180             return ::ceil(d);
00181         }
00182 
00184         template <class T>
00185             static inline T sign(T a)
00186             {
00187                 if (a==0)
00188                     return 0;
00189                 else return (a<0) ? (-1) : (+1);
00190             }
00191 
00193         template <class T>
00194             static inline void swap(T & a,T &b)
00195             {
00196                 T c=a;
00197                 a=b;
00198                 b=c;
00199             }
00200 
00202         template <class T>
00203             static inline T sq(T x)
00204             {
00205                 return x*x;
00206             }
00207 
00209         static inline float32_t sqrt(float32_t x)
00210         {
00211             return ::sqrtf(x);
00212         }
00213 
00215         static inline float64_t sqrt(float64_t x)
00216         {
00217             return ::sqrt(x);
00218         }
00219 
00221         static inline float128_t sqrt(float128_t x)
00222         {
00223             //fall back to double precision sqrt if sqrtl is not
00224             //available
00225 #ifdef HAVE_SQRTL
00226             return ::sqrtl(x);
00227 #else
00228             return ::sqrt(x);
00229 #endif
00230         }
00231 
00232 
00234         static inline float128_t powl(float128_t x, float128_t n)
00235         {
00236             //fall back to double precision pow if powl is not
00237             //available
00238 #ifdef HAVE_POWL
00239             return ::powl((long double) x, (long double) n);
00240 #else
00241             return ::pow((double) x, (double) n);
00242 #endif
00243         }
00244 
00245         static inline int32_t pow(int32_t x, int32_t n)
00246         {
00247             ASSERT(n>=0);
00248             int32_t result=1;
00249             while (n--)
00250                 result*=x;
00251 
00252             return result;
00253         }
00254 
00255         static inline float64_t pow(float64_t x, int32_t n)
00256         {
00257             ASSERT(n>=0);
00258             float64_t result=1;
00259             while (n--)
00260                 result*=x;
00261 
00262             return result;
00263         }
00264 
00265         static inline float64_t pow(float64_t x, float64_t n)
00266         {
00267             return ::pow((double) x, (double) n);
00268         }
00269 
00270         static inline float64_t log10(float64_t v)
00271         {
00272             return ::log(v)/::log(10.0);
00273         }
00274 
00275 #ifdef HAVE_LOG2
00276         static inline float64_t log2(float64_t v)
00277         {
00278             return ::log2(v);
00279         }
00280 #else
00281         static inline float64_t log2(float64_t v)
00282         {
00283             return ::log(v)/::log(2.0);
00284         }
00285 #endif //HAVE_LOG2
00286 
00287         static inline float64_t log(float64_t v)
00288         {
00289             return ::log(v);
00290         }
00291 
00292         template <class T>
00293         static void transpose_matrix(
00294             T*& matrix, int32_t& num_feat, int32_t& num_vec)
00295         {
00296             T* transposed=new T[num_vec*num_feat];
00297             for (int32_t i=0; i<num_vec; i++)
00298             {
00299                 for (int32_t j=0; j<num_feat; j++)
00300                     transposed[i+j*num_vec]=matrix[i*num_feat+j];
00301             }
00302 
00303             delete[] matrix;
00304             matrix=transposed;
00305 
00306             CMath::swap(num_feat, num_vec);
00307         }
00308 
00309 #ifdef HAVE_LAPACK
00312         static float64_t* pinv(
00313             float64_t* matrix, int32_t rows, int32_t cols,
00314             float64_t* target=NULL);
00315 
00316 
00317         //C := alpha*op( A )*op( B ) + beta*C
00318         //op( X ) = X   or   op( X ) = X',
00319         static inline void dgemm(
00320             double alpha, const double* A, int rows, int cols,
00321             CBLAS_TRANSPOSE transposeA, double *B, int cols_B,
00322             CBLAS_TRANSPOSE transposeB, double beta, double *C)
00323         {
00324             cblas_dgemm(CblasColMajor, transposeA, transposeB, rows, cols, cols_B,
00325                     alpha, A, cols, B, cols_B, beta, C, cols);
00326         }
00327 
00328         //y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
00329         static inline void dgemv(
00330             double alpha, const double *A, int rows, int cols,
00331             const CBLAS_TRANSPOSE transposeA, const double* X, double beta,
00332             double* Y)
00333         {
00334             cblas_dgemv(CblasColMajor, transposeA,
00335                     rows, cols, alpha, A, cols,
00336                     X, 1, beta, Y, 1);
00337         }
00338 #endif
00339 
00340         static inline int64_t factorial(int32_t n)
00341         {
00342             int64_t res=1;
00343             for (int i=2; i<=n; i++)
00344                 res*=i ;
00345             return res ;
00346         }
00347 
00348         static void init_random(uint32_t initseed=0)
00349         {
00350             if (initseed==0)
00351             {
00352                 struct timeval tv;
00353                 gettimeofday(&tv, NULL);
00354                 seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
00355             }
00356             else
00357                 seed=initseed;
00358 #if !defined(CYGWIN) && !defined(__INTERIX)
00359             //seed=42
00360             //SG_SPRINT("initializing random number generator with %d (seed size %d)\n", seed, RNG_SEED_SIZE);
00361             initstate(seed, CMath::rand_state, RNG_SEED_SIZE);
00362 #endif
00363         }
00364 
00365         static inline int64_t random()
00366         {
00367 #if defined(CYGWIN) || defined(__INTERIX)
00368             return rand();
00369 #else
00370             return ::random();
00371 #endif
00372         }
00373 
00374         static inline int32_t random(int32_t min_value, int32_t max_value)
00375         {
00376             int32_t ret = min_value + (int32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00377             ASSERT(ret>=min_value && ret<=max_value);
00378             return ret ;
00379         }
00380 
00381         static inline float32_t random(float32_t min_value, float32_t max_value)
00382         {
00383             float32_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00384 
00385             if (ret<min_value || ret>max_value)
00386                 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00387             ASSERT(ret>=min_value && ret<=max_value);
00388             return ret;
00389         }
00390 
00391         static inline float64_t random(float64_t min_value, float64_t max_value)
00392         {
00393             float64_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00394 
00395             if (ret<min_value || ret>max_value)
00396                 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00397             ASSERT(ret>=min_value && ret<=max_value);
00398             return ret;
00399         }
00400 
00401         template <class T>
00402             static void fill_vector(T* vec, int32_t len, T value)
00403             {
00404                 for (int32_t i=0; i<len; i++)
00405                     vec[i]=value;
00406             }
00407         template <class T>
00408             static void range_fill_vector(T* vec, int32_t len, T start=0)
00409             {
00410                 for (int32_t i=0; i<len; i++)
00411                     vec[i]=i+start;
00412             }
00413 
00414         template <class T>
00415             static void random_vector(T* vec, int32_t len, T min_value, T max_value)
00416             {
00417                 for (int32_t i=0; i<len; i++)
00418                     vec[i]=CMath::random(min_value, max_value);
00419             }
00420 
00421         static inline int32_t* randperm(int32_t n)
00422         {
00423             int32_t* perm = new int32_t[n];
00424 
00425             if (!perm)
00426                 return NULL;
00427             for (int32_t i = 0; i < n; i++)
00428                 perm[i] = i;
00429             for (int32_t i = 0; i < n; i++)
00430                 swap(perm[random(0, n - 1)], perm[i]);
00431             return perm;
00432         }
00433 
00434         static inline int64_t nchoosek(int32_t n, int32_t k)
00435         {
00436             int64_t res=1;
00437 
00438             for (int32_t i=n-k+1; i<=n; i++)
00439                 res*=i;
00440 
00441             return res/factorial(k);
00442         }
00443 
00445         template <class T>
00446             static inline void vec1_plus_scalar_times_vec2(T* vec1,
00447                     T scalar, const T* vec2, int32_t n)
00448             {
00449                 for (int32_t i=0; i<n; i++)
00450                     vec1[i]+=scalar*vec2[i];
00451             }
00452 
00454         static inline float64_t dot(const float64_t* v1, const float64_t* v2, int32_t n)
00455         {
00456             float64_t r=0;
00457 #ifdef HAVE_LAPACK
00458             int32_t skip=1;
00459             r = cblas_ddot(n, v1, skip, v2, skip);
00460 #else
00461             for (int32_t i=0; i<n; i++)
00462                 r+=v1[i]*v2[i];
00463 #endif
00464             return r;
00465         }
00466 
00468         static inline float32_t dot(
00469             const float32_t* v1, const float32_t* v2, int32_t n)
00470         {
00471             float64_t r=0;
00472 #ifdef HAVE_LAPACK
00473             int32_t skip=1;
00474             r = cblas_sdot(n, v1, skip, v2, skip);
00475 #else
00476             for (int32_t i=0; i<n; i++)
00477                 r+=v1[i]*v2[i];
00478 #endif
00479             return r;
00480         }
00481 
00483         static inline float64_t dot(
00484             const int64_t* v1, const int64_t* v2, int32_t n)
00485         {
00486             float64_t r=0;
00487             for (int32_t i=0; i<n; i++)
00488                 r+=((float64_t) v1[i])*v2[i];
00489 
00490             return r;
00491         }
00492 
00494         static inline float64_t dot(
00495             const uint16_t* v1, const uint16_t* v2, int32_t n)
00496         {
00497             float64_t r=0;
00498             for (int32_t i=0; i<n; i++)
00499                 r+=((float64_t) v1[i])*v2[i];
00500 
00501             return r;
00502         }
00503 
00505         static inline float64_t dot(
00506             const char* v1, const char* v2, int32_t n)
00507         {
00508             float64_t r=0;
00509             for (int32_t i=0; i<n; i++)
00510                 r+=((float64_t) v1[i])*v2[i];
00511 
00512             return r;
00513         }
00514 
00516         static inline float64_t dot(
00517             const uint8_t* v1, const uint8_t* v2, int32_t n)
00518         {
00519             float64_t r=0;
00520             for (int32_t i=0; i<n; i++)
00521                 r+=((float64_t) v1[i])*v2[i];
00522 
00523             return r;
00524         }
00525 
00527         static inline float64_t dot(
00528             const float64_t* v1, const char* v2, int32_t n)
00529         {
00530             float64_t r=0;
00531             for (int32_t i=0; i<n; i++)
00532                 r+=((float64_t) v1[i])*v2[i];
00533 
00534             return r;
00535         }
00536 
00538         template <class T>
00539             static inline void add(
00540                 T* target, T alpha, const T* v1, T beta, const T* v2,
00541                 int32_t len)
00542             {
00543                 for (int32_t i=0; i<len; i++)
00544                     target[i]=alpha*v1[i]+beta*v2[i];
00545             }
00546 
00548         template <class T>
00549             static inline void add_scalar(T alpha, T* vec, int32_t len)
00550             {
00551                 for (int32_t i=0; i<len; i++)
00552                     vec[i]+=alpha;
00553             }
00554 
00556         template <class T>
00557             static inline void scale_vector(T alpha, T* vec, int32_t len)
00558             {
00559                 for (int32_t i=0; i<len; i++)
00560                     vec[i]*=alpha;
00561             }
00562 
00564         template <class T>
00565             static inline T sum(T* vec, int32_t len)
00566             {
00567                 T result=0;
00568                 for (int32_t i=0; i<len; i++)
00569                     result+=vec[i];
00570 
00571                 return result;
00572             }
00573 
00575         template <class T>
00576             static inline T max(T* vec, int32_t len)
00577             {
00578                 ASSERT(len>0);
00579                 T maxv=vec[0];
00580 
00581                 for (int32_t i=1; i<len; i++)
00582                     maxv=CMath::max(vec[i], maxv);
00583 
00584                 return maxv;
00585             }
00586 
00588         template <class T>
00589             static inline T sum_abs(T* vec, int32_t len)
00590             {
00591                 T result=0;
00592                 for (int32_t i=0; i<len; i++)
00593                     result+=CMath::abs(vec[i]);
00594 
00595                 return result;
00596             }
00597 
00598         static inline float64_t mean(float64_t* vec, int32_t len)
00599         {
00600             ASSERT(vec);
00601             ASSERT(len>0);
00602 
00603             float64_t mean=0;
00604             for (int32_t i=0; i<len; i++)
00605                 mean+=vec[i];
00606             return mean/len;
00607         }
00608 
00609         static inline float64_t trace(
00610             float64_t* mat, int32_t cols, int32_t rows)
00611         {
00612             float64_t trace=0;
00613             for (int32_t i=0; i<rows; i++)
00614                 trace+=mat[i*cols+i];
00615             return trace;
00616         }
00617 
00621         static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
00622         static void sort(float64_t *a, int32_t*idx, int32_t N);
00623 
00624         /*
00625          * Inline function to extract the byte at position p (from left)
00626          * of an 64 bit integer. The function is somewhat identical to 
00627          * accessing an array of characters via [].
00628          */
00629 
00631         template <class T>
00632             inline static void radix_sort(T* array, size_t size)
00633             {
00634                 radix_sort_helper(array,size,0);
00635             }
00636 
00637         template <class T>
00638             static inline uint8_t byte(T word, uint16_t p)
00639             {
00640                 return (word >> (sizeof(T)-p-1) * 8) & 0xff;
00641             }
00642 
00643         template <class T>
00644             static void radix_sort_helper(T* array, size_t size, uint16_t i)
00645             {
00646                 static size_t count[256], nc, cmin;
00647                 T *ak;
00648                 uint8_t c=0;
00649                 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
00650                 T *an, *aj, *pile[256];
00651                 size_t *cp, cmax;
00652 
00653                 /* Push initial array to stack */
00654                 sp = s;
00655                 radix_push(array, size, i);
00656 
00657                 /* Loop until all digits have been sorted */
00658                 while (sp>s) {
00659                     radix_pop(array, size, i);
00660                     an = array + size;
00661 
00662                     /* Make character histogram */
00663                     if (nc == 0) {
00664                         cmin = 0xff;
00665                         for (ak = array; ak < an; ak++) {
00666                             c = byte(*ak, i);
00667                             count[c]++;
00668                             if (count[c] == 1) {
00669                                 /* Determine smallest character */
00670                                 if (c < cmin)
00671                                     cmin = c;
00672                                 nc++;
00673                             }
00674                         }
00675 
00676                         /* Sort recursively if stack size too small */
00677                         if (sp + nc > s + RADIX_STACK_SIZE) {
00678                             radix_sort_helper(array, size, i);
00679                             continue;
00680                         }
00681                     }
00682 
00683                     /*
00684                      * Set pile[]; push incompletely sorted bins onto stack.
00685                      * pile[] = pointers to last out-of-place element in bins.
00686                      * Before permuting: pile[c-1] + count[c] = pile[c];
00687                      * during deal: pile[c] counts down to pile[c-1].
00688                      */
00689                     olds = bigs = sp;
00690                     cmax = 2;
00691                     ak = array;
00692                     pile[0xff] = an;
00693                     for (cp = count + cmin; nc > 0; cp++) {
00694                         /* Find next non-empty pile */
00695                         while (*cp == 0)
00696                             cp++;
00697                         /* Pile with several entries */
00698                         if (*cp > 1) {
00699                             /* Determine biggest pile */
00700                             if (*cp > cmax) {
00701                                 cmax = *cp;
00702                                 bigs = sp;
00703                             }
00704 
00705                             if (i < sizeof(T)-1)
00706                                 radix_push(ak, *cp, (uint16_t) (i + 1));
00707                         }
00708                         pile[cp - count] = ak += *cp;
00709                         nc--;
00710                     }
00711 
00712                     /* Play it safe -- biggest bin last. */
00713                     swap(*olds, *bigs);
00714 
00715                     /*
00716                      * Permute misplacements home. Already home: everything
00717                      * before aj, and in pile[c], items from pile[c] on.
00718                      * Inner loop:
00719                      *      r = next element to put in place;
00720                      *      ak = pile[r[i]] = location to put the next element.
00721                      *      aj = bottom of 1st disordered bin.
00722                      * Outer loop:
00723                      *      Once the 1st disordered bin is done, ie. aj >= ak,
00724                      *      aj<-aj + count[c] connects the bins in array linked list;
00725                      *      reset count[c].
00726                      */
00727                     aj = array;
00728                     while (aj<an)
00729                     {
00730                         T r;
00731 
00732                         for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
00733                             swap(*ak, r);
00734 
00735                         *aj = r;
00736                         aj += count[c];
00737                         count[c] = 0;
00738                     }
00739                 }
00740             }
00743         template <class T>
00744             static void insertion_sort(T* output, int32_t size)
00745             {
00746                 for (int32_t i=0; i<size-1; i++)
00747                 {
00748                     int32_t j=i-1;
00749                     T value=output[i];
00750                     while (j >= 0 && output[j] > value)
00751                     {
00752                         output[j+1] = output[j];
00753                         j--;
00754                     }
00755                     output[j+1]=value;
00756                 }
00757             }
00758 
00759 
00762         //template <class T>
00763         //static void qsort(T* output, int32_t size) ;
00764         template <class T>
00765             static void qsort(T* output, int32_t size)
00766             {
00767                 if (size==2)
00768                 {
00769                     if (output[0] > output [1])
00770                         swap(output[0],output[1]);
00771                     return;
00772                 }
00773                 T split=output[random(0,size-1)];
00774 
00775                 int32_t left=0;
00776                 int32_t right=size-1;
00777 
00778                 while (left<=right)
00779                 {
00780                     while (output[left] < split)
00781                         left++;
00782                     while (output[right] > split)
00783                         right--;
00784 
00785                     if (left<=right)
00786                     {
00787                         swap(output[left],output[right]);
00788                         left++;
00789                         right--;
00790                     }
00791                 }
00792 
00793                 if (right+1> 1)
00794                     qsort(output,right+1);
00795 
00796                 if (size-left> 1)
00797                     qsort(&output[left],size-left);
00798             }
00799 
00801         template <class T> static void display_vector(
00802             T* vector, int32_t n, const char* name="vector");
00803 
00805         template <class T> static void display_matrix(
00806             T* matrix, int32_t rows, int32_t cols, const char* name="matrix");
00807 
00813         template <class T1,class T2>
00814             static void qsort_index(T1* output, T2* index, uint32_t size);
00815 
00816         template <class T1,class T2>
00817             static void* parallel_qsort_index(void* p);
00818 
00824         template <class T1,class T2>
00825             static void qsort_backward_index(
00826                 T1* output, T2* index, int32_t size);
00827 
00828         /* finds the smallest element in output and puts that element as the 
00829            first element  */
00830         template <class T>
00831             static void min(float64_t* output, T* index, int32_t size);
00832 
00833         /* finds the n smallest elements in output and puts these elements as the 
00834            first n elements  */
00835         template <class T>
00836             static void nmin(
00837                 float64_t* output, T* index, int32_t size, int32_t n);
00838 
00839         /* performs a inplace unique of a vector of type T using quicksort
00840          * returns the new number of elements */
00841         template <class T>
00842             static int32_t unique(T* output, int32_t size)
00843             {
00844                 qsort(output, size);
00845                 int32_t i,j=0 ;
00846                 for (i=0; i<size; i++)
00847                     if (i==0 || output[i]!=output[i-1])
00848                         output[j++]=output[i];
00849                 return j ;
00850             }
00851 
00852         /* finds an element in a sorted array via binary search
00853          * returns -1 if not found */
00854         template <class T>
00855             static int32_t binary_search_helper(T* output, int32_t size, T elem)
00856             {
00857                 int32_t start=0;
00858                 int32_t end=size-1;
00859 
00860                 if (size<1)
00861                     return -1;
00862 
00863                 while (start<end)
00864                 {
00865                     int32_t middle=(start+end)/2; 
00866 
00867                     if (output[middle]>elem)
00868                         end=middle-1;
00869                     else if (output[middle]<elem)
00870                         start=middle+1;
00871                     else
00872                         return middle;
00873                 }
00874 
00875                 if (output[start]==elem)
00876                     return start;
00877                 else
00878                     return -1;
00879             }
00880 
00881         /* finds an element in a sorted array via binary search
00882          *     * returns -1 if not found */
00883         template <class T>
00884             static inline int32_t binary_search(T* output, int32_t size, T elem)
00885             {
00886                 int32_t ind = binary_search_helper(output, size, elem);
00887                 if (ind >= 0 && output[ind] == elem)
00888                     return ind;
00889                 return -1;
00890             }
00891 
00892         /* finds an element in a sorted array via binary search 
00893          * if it exists, else the index the largest smaller element
00894          * is returned
00895          * note: a successor is not mandatory */
00896         template <class T>
00897             static int32_t binary_search_max_lower_equal(
00898                 T* output, int32_t size, T elem)
00899             {
00900                 int32_t ind = binary_search_helper(output, size, elem);
00901 
00902                 if (output[ind]<=elem)
00903                     return ind;
00904                 if (ind>0 && output[ind-1] <= elem)
00905                     return ind-1;
00906                 return -1;
00907             }
00908 
00911         static float64_t Align(
00912             char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
00913 
00918         static int32_t calcroc(
00919             float64_t* fp, float64_t* tp, float64_t* output, int32_t* label,
00920             int32_t& size, int32_t& possize, int32_t& negsize,
00921             float64_t& tresh, FILE* rocfile);
00923 
00926         static float64_t mutual_info(float64_t* p1, float64_t* p2, int32_t len);
00927 
00930         static float64_t relative_entropy(
00931             float64_t* p, float64_t* q, int32_t len);
00932 
00934         static float64_t entropy(float64_t* p, int32_t len);
00935 
00937         inline static uint32_t get_seed()
00938         {
00939             return CMath::seed;
00940         }
00941 
00942 
00953 #ifdef USE_LOGCACHE
00954         static inline float64_t logarithmic_sum(float64_t p, float64_t q)
00955         {
00956             float64_t diff;
00957 
00958             if (!finite(p))
00959                 return q;
00960             if (!finite(q))
00961             {
00962                 SG_SWARNING(("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
00963                         return NAN;
00964                         }
00965                         diff = p - q;
00966                         if (diff > 0)
00967                         return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
00968                         return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
00969                         }
00970 
00972                         static void init_log_table();
00973 
00975                         static int32_t determine_logrange();
00976 
00978                         static int32_t determine_logaccuracy(int32_t range);
00979 #else
00980                         static inline float64_t logarithmic_sum(
00981                             float64_t p, float64_t q)
00982                         {
00983                         float64_t diff;
00984 
00985                         if (!finite(p))
00986                             return q;
00987                         if (!finite(q))
00988                             return p;
00989                         diff = p - q;
00990                         if (diff > 0)
00991                             return diff > LOGRANGE? p : p + log(1 + exp(-diff));
00992                         return -diff > LOGRANGE? q : q + log(1 + exp(diff));
00993                         }
00994 #endif
00995 #ifdef LOG_SUM_ARRAY
00996 
01001                 static inline float64_t logarithmic_sum_array(
01002                     float64_t *p, int32_t len)
01003                 {
01004                     if (len<=2)
01005                     {
01006                         if (len==2)
01007                             return logarithmic_sum(p[0],p[1]) ;
01008                         if (len==1)
01009                             return p[0];
01010                         return -INFTY ;
01011                     }
01012                     else
01013                     {
01014                         register float64_t *pp=p ;
01015                         if (len%2==1) pp++ ;
01016                         for (register int32_t j=0; j < len>>1; j++)
01017                             pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
01018                     }
01019                     return logarithmic_sum_array(p,len%2+len>>1) ;
01020                 } 
01021 #endif
01022 
01023     public:
01026 
01027                 static const float64_t INFTY;
01028                 static const float64_t ALMOST_INFTY;
01029 
01031                 static const float64_t ALMOST_NEG_INFTY;
01032 
01034                 static int32_t LOGRANGE;
01035 
01037                 static uint32_t seed;
01038 
01039 #ifdef USE_LOGCACHE 
01040 
01042                 static int32_t LOGACCURACY;
01044     protected:
01046                 static float64_t* logtable; 
01047 #endif
01048                 static char* rand_state;
01049 };
01050 
01051     template <class T1,class T2>
01052 void* CMath::parallel_qsort_index(void* p)
01053 {
01054     struct thread_qsort* ps=(thread_qsort*) p;
01055     T1* output=ps->output;
01056     T2* index=ps->index;
01057     int32_t size=ps->size;
01058     int32_t* qsort_threads=ps->qsort_threads;
01059     int32_t sort_limit=ps->sort_limit;
01060     int32_t num_threads=ps->num_threads;
01061 
01062     if (size==2)
01063     {
01064         if (output[0] > output [1])
01065         {
01066             swap(output[0], output[1]);
01067             swap(index[0], index[1]);
01068         }
01069         return NULL;
01070     }
01071     /*float64_t split=output[(((uint64_t) size)*rand())/(((uint64_t)RAND_MAX)+1)];*/
01072     float64_t split=output[size/2];
01073 
01074     int32_t left=0;
01075     int32_t right=size-1;
01076 
01077     while (left<=right)
01078     {
01079         while (output[left] < split)
01080             left++;
01081         while (output[right] > split)
01082             right--;
01083 
01084         if (left<=right)
01085         {
01086             swap(output[left], output[right]);
01087             swap(index[left], index[right]);
01088             left++;
01089             right--;
01090         }
01091     }
01092     bool lthread_start=false;
01093     bool rthread_start=false;
01094     pthread_t lthread;
01095     pthread_t rthread;
01096     struct thread_qsort t1;
01097     struct thread_qsort t2;
01098 
01099     if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
01100         qsort_index(output,index,right+1);
01101     else if (right+1> 1)
01102     {
01103         *qsort_threads++;
01104         lthread_start=true;
01105         t1.output=output;
01106         t1.index=index;
01107         t1.size=right+1;
01108         if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, (void*) &t1) != 0)
01109         {
01110             lthread_start=false;
01111             *qsort_threads--;
01112             qsort_index(output,index,right+1);
01113         }
01114     }
01115 
01116 
01117     if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
01118         qsort_index(&output[left],&index[left], size-left);
01119     else if (size-left> 1)
01120     {
01121         *qsort_threads++;
01122         rthread_start=true;
01123         t2.output=&output[left];
01124         t2.index=&index[left];
01125         t2.size=size-left;
01126         if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, (void*)&t2) != 0)
01127         {
01128             rthread_start=false;
01129             *qsort_threads--;
01130             qsort_index(&output[left],&index[left], size-left);
01131         }
01132     }
01133 
01134     if (lthread_start)
01135     {
01136         pthread_join(lthread, NULL);
01137         *qsort_threads--;
01138     }
01139 
01140     if (rthread_start)
01141     {
01142         pthread_join(rthread, NULL);
01143         *qsort_threads--;
01144     }
01145 
01146     return NULL;
01147 }
01148 
01149 
01150 //implementations of template functions
01151     template <class T1,class T2>
01152 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
01153 {
01154     if (size==2)
01155     {
01156         if (output[0] > output [1])
01157         {
01158             swap(output[0],output[1]);
01159             swap(index[0],index[1]);
01160         }
01161         return;
01162     }
01163     T1 split=output[random(0,size-1)];
01164 
01165     int32_t left=0;
01166     int32_t right=size-1;
01167 
01168     while (left<=right)
01169     {
01170         while (output[left] < split)
01171             left++;
01172         while (output[right] > split)
01173             right--;
01174 
01175         if (left<=right)
01176         {
01177             swap(output[left],output[right]);
01178             swap(index[left],index[right]);
01179             left++;
01180             right--;
01181         }
01182     }
01183 
01184     if (right+1> 1)
01185         qsort_index(output,index,right+1);
01186 
01187     if (size-left> 1)
01188         qsort_index(&output[left],&index[left], size-left);
01189 }
01190 
01191     template <class T1,class T2>
01192 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
01193 {
01194     if (size==2)
01195     {
01196         if (output[0] < output [1])
01197         {
01198             swap(output[0],output[1]);
01199             swap(index[0],index[1]);
01200         }
01201         return;
01202     }
01203 
01204     T1 split=output[random(0,size-1)];
01205 
01206     int32_t left=0;
01207     int32_t right=size-1;
01208 
01209     while (left<=right)
01210     {
01211         while (output[left] > split)
01212             left++;
01213         while (output[right] < split)
01214             right--;
01215 
01216         if (left<=right)
01217         {
01218             swap(output[left],output[right]);
01219             swap(index[left],index[right]);
01220             left++;
01221             right--;
01222         }
01223     }
01224 
01225     if (right+1> 1)
01226         qsort_backward_index(output,index,right+1);
01227 
01228     if (size-left> 1)
01229         qsort_backward_index(&output[left],&index[left], size-left);
01230 }
01231 
01232     template <class T> 
01233 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
01234 {
01235     if (6*n*size<13*size*CMath::log(size))
01236         for (int32_t i=0; i<n; i++)
01237             min(&output[i], &index[i], size-i) ;
01238     else
01239         qsort_index(output, index, size) ;
01240 }
01241 
01242 /* move the smallest entry in the array to the beginning */
01243     template <class T>
01244 void CMath::min(float64_t* output, T* index, int32_t size)
01245 {
01246     if (size<=1)
01247         return;
01248     float64_t min_elem=output[0];
01249     int32_t min_index=0;
01250     for (int32_t i=1; i<size; i++)
01251     {
01252         if (output[i]<min_elem)
01253         {
01254             min_index=i;
01255             min_elem=output[i];
01256         }
01257     }
01258     swap(output[0], output[min_index]);
01259     swap(index[0], index[min_index]);
01260 }
01261 #endif

SHOGUN Machine Learning Toolbox - Documentation