SHOGUN  4.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Math.h
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Thoralf Klein
8  * Written (W) 2013 Soumyajit De
9  * Written (W) 2012 Fernando Jose Iglesias Garcia
10  * Written (W) 2011 Siddharth Kherada
11  * Written (W) 2011 Justin Patera
12  * Written (W) 2011 Alesis Novik
13  * Written (W) 2011-2013 Heiko Strathmann
14  * Written (W) 1999-2009 Soeren Sonnenburg
15  * Written (W) 1999-2008 Gunnar Raetsch
16  * Written (W) 2007 Konrad Rieck
17  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
18  */
19 
20 #ifndef __MATHEMATICS_H_
21 #define __MATHEMATICS_H_
22 
23 #include <shogun/lib/config.h>
24 
25 #include <shogun/base/SGObject.h>
26 #include <shogun/lib/common.h>
27 #include <shogun/io/SGIO.h>
28 #include <shogun/base/Parallel.h>
30 #include <shogun/lib/SGVector.h>
31 #include <algorithm>
32 
33 #ifndef _USE_MATH_DEFINES
34 #define _USE_MATH_DEFINES
35 #endif
36 
37 #include <math.h>
38 #include <float.h>
39 #include <sys/types.h>
40 #ifndef _WIN32
41 #include <unistd.h>
42 #endif
43 
44 #ifdef HAVE_PTHREAD
45 #include <pthread.h>
46 #endif
47 
48 #ifdef SUNOS
49 #include <ieeefp.h>
50 #endif
51 
53 #ifdef log2
54 #define cygwin_log2 log2
55 #undef log2
56 #endif
57 
58 #ifndef M_PI
59 #define M_PI 3.14159265358979323846
60 #endif
61 
62 #ifdef _WIN32
63 #ifndef isnan
64 #define isnan _isnan
65 #endif
66 
67 #ifndef isfinite
68 #define isfinite _isfinite
69 #endif
70 #endif //_WIN32
71 
72 /* Size of RNG seed */
73 #define RNG_SEED_SIZE 256
74 
75 /* Maximum stack size */
76 #define RADIX_STACK_SIZE 512
77 
78 /* Stack macros */
79 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
80 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
81 
82 #ifndef DOXYGEN_SHOULD_SKIP_THIS
83 
84 template <class T> struct radix_stack_t
85 {
87  T *sa;
89  size_t sn;
91  uint16_t si;
92 };
93 
95 template <class T1, class T2> struct thread_qsort
96 {
98  T1* output;
100  T2* index;
102  uint32_t size;
103 
105  int32_t* qsort_threads;
107  int32_t sort_limit;
109  int32_t num_threads;
110 };
111 #endif // DOXYGEN_SHOULD_SKIP_THIS
112 
113 #define COMPLEX128_ERROR_ONEARG(function) \
114 static inline complex128_t function(complex128_t a) \
115 { \
116  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
117  #function);\
118  return complex128_t(0.0, 0.0); \
119 }
120 
121 #define COMPLEX128_STDMATH(function) \
122 static inline complex128_t function(complex128_t a) \
123 { \
124  return std::function(a); \
125 }
126 
127 namespace shogun
128 {
130  extern CRandom* sg_rand;
131  class CSGObject;
134 class CMath : public CSGObject
135 {
136  public:
140  CMath();
142 
144  virtual ~CMath();
146 
150 
156  template <class T>
157  static inline T min(T a, T b)
158  {
159  return (a<=b) ? a : b;
160  }
161 
167  template <class T>
168  static inline T max(T a, T b)
169  {
170  return (a>=b) ? a : b;
171  }
172 
178  template <class T>
179  static inline T abs(T a)
180  {
181  // can't be a>=0?(a):(-a), because compiler complains about
182  // 'comparison always true' when T is unsigned
183  if (a==0)
184  return 0;
185  else if (a>0)
186  return a;
187  else
188  return -a;
189  }
190 
195  static inline float64_t abs(complex128_t a)
196  {
197  float64_t a_real=a.real();
198  float64_t a_imag=a.imag();
199  return (CMath::sqrt(a_real*a_real+a_imag*a_imag));
200  }
202 
208  template <class T>
209  static T min(T* vec, int32_t len)
210  {
211  ASSERT(len>0)
212  T minv=vec[0];
213 
214  for (int32_t i=1; i<len; i++)
215  minv=min(vec[i], minv);
216 
217  return minv;
218  }
219 
225  template <class T>
226  static T max(T* vec, int32_t len)
227  {
228  ASSERT(len>0)
229  T maxv=vec[0];
230 
231  for (int32_t i=1; i<len; i++)
232  maxv=max(vec[i], maxv);
233 
234  return maxv;
235  }
236 
243  template <class T>
244  static inline T clamp(T value, T lb, T ub)
245  {
246  if (value<=lb)
247  return lb;
248  else if (value>=ub)
249  return ub;
250  else
251  return value;
252  }
253 
261  template <class T>
262  static int32_t arg_max(T * vec, int32_t inc, int32_t len, T * maxv_ptr = NULL)
263  {
264  ASSERT(len > 0 || inc > 0)
265 
266  T maxv = vec[0];
267  int32_t maxIdx = 0;
268 
269  for (int32_t i = 1, j = inc ; i < len ; i++, j += inc)
270  {
271  if (vec[j] > maxv)
272  maxv = vec[j], maxIdx = i;
273  }
274 
275  if (maxv_ptr != NULL)
276  *maxv_ptr = maxv;
277 
278  return maxIdx;
279  }
280 
288  template <class T>
289  static int32_t arg_min(T * vec, int32_t inc, int32_t len, T * minv_ptr = NULL)
290  {
291  ASSERT(len > 0 || inc > 0)
292 
293  T minv = vec[0];
294  int32_t minIdx = 0;
295 
296  for (int32_t i = 1, j = inc ; i < len ; i++, j += inc)
297  {
298  if (vec[j] < minv)
299  minv = vec[j], minIdx = i;
300  }
301 
302  if (minv_ptr != NULL)
303  *minv_ptr = minv;
304 
305  return minIdx;
306  }
307 
310 
317  template <class T>
318  static inline bool fequals_abs(const T& a, const T& b,
319  const float64_t eps)
320  {
321  const T diff = CMath::abs<T>((a-b));
322  return (diff < eps);
323  }
324 
334  template <class T>
335  static inline bool fequals(const T& a, const T& b,
336  const float64_t eps, bool tolerant=false)
337  {
338  const T absA = CMath::abs<T>(a);
339  const T absB = CMath::abs<T>(b);
340  const T diff = CMath::abs<T>((a-b));
341  T comp;
342 
343  // Handle this separately since NAN is unordered
345  return true;
346 
347  // Required for JSON Serialization Tests
348  if (tolerant)
349  return CMath::fequals_abs<T>(a, b, eps);
350 
351  // handles float32_t and float64_t separately
352  if (sizeof(T) == 4)
354 
355  else
357 
358  if (a==b)
359  return true;
360 
361  // both a and b are 0 and relative error is less meaningful
362  else if ( (a==0) || (b==0) || (diff < comp) )
363  return (diff<(eps * comp));
364 
365  // use max(relative error, diff) to handle large eps
366  else
367  {
368  T check = ((diff/(absA + absB)) > diff)?
369  (diff/(absA + absB)):diff;
370  return (check < eps);
371  }
372  }
373 
374  /* Get the corresponding absolute tolerance for unit test given a relative tolerance
375  *
376  * Note that a unit test will be passed only when
377  * \f[
378  * |v_\text{true} - v_\text{predict}| \leq tol_\text{relative} * |v_\text{true}|
379  * \f] which is equivalent to
380  * \f[
381  * |v_\text{true} - v_\text{predict}| \leq tol_\text{absolute}
382  * \f] where
383  * \f[
384  * tol_\text{absolute} = tol_\text{relative} * |v_\text{true}|
385  * \f]
386  *
387  * @param true_value true value should be finite (neither NAN nor INF)
388  * @param rel_tolerance relative tolerance should be positive and less than 1.0
389  *
390  * @return the corresponding absolute tolerance
391  */
392  static float64_t get_abs_tolerance(float64_t true_value, float64_t rel_tolerance);
393 
398  static inline float64_t round(float64_t d)
399  {
400  return ::floor(d+0.5);
401  }
402 
407  static inline float64_t floor(float64_t d)
408  {
409  return ::floor(d);
410  }
411 
416  static inline float64_t ceil(float64_t d)
417  {
418  return ::ceil(d);
419  }
420 
425  template <class T>
426  static inline T sign(T a)
427  {
428  if (a==0)
429  return 0;
430  else return (a<0) ? (-1) : (+1);
431  }
432 
437  template <class T>
438  static inline void swap(T &a,T &b)
439  {
440  T c=a;
441  a=b;
442  b=c;
443  }
444 
449  template <class T>
450  static inline T sq(T x)
451  {
452  return x*x;
453  }
454 
459  static inline float32_t sqrt(float32_t x)
460  {
461  return ::sqrtf(x);
462  }
463 
468  static inline float64_t sqrt(float64_t x)
469  {
470  return ::sqrt(x);
471  }
472 
477  static inline floatmax_t sqrt(floatmax_t x)
478  {
479  //fall back to double precision sqrt if sqrtl is not
480  //available
481 #ifdef HAVE_SQRTL
482  return ::sqrtl(x);
483 #else
484  return ::sqrt(x);
485 #endif
486  }
487 
490 
491 
495  static inline float32_t invsqrt(float32_t x)
496  {
497  union float_to_int
498  {
499  float32_t f;
500  int32_t i;
501  };
502 
503  float_to_int tmp;
504  tmp.f=x;
505 
506  float32_t xhalf = 0.5f * x;
507  // store floating-point bits in integer tmp.i
508  // and do initial guess for Newton's method
509  tmp.i = 0x5f3759d5 - (tmp.i >> 1);
510  x = tmp.f; // convert new bits into float
511  x = x*(1.5f - xhalf*x*x); // One round of Newton's method
512  return x;
513  }
514 
524  static inline floatmax_t powl(floatmax_t x, floatmax_t n)
525  {
526  //fall back to double precision pow if powl is not
527  //available
528 #ifdef HAVE_POWL
529  return ::powl((long double) x, (long double) n);
530 #else
531  return ::pow((double) x, (double) n);
532 #endif
533  }
534 
535  static inline int32_t pow(bool x, int32_t n)
536  {
537  return 0;
538  }
539 
544  static inline int32_t pow(int32_t x, int32_t n)
545  {
546  ASSERT(n>=0)
547  int32_t result=1;
548  while (n--)
549  result*=x;
550 
551  return result;
552  }
553 
558  static inline float64_t pow(float64_t x, int32_t n)
559  {
560  if (n>=0)
561  {
562  float64_t result=1;
563  while (n--)
564  result*=x;
565 
566  return result;
567  }
568  else
569  return ::pow((double)x, (double)n);
570  }
571 
576  static inline float64_t pow(float64_t x, float64_t n)
577  {
578  return ::pow((double) x, (double) n);
579  }
580 
585  static inline complex128_t pow(complex128_t x, int32_t n)
586  {
587  return std::pow(x, n);
588  }
589 
595  {
596  return std::pow(x, n);
597  }
598 
604  {
605  return std::pow(x, n);
606  }
607 
613  {
614  return std::pow(x, n);
615  }
617 
621  static inline float64_t exp(float64_t x)
622  {
623  return ::exp((double) x);
624  }
625 
627  static inline float64_t dot(const bool* v1, const bool* v2, int32_t n)
628  {
629  float64_t r=0;
630  for (int32_t i=0; i<n; i++)
631  r+=((v1[i]) ? 1 : 0) * ((v2[i]) ? 1 : 0);
632  return r;
633  }
634 
636  static inline floatmax_t dot(const floatmax_t* v1, const floatmax_t* v2, int32_t n)
637  {
638  floatmax_t r=0;
639  for (int32_t i=0; i<n; i++)
640  r+=v1[i]*v2[i];
641  return r;
642  }
643 
644 
646  static float64_t dot(const float64_t* v1, const float64_t* v2, int32_t n);
647 
649  static float32_t dot(const float32_t* v1, const float32_t* v2, int32_t n);
650 
652  static inline float64_t dot(
653  const uint64_t* v1, const uint64_t* v2, int32_t n)
654  {
655  float64_t r=0;
656  for (int32_t i=0; i<n; i++)
657  r+=((float64_t) v1[i])*v2[i];
658 
659  return r;
660  }
662  static inline float64_t dot(
663  const int64_t* v1, const int64_t* v2, int32_t n)
664  {
665  float64_t r=0;
666  for (int32_t i=0; i<n; i++)
667  r+=((float64_t) v1[i])*v2[i];
668 
669  return r;
670  }
671 
673  static inline float64_t dot(
674  const int32_t* v1, const int32_t* v2, int32_t n)
675  {
676  float64_t r=0;
677  for (int32_t i=0; i<n; i++)
678  r+=((float64_t) v1[i])*v2[i];
679 
680  return r;
681  }
682 
684  static inline float64_t dot(
685  const uint32_t* v1, const uint32_t* v2, int32_t n)
686  {
687  float64_t r=0;
688  for (int32_t i=0; i<n; i++)
689  r+=((float64_t) v1[i])*v2[i];
690 
691  return r;
692  }
693 
695  static inline float64_t dot(
696  const uint16_t* v1, const uint16_t* v2, int32_t n)
697  {
698  float64_t r=0;
699  for (int32_t i=0; i<n; i++)
700  r+=((float64_t) v1[i])*v2[i];
701 
702  return r;
703  }
704 
706  static inline float64_t dot(
707  const int16_t* v1, const int16_t* v2, int32_t n)
708  {
709  float64_t r=0;
710  for (int32_t i=0; i<n; i++)
711  r+=((float64_t) v1[i])*v2[i];
712 
713  return r;
714  }
715 
717  static inline float64_t dot(
718  const char* v1, const char* v2, int32_t n)
719  {
720  float64_t r=0;
721  for (int32_t i=0; i<n; i++)
722  r+=((float64_t) v1[i])*v2[i];
723 
724  return r;
725  }
726 
728  static inline float64_t dot(
729  const uint8_t* v1, const uint8_t* v2, int32_t n)
730  {
731  float64_t r=0;
732  for (int32_t i=0; i<n; i++)
733  r+=((float64_t) v1[i])*v2[i];
734 
735  return r;
736  }
737 
739  static inline float64_t dot(
740  const int8_t* v1, const int8_t* v2, int32_t n)
741  {
742  float64_t r=0;
743  for (int32_t i=0; i<n; i++)
744  r+=((float64_t) v1[i])*v2[i];
745 
746  return r;
747  }
748 
750  static inline float64_t dot(
751  const float64_t* v1, const char* v2, int32_t n)
752  {
753  float64_t r=0;
754  for (int32_t i=0; i<n; i++)
755  r+=((float64_t) v1[i])*v2[i];
756 
757  return r;
758  }
759 
762 
763 
771  static inline float64_t tan(float64_t x)
772  {
773  return ::tan((double) x);
774  }
775 
778 
779 
783  static inline float64_t atan(float64_t x)
784  {
785  return ::atan((double) x);
786  }
787 
790 
791 
796  static inline float64_t atan2(float64_t y, float64_t x)
797  {
798  return ::atan2((double) y, (double) x);
799  }
800 
803 
804 
808  static inline float64_t tanh(float64_t x)
809  {
810  return ::tanh((double) x);
811  }
812 
815 
816 
820  static inline float64_t sin(float64_t x)
821  {
822  return ::sin(x);
823  }
824 
827 
828 
832  static inline float64_t asin(float64_t x)
833  {
834  return ::asin(x);
835  }
836 
839 
840 
844  static inline float64_t sinh(float64_t x)
845  {
846  return ::sinh(x);
847  }
848 
851 
852 
856  static inline float64_t cos(float64_t x)
857  {
858  return ::cos(x);
859  }
860 
863 
864 
868  static inline float64_t acos(float64_t x)
869  {
870  return ::acos(x);
871  }
872 
875 
876 
880  static inline float64_t cosh(float64_t x)
881  {
882  return ::cosh(x);
883  }
884 
888 
897  static inline float64_t log10(float64_t v)
898  {
899  return ::log(v)/::log(10.0);
900  }
901 
904 
905 
909  static inline float64_t log2(float64_t v)
910  {
911 #ifdef HAVE_LOG2
912  return ::log2(v);
913 #else
914  return ::log(v)/::log(2.0);
915 #endif //HAVE_LOG2
916  }
917 
922  static inline float64_t log(float64_t v)
923  {
924  return ::log(v);
925  }
926 
929 
930  static inline index_t floor_log(index_t n)
931  {
932  index_t i;
933  for (i = 0; n != 0; i++)
934  n >>= 1;
935 
936  return i;
937  }
939 
946  static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
947  {
948  ASSERT(len>0 && xy)
949 
950  float64_t area = 0.0;
951 
952  if (!reversed)
953  {
954  for (int i=1; i<len; i++)
955  area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
956  }
957  else
958  {
959  for (int i=1; i<len; i++)
960  area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
961  }
962 
963  return area;
964  }
965 
971  static bool strtof(const char* str, float32_t* float_result);
972 
978  static bool strtod(const char* str, float64_t* double_result);
979 
985  static bool strtold(const char* str, floatmax_t* long_double_result);
986 
991  static inline int64_t factorial(int32_t n)
992  {
993  int64_t res=1;
994  for (int i=2; i<=n; i++)
995  res*=i ;
996  return res ;
997  }
998 
1006  static void init_random(uint32_t initseed=0)
1007  {
1008  if (initseed==0)
1010  else
1011  seed=initseed;
1012 
1013  sg_rand->set_seed(seed);
1014  }
1015 
1019  static inline uint64_t random()
1020  {
1021  return sg_rand->random_64();
1022  }
1023 
1027  static inline uint64_t random(uint64_t min_value, uint64_t max_value)
1028  {
1029  return sg_rand->random(min_value, max_value);
1030  }
1031 
1037  static inline int64_t random(int64_t min_value, int64_t max_value)
1038  {
1039  return sg_rand->random(min_value, max_value);
1040  }
1041 
1047  static inline uint32_t random(uint32_t min_value, uint32_t max_value)
1048  {
1049  return sg_rand->random(min_value, max_value);
1050  }
1051 
1057  static inline int32_t random(int32_t min_value, int32_t max_value)
1058  {
1059  return sg_rand->random(min_value, max_value);
1060  }
1061 
1067  static inline float32_t random(float32_t min_value, float32_t max_value)
1068  {
1069  return sg_rand->random(min_value, max_value);
1070  }
1071 
1077  static inline float64_t random(float64_t min_value, float64_t max_value)
1078  {
1079  return sg_rand->random(min_value, max_value);
1080  }
1081 
1087  static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
1088  {
1089  return sg_rand->random(min_value, max_value);
1090  }
1091 
1095  static inline float32_t normal_random(float32_t mean, float32_t std_dev)
1096  {
1097  // sets up variables & makes sure rand_s.range == (0,1)
1098  float32_t ret;
1099  float32_t rand_u;
1100  float32_t rand_v;
1101  float32_t rand_s;
1102  do
1103  {
1104  rand_u = CMath::random(-1.0, 1.0);
1105  rand_v = CMath::random(-1.0, 1.0);
1106  rand_s = rand_u*rand_u + rand_v*rand_v;
1107  } while ((rand_s == 0) || (rand_s >= 1));
1108 
1109  // the meat & potatos, and then the mean & standard deviation shifting...
1110  ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
1111  ret = std_dev*ret + mean;
1112  return ret;
1113  }
1114 
1118  static inline float64_t normal_random(float64_t mean, float64_t std_dev)
1119  {
1120  return sg_rand->normal_distrib(mean, std_dev);
1121  }
1122 
1125  static inline float32_t randn_float()
1126  {
1127  return normal_random(0.0, 1.0);
1128  }
1129 
1132  static inline float64_t randn_double()
1133  {
1134  return sg_rand->std_normal_distrib();
1135  }
1137 
1143  template <class T>
1144  static void permute(SGVector<T> v, CRandom* rand=NULL)
1145  {
1146  if (rand)
1147  {
1148  for (index_t i=0; i<v.vlen; ++i)
1149  swap(v[i], v[rand->random(i, v.vlen-1)]);
1150  }
1151  else
1152  {
1153  for (index_t i=0; i<v.vlen; ++i)
1154  swap(v[i], v[random(i, v.vlen-1)]);
1155  }
1156  }
1157 
1163  template <class T>
1164  static int32_t get_num_nonzero(T* vec, int32_t len)
1165  {
1166  int32_t nnz = 0;
1167  for (index_t i=0; i<len; ++i)
1168  nnz += vec[i] != 0;
1169 
1170  return nnz;
1171  }
1172 
1178  static int32_t get_num_nonzero(complex128_t* vec, int32_t len)
1179  {
1180  int32_t nnz=0;
1181  for (index_t i=0; i<len; ++i)
1182  nnz+=vec[i]!=0.0;
1183  return nnz;
1184  }
1185 
1191  static inline int64_t nchoosek(int32_t n, int32_t k)
1192  {
1193  int64_t res=1;
1194 
1195  for (int32_t i=n-k+1; i<=n; i++)
1196  res*=i;
1197 
1198  return res/factorial(k);
1199  }
1200 
1207  static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100);
1208 
1215  template <class T>
1216  static float64_t* linspace(T start, T end, int32_t n)
1217  {
1218  float64_t* output = SG_MALLOC(float64_t, n);
1219  linspace(output, start, end, n);
1220 
1221  return output;
1222  }
1223 
1230  template <class T>
1231  static SGVector<float64_t> linspace_vec(T start, T end, int32_t n)
1232  {
1233  return SGVector<float64_t>(linspace(start, end, n), n);
1234  }
1235 
1241  template <class T>
1242  static T log_sum_exp(SGVector<T> values)
1243  {
1244  REQUIRE(values.vector, "Values are empty");
1245  REQUIRE(values.vlen>0,"number of values supplied is 0\n");
1246 
1247  /* find minimum element index */
1248  index_t min_index=0;
1249  T X0=values[0];
1250  if (values.vlen>1)
1251  {
1252  for (index_t i=1; i<values.vlen; ++i)
1253  {
1254  if (values[i]<X0)
1255  {
1256  X0=values[i];
1257  min_index=i;
1258  }
1259  }
1260  }
1261 
1262  /* remove element from vector copy and compute log sum exp */
1263  SGVector<T> values_without_X0(values.vlen-1);
1264  index_t from_idx=0;
1265  index_t to_idx=0;
1266  for (from_idx=0; from_idx<values.vlen; ++from_idx)
1267  {
1268  if (from_idx!=min_index)
1269  {
1270  values_without_X0[to_idx]=exp(values[from_idx]-X0);
1271  to_idx++;
1272  }
1273  }
1274 
1275  return X0+log(SGVector<T>::sum(values_without_X0)+1);
1276  }
1277 
1284  template <class T>
1285  static T log_mean_exp(SGVector<T> values)
1286  {
1287  return log_sum_exp(values) - log(values.vlen);
1288  }
1289 
1297  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
1298 
1305  static void sort(float64_t *a, int32_t*idx, int32_t N);
1306 
1312  template <class T>
1313  static void qsort(T* output, int32_t size)
1314  {
1315  if (size<=1)
1316  return;
1317 
1318  if (size==2)
1319  {
1320  if (output[0] > output [1])
1321  CMath::swap(output[0],output[1]);
1322  return;
1323  }
1324  //T split=output[random(0,size-1)];
1325  T split=output[size/2];
1326 
1327  int32_t left=0;
1328  int32_t right=size-1;
1329 
1330  while (left<=right)
1331  {
1332  while (output[left] < split)
1333  left++;
1334  while (output[right] > split)
1335  right--;
1336 
1337  if (left<=right)
1338  {
1339  CMath::swap(output[left],output[right]);
1340  left++;
1341  right--;
1342  }
1343  }
1344 
1345  if (right+1> 1)
1346  qsort(output,right+1);
1347 
1348  if (size-left> 1)
1349  qsort(&output[left],size-left);
1350  }
1351 
1357  template <class T>
1358  static void insertion_sort(T* output, int32_t size)
1359  {
1360  for (int32_t i=0; i<size-1; i++)
1361  {
1362  int32_t j=i-1;
1363  T value=output[i];
1364  while (j >= 0 && output[j] > value)
1365  {
1366  output[j+1] = output[j];
1367  j--;
1368  }
1369  output[j+1]=value;
1370  }
1371  }
1372 
1377  template <class T>
1378  inline static void radix_sort(T* array, int32_t size)
1379  {
1380  radix_sort_helper(array,size,0);
1381  }
1382 
1390  template <class T>
1391  static inline uint8_t byte(T word, uint16_t p)
1392  {
1393  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
1394  }
1395 
1397  static inline uint8_t byte(complex128_t word, uint16_t p)
1398  {
1399  SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
1400  return uint8_t(0);
1401  }
1402 
1408  template <class T>
1409  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
1410  {
1411  static size_t count[256], nc, cmin;
1412  T *ak;
1413  uint8_t c=0;
1414  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
1415  T *an, *aj, *pile[256];
1416  size_t *cp, cmax;
1417 
1418  /* Push initial array to stack */
1419  sp = s;
1420  radix_push(array, size, i);
1421 
1422  /* Loop until all digits have been sorted */
1423  while (sp>s) {
1424  radix_pop(array, size, i);
1425  an = array + size;
1426 
1427  /* Make character histogram */
1428  if (nc == 0) {
1429  cmin = 0xff;
1430  for (ak = array; ak < an; ak++) {
1431  c = byte(*ak, i);
1432  count[c]++;
1433  if (count[c] == 1) {
1434  /* Determine smallest character */
1435  if (c < cmin)
1436  cmin = c;
1437  nc++;
1438  }
1439  }
1440 
1441  /* Sort recursively if stack size too small */
1442  if (sp + nc > s + RADIX_STACK_SIZE) {
1443  radix_sort_helper(array, size, i);
1444  continue;
1445  }
1446  }
1447 
1448  /*
1449  * Set pile[]; push incompletely sorted bins onto stack.
1450  * pile[] = pointers to last out-of-place element in bins.
1451  * Before permuting: pile[c-1] + count[c] = pile[c];
1452  * during deal: pile[c] counts down to pile[c-1].
1453  */
1454  olds = bigs = sp;
1455  cmax = 2;
1456  ak = array;
1457  pile[0xff] = an;
1458  for (cp = count + cmin; nc > 0; cp++) {
1459  /* Find next non-empty pile */
1460  while (*cp == 0)
1461  cp++;
1462  /* Pile with several entries */
1463  if (*cp > 1) {
1464  /* Determine biggest pile */
1465  if (*cp > cmax) {
1466  cmax = *cp;
1467  bigs = sp;
1468  }
1469 
1470  if (i < sizeof(T)-1)
1471  radix_push(ak, *cp, (uint16_t) (i + 1));
1472  }
1473  pile[cp - count] = ak += *cp;
1474  nc--;
1475  }
1476 
1477  /* Play it safe -- biggest bin last. */
1478  swap(*olds, *bigs);
1479 
1480  /*
1481  * Permute misplacements home. Already home: everything
1482  * before aj, and in pile[c], items from pile[c] on.
1483  * Inner loop:
1484  * r = next element to put in place;
1485  * ak = pile[r[i]] = location to put the next element.
1486  * aj = bottom of 1st disordered bin.
1487  * Outer loop:
1488  * Once the 1st disordered bin is done, ie. aj >= ak,
1489  * aj<-aj + count[c] connects the bins in array linked list;
1490  * reset count[c].
1491  */
1492  aj = array;
1493  while (aj<an)
1494  {
1495  T r;
1496 
1497  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
1498  swap(*ak, r);
1499 
1500  *aj = r;
1501  aj += count[c];
1502  count[c] = 0;
1503  }
1504  }
1505  }
1506 
1508  static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
1509  {
1510  SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
1511  }
1512 
1519  template <class T>
1520  static void qsort(T** vector, index_t length)
1521  {
1522  if (length<=1)
1523  return;
1524 
1525  if (length==2)
1526  {
1527  if (*vector[0]>*vector[1])
1528  swap(vector[0],vector[1]);
1529  return;
1530  }
1531  T* split=vector[length/2];
1532 
1533  int32_t left=0;
1534  int32_t right=length-1;
1535 
1536  while (left<=right)
1537  {
1538  while (*vector[left]<*split)
1539  ++left;
1540  while (*vector[right]>*split)
1541  --right;
1542 
1543  if (left<=right)
1544  {
1545  swap(vector[left],vector[right]);
1546  ++left;
1547  --right;
1548  }
1549  }
1550 
1551  if (right+1>1)
1552  qsort(vector,right+1);
1553 
1554  if (length-left>1)
1555  qsort(&vector[left],length-left);
1556  }
1557 
1559  static void qsort(complex128_t** vector, index_t length)
1560  {
1561  SG_SERROR("CMath::qsort():: Not supported for complex128_t\n");
1562  }
1563 
1567  template <class T>
1568  static void qsort(SGVector<T> vector)
1569  {
1570  qsort<T>(vector, vector.size());
1571  }
1572 
1574  template<class T>
1576  {
1578  IndexSorter(const SGVector<T> *vec) { data = vec->vector; }
1579 
1581  bool operator() (index_t i, index_t j) const
1582  {
1583  return abs(data[i]-data[j])>std::numeric_limits<T>::epsilon()
1584  && data[i]<data[j];
1585  }
1586 
1588  const T* data;
1589  };
1590 
1598  template<class T>
1600  {
1601  IndexSorter<T> cmp(&vector);
1602  SGVector<index_t> idx(vector.size());
1603  for (index_t i=0; i < vector.size(); ++i)
1604  idx[i] = i;
1605 
1606  std::sort(idx.vector, idx.vector+vector.size(), cmp);
1607 
1608  return idx;
1609  }
1610 
1616  template <class T>
1617  static bool is_sorted(SGVector<T> vector)
1618  {
1619  if (vector.size() < 2)
1620  return true;
1621 
1622  for(int32_t i=1; i<vector.size(); i++)
1623  {
1624  if (vector[i-1] > vector[i])
1625  return false;
1626  }
1627 
1628  return true;
1629  }
1630 
1635  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1636  {
1637  ASSERT(width>=0)
1638  for (int i=0; i<width; i++)
1639  {
1640  T mask = ((T) 1)<<(sizeof(T)*8-1);
1641  while (mask)
1642  {
1643  if (mask & word)
1644  SG_SPRINT("1")
1645  else
1646  SG_SPRINT("0")
1647 
1648  mask>>=1;
1649  }
1650  }
1651  }
1652 
1654  static void display_bits(complex128_t word,
1655  int32_t width=8*sizeof(complex128_t))
1656  {
1657  SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
1658  }
1659 
1668  template <class T1,class T2>
1669  static void qsort_index(T1* output, T2* index, uint32_t size);
1670 
1672  template <class T>
1673  static void qsort_index(complex128_t* output, T* index, uint32_t size)
1674  {
1675  SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
1676  }
1677 
1686  template <class T1,class T2>
1687  static void qsort_backward_index(
1688  T1* output, T2* index, int32_t size);
1689 
1691  template <class T>
1693  complex128_t* output, T* index, uint32_t size)
1694  {
1695  SG_SERROR("CMath::qsort_backword_index():: \
1696  Not supported for complex128_t\n");
1697  }
1698 
1710  template <class T1,class T2>
1711  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1712  {
1713  int32_t n=0;
1714  thread_qsort<T1,T2> t;
1715  t.output=output;
1716  t.index=index;
1717  t.size=size;
1718  t.qsort_threads=&n;
1719  t.sort_limit=limit;
1720  t.num_threads=n_threads;
1721  parallel_qsort_index<T1,T2>(&t);
1722  }
1723 
1725  template <class T>
1726  inline static void parallel_qsort_index(complex128_t* output, T* index,
1727  uint32_t size, int32_t n_threads, int32_t limit=0)
1728  {
1729  SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1730  }
1731 
1733  template <class T1,class T2>
1734  static void* parallel_qsort_index(void* p);
1735 
1736 
1743  template <class T>
1744  static void min(float64_t* output, T* index, int32_t size);
1745 
1747  static void min(float64_t* output, complex128_t* index, int32_t size)
1748  {
1749  SG_SERROR("CMath::min():: Not supported for complex128_t\n");
1750  }
1751 
1755  template <class T>
1756  static void nmin(
1757  float64_t* output, T* index, int32_t size, int32_t n);
1758 
1760  static void nmin(float64_t* output, complex128_t* index,
1761  int32_t size, int32_t n)
1762  {
1763  SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
1764  }
1765 
1766 
1767 
1771  template <class T>
1772  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1773  {
1774  int32_t start=0;
1775  int32_t end=size-1;
1776 
1777  if (size<1)
1778  return 0;
1779 
1780  while (start<end)
1781  {
1782  int32_t middle=(start+end)/2;
1783 
1784  if (output[middle]>elem)
1785  end=middle-1;
1786  else if (output[middle]<elem)
1787  start=middle+1;
1788  else
1789  return middle;
1790  }
1791 
1792  return start;
1793  }
1794 
1796  static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
1797  {
1798  SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
1799  return int32_t(0);
1800  }
1801 
1808  template <class T>
1809  static inline int32_t binary_search(T* output, int32_t size, T elem)
1810  {
1811  int32_t ind = binary_search_helper(output, size, elem);
1812  if (ind >= 0 && output[ind] == elem)
1813  return ind;
1814  return -1;
1815  }
1816 
1818  static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
1819  {
1820  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1821  return int32_t(-1);
1822  }
1823 
1824 
1832  template<class T>
1833  static inline int32_t binary_search(T** vector, index_t length,
1834  T* elem)
1835  {
1836  int32_t start=0;
1837  int32_t end=length-1;
1838 
1839  if (length<1)
1840  return -1;
1841 
1842  while (start<end)
1843  {
1844  int32_t middle=(start+end)/2;
1845 
1846  if (*vector[middle]>*elem)
1847  end=middle-1;
1848  else if (*vector[middle]<*elem)
1849  start=middle+1;
1850  else
1851  {
1852  start=middle;
1853  break;
1854  }
1855  }
1856 
1857  if (start>=0&&*vector[start]==*elem)
1858  return start;
1859 
1860  return -1;
1861  }
1862 
1864  static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
1865  {
1866  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1867  return int32_t(-1);
1868  }
1869 
1876  template <class T>
1878  T* output, int32_t size, T elem)
1879  {
1880  int32_t ind = binary_search_helper(output, size, elem);
1881 
1882  if (output[ind]<=elem)
1883  return ind;
1884  if (ind>0 && output[ind-1] <= elem)
1885  return ind-1;
1886  return -1;
1887  }
1888 
1891  int32_t size, complex128_t elem)
1892  {
1893  SG_SERROR("CMath::binary_search_max_lower_equal():: \
1894  Not supported for complex128_t\n");
1895  return int32_t(-1);
1896  }
1897 
1906  static float64_t Align(
1907  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1908 
1909 
1911 
1917  {
1918  return c.real();
1919  }
1920 
1926  {
1927  return c.imag();
1928  }
1929 
1931  inline static uint32_t get_seed()
1932  {
1933  return CMath::seed;
1934  }
1935 
1937  inline static uint32_t get_log_range()
1938  {
1939  return CMath::LOGRANGE;
1940  }
1941 
1942 #ifdef USE_LOGCACHE
1943  inline static uint32_t get_log_accuracy()
1945  {
1946  return CMath::LOGACCURACY;
1947  }
1948 #endif
1949 
1951  static int is_finite(double f);
1952 
1954  static int is_infinity(double f);
1955 
1957  static int is_nan(double f);
1958 
1968 #ifdef USE_LOGCACHE
1969  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1970  {
1971  float64_t diff;
1972 
1973  if (!CMath::is_finite(p))
1974  return q;
1975 
1976  if (!CMath::is_finite(q))
1977  {
1978  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
1979  return NOT_A_NUMBER;
1980  }
1981  diff = p - q;
1982  if (diff > 0)
1983  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1984  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1985  }
1986 
1988  static void init_log_table();
1989 
1991  static int32_t determine_logrange();
1992 
1994  static int32_t determine_logaccuracy(int32_t range);
1995 #else
1997  float64_t p, float64_t q)
1998  {
1999  float64_t diff;
2000 
2001  if (!CMath::is_finite(p))
2002  return q;
2003  if (!CMath::is_finite(q))
2004  return p;
2005  diff = p - q;
2006  if (diff > 0)
2007  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
2008  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
2009  }
2010 #endif
2011 #ifdef USE_LOGSUMARRAY
2012 
2017  static inline float64_t logarithmic_sum_array(
2018  float64_t *p, int32_t len)
2019  {
2020  if (len<=2)
2021  {
2022  if (len==2)
2023  return logarithmic_sum(p[0],p[1]) ;
2024  if (len==1)
2025  return p[0];
2026  return -INFTY ;
2027  }
2028  else
2029  {
2030  float64_t *pp=p ;
2031  if (len%2==1) pp++ ;
2032  for (int32_t j=0; j < len>>1; j++)
2033  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
2034  }
2035  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
2036  }
2037 #endif //USE_LOGSUMARRAY
2038 
2039 
2041  virtual const char* get_name() const { return "Math"; }
2042  public:
2045  static const float64_t NOT_A_NUMBER;
2048  static const float64_t INFTY;
2049  static const float64_t ALMOST_INFTY;
2050 
2053 
2055  static const float64_t PI;
2056 
2059 
2060  /* Largest and smallest possible float64_t */
2063 
2064  /* Floating point Limits, Normalized */
2065  static const float32_t F_MAX_VAL32;
2067  static const float64_t F_MAX_VAL64;
2069 
2070  /* Floating point limits, Denormalized */
2071  static const float32_t F_MIN_VAL32;
2072  static const float64_t F_MIN_VAL64;
2073 
2074  protected:
2076  static int32_t LOGRANGE;
2077 
2079  static uint32_t seed;
2080 
2081 #ifdef USE_LOGCACHE
2082 
2084  static int32_t LOGACCURACY;
2086  static float64_t* logtable;
2088 #endif
2089 };
2090 
2091 //implementations of template functions
2092 template <class T1,class T2>
2094  {
2095  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
2096  T1* output=ps->output;
2097  T2* index=ps->index;
2098  int32_t size=ps->size;
2099  int32_t* qsort_threads=ps->qsort_threads;
2100  int32_t sort_limit=ps->sort_limit;
2101  int32_t num_threads=ps->num_threads;
2102 
2103  if (size<2)
2104  {
2105  return NULL;
2106  }
2107 
2108  if (size==2)
2109  {
2110  if (output[0] > output [1])
2111  {
2112  swap(output[0], output[1]);
2113  swap(index[0], index[1]);
2114  }
2115  return NULL;
2116  }
2117  //T1 split=output[random(0,size-1)];
2118  T1 split=output[size/2];
2119 
2120  int32_t left=0;
2121  int32_t right=size-1;
2122 
2123  while (left<=right)
2124  {
2125  while (output[left] < split)
2126  left++;
2127  while (output[right] > split)
2128  right--;
2129 
2130  if (left<=right)
2131  {
2132  swap(output[left], output[right]);
2133  swap(index[left], index[right]);
2134  left++;
2135  right--;
2136  }
2137  }
2138  bool lthread_start=false;
2139  bool rthread_start=false;
2140  pthread_t lthread;
2141  pthread_t rthread;
2142  struct thread_qsort<T1,T2> t1;
2143  struct thread_qsort<T1,T2> t2;
2144 
2145  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
2146  qsort_index(output,index,right+1);
2147  else if (right+1> 1)
2148  {
2149  (*qsort_threads)++;
2150  lthread_start=true;
2151  t1.output=output;
2152  t1.index=index;
2153  t1.size=right+1;
2154  t1.qsort_threads=qsort_threads;
2155  t1.sort_limit=sort_limit;
2156  t1.num_threads=num_threads;
2157  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
2158  {
2159  lthread_start=false;
2160  (*qsort_threads)--;
2161  qsort_index(output,index,right+1);
2162  }
2163  }
2164 
2165 
2166  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
2167  qsort_index(&output[left],&index[left], size-left);
2168  else if (size-left> 1)
2169  {
2170  (*qsort_threads)++;
2171  rthread_start=true;
2172  t2.output=&output[left];
2173  t2.index=&index[left];
2174  t2.size=size-left;
2175  t2.qsort_threads=qsort_threads;
2176  t2.sort_limit=sort_limit;
2177  t2.num_threads=num_threads;
2178  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
2179  {
2180  rthread_start=false;
2181  (*qsort_threads)--;
2182  qsort_index(&output[left],&index[left], size-left);
2183  }
2184  }
2185 
2186  if (lthread_start)
2187  {
2188  pthread_join(lthread, NULL);
2189  (*qsort_threads)--;
2190  }
2191 
2192  if (rthread_start)
2193  {
2194  pthread_join(rthread, NULL);
2195  (*qsort_threads)--;
2196  }
2197 
2198  return NULL;
2199  }
2200 
2201  template <class T1,class T2>
2202 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
2203 {
2204  if (size<=1)
2205  return;
2206 
2207  if (size==2)
2208  {
2209  if (output[0] > output [1])
2210  {
2211  swap(output[0],output[1]);
2212  swap(index[0],index[1]);
2213  }
2214  return;
2215  }
2216  //T1 split=output[random(0,size-1)];
2217  T1 split=output[size/2];
2218 
2219  int32_t left=0;
2220  int32_t right=size-1;
2221 
2222  while (left<=right)
2223  {
2224  while (output[left] < split)
2225  left++;
2226  while (output[right] > split)
2227  right--;
2228 
2229  if (left<=right)
2230  {
2231  swap(output[left],output[right]);
2232  swap(index[left],index[right]);
2233  left++;
2234  right--;
2235  }
2236  }
2237 
2238  if (right+1> 1)
2239  qsort_index(output,index,right+1);
2240 
2241  if (size-left> 1)
2242  qsort_index(&output[left],&index[left], size-left);
2243 }
2244 
2245  template <class T1,class T2>
2246 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
2247 {
2248  if (size<=1)
2249  return;
2250 
2251  if (size==2)
2252  {
2253  if (output[0] < output [1])
2254  {
2255  swap(output[0],output[1]);
2256  swap(index[0],index[1]);
2257  }
2258  return;
2259  }
2260 
2261  //T1 split=output[random(0,size-1)];
2262  T1 split=output[size/2];
2263 
2264  int32_t left=0;
2265  int32_t right=size-1;
2266 
2267  while (left<=right)
2268  {
2269  while (output[left] > split)
2270  left++;
2271  while (output[right] < split)
2272  right--;
2273 
2274  if (left<=right)
2275  {
2276  swap(output[left],output[right]);
2277  swap(index[left],index[right]);
2278  left++;
2279  right--;
2280  }
2281  }
2282 
2283  if (right+1> 1)
2284  qsort_backward_index(output,index,right+1);
2285 
2286  if (size-left> 1)
2287  qsort_backward_index(&output[left],&index[left], size-left);
2288 }
2289 
2290  template <class T>
2291 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
2292 {
2293  if (6*n*size<13*size*CMath::log(size))
2294  for (int32_t i=0; i<n; i++)
2295  min(&output[i], &index[i], size-i);
2296  else
2297  qsort_index(output, index, size);
2298 }
2299 
2300 /* move the smallest entry in the array to the beginning */
2301  template <class T>
2302 void CMath::min(float64_t* output, T* index, int32_t size)
2303 {
2304  if (size<=1)
2305  return;
2306  float64_t min_elem=output[0];
2307  int32_t min_index=0;
2308  for (int32_t i=1; i<size; i++)
2309  {
2310  if (output[i]<min_elem)
2311  {
2312  min_index=i;
2313  min_elem=output[i];
2314  }
2315  }
2316  swap(output[0], output[min_index]);
2317  swap(index[0], index[min_index]);
2318 }
2319 
2321 template <>
2322 inline float64_t* CMath::linspace<complex128_t>(complex128_t start, complex128_t end, int32_t n)
2323 {
2324  SG_SERROR("SGVector::linspace():: Not supported for complex128_t\n");
2325  return NULL;
2326 }
2327 
2328 #define COMPLEX128_ERROR_ONEVECARG_RETURNS_T(function, return_type, return_statement) \
2329 template <> \
2330 inline return_type CMath::function<complex128_t>(SGVector<complex128_t> vector) \
2331 { \
2332  SG_SERROR("CMath::%s():: Not supported for complex128_t\n", \
2333  #function); \
2334  return_statement; \
2335 }
2336 
2337 #define COMPLEX128_ERROR_ONEARG_T(function) \
2338 template <> \
2339 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
2340 { \
2341  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2342  #function);\
2343  return complex128_t(0.0, 0.0); \
2344 }
2345 
2346 #define COMPLEX128_ERROR_TWOARGS_T(function) \
2347 template <> \
2348 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \
2349 { \
2350  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2351  #function);\
2352  return complex128_t(0.0, 0.0); \
2353 }
2354 
2355 #define COMPLEX128_ERROR_THREEARGS_T(function) \
2356 template <> \
2357 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \
2358 { \
2359  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2360  #function);\
2361  return complex128_t(0.0, 0.0); \
2362 }
2363 
2364 #define COMPLEX128_ERROR_SORT_T(function) \
2365 template <> \
2366 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \
2367 { \
2368  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2369  #function);\
2370 }
2371 
2372 #define COMPLEX128_ERROR_ARG_MAX_MIN(function) \
2373 template <> \
2374 inline int32_t CMath::function<complex128_t>(complex128_t * a, int32_t b, int32_t c, complex128_t * d) \
2375 { \
2376  int32_t maxIdx=0; \
2377  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2378  #function);\
2379  return maxIdx; \
2380 }
2381 
2384 
2386 COMPLEX128_ERROR_ONEVECARG_RETURNS_T(argsort, SGVector<index_t>, SGVector<index_t> idx(vector.size());return idx;)
2387 
2389 COMPLEX128_ERROR_ONEVECARG_RETURNS_T(is_sorted, bool, return false;)
2390 
2393 
2396 
2399 
2401 // COMPLEX128_ERROR_ONEARG_T(sign)
2402 
2405 
2408 
2411 
2414 
2417 
2418 }
2419 #undef COMPLEX128_ERROR_ONEARG
2420 #undef COMPLEX128_ERROR_ONEARG_T
2421 #undef COMPLEX128_ERROR_TWOARGS_T
2422 #undef COMPLEX128_ERROR_THREEARGS_T
2423 #undef COMPLEX128_STDMATH
2424 #undef COMPLEX128_ERROR_SORT_T
2425 #endif
static float64_t sin(float64_t x)
tanh(x), x being a complex128_t
Definition: Math.h:820
static float64_t normal_random(float64_t mean, float64_t std_dev)
Definition: Math.h:1118
static const float32_t F_MAX_VAL32
Definition: Math.h:2065
static const float64_t MACHINE_EPSILON
Definition: Math.h:2058
static bool strtof(const char *str, float32_t *float_result)
Definition: Math.cpp:281
static void permute(SGVector< T > v, CRandom *rand=NULL)
Definition: Math.h:1144
static int32_t binary_search(complex128_t *output, int32_t size, complex128_t elem)
binary_search not implemented for complex128_t
Definition: Math.h:1818
static void parallel_qsort_index(T1 *output, T2 *index, uint32_t size, int32_t n_threads, int32_t limit=262144)
Definition: Math.h:1711
static floatmax_t dot(const floatmax_t *v1, const floatmax_t *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
Definition: Math.h:636
static uint32_t seed
random generator seed
Definition: Math.h:2079
std::complex< float64_t > complex128_t
Definition: common.h:67
static int is_finite(double f)
checks whether a float is finite
Definition: Math.cpp:266
static float64_t Align(char *seq1, char *seq2, int32_t l1, int32_t l2, float64_t gapCost)
Definition: Math.cpp:178
static int32_t arg_max(T *vec, int32_t inc, int32_t len, T *maxv_ptr=NULL)
Definition: Math.h:262
static float64_t sqrt(float64_t x)
Definition: Math.h:468
static floatmax_t random(floatmax_t min_value, floatmax_t max_value)
Definition: Math.h:1087
static floatmax_t sqrt(floatmax_t x)
Definition: Math.h:477
static void linspace(float64_t *output, float64_t start, float64_t end, int32_t n=100)
Definition: Math.cpp:221
int32_t index_t
Definition: common.h:62
static float64_t ceil(float64_t d)
Definition: Math.h:416
static bool strtod(const char *str, float64_t *double_result)
Definition: Math.cpp:312
static complex128_t pow(complex128_t x, int32_t n)
Definition: Math.h:585
virtual ~CMath()
Destructor - frees logtable.
Definition: Math.cpp:82
static const float64_t INFTY
infinity
Definition: Math.h:2048
static SGVector< float64_t > linspace_vec(T start, T end, int32_t n)
Definition: Math.h:1231
static int32_t binary_search_helper(T *output, int32_t size, T elem)
Definition: Math.h:1772
#define COMPLEX128_ERROR_THREEARGS_T(function)
Definition: Math.h:2355
static void nmin(float64_t *output, T *index, int32_t size, int32_t n)
Definition: Math.h:2291
static void qsort_index(T1 *output, T2 *index, uint32_t size)
Definition: Math.h:2202
static void qsort_index(complex128_t *output, T *index, uint32_t size)
qsort_index not implemented for complex128_t
Definition: Math.h:1673
static float64_t log10(float64_t v)
Definition: Math.h:897
#define COMPLEX128_ERROR_TWOARGS_T(function)
Definition: Math.h:2346
#define SG_SWARNING(...)
Definition: SGIO.h:178
static void qsort(SGVector< T > vector)
Definition: Math.h:1568
#define COMPLEX128_ERROR_ONEVECARG_RETURNS_T(function, return_type, return_statement)
Definition: Math.h:2328
static float32_t normal_random(float32_t mean, float32_t std_dev)
Definition: Math.h:1095
static float64_t random(float64_t min_value, float64_t max_value)
Definition: Math.h:1077
static int32_t binary_search(T *output, int32_t size, T elem)
Definition: Math.h:1809
static T sq(T x)
Definition: Math.h:450
static uint32_t get_seed()
returns number generator seed
Definition: Math.h:1931
#define COMPLEX128_ERROR_ARG_MAX_MIN(function)
Definition: Math.h:2372
static float32_t randn_float()
Definition: Math.h:1125
static float64_t atan2(float64_t y, float64_t x)
atan(x), x being a complex128_t not implemented
Definition: Math.h:796
static const float64_t MIN_REAL_NUMBER
Definition: Math.h:2062
static float64_t dot(const char *v1, const char *v2, int32_t n)
Compute dot product between v1 and v2 (for 8bit (un)signed ints)
Definition: Math.h:717
static float64_t randn_double()
Definition: Math.h:1132
static int32_t binary_search_max_lower_equal(T *output, int32_t size, T elem)
Definition: Math.h:1877
static float64_t atan(float64_t x)
tan(x), x being a complex128_t
Definition: Math.h:783
#define REQUIRE(x,...)
Definition: SGIO.h:206
static const float64_t F_MAX_VAL64
Definition: Math.h:2067
static float64_t cosh(float64_t x)
acos(x), x being a complex128_t not implemented
Definition: Math.h:880
static uint32_t get_log_range()
returns range of logtable
Definition: Math.h:1937
void split(v_array< ds_node< P > > &point_set, v_array< ds_node< P > > &far_set, int max_scale)
Definition: JLCoverTree.h:150
static float32_t random(float32_t min_value, float32_t max_value)
Definition: Math.h:1067
static const float32_t F_MIN_VAL32
Definition: Math.h:2071
CRandom * sg_rand
Definition: init.cpp:39
static float64_t dot(const int8_t *v1, const int8_t *v2, int32_t n)
Compute dot product between v1 and v2 (for 8bit (un)signed ints)
Definition: Math.h:739
IndexSorter(const SGVector< T > *vec)
Definition: Math.h:1578
static bool fequals_abs(const T &a, const T &b, const float64_t eps)
Definition: Math.h:318
#define RADIX_STACK_SIZE
Definition: Math.h:76
static float64_t imag(complex128_t c)
Definition: Math.h:1925
static float64_t floor(float64_t d)
Definition: Math.h:407
static uint64_t random()
Definition: Math.h:1019
static float64_t * linspace(T start, T end, int32_t n)
Definition: Math.h:1216
static int32_t get_num_nonzero(complex128_t *vec, int32_t len)
Definition: Math.h:1178
static const float32_t F_MIN_NORM_VAL32
Definition: Math.h:2066
static float64_t real(complex128_t c)
Definition: Math.h:1916
static uint8_t byte(complex128_t word, uint16_t p)
byte not implemented for complex128_t
Definition: Math.h:1397
static void qsort(T *output, int32_t size)
Definition: Math.h:1313
static int32_t LOGRANGE
range for logtable: log(1+exp(x)) -LOGRANGE <= x <= 0
Definition: Math.h:2076
static const float64_t ALMOST_NEG_INFTY
almost neg (log) infinity
Definition: Math.h:2052
static bool fequals(const T &a, const T &b, const float64_t eps, bool tolerant=false)
Definition: Math.h:335
static complex128_t pow(complex128_t x, complex128_t n)
Definition: Math.h:594
static float32_t invsqrt(float32_t x)
x^0.5, x being a complex128_t
Definition: Math.h:495
int32_t size() const
Definition: SGVector.h:113
index_t vlen
Definition: SGVector.h:494
static uint8_t byte(T word, uint16_t p)
Definition: Math.h:1391
static complex128_t pow(float64_t x, complex128_t n)
Definition: Math.h:612
#define SG_SPRINT(...)
Definition: SGIO.h:180
CMath()
Constructor - initializes log-table.
Definition: Math.cpp:65
static float64_t pow(float64_t x, int32_t n)
Definition: Math.h:558
#define ASSERT(x)
Definition: SGIO.h:201
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:115
static int32_t random(int32_t min_value, int32_t max_value)
Definition: Math.h:1057
static float64_t pow(float64_t x, float64_t n)
Definition: Math.h:576
shogun vector
static float64_t dot(const uint8_t *v1, const uint8_t *v2, int32_t n)
Compute dot product between v1 and v2 (for 8bit (un)signed ints)
Definition: Math.h:728
bool operator()(index_t i, index_t j) const
Definition: Math.h:1581
static void qsort(T **vector, index_t length)
Definition: Math.h:1520
static void init_random(uint32_t initseed=0)
Definition: Math.h:1006
static int32_t pow(int32_t x, int32_t n)
Definition: Math.h:544
#define radix_pop(a, n, i)
Definition: Math.h:80
static float64_t dot(const uint16_t *v1, const uint16_t *v2, int32_t n)
Compute dot product between v1 and v2 (for 16bit unsigned ints)
Definition: Math.h:695
double float64_t
Definition: common.h:50
static void parallel_qsort_index(complex128_t *output, T *index, uint32_t size, int32_t n_threads, int32_t limit=0)
parallel_qsort_index not implemented for complex128_t
Definition: Math.h:1726
static void min(float64_t *output, complex128_t *index, int32_t size)
complex128_t cannot be used as index
Definition: Math.h:1747
long double floatmax_t
Definition: common.h:51
static float64_t dot(const int16_t *v1, const int16_t *v2, int32_t n)
Compute dot product between v1 and v2 (for 16bit unsigned ints)
Definition: Math.h:706
#define COMPLEX128_ERROR_ONEARG(function)
Definition: Math.h:113
static SGVector< index_t > argsort(SGVector< T > vector)
Definition: Math.h:1599
static float64_t cos(float64_t x)
sinh(x), x being a complex128_t
Definition: Math.h:856
static float64_t dot(const uint32_t *v1, const uint32_t *v2, int32_t n)
Compute dot product between v1 and v2 (for 32bit unsigned ints)
Definition: Math.h:684
static void qsort_backword_index(complex128_t *output, T *index, uint32_t size)
qsort_backword_index not implemented for complex128_t
Definition: Math.h:1692
static int32_t get_num_nonzero(T *vec, int32_t len)
Definition: Math.h:1164
static T log_sum_exp(SGVector< T > values)
Definition: Math.h:1242
static T max(T a, T b)
Definition: Math.h:168
static float64_t area_under_curve(float64_t *xy, int32_t len, bool reversed)
Definition: Math.h:946
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
Definition: Math.h:627
static void display_bits(T word, int32_t width=8 *sizeof(T))
Definition: Math.h:1635
static uint64_t random(uint64_t min_value, uint64_t max_value)
Definition: Math.h:1027
static int64_t factorial(int32_t n)
Definition: Math.h:991
static void qsort(complex128_t **vector, index_t length)
qsort not implemented for complex128_t
Definition: Math.h:1559
static int32_t arg_min(T *vec, int32_t inc, int32_t len, T *minv_ptr=NULL)
Definition: Math.h:289
static float64_t tan(float64_t x)
Definition: Math.h:771
static int32_t binary_search_max_lower_equal(complex128_t *output, int32_t size, complex128_t elem)
binary_search_max_lower_equal not implemented for complex128_t
Definition: Math.h:1890
float float32_t
Definition: common.h:49
static float64_t sinh(float64_t x)
asin(x), x being a complex128_t not implemented
Definition: Math.h:844
: Pseudo random number geneartor
Definition: Random.h:34
static float64_t get_abs_tolerance(float64_t true_value, float64_t rel_tolerance)
Definition: Math.cpp:381
static int64_t nchoosek(int32_t n, int32_t k)
Definition: Math.h:1191
static T min(T *vec, int32_t len)
Definition: Math.h:209
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
static float64_t dot(const int32_t *v1, const int32_t *v2, int32_t n)
Compute dot product between v1 and v2 (for 32bit ints)
Definition: Math.h:673
static int32_t binary_search(T **vector, index_t length, T *elem)
Definition: Math.h:1833
static int is_infinity(double f)
checks whether a float is infinity
Definition: Math.cpp:247
static float64_t acos(float64_t x)
cos(x), x being a complex128_t
Definition: Math.h:868
static float64_t abs(complex128_t a)
Definition: Math.h:195
static float64_t log2(float64_t v)
log10(x), x being a complex128_t
Definition: Math.h:909
static void display_bits(complex128_t word, int32_t width=8 *sizeof(complex128_t))
disply_bits not implemented for complex128_t
Definition: Math.h:1654
#define COMPLEX128_STDMATH(function)
Definition: Math.h:121
static void radix_sort_helper(T *array, int32_t size, uint16_t i)
Definition: Math.h:1409
static T sign(T a)
Definition: Math.h:426
static float64_t tanh(float64_t x)
atan2(x), x being a complex128_t not implemented
Definition: Math.h:808
static int is_nan(double f)
checks whether a float is nan
Definition: Math.cpp:234
static float64_t dot(const float64_t *v1, const char *v2, int32_t n)
Compute dot product between v1 and v2.
Definition: Math.h:750
virtual const char * get_name() const
Definition: Math.h:2041
static float64_t asin(float64_t x)
sin(x), x being a complex128_t
Definition: Math.h:832
static T min(T a, T b)
Definition: Math.h:157
#define SG_SERROR(...)
Definition: SGIO.h:179
static float64_t exp(float64_t x)
Definition: Math.h:621
static const float64_t F_MIN_VAL64
Definition: Math.h:2072
static const float64_t F_MIN_NORM_VAL64
Definition: Math.h:2068
static float64_t log(float64_t v)
Definition: Math.h:922
static void radix_sort_helper(complex128_t *array, int32_t size, uint16_t i)
radix_sort_helper not implemented for complex128_t
Definition: Math.h:1508
static const float64_t ALMOST_INFTY
Definition: Math.h:2049
Class which collects generic mathematical functions.
Definition: Math.h:134
static T log_mean_exp(SGVector< T > values)
Definition: Math.h:1285
static void insertion_sort(T *output, int32_t size)
Definition: Math.h:1358
static void swap(T &a, T &b)
Definition: Math.h:438
static complex128_t pow(complex128_t x, float64_t n)
Definition: Math.h:603
static void sort(int32_t *a, int32_t cols, int32_t sort_col=0)
Definition: Math.cpp:139
static T max(T *vec, int32_t len)
Definition: Math.h:226
static int32_t binary_search_helper(complex128_t *output, int32_t size, complex128_t elem)
binary_search_helper not implemented for complex128_t
Definition: Math.h:1796
static float64_t round(float64_t d)
Definition: Math.h:398
static bool is_sorted(SGVector< T > vector)
Definition: Math.h:1617
static float32_t sqrt(float32_t x)
Definition: Math.h:459
static uint32_t generate_seed()
Definition: Random.cpp:315
static index_t floor_log(index_t n)
log(x), x being a complex128_t
Definition: Math.h:930
static float64_t dot(const int64_t *v1, const int64_t *v2, int32_t n)
Compute dot product between v1 and v2 (for 64bit ints)
Definition: Math.h:662
static void radix_sort(T *array, int32_t size)
Definition: Math.h:1378
static uint32_t random(uint32_t min_value, uint32_t max_value)
Definition: Math.h:1047
static floatmax_t powl(floatmax_t x, floatmax_t n)
Definition: Math.h:524
static float64_t dot(const uint64_t *v1, const uint64_t *v2, int32_t n)
compute dot product between v1 and v2 (for 64bit unsigned ints)
Definition: Math.h:652
static float64_t logarithmic_sum(float64_t p, float64_t q)
Definition: Math.h:1996
#define radix_push(a, n, i)
Definition: Math.h:79
static T clamp(T value, T lb, T ub)
Definition: Math.h:244
static int32_t binary_search(complex128_t **vector, index_t length, complex128_t *elem)
binary_search not implemented for complex128_t
Definition: Math.h:1864
#define COMPLEX128_ERROR_SORT_T(function)
Definition: Math.h:2364
static bool strtold(const char *str, floatmax_t *long_double_result)
Definition: Math.cpp:343
static int32_t pow(bool x, int32_t n)
Definition: Math.h:535
static const float64_t NOT_A_NUMBER
not a number
Definition: Math.h:2046
static void qsort_backward_index(T1 *output, T2 *index, int32_t size)
Definition: Math.h:2246
static const float64_t MAX_REAL_NUMBER
Definition: Math.h:2061
static T abs(T a)
Definition: Math.h:179
static void nmin(float64_t *output, complex128_t *index, int32_t size, int32_t n)
complex128_t cannot be used as index
Definition: Math.h:1760
static int64_t random(int64_t min_value, int64_t max_value)
Definition: Math.h:1037
static const float64_t PI
Definition: Math.h:2055

SHOGUN Machine Learning Toolbox - Documentation