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 /* Maximum stack size */
00060 #define RADIX_STACK_SIZE        512
00061 
00062 /* Stack macros */
00063 #define radix_push(a, n, i)     sp->sa = a, sp->sn = n, (sp++)->si = i
00064 #define radix_pop(a, n, i)      a = (--sp)->sa, n = sp->sn, i = sp->si
00065 
00067 template <class T> struct radix_stack_t
00068 {
00070     T *sa;
00072     size_t sn;
00074     unsigned short si;
00075 };
00076 
00078 struct pair
00079 {
00081     int idx1;
00083     int idx2;
00084 };
00085 
00087 struct thread_qsort
00088 {
00090     DREAL* output;
00092     INT* index;
00094     INT size;
00095 
00097     INT* qsort_threads;
00099     INT sort_limit;
00101     INT num_threads;
00102 };
00103 
00107 class CMath : public CSGObject
00108 {
00109     public:
00113 
00114         CMath();
00115 
00117         virtual ~CMath();
00119 
00123 
00125         //
00126         template <class T>
00127             static inline T min(T a, T b)
00128             {
00129                 return ((a)<=(b))?(a):(b);
00130             }
00131 
00133         template <class T>
00134             static inline T max(T a, T b) 
00135             {
00136                 return ((a)>=(b))?(a):(b);
00137             }
00138 
00140         template <class T>
00141             static inline T clamp(T value, T lb, T ub) 
00142             {
00143                 if (value<=lb)
00144                     return lb;
00145                 else if (value>=ub)
00146                     return ub;
00147                 else
00148                     return value;
00149             }
00150 
00152         template <class T>
00153             static inline T abs(T a)
00154             {
00155                 // can't be a>=0?(a):(-a), because compiler complains about
00156                 // 'comparison always true' when T is unsigned
00157                 if (a==0)
00158                     return 0;
00159                 else if (a>0)
00160                     return a;
00161                 else
00162                     return -a;
00163             }
00165 
00168 
00170         static UINT crc32(BYTE *data, INT len);
00171 
00172         static inline DREAL round(DREAL d)
00173         {
00174             return ::floor(d+0.5);
00175         }
00176 
00177         static inline DREAL floor(DREAL d)
00178         {
00179             return ::floor(d);
00180         }
00181 
00182         static inline DREAL ceil(DREAL d)
00183         {
00184             return ::ceil(d);
00185         }
00186 
00188         template <class T>
00189             static inline T sign(T a)
00190             {
00191                 if (a==0)
00192                     return 0;
00193                 else return (a<0) ? (-1) : (+1);
00194             }
00195 
00197         template <class T>
00198             static inline void swap(T & a,T &b)
00199             {
00200                 T c=a;
00201                 a=b;
00202                 b=c;
00203             }
00204 
00206         template <class T>
00207             static inline T sq(T x)
00208             {
00209                 return x*x;
00210             }
00211 
00213         static inline SHORTREAL sqrt(SHORTREAL x)
00214         {
00215             return ::sqrtf(x);
00216         }
00217 
00219         static inline DREAL sqrt(DREAL x)
00220         {
00221             return ::sqrt(x);
00222         }
00223 
00225         static inline LONGREAL sqrt(LONGREAL x)
00226         {
00227             //fall back to double precision sqrt if sqrtl is not
00228             //available
00229 #ifdef HAVE_SQRTL
00230             return ::sqrtl(x);
00231 #else
00232             return ::sqrt(x);
00233 #endif
00234         }
00235 
00236 
00238         static inline long double powl(long double x, long double n)
00239         {
00240             //fall back to double precision pow if powl is not
00241             //available
00242 #ifdef HAVE_POWL
00243             return ::powl(x, n);
00244 #else
00245             return ::pow(x, n);
00246 #endif
00247         }
00248 
00249         static inline INT pow(INT x, INT n)
00250         {
00251             ASSERT(n>=0);
00252             INT result=1;
00253             while (n--)
00254                 result*=x;
00255 
00256             return result;
00257         }
00258 
00259         static inline DREAL pow(DREAL x, DREAL n)
00260         {
00261             return ::pow(x, n);
00262         }
00263 
00264         static inline DREAL log10(DREAL v)
00265         {
00266             return ::log(v)/::log(10.0);
00267         }
00268 
00269 #ifdef HAVE_LOG2
00270         static inline DREAL log2(DREAL v)
00271         {
00272             return ::log2(v);
00273         }
00274 #else
00275         static inline DREAL log2(DREAL v)
00276         {
00277             return ::log(v)/::log(2.0);
00278         }
00279 #endif //HAVE_LOG2
00280 
00281         static inline DREAL log(DREAL v)
00282         {
00283             return ::log(v);
00284         }
00285 
00286         template <class T> static void transpose_matrix(T*& matrix, INT& num_feat, INT& num_vec)
00287         {
00288             T* transposed=new T[num_vec*num_feat];
00289             for (INT i=0; i<num_vec; i++)
00290             {
00291                 for (INT j=0; j<num_feat; j++)
00292                     transposed[j+i*num_feat]=matrix[i+j*num_vec];
00293             }
00294 
00295             delete[] matrix;
00296             matrix=transposed;
00297 
00298             INT tmp=num_feat;
00299             num_feat=num_vec;
00300             num_vec=tmp;
00301         }
00302 
00303 #ifdef HAVE_LAPACK
00306         static DREAL* pinv(DREAL* matrix, INT rows, INT cols, DREAL* target=NULL);
00307 
00308 
00309         //C := alpha*op( A )*op( B ) + beta*C
00310         //op( X ) = X   or   op( X ) = X',
00311         static inline void dgemm(double alpha, const double* A, int rows, int cols, CBLAS_TRANSPOSE transposeA,
00312                 double *B, int cols_B, CBLAS_TRANSPOSE transposeB,
00313                 double beta, double *C)
00314         {
00315             cblas_dgemm(CblasColMajor, transposeA, transposeB, rows, cols, cols_B,
00316                     alpha, A, cols, B, cols_B, beta, C, cols);
00317         }
00318 
00319         //y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
00320         static inline void dgemv(double alpha, const double *A, int rows, int cols, const CBLAS_TRANSPOSE transposeA,
00321                 const double* X, double beta, double* Y)
00322         {
00323             cblas_dgemv(CblasColMajor, transposeA,
00324                     rows, cols, alpha, A, cols,
00325                     X, 1, beta, Y, 1);
00326         }
00327 #endif
00328 
00329         static inline LONG factorial(INT n)
00330         {
00331             LONG res=1 ;
00332             for (int i=2; i<=n; i++)
00333                 res*=i ;
00334             return res ;
00335         }
00336 
00337         static void init_random(UINT initseed=0)
00338         {
00339             if (initseed==0)
00340             {
00341                 struct timeval tv;
00342                 gettimeofday(&tv, NULL);
00343                 seed=(UINT) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
00344             }
00345             else
00346                 seed=initseed;
00347 #if !defined(CYGWIN) && !defined(__INTERIX)
00348             //seed=42;
00349             //SG_SPRINT("initializing random number generator with %d\n", seed);
00350             initstate(seed, CMath::rand_state, sizeof(CMath::rand_state));
00351 #endif
00352         }
00353 
00354         static inline LONG random()
00355         {
00356 #if defined(CYGWIN) || defined(__INTERIX)
00357             return rand();
00358 #else
00359             return ::random();
00360 #endif
00361         }
00362 
00363         static inline INT random(INT min_value, INT max_value)
00364         {
00365             INT ret = min_value + (INT) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
00366             ASSERT(ret>=min_value && ret<=max_value);
00367             return ret ;
00368         }
00369 
00370         static inline SHORTREAL random(SHORTREAL min_value, SHORTREAL max_value)
00371         {
00372             SHORTREAL ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00373 
00374             if (ret<min_value || ret>max_value)
00375                 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00376             ASSERT(ret>=min_value && ret<=max_value);
00377             return ret;
00378         }
00379 
00380         static inline DREAL random(DREAL min_value, DREAL max_value)
00381         {
00382             DREAL ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
00383 
00384             if (ret<min_value || ret>max_value)
00385                 SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
00386             ASSERT(ret>=min_value && ret<=max_value);
00387             return ret;
00388         }
00389 
00390         template <class T>
00391             static void fill_vector(T* vec, INT len, T value)
00392             {
00393                 for (INT i=0; i<len; i++)
00394                     vec[i]=value;
00395             }
00396         template <class T>
00397             static void range_fill_vector(T* vec, INT len, T start=0)
00398             {
00399                 for (INT i=0; i<len; i++)
00400                     vec[i]=i+start;
00401             }
00402 
00403         template <class T>
00404             static void random_vector(T* vec, INT len, T min_value, T max_value)
00405             {
00406                 for (INT i=0; i<len; i++)
00407                     vec[i]=CMath::random(min_value, max_value);
00408             }
00409 
00410         static inline INT* randperm(INT n)
00411         {
00412             INT* perm = new INT[n];
00413 
00414             if (!perm)
00415                 return NULL;
00416             for (INT i = 0; i < n; i++)
00417                 perm[i] = i;
00418             for (INT i = 0; i < n; i++)
00419                 swap(perm[random(0, n - 1)], perm[i]);
00420             return perm;
00421         }
00422 
00423         static inline LONG nchoosek(INT n, INT k)
00424         {
00425             long res=1 ;
00426 
00427             for (INT i=n-k+1; i<=n; i++)
00428                 res*=i ;
00429 
00430             return res/factorial(k) ;
00431         }
00432 
00434         template <class T>
00435             static inline void vec1_plus_scalar_times_vec2(T* vec1,
00436                     T scalar, const T* vec2, INT n)
00437             {
00438                 for (INT i=0; i<n; i++)
00439                     vec1[i]+=scalar*vec2[i];
00440             }
00441 
00443         static inline DREAL dot(const DREAL* v1, const DREAL* v2, INT n)
00444         {
00445             DREAL r=0;
00446 #ifdef HAVE_LAPACK
00447             INT skip=1;
00448             r = cblas_ddot(n, v1, skip, v2, skip);
00449 #else
00450             for (INT i=0; i<n; i++)
00451                 r+=v1[i]*v2[i];
00452 #endif
00453             return r;
00454         }
00455 
00457         static inline SHORTREAL dot(const SHORTREAL* v1, const SHORTREAL* v2, INT n)
00458         {
00459             DREAL r=0;
00460 #ifdef HAVE_LAPACK
00461             INT skip=1;
00462             r = cblas_sdot(n, v1, skip, v2, skip);
00463 #else
00464             for (INT i=0; i<n; i++)
00465                 r+=v1[i]*v2[i];
00466 #endif
00467             return r;
00468         }
00469 
00471         template <class T>
00472             static inline void add(T* target, T alpha, const T* v1, T beta, const T* v2, INT len)
00473             {
00474                 for (INT i=0; i<len; i++)
00475                     target[i]=alpha*v1[i]+beta*v2[i];
00476             }
00477 
00479         template <class T>
00480             static inline void add_scalar(T alpha, T* vec, INT len)
00481             {
00482                 for (INT i=0; i<len; i++)
00483                     vec[i]+=alpha;
00484             }
00485 
00487         template <class T>
00488             static inline void scale_vector(T alpha, T* vec, INT len)
00489             {
00490                 for (INT i=0; i<len; i++)
00491                     vec[i]*=alpha;
00492             }
00493 
00495         template <class T>
00496             static inline T sum(T* vec, INT len)
00497             {
00498                 T result=0;
00499                 for (INT i=0; i<len; i++)
00500                     result+=vec[i];
00501 
00502                 return result;
00503             }
00504 
00506         template <class T>
00507             static inline T max(T* vec, INT len)
00508             {
00509                 ASSERT(len>0);
00510                 T maxv=vec[0];
00511 
00512                 for (INT i=1; i<len; i++)
00513                     maxv=CMath::max(vec[i], maxv);
00514 
00515                 return maxv;
00516             }
00517 
00519         template <class T>
00520             static inline T sum_abs(T* vec, INT len)
00521             {
00522                 T result=0;
00523                 for (INT i=0; i<len; i++)
00524                     result+=CMath::abs(vec[i]);
00525 
00526                 return result;
00527             }
00528 
00529         static inline DREAL mean(DREAL* vec, INT len)
00530         {
00531             ASSERT(vec);
00532             ASSERT(len>0);
00533 
00534             DREAL mean=0;
00535             for (INT i=0; i<len; i++)
00536                 mean+=vec[i];
00537             return mean/len;
00538         }
00539 
00540         static inline DREAL trace(DREAL* mat, INT cols, INT rows)
00541         {
00542             DREAL trace=0;
00543             for (INT i=0; i<rows; i++)
00544                 trace+=mat[i*cols+i];
00545             return trace;
00546         }
00547 
00551         static void sort(INT *a, INT cols, INT sort_col=0) ;
00552         static void sort(DREAL *a, INT*idx, INT N) ;
00553 
00554         /*
00555          * Inline function to extract the byte at position p (from left)
00556          * of an 64 bit integer. The function is somewhat identical to 
00557          * accessing an array of characters via [].
00558          */
00559 
00561         template <class T>
00562             inline static void radix_sort(T* array, size_t size)
00563             {
00564                 radix_sort_helper(array,size,0);
00565             }
00566 
00567         template <class T>
00568             static inline BYTE byte(T word, unsigned short p)
00569             {
00570                 return (word >> (sizeof(T)-p-1) * 8) & 0xff;
00571             }
00572 
00573         template <class T>
00574             static void radix_sort_helper(T* array, size_t size, unsigned short i)
00575             {
00576                 static size_t count[256], nc, cmin;
00577                 T *ak;
00578                 BYTE c=0;
00579                 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
00580                 T *an, *aj, *pile[256];
00581                 size_t *cp, cmax;
00582 
00583                 /* Push initial array to stack */
00584                 sp = s;
00585                 radix_push(array, size, i);
00586 
00587                 /* Loop until all digits have been sorted */
00588                 while (sp>s) {
00589                     radix_pop(array, size, i);
00590                     an = array + size;
00591 
00592                     /* Make character histogram */
00593                     if (nc == 0) {
00594                         cmin = 0xff;
00595                         for (ak = array; ak < an; ak++) {
00596                             c = byte(*ak, i);
00597                             count[c]++;
00598                             if (count[c] == 1) {
00599                                 /* Determine smallest character */
00600                                 if (c < cmin)
00601                                     cmin = c;
00602                                 nc++;
00603                             }
00604                         }
00605 
00606                         /* Sort recursively if stack size too small */
00607                         if (sp + nc > s + RADIX_STACK_SIZE) {
00608                             radix_sort_helper(array, size, i);
00609                             continue;
00610                         }
00611                     }
00612 
00613                     /*
00614                      * Set pile[]; push incompletely sorted bins onto stack.
00615                      * pile[] = pointers to last out-of-place element in bins.
00616                      * Before permuting: pile[c-1] + count[c] = pile[c];
00617                      * during deal: pile[c] counts down to pile[c-1].
00618                      */
00619                     olds = bigs = sp;
00620                     cmax = 2;
00621                     ak = array;
00622                     pile[0xff] = an;
00623                     for (cp = count + cmin; nc > 0; cp++) {
00624                         /* Find next non-empty pile */
00625                         while (*cp == 0)
00626                             cp++;
00627                         /* Pile with several entries */
00628                         if (*cp > 1) {
00629                             /* Determine biggest pile */
00630                             if (*cp > cmax) {
00631                                 cmax = *cp;
00632                                 bigs = sp;
00633                             }
00634 
00635                             if (i < sizeof(T)-1)
00636                                 radix_push(ak, *cp, (unsigned short) (i + 1));
00637                         }
00638                         pile[cp - count] = ak += *cp;
00639                         nc--;
00640                     }
00641 
00642                     /* Play it safe -- biggest bin last. */
00643                     swap(*olds, *bigs);
00644 
00645                     /*
00646                      * Permute misplacements home. Already home: everything
00647                      * before aj, and in pile[c], items from pile[c] on.
00648                      * Inner loop:
00649                      *      r = next element to put in place;
00650                      *      ak = pile[r[i]] = location to put the next element.
00651                      *      aj = bottom of 1st disordered bin.
00652                      * Outer loop:
00653                      *      Once the 1st disordered bin is done, ie. aj >= ak,
00654                      *      aj<-aj + count[c] connects the bins in array linked list;
00655                      *      reset count[c].
00656                      */
00657                     aj = array;
00658                     while (aj<an)
00659                     {
00660                         T r;
00661 
00662                         for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
00663                             swap(*ak, r);
00664 
00665                         *aj = r;
00666                         aj += count[c];
00667                         count[c] = 0;
00668                     }
00669                 }
00670             }
00673         template <class T>
00674             static void insertion_sort(T* output, INT size)
00675             {
00676                 for (INT i=0; i<size-1; i++)
00677                 {
00678                     INT j=i-1;
00679                     T value=output[i];
00680                     while (j >= 0 && output[j] > value)
00681                     {
00682                         output[j+1] = output[j];
00683                         j--;
00684                     }
00685                     output[j+1]=value;
00686                 }
00687             }
00688 
00689 
00692         //template <class T>
00693         //static void qsort(T* output, INT size) ;
00694         template <class T>
00695             static void qsort(T* output, INT size)
00696             {
00697                 if (size==2)
00698                 {
00699                     if (output[0] > output [1])
00700                         swap(output[0],output[1]);
00701                     return;
00702                 }
00703                 T split=output[random(0,size-1)];
00704 
00705                 INT left=0;
00706                 INT right=size-1;
00707 
00708                 while (left<=right)
00709                 {
00710                     while (output[left] < split)
00711                         left++;
00712                     while (output[right] > split)
00713                         right--;
00714 
00715                     if (left<=right)
00716                     {
00717                         swap(output[left],output[right]);
00718                         left++;
00719                         right--;
00720                     }
00721                 }
00722 
00723                 if (right+1> 1)
00724                     qsort(output,right+1);
00725 
00726                 if (size-left> 1)
00727                     qsort(&output[left],size-left);
00728             }
00729 
00731         template <class T> static void display_vector(T* vector, INT n, const char* name="vector");
00732 
00734         template <class T> static void display_matrix(T* matrix, INT rows, INT cols, const char* name="matrix");
00735 
00741         template <class T1,class T2>
00742             static void qsort_index(T1* output, T2* index, UINT size);
00743 
00744         template <class T1,class T2>
00745             static void* parallel_qsort_index(void* p);
00746 
00752         template <class T1,class T2>
00753             static void qsort_backward_index(T1* output, T2* index, INT size);
00754 
00755         /* finds the smallest element in output and puts that element as the 
00756            first element  */
00757         template <class T>
00758             static void min(DREAL* output, T* index, INT size);
00759 
00760         /* finds the n smallest elements in output and puts these elements as the 
00761            first n elements  */
00762         template <class T>
00763             static void nmin(DREAL* output, T* index, INT size, INT n);
00764 
00765         /* performs a inplace unique of a vector of type T using quicksort 
00766          * returns the new number of elements */
00767         template <class T>
00768             static INT unique(T* output, INT size) 
00769             {
00770                 qsort(output, size);
00771                 INT i,j=0 ;
00772                 for (i=0; i<size; i++)
00773                     if (i==0 || output[i]!=output[i-1])
00774                         output[j++]=output[i];
00775                 return j ;
00776             }
00777 
00778         /* finds an element in a sorted array via binary search
00779          * returns -1 if not found */
00780         template <class T>
00781             static INT binary_search_helper(T* output, INT size, T elem)
00782             {
00783                 INT start=0;
00784                 INT end=size-1;
00785 
00786                 if (size<1)
00787                     return -1;
00788 
00789                 while (start<end)
00790                 {
00791                     INT middle=(start+end)/2; 
00792 
00793                     if (output[middle]>elem)
00794                         end=middle-1;
00795                     else if (output[middle]<elem)
00796                         start=middle+1;
00797                     else
00798                         return middle;
00799                 }
00800 
00801                 if (output[start]==elem)
00802                     return start;
00803                 else
00804                     return -1;
00805             }
00806 
00807         /* finds an element in a sorted array via binary search
00808          *     * returns -1 if not found */
00809         template <class T>
00810             static inline INT binary_search(T* output, INT size, T elem)
00811             {
00812                 INT ind = binary_search_helper(output, size, elem);
00813                 if (ind >= 0 && output[ind] == elem)
00814                     return ind;
00815                 return -1;
00816             }
00817 
00818         /* finds an element in a sorted array via binary search 
00819          * if it exists, else the index the largest smaller element
00820          * is returned
00821          * note: a successor is not mandatory */
00822         template <class T>
00823             static INT binary_search_max_lower_equal(T* output, INT size, T elem)
00824             {
00825                 INT ind = binary_search_helper(output, size, elem);
00826 
00827                 if (output[ind]<=elem)
00828                     return ind;
00829                 if (ind>0 && output[ind-1] <= elem)
00830                     return ind-1;
00831                 return -1;
00832             }
00833 
00836         static DREAL Align(CHAR * seq1, CHAR* seq2, INT l1, INT l2, DREAL gapCost);
00837 
00842         static INT calcroc(DREAL* fp, DREAL* tp, DREAL* output, INT* label, INT& size, INT& possize, INT& negsize, DREAL& tresh, FILE* rocfile);
00844 
00847         static double mutual_info(DREAL* p1, DREAL* p2, INT len);
00848 
00851         static double relative_entropy(DREAL* p, DREAL* q, INT len);
00852 
00854         static double entropy(DREAL* p, INT len);
00855 
00857         inline static UINT get_seed()
00858         {
00859             return CMath::seed;
00860         }
00861 
00862 
00873 #ifdef USE_LOGCACHE
00874         static inline DREAL logarithmic_sum(DREAL p, DREAL q)
00875         {
00876             DREAL diff;
00877 
00878             if (!finite(p))
00879                 return q;
00880             if (!finite(q))
00881             {
00882                 SG_SWARNING(("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
00883                         return NAN;
00884                         }
00885                         diff = p - q;
00886                         if (diff > 0)
00887                         return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
00888                         return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
00889                         }
00890 
00892                         static void init_log_table();
00893 
00895                         static INT determine_logrange();
00896 
00898                         static INT determine_logaccuracy(INT range);
00899 #else
00900                         static inline DREAL logarithmic_sum(DREAL p, DREAL q)
00901                         {
00902                         DREAL diff;
00903 
00904                         if (!finite(p))
00905                             return q;
00906                         if (!finite(q))
00907                             return p;
00908                         diff = p - q;
00909                         if (diff > 0)
00910                             return diff > LOGRANGE? p : p + log(1 + exp(-diff));
00911                         return -diff > LOGRANGE? q : q + log(1 + exp(diff));
00912                         }
00913 #endif
00914 #ifdef LOG_SUM_ARRAY
00915 
00920                 static inline DREAL logarithmic_sum_array(DREAL *p, INT len)
00921                 {
00922                     if (len<=2)
00923                     {
00924                         if (len==2)
00925                             return logarithmic_sum(p[0],p[1]) ;
00926                         if (len==1)
00927                             return p[0];
00928                         return -INFTY ;
00929                     }
00930                     else
00931                     {
00932                         register DREAL *pp=p ;
00933                         if (len%2==1) pp++ ;
00934                         for (register INT j=0; j < len>>1; j++)
00935                             pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
00936                     }
00937                     return logarithmic_sum_array(p,len%2+len>>1) ;
00938                 } 
00939 #endif
00940 
00941     public:
00944 
00945                 static const DREAL INFTY;
00946                 static const DREAL ALMOST_INFTY;
00947 
00949                 static const DREAL ALMOST_NEG_INFTY;
00950 
00952                 static INT LOGRANGE;
00953 
00955                 static UINT seed;
00956 
00957 #ifdef USE_LOGCACHE 
00958 
00960                 static INT LOGACCURACY;
00962     protected:
00964                 static DREAL* logtable; 
00965 #endif
00966                 static CHAR* rand_state;
00967 };
00968 
00969     template <class T1,class T2>
00970 void* CMath::parallel_qsort_index(void* p)
00971 {
00972     struct thread_qsort* ps=(thread_qsort*) p;
00973     T1* output=ps->output;
00974     T2* index=ps->index;
00975     INT size=ps->size;
00976     INT* qsort_threads=ps->qsort_threads;
00977     INT sort_limit=ps->sort_limit;
00978     INT num_threads=ps->num_threads;
00979 
00980     if (size==2)
00981     {
00982         if (output[0] > output [1])
00983         {
00984             swap(output[0], output[1]);
00985             swap(index[0], index[1]);
00986         }
00987         return NULL;
00988     }
00989     /*double split=output[(((uint64_t) size)*rand())/(((uint64_t)RAND_MAX)+1)];*/
00990     double split=output[size/2];
00991 
00992     INT left=0;
00993     INT right=size-1;
00994 
00995     while (left<=right)
00996     {
00997         while (output[left] < split)
00998             left++;
00999         while (output[right] > split)
01000             right--;
01001 
01002         if (left<=right)
01003         {
01004             swap(output[left], output[right]);
01005             swap(index[left], index[right]);
01006             left++;
01007             right--;
01008         }
01009     }
01010     bool lthread_start=false;
01011     bool rthread_start=false;
01012     pthread_t lthread;
01013     pthread_t rthread;
01014     struct thread_qsort t1;
01015     struct thread_qsort t2;
01016 
01017     if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
01018         qsort_index(output,index,right+1);
01019     else if (right+1> 1)
01020     {
01021         *qsort_threads++;
01022         lthread_start=true;
01023         t1.output=output;
01024         t1.index=index;
01025         t1.size=right+1;
01026         if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, (void*) &t1) != 0)
01027         {
01028             lthread_start=false;
01029             *qsort_threads--;
01030             qsort_index(output,index,right+1);
01031         }
01032     }
01033 
01034 
01035     if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
01036         qsort_index(&output[left],&index[left], size-left);
01037     else if (size-left> 1)
01038     {
01039         *qsort_threads++;
01040         rthread_start=true;
01041         t2.output=&output[left];
01042         t2.index=&index[left];
01043         t2.size=size-left;
01044         if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, (void*)&t2) != 0)
01045         {
01046             rthread_start=false;
01047             *qsort_threads--;
01048             qsort_index(&output[left],&index[left], size-left);
01049         }
01050     }
01051 
01052     if (lthread_start)
01053     {
01054         pthread_join(lthread, NULL);
01055         *qsort_threads--;
01056     }
01057 
01058     if (rthread_start)
01059     {
01060         pthread_join(rthread, NULL);
01061         *qsort_threads--;
01062     }
01063 
01064     return NULL;
01065 }
01066 
01067 
01068 //implementations of template functions
01069     template <class T1,class T2>
01070 void CMath::qsort_index(T1* output, T2* index, UINT size)
01071 {
01072     if (size==2)
01073     {
01074         if (output[0] > output [1])
01075         {
01076             swap(output[0],output[1]);
01077             swap(index[0],index[1]);
01078         }
01079         return;
01080     }
01081     T1 split=output[random(0,size-1)];
01082 
01083     INT left=0;
01084     INT right=size-1;
01085 
01086     while (left<=right)
01087     {
01088         while (output[left] < split)
01089             left++;
01090         while (output[right] > split)
01091             right--;
01092 
01093         if (left<=right)
01094         {
01095             swap(output[left],output[right]);
01096             swap(index[left],index[right]);
01097             left++;
01098             right--;
01099         }
01100     }
01101 
01102     if (right+1> 1)
01103         qsort_index(output,index,right+1);
01104 
01105     if (size-left> 1)
01106         qsort_index(&output[left],&index[left], size-left);
01107 }
01108 
01109     template <class T1,class T2>
01110 void CMath::qsort_backward_index(T1* output, T2* index, INT size)
01111 {
01112     if (size==2)
01113     {
01114         if (output[0] < output [1])
01115         {
01116             swap(output[0],output[1]);
01117             swap(index[0],index[1]);
01118         }
01119         return;
01120     }
01121 
01122     T1 split=output[random(0,size-1)];
01123 
01124     INT left=0;
01125     INT right=size-1;
01126 
01127     while (left<=right)
01128     {
01129         while (output[left] > split)
01130             left++;
01131         while (output[right] < split)
01132             right--;
01133 
01134         if (left<=right)
01135         {
01136             swap(output[left],output[right]);
01137             swap(index[left],index[right]);
01138             left++;
01139             right--;
01140         }
01141     }
01142 
01143     if (right+1> 1)
01144         qsort_backward_index(output,index,right+1);
01145 
01146     if (size-left> 1)
01147         qsort_backward_index(&output[left],&index[left], size-left);
01148 }
01149 
01150     template <class T> 
01151 void CMath::nmin(DREAL* output, T* index, INT size, INT n)
01152 {
01153     if (6*n*size<13*size*CMath::log(size))
01154         for (INT i=0; i<n; i++)
01155             min(&output[i], &index[i], size-i) ;
01156     else
01157         qsort_index(output, index, size) ;
01158 }
01159 
01160 /* move the smallest entry in the array to the beginning */
01161     template <class T>
01162 void CMath::min(DREAL* output, T* index, INT size)
01163 {
01164     if (size<=1)
01165         return;
01166     DREAL min_elem=output[0];
01167     INT min_index=0;
01168     for (INT i=1; i<size; i++)
01169     {
01170         if (output[i]<min_elem)
01171         {
01172             min_index=i;
01173             min_elem=output[i];
01174         }
01175     }
01176     swap(output[0], output[min_index]);
01177     swap(index[0], index[min_index]);
01178 }
01179 #endif

SHOGUN Machine Learning Toolbox - Documentation