SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Math.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) 2013 Thoralf Klein
00008  * Written (W) 2013 Soumyajit De
00009  * Written (W) 2012 Fernando Jose Iglesias Garcia
00010  * Written (W) 2011 Siddharth Kherada
00011  * Written (W) 2011 Justin Patera
00012  * Written (W) 2011 Alesis Novik
00013  * Written (W) 2011-2013 Heiko Strathmann
00014  * Written (W) 1999-2009 Soeren Sonnenburg
00015  * Written (W) 1999-2008 Gunnar Raetsch
00016  * Written (W) 2007 Konrad Rieck
00017  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
00018  */
00019 
00020 #ifndef __MATHEMATICS_H_
00021 #define __MATHEMATICS_H_
00022 
00023 #include <shogun/base/SGObject.h>
00024 #include <shogun/lib/common.h>
00025 #include <shogun/io/SGIO.h>
00026 #include <shogun/base/Parallel.h>
00027 #include <shogun/mathematics/Random.h>
00028 
00029 #ifndef _USE_MATH_DEFINES
00030 #define _USE_MATH_DEFINES
00031 #endif
00032 
00033 #include <math.h>
00034 #include <stdio.h>
00035 #include <float.h>
00036 #include <sys/types.h>
00037 #ifndef _WIN32
00038 #include <unistd.h>
00039 #endif
00040 
00041 #ifdef HAVE_PTHREAD
00042 #include <pthread.h>
00043 #endif
00044 
00045 #ifdef SUNOS
00046 #include <ieeefp.h>
00047 #endif
00048 
00050 #ifdef log2
00051 #define cygwin_log2 log2
00052 #undef log2
00053 #endif
00054 
00055 #ifndef M_PI
00056 #define M_PI 3.14159265358979323846
00057 #endif
00058 
00059 #ifdef _WIN32
00060 #ifndef isnan
00061 #define isnan _isnan
00062 #endif
00063 
00064 #ifndef isfinite
00065 #define isfinite _isfinite
00066 #endif
00067 #endif //_WIN32
00068 
00069 /* Size of RNG seed */
00070 #define RNG_SEED_SIZE 256
00071 
00072 /* Maximum stack size */
00073 #define RADIX_STACK_SIZE        512
00074 
00075 /* Stack macros */
00076 #define radix_push(a, n, i)     sp->sa = a, sp->sn = n, (sp++)->si = i
00077 #define radix_pop(a, n, i)      a = (--sp)->sa, n = sp->sn, i = sp->si
00078 
00079 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00080 
00081 template <class T> struct radix_stack_t
00082 {
00084     T *sa;
00086     size_t sn;
00088     uint16_t si;
00089 };
00090 
00092 
00094 template <class T1, class T2> struct thread_qsort
00095 {
00097     T1* output;
00099     T2* index;
00101     uint32_t size;
00102 
00104     int32_t* qsort_threads;
00106     int32_t sort_limit;
00108     int32_t num_threads;
00109 };
00110 #endif // DOXYGEN_SHOULD_SKIP_THIS
00111 
00112 #define COMPLEX128_ERROR_ONEARG(function)   \
00113 static inline complex128_t function(complex128_t a) \
00114 {   \
00115     SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
00116         #function);\
00117     return complex128_t(0.0, 0.0);  \
00118 }
00119 
00120 #define COMPLEX128_STDMATH(function)    \
00121 static inline complex128_t function(complex128_t a) \
00122 {   \
00123     return std::function(a);    \
00124 }
00125 
00126 namespace shogun
00127 {
00129     extern CRandom* sg_rand;
00130     class CSGObject;
00133 class CMath : public CSGObject
00134 {
00135     public:
00139 
00140         CMath();
00141 
00143         virtual ~CMath();
00145 
00149 
00151         //
00152         template <class T>
00153             static inline T min(T a, T b)
00154             {
00155                 return (a<=b) ? a : b;
00156             }
00157 
00159         template <class T>
00160             static inline T max(T a, T b)
00161             {
00162                 return (a>=b) ? a : b;
00163             }
00164 
00166         template <class T>
00167             static inline T clamp(T value, T lb, T ub)
00168             {
00169                 if (value<=lb)
00170                     return lb;
00171                 else if (value>=ub)
00172                     return ub;
00173                 else
00174                     return value;
00175             }
00176 
00178         template <class T>
00179             static inline T abs(T a)
00180             {
00181                 // can't be a>=0?(a):(-a), because compiler complains about
00182                 // 'comparison always true' when T is unsigned
00183                 if (a==0)
00184                     return 0;
00185                 else if (a>0)
00186                     return a;
00187                 else
00188                     return -a;
00189             }
00190 
00192         static inline float64_t abs(complex128_t a)
00193         {
00194             float64_t a_real=a.real();
00195             float64_t a_imag=a.imag();
00196             return (CMath::sqrt(a_real*a_real+a_imag*a_imag));
00197         }
00199 
00202 
00203         static inline float64_t round(float64_t d)
00204         {
00205             return ::floor(d+0.5);
00206         }
00207 
00208         static inline float64_t floor(float64_t d)
00209         {
00210             return ::floor(d);
00211         }
00212 
00213         static inline float64_t ceil(float64_t d)
00214         {
00215             return ::ceil(d);
00216         }
00217 
00219         template <class T>
00220             static inline T sign(T a)
00221             {
00222                 if (a==0)
00223                     return 0;
00224                 else return (a<0) ? (-1) : (+1);
00225             }
00226 
00228         template <class T>
00229             static inline void swap(T &a,T &b)
00230             {
00231             /* register for fast swaps */
00232                 register T c=a;
00233                 a=b;
00234                 b=c;
00235             }
00236 
00238         template <class T>
00239             static inline T sq(T x)
00240             {
00241                 return x*x;
00242             }
00243 
00245         static inline float32_t sqrt(float32_t x)
00246         {
00247             return ::sqrtf(x);
00248         }
00249 
00251         static inline float64_t sqrt(float64_t x)
00252         {
00253             return ::sqrt(x);
00254         }
00255 
00257         static inline floatmax_t sqrt(floatmax_t x)
00258         {
00259             //fall back to double precision sqrt if sqrtl is not
00260             //available
00261 #ifdef HAVE_SQRTL
00262             return ::sqrtl(x);
00263 #else
00264             return ::sqrt(x);
00265 #endif
00266         }
00267 
00269         COMPLEX128_STDMATH(sqrt)
00270 
00271         
00272         static inline float32_t invsqrt(float32_t x)
00273         {
00274             union float_to_int
00275             {
00276                 float32_t f;
00277                 int32_t i;
00278             };
00279 
00280             float_to_int tmp;
00281             tmp.f=x;
00282 
00283             float32_t xhalf = 0.5f * x;
00284             // store floating-point bits in integer tmp.i
00285             // and do initial guess for Newton's method
00286             tmp.i = 0x5f3759d5 - (tmp.i >> 1);
00287             x = tmp.f; // convert new bits into float
00288             x = x*(1.5f - xhalf*x*x); // One round of Newton's method
00289             return x;
00290         }
00291 
00293         static inline floatmax_t powl(floatmax_t x, floatmax_t n)
00294         {
00295             //fall back to double precision pow if powl is not
00296             //available
00297 #ifdef HAVE_POWL
00298             return ::powl((long double) x, (long double) n);
00299 #else
00300             return ::pow((double) x, (double) n);
00301 #endif
00302         }
00303 
00304         static inline int32_t pow(bool x, int32_t n)
00305         {
00306             return 0;
00307         }
00308 
00309         static inline int32_t pow(int32_t x, int32_t n)
00310         {
00311             ASSERT(n>=0)
00312             int32_t result=1;
00313             while (n--)
00314                 result*=x;
00315 
00316             return result;
00317         }
00318 
00319         static inline float64_t pow(float64_t x, int32_t n)
00320         {
00321             if (n>=0)
00322             {
00323                 float64_t result=1;
00324                 while (n--)
00325                     result*=x;
00326 
00327                 return result;
00328             }
00329             else
00330                 return ::pow((double)x, (double)n);
00331         }
00332 
00333         static inline float64_t pow(float64_t x, float64_t n)
00334         {
00335             return ::pow((double) x, (double) n);
00336         }
00337 
00339         static inline complex128_t pow(complex128_t x, int32_t n)
00340         {
00341             return std::pow(x, n);
00342         }
00343 
00344         static inline complex128_t pow(complex128_t x, complex128_t n)
00345         {
00346             return std::pow(x, n);
00347         }
00348 
00349         static inline complex128_t pow(complex128_t x, float64_t n)
00350         {
00351             return std::pow(x, n);
00352         }
00353 
00354         static inline complex128_t pow(float64_t x, complex128_t n)
00355         {
00356             return std::pow(x, n);
00357         }
00358 
00359         static inline float64_t exp(float64_t x)
00360         {
00361             return ::exp((double) x);
00362         }
00363 
00365         COMPLEX128_STDMATH(exp)
00366 
00367         
00368         static inline float64_t tan(float64_t x)
00369         {
00370             return ::tan((double) x);
00371         }
00372 
00374         COMPLEX128_STDMATH(tan)
00375 
00376         
00377         static inline float64_t atan(float64_t x)
00378         {
00379             return ::atan((double) x);
00380         }
00381 
00383         COMPLEX128_ERROR_ONEARG(atan)
00384 
00385         
00386         static inline float64_t atan2(float64_t x, float64_t y)
00387         {
00388             return ::atan2((double) x, (double) y);
00389         }
00390 
00392         COMPLEX128_ERROR_ONEARG(atan2)
00393 
00394         
00395         static inline float64_t tanh(float64_t x)
00396         {
00397             return ::tanh((double) x);
00398         }
00399 
00401         COMPLEX128_STDMATH(tanh)
00402 
00403         static inline float64_t log10(float64_t v)
00404         {
00405             return ::log(v)/::log(10.0);
00406         }
00407 
00409         COMPLEX128_STDMATH(log10)
00410 
00411         static inline float64_t log2(float64_t v)
00412         {
00413 #ifdef HAVE_LOG2
00414             return ::log2(v);
00415 #else
00416             return ::log(v)/::log(2.0);
00417 #endif //HAVE_LOG2
00418         }
00419 
00420         static inline float64_t log(float64_t v)
00421         {
00422             return ::log(v);
00423         }
00424 
00426         COMPLEX128_STDMATH(log)
00427 
00428         static inline index_t floor_log(index_t n)
00429         {
00430             index_t i;
00431             for (i = 0; n != 0; i++)
00432                 n >>= 1;
00433 
00434             return i;
00435         }
00436 
00437         static inline float64_t sin(float64_t x)
00438         {
00439             return ::sin(x);
00440         }
00441 
00443         COMPLEX128_STDMATH(sin)
00444 
00445         static inline float64_t asin(float64_t x)
00446         {
00447             return ::asin(x);
00448         }
00449 
00451         COMPLEX128_ERROR_ONEARG(asin)
00452 
00453         static inline float64_t sinh(float64_t x)
00454         {
00455             return ::sinh(x);
00456         }
00457 
00459         COMPLEX128_STDMATH(sinh)
00460 
00461         static inline float64_t cos(float64_t x)
00462         {
00463             return ::cos(x);
00464         }
00465 
00467         COMPLEX128_STDMATH(cos)
00468 
00469         static inline float64_t acos(float64_t x)
00470         {
00471             return ::acos(x);
00472         }
00473 
00475         COMPLEX128_ERROR_ONEARG(acos)
00476 
00477         static inline float64_t cosh(float64_t x)
00478         {
00479             return ::cosh(x);
00480         }
00481 
00483         COMPLEX128_STDMATH(cosh)
00484 
00485         static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
00486         {
00487             ASSERT(len>0 && xy)
00488 
00489             float64_t area = 0.0;
00490 
00491             if (!reversed)
00492             {
00493                 for (int i=1; i<len; i++)
00494                     area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
00495             }
00496             else
00497             {
00498                 for (int i=1; i<len; i++)
00499                     area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
00500             }
00501 
00502             return area;
00503         }
00504 
00505         static bool strtof(const char* str, float32_t* float_result);
00506         static bool strtod(const char* str, float64_t* double_result);
00507         static bool strtold(const char* str, floatmax_t* long_double_result);
00508 
00509         static inline int64_t factorial(int32_t n)
00510         {
00511             int64_t res=1;
00512             for (int i=2; i<=n; i++)
00513                 res*=i ;
00514             return res ;
00515         }
00516 
00517         static void init_random(uint32_t initseed=0)
00518         {
00519             if (initseed==0)
00520                 seed = CRandom::generate_seed();
00521             else
00522                 seed=initseed;
00523 
00524             sg_rand->set_seed(seed);
00525         }
00526 
00527         static inline uint64_t random()
00528         {
00529             return sg_rand->random_64();
00530         }
00531 
00532         static inline uint64_t random(uint64_t min_value, uint64_t max_value)
00533         {
00534             return sg_rand->random(min_value, max_value);
00535         }
00536 
00537         static inline int64_t random(int64_t min_value, int64_t max_value)
00538         {
00539             return sg_rand->random(min_value, max_value);
00540         }
00541 
00542         static inline uint32_t random(uint32_t min_value, uint32_t max_value)
00543         {
00544             return sg_rand->random(min_value, max_value);
00545         }
00546 
00547         static inline int32_t random(int32_t min_value, int32_t max_value)
00548         {
00549             return sg_rand->random(min_value, max_value);
00550         }
00551 
00552         static inline float32_t random(float32_t min_value, float32_t max_value)
00553         {
00554             return sg_rand->random(min_value, max_value);
00555         }
00556 
00557         static inline float64_t random(float64_t min_value, float64_t max_value)
00558         {
00559             return sg_rand->random(min_value, max_value);
00560         }
00561 
00562         static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
00563         {
00564             return sg_rand->random(min_value, max_value);
00565         }
00566 
00570         static inline float32_t normal_random(float32_t mean, float32_t std_dev)
00571         {
00572             // sets up variables & makes sure rand_s.range == (0,1)
00573             float32_t ret;
00574             float32_t rand_u;
00575             float32_t rand_v;
00576             float32_t rand_s;
00577             do
00578             {
00579                 rand_u = CMath::random(-1.0, 1.0);
00580                 rand_v = CMath::random(-1.0, 1.0);
00581                 rand_s = rand_u*rand_u + rand_v*rand_v;
00582             } while ((rand_s == 0) || (rand_s >= 1));
00583 
00584             // the meat & potatos, and then the mean & standard deviation shifting...
00585             ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
00586             ret = std_dev*ret + mean;
00587             return ret;
00588         }
00589 
00593         static inline float64_t normal_random(float64_t mean, float64_t std_dev)
00594         {
00595             return sg_rand->normal_distrib(mean, std_dev);
00596         }
00597 
00600         static inline float32_t randn_float()
00601         {
00602             return normal_random(0.0, 1.0);
00603         }
00604 
00607         static inline float64_t randn_double()
00608         {
00609             return sg_rand->std_normal_distrib();
00610         }
00611 
00612         template <class T>
00613             static int32_t get_num_nonzero(T* vec, int32_t len)
00614             {
00615                 int32_t nnz = 0;
00616                 for (index_t i=0; i<len; ++i)
00617                     nnz += vec[i] != 0;
00618 
00619                 return nnz;
00620             }
00621 
00622         static int32_t get_num_nonzero(complex128_t* vec, int32_t len)
00623         {
00624                 int32_t nnz=0;
00625                 for (index_t i=0; i<len; ++i)
00626                     nnz+=vec[i]!=0.0;
00627                 return nnz;
00628         }
00629 
00630         static inline int64_t nchoosek(int32_t n, int32_t k)
00631         {
00632             int64_t res=1;
00633 
00634             for (int32_t i=n-k+1; i<=n; i++)
00635                 res*=i;
00636 
00637             return res/factorial(k);
00638         }
00639 
00647         static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100);
00648 
00655         template <class T>
00656         static T log_sum_exp(SGVector<T> values)
00657         {
00658             REQUIRE(values.vector, "Values are empty");
00659 
00660             /* find minimum element index */
00661             index_t min_index=0;
00662             T X0=values[0];
00663             if (values.vlen>1)
00664             {
00665                 for (index_t i=1; i<values.vlen; ++i)
00666                 {
00667                     if (values[i]<X0)
00668                     {
00669                         X0=values[i];
00670                         min_index=i;
00671                     }
00672                 }
00673             }
00674 
00675             /* remove element from vector copy and compute log sum exp */
00676             SGVector<T> values_without_X0(values.vlen-1);
00677             index_t from_idx=0;
00678             index_t to_idx=0;
00679             for (from_idx=0; from_idx<values.vlen; ++from_idx)
00680             {
00681                 if (from_idx!=min_index)
00682                 {
00683                     values_without_X0[to_idx]=exp(values[from_idx]-X0);
00684                     to_idx++;
00685                 }
00686             }
00687 
00688             return X0+log(SGVector<T>::sum(values_without_X0)+1);
00689         }
00690 
00697         template <class T>
00698         static T log_mean_exp(SGVector<T> values)
00699         {
00700             return log_sum_exp(values) - log(values.vlen);
00701         }
00702 
00706         static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
00707         static void sort(float64_t *a, int32_t*idx, int32_t N);
00708 
00711         template <class T>
00712             static void qsort(T* output, int32_t size)
00713             {
00714                 if (size<=1)
00715                     return;
00716 
00717                 if (size==2)
00718                 {
00719                     if (output[0] > output [1])
00720                         CMath::swap(output[0],output[1]);
00721                     return;
00722                 }
00723                 //T split=output[random(0,size-1)];
00724                 T split=output[size/2];
00725 
00726                 int32_t left=0;
00727                 int32_t right=size-1;
00728 
00729                 while (left<=right)
00730                 {
00731                     while (output[left] < split)
00732                         left++;
00733                     while (output[right] > split)
00734                         right--;
00735 
00736                     if (left<=right)
00737                     {
00738                         CMath::swap(output[left],output[right]);
00739                         left++;
00740                         right--;
00741                     }
00742                 }
00743 
00744                 if (right+1> 1)
00745                     qsort(output,right+1);
00746 
00747                 if (size-left> 1)
00748                     qsort(&output[left],size-left);
00749             }
00750 
00753         template <class T>
00754             static void insertion_sort(T* output, int32_t size)
00755             {
00756                 for (int32_t i=0; i<size-1; i++)
00757                 {
00758                     int32_t j=i-1;
00759                     T value=output[i];
00760                     while (j >= 0 && output[j] > value)
00761                     {
00762                         output[j+1] = output[j];
00763                         j--;
00764                     }
00765                     output[j+1]=value;
00766                 }
00767             }
00768 
00770         template <class T>
00771             inline static void radix_sort(T* array, int32_t size)
00772             {
00773                 radix_sort_helper(array,size,0);
00774             }
00775 
00776         /*
00777          * Inline function to extract the byte at position p (from left)
00778          * of an 64 bit integer. The function is somewhat identical to
00779          * accessing an array of characters via [].
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 
00788         static inline uint8_t byte(complex128_t word, uint16_t p)
00789         {
00790             SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
00791             return uint8_t(0);
00792         }
00793 
00794         template <class T>
00795             static void radix_sort_helper(T* array, int32_t size, uint16_t i)
00796             {
00797                 static size_t count[256], nc, cmin;
00798                 T *ak;
00799                 uint8_t c=0;
00800                 radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
00801                 T *an, *aj, *pile[256];
00802                 size_t *cp, cmax;
00803 
00804                 /* Push initial array to stack */
00805                 sp = s;
00806                 radix_push(array, size, i);
00807 
00808                 /* Loop until all digits have been sorted */
00809                 while (sp>s) {
00810                     radix_pop(array, size, i);
00811                     an = array + size;
00812 
00813                     /* Make character histogram */
00814                     if (nc == 0) {
00815                         cmin = 0xff;
00816                         for (ak = array; ak < an; ak++) {
00817                             c = byte(*ak, i);
00818                             count[c]++;
00819                             if (count[c] == 1) {
00820                                 /* Determine smallest character */
00821                                 if (c < cmin)
00822                                     cmin = c;
00823                                 nc++;
00824                             }
00825                         }
00826 
00827                         /* Sort recursively if stack size too small */
00828                         if (sp + nc > s + RADIX_STACK_SIZE) {
00829                             radix_sort_helper(array, size, i);
00830                             continue;
00831                         }
00832                     }
00833 
00834                     /*
00835                      * Set pile[]; push incompletely sorted bins onto stack.
00836                      * pile[] = pointers to last out-of-place element in bins.
00837                      * Before permuting: pile[c-1] + count[c] = pile[c];
00838                      * during deal: pile[c] counts down to pile[c-1].
00839                      */
00840                     olds = bigs = sp;
00841                     cmax = 2;
00842                     ak = array;
00843                     pile[0xff] = an;
00844                     for (cp = count + cmin; nc > 0; cp++) {
00845                         /* Find next non-empty pile */
00846                         while (*cp == 0)
00847                             cp++;
00848                         /* Pile with several entries */
00849                         if (*cp > 1) {
00850                             /* Determine biggest pile */
00851                             if (*cp > cmax) {
00852                                 cmax = *cp;
00853                                 bigs = sp;
00854                             }
00855 
00856                             if (i < sizeof(T)-1)
00857                                 radix_push(ak, *cp, (uint16_t) (i + 1));
00858                         }
00859                         pile[cp - count] = ak += *cp;
00860                         nc--;
00861                     }
00862 
00863                     /* Play it safe -- biggest bin last. */
00864                     swap(*olds, *bigs);
00865 
00866                     /*
00867                      * Permute misplacements home. Already home: everything
00868                      * before aj, and in pile[c], items from pile[c] on.
00869                      * Inner loop:
00870                      *      r = next element to put in place;
00871                      *      ak = pile[r[i]] = location to put the next element.
00872                      *      aj = bottom of 1st disordered bin.
00873                      * Outer loop:
00874                      *      Once the 1st disordered bin is done, ie. aj >= ak,
00875                      *      aj<-aj + count[c] connects the bins in array linked list;
00876                      *      reset count[c].
00877                      */
00878                     aj = array;
00879                     while (aj<an)
00880                     {
00881                         T r;
00882 
00883                         for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
00884                             swap(*ak, r);
00885 
00886                         *aj = r;
00887                         aj += count[c];
00888                         count[c] = 0;
00889                     }
00890                 }
00891             }
00892 
00894         static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
00895         {
00896             SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
00897         }
00898 
00908         template <class T>
00909             static void qsort(T** vector, index_t length)
00910             {
00911                 if (length<=1)
00912                     return;
00913 
00914                 if (length==2)
00915                 {
00916                     if (*vector[0]>*vector[1])
00917                         swap(vector[0],vector[1]);
00918                     return;
00919                 }
00920                 T* split=vector[length/2];
00921 
00922                 int32_t left=0;
00923                 int32_t right=length-1;
00924 
00925                 while (left<=right)
00926                 {
00927                     while (*vector[left]<*split)
00928                         ++left;
00929                     while (*vector[right]>*split)
00930                         --right;
00931 
00932                     if (left<=right)
00933                     {
00934                         swap(vector[left],vector[right]);
00935                         ++left;
00936                         --right;
00937                     }
00938                 }
00939 
00940                 if (right+1>1)
00941                     qsort(vector,right+1);
00942 
00943                 if (length-left>1)
00944                     qsort(&vector[left],length-left);
00945             }
00946 
00948         static void qsort(complex128_t** vector, index_t length)
00949         {
00950             SG_SERROR("CMath::qsort():: Not supported for complex128_t\n");
00951         }
00952 
00954         template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
00955         {
00956             ASSERT(width>=0)
00957             for (int i=0; i<width; i++)
00958             {
00959                 T mask = ((T) 1)<<(sizeof(T)*8-1);
00960                 while (mask)
00961                 {
00962                     if (mask & word)
00963                         SG_SPRINT("1")
00964                     else
00965                         SG_SPRINT("0")
00966 
00967                     mask>>=1;
00968                 }
00969             }
00970         }
00971 
00973         static void display_bits(complex128_t word,
00974             int32_t width=8*sizeof(complex128_t))
00975         {
00976             SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
00977         }
00978 
00984         template <class T1,class T2>
00985             static void qsort_index(T1* output, T2* index, uint32_t size);
00986 
00988         template <class T>
00989             static void qsort_index(complex128_t* output, T* index, uint32_t size)
00990             {
00991                 SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
00992             }
00993 
00999         template <class T1,class T2>
01000             static void qsort_backward_index(
01001                 T1* output, T2* index, int32_t size);
01002 
01004         template <class T>
01005             static void qsort_backword_index(
01006                 complex128_t* output, T* index, uint32_t size)
01007             {
01008                 SG_SERROR("CMath::qsort_backword_index():: \
01009                     Not supported for complex128_t\n");
01010             }
01011 
01019         template <class T1,class T2>
01020             inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
01021             {
01022                 int32_t n=0;
01023                 thread_qsort<T1,T2> t;
01024                 t.output=output;
01025                 t.index=index;
01026                 t.size=size;
01027                 t.qsort_threads=&n;
01028                 t.sort_limit=limit;
01029                 t.num_threads=n_threads;
01030                 parallel_qsort_index<T1,T2>(&t);
01031             }
01032 
01034         template <class T>
01035             inline static void parallel_qsort_index(complex128_t* output, T* index,
01036                 uint32_t size, int32_t n_threads, int32_t limit=0)
01037             {
01038                 SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
01039             }
01040 
01041 
01042         template <class T1,class T2>
01043             static void* parallel_qsort_index(void* p);
01044 
01045 
01046         /* finds the smallest element in output and puts that element as the
01047            first element  */
01048         template <class T>
01049             static void min(float64_t* output, T* index, int32_t size);
01050 
01052         static void min(float64_t* output, complex128_t* index, int32_t size)
01053         {
01054             SG_SERROR("CMath::min():: Not supported for complex128_t\n");
01055         }
01056 
01057         /* finds the n smallest elements in output and puts these elements as the
01058            first n elements  */
01059         template <class T>
01060             static void nmin(
01061                 float64_t* output, T* index, int32_t size, int32_t n);
01062 
01064         static void nmin(float64_t* output, complex128_t* index,
01065             int32_t size, int32_t n)
01066         {
01067             SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
01068         }
01069 
01070 
01071 
01072         /* finds an element in a sorted array via binary search
01073          * returns -1 if not found */
01074         template <class T>
01075             static int32_t binary_search_helper(T* output, int32_t size, T elem)
01076             {
01077                 int32_t start=0;
01078                 int32_t end=size-1;
01079 
01080                 if (size<1)
01081                     return 0;
01082 
01083                 while (start<end)
01084                 {
01085                     int32_t middle=(start+end)/2;
01086 
01087                     if (output[middle]>elem)
01088                         end=middle-1;
01089                     else if (output[middle]<elem)
01090                         start=middle+1;
01091                     else
01092                         return middle;
01093                 }
01094 
01095                 return start;
01096             }
01097 
01099         static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
01100         {
01101             SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
01102             return int32_t(0);
01103         }
01104 
01105         /* finds an element in a sorted array via binary search
01106          *     * returns -1 if not found */
01107         template <class T>
01108             static inline int32_t binary_search(T* output, int32_t size, T elem)
01109             {
01110                 int32_t ind = binary_search_helper(output, size, elem);
01111                 if (ind >= 0 && output[ind] == elem)
01112                     return ind;
01113                 return -1;
01114             }
01115 
01117         static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
01118         {
01119             SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
01120             return int32_t(-1);
01121         }
01122 
01123 
01124         /* Finds an element in a sorted array of pointers via binary search
01125          * Every element is dereferenced once before being compared
01126          *
01127          * @param array array of pointers to search in (assumed being sorted)
01128          * @param length length of array
01129          * @param elem pointer to element to search for
01130          * @return index of elem, -1 if not found */
01131         template<class T>
01132             static inline int32_t binary_search(T** vector, index_t length,
01133                     T* elem)
01134             {
01135                 int32_t start=0;
01136                 int32_t end=length-1;
01137 
01138                 if (length<1)
01139                     return -1;
01140 
01141                 while (start<end)
01142                 {
01143                     int32_t middle=(start+end)/2;
01144 
01145                     if (*vector[middle]>*elem)
01146                         end=middle-1;
01147                     else if (*vector[middle]<*elem)
01148                         start=middle+1;
01149                     else
01150                     {
01151                         start=middle;
01152                         break;
01153                     }
01154                 }
01155 
01156                 if (start>=0&&*vector[start]==*elem)
01157                     return start;
01158 
01159                 return -1;
01160             }
01161 
01163         static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
01164         {
01165             SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
01166             return int32_t(-1);
01167         }
01168 
01169         template <class T>
01170             static int32_t binary_search_max_lower_equal(
01171                 T* output, int32_t size, T elem)
01172             {
01173                 int32_t ind = binary_search_helper(output, size, elem);
01174 
01175                 if (output[ind]<=elem)
01176                     return ind;
01177                 if (ind>0 && output[ind-1] <= elem)
01178                     return ind-1;
01179                 return -1;
01180             }
01181 
01183         static int32_t binary_search_max_lower_equal(complex128_t* output,
01184             int32_t size, complex128_t elem)
01185         {
01186             SG_SERROR("CMath::binary_search_max_lower_equal():: \
01187                 Not supported for complex128_t\n");
01188             return int32_t(-1);
01189         }
01190 
01193         static float64_t Align(
01194             char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
01195 
01196 
01198 
01200         static float64_t real(complex128_t c)
01201         {
01202             return c.real();
01203         }
01204 
01206         static float64_t imag(complex128_t c)
01207         {
01208             return c.imag();
01209         }
01210 
01212         inline static uint32_t get_seed()
01213         {
01214             return CMath::seed;
01215         }
01216 
01218         inline static uint32_t get_log_range()
01219         {
01220             return CMath::LOGRANGE;
01221         }
01222 
01223 #ifdef USE_LOGCACHE
01224 
01225         inline static uint32_t get_log_accuracy()
01226         {
01227             return CMath::LOGACCURACY;
01228         }
01229 #endif
01230 
01232         static int is_finite(double f);
01233 
01235         static int is_infinity(double f);
01236 
01238         static int is_nan(double f);
01239 
01250 #ifdef USE_LOGCACHE
01251         static inline float64_t logarithmic_sum(float64_t p, float64_t q)
01252         {
01253             float64_t diff;
01254 
01255             if (!CMath::is_finite(p))
01256                 return q;
01257 
01258             if (!CMath::is_finite(q))
01259             {
01260                 SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
01261                 return NOT_A_NUMBER;
01262             }
01263             diff = p - q;
01264             if (diff > 0)
01265                 return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
01266             return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
01267         }
01268 
01270         static void init_log_table();
01271 
01273         static int32_t determine_logrange();
01274 
01276         static int32_t determine_logaccuracy(int32_t range);
01277 #else
01278         static inline float64_t logarithmic_sum(
01279                 float64_t p, float64_t q)
01280         {
01281             float64_t diff;
01282 
01283             if (!CMath::is_finite(p))
01284                 return q;
01285             if (!CMath::is_finite(q))
01286                 return p;
01287             diff = p - q;
01288             if (diff > 0)
01289                 return diff > LOGRANGE? p : p + log(1 + exp(-diff));
01290             return -diff > LOGRANGE? q : q + log(1 + exp(diff));
01291         }
01292 #endif
01293 #ifdef USE_LOGSUMARRAY
01294 
01299                 static inline float64_t logarithmic_sum_array(
01300                     float64_t *p, int32_t len)
01301                 {
01302                     if (len<=2)
01303                     {
01304                         if (len==2)
01305                             return logarithmic_sum(p[0],p[1]) ;
01306                         if (len==1)
01307                             return p[0];
01308                         return -INFTY ;
01309                     }
01310                     else
01311                     {
01312                         register float64_t *pp=p ;
01313                         if (len%2==1) pp++ ;
01314                         for (register int32_t j=0; j < len>>1; j++)
01315                             pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
01316                     }
01317                     return logarithmic_sum_array(p,len%2 + (len>>1)) ;
01318                 }
01319 #endif //USE_LOGSUMARRAY
01320 
01321 
01323                 virtual const char* get_name() const { return "Math"; }
01324     public:
01327 
01328                 static const float64_t NOT_A_NUMBER;
01330                 static const float64_t INFTY;
01331                 static const float64_t ALMOST_INFTY;
01332 
01334                 static const float64_t ALMOST_NEG_INFTY;
01335 
01337                 static const float64_t PI;
01338 
01340                 static const float64_t MACHINE_EPSILON;
01341 
01342                 /* largest and smallest possible float64_t */
01343                 static const float64_t MAX_REAL_NUMBER;
01344                 static const float64_t MIN_REAL_NUMBER;
01345 
01346     protected:
01348                 static int32_t LOGRANGE;
01349 
01351                 static uint32_t seed;
01352 
01353 #ifdef USE_LOGCACHE
01354 
01356                 static int32_t LOGACCURACY;
01358 
01359                 static float64_t* logtable;
01360 #endif
01361 };
01362 
01363 //implementations of template functions
01364 template <class T1,class T2>
01365 void* CMath::parallel_qsort_index(void* p)
01366     {
01367         struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
01368         T1* output=ps->output;
01369         T2* index=ps->index;
01370         int32_t size=ps->size;
01371         int32_t* qsort_threads=ps->qsort_threads;
01372         int32_t sort_limit=ps->sort_limit;
01373         int32_t num_threads=ps->num_threads;
01374 
01375         if (size<2)
01376         {
01377             return NULL;
01378         }
01379 
01380         if (size==2)
01381         {
01382             if (output[0] > output [1])
01383             {
01384                 swap(output[0], output[1]);
01385                 swap(index[0], index[1]);
01386             }
01387             return NULL;
01388         }
01389         //T1 split=output[random(0,size-1)];
01390         T1 split=output[size/2];
01391 
01392         int32_t left=0;
01393         int32_t right=size-1;
01394 
01395         while (left<=right)
01396         {
01397             while (output[left] < split)
01398                 left++;
01399             while (output[right] > split)
01400                 right--;
01401 
01402             if (left<=right)
01403             {
01404                 swap(output[left], output[right]);
01405                 swap(index[left], index[right]);
01406                 left++;
01407                 right--;
01408             }
01409         }
01410         bool lthread_start=false;
01411         bool rthread_start=false;
01412         pthread_t lthread;
01413         pthread_t rthread;
01414         struct thread_qsort<T1,T2> t1;
01415         struct thread_qsort<T1,T2> t2;
01416 
01417         if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
01418             qsort_index(output,index,right+1);
01419         else if (right+1> 1)
01420         {
01421             (*qsort_threads)++;
01422             lthread_start=true;
01423             t1.output=output;
01424             t1.index=index;
01425             t1.size=right+1;
01426             t1.qsort_threads=qsort_threads;
01427             t1.sort_limit=sort_limit;
01428             t1.num_threads=num_threads;
01429             if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
01430             {
01431                 lthread_start=false;
01432                 (*qsort_threads)--;
01433                 qsort_index(output,index,right+1);
01434             }
01435         }
01436 
01437 
01438         if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
01439             qsort_index(&output[left],&index[left], size-left);
01440         else if (size-left> 1)
01441         {
01442             (*qsort_threads)++;
01443             rthread_start=true;
01444             t2.output=&output[left];
01445             t2.index=&index[left];
01446             t2.size=size-left;
01447             t2.qsort_threads=qsort_threads;
01448             t2.sort_limit=sort_limit;
01449             t2.num_threads=num_threads;
01450             if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
01451             {
01452                 rthread_start=false;
01453                 (*qsort_threads)--;
01454                 qsort_index(&output[left],&index[left], size-left);
01455             }
01456         }
01457 
01458         if (lthread_start)
01459         {
01460             pthread_join(lthread, NULL);
01461             (*qsort_threads)--;
01462         }
01463 
01464         if (rthread_start)
01465         {
01466             pthread_join(rthread, NULL);
01467             (*qsort_threads)--;
01468         }
01469 
01470         return NULL;
01471     }
01472 
01473     template <class T1,class T2>
01474 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
01475 {
01476     if (size<=1)
01477         return;
01478 
01479     if (size==2)
01480     {
01481         if (output[0] > output [1])
01482         {
01483             swap(output[0],output[1]);
01484             swap(index[0],index[1]);
01485         }
01486         return;
01487     }
01488     //T1 split=output[random(0,size-1)];
01489     T1 split=output[size/2];
01490 
01491     int32_t left=0;
01492     int32_t right=size-1;
01493 
01494     while (left<=right)
01495     {
01496         while (output[left] < split)
01497             left++;
01498         while (output[right] > split)
01499             right--;
01500 
01501         if (left<=right)
01502         {
01503             swap(output[left],output[right]);
01504             swap(index[left],index[right]);
01505             left++;
01506             right--;
01507         }
01508     }
01509 
01510     if (right+1> 1)
01511         qsort_index(output,index,right+1);
01512 
01513     if (size-left> 1)
01514         qsort_index(&output[left],&index[left], size-left);
01515 }
01516 
01517     template <class T1,class T2>
01518 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
01519 {
01520     if (size<=1)
01521         return;
01522 
01523     if (size==2)
01524     {
01525         if (output[0] < output [1])
01526         {
01527             swap(output[0],output[1]);
01528             swap(index[0],index[1]);
01529         }
01530         return;
01531     }
01532 
01533     //T1 split=output[random(0,size-1)];
01534     T1 split=output[size/2];
01535 
01536     int32_t left=0;
01537     int32_t right=size-1;
01538 
01539     while (left<=right)
01540     {
01541         while (output[left] > split)
01542             left++;
01543         while (output[right] < split)
01544             right--;
01545 
01546         if (left<=right)
01547         {
01548             swap(output[left],output[right]);
01549             swap(index[left],index[right]);
01550             left++;
01551             right--;
01552         }
01553     }
01554 
01555     if (right+1> 1)
01556         qsort_backward_index(output,index,right+1);
01557 
01558     if (size-left> 1)
01559         qsort_backward_index(&output[left],&index[left], size-left);
01560 }
01561 
01562     template <class T>
01563 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
01564 {
01565     if (6*n*size<13*size*CMath::log(size))
01566         for (int32_t i=0; i<n; i++)
01567             min(&output[i], &index[i], size-i) ;
01568     else
01569         qsort_index(output, index, size) ;
01570 }
01571 
01572 /* move the smallest entry in the array to the beginning */
01573     template <class T>
01574 void CMath::min(float64_t* output, T* index, int32_t size)
01575 {
01576     if (size<=1)
01577         return;
01578     float64_t min_elem=output[0];
01579     int32_t min_index=0;
01580     for (int32_t i=1; i<size; i++)
01581     {
01582         if (output[i]<min_elem)
01583         {
01584             min_index=i;
01585             min_elem=output[i];
01586         }
01587     }
01588     swap(output[0], output[min_index]);
01589     swap(index[0], index[min_index]);
01590 }
01591 
01592 #define COMPLEX128_ERROR_ONEARG_T(function) \
01593 template <> \
01594 inline complex128_t CMath::function<complex128_t>(complex128_t a)   \
01595 {   \
01596     SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
01597         #function);\
01598     return complex128_t(0.0, 0.0);  \
01599 }
01600 
01601 #define COMPLEX128_ERROR_TWOARGS_T(function) \
01602 template <> \
01603 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b)   \
01604 {   \
01605     SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
01606         #function);\
01607     return complex128_t(0.0, 0.0);  \
01608 }
01609 
01610 #define COMPLEX128_ERROR_THREEARGS_T(function) \
01611 template <> \
01612 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c)   \
01613 {   \
01614     SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
01615         #function);\
01616     return complex128_t(0.0, 0.0);  \
01617 }
01618 
01619 #define COMPLEX128_ERROR_SORT_T(function)   \
01620 template <> \
01621 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b)  \
01622 {   \
01623     SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
01624         #function);\
01625 }
01626 
01628 COMPLEX128_ERROR_TWOARGS_T(min)
01629 
01630 
01631 COMPLEX128_ERROR_TWOARGS_T(max)
01632 
01634 COMPLEX128_ERROR_THREEARGS_T(clamp)
01635 
01637 // COMPLEX128_ERROR_ONEARG_T(sign)
01638 
01640 COMPLEX128_ERROR_SORT_T(qsort)
01641 
01643 COMPLEX128_ERROR_SORT_T(insertion_sort)
01644 
01646 COMPLEX128_ERROR_SORT_T(radix_sort)
01647 
01648 }
01649 #undef COMPLEX128_ERROR_ONEARG
01650 #undef COMPLEX128_ERROR_ONEARG_T
01651 #undef COMPLEX128_ERROR_TWOARGS_T
01652 #undef COMPLEX128_ERROR_THREEARGS_T
01653 #undef COMPLEX128_STDMATH
01654 #undef COMPLEX128_ERROR_SORT_T
01655 #endif 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation