SHOGUN  6.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
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/lib/common.h>
26 #include <shogun/base/Parallel.h>
28 #include <shogun/lib/SGVector.h>
29 #include <algorithm>
30 
31 #ifndef _USE_MATH_DEFINES
32 #define _USE_MATH_DEFINES
33 #endif
34 #include <math.h>
35 #include <cfloat>
36 
37 #ifndef _WIN32
38 #include <unistd.h>
39 #endif
40 
41 #ifdef HAVE_PTHREAD
42 #include <pthread.h>
43 #endif
44 
45 #ifdef SUNOS
46 #include <ieeefp.h>
47 #endif
48 
50 #ifdef log2
51 #define cygwin_log2 log2
52 #undef log2
53 #endif
54 
55 #ifndef M_PI
56 #define M_PI 3.14159265358979323846
57 #endif
58 
59 #ifdef _WIN32
60 #ifndef isnan
61 #define isnan _isnan
62 #endif
63 
64 #ifndef isfinite
65 #define isfinite _isfinite
66 #endif
67 #endif //_WIN32
68 
69 /* Size of RNG seed */
70 #define RNG_SEED_SIZE 256
71 
72 /* Maximum stack size */
73 #define RADIX_STACK_SIZE 512
74 
75 /* Stack macros */
76 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
77 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
78 
79 #ifndef DOXYGEN_SHOULD_SKIP_THIS
80 
81 template <class T> struct radix_stack_t
82 {
84  T *sa;
86  size_t sn;
88  uint16_t si;
89 };
90 
92 template <class T1, class T2> struct thread_qsort
93 {
95  T1* output;
97  T2* index;
99  uint32_t size;
100 
102  int32_t* qsort_threads;
104  int32_t sort_limit;
106  int32_t num_threads;
107 };
108 #endif // DOXYGEN_SHOULD_SKIP_THIS
109 
110 #define COMPLEX128_ERROR_ONEARG(function) \
111 static inline complex128_t function(complex128_t a) \
112 { \
113  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
114  #function);\
115  return complex128_t(0.0, 0.0); \
116 }
117 
118 #define COMPLEX128_STDMATH(function) \
119 static inline complex128_t function(complex128_t a) \
120 { \
121  return std::function(a); \
122 }
123 
124 namespace shogun
125 {
127  extern CRandom* sg_rand;
130 class CMath : public CSGObject
131 {
132  public:
136  CMath();
138 
140  virtual ~CMath();
142 
146 
152  template <class T>
153  static inline T min(T a, T b)
154  {
155  return (a<=b) ? a : b;
156  }
157 
163  template <class T>
164  static inline T max(T a, T b)
165  {
166  return (a>=b) ? a : b;
167  }
168 
174  template <class T>
175  static inline T abs(T a)
176  {
177  // can't be a>=0?(a):(-a), because compiler complains about
178  // 'comparison always true' when T is unsigned
179  if (a==0)
180  return 0;
181  else if (a>0)
182  return a;
183  else
184  return -a;
185  }
186 
191  static inline float64_t abs(complex128_t a)
192  {
193  float64_t a_real=a.real();
194  float64_t a_imag=a.imag();
195  return (CMath::sqrt(a_real*a_real+a_imag*a_imag));
196  }
198 
204  template <class T>
205  static T min(T* vec, int32_t len)
206  {
207  ASSERT(len>0)
208  T minv=vec[0];
209 
210  for (int32_t i=1; i<len; i++)
211  minv=min(vec[i], minv);
212 
213  return minv;
214  }
215 
221  template <class T>
222  static T max(T* vec, int32_t len)
223  {
224  ASSERT(len>0)
225  T maxv=vec[0];
226 
227  for (int32_t i=1; i<len; i++)
228  maxv=max(vec[i], maxv);
229 
230  return maxv;
231  }
232 
239  template <class T>
240  static inline T clamp(T value, T lb, T ub)
241  {
242  if (value<=lb)
243  return lb;
244  else if (value>=ub)
245  return ub;
246  else
247  return value;
248  }
249 
257  template <class T>
258  static int32_t arg_max(T * vec, int32_t inc, int32_t len, T * maxv_ptr = NULL)
259  {
260  ASSERT(len > 0 || inc > 0)
261 
262  T maxv = vec[0];
263  int32_t maxIdx = 0;
264 
265  for (int32_t i = 1, j = inc ; i < len ; i++, j += inc)
266  {
267  if (vec[j] > maxv)
268  maxv = vec[j], maxIdx = i;
269  }
270 
271  if (maxv_ptr != NULL)
272  *maxv_ptr = maxv;
273 
274  return maxIdx;
275  }
276 
284  template <class T>
285  static int32_t arg_min(T * vec, int32_t inc, int32_t len, T * minv_ptr = NULL)
286  {
287  ASSERT(len > 0 || inc > 0)
288 
289  T minv = vec[0];
290  int32_t minIdx = 0;
291 
292  for (int32_t i = 1, j = inc ; i < len ; i++, j += inc)
293  {
294  if (vec[j] < minv)
295  minv = vec[j], minIdx = i;
296  }
297 
298  if (minv_ptr != NULL)
299  *minv_ptr = minv;
300 
301  return minIdx;
302  }
303 
306 
313  template <class T>
314  static inline bool fequals_abs(const T& a, const T& b,
315  const float64_t eps)
316  {
317  const T diff = CMath::abs<T>((a-b));
318  return (diff < eps);
319  }
320 
330  template <class T>
331  static inline bool fequals(const T& a, const T& b,
332  const float64_t eps, bool tolerant=false)
333  {
334  const T absA = CMath::abs<T>(a);
335  const T absB = CMath::abs<T>(b);
336  const T diff = CMath::abs<T>((a-b));
337  T comp;
338 
339  // Handle this separately since NAN is unordered
341  return true;
342 
343  // Required for JSON Serialization Tests
344  if (tolerant)
345  return CMath::fequals_abs<T>(a, b, eps);
346 
347  // handles float32_t and float64_t separately
348  if (sizeof(T) == 4)
350 
351  else
353 
354  if (a == b)
355  return true;
356 
357  // both a and b are 0 and relative error is less meaningful
358  else if ((a == 0) || (b == 0) || (diff < comp))
359  return (diff < (eps * comp));
360  // use max(relative error, diff) to handle large eps
361  else
362  {
363  T check = ((diff/(absA + absB)) > diff)?
364  (diff/(absA + absB)):diff;
365  return (check < eps);
366  }
367  }
368 
369  /* Get the corresponding absolute tolerance for unit test given a relative tolerance
370  *
371  * Note that a unit test will be passed only when
372  * \f[
373  * |v_\text{true} - v_\text{predict}| \leq tol_\text{relative} * |v_\text{true}|
374  * \f] which is equivalent to
375  * \f[
376  * |v_\text{true} - v_\text{predict}| \leq tol_\text{absolute}
377  * \f] where
378  * \f[
379  * tol_\text{absolute} = tol_\text{relative} * |v_\text{true}|
380  * \f]
381  *
382  * @param true_value true value should be finite (neither NAN nor INF)
383  * @param rel_tolerance relative tolerance should be positive and less than 1.0
384  *
385  * @return the corresponding absolute tolerance
386  */
387  static float64_t get_abs_tolerance(float64_t true_value, float64_t rel_tolerance);
388 
393  static inline float64_t round(float64_t d)
394  {
395  return ::floor(d+0.5);
396  }
397 
402  static inline float64_t floor(float64_t d)
403  {
404  return ::floor(d);
405  }
406 
411  static inline float64_t ceil(float64_t d)
412  {
413  return ::ceil(d);
414  }
415 
420  template <class T>
421  static inline T sign(T a)
422  {
423  if (a==0)
424  return 0;
425  else return (a<0) ? (-1) : (+1);
426  }
427 
432  template <class T>
433  static inline void swap(T &a,T &b)
434  {
435  T c=a;
436  a=b;
437  b=c;
438  }
439 
444  template <class T>
445  static inline T sq(T x)
446  {
447  return x*x;
448  }
449 
454  static inline float32_t sqrt(float32_t x)
455  {
456  return ::sqrtf(x);
457  }
458 
463  static inline float64_t sqrt(float64_t x)
464  {
465  return ::sqrt(x);
466  }
467 
472  static inline floatmax_t sqrt(floatmax_t x)
473  {
474  //fall back to double precision sqrt if sqrtl is not
475  //available
476 #ifdef HAVE_SQRTL
477  return ::sqrtl(x);
478 #else
479  return ::sqrt(x);
480 #endif
481  }
482 
485 
486 
490  static inline float32_t invsqrt(float32_t x)
491  {
492  union float_to_int
493  {
494  float32_t f;
495  int32_t i;
496  };
497 
498  float_to_int tmp;
499  tmp.f=x;
500 
501  float32_t xhalf = 0.5f * x;
502  // store floating-point bits in integer tmp.i
503  // and do initial guess for Newton's method
504  tmp.i = 0x5f3759d5 - (tmp.i >> 1);
505  x = tmp.f; // convert new bits into float
506  x = x*(1.5f - xhalf*x*x); // One round of Newton's method
507  return x;
508  }
509 
519  static inline floatmax_t powl(floatmax_t x, floatmax_t n)
520  {
521  //fall back to double precision pow if powl is not
522  //available
523 #ifdef HAVE_POWL
524  return ::powl((long double) x, (long double) n);
525 #else
526  return ::pow((double) x, (double) n);
527 #endif
528  }
529 
530  static inline int32_t pow(bool x, int32_t n)
531  {
532  return 0;
533  }
534 
539  static inline int32_t pow(int32_t x, int32_t n)
540  {
541  ASSERT(n>=0)
542  int32_t result=1;
543  while (n--)
544  result*=x;
545 
546  return result;
547  }
548 
553  static inline float64_t pow(float64_t x, int32_t n)
554  {
555  if (n>=0)
556  {
557  float64_t result=1;
558  while (n--)
559  result*=x;
560 
561  return result;
562  }
563  else
564  return ::pow((double)x, (double)n);
565  }
566 
571  static inline float64_t pow(float64_t x, float64_t n)
572  {
573  return ::pow((double) x, (double) n);
574  }
575 
580  static inline complex128_t pow(complex128_t x, int32_t n)
581  {
582  return std::pow(x, n);
583  }
584 
590  {
591  return std::pow(x, n);
592  }
593 
599  {
600  return std::pow(x, n);
601  }
602 
608  {
609  return std::pow(x, n);
610  }
612 
616  static inline float64_t exp(float64_t x)
617  {
618  return ::exp((double) x);
619  }
620 
622  static inline float64_t dot(const bool* v1, const bool* v2, int32_t n)
623  {
624  float64_t r=0;
625  for (int32_t i=0; i<n; i++)
626  r+=((v1[i]) ? 1 : 0) * ((v2[i]) ? 1 : 0);
627  return r;
628  }
629 
631  static inline floatmax_t dot(const floatmax_t* v1, const floatmax_t* v2, int32_t n)
632  {
633  floatmax_t r=0;
634  for (int32_t i=0; i<n; i++)
635  r+=v1[i]*v2[i];
636  return r;
637  }
638 
639 
641  static float64_t dot(const float64_t* v1, const float64_t* v2, int32_t n);
642 
644  static float32_t dot(const float32_t* v1, const float32_t* v2, int32_t n);
645 
647  static inline float64_t dot(
648  const uint64_t* v1, const uint64_t* v2, int32_t n)
649  {
650  float64_t r=0;
651  for (int32_t i=0; i<n; i++)
652  r+=((float64_t) v1[i])*v2[i];
653 
654  return r;
655  }
657  static inline float64_t dot(
658  const int64_t* v1, const int64_t* v2, int32_t n)
659  {
660  float64_t r=0;
661  for (int32_t i=0; i<n; i++)
662  r+=((float64_t) v1[i])*v2[i];
663 
664  return r;
665  }
666 
668  static inline float64_t dot(
669  const int32_t* v1, const int32_t* v2, int32_t n)
670  {
671  float64_t r=0;
672  for (int32_t i=0; i<n; i++)
673  r+=((float64_t) v1[i])*v2[i];
674 
675  return r;
676  }
677 
679  static inline float64_t dot(
680  const uint32_t* v1, const uint32_t* v2, int32_t n)
681  {
682  float64_t r=0;
683  for (int32_t i=0; i<n; i++)
684  r+=((float64_t) v1[i])*v2[i];
685 
686  return r;
687  }
688 
690  static inline float64_t dot(
691  const uint16_t* v1, const uint16_t* v2, int32_t n)
692  {
693  float64_t r=0;
694  for (int32_t i=0; i<n; i++)
695  r+=((float64_t) v1[i])*v2[i];
696 
697  return r;
698  }
699 
701  static inline float64_t dot(
702  const int16_t* v1, const int16_t* v2, int32_t n)
703  {
704  float64_t r=0;
705  for (int32_t i=0; i<n; i++)
706  r+=((float64_t) v1[i])*v2[i];
707 
708  return r;
709  }
710 
712  static inline float64_t dot(
713  const char* v1, const char* v2, int32_t n)
714  {
715  float64_t r=0;
716  for (int32_t i=0; i<n; i++)
717  r+=((float64_t) v1[i])*v2[i];
718 
719  return r;
720  }
721 
723  static inline float64_t dot(
724  const uint8_t* v1, const uint8_t* v2, int32_t n)
725  {
726  float64_t r=0;
727  for (int32_t i=0; i<n; i++)
728  r+=((float64_t) v1[i])*v2[i];
729 
730  return r;
731  }
732 
734  static inline float64_t dot(
735  const int8_t* v1, const int8_t* v2, int32_t n)
736  {
737  float64_t r=0;
738  for (int32_t i=0; i<n; i++)
739  r+=((float64_t) v1[i])*v2[i];
740 
741  return r;
742  }
743 
745  static inline float64_t dot(
746  const float64_t* v1, const char* v2, int32_t n)
747  {
748  float64_t r=0;
749  for (int32_t i=0; i<n; i++)
750  r+=((float64_t) v1[i])*v2[i];
751 
752  return r;
753  }
754 
757 
758 
766  static inline float64_t tan(float64_t x)
767  {
768  return ::tan((double) x);
769  }
770 
773 
774 
778  static inline float64_t atan(float64_t x)
779  {
780  return ::atan((double) x);
781  }
782 
785 
786 
791  static inline float64_t atan2(float64_t y, float64_t x)
792  {
793  return ::atan2((double) y, (double) x);
794  }
795 
798 
799 
803  static inline float64_t tanh(float64_t x)
804  {
805  return ::tanh((double) x);
806  }
807 
810 
811 
815  static inline float64_t sin(float64_t x)
816  {
817  return ::sin(x);
818  }
819 
822 
823 
827  static inline float64_t asin(float64_t x)
828  {
829  return ::asin(x);
830  }
831 
834 
835 
839  static inline float64_t sinh(float64_t x)
840  {
841  return ::sinh(x);
842  }
843 
846 
847 
851  static inline float64_t cos(float64_t x)
852  {
853  return ::cos(x);
854  }
855 
858 
859 
863  static inline float64_t acos(float64_t x)
864  {
865  return ::acos(x);
866  }
867 
870 
871 
875  static inline float64_t cosh(float64_t x)
876  {
877  return ::cosh(x);
878  }
879 
883 
892  static inline float64_t log10(float64_t v)
893  {
894  return ::log(v)/::log(10.0);
895  }
896 
899 
900 
904  static inline float64_t log2(float64_t v)
905  {
906 #ifdef HAVE_LOG2
907  return ::log2(v);
908 #else
909  return ::log(v)/::log(2.0);
910 #endif //HAVE_LOG2
911  }
912 
917  static inline float64_t log(float64_t v)
918  {
919  return ::log(v);
920  }
921 
924 
925  static inline index_t floor_log(index_t n)
926  {
927  index_t i;
928  for (i = 0; n != 0; i++)
929  n >>= 1;
930 
931  return i;
932  }
934 
941  static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
942  {
943  ASSERT(len>0 && xy)
944 
945  float64_t area = 0.0;
946 
947  if (!reversed)
948  {
949  for (int i=1; i<len; i++)
950  area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
951  }
952  else
953  {
954  for (int i=1; i<len; i++)
955  area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
956  }
957 
958  return area;
959  }
960 
966  static bool strtof(const char* str, float32_t* float_result);
967 
973  static bool strtod(const char* str, float64_t* double_result);
974 
980  static bool strtold(const char* str, floatmax_t* long_double_result);
981 
986  static inline int64_t factorial(int32_t n)
987  {
988  int64_t res=1;
989  for (int i=2; i<=n; i++)
990  res*=i ;
991  return res ;
992  }
993 
1001  static void init_random(uint32_t initseed=0)
1002  {
1003  if (initseed==0)
1005  else
1006  seed=initseed;
1007 
1008  sg_rand->set_seed(seed);
1009  }
1010 
1014  static inline uint64_t random()
1015  {
1016  return sg_rand->random_64();
1017  }
1018 
1022  static inline uint64_t random(uint64_t min_value, uint64_t max_value)
1023  {
1024  return sg_rand->random(min_value, max_value);
1025  }
1026 
1032  static inline int64_t random(int64_t min_value, int64_t max_value)
1033  {
1034  return sg_rand->random(min_value, max_value);
1035  }
1036 
1042  static inline uint32_t random(uint32_t min_value, uint32_t max_value)
1043  {
1044  return sg_rand->random(min_value, max_value);
1045  }
1046 
1052  static inline int32_t random(int32_t min_value, int32_t max_value)
1053  {
1054  return sg_rand->random(min_value, max_value);
1055  }
1056 
1062  static inline float32_t random(float32_t min_value, float32_t max_value)
1063  {
1064  return sg_rand->random(min_value, max_value);
1065  }
1066 
1072  static inline float64_t random(float64_t min_value, float64_t max_value)
1073  {
1074  return sg_rand->random(min_value, max_value);
1075  }
1076 
1082  static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
1083  {
1084  return sg_rand->random(min_value, max_value);
1085  }
1086 
1091  {
1092  // sets up variables & makes sure rand_s.range == (0,1)
1093  float32_t ret;
1094  float32_t rand_u;
1095  float32_t rand_v;
1096  float32_t rand_s;
1097  do
1098  {
1099  rand_u = static_cast<float32_t>(CMath::random(-1.0, 1.0));
1100  rand_v = static_cast<float32_t>(CMath::random(-1.0, 1.0));
1101  rand_s = rand_u*rand_u + rand_v*rand_v;
1102  } while ((rand_s == 0) || (rand_s >= 1));
1103 
1104  // the meat & potatos, and then the mean & standard deviation shifting...
1105  ret = static_cast<float32_t>(rand_u*CMath::sqrt(-2.0*CMath::log(rand_s)/rand_s));
1106  ret = std_dev*ret + mean;
1107  return ret;
1108  }
1109 
1114  {
1115  return sg_rand->normal_distrib(mean, std_dev);
1116  }
1117 
1120  static inline float32_t randn_float()
1121  {
1122  return static_cast<float32_t>(normal_random(0.0, 1.0));
1123  }
1124 
1127  static inline float64_t randn_double()
1128  {
1129  return sg_rand->std_normal_distrib();
1130  }
1132 
1140  static int32_t gcd(int32_t a, int32_t b)
1141  {
1142  REQUIRE((a>=0 && b>0) || (b>=0 && a>0), "gcd(%d,%d) is not defined.\n",
1143  a, b);
1144 
1145  if (1 == a || 1 == b)
1146  return 1;
1147 
1148  while (0 < a && 0 < b)
1149  {
1150  if (a > b)
1151  a %= b;
1152  else
1153  b %= a;
1154  }
1155 
1156  return 0 == a ? b : a;
1157  }
1158 
1164  template <class T>
1165  static void permute(SGVector<T> v, CRandom* rand=NULL)
1166  {
1167  if (rand)
1168  {
1169  for (index_t i=0; i<v.vlen; ++i)
1170  swap(v[i], v[rand->random(i, v.vlen-1)]);
1171  }
1172  else
1173  {
1174  for (index_t i=0; i<v.vlen; ++i)
1175  swap(v[i], v[random(i, v.vlen-1)]);
1176  }
1177  }
1178 
1184  template <class T>
1185  static int32_t get_num_nonzero(T* vec, int32_t len)
1186  {
1187  int32_t nnz = 0;
1188  for (index_t i=0; i<len; ++i)
1189  nnz += vec[i] != 0;
1190 
1191  return nnz;
1192  }
1193 
1199  static int32_t get_num_nonzero(complex128_t* vec, int32_t len)
1200  {
1201  int32_t nnz=0;
1202  for (index_t i=0; i<len; ++i)
1203  nnz+=vec[i]!=0.0;
1204  return nnz;
1205  }
1206 
1212  static inline int64_t nchoosek(int32_t n, int32_t k)
1213  {
1214  int64_t res=1;
1215 
1216  for (int32_t i=n-k+1; i<=n; i++)
1217  res*=i;
1218 
1219  return res/factorial(k);
1220  }
1221 
1228  static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100);
1229 
1236  template <class T>
1237  static float64_t* linspace(T start, T end, int32_t n)
1238  {
1239  float64_t* output = SG_MALLOC(float64_t, n);
1240  linspace(output, start, end, n);
1241 
1242  return output;
1243  }
1244 
1251  template <class T>
1252  static SGVector<float64_t> linspace_vec(T start, T end, int32_t n)
1253  {
1254  return SGVector<float64_t>(linspace(start, end, n), n);
1255  }
1256 
1262  template <class T>
1263  static T log_sum_exp(SGVector<T> values)
1264  {
1265  REQUIRE(values.vector, "Values are empty");
1266  REQUIRE(values.vlen>0,"number of values supplied is 0\n");
1267 
1268  /* find minimum element index */
1269  index_t min_index=0;
1270  T X0=values[0];
1271  if (values.vlen>1)
1272  {
1273  for (index_t i=1; i<values.vlen; ++i)
1274  {
1275  if (values[i]<X0)
1276  {
1277  X0=values[i];
1278  min_index=i;
1279  }
1280  }
1281  }
1282 
1283  /* remove element from vector copy and compute log sum exp */
1284  SGVector<T> values_without_X0(values.vlen-1);
1285  index_t from_idx=0;
1286  index_t to_idx=0;
1287  for (from_idx=0; from_idx<values.vlen; ++from_idx)
1288  {
1289  if (from_idx!=min_index)
1290  {
1291  values_without_X0[to_idx]=exp(values[from_idx]-X0);
1292  to_idx++;
1293  }
1294  }
1295 
1296  return X0+log(SGVector<T>::sum(values_without_X0)+1);
1297  }
1298 
1305  template <class T>
1306  static T log_mean_exp(SGVector<T> values)
1307  {
1308  return log_sum_exp(values) - log(values.vlen);
1309  }
1310 
1318  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
1319 
1326  static void sort(float64_t *a, int32_t*idx, int32_t N);
1327 
1333  template <class T>
1334  static void qsort(T* output, int32_t size)
1335  {
1336  if (size<=1)
1337  return;
1338 
1339  if (size==2)
1340  {
1341  if (output[0] > output [1])
1342  CMath::swap(output[0],output[1]);
1343  return;
1344  }
1345  //T split=output[random(0,size-1)];
1346  T split=output[size/2];
1347 
1348  int32_t left=0;
1349  int32_t right=size-1;
1350 
1351  while (left<=right)
1352  {
1353  while (output[left] < split)
1354  left++;
1355  while (output[right] > split)
1356  right--;
1357 
1358  if (left<=right)
1359  {
1360  CMath::swap(output[left],output[right]);
1361  left++;
1362  right--;
1363  }
1364  }
1365 
1366  if (right+1> 1)
1367  qsort(output,right+1);
1368 
1369  if (size-left> 1)
1370  qsort(&output[left],size-left);
1371  }
1372 
1378  template <class T>
1379  static void insertion_sort(T* output, int32_t size)
1380  {
1381  for (int32_t i=0; i<size-1; i++)
1382  {
1383  int32_t j=i-1;
1384  T value=output[i];
1385  while (j >= 0 && output[j] > value)
1386  {
1387  output[j+1] = output[j];
1388  j--;
1389  }
1390  output[j+1]=value;
1391  }
1392  }
1393 
1398  template <class T>
1399  inline static void radix_sort(T* array, int32_t size)
1400  {
1401  radix_sort_helper(array,size,0);
1402  }
1403 
1411  template <class T>
1412  static inline uint8_t byte(T word, uint16_t p)
1413  {
1414  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
1415  }
1416 
1418  static inline uint8_t byte(complex128_t word, uint16_t p)
1419  {
1420  SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
1421  return uint8_t(0);
1422  }
1423 
1429  template <class T>
1430  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
1431  {
1432  static size_t count[256], nc, cmin;
1433  T *ak;
1434  uint8_t c=0;
1435  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
1436  T *an, *aj, *pile[256];
1437  size_t *cp, cmax;
1438 
1439  /* Push initial array to stack */
1440  sp = s;
1441  radix_push(array, size, i);
1442 
1443  /* Loop until all digits have been sorted */
1444  while (sp>s) {
1445  radix_pop(array, size, i);
1446  an = array + size;
1447 
1448  /* Make character histogram */
1449  if (nc == 0) {
1450  cmin = 0xff;
1451  for (ak = array; ak < an; ak++) {
1452  c = byte(*ak, i);
1453  count[c]++;
1454  if (count[c] == 1) {
1455  /* Determine smallest character */
1456  if (c < cmin)
1457  cmin = c;
1458  nc++;
1459  }
1460  }
1461 
1462  /* Sort recursively if stack size too small */
1463  if (sp + nc > s + RADIX_STACK_SIZE) {
1464  radix_sort_helper(array, size, i);
1465  continue;
1466  }
1467  }
1468 
1469  /*
1470  * Set pile[]; push incompletely sorted bins onto stack.
1471  * pile[] = pointers to last out-of-place element in bins.
1472  * Before permuting: pile[c-1] + count[c] = pile[c];
1473  * during deal: pile[c] counts down to pile[c-1].
1474  */
1475  olds = bigs = sp;
1476  cmax = 2;
1477  ak = array;
1478  pile[0xff] = an;
1479  for (cp = count + cmin; nc > 0; cp++) {
1480  /* Find next non-empty pile */
1481  while (*cp == 0)
1482  cp++;
1483  /* Pile with several entries */
1484  if (*cp > 1) {
1485  /* Determine biggest pile */
1486  if (*cp > cmax) {
1487  cmax = *cp;
1488  bigs = sp;
1489  }
1490 
1491  if (i < sizeof(T)-1)
1492  radix_push(ak, *cp, (uint16_t) (i + 1));
1493  }
1494  pile[cp - count] = ak += *cp;
1495  nc--;
1496  }
1497 
1498  /* Play it safe -- biggest bin last. */
1499  swap(*olds, *bigs);
1500 
1501  /*
1502  * Permute misplacements home. Already home: everything
1503  * before aj, and in pile[c], items from pile[c] on.
1504  * Inner loop:
1505  * r = next element to put in place;
1506  * ak = pile[r[i]] = location to put the next element.
1507  * aj = bottom of 1st disordered bin.
1508  * Outer loop:
1509  * Once the 1st disordered bin is done, ie. aj >= ak,
1510  * aj<-aj + count[c] connects the bins in array linked list;
1511  * reset count[c].
1512  */
1513  aj = array;
1514  while (aj<an)
1515  {
1516  T r;
1517 
1518  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
1519  swap(*ak, r);
1520 
1521  *aj = r;
1522  aj += count[c];
1523  count[c] = 0;
1524  }
1525  }
1526  }
1527 
1529  static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
1530  {
1531  SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
1532  }
1533 
1540  template <class T>
1541  static void qsort(T** vector, index_t length)
1542  {
1543  if (length<=1)
1544  return;
1545 
1546  if (length==2)
1547  {
1548  if (*vector[0]>*vector[1])
1549  swap(vector[0],vector[1]);
1550  return;
1551  }
1552  T* split=vector[length/2];
1553 
1554  int32_t left=0;
1555  int32_t right=length-1;
1556 
1557  while (left<=right)
1558  {
1559  while (*vector[left]<*split)
1560  ++left;
1561  while (*vector[right]>*split)
1562  --right;
1563 
1564  if (left<=right)
1565  {
1566  swap(vector[left],vector[right]);
1567  ++left;
1568  --right;
1569  }
1570  }
1571 
1572  if (right+1>1)
1573  qsort(vector,right+1);
1574 
1575  if (length-left>1)
1576  qsort(&vector[left],length-left);
1577  }
1578 
1580  static void qsort(complex128_t** vector, index_t length)
1581  {
1582  SG_SERROR("CMath::qsort():: Not supported for complex128_t\n");
1583  }
1584 
1588  template <class T>
1589  static void qsort(SGVector<T> vector)
1590  {
1591  qsort<T>(vector, vector.size());
1592  }
1593 
1595  template<class T>
1597  {
1599  IndexSorter(const SGVector<T> *vec) { data = vec->vector; }
1600 
1602  bool operator() (index_t i, index_t j) const
1603  {
1604  return abs(data[i]-data[j])>std::numeric_limits<T>::epsilon()
1605  && data[i]<data[j];
1606  }
1607 
1609  const T* data;
1610  };
1611 
1619  template<class T>
1621  {
1622  IndexSorter<T> cmp(&vector);
1623  SGVector<index_t> idx(vector.size());
1624  for (index_t i=0; i < vector.size(); ++i)
1625  idx[i] = i;
1626 
1627  std::sort(idx.vector, idx.vector+vector.size(), cmp);
1628 
1629  return idx;
1630  }
1631 
1637  template <class T>
1638  static bool is_sorted(SGVector<T> vector)
1639  {
1640  if (vector.size() < 2)
1641  return true;
1642 
1643  for(int32_t i=1; i<vector.size(); i++)
1644  {
1645  if (vector[i-1] > vector[i])
1646  return false;
1647  }
1648 
1649  return true;
1650  }
1651 
1656  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1657  {
1658  ASSERT(width>=0)
1659  for (int i=0; i<width; i++)
1660  {
1661  T mask = ((T) 1)<<(sizeof(T)*8-1);
1662  while (mask)
1663  {
1664  if (mask & word)
1665  SG_SPRINT("1")
1666  else
1667  SG_SPRINT("0")
1668 
1669  mask>>=1;
1670  }
1671  }
1672  }
1673 
1675  static void display_bits(complex128_t word,
1676  int32_t width=8*sizeof(complex128_t))
1677  {
1678  SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
1679  }
1680 
1689  template <class T1,class T2>
1690  static void qsort_index(T1* output, T2* index, uint32_t size);
1691 
1693  template <class T>
1694  static void qsort_index(complex128_t* output, T* index, uint32_t size)
1695  {
1696  SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
1697  }
1698 
1707  template <class T1,class T2>
1708  static void qsort_backward_index(
1709  T1* output, T2* index, int32_t size);
1710 
1712  template <class T>
1714  complex128_t* output, T* index, uint32_t size)
1715  {
1716  SG_SERROR("CMath::qsort_backword_index():: \
1717  Not supported for complex128_t\n");
1718  }
1719 
1731  template <class T1,class T2>
1732  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1733  {
1734  int32_t n=0;
1735  thread_qsort<T1,T2> t;
1736  t.output=output;
1737  t.index=index;
1738  t.size=size;
1739  t.qsort_threads=&n;
1740  t.sort_limit=limit;
1741  t.num_threads=n_threads;
1742  parallel_qsort_index<T1,T2>(&t);
1743  }
1744 
1746  template <class T>
1747  inline static void parallel_qsort_index(complex128_t* output, T* index,
1748  uint32_t size, int32_t n_threads, int32_t limit=0)
1749  {
1750  SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1751  }
1752 
1754  template <class T1,class T2>
1755  static void* parallel_qsort_index(void* p);
1756 
1757 
1764  template <class T>
1765  static void min(float64_t* output, T* index, int32_t size);
1766 
1768  static void min(float64_t* output, complex128_t* index, int32_t size)
1769  {
1770  SG_SERROR("CMath::min():: Not supported for complex128_t\n");
1771  }
1772 
1776  template <class T>
1777  static void nmin(
1778  float64_t* output, T* index, int32_t size, int32_t n);
1779 
1781  static void nmin(float64_t* output, complex128_t* index,
1782  int32_t size, int32_t n)
1783  {
1784  SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
1785  }
1786 
1787 
1788 
1792  template <class T>
1793  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1794  {
1795  int32_t start=0;
1796  int32_t end=size-1;
1797 
1798  if (size<1)
1799  return 0;
1800 
1801  while (start<end)
1802  {
1803  int32_t middle=(start+end)/2;
1804 
1805  if (output[middle]>elem)
1806  end=middle-1;
1807  else if (output[middle]<elem)
1808  start=middle+1;
1809  else
1810  return middle;
1811  }
1812 
1813  return start;
1814  }
1815 
1817  static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
1818  {
1819  SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
1820  return int32_t(0);
1821  }
1822 
1829  template <class T>
1830  static inline int32_t binary_search(T* output, int32_t size, T elem)
1831  {
1832  int32_t ind = binary_search_helper(output, size, elem);
1833  if (ind >= 0 && output[ind] == elem)
1834  return ind;
1835  return -1;
1836  }
1837 
1839  static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
1840  {
1841  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1842  return int32_t(-1);
1843  }
1844 
1845 
1853  template<class T>
1854  static inline int32_t binary_search(T** vector, index_t length,
1855  T* elem)
1856  {
1857  int32_t start=0;
1858  int32_t end=length-1;
1859 
1860  if (length<1)
1861  return -1;
1862 
1863  while (start<end)
1864  {
1865  int32_t middle=(start+end)/2;
1866 
1867  if (*vector[middle]>*elem)
1868  end=middle-1;
1869  else if (*vector[middle]<*elem)
1870  start=middle+1;
1871  else
1872  {
1873  start=middle;
1874  break;
1875  }
1876  }
1877 
1878  if (start>=0&&*vector[start]==*elem)
1879  return start;
1880 
1881  return -1;
1882  }
1883 
1885  static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
1886  {
1887  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1888  return int32_t(-1);
1889  }
1890 
1897  template <class T>
1899  T* output, int32_t size, T elem)
1900  {
1901  int32_t ind = binary_search_helper(output, size, elem);
1902 
1903  if (output[ind]<=elem)
1904  return ind;
1905  if (ind>0 && output[ind-1] <= elem)
1906  return ind-1;
1907  return -1;
1908  }
1909 
1912  int32_t size, complex128_t elem)
1913  {
1914  SG_SERROR("CMath::binary_search_max_lower_equal():: \
1915  Not supported for complex128_t\n");
1916  return int32_t(-1);
1917  }
1918 
1927  static float64_t Align(
1928  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1929 
1930 
1932 
1938  {
1939  return c.real();
1940  }
1941 
1947  {
1948  return c.imag();
1949  }
1950 
1952  inline static uint32_t get_seed()
1953  {
1954  return CMath::seed;
1955  }
1956 
1958  inline static uint32_t get_log_range()
1959  {
1960  return CMath::LOGRANGE;
1961  }
1962 
1963 #ifdef USE_LOGCACHE
1964  inline static uint32_t get_log_accuracy()
1966  {
1967  return CMath::LOGACCURACY;
1968  }
1969 #endif
1970 
1972  static int is_finite(double f);
1973 
1975  static int is_infinity(double f);
1976 
1978  static int is_nan(double f);
1979 
1989 #ifdef USE_LOGCACHE
1990  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1991  {
1992  float64_t diff;
1993 
1994  if (!CMath::is_finite(p))
1995  return q;
1996 
1997  if (!CMath::is_finite(q))
1998  {
1999  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
2000  return NOT_A_NUMBER;
2001  }
2002  diff = p - q;
2003  if (diff > 0)
2004  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
2005  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
2006  }
2007 
2009  static void init_log_table();
2010 
2012  static int32_t determine_logrange();
2013 
2015  static int32_t determine_logaccuracy(int32_t range);
2016 #else
2018  float64_t p, float64_t q)
2019  {
2020  float64_t diff;
2021 
2022  if (!CMath::is_finite(p))
2023  return q;
2024  if (!CMath::is_finite(q))
2025  return p;
2026  diff = p - q;
2027  if (diff > 0)
2028  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
2029  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
2030  }
2031 #endif
2032 #ifdef USE_LOGSUMARRAY
2033 
2038  static inline float64_t logarithmic_sum_array(
2039  float64_t *p, int32_t len)
2040  {
2041  if (len<=2)
2042  {
2043  if (len==2)
2044  return logarithmic_sum(p[0],p[1]) ;
2045  if (len==1)
2046  return p[0];
2047  return -INFTY ;
2048  }
2049  else
2050  {
2051  float64_t *pp=p ;
2052  if (len%2==1) pp++ ;
2053  for (int32_t j=0; j < len>>1; j++)
2054  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
2055  }
2056  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
2057  }
2058 #endif //USE_LOGSUMARRAY
2059 
2060 
2062  virtual const char* get_name() const { return "Math"; }
2063  public:
2066  static const float64_t NOT_A_NUMBER;
2069  static const float64_t INFTY;
2070  static const float64_t ALMOST_INFTY;
2071 
2074 
2076  static const float64_t PI;
2077 
2080 
2081  /* Largest and smallest possible float64_t */
2084 
2085  /* Floating point Limits, Normalized */
2086  static const float32_t F_MAX_VAL32;
2088  static const float64_t F_MAX_VAL64;
2090 
2091  /* Floating point limits, Denormalized */
2092  static const float32_t F_MIN_VAL32;
2093  static const float64_t F_MIN_VAL64;
2094 
2095  protected:
2097  static int32_t LOGRANGE;
2098 
2100  static uint32_t seed;
2101 
2102 #ifdef USE_LOGCACHE
2103 
2105  static int32_t LOGACCURACY;
2107  static float64_t* logtable;
2109 #endif
2110 };
2111 
2112 //implementations of template functions
2113 template <class T1,class T2>
2115  {
2116  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
2117  T1* output=ps->output;
2118  T2* index=ps->index;
2119  int32_t size=ps->size;
2120  int32_t* qsort_threads=ps->qsort_threads;
2121  int32_t sort_limit=ps->sort_limit;
2122  int32_t num_threads=ps->num_threads;
2123 
2124  if (size<2)
2125  {
2126  return NULL;
2127  }
2128 
2129  if (size==2)
2130  {
2131  if (output[0] > output [1])
2132  {
2133  swap(output[0], output[1]);
2134  swap(index[0], index[1]);
2135  }
2136  return NULL;
2137  }
2138  //T1 split=output[random(0,size-1)];
2139  T1 split=output[size/2];
2140 
2141  int32_t left=0;
2142  int32_t right=size-1;
2143 
2144  while (left<=right)
2145  {
2146  while (output[left] < split)
2147  left++;
2148  while (output[right] > split)
2149  right--;
2150 
2151  if (left<=right)
2152  {
2153  swap(output[left], output[right]);
2154  swap(index[left], index[right]);
2155  left++;
2156  right--;
2157  }
2158  }
2159  bool lthread_start=false;
2160  bool rthread_start=false;
2161  pthread_t lthread;
2162  pthread_t rthread;
2163  struct thread_qsort<T1,T2> t1;
2164  struct thread_qsort<T1,T2> t2;
2165 
2166  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
2167  qsort_index(output,index,right+1);
2168  else if (right+1> 1)
2169  {
2170  (*qsort_threads)++;
2171  lthread_start=true;
2172  t1.output=output;
2173  t1.index=index;
2174  t1.size=right+1;
2175  t1.qsort_threads=qsort_threads;
2176  t1.sort_limit=sort_limit;
2177  t1.num_threads=num_threads;
2178  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
2179  {
2180  lthread_start=false;
2181  (*qsort_threads)--;
2182  qsort_index(output,index,right+1);
2183  }
2184  }
2185 
2186 
2187  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
2188  qsort_index(&output[left],&index[left], size-left);
2189  else if (size-left> 1)
2190  {
2191  (*qsort_threads)++;
2192  rthread_start=true;
2193  t2.output=&output[left];
2194  t2.index=&index[left];
2195  t2.size=size-left;
2196  t2.qsort_threads=qsort_threads;
2197  t2.sort_limit=sort_limit;
2198  t2.num_threads=num_threads;
2199  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
2200  {
2201  rthread_start=false;
2202  (*qsort_threads)--;
2203  qsort_index(&output[left],&index[left], size-left);
2204  }
2205  }
2206 
2207  if (lthread_start)
2208  {
2209  pthread_join(lthread, NULL);
2210  (*qsort_threads)--;
2211  }
2212 
2213  if (rthread_start)
2214  {
2215  pthread_join(rthread, NULL);
2216  (*qsort_threads)--;
2217  }
2218 
2219  return NULL;
2220  }
2221 
2222  template <class T1,class T2>
2223 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
2224 {
2225  if (size<=1)
2226  return;
2227 
2228  if (size==2)
2229  {
2230  if (output[0] > output [1])
2231  {
2232  swap(output[0],output[1]);
2233  swap(index[0],index[1]);
2234  }
2235  return;
2236  }
2237  //T1 split=output[random(0,size-1)];
2238  T1 split=output[size/2];
2239 
2240  int32_t left=0;
2241  int32_t right=size-1;
2242 
2243  while (left<=right)
2244  {
2245  while (output[left] < split)
2246  left++;
2247  while (output[right] > split)
2248  right--;
2249 
2250  if (left<=right)
2251  {
2252  swap(output[left],output[right]);
2253  swap(index[left],index[right]);
2254  left++;
2255  right--;
2256  }
2257  }
2258 
2259  if (right+1> 1)
2260  qsort_index(output,index,right+1);
2261 
2262  if (size-left> 1)
2263  qsort_index(&output[left],&index[left], size-left);
2264 }
2265 
2266  template <class T1,class T2>
2267 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
2268 {
2269  if (size<=1)
2270  return;
2271 
2272  if (size==2)
2273  {
2274  if (output[0] < output [1])
2275  {
2276  swap(output[0],output[1]);
2277  swap(index[0],index[1]);
2278  }
2279  return;
2280  }
2281 
2282  //T1 split=output[random(0,size-1)];
2283  T1 split=output[size/2];
2284 
2285  int32_t left=0;
2286  int32_t right=size-1;
2287 
2288  while (left<=right)
2289  {
2290  while (output[left] > split)
2291  left++;
2292  while (output[right] < split)
2293  right--;
2294 
2295  if (left<=right)
2296  {
2297  swap(output[left],output[right]);
2298  swap(index[left],index[right]);
2299  left++;
2300  right--;
2301  }
2302  }
2303 
2304  if (right+1> 1)
2305  qsort_backward_index(output,index,right+1);
2306 
2307  if (size-left> 1)
2308  qsort_backward_index(&output[left],&index[left], size-left);
2309 }
2310 
2311  template <class T>
2312 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
2313 {
2314  if (6*n*size<13*size*CMath::log(size))
2315  for (int32_t i=0; i<n; i++)
2316  min(&output[i], &index[i], size-i);
2317  else
2318  qsort_index(output, index, size);
2319 }
2320 
2321 /* move the smallest entry in the array to the beginning */
2322  template <class T>
2323 void CMath::min(float64_t* output, T* index, int32_t size)
2324 {
2325  if (size<=1)
2326  return;
2327  float64_t min_elem=output[0];
2328  int32_t min_index=0;
2329  for (int32_t i=1; i<size; i++)
2330  {
2331  if (output[i]<min_elem)
2332  {
2333  min_index=i;
2334  min_elem=output[i];
2335  }
2336  }
2337  swap(output[0], output[min_index]);
2338  swap(index[0], index[min_index]);
2339 }
2340 
2342 template <>
2343 inline float64_t* CMath::linspace<complex128_t>(complex128_t start, complex128_t end, int32_t n)
2344 {
2345  SG_SERROR("SGVector::linspace():: Not supported for complex128_t\n");
2346  return NULL;
2347 }
2348 
2349 #define COMPLEX128_ERROR_ONEVECARG_RETURNS_T(function, return_type, return_statement) \
2350 template <> \
2351 inline return_type CMath::function<complex128_t>(SGVector<complex128_t> vector) \
2352 { \
2353  SG_SERROR("CMath::%s():: Not supported for complex128_t\n", \
2354  #function); \
2355  return_statement; \
2356 }
2357 
2358 #define COMPLEX128_ERROR_ONEARG_T(function) \
2359 template <> \
2360 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
2361 { \
2362  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2363  #function);\
2364  return complex128_t(0.0, 0.0); \
2365 }
2366 
2367 #define COMPLEX128_ERROR_TWOARGS_T(function) \
2368 template <> \
2369 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \
2370 { \
2371  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2372  #function);\
2373  return complex128_t(0.0, 0.0); \
2374 }
2375 
2376 #define COMPLEX128_ERROR_THREEARGS_T(function) \
2377 template <> \
2378 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \
2379 { \
2380  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2381  #function);\
2382  return complex128_t(0.0, 0.0); \
2383 }
2384 
2385 #define COMPLEX128_ERROR_SORT_T(function) \
2386 template <> \
2387 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \
2388 { \
2389  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2390  #function);\
2391 }
2392 
2393 #define COMPLEX128_ERROR_ARG_MAX_MIN(function) \
2394 template <> \
2395 inline int32_t CMath::function<complex128_t>(complex128_t * a, int32_t b, int32_t c, complex128_t * d) \
2396 { \
2397  int32_t maxIdx=0; \
2398  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2399  #function);\
2400  return maxIdx; \
2401 }
2402 
2405 
2407 COMPLEX128_ERROR_ONEVECARG_RETURNS_T(argsort, SGVector<index_t>, SGVector<index_t> idx(vector.size());return idx;)
2408 
2410 COMPLEX128_ERROR_ONEVECARG_RETURNS_T(is_sorted, bool, return false;)
2411 
2414 
2417 
2420 
2422 // COMPLEX128_ERROR_ONEARG_T(sign)
2423 
2426 
2429 
2432 
2435 
2438 
2439 }
2440 #undef COMPLEX128_ERROR_ONEARG
2441 #undef COMPLEX128_ERROR_ONEARG_T
2442 #undef COMPLEX128_ERROR_TWOARGS_T
2443 #undef COMPLEX128_ERROR_THREEARGS_T
2444 #undef COMPLEX128_STDMATH
2445 #undef COMPLEX128_ERROR_SORT_T
2446 #endif
static float64_t sin(float64_t x)
tanh(x), x being a complex128_t
Definition: Math.h:815
static float64_t normal_random(float64_t mean, float64_t std_dev)
Definition: Math.h:1113
static const float32_t F_MAX_VAL32
Definition: Math.h:2086
static const float64_t MACHINE_EPSILON
Definition: Math.h:2079
static bool strtof(const char *str, float32_t *float_result)
Definition: Math.cpp:274
static int32_t gcd(int32_t a, int32_t b)
Definition: Math.h:1140
static void permute(SGVector< T > v, CRandom *rand=NULL)
Definition: Math.h:1165
static int32_t binary_search(complex128_t *output, int32_t size, complex128_t elem)
binary_search not implemented for complex128_t
Definition: Math.h:1839
static void parallel_qsort_index(T1 *output, T2 *index, uint32_t size, int32_t n_threads, int32_t limit=262144)
Definition: Math.h:1732
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:631
static uint32_t seed
random generator seed
Definition: Math.h:2100
std::complex< float64_t > complex128_t
Definition: common.h:77
static int is_finite(double f)
checks whether a float is finite
Definition: Math.cpp:259
static float64_t Align(char *seq1, char *seq2, int32_t l1, int32_t l2, float64_t gapCost)
Definition: Math.cpp:172
static int32_t arg_max(T *vec, int32_t inc, int32_t len, T *maxv_ptr=NULL)
Definition: Math.h:258
static float64_t sqrt(float64_t x)
Definition: Math.h:463
static floatmax_t random(floatmax_t min_value, floatmax_t max_value)
Definition: Math.h:1082
static floatmax_t sqrt(floatmax_t x)
Definition: Math.h:472
static void linspace(float64_t *output, float64_t start, float64_t end, int32_t n=100)
Definition: Math.cpp:215
int32_t index_t
Definition: common.h:72
static float64_t ceil(float64_t d)
Definition: Math.h:411
static bool strtod(const char *str, float64_t *double_result)
Definition: Math.cpp:305
static complex128_t pow(complex128_t x, int32_t n)
Definition: Math.h:580
virtual ~CMath()
Destructor - frees logtable.
Definition: Math.cpp:76
static const float64_t INFTY
infinity
Definition: Math.h:2069
static SGVector< float64_t > linspace_vec(T start, T end, int32_t n)
Definition: Math.h:1252
static int32_t binary_search_helper(T *output, int32_t size, T elem)
Definition: Math.h:1793
#define COMPLEX128_ERROR_THREEARGS_T(function)
Definition: Math.h:2376
static void nmin(float64_t *output, T *index, int32_t size, int32_t n)
Definition: Math.h:2312
static void qsort_index(T1 *output, T2 *index, uint32_t size)
Definition: Math.h:2223
static void qsort_index(complex128_t *output, T *index, uint32_t size)
qsort_index not implemented for complex128_t
Definition: Math.h:1694
static float64_t log10(float64_t v)
Definition: Math.h:892
#define COMPLEX128_ERROR_TWOARGS_T(function)
Definition: Math.h:2367
#define SG_SWARNING(...)
Definition: SGIO.h:177
static void qsort(SGVector< T > vector)
Definition: Math.h:1589
#define COMPLEX128_ERROR_ONEVECARG_RETURNS_T(function, return_type, return_statement)
Definition: Math.h:2349
static float32_t normal_random(float32_t mean, float32_t std_dev)
Definition: Math.h:1090
static float64_t random(float64_t min_value, float64_t max_value)
Definition: Math.h:1072
static int32_t binary_search(T *output, int32_t size, T elem)
Definition: Math.h:1830
static T sq(T x)
Definition: Math.h:445
static uint32_t get_seed()
returns number generator seed
Definition: Math.h:1952
#define COMPLEX128_ERROR_ARG_MAX_MIN(function)
Definition: Math.h:2393
static float32_t randn_float()
Definition: Math.h:1120
static float64_t atan2(float64_t y, float64_t x)
atan(x), x being a complex128_t not implemented
Definition: Math.h:791
static const float64_t MIN_REAL_NUMBER
Definition: Math.h:2083
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:712
static float64_t randn_double()
Definition: Math.h:1127
static int32_t binary_search_max_lower_equal(T *output, int32_t size, T elem)
Definition: Math.h:1898
static float64_t atan(float64_t x)
tan(x), x being a complex128_t
Definition: Math.h:778
#define REQUIRE(x,...)
Definition: SGIO.h:205
static const float64_t F_MAX_VAL64
Definition: Math.h:2088
static float64_t cosh(float64_t x)
acos(x), x being a complex128_t not implemented
Definition: Math.h:875
static uint32_t get_log_range()
returns range of logtable
Definition: Math.h:1958
void split(v_array< ds_node< P > > &point_set, v_array< ds_node< P > > &far_set, int max_scale)
Definition: JLCoverTree.h:154
static float32_t random(float32_t min_value, float32_t max_value)
Definition: Math.h:1062
static const float32_t F_MIN_VAL32
Definition: Math.h:2092
CRandom * sg_rand
Definition: init.cpp:41
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:734
std::enable_if<!std::is_same< T, complex128_t >::value, float64_t >::type mean(const Container< T > &a)
IndexSorter(const SGVector< T > *vec)
Definition: Math.h:1599
static bool fequals_abs(const T &a, const T &b, const float64_t eps)
Definition: Math.h:314
#define RADIX_STACK_SIZE
Definition: Math.h:73
static float64_t imag(complex128_t c)
Definition: Math.h:1946
static float64_t floor(float64_t d)
Definition: Math.h:402
static uint64_t random()
Definition: Math.h:1014
static float64_t * linspace(T start, T end, int32_t n)
Definition: Math.h:1237
static int32_t get_num_nonzero(complex128_t *vec, int32_t len)
Definition: Math.h:1199
static const float32_t F_MIN_NORM_VAL32
Definition: Math.h:2087
static float64_t real(complex128_t c)
Definition: Math.h:1937
static uint8_t byte(complex128_t word, uint16_t p)
byte not implemented for complex128_t
Definition: Math.h:1418
static void qsort(T *output, int32_t size)
Definition: Math.h:1334
static int32_t LOGRANGE
range for logtable: log(1+exp(x)) -LOGRANGE <= x <= 0
Definition: Math.h:2097
static const float64_t ALMOST_NEG_INFTY
almost neg (log) infinity
Definition: Math.h:2073
static bool fequals(const T &a, const T &b, const float64_t eps, bool tolerant=false)
Definition: Math.h:331
static complex128_t pow(complex128_t x, complex128_t n)
Definition: Math.h:589
static float32_t invsqrt(float32_t x)
x^0.5, x being a complex128_t
Definition: Math.h:490
int32_t size() const
Definition: SGVector.h:136
index_t vlen
Definition: SGVector.h:545
static uint8_t byte(T word, uint16_t p)
Definition: Math.h:1412
static complex128_t pow(float64_t x, complex128_t n)
Definition: Math.h:607
#define SG_SPRINT(...)
Definition: SGIO.h:179
CMath()
Constructor - initializes log-table.
Definition: Math.cpp:59
static float64_t pow(float64_t x, int32_t n)
Definition: Math.h:553
#define ASSERT(x)
Definition: SGIO.h:200
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:125
static int32_t random(int32_t min_value, int32_t max_value)
Definition: Math.h:1052
static float64_t pow(float64_t x, float64_t n)
Definition: Math.h:571
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:723
bool operator()(index_t i, index_t j) const
Definition: Math.h:1602
static void qsort(T **vector, index_t length)
Definition: Math.h:1541
static void init_random(uint32_t initseed=0)
Definition: Math.h:1001
static int32_t pow(int32_t x, int32_t n)
Definition: Math.h:539
#define radix_pop(a, n, i)
Definition: Math.h:77
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:690
double float64_t
Definition: common.h:60
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:1747
static void min(float64_t *output, complex128_t *index, int32_t size)
complex128_t cannot be used as index
Definition: Math.h:1768
long double floatmax_t
Definition: common.h:61
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:701
#define COMPLEX128_ERROR_ONEARG(function)
Definition: Math.h:110
static SGVector< index_t > argsort(SGVector< T > vector)
Definition: Math.h:1620
static float64_t cos(float64_t x)
sinh(x), x being a complex128_t
Definition: Math.h:851
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:679
static void qsort_backword_index(complex128_t *output, T *index, uint32_t size)
qsort_backword_index not implemented for complex128_t
Definition: Math.h:1713
static int32_t get_num_nonzero(T *vec, int32_t len)
Definition: Math.h:1185
static T log_sum_exp(SGVector< T > values)
Definition: Math.h:1263
static T max(T a, T b)
Definition: Math.h:164
static float64_t area_under_curve(float64_t *xy, int32_t len, bool reversed)
Definition: Math.h:941
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:622
static void display_bits(T word, int32_t width=8 *sizeof(T))
Definition: Math.h:1656
static uint64_t random(uint64_t min_value, uint64_t max_value)
Definition: Math.h:1022
static int64_t factorial(int32_t n)
Definition: Math.h:986
static void qsort(complex128_t **vector, index_t length)
qsort not implemented for complex128_t
Definition: Math.h:1580
static int32_t arg_min(T *vec, int32_t inc, int32_t len, T *minv_ptr=NULL)
Definition: Math.h:285
static float64_t tan(float64_t x)
Definition: Math.h:766
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:1911
float float32_t
Definition: common.h:59
static float64_t sinh(float64_t x)
asin(x), x being a complex128_t not implemented
Definition: Math.h:839
: 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:374
static int64_t nchoosek(int32_t n, int32_t k)
Definition: Math.h:1212
static T min(T *vec, int32_t len)
Definition: Math.h:205
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:668
static int32_t binary_search(T **vector, index_t length, T *elem)
Definition: Math.h:1854
static int is_infinity(double f)
checks whether a float is infinity
Definition: Math.cpp:241
static float64_t acos(float64_t x)
cos(x), x being a complex128_t
Definition: Math.h:863
static float64_t abs(complex128_t a)
Definition: Math.h:191
static float64_t log2(float64_t v)
log10(x), x being a complex128_t
Definition: Math.h:904
static void display_bits(complex128_t word, int32_t width=8 *sizeof(complex128_t))
disply_bits not implemented for complex128_t
Definition: Math.h:1675
#define COMPLEX128_STDMATH(function)
Definition: Math.h:118
static void radix_sort_helper(T *array, int32_t size, uint16_t i)
Definition: Math.h:1430
static T sign(T a)
Definition: Math.h:421
static float64_t tanh(float64_t x)
atan2(x), x being a complex128_t not implemented
Definition: Math.h:803
static int is_nan(double f)
checks whether a float is nan
Definition: Math.cpp:228
static float64_t dot(const float64_t *v1, const char *v2, int32_t n)
Compute dot product between v1 and v2.
Definition: Math.h:745
virtual const char * get_name() const
Definition: Math.h:2062
static float64_t asin(float64_t x)
sin(x), x being a complex128_t
Definition: Math.h:827
static T min(T a, T b)
Definition: Math.h:153
#define SG_SERROR(...)
Definition: SGIO.h:178
static float64_t exp(float64_t x)
Definition: Math.h:616
static const float64_t F_MIN_VAL64
Definition: Math.h:2093
static const float64_t F_MIN_NORM_VAL64
Definition: Math.h:2089
static float64_t log(float64_t v)
Definition: Math.h:917
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:1529
static const float64_t ALMOST_INFTY
Definition: Math.h:2070
Class which collects generic mathematical functions.
Definition: Math.h:130
static T log_mean_exp(SGVector< T > values)
Definition: Math.h:1306
static void insertion_sort(T *output, int32_t size)
Definition: Math.h:1379
static void swap(T &a, T &b)
Definition: Math.h:433
static complex128_t pow(complex128_t x, float64_t n)
Definition: Math.h:598
static void sort(int32_t *a, int32_t cols, int32_t sort_col=0)
Definition: Math.cpp:133
static T max(T *vec, int32_t len)
Definition: Math.h:222
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:1817
static float64_t round(float64_t d)
Definition: Math.h:393
static bool is_sorted(SGVector< T > vector)
Definition: Math.h:1638
static float32_t sqrt(float32_t x)
Definition: Math.h:454
static uint32_t generate_seed()
Definition: Random.cpp:346
static index_t floor_log(index_t n)
log(x), x being a complex128_t
Definition: Math.h:925
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:657
static void radix_sort(T *array, int32_t size)
Definition: Math.h:1399
static uint32_t random(uint32_t min_value, uint32_t max_value)
Definition: Math.h:1042
static floatmax_t powl(floatmax_t x, floatmax_t n)
Definition: Math.h:519
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:647
static float64_t logarithmic_sum(float64_t p, float64_t q)
Definition: Math.h:2017
#define radix_push(a, n, i)
Definition: Math.h:76
static T clamp(T value, T lb, T ub)
Definition: Math.h:240
static int32_t binary_search(complex128_t **vector, index_t length, complex128_t *elem)
binary_search not implemented for complex128_t
Definition: Math.h:1885
#define COMPLEX128_ERROR_SORT_T(function)
Definition: Math.h:2385
static bool strtold(const char *str, floatmax_t *long_double_result)
Definition: Math.cpp:336
static int32_t pow(bool x, int32_t n)
Definition: Math.h:530
static const float64_t NOT_A_NUMBER
not a number
Definition: Math.h:2067
static void qsort_backward_index(T1 *output, T2 *index, int32_t size)
Definition: Math.h:2267
static const float64_t MAX_REAL_NUMBER
Definition: Math.h:2082
static T abs(T a)
Definition: Math.h:175
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:1781
static int64_t random(int64_t min_value, int64_t max_value)
Definition: Math.h:1032
static const float64_t PI
Definition: Math.h:2076

SHOGUN Machine Learning Toolbox - Documentation