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

SHOGUN Machine Learning Toolbox - Documentation