SHOGUN  6.1.3
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 <cmath>
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 
49 #ifndef M_PI
50 #define M_PI 3.14159265358979323846
51 #endif
52 
53 /* Size of RNG seed */
54 #define RNG_SEED_SIZE 256
55 
56 /* Maximum stack size */
57 #define RADIX_STACK_SIZE 512
58 
59 /* Stack macros */
60 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
61 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
62 
63 #ifndef DOXYGEN_SHOULD_SKIP_THIS
64 
65 template <class T> struct radix_stack_t
66 {
68  T *sa;
70  size_t sn;
72  uint16_t si;
73 };
74 
76 template <class T1, class T2> struct thread_qsort
77 {
79  T1* output;
81  T2* index;
83  uint32_t size;
84 
86  int32_t* qsort_threads;
88  int32_t sort_limit;
90  int32_t num_threads;
91 };
92 #endif // DOXYGEN_SHOULD_SKIP_THIS
93 
94 #define COMPLEX128_ERROR_ONEARG(function) \
95 static inline complex128_t function(complex128_t a) \
96 { \
97  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
98  #function);\
99  return complex128_t(0.0, 0.0); \
100 }
101 
102 #define COMPLEX128_STDMATH(function) \
103 static inline complex128_t function(complex128_t a) \
104 { \
105  return std::function(a); \
106 }
107 
108 namespace shogun
109 {
111  extern CRandom* sg_rand;
114 class CMath : public CSGObject
115 {
116  public:
120  CMath();
122 
124  virtual ~CMath();
126 
127 #ifndef SWIG // SWIG should skip this part
128 
131 
137  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
138  static inline T min(T a, T b)
139  {
140  return std::min(a, b);
141  }
142 
148  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
149  static inline T max(T a, T b)
150  {
151  return std::max(a, b);
152  }
153 #endif
154 
160  template <class T>
161  static inline T abs(T a)
162  {
163  // can't be a>=0?(a):(-a), because compiler complains about
164  // 'comparison always true' when T is unsigned
165  if (a==0)
166  return 0;
167  else if (a>0)
168  return a;
169  else
170  return -a;
171  }
172 
177  static inline float64_t abs(complex128_t a)
178  {
179  float64_t a_real=a.real();
180  float64_t a_imag=a.imag();
181  return (CMath::sqrt(a_real*a_real+a_imag*a_imag));
182  }
184 
190  template <class T>
191  static T min(T* vec, int32_t len)
192  {
193  ASSERT(len>0)
194  return *std::min_element(vec, vec+len);
195  }
196 
202  template <class T>
203  static T max(T* vec, int32_t len)
204  {
205  ASSERT(len>0)
206  return *std::max_element(vec, vec+len);
207  }
208 
209 #ifndef SWIG // SWIG should skip this part
210 
216  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
217  static inline T clamp(T value, T lb, T ub)
218  {
219  if (value<=lb)
220  return lb;
221  else if (value>=ub)
222  return ub;
223  else
224  return value;
225  }
226 
234  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
235  static int32_t arg_max(T * vec, int32_t inc, int32_t len, T * maxv_ptr = NULL)
236  {
237  ASSERT(len > 0 || inc > 0)
238 
239  T maxv = vec[0];
240  int32_t maxIdx = 0;
241 
242  for (int32_t i = 1, j = inc ; i < len ; i++, j += inc)
243  {
244  if (vec[j] > maxv)
245  maxv = vec[j], maxIdx = i;
246  }
247 
248  if (maxv_ptr != NULL)
249  *maxv_ptr = maxv;
250 
251  return maxIdx;
252  }
253 
261  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
262  static int32_t arg_min(T * vec, int32_t inc, int32_t len, T * minv_ptr = NULL)
263  {
264  ASSERT(len > 0 || inc > 0)
265 
266  T minv = vec[0];
267  int32_t minIdx = 0;
268 
269  for (int32_t i = 1, j = inc ; i < len ; i++, j += inc)
270  {
271  if (vec[j] < minv)
272  minv = vec[j], minIdx = i;
273  }
274 
275  if (minv_ptr != NULL)
276  *minv_ptr = minv;
277 
278  return minIdx;
279  }
280 
283 
290  template <class T, class = typename std::enable_if<std::is_floating_point<T>::value>::type>
291  static inline bool fequals_abs(const T& a, const T& b,
292  const float64_t eps)
293  {
294  const T diff = CMath::abs<T>((a-b));
295  return (diff < eps);
296  }
297 
307  template <class T, class = typename std::enable_if<std::is_floating_point<T>::value>::type>
308  static inline bool fequals(const T& a, const T& b,
309  const float64_t eps, bool tolerant=false)
310  {
311  const T absA = CMath::abs<T>(a);
312  const T absB = CMath::abs<T>(b);
313  const T diff = CMath::abs<T>((a-b));
314 
315  // Handle this separately since NAN is unordered
316  if (CMath::is_nan((float64_t)a) && CMath::is_nan((float64_t)b))
317  return true;
318 
319  // Required for JSON Serialization Tests
320  if (tolerant)
321  return CMath::fequals_abs<T>(a, b, eps);
322 
323  // handles float32_t and float64_t separately
324  T comp = (std::is_same<float32_t, T>::value) ? CMath::F_MIN_NORM_VAL32 : CMath::F_MIN_NORM_VAL64;
325 
326  if (a == b)
327  return true;
328 
329  // both a and b are 0 and relative error is less meaningful
330  else if ((a == 0) || (b == 0) || (diff < comp))
331  return (diff < (eps * comp));
332  // use max(relative error, diff) to handle large eps
333  else
334  {
335  T check = ((diff/(absA + absB)) > diff)?
336  (diff/(absA + absB)):diff;
337  return (check < eps);
338  }
339  }
340 #endif
341 
342  /* Get the corresponding absolute tolerance for unit test given a relative tolerance
343  *
344  * Note that a unit test will be passed only when
345  * \f[
346  * |v_\text{true} - v_\text{predict}| \leq tol_\text{relative} * |v_\text{true}|
347  * \f] which is equivalent to
348  * \f[
349  * |v_\text{true} - v_\text{predict}| \leq tol_\text{absolute}
350  * \f] where
351  * \f[
352  * tol_\text{absolute} = tol_\text{relative} * |v_\text{true}|
353  * \f]
354  *
355  * @param true_value true value should be finite (neither NAN nor INF)
356  * @param rel_tolerance relative tolerance should be positive and less than 1.0
357  *
358  * @return the corresponding absolute tolerance
359  */
360  static float64_t get_abs_tolerance(float64_t true_value, float64_t rel_tolerance);
361 
366  static inline float64_t round(float64_t d)
367  {
368  return std::floor(d+0.5);
369  }
370 
375  static inline float64_t floor(float64_t d)
376  {
377  return std::floor(d);
378  }
379 
384  static inline float64_t ceil(float64_t d)
385  {
386  return std::ceil(d);
387  }
388 
393  template <class T>
394  static inline T sign(T a)
395  {
396  if (a==0)
397  return 0;
398  else return (a<0) ? (-1) : (+1);
399  }
400 
405  template <class T>
406  static inline void swap(T &a,T &b)
407  {
408  T c=a;
409  a=b;
410  b=c;
411  }
412 
417  template <class T>
418  static inline T sq(T x)
419  {
420  return x*x;
421  }
422 
427  template <class T>
428  static inline T sqrt(T x)
429  {
430  return std::sqrt(x);
431  }
432 
434  COMPLEX128_STDMATH(sqrt)
435 
436 
440  static inline float32_t invsqrt(float32_t x)
441  {
442  union float_to_int
443  {
444  float32_t f;
445  int32_t i;
446  };
447 
448  float_to_int tmp;
449  tmp.f=x;
450 
451  float32_t xhalf = 0.5f * x;
452  // store floating-point bits in integer tmp.i
453  // and do initial guess for Newton's method
454  tmp.i = 0x5f3759d5 - (tmp.i >> 1);
455  x = tmp.f; // convert new bits into float
456  x = x*(1.5f - xhalf*x*x); // One round of Newton's method
457  return x;
458  }
459 
469  static inline floatmax_t powl(floatmax_t x, floatmax_t n)
470  {
471  return std::pow(x, n);
472  }
473 
474  static inline int32_t pow(bool x, int32_t n)
475  {
476  return 0;
477  }
478 
483  static inline int32_t pow(int32_t x, int32_t n)
484  {
485  ASSERT(n>=0)
486  int32_t result=1;
487  while (n--)
488  result*=x;
489 
490  return result;
491  }
492 
497  static inline float64_t pow(float64_t x, int32_t n)
498  {
499  return std::pow(x, n);
500  }
501 
506  static inline float64_t pow(float64_t x, float64_t n)
507  {
508  return std::pow(x, n);
509  }
510 
515  static inline complex128_t pow(complex128_t x, int32_t n)
516  {
517  return std::pow(x, n);
518  }
519 
525  {
526  return std::pow(x, n);
527  }
528 
534  {
535  return std::pow(x, n);
536  }
537 
543  {
544  return std::pow(x, n);
545  }
547 
551  static inline float64_t exp(float64_t x)
552  {
553  return std::exp(x);
554  }
555 
557  COMPLEX128_STDMATH(exp)
558 
559 
567  static inline float64_t tan(float64_t x)
568  {
569  return std::tan(x);
570  }
571 
573  COMPLEX128_STDMATH(tan)
574 
575 
579  static inline float64_t atan(float64_t x)
580  {
581  return std::atan(x);
582  }
583 
586 
587 
592  static inline float64_t atan2(float64_t y, float64_t x)
593  {
594  return std::atan2(y, x);
595  }
596 
599 
600 
604  static inline float64_t tanh(float64_t x)
605  {
606  return std::tanh(x);
607  }
608 
610  COMPLEX128_STDMATH(tanh)
611 
612 
616  static inline float64_t sin(float64_t x)
617  {
618  return std::sin(x);
619  }
620 
622  COMPLEX128_STDMATH(sin)
623 
624 
628  static inline float64_t asin(float64_t x)
629  {
630  return std::asin(x);
631  }
632 
635 
636 
640  static inline float64_t sinh(float64_t x)
641  {
642  return std::sinh(x);
643  }
644 
646  COMPLEX128_STDMATH(sinh)
647 
648 
652  static inline float64_t cos(float64_t x)
653  {
654  return std::cos(x);
655  }
656 
658  COMPLEX128_STDMATH(cos)
659 
660 
664  static inline float64_t acos(float64_t x)
665  {
666  return std::acos(x);
667  }
668 
671 
672 
676  static inline float64_t cosh(float64_t x)
677  {
678  return std::cosh(x);
679  }
680 
682  COMPLEX128_STDMATH(cosh)
684 
693  static inline float64_t log10(float64_t v)
694  {
695  return std::log10(v);
696  }
697 
699  COMPLEX128_STDMATH(log10)
700 
701 
705  static inline float64_t log2(float64_t v)
706  {
707  return std::log2(v);
708  }
709 
714  static inline float64_t log(float64_t v)
715  {
716  return std::log(v);
717  }
718 
720  COMPLEX128_STDMATH(log)
721 
722  static inline index_t floor_log(index_t n)
723  {
724  index_t i;
725  for (i = 0; n != 0; i++)
726  n >>= 1;
727 
728  return i;
729  }
731 
738  static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
739  {
740  ASSERT(len>0 && xy)
741 
742  float64_t area = 0.0;
743 
744  if (!reversed)
745  {
746  for (int i=1; i<len; i++)
747  area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
748  }
749  else
750  {
751  for (int i=1; i<len; i++)
752  area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
753  }
754 
755  return area;
756  }
757 
763  static bool strtof(const char* str, float32_t* float_result);
764 
770  static bool strtod(const char* str, float64_t* double_result);
771 
777  static bool strtold(const char* str, floatmax_t* long_double_result);
778 
783  static inline int64_t factorial(int32_t n)
784  {
785  int64_t res=1;
786  for (int i=2; i<=n; i++)
787  res*=i ;
788  return res ;
789  }
790 
798  static void init_random(uint32_t initseed=0)
799  {
800  if (initseed==0)
801  seed = CRandom::generate_seed();
802  else
803  seed=initseed;
804 
805  sg_rand->set_seed(seed);
806  }
807 
811  static inline uint64_t random()
812  {
813  return sg_rand->random_64();
814  }
815 
819  static inline uint64_t random(uint64_t min_value, uint64_t max_value)
820  {
821  return sg_rand->random(min_value, max_value);
822  }
823 
829  static inline int64_t random(int64_t min_value, int64_t max_value)
830  {
831  return sg_rand->random(min_value, max_value);
832  }
833 
839  static inline uint32_t random(uint32_t min_value, uint32_t max_value)
840  {
841  return sg_rand->random(min_value, max_value);
842  }
843 
849  static inline int32_t random(int32_t min_value, int32_t max_value)
850  {
851  return sg_rand->random(min_value, max_value);
852  }
853 
859  static inline float32_t random(float32_t min_value, float32_t max_value)
860  {
861  return sg_rand->random(min_value, max_value);
862  }
863 
869  static inline float64_t random(float64_t min_value, float64_t max_value)
870  {
871  return sg_rand->random(min_value, max_value);
872  }
873 
879  static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
880  {
881  return sg_rand->random(min_value, max_value);
882  }
883 
888  {
889  // sets up variables & makes sure rand_s.range == (0,1)
890  float32_t ret;
891  float32_t rand_u;
892  float32_t rand_v;
893  float32_t rand_s;
894  do
895  {
896  rand_u = static_cast<float32_t>(CMath::random(-1.0, 1.0));
897  rand_v = static_cast<float32_t>(CMath::random(-1.0, 1.0));
898  rand_s = rand_u*rand_u + rand_v*rand_v;
899  } while ((rand_s == 0) || (rand_s >= 1));
900 
901  // the meat & potatos, and then the mean & standard deviation shifting...
902  ret = static_cast<float32_t>(rand_u*CMath::sqrt(-2.0*CMath::log(rand_s)/rand_s));
903  ret = std_dev*ret + mean;
904  return ret;
905  }
906 
911  {
912  return sg_rand->normal_distrib(mean, std_dev);
913  }
914 
917  static inline float32_t randn_float()
918  {
919  return static_cast<float32_t>(normal_random(0.0, 1.0));
920  }
921 
924  static inline float64_t randn_double()
925  {
926  return sg_rand->std_normal_distrib();
927  }
929 
937  static int32_t gcd(int32_t a, int32_t b)
938  {
939  REQUIRE((a>=0 && b>0) || (b>=0 && a>0), "gcd(%d,%d) is not defined.\n",
940  a, b);
941 
942  if (1 == a || 1 == b)
943  return 1;
944 
945  while (0 < a && 0 < b)
946  {
947  if (a > b)
948  a %= b;
949  else
950  b %= a;
951  }
952 
953  return 0 == a ? b : a;
954  }
955 
961  template <class T>
962  static void permute(SGVector<T> v, CRandom* rand=NULL)
963  {
964  if (rand)
965  {
966  for (index_t i=0; i<v.vlen; ++i)
967  swap(v[i], v[rand->random(i, v.vlen-1)]);
968  }
969  else
970  {
971  for (index_t i=0; i<v.vlen; ++i)
972  swap(v[i], v[random(i, v.vlen-1)]);
973  }
974  }
975 
981  template <class T>
982  static int32_t get_num_nonzero(T* vec, int32_t len)
983  {
984  int32_t nnz = 0;
985  for (index_t i=0; i<len; ++i)
986  nnz += vec[i] != 0;
987 
988  return nnz;
989  }
990 
996  static int32_t get_num_nonzero(complex128_t* vec, int32_t len)
997  {
998  int32_t nnz=0;
999  for (index_t i=0; i<len; ++i)
1000  nnz+=vec[i]!=0.0;
1001  return nnz;
1002  }
1003 
1009  static inline int64_t nchoosek(int32_t n, int32_t k)
1010  {
1011  int64_t res=1;
1012 
1013  for (int32_t i=n-k+1; i<=n; i++)
1014  res*=i;
1015 
1016  return res/factorial(k);
1017  }
1018 
1025  static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100);
1026 
1027 #ifndef SWIG // SWIG should skip this part
1028 
1034  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1035  static float64_t* linspace(T start, T end, int32_t n)
1036  {
1037  float64_t* output = SG_MALLOC(float64_t, n);
1038  linspace(output, start, end, n);
1039 
1040  return output;
1041  }
1042 #endif
1043 
1050  template <class T>
1051  static SGVector<float64_t> linspace_vec(T start, T end, int32_t n)
1052  {
1053  return SGVector<float64_t>(linspace(start, end, n), n);
1054  }
1055 
1061  template <class T>
1062  static T log_sum_exp(SGVector<T> values)
1063  {
1064  REQUIRE(values.vector, "Values are empty");
1065  REQUIRE(values.vlen>0,"number of values supplied is 0\n");
1066 
1067  /* find minimum element index */
1068  index_t min_index=0;
1069  T X0=values[0];
1070  if (values.vlen>1)
1071  {
1072  for (index_t i=1; i<values.vlen; ++i)
1073  {
1074  if (values[i]<X0)
1075  {
1076  X0=values[i];
1077  min_index=i;
1078  }
1079  }
1080  }
1081 
1082  /* remove element from vector copy and compute log sum exp */
1083  SGVector<T> values_without_X0(values.vlen-1);
1084  index_t from_idx=0;
1085  index_t to_idx=0;
1086  for (from_idx=0; from_idx<values.vlen; ++from_idx)
1087  {
1088  if (from_idx!=min_index)
1089  {
1090  values_without_X0[to_idx]=exp(values[from_idx]-X0);
1091  to_idx++;
1092  }
1093  }
1094 
1095  return X0+log(SGVector<T>::sum(values_without_X0)+1);
1096  }
1097 
1104  template <class T>
1105  static T log_mean_exp(SGVector<T> values)
1106  {
1107  return log_sum_exp(values) - log(values.vlen);
1108  }
1109 
1117  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
1118 
1125  static void sort(float64_t *a, int32_t*idx, int32_t N);
1126 
1127 #ifndef SWIG // SWIG should skip this part
1128 
1133  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1134  static void qsort(T* output, int32_t size)
1135  {
1136  if (size<=1)
1137  return;
1138 
1139  if (size==2)
1140  {
1141  if (output[0] > output [1])
1142  CMath::swap(output[0],output[1]);
1143  return;
1144  }
1145  //T split=output[random(0,size-1)];
1146  T split=output[size/2];
1147 
1148  int32_t left=0;
1149  int32_t right=size-1;
1150 
1151  while (left<=right)
1152  {
1153  while (output[left] < split)
1154  left++;
1155  while (output[right] > split)
1156  right--;
1157 
1158  if (left<=right)
1159  {
1160  CMath::swap(output[left],output[right]);
1161  left++;
1162  right--;
1163  }
1164  }
1165 
1166  if (right+1> 1)
1167  qsort(output,right+1);
1168 
1169  if (size-left> 1)
1170  qsort(&output[left],size-left);
1171  }
1172 
1178  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1179  static void insertion_sort(T* output, int32_t size)
1180  {
1181  for (int32_t i=0; i<size-1; i++)
1182  {
1183  int32_t j=i-1;
1184  T value=output[i];
1185  while (j >= 0 && output[j] > value)
1186  {
1187  output[j+1] = output[j];
1188  j--;
1189  }
1190  output[j+1]=value;
1191  }
1192  }
1193 
1198  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1199  inline static void radix_sort(T* array, int32_t size)
1200  {
1201  radix_sort_helper(array,size,0);
1202  }
1203 #endif
1204 
1212  template <class T>
1213  static inline uint8_t byte(T word, uint16_t p)
1214  {
1215  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
1216  }
1217 
1219  static inline uint8_t byte(complex128_t word, uint16_t p)
1220  {
1221  SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
1222  return uint8_t(0);
1223  }
1224 
1230  template <class T>
1231  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
1232  {
1233  static size_t count[256], nc, cmin;
1234  T *ak;
1235  uint8_t c=0;
1236  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
1237  T *an, *aj, *pile[256];
1238  size_t *cp, cmax;
1239 
1240  /* Push initial array to stack */
1241  sp = s;
1242  radix_push(array, size, i);
1243 
1244  /* Loop until all digits have been sorted */
1245  while (sp>s) {
1246  radix_pop(array, size, i);
1247  an = array + size;
1248 
1249  /* Make character histogram */
1250  if (nc == 0) {
1251  cmin = 0xff;
1252  for (ak = array; ak < an; ak++) {
1253  c = byte(*ak, i);
1254  count[c]++;
1255  if (count[c] == 1) {
1256  /* Determine smallest character */
1257  if (c < cmin)
1258  cmin = c;
1259  nc++;
1260  }
1261  }
1262 
1263  /* Sort recursively if stack size too small */
1264  if (sp + nc > s + RADIX_STACK_SIZE) {
1265  radix_sort_helper(array, size, i);
1266  continue;
1267  }
1268  }
1269 
1270  /*
1271  * Set pile[]; push incompletely sorted bins onto stack.
1272  * pile[] = pointers to last out-of-place element in bins.
1273  * Before permuting: pile[c-1] + count[c] = pile[c];
1274  * during deal: pile[c] counts down to pile[c-1].
1275  */
1276  olds = bigs = sp;
1277  cmax = 2;
1278  ak = array;
1279  pile[0xff] = an;
1280  for (cp = count + cmin; nc > 0; cp++) {
1281  /* Find next non-empty pile */
1282  while (*cp == 0)
1283  cp++;
1284  /* Pile with several entries */
1285  if (*cp > 1) {
1286  /* Determine biggest pile */
1287  if (*cp > cmax) {
1288  cmax = *cp;
1289  bigs = sp;
1290  }
1291 
1292  if (i < sizeof(T)-1)
1293  radix_push(ak, *cp, (uint16_t) (i + 1));
1294  }
1295  pile[cp - count] = ak += *cp;
1296  nc--;
1297  }
1298 
1299  /* Play it safe -- biggest bin last. */
1300  swap(*olds, *bigs);
1301 
1302  /*
1303  * Permute misplacements home. Already home: everything
1304  * before aj, and in pile[c], items from pile[c] on.
1305  * Inner loop:
1306  * r = next element to put in place;
1307  * ak = pile[r[i]] = location to put the next element.
1308  * aj = bottom of 1st disordered bin.
1309  * Outer loop:
1310  * Once the 1st disordered bin is done, ie. aj >= ak,
1311  * aj<-aj + count[c] connects the bins in array linked list;
1312  * reset count[c].
1313  */
1314  aj = array;
1315  while (aj<an)
1316  {
1317  T r;
1318 
1319  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
1320  swap(*ak, r);
1321 
1322  *aj = r;
1323  aj += count[c];
1324  count[c] = 0;
1325  }
1326  }
1327  }
1328 
1330  static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
1331  {
1332  SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
1333  }
1334 
1335 #ifndef SWIG // SWIG should skip this part
1336 
1342  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1343  static void qsort(T** vector, index_t length)
1344  {
1345  if (length<=1)
1346  return;
1347 
1348  if (length==2)
1349  {
1350  if (*vector[0]>*vector[1])
1351  swap(vector[0],vector[1]);
1352  return;
1353  }
1354  T* split=vector[length/2];
1355 
1356  int32_t left=0;
1357  int32_t right=length-1;
1358 
1359  while (left<=right)
1360  {
1361  while (*vector[left]<*split)
1362  ++left;
1363  while (*vector[right]>*split)
1364  --right;
1365 
1366  if (left<=right)
1367  {
1368  swap(vector[left],vector[right]);
1369  ++left;
1370  --right;
1371  }
1372  }
1373 
1374  if (right+1>1)
1375  qsort(vector,right+1);
1376 
1377  if (length-left>1)
1378  qsort(&vector[left],length-left);
1379  }
1380 
1384  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1385  static void qsort(SGVector<T> vector)
1386  {
1387  qsort<T>(vector, vector.size());
1388  }
1389 #endif
1390 
1392  template<class T>
1394  {
1396  IndexSorter(const SGVector<T> *vec) { data = vec->vector; }
1397 
1399  bool operator() (index_t i, index_t j) const
1400  {
1401  return abs(data[i]-data[j])>std::numeric_limits<T>::epsilon()
1402  && data[i]<data[j];
1403  }
1404 
1406  const T* data;
1407  };
1408 
1409 #ifndef SWIG // SWIG should skip this part
1410 
1417  template<class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1419  {
1420  IndexSorter<T> cmp(&vector);
1421  SGVector<index_t> idx(vector.size());
1422  for (index_t i=0; i < vector.size(); ++i)
1423  idx[i] = i;
1424 
1425  std::sort(idx.vector, idx.vector+vector.size(), cmp);
1426 
1427  return idx;
1428  }
1429 
1435  template <class T, class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1436  static bool is_sorted(SGVector<T> vector)
1437  {
1438  if (vector.size() < 2)
1439  return true;
1440 
1441  for(int32_t i=1; i<vector.size(); i++)
1442  {
1443  if (vector[i-1] > vector[i])
1444  return false;
1445  }
1446 
1447  return true;
1448  }
1449 #endif
1450 
1455  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1456  {
1457  ASSERT(width>=0)
1458  for (int i=0; i<width; i++)
1459  {
1460  T mask = ((T) 1)<<(sizeof(T)*8-1);
1461  while (mask)
1462  {
1463  if (mask & word)
1464  SG_SPRINT("1")
1465  else
1466  SG_SPRINT("0")
1467 
1468  mask>>=1;
1469  }
1470  }
1471  }
1472 
1474  static void display_bits(complex128_t word,
1475  int32_t width=8*sizeof(complex128_t))
1476  {
1477  SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
1478  }
1479 
1488  template <class T1,class T2>
1489  static void qsort_index(T1* output, T2* index, uint32_t size);
1490 
1492  template <class T>
1493  static void qsort_index(complex128_t* output, T* index, uint32_t size)
1494  {
1495  SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
1496  }
1497 
1506  template <class T1,class T2>
1507  static void qsort_backward_index(
1508  T1* output, T2* index, int32_t size);
1509 
1511  template <class T>
1513  complex128_t* output, T* index, uint32_t size)
1514  {
1515  SG_SERROR("CMath::qsort_backword_index():: \
1516  Not supported for complex128_t\n");
1517  }
1518 
1530  template <class T1,class T2>
1531  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1532  {
1533  int32_t n=0;
1534  thread_qsort<T1,T2> t;
1535  t.output=output;
1536  t.index=index;
1537  t.size=size;
1538  t.qsort_threads=&n;
1539  t.sort_limit=limit;
1540  t.num_threads=n_threads;
1541  parallel_qsort_index<T1,T2>(&t);
1542  }
1543 
1545  template <class T>
1546  inline static void parallel_qsort_index(complex128_t* output, T* index,
1547  uint32_t size, int32_t n_threads, int32_t limit=0)
1548  {
1549  SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1550  }
1551 
1553  template <class T1,class T2>
1554  static void* parallel_qsort_index(void* p);
1555 
1556 
1563  template <class T>
1564  static void min(float64_t* output, T* index, int32_t size);
1565 
1567  static void min(float64_t* output, complex128_t* index, int32_t size)
1568  {
1569  SG_SERROR("CMath::min():: Not supported for complex128_t\n");
1570  }
1571 
1575  template <class T>
1576  static void nmin(
1577  float64_t* output, T* index, int32_t size, int32_t n);
1578 
1580  static void nmin(float64_t* output, complex128_t* index,
1581  int32_t size, int32_t n)
1582  {
1583  SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
1584  }
1585 
1586 
1587 
1591  template <class T>
1592  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1593  {
1594  int32_t start=0;
1595  int32_t end=size-1;
1596 
1597  if (size<1)
1598  return 0;
1599 
1600  while (start<end)
1601  {
1602  int32_t middle=(start+end)/2;
1603 
1604  if (output[middle]>elem)
1605  end=middle-1;
1606  else if (output[middle]<elem)
1607  start=middle+1;
1608  else
1609  return middle;
1610  }
1611 
1612  return start;
1613  }
1614 
1616  static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
1617  {
1618  SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
1619  return int32_t(0);
1620  }
1621 
1628  template <class T>
1629  static inline int32_t binary_search(T* output, int32_t size, T elem)
1630  {
1631  int32_t ind = binary_search_helper(output, size, elem);
1632  if (ind >= 0 && output[ind] == elem)
1633  return ind;
1634  return -1;
1635  }
1636 
1638  static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
1639  {
1640  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1641  return int32_t(-1);
1642  }
1643 
1644 
1652  template<class T>
1653  static inline int32_t binary_search(T** vector, index_t length,
1654  T* elem)
1655  {
1656  int32_t start=0;
1657  int32_t end=length-1;
1658 
1659  if (length<1)
1660  return -1;
1661 
1662  while (start<end)
1663  {
1664  int32_t middle=(start+end)/2;
1665 
1666  if (*vector[middle]>*elem)
1667  end=middle-1;
1668  else if (*vector[middle]<*elem)
1669  start=middle+1;
1670  else
1671  {
1672  start=middle;
1673  break;
1674  }
1675  }
1676 
1677  if (start>=0&&*vector[start]==*elem)
1678  return start;
1679 
1680  return -1;
1681  }
1682 
1684  static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
1685  {
1686  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1687  return int32_t(-1);
1688  }
1689 
1696  template <class T>
1698  T* output, int32_t size, T elem)
1699  {
1700  int32_t ind = binary_search_helper(output, size, elem);
1701 
1702  if (output[ind]<=elem)
1703  return ind;
1704  if (ind>0 && output[ind-1] <= elem)
1705  return ind-1;
1706  return -1;
1707  }
1708 
1711  int32_t size, complex128_t elem)
1712  {
1713  SG_SERROR("CMath::binary_search_max_lower_equal():: \
1714  Not supported for complex128_t\n");
1715  return int32_t(-1);
1716  }
1717 
1726  static float64_t Align(
1727  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1728 
1729 
1731 
1737  {
1738  return c.real();
1739  }
1740 
1746  {
1747  return c.imag();
1748  }
1749 
1751  inline static uint32_t get_seed()
1752  {
1753  return CMath::seed;
1754  }
1755 
1757  inline static uint32_t get_log_range()
1758  {
1759  return CMath::LOGRANGE;
1760  }
1761 
1762 #ifdef USE_LOGCACHE
1763  inline static uint32_t get_log_accuracy()
1765  {
1766  return CMath::LOGACCURACY;
1767  }
1768 #endif
1769 
1771  static int is_finite(double f);
1772 
1774  static int is_infinity(double f);
1775 
1777  static int is_nan(double f);
1778 
1788 #ifdef USE_LOGCACHE
1789  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1790  {
1791  float64_t diff;
1792 
1793  if (!CMath::is_finite(p))
1794  return q;
1795 
1796  if (!CMath::is_finite(q))
1797  {
1798  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
1799  return NOT_A_NUMBER;
1800  }
1801  diff = p - q;
1802  if (diff > 0)
1803  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1804  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1805  }
1806 
1808  static void init_log_table();
1809 
1811  static int32_t determine_logrange();
1812 
1814  static int32_t determine_logaccuracy(int32_t range);
1815 #else
1817  float64_t p, float64_t q)
1818  {
1819  float64_t diff;
1820 
1821  if (!CMath::is_finite(p))
1822  return q;
1823  if (!CMath::is_finite(q))
1824  return p;
1825  diff = p - q;
1826  if (diff > 0)
1827  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
1828  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
1829  }
1830 #endif
1831 #ifdef USE_LOGSUMARRAY
1832 
1837  static inline float64_t logarithmic_sum_array(
1838  float64_t *p, int32_t len)
1839  {
1840  if (len<=2)
1841  {
1842  if (len==2)
1843  return logarithmic_sum(p[0],p[1]) ;
1844  if (len==1)
1845  return p[0];
1846  return -INFTY ;
1847  }
1848  else
1849  {
1850  float64_t *pp=p ;
1851  if (len%2==1) pp++ ;
1852  for (int32_t j=0; j < len>>1; j++)
1853  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
1854  }
1855  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1856  }
1857 #endif //USE_LOGSUMARRAY
1858 
1859 
1861  virtual const char* get_name() const { return "Math"; }
1862  public:
1865  static const float64_t NOT_A_NUMBER;
1868  static const float64_t INFTY;
1869  static const float64_t ALMOST_INFTY;
1870 
1873 
1875  static const float64_t PI;
1876 
1879 
1880  /* Largest and smallest possible float64_t */
1883 
1884  /* Floating point Limits, Normalized */
1885  static const float32_t F_MAX_VAL32;
1887  static const float64_t F_MAX_VAL64;
1889 
1890  /* Floating point limits, Denormalized */
1891  static const float32_t F_MIN_VAL32;
1892  static const float64_t F_MIN_VAL64;
1893 
1894  protected:
1896  static int32_t LOGRANGE;
1897 
1899  static uint32_t seed;
1900 
1901 #ifdef USE_LOGCACHE
1902 
1904  static int32_t LOGACCURACY;
1906  static float64_t* logtable;
1908 #endif
1909 };
1910 
1911 //implementations of template functions
1912 template <class T1,class T2>
1913 void* CMath::parallel_qsort_index(void* p)
1914  {
1915  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1916  T1* output=ps->output;
1917  T2* index=ps->index;
1918  int32_t size=ps->size;
1919  int32_t* qsort_threads=ps->qsort_threads;
1920  int32_t sort_limit=ps->sort_limit;
1921  int32_t num_threads=ps->num_threads;
1922 
1923  if (size<2)
1924  {
1925  return NULL;
1926  }
1927 
1928  if (size==2)
1929  {
1930  if (output[0] > output [1])
1931  {
1932  swap(output[0], output[1]);
1933  swap(index[0], index[1]);
1934  }
1935  return NULL;
1936  }
1937  //T1 split=output[random(0,size-1)];
1938  T1 split=output[size/2];
1939 
1940  int32_t left=0;
1941  int32_t right=size-1;
1942 
1943  while (left<=right)
1944  {
1945  while (output[left] < split)
1946  left++;
1947  while (output[right] > split)
1948  right--;
1949 
1950  if (left<=right)
1951  {
1952  swap(output[left], output[right]);
1953  swap(index[left], index[right]);
1954  left++;
1955  right--;
1956  }
1957  }
1958  bool lthread_start=false;
1959  bool rthread_start=false;
1960  pthread_t lthread;
1961  pthread_t rthread;
1962  struct thread_qsort<T1,T2> t1;
1963  struct thread_qsort<T1,T2> t2;
1964 
1965  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1966  qsort_index(output,index,right+1);
1967  else if (right+1> 1)
1968  {
1969  (*qsort_threads)++;
1970  lthread_start=true;
1971  t1.output=output;
1972  t1.index=index;
1973  t1.size=right+1;
1974  t1.qsort_threads=qsort_threads;
1975  t1.sort_limit=sort_limit;
1976  t1.num_threads=num_threads;
1977  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1978  {
1979  lthread_start=false;
1980  (*qsort_threads)--;
1981  qsort_index(output,index,right+1);
1982  }
1983  }
1984 
1985 
1986  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1987  qsort_index(&output[left],&index[left], size-left);
1988  else if (size-left> 1)
1989  {
1990  (*qsort_threads)++;
1991  rthread_start=true;
1992  t2.output=&output[left];
1993  t2.index=&index[left];
1994  t2.size=size-left;
1995  t2.qsort_threads=qsort_threads;
1996  t2.sort_limit=sort_limit;
1997  t2.num_threads=num_threads;
1998  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1999  {
2000  rthread_start=false;
2001  (*qsort_threads)--;
2002  qsort_index(&output[left],&index[left], size-left);
2003  }
2004  }
2005 
2006  if (lthread_start)
2007  {
2008  pthread_join(lthread, NULL);
2009  (*qsort_threads)--;
2010  }
2011 
2012  if (rthread_start)
2013  {
2014  pthread_join(rthread, NULL);
2015  (*qsort_threads)--;
2016  }
2017 
2018  return NULL;
2019  }
2020 
2021  template <class T1,class T2>
2022 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
2023 {
2024  if (size<=1)
2025  return;
2026 
2027  if (size==2)
2028  {
2029  if (output[0] > output [1])
2030  {
2031  swap(output[0],output[1]);
2032  swap(index[0],index[1]);
2033  }
2034  return;
2035  }
2036  //T1 split=output[random(0,size-1)];
2037  T1 split=output[size/2];
2038 
2039  int32_t left=0;
2040  int32_t right=size-1;
2041 
2042  while (left<=right)
2043  {
2044  while (output[left] < split)
2045  left++;
2046  while (output[right] > split)
2047  right--;
2048 
2049  if (left<=right)
2050  {
2051  swap(output[left],output[right]);
2052  swap(index[left],index[right]);
2053  left++;
2054  right--;
2055  }
2056  }
2057 
2058  if (right+1> 1)
2059  qsort_index(output,index,right+1);
2060 
2061  if (size-left> 1)
2062  qsort_index(&output[left],&index[left], size-left);
2063 }
2064 
2065  template <class T1,class T2>
2066 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
2067 {
2068  if (size<=1)
2069  return;
2070 
2071  if (size==2)
2072  {
2073  if (output[0] < output [1])
2074  {
2075  swap(output[0],output[1]);
2076  swap(index[0],index[1]);
2077  }
2078  return;
2079  }
2080 
2081  //T1 split=output[random(0,size-1)];
2082  T1 split=output[size/2];
2083 
2084  int32_t left=0;
2085  int32_t right=size-1;
2086 
2087  while (left<=right)
2088  {
2089  while (output[left] > split)
2090  left++;
2091  while (output[right] < split)
2092  right--;
2093 
2094  if (left<=right)
2095  {
2096  swap(output[left],output[right]);
2097  swap(index[left],index[right]);
2098  left++;
2099  right--;
2100  }
2101  }
2102 
2103  if (right+1> 1)
2104  qsort_backward_index(output,index,right+1);
2105 
2106  if (size-left> 1)
2107  qsort_backward_index(&output[left],&index[left], size-left);
2108 }
2109 
2110  template <class T>
2111 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
2112 {
2113  if (6*n*size<13*size*CMath::log(size))
2114  for (int32_t i=0; i<n; i++)
2115  min(&output[i], &index[i], size-i);
2116  else
2117  qsort_index(output, index, size);
2118 }
2119 
2120 /* move the smallest entry in the array to the beginning */
2121  template <class T>
2122 void CMath::min(float64_t* output, T* index, int32_t size)
2123 {
2124  if (size<=1)
2125  return;
2126  float64_t min_elem=output[0];
2127  int32_t min_index=0;
2128  for (int32_t i=1; i<size; i++)
2129  {
2130  if (output[i]<min_elem)
2131  {
2132  min_index=i;
2133  min_elem=output[i];
2134  }
2135  }
2136  swap(output[0], output[min_index]);
2137  swap(index[0], index[min_index]);
2138 }
2139 
2140 #define COMPLEX128_ERROR_ONEARG_T(function) \
2141 template <> \
2142 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
2143 { \
2144  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
2145  #function);\
2146  return complex128_t(0.0, 0.0); \
2147 }
2148 
2150 // COMPLEX128_ERROR_ONEARG_T(sign)
2151 
2152 }
2153 #undef COMPLEX128_ERROR_ONEARG
2154 #undef COMPLEX128_ERROR_ONEARG_T
2155 #undef COMPLEX128_STDMATH
2156 #endif
static float64_t normal_random(float64_t mean, float64_t std_dev)
Definition: Math.h:910
static const float32_t F_MAX_VAL32
Definition: Math.h:1885
static const float64_t MACHINE_EPSILON
Definition: Math.h:1878
static int32_t gcd(int32_t a, int32_t b)
Definition: Math.h:937
static void permute(SGVector< T > v, CRandom *rand=NULL)
Definition: Math.h:962
static int32_t binary_search(complex128_t *output, int32_t size, complex128_t elem)
binary_search not implemented for complex128_t
Definition: Math.h:1638
static void parallel_qsort_index(T1 *output, T2 *index, uint32_t size, int32_t n_threads, int32_t limit=262144)
Definition: Math.h:1531
static uint32_t seed
random generator seed
Definition: Math.h:1899
static void qsort(SGVector< T > vector)
Definition: Math.h:1385
std::complex< float64_t > complex128_t
Definition: common.h:77
uint64_t random(uint64_t min_value, uint64_t max_value)
Definition: Random.h:101
float64_t std_normal_distrib() const
Definition: Random.cpp:268
static T max(T a, T b)
Definition: Math.h:149
static floatmax_t random(floatmax_t min_value, floatmax_t max_value)
Definition: Math.h:879
int32_t index_t
Definition: common.h:72
static float64_t ceil(float64_t d)
Definition: Math.h:384
static complex128_t pow(complex128_t x, int32_t n)
Definition: Math.h:515
static const float64_t INFTY
infinity
Definition: Math.h:1868
static SGVector< float64_t > linspace_vec(T start, T end, int32_t n)
Definition: Math.h:1051
static int32_t binary_search_helper(T *output, int32_t size, T elem)
Definition: Math.h:1592
static void qsort_index(complex128_t *output, T *index, uint32_t size)
qsort_index not implemented for complex128_t
Definition: Math.h:1493
static T sqrt(T x)
Definition: Math.h:428
#define SG_SWARNING(...)
Definition: SGIO.h:163
static float32_t normal_random(float32_t mean, float32_t std_dev)
Definition: Math.h:887
uint64_t random_64() const
Definition: Random.cpp:137
static float64_t random(float64_t min_value, float64_t max_value)
Definition: Math.h:869
static int32_t binary_search(T *output, int32_t size, T elem)
Definition: Math.h:1629
static T sq(T x)
Definition: Math.h:418
static uint32_t get_seed()
returns number generator seed
Definition: Math.h:1751
static float32_t randn_float()
Definition: Math.h:917
static const float64_t MIN_REAL_NUMBER
Definition: Math.h:1882
static float64_t randn_double()
Definition: Math.h:924
static int32_t binary_search_max_lower_equal(T *output, int32_t size, T elem)
Definition: Math.h:1697
#define REQUIRE(x,...)
Definition: SGIO.h:181
static const float64_t F_MAX_VAL64
Definition: Math.h:1887
static T clamp(T value, T lb, T ub)
Definition: Math.h:217
static uint32_t get_log_range()
returns range of logtable
Definition: Math.h:1757
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:859
static const float32_t F_MIN_VAL32
Definition: Math.h:1891
CRandom * sg_rand
Definition: init.cpp:47
Range< T > range(T rend)
Definition: range.h:136
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:1396
#define RADIX_STACK_SIZE
Definition: Math.h:57
static void qsort(T **vector, index_t length)
Definition: Math.h:1343
static float64_t imag(complex128_t c)
Definition: Math.h:1745
static float64_t floor(float64_t d)
Definition: Math.h:375
static void insertion_sort(T *output, int32_t size)
Definition: Math.h:1179
static SGVector< index_t > argsort(SGVector< T > vector)
Definition: Math.h:1418
static uint64_t random()
Definition: Math.h:811
static int32_t arg_min(T *vec, int32_t inc, int32_t len, T *minv_ptr=NULL)
Definition: Math.h:262
static int32_t get_num_nonzero(complex128_t *vec, int32_t len)
Definition: Math.h:996
static const float32_t F_MIN_NORM_VAL32
Definition: Math.h:1886
static float64_t real(complex128_t c)
Definition: Math.h:1736
static uint8_t byte(complex128_t word, uint16_t p)
byte not implemented for complex128_t
Definition: Math.h:1219
static int32_t LOGRANGE
range for logtable: log(1+exp(x)) -LOGRANGE <= x <= 0
Definition: Math.h:1896
static const float64_t ALMOST_NEG_INFTY
almost neg (log) infinity
Definition: Math.h:1872
static complex128_t pow(complex128_t x, complex128_t n)
Definition: Math.h:524
static bool fequals(const T &a, const T &b, const float64_t eps, bool tolerant=false)
Definition: Math.h:308
static uint8_t byte(T word, uint16_t p)
Definition: Math.h:1213
float64_t normal_distrib(float64_t mu, float64_t sigma) const
Definition: Random.cpp:263
static complex128_t pow(float64_t x, complex128_t n)
Definition: Math.h:542
#define SG_SPRINT(...)
Definition: SGIO.h:165
static float64_t pow(float64_t x, int32_t n)
Definition: Math.h:497
#define ASSERT(x)
Definition: SGIO.h:176
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:124
static int32_t random(int32_t min_value, int32_t max_value)
Definition: Math.h:849
static float64_t pow(float64_t x, float64_t n)
Definition: Math.h:506
static bool is_sorted(SGVector< T > vector)
Definition: Math.h:1436
int32_t size() const
Definition: SGVector.h:156
static void init_random(uint32_t initseed=0)
Definition: Math.h:798
static int32_t pow(int32_t x, int32_t n)
Definition: Math.h:483
#define radix_pop(a, n, i)
Definition: Math.h:61
double float64_t
Definition: common.h:60
static float64_t * linspace(T start, T end, int32_t n)
Definition: Math.h:1035
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:1546
static void min(float64_t *output, complex128_t *index, int32_t size)
complex128_t cannot be used as index
Definition: Math.h:1567
long double floatmax_t
Definition: common.h:61
#define COMPLEX128_ERROR_ONEARG(function)
Definition: Math.h:94
static void radix_sort(T *array, int32_t size)
Definition: Math.h:1199
shogun vector
static void qsort_backword_index(complex128_t *output, T *index, uint32_t size)
qsort_backword_index not implemented for complex128_t
Definition: Math.h:1512
static int32_t get_num_nonzero(T *vec, int32_t len)
Definition: Math.h:982
static T log_sum_exp(SGVector< T > values)
Definition: Math.h:1062
static float64_t area_under_curve(float64_t *xy, int32_t len, bool reversed)
Definition: Math.h:738
static void display_bits(T word, int32_t width=8 *sizeof(T))
Definition: Math.h:1455
static uint64_t random(uint64_t min_value, uint64_t max_value)
Definition: Math.h:819
static bool fequals_abs(const T &a, const T &b, const float64_t eps)
Definition: Math.h:291
static int64_t factorial(int32_t n)
Definition: Math.h:783
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:1710
float float32_t
Definition: common.h:59
: Pseudo random number geneartor
Definition: Random.h:34
static int64_t nchoosek(int32_t n, int32_t k)
Definition: Math.h:1009
static T min(T *vec, int32_t len)
Definition: Math.h:191
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
static int32_t binary_search(T **vector, index_t length, T *elem)
Definition: Math.h:1653
void set_seed(uint32_t seed)
Definition: Random.cpp:56
static float64_t abs(complex128_t a)
Definition: Math.h:177
static void display_bits(complex128_t word, int32_t width=8 *sizeof(complex128_t))
disply_bits not implemented for complex128_t
Definition: Math.h:1474
#define COMPLEX128_STDMATH(function)
Definition: Math.h:102
static void radix_sort_helper(T *array, int32_t size, uint16_t i)
Definition: Math.h:1231
static T sign(T a)
Definition: Math.h:394
virtual const char * get_name() const
Definition: Math.h:1861
#define SG_SERROR(...)
Definition: SGIO.h:164
static float64_t exp(float64_t x)
Definition: Math.h:551
static const float64_t F_MIN_VAL64
Definition: Math.h:1892
static const float64_t F_MIN_NORM_VAL64
Definition: Math.h:1888
static float64_t log(float64_t v)
Definition: Math.h:714
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:1330
static const float64_t ALMOST_INFTY
Definition: Math.h:1869
Class which collects generic mathematical functions.
Definition: Math.h:114
static T log_mean_exp(SGVector< T > values)
Definition: Math.h:1105
static void swap(T &a, T &b)
Definition: Math.h:406
static complex128_t pow(complex128_t x, float64_t n)
Definition: Math.h:533
static T max(T *vec, int32_t len)
Definition: Math.h:203
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:1616
static float64_t round(float64_t d)
Definition: Math.h:366
T max(const Container< T > &a)
static uint32_t random(uint32_t min_value, uint32_t max_value)
Definition: Math.h:839
static floatmax_t powl(floatmax_t x, floatmax_t n)
Definition: Math.h:469
static float64_t logarithmic_sum(float64_t p, float64_t q)
Definition: Math.h:1816
#define radix_push(a, n, i)
Definition: Math.h:60
static T min(T a, T b)
Definition: Math.h:138
static int32_t binary_search(complex128_t **vector, index_t length, complex128_t *elem)
binary_search not implemented for complex128_t
Definition: Math.h:1684
static int32_t pow(bool x, int32_t n)
Definition: Math.h:474
static const float64_t MAX_REAL_NUMBER
Definition: Math.h:1881
index_t vlen
Definition: SGVector.h:571
static int32_t arg_max(T *vec, int32_t inc, int32_t len, T *maxv_ptr=NULL)
Definition: Math.h:235
static T abs(T a)
Definition: Math.h:161
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:1580
static int64_t random(int64_t min_value, int64_t max_value)
Definition: Math.h:829
static void qsort(T *output, int32_t size)
Definition: Math.h:1134
static const float64_t PI
Definition: Math.h:1875

SHOGUN Machine Learning Toolbox - Documentation