SHOGUN  3.2.1
 全部  命名空间 文件 函数 变量 类型定义 枚举 枚举值 友元 宏定义 
Math.h
浏览该文件的文档.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Thoralf Klein
8  * Written (W) 2013 Soumyajit De
9  * Written (W) 2012 Fernando Jose Iglesias Garcia
10  * Written (W) 2011 Siddharth Kherada
11  * Written (W) 2011 Justin Patera
12  * Written (W) 2011 Alesis Novik
13  * Written (W) 2011-2013 Heiko Strathmann
14  * Written (W) 1999-2009 Soeren Sonnenburg
15  * Written (W) 1999-2008 Gunnar Raetsch
16  * Written (W) 2007 Konrad Rieck
17  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
18  */
19 
20 #ifndef __MATHEMATICS_H_
21 #define __MATHEMATICS_H_
22 
23 #include <shogun/lib/config.h>
24 
25 #include <shogun/base/SGObject.h>
26 #include <shogun/lib/common.h>
27 #include <shogun/io/SGIO.h>
28 #include <shogun/base/Parallel.h>
30 #include <shogun/lib/SGVector.h>
31 
32 #ifndef _USE_MATH_DEFINES
33 #define _USE_MATH_DEFINES
34 #endif
35 
36 #include <math.h>
37 #include <stdio.h>
38 #include <float.h>
39 #include <sys/types.h>
40 #ifndef _WIN32
41 #include <unistd.h>
42 #endif
43 
44 #ifdef HAVE_PTHREAD
45 #include <pthread.h>
46 #endif
47 
48 #ifdef SUNOS
49 #include <ieeefp.h>
50 #endif
51 
53 #ifdef log2
54 #define cygwin_log2 log2
55 #undef log2
56 #endif
57 
58 #ifndef M_PI
59 #define M_PI 3.14159265358979323846
60 #endif
61 
62 #ifdef _WIN32
63 #ifndef isnan
64 #define isnan _isnan
65 #endif
66 
67 #ifndef isfinite
68 #define isfinite _isfinite
69 #endif
70 #endif //_WIN32
71 
72 /* Size of RNG seed */
73 #define RNG_SEED_SIZE 256
74 
75 /* Maximum stack size */
76 #define RADIX_STACK_SIZE 512
77 
78 /* Stack macros */
79 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
80 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
81 
82 #ifndef DOXYGEN_SHOULD_SKIP_THIS
83 
84 template <class T> struct radix_stack_t
85 {
87  T *sa;
89  size_t sn;
91  uint16_t si;
92 };
93 
95 
97 template <class T1, class T2> struct thread_qsort
98 {
100  T1* output;
102  T2* index;
104  uint32_t size;
105 
107  int32_t* qsort_threads;
109  int32_t sort_limit;
111  int32_t num_threads;
112 };
113 #endif // DOXYGEN_SHOULD_SKIP_THIS
114 
115 #define COMPLEX128_ERROR_ONEARG(function) \
116 static inline complex128_t function(complex128_t a) \
117 { \
118  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
119  #function);\
120  return complex128_t(0.0, 0.0); \
121 }
122 
123 #define COMPLEX128_STDMATH(function) \
124 static inline complex128_t function(complex128_t a) \
125 { \
126  return std::function(a); \
127 }
128 
129 namespace shogun
130 {
132  extern CRandom* sg_rand;
133  class CSGObject;
136 class CMath : public CSGObject
137 {
138  public:
142  CMath();
144 
146  virtual ~CMath();
148 
152 
154  //
155  template <class T>
156  static inline T min(T a, T b)
157  {
158  return (a<=b) ? a : b;
159  }
160 
162  template <class T>
163  static inline T max(T a, T b)
164  {
165  return (a>=b) ? a : b;
166  }
167 
169 
170  template <class T>
171  static T min(T* vec, int32_t len)
172  {
173  ASSERT(len>0)
174  T minv=vec[0];
175 
176  for (int32_t i=1; i<len; i++)
177  minv=min(vec[i], minv);
178 
179  return minv;
180  }
181 
183 
184  template <class T>
185  static T max(T* vec, int32_t len)
186  {
187  ASSERT(len>0)
188  T maxv=vec[0];
189 
190  for (int32_t i=1; i<len; i++)
191  maxv=max(vec[i], maxv);
192 
193  return maxv;
194  }
195 
197  template <class T>
198  static inline T clamp(T value, T lb, T ub)
199  {
200  if (value<=lb)
201  return lb;
202  else if (value>=ub)
203  return ub;
204  else
205  return value;
206  }
207 
209  template <class T>
210  static inline T abs(T a)
211  {
212  // can't be a>=0?(a):(-a), because compiler complains about
213  // 'comparison always true' when T is unsigned
214  if (a==0)
215  return 0;
216  else if (a>0)
217  return a;
218  else
219  return -a;
220  }
221 
223  static inline float64_t abs(complex128_t a)
224  {
225  float64_t a_real=a.real();
226  float64_t a_imag=a.imag();
227  return (CMath::sqrt(a_real*a_real+a_imag*a_imag));
228  }
230 
232  template <class T>
233  static int32_t arg_max(T * vec, int32_t inc, int32_t len, T * maxv_ptr = NULL)
234  {
235  ASSERT(len > 0 || inc > 0)
236 
237  T maxv = vec[0];
238  int32_t maxIdx = 0;
239 
240  for (int32_t i = 1, j = inc ; i < len ; i++, j += inc)
241  {
242  if (vec[j] > maxv)
243  maxv = vec[j], maxIdx = i;
244  }
245 
246  if (maxv_ptr != NULL)
247  *maxv_ptr = maxv;
248 
249  return maxIdx;
250  }
251 
253 
254  template <class T>
255  static int32_t arg_min(T * vec, int32_t inc, int32_t len, T * minv_ptr = NULL)
256  {
257  ASSERT(len > 0 || inc > 0)
258 
259  T minv = vec[0];
260  int32_t minIdx = 0;
261 
262  for (int32_t i = 1, j = inc ; i < len ; i++, j += inc)
263  {
264  if (vec[j] < minv)
265  minv = vec[j], minIdx = i;
266  }
267 
268  if (minv_ptr != NULL)
269  *minv_ptr = minv;
270 
271  return minIdx;
272  }
273 
276 
283  template <class T>
284  static inline bool fequals_abs(const T& a, const T& b,
285  const float64_t eps)
286  {
287  const T diff = CMath::abs<T>((a-b));
288  return (diff < eps);
289  }
290 
300  template <class T>
301  static inline bool fequals(const T& a, const T& b,
302  const float64_t eps, bool tolerant=false)
303  {
304  const T absA = CMath::abs<T>(a);
305  const T absB = CMath::abs<T>(b);
306  const T diff = CMath::abs<T>((a-b));
307  T comp;
308 
309  // Handle this separately since NAN is unordered
311  return true;
312 
313  // Required for JSON Serialization Tests
314  if (tolerant)
315  return CMath::fequals_abs<T>(a, b, eps);
316 
317  // handles float32_t and float64_t separately
318  if (sizeof(T) == 4)
320 
321  else
323 
324  if (a==b)
325  return true;
326 
327  // both a and b are 0 and relative error is less meaningful
328  else if ( (a==0) || (b==0) || (diff < comp) )
329  return (diff<(eps * comp));
330 
331  // use max(relative error, diff) to handle large eps
332  else
333  {
334  T check = ((diff/(absA + absB)) > diff)?
335  (diff/(absA + absB)):diff;
336  return (check < eps);
337  }
338  }
339 
340  /* Get the corresponding absolute tolerance for unit test given a relative tolerance
341  *
342  * Note that a unit test will be passed only when
343  * \f[
344  * |v_\text{true} - v_\text{predict}| \leq tol_\text{relative} * |v_\text{true}|
345  * \f] which is equivalent to
346  * \f[
347  * |v_\text{true} - v_\text{predict}| \leq tol_\text{absolute}
348  * \f] where
349  * \f[
350  * tol_\text{absolute} = tol_\text{relative} * |v_\text{true}|
351  * \f]
352  *
353  * @param true_value true value should be finite (neither NAN nor INF)
354  * @param rel_tolorance relative tolerance should be positive and less than 1.0
355  *
356  * @return the corresponding absolute tolerance
357  */
358  static float64_t get_abs_tolorance(float64_t true_value, float64_t rel_tolorance);
359 
360  static inline float64_t round(float64_t d)
361  {
362  return ::floor(d+0.5);
363  }
364 
365  static inline float64_t floor(float64_t d)
366  {
367  return ::floor(d);
368  }
369 
370  static inline float64_t ceil(float64_t d)
371  {
372  return ::ceil(d);
373  }
374 
376  template <class T>
377  static inline T sign(T a)
378  {
379  if (a==0)
380  return 0;
381  else return (a<0) ? (-1) : (+1);
382  }
383 
385  template <class T>
386  static inline void swap(T &a,T &b)
387  {
388  T c=a;
389  a=b;
390  b=c;
391  }
392 
394  template <class T>
395  static inline T sq(T x)
396  {
397  return x*x;
398  }
399 
401  static inline float32_t sqrt(float32_t x)
402  {
403  return ::sqrtf(x);
404  }
405 
407  static inline float64_t sqrt(float64_t x)
408  {
409  return ::sqrt(x);
410  }
411 
413  static inline floatmax_t sqrt(floatmax_t x)
414  {
415  //fall back to double precision sqrt if sqrtl is not
416  //available
417 #ifdef HAVE_SQRTL
418  return ::sqrtl(x);
419 #else
420  return ::sqrt(x);
421 #endif
422  }
423 
426 
427 
428  static inline float32_t invsqrt(float32_t x)
429  {
430  union float_to_int
431  {
432  float32_t f;
433  int32_t i;
434  };
435 
436  float_to_int tmp;
437  tmp.f=x;
438 
439  float32_t xhalf = 0.5f * x;
440  // store floating-point bits in integer tmp.i
441  // and do initial guess for Newton's method
442  tmp.i = 0x5f3759d5 - (tmp.i >> 1);
443  x = tmp.f; // convert new bits into float
444  x = x*(1.5f - xhalf*x*x); // One round of Newton's method
445  return x;
446  }
447 
449  static inline floatmax_t powl(floatmax_t x, floatmax_t n)
450  {
451  //fall back to double precision pow if powl is not
452  //available
453 #ifdef HAVE_POWL
454  return ::powl((long double) x, (long double) n);
455 #else
456  return ::pow((double) x, (double) n);
457 #endif
458  }
459 
460  static inline int32_t pow(bool x, int32_t n)
461  {
462  return 0;
463  }
464 
465  static inline int32_t pow(int32_t x, int32_t n)
466  {
467  ASSERT(n>=0)
468  int32_t result=1;
469  while (n--)
470  result*=x;
471 
472  return result;
473  }
474 
475  static inline float64_t pow(float64_t x, int32_t n)
476  {
477  if (n>=0)
478  {
479  float64_t result=1;
480  while (n--)
481  result*=x;
482 
483  return result;
484  }
485  else
486  return ::pow((double)x, (double)n);
487  }
488 
489  static inline float64_t pow(float64_t x, float64_t n)
490  {
491  return ::pow((double) x, (double) n);
492  }
493 
495  static inline complex128_t pow(complex128_t x, int32_t n)
496  {
497  return std::pow(x, n);
498  }
499 
501  {
502  return std::pow(x, n);
503  }
504 
506  {
507  return std::pow(x, n);
508  }
509 
511  {
512  return std::pow(x, n);
513  }
514 
515  static inline float64_t exp(float64_t x)
516  {
517  return ::exp((double) x);
518  }
519 
522 
523 
524  static inline float64_t tan(float64_t x)
525  {
526  return ::tan((double) x);
527  }
528 
531 
532 
533  static inline float64_t atan(float64_t x)
534  {
535  return ::atan((double) x);
536  }
537 
540 
541 
542  static inline float64_t atan2(float64_t x, float64_t y)
543  {
544  return ::atan2((double) x, (double) y);
545  }
546 
549 
550 
551  static inline float64_t tanh(float64_t x)
552  {
553  return ::tanh((double) x);
554  }
555 
558 
559  static inline float64_t log10(float64_t v)
560  {
561  return ::log(v)/::log(10.0);
562  }
563 
566 
567  static inline float64_t log2(float64_t v)
568  {
569 #ifdef HAVE_LOG2
570  return ::log2(v);
571 #else
572  return ::log(v)/::log(2.0);
573 #endif //HAVE_LOG2
574  }
575 
576  static inline float64_t log(float64_t v)
577  {
578  return ::log(v);
579  }
580 
583 
584  static inline index_t floor_log(index_t n)
585  {
586  index_t i;
587  for (i = 0; n != 0; i++)
588  n >>= 1;
589 
590  return i;
591  }
592 
593  static inline float64_t sin(float64_t x)
594  {
595  return ::sin(x);
596  }
597 
600 
601  static inline float64_t asin(float64_t x)
602  {
603  return ::asin(x);
604  }
605 
608 
609  static inline float64_t sinh(float64_t x)
610  {
611  return ::sinh(x);
612  }
613 
616 
617  static inline float64_t cos(float64_t x)
618  {
619  return ::cos(x);
620  }
621 
624 
625  static inline float64_t acos(float64_t x)
626  {
627  return ::acos(x);
628  }
629 
632 
633  static inline float64_t cosh(float64_t x)
634  {
635  return ::cosh(x);
636  }
637 
640 
641  static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
642  {
643  ASSERT(len>0 && xy)
644 
645  float64_t area = 0.0;
646 
647  if (!reversed)
648  {
649  for (int i=1; i<len; i++)
650  area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
651  }
652  else
653  {
654  for (int i=1; i<len; i++)
655  area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
656  }
657 
658  return area;
659  }
660 
661  static bool strtof(const char* str, float32_t* float_result);
662  static bool strtod(const char* str, float64_t* double_result);
663  static bool strtold(const char* str, floatmax_t* long_double_result);
664 
665  static inline int64_t factorial(int32_t n)
666  {
667  int64_t res=1;
668  for (int i=2; i<=n; i++)
669  res*=i ;
670  return res ;
671  }
672 
673  static void init_random(uint32_t initseed=0)
674  {
675  if (initseed==0)
677  else
678  seed=initseed;
679 
680  sg_rand->set_seed(seed);
681  }
682 
683  static inline uint64_t random()
684  {
685  return sg_rand->random_64();
686  }
687 
688  static inline uint64_t random(uint64_t min_value, uint64_t max_value)
689  {
690  return sg_rand->random(min_value, max_value);
691  }
692 
693  static inline int64_t random(int64_t min_value, int64_t max_value)
694  {
695  return sg_rand->random(min_value, max_value);
696  }
697 
698  static inline uint32_t random(uint32_t min_value, uint32_t max_value)
699  {
700  return sg_rand->random(min_value, max_value);
701  }
702 
703  static inline int32_t random(int32_t min_value, int32_t max_value)
704  {
705  return sg_rand->random(min_value, max_value);
706  }
707 
708  static inline float32_t random(float32_t min_value, float32_t max_value)
709  {
710  return sg_rand->random(min_value, max_value);
711  }
712 
713  static inline float64_t random(float64_t min_value, float64_t max_value)
714  {
715  return sg_rand->random(min_value, max_value);
716  }
717 
718  static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
719  {
720  return sg_rand->random(min_value, max_value);
721  }
722 
726  static inline float32_t normal_random(float32_t mean, float32_t std_dev)
727  {
728  // sets up variables & makes sure rand_s.range == (0,1)
729  float32_t ret;
730  float32_t rand_u;
731  float32_t rand_v;
732  float32_t rand_s;
733  do
734  {
735  rand_u = CMath::random(-1.0, 1.0);
736  rand_v = CMath::random(-1.0, 1.0);
737  rand_s = rand_u*rand_u + rand_v*rand_v;
738  } while ((rand_s == 0) || (rand_s >= 1));
739 
740  // the meat & potatos, and then the mean & standard deviation shifting...
741  ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
742  ret = std_dev*ret + mean;
743  return ret;
744  }
745 
749  static inline float64_t normal_random(float64_t mean, float64_t std_dev)
750  {
751  return sg_rand->normal_distrib(mean, std_dev);
752  }
753 
756  static inline float32_t randn_float()
757  {
758  return normal_random(0.0, 1.0);
759  }
760 
763  static inline float64_t randn_double()
764  {
765  return sg_rand->std_normal_distrib();
766  }
767 
774  template <class T>
775  static void permute(SGVector<T> v, CRandom* rand=NULL)
776  {
777  if (rand)
778  {
779  for (index_t i=0; i<v.vlen; ++i)
780  swap(v[i], v[rand->random(i, v.vlen-1)]);
781  }
782  else
783  {
784  for (index_t i=0; i<v.vlen; ++i)
785  swap(v[i], v[random(i, v.vlen-1)]);
786  }
787  }
788 
789  template <class T>
790  static int32_t get_num_nonzero(T* vec, int32_t len)
791  {
792  int32_t nnz = 0;
793  for (index_t i=0; i<len; ++i)
794  nnz += vec[i] != 0;
795 
796  return nnz;
797  }
798 
799  static int32_t get_num_nonzero(complex128_t* vec, int32_t len)
800  {
801  int32_t nnz=0;
802  for (index_t i=0; i<len; ++i)
803  nnz+=vec[i]!=0.0;
804  return nnz;
805  }
806 
807  static inline int64_t nchoosek(int32_t n, int32_t k)
808  {
809  int64_t res=1;
810 
811  for (int32_t i=n-k+1; i<=n; i++)
812  res*=i;
813 
814  return res/factorial(k);
815  }
816 
824  static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100);
825 
833  template <class T>
834  static float64_t* linspace(T start, T end, int32_t n)
835  {
836  float64_t* output = SG_MALLOC(float64_t, n);
837  linspace(output, start, end, n);
838 
839  return output;
840  }
841 
849  template <class T>
850  static SGVector<float64_t> linspace_vec(T start, T end, int32_t n)
851  {
852  return SGVector<float64_t>(linspace(start, end, n), n);
853  }
854 
861  template <class T>
862  static T log_sum_exp(SGVector<T> values)
863  {
864  REQUIRE(values.vector, "Values are empty");
865  REQUIRE(values.vlen>0,"number of values supplied is 0\n");
866 
867  /* find minimum element index */
868  index_t min_index=0;
869  T X0=values[0];
870  if (values.vlen>1)
871  {
872  for (index_t i=1; i<values.vlen; ++i)
873  {
874  if (values[i]<X0)
875  {
876  X0=values[i];
877  min_index=i;
878  }
879  }
880  }
881 
882  /* remove element from vector copy and compute log sum exp */
883  SGVector<T> values_without_X0(values.vlen-1);
884  index_t from_idx=0;
885  index_t to_idx=0;
886  for (from_idx=0; from_idx<values.vlen; ++from_idx)
887  {
888  if (from_idx!=min_index)
889  {
890  values_without_X0[to_idx]=exp(values[from_idx]-X0);
891  to_idx++;
892  }
893  }
894 
895  return X0+log(SGVector<T>::sum(values_without_X0)+1);
896  }
897 
904  template <class T>
905  static T log_mean_exp(SGVector<T> values)
906  {
907  return log_sum_exp(values) - log(values.vlen);
908  }
909 
913  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
914  static void sort(float64_t *a, int32_t*idx, int32_t N);
915 
918  template <class T>
919  static void qsort(T* output, int32_t size)
920  {
921  if (size<=1)
922  return;
923 
924  if (size==2)
925  {
926  if (output[0] > output [1])
927  CMath::swap(output[0],output[1]);
928  return;
929  }
930  //T split=output[random(0,size-1)];
931  T split=output[size/2];
932 
933  int32_t left=0;
934  int32_t right=size-1;
935 
936  while (left<=right)
937  {
938  while (output[left] < split)
939  left++;
940  while (output[right] > split)
941  right--;
942 
943  if (left<=right)
944  {
945  CMath::swap(output[left],output[right]);
946  left++;
947  right--;
948  }
949  }
950 
951  if (right+1> 1)
952  qsort(output,right+1);
953 
954  if (size-left> 1)
955  qsort(&output[left],size-left);
956  }
957 
960  template <class T>
961  static void insertion_sort(T* output, int32_t size)
962  {
963  for (int32_t i=0; i<size-1; i++)
964  {
965  int32_t j=i-1;
966  T value=output[i];
967  while (j >= 0 && output[j] > value)
968  {
969  output[j+1] = output[j];
970  j--;
971  }
972  output[j+1]=value;
973  }
974  }
975 
977  template <class T>
978  inline static void radix_sort(T* array, int32_t size)
979  {
980  radix_sort_helper(array,size,0);
981  }
982 
983  /*
984  * Inline function to extract the byte at position p (from left)
985  * of an 64 bit integer. The function is somewhat identical to
986  * accessing an array of characters via [].
987  */
988  template <class T>
989  static inline uint8_t byte(T word, uint16_t p)
990  {
991  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
992  }
993 
995  static inline uint8_t byte(complex128_t word, uint16_t p)
996  {
997  SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
998  return uint8_t(0);
999  }
1000 
1001  template <class T>
1002  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
1003  {
1004  static size_t count[256], nc, cmin;
1005  T *ak;
1006  uint8_t c=0;
1007  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
1008  T *an, *aj, *pile[256];
1009  size_t *cp, cmax;
1010 
1011  /* Push initial array to stack */
1012  sp = s;
1013  radix_push(array, size, i);
1014 
1015  /* Loop until all digits have been sorted */
1016  while (sp>s) {
1017  radix_pop(array, size, i);
1018  an = array + size;
1019 
1020  /* Make character histogram */
1021  if (nc == 0) {
1022  cmin = 0xff;
1023  for (ak = array; ak < an; ak++) {
1024  c = byte(*ak, i);
1025  count[c]++;
1026  if (count[c] == 1) {
1027  /* Determine smallest character */
1028  if (c < cmin)
1029  cmin = c;
1030  nc++;
1031  }
1032  }
1033 
1034  /* Sort recursively if stack size too small */
1035  if (sp + nc > s + RADIX_STACK_SIZE) {
1036  radix_sort_helper(array, size, i);
1037  continue;
1038  }
1039  }
1040 
1041  /*
1042  * Set pile[]; push incompletely sorted bins onto stack.
1043  * pile[] = pointers to last out-of-place element in bins.
1044  * Before permuting: pile[c-1] + count[c] = pile[c];
1045  * during deal: pile[c] counts down to pile[c-1].
1046  */
1047  olds = bigs = sp;
1048  cmax = 2;
1049  ak = array;
1050  pile[0xff] = an;
1051  for (cp = count + cmin; nc > 0; cp++) {
1052  /* Find next non-empty pile */
1053  while (*cp == 0)
1054  cp++;
1055  /* Pile with several entries */
1056  if (*cp > 1) {
1057  /* Determine biggest pile */
1058  if (*cp > cmax) {
1059  cmax = *cp;
1060  bigs = sp;
1061  }
1062 
1063  if (i < sizeof(T)-1)
1064  radix_push(ak, *cp, (uint16_t) (i + 1));
1065  }
1066  pile[cp - count] = ak += *cp;
1067  nc--;
1068  }
1069 
1070  /* Play it safe -- biggest bin last. */
1071  swap(*olds, *bigs);
1072 
1073  /*
1074  * Permute misplacements home. Already home: everything
1075  * before aj, and in pile[c], items from pile[c] on.
1076  * Inner loop:
1077  * r = next element to put in place;
1078  * ak = pile[r[i]] = location to put the next element.
1079  * aj = bottom of 1st disordered bin.
1080  * Outer loop:
1081  * Once the 1st disordered bin is done, ie. aj >= ak,
1082  * aj<-aj + count[c] connects the bins in array linked list;
1083  * reset count[c].
1084  */
1085  aj = array;
1086  while (aj<an)
1087  {
1088  T r;
1089 
1090  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
1091  swap(*ak, r);
1092 
1093  *aj = r;
1094  aj += count[c];
1095  count[c] = 0;
1096  }
1097  }
1098  }
1099 
1101  static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
1102  {
1103  SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
1104  }
1105 
1115  template <class T>
1116  static void qsort(T** vector, index_t length)
1117  {
1118  if (length<=1)
1119  return;
1120 
1121  if (length==2)
1122  {
1123  if (*vector[0]>*vector[1])
1124  swap(vector[0],vector[1]);
1125  return;
1126  }
1127  T* split=vector[length/2];
1128 
1129  int32_t left=0;
1130  int32_t right=length-1;
1131 
1132  while (left<=right)
1133  {
1134  while (*vector[left]<*split)
1135  ++left;
1136  while (*vector[right]>*split)
1137  --right;
1138 
1139  if (left<=right)
1140  {
1141  swap(vector[left],vector[right]);
1142  ++left;
1143  --right;
1144  }
1145  }
1146 
1147  if (right+1>1)
1148  qsort(vector,right+1);
1149 
1150  if (length-left>1)
1151  qsort(&vector[left],length-left);
1152  }
1153 
1155  static void qsort(complex128_t** vector, index_t length)
1156  {
1157  SG_SERROR("CMath::qsort():: Not supported for complex128_t\n");
1158  }
1159 
1161  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1162  {
1163  ASSERT(width>=0)
1164  for (int i=0; i<width; i++)
1165  {
1166  T mask = ((T) 1)<<(sizeof(T)*8-1);
1167  while (mask)
1168  {
1169  if (mask & word)
1170  SG_SPRINT("1")
1171  else
1172  SG_SPRINT("0")
1173 
1174  mask>>=1;
1175  }
1176  }
1177  }
1178 
1180  static void display_bits(complex128_t word,
1181  int32_t width=8*sizeof(complex128_t))
1182  {
1183  SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
1184  }
1185 
1191  template <class T1,class T2>
1192  static void qsort_index(T1* output, T2* index, uint32_t size);
1193 
1195  template <class T>
1196  static void qsort_index(complex128_t* output, T* index, uint32_t size)
1197  {
1198  SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
1199  }
1200 
1206  template <class T1,class T2>
1207  static void qsort_backward_index(
1208  T1* output, T2* index, int32_t size);
1209 
1211  template <class T>
1213  complex128_t* output, T* index, uint32_t size)
1214  {
1215  SG_SERROR("CMath::qsort_backword_index():: \
1216  Not supported for complex128_t\n");
1217  }
1218 
1226  template <class T1,class T2>
1227  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1228  {
1229  int32_t n=0;
1230  thread_qsort<T1,T2> t;
1231  t.output=output;
1232  t.index=index;
1233  t.size=size;
1234  t.qsort_threads=&n;
1235  t.sort_limit=limit;
1236  t.num_threads=n_threads;
1237  parallel_qsort_index<T1,T2>(&t);
1238  }
1239 
1241  template <class T>
1242  inline static void parallel_qsort_index(complex128_t* output, T* index,
1243  uint32_t size, int32_t n_threads, int32_t limit=0)
1244  {
1245  SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1246  }
1247 
1248 
1249  template <class T1,class T2>
1250  static void* parallel_qsort_index(void* p);
1251 
1252 
1253  /* finds the smallest element in output and puts that element as the
1254  first element */
1255  template <class T>
1256  static void min(float64_t* output, T* index, int32_t size);
1257 
1259  static void min(float64_t* output, complex128_t* index, int32_t size)
1260  {
1261  SG_SERROR("CMath::min():: Not supported for complex128_t\n");
1262  }
1263 
1264  /* finds the n smallest elements in output and puts these elements as the
1265  first n elements */
1266  template <class T>
1267  static void nmin(
1268  float64_t* output, T* index, int32_t size, int32_t n);
1269 
1271  static void nmin(float64_t* output, complex128_t* index,
1272  int32_t size, int32_t n)
1273  {
1274  SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
1275  }
1276 
1277 
1278 
1279  /* finds an element in a sorted array via binary search
1280  * returns -1 if not found */
1281  template <class T>
1282  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1283  {
1284  int32_t start=0;
1285  int32_t end=size-1;
1286 
1287  if (size<1)
1288  return 0;
1289 
1290  while (start<end)
1291  {
1292  int32_t middle=(start+end)/2;
1293 
1294  if (output[middle]>elem)
1295  end=middle-1;
1296  else if (output[middle]<elem)
1297  start=middle+1;
1298  else
1299  return middle;
1300  }
1301 
1302  return start;
1303  }
1304 
1306  static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
1307  {
1308  SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
1309  return int32_t(0);
1310  }
1311 
1312  /* finds an element in a sorted array via binary search
1313  * * returns -1 if not found */
1314  template <class T>
1315  static inline int32_t binary_search(T* output, int32_t size, T elem)
1316  {
1317  int32_t ind = binary_search_helper(output, size, elem);
1318  if (ind >= 0 && output[ind] == elem)
1319  return ind;
1320  return -1;
1321  }
1322 
1324  static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
1325  {
1326  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1327  return int32_t(-1);
1328  }
1329 
1330 
1331  /* Finds an element in a sorted array of pointers via binary search
1332  * Every element is dereferenced once before being compared
1333  *
1334  * @param array array of pointers to search in (assumed being sorted)
1335  * @param length length of array
1336  * @param elem pointer to element to search for
1337  * @return index of elem, -1 if not found */
1338  template<class T>
1339  static inline int32_t binary_search(T** vector, index_t length,
1340  T* elem)
1341  {
1342  int32_t start=0;
1343  int32_t end=length-1;
1344 
1345  if (length<1)
1346  return -1;
1347 
1348  while (start<end)
1349  {
1350  int32_t middle=(start+end)/2;
1351 
1352  if (*vector[middle]>*elem)
1353  end=middle-1;
1354  else if (*vector[middle]<*elem)
1355  start=middle+1;
1356  else
1357  {
1358  start=middle;
1359  break;
1360  }
1361  }
1362 
1363  if (start>=0&&*vector[start]==*elem)
1364  return start;
1365 
1366  return -1;
1367  }
1368 
1370  static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
1371  {
1372  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1373  return int32_t(-1);
1374  }
1375 
1376  template <class T>
1378  T* output, int32_t size, T elem)
1379  {
1380  int32_t ind = binary_search_helper(output, size, elem);
1381 
1382  if (output[ind]<=elem)
1383  return ind;
1384  if (ind>0 && output[ind-1] <= elem)
1385  return ind-1;
1386  return -1;
1387  }
1388 
1391  int32_t size, complex128_t elem)
1392  {
1393  SG_SERROR("CMath::binary_search_max_lower_equal():: \
1394  Not supported for complex128_t\n");
1395  return int32_t(-1);
1396  }
1397 
1400  static float64_t Align(
1401  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1402 
1403 
1405 
1408  {
1409  return c.real();
1410  }
1411 
1414  {
1415  return c.imag();
1416  }
1417 
1419  inline static uint32_t get_seed()
1420  {
1421  return CMath::seed;
1422  }
1423 
1425  inline static uint32_t get_log_range()
1426  {
1427  return CMath::LOGRANGE;
1428  }
1429 
1430 #ifdef USE_LOGCACHE
1431  inline static uint32_t get_log_accuracy()
1433  {
1434  return CMath::LOGACCURACY;
1435  }
1436 #endif
1437 
1439  static int is_finite(double f);
1440 
1442  static int is_infinity(double f);
1443 
1445  static int is_nan(double f);
1446 
1457 #ifdef USE_LOGCACHE
1458  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1459  {
1460  float64_t diff;
1461 
1462  if (!CMath::is_finite(p))
1463  return q;
1464 
1465  if (!CMath::is_finite(q))
1466  {
1467  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
1468  return NOT_A_NUMBER;
1469  }
1470  diff = p - q;
1471  if (diff > 0)
1472  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1473  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1474  }
1475 
1477  static void init_log_table();
1478 
1480  static int32_t determine_logrange();
1481 
1483  static int32_t determine_logaccuracy(int32_t range);
1484 #else
1486  float64_t p, float64_t q)
1487  {
1488  float64_t diff;
1489 
1490  if (!CMath::is_finite(p))
1491  return q;
1492  if (!CMath::is_finite(q))
1493  return p;
1494  diff = p - q;
1495  if (diff > 0)
1496  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
1497  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
1498  }
1499 #endif
1500 #ifdef USE_LOGSUMARRAY
1501 
1506  static inline float64_t logarithmic_sum_array(
1507  float64_t *p, int32_t len)
1508  {
1509  if (len<=2)
1510  {
1511  if (len==2)
1512  return logarithmic_sum(p[0],p[1]) ;
1513  if (len==1)
1514  return p[0];
1515  return -INFTY ;
1516  }
1517  else
1518  {
1519  float64_t *pp=p ;
1520  if (len%2==1) pp++ ;
1521  for (int32_t j=0; j < len>>1; j++)
1522  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
1523  }
1524  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1525  }
1526 #endif //USE_LOGSUMARRAY
1527 
1528 
1530  virtual const char* get_name() const { return "Math"; }
1531  public:
1534  static const float64_t NOT_A_NUMBER;
1537  static const float64_t INFTY;
1538  static const float64_t ALMOST_INFTY;
1539 
1542 
1544  static const float64_t PI;
1545 
1548 
1549  /* largest and smallest possible float64_t */
1552 
1553  /* Floating point Limits, Normalized */
1554  static const float32_t F_MAX_VAL32;
1556  static const float64_t F_MAX_VAL64;
1558 
1559  /* Floating point limits, Denormalized */
1560  static const float32_t F_MIN_VAL32;
1561  static const float64_t F_MIN_VAL64;
1562 
1563  protected:
1565  static int32_t LOGRANGE;
1566 
1568  static uint32_t seed;
1569 
1570 #ifdef USE_LOGCACHE
1571 
1573  static int32_t LOGACCURACY;
1575  static float64_t* logtable;
1577 #endif
1578 };
1579 
1580 //implementations of template functions
1581 template <class T1,class T2>
1583  {
1584  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1585  T1* output=ps->output;
1586  T2* index=ps->index;
1587  int32_t size=ps->size;
1588  int32_t* qsort_threads=ps->qsort_threads;
1589  int32_t sort_limit=ps->sort_limit;
1590  int32_t num_threads=ps->num_threads;
1591 
1592  if (size<2)
1593  {
1594  return NULL;
1595  }
1596 
1597  if (size==2)
1598  {
1599  if (output[0] > output [1])
1600  {
1601  swap(output[0], output[1]);
1602  swap(index[0], index[1]);
1603  }
1604  return NULL;
1605  }
1606  //T1 split=output[random(0,size-1)];
1607  T1 split=output[size/2];
1608 
1609  int32_t left=0;
1610  int32_t right=size-1;
1611 
1612  while (left<=right)
1613  {
1614  while (output[left] < split)
1615  left++;
1616  while (output[right] > split)
1617  right--;
1618 
1619  if (left<=right)
1620  {
1621  swap(output[left], output[right]);
1622  swap(index[left], index[right]);
1623  left++;
1624  right--;
1625  }
1626  }
1627  bool lthread_start=false;
1628  bool rthread_start=false;
1629  pthread_t lthread;
1630  pthread_t rthread;
1631  struct thread_qsort<T1,T2> t1;
1632  struct thread_qsort<T1,T2> t2;
1633 
1634  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1635  qsort_index(output,index,right+1);
1636  else if (right+1> 1)
1637  {
1638  (*qsort_threads)++;
1639  lthread_start=true;
1640  t1.output=output;
1641  t1.index=index;
1642  t1.size=right+1;
1643  t1.qsort_threads=qsort_threads;
1644  t1.sort_limit=sort_limit;
1645  t1.num_threads=num_threads;
1646  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1647  {
1648  lthread_start=false;
1649  (*qsort_threads)--;
1650  qsort_index(output,index,right+1);
1651  }
1652  }
1653 
1654 
1655  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1656  qsort_index(&output[left],&index[left], size-left);
1657  else if (size-left> 1)
1658  {
1659  (*qsort_threads)++;
1660  rthread_start=true;
1661  t2.output=&output[left];
1662  t2.index=&index[left];
1663  t2.size=size-left;
1664  t2.qsort_threads=qsort_threads;
1665  t2.sort_limit=sort_limit;
1666  t2.num_threads=num_threads;
1667  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1668  {
1669  rthread_start=false;
1670  (*qsort_threads)--;
1671  qsort_index(&output[left],&index[left], size-left);
1672  }
1673  }
1674 
1675  if (lthread_start)
1676  {
1677  pthread_join(lthread, NULL);
1678  (*qsort_threads)--;
1679  }
1680 
1681  if (rthread_start)
1682  {
1683  pthread_join(rthread, NULL);
1684  (*qsort_threads)--;
1685  }
1686 
1687  return NULL;
1688  }
1689 
1690  template <class T1,class T2>
1691 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
1692 {
1693  if (size<=1)
1694  return;
1695 
1696  if (size==2)
1697  {
1698  if (output[0] > output [1])
1699  {
1700  swap(output[0],output[1]);
1701  swap(index[0],index[1]);
1702  }
1703  return;
1704  }
1705  //T1 split=output[random(0,size-1)];
1706  T1 split=output[size/2];
1707 
1708  int32_t left=0;
1709  int32_t right=size-1;
1710 
1711  while (left<=right)
1712  {
1713  while (output[left] < split)
1714  left++;
1715  while (output[right] > split)
1716  right--;
1717 
1718  if (left<=right)
1719  {
1720  swap(output[left],output[right]);
1721  swap(index[left],index[right]);
1722  left++;
1723  right--;
1724  }
1725  }
1726 
1727  if (right+1> 1)
1728  qsort_index(output,index,right+1);
1729 
1730  if (size-left> 1)
1731  qsort_index(&output[left],&index[left], size-left);
1732 }
1733 
1734  template <class T1,class T2>
1735 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
1736 {
1737  if (size<=1)
1738  return;
1739 
1740  if (size==2)
1741  {
1742  if (output[0] < output [1])
1743  {
1744  swap(output[0],output[1]);
1745  swap(index[0],index[1]);
1746  }
1747  return;
1748  }
1749 
1750  //T1 split=output[random(0,size-1)];
1751  T1 split=output[size/2];
1752 
1753  int32_t left=0;
1754  int32_t right=size-1;
1755 
1756  while (left<=right)
1757  {
1758  while (output[left] > split)
1759  left++;
1760  while (output[right] < split)
1761  right--;
1762 
1763  if (left<=right)
1764  {
1765  swap(output[left],output[right]);
1766  swap(index[left],index[right]);
1767  left++;
1768  right--;
1769  }
1770  }
1771 
1772  if (right+1> 1)
1773  qsort_backward_index(output,index,right+1);
1774 
1775  if (size-left> 1)
1776  qsort_backward_index(&output[left],&index[left], size-left);
1777 }
1778 
1779  template <class T>
1780 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
1781 {
1782  if (6*n*size<13*size*CMath::log(size))
1783  for (int32_t i=0; i<n; i++)
1784  min(&output[i], &index[i], size-i);
1785  else
1786  qsort_index(output, index, size);
1787 }
1788 
1789 /* move the smallest entry in the array to the beginning */
1790  template <class T>
1791 void CMath::min(float64_t* output, T* index, int32_t size)
1792 {
1793  if (size<=1)
1794  return;
1795  float64_t min_elem=output[0];
1796  int32_t min_index=0;
1797  for (int32_t i=1; i<size; i++)
1798  {
1799  if (output[i]<min_elem)
1800  {
1801  min_index=i;
1802  min_elem=output[i];
1803  }
1804  }
1805  swap(output[0], output[min_index]);
1806  swap(index[0], index[min_index]);
1807 }
1808 
1809 template <>
1810 inline float64_t* CMath::linspace<complex128_t>(complex128_t start, complex128_t end, int32_t n)
1811 {
1812  SG_SERROR("SGVector::linspace():: Not supported for complex128_t\n");
1813  return NULL;
1814 }
1815 
1816 #define COMPLEX128_ERROR_ONEARG_T(function) \
1817 template <> \
1818 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
1819 { \
1820  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1821  #function);\
1822  return complex128_t(0.0, 0.0); \
1823 }
1824 
1825 #define COMPLEX128_ERROR_TWOARGS_T(function) \
1826 template <> \
1827 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \
1828 { \
1829  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1830  #function);\
1831  return complex128_t(0.0, 0.0); \
1832 }
1833 
1834 #define COMPLEX128_ERROR_THREEARGS_T(function) \
1835 template <> \
1836 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \
1837 { \
1838  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1839  #function);\
1840  return complex128_t(0.0, 0.0); \
1841 }
1842 
1843 #define COMPLEX128_ERROR_SORT_T(function) \
1844 template <> \
1845 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \
1846 { \
1847  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1848  #function);\
1849 }
1850 
1851 #define COMPLEX128_ERROR_ARG_MAX_MIN(function) \
1852 template <> \
1853 inline int32_t CMath::function<complex128_t>(complex128_t * a, int32_t b, int32_t c, complex128_t * d) \
1854 { \
1855  int32_t maxIdx=0; \
1856  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1857  #function);\
1858  return maxIdx; \
1859 }
1860 
1863 
1866 
1869 
1871 // COMPLEX128_ERROR_ONEARG_T(sign)
1872 
1875 
1878 
1881 
1884 
1887 
1888 }
1889 #undef COMPLEX128_ERROR_ONEARG
1890 #undef COMPLEX128_ERROR_ONEARG_T
1891 #undef COMPLEX128_ERROR_TWOARGS_T
1892 #undef COMPLEX128_ERROR_THREEARGS_T
1893 #undef COMPLEX128_STDMATH
1894 #undef COMPLEX128_ERROR_SORT_T
1895 #endif
static float64_t sin(float64_t x)
Definition: Math.h:593
static float64_t normal_random(float64_t mean, float64_t std_dev)
Definition: Math.h:749
static const float32_t F_MAX_VAL32
Definition: Math.h:1554
static const float64_t MACHINE_EPSILON
Definition: Math.h:1547
static bool strtof(const char *str, float32_t *float_result)
Definition: Math.cpp:264
static void permute(SGVector< T > v, CRandom *rand=NULL)
Definition: Math.h:775
static int32_t binary_search(complex128_t *output, int32_t size, complex128_t elem)
binary_search not implemented for complex128_t
Definition: Math.h:1324
static void parallel_qsort_index(T1 *output, T2 *index, uint32_t size, int32_t n_threads, int32_t limit=262144)
Definition: Math.h:1227
static uint32_t seed
random generator seed
Definition: Math.h:1568
std::complex< float64_t > complex128_t
Definition: common.h:67
static int is_finite(double f)
checks whether a float is finite
Definition: Math.cpp:249
static float64_t Align(char *seq1, char *seq2, int32_t l1, int32_t l2, float64_t gapCost)
Definition: Math.cpp:161
static int32_t arg_max(T *vec, int32_t inc, int32_t len, T *maxv_ptr=NULL)
Definition: Math.h:233
static float64_t sqrt(float64_t x)
x^0.5
Definition: Math.h:407
static floatmax_t random(floatmax_t min_value, floatmax_t max_value)
Definition: Math.h:718
static floatmax_t sqrt(floatmax_t x)
x^0.5
Definition: Math.h:413
static void linspace(float64_t *output, float64_t start, float64_t end, int32_t n=100)
Definition: Math.cpp:204
int32_t index_t
Definition: common.h:62
static float64_t ceil(float64_t d)
Definition: Math.h:370
static bool strtod(const char *str, float64_t *double_result)
Definition: Math.cpp:295
static complex128_t pow(complex128_t x, int32_t n)
x^n, x or n being a complex128_t
Definition: Math.h:495
virtual ~CMath()
Destructor - frees logtable.
Definition: Math.cpp:83
static const float64_t INFTY
infinity
Definition: Math.h:1537
static SGVector< float64_t > linspace_vec(T start, T end, int32_t n)
Definition: Math.h:850
static int32_t binary_search_helper(T *output, int32_t size, T elem)
Definition: Math.h:1282
#define COMPLEX128_ERROR_THREEARGS_T(function)
Definition: Math.h:1834
static void nmin(float64_t *output, T *index, int32_t size, int32_t n)
Definition: Math.h:1780
static void qsort_index(T1 *output, T2 *index, uint32_t size)
Definition: Math.h:1691
static void qsort_index(complex128_t *output, T *index, uint32_t size)
qsort_index not implemented for complex128_t
Definition: Math.h:1196
static float64_t log10(float64_t v)
tanh(x), x being a complex128_t
Definition: Math.h:559
#define COMPLEX128_ERROR_TWOARGS_T(function)
Definition: Math.h:1825
#define SG_SWARNING(...)
Definition: SGIO.h:178
static float32_t normal_random(float32_t mean, float32_t std_dev)
Definition: Math.h:726
static float64_t random(float64_t min_value, float64_t max_value)
Definition: Math.h:713
static int32_t binary_search(T *output, int32_t size, T elem)
Definition: Math.h:1315
static T sq(T x)
x^2
Definition: Math.h:395
static uint32_t get_seed()
returns number generator seed
Definition: Math.h:1419
#define COMPLEX128_ERROR_ARG_MAX_MIN(function)
Definition: Math.h:1851
static float32_t randn_float()
Definition: Math.h:756
static const float64_t MIN_REAL_NUMBER
Definition: Math.h:1551
static float64_t randn_double()
Definition: Math.h:763
static int32_t binary_search_max_lower_equal(T *output, int32_t size, T elem)
Definition: Math.h:1377
static float64_t atan(float64_t x)
tan(x), x being a complex128_t
Definition: Math.h:533
#define REQUIRE(x,...)
Definition: SGIO.h:206
static const float64_t F_MAX_VAL64
Definition: Math.h:1556
static float64_t cosh(float64_t x)
acos(x), x being a complex128_t not implemented
Definition: Math.h:633
static uint32_t get_log_range()
returns range of logtable
Definition: Math.h:1425
void split(v_array< ds_node< P > > &point_set, v_array< ds_node< P > > &far_set, int max_scale)
Definition: JLCoverTree.h:149
static float32_t random(float32_t min_value, float32_t max_value)
Definition: Math.h:708
static const float32_t F_MIN_VAL32
Definition: Math.h:1560
CRandom * sg_rand
Definition: init.cpp:39
static bool fequals_abs(const T &a, const T &b, const float64_t eps)
Definition: Math.h:284
static float64_t get_abs_tolorance(float64_t true_value, float64_t rel_tolorance)
Definition: Math.cpp:364
#define RADIX_STACK_SIZE
Definition: Math.h:76
static float64_t imag(complex128_t c)
returns imag part of a complex128_t number
Definition: Math.h:1413
static float64_t floor(float64_t d)
Definition: Math.h:365
static uint64_t random()
Definition: Math.h:683
static float64_t * linspace(T start, T end, int32_t n)
Definition: Math.h:834
static int32_t get_num_nonzero(complex128_t *vec, int32_t len)
Definition: Math.h:799
static const float32_t F_MIN_NORM_VAL32
Definition: Math.h:1555
static float64_t real(complex128_t c)
returns real part of a complex128_t number
Definition: Math.h:1407
static uint8_t byte(complex128_t word, uint16_t p)
byte not implemented for complex128_t
Definition: Math.h:995
static void qsort(T *output, int32_t size)
Definition: Math.h:919
static int32_t LOGRANGE
range for logtable: log(1+exp(x)) -LOGRANGE <= x <= 0
Definition: Math.h:1565
static const float64_t ALMOST_NEG_INFTY
almost neg (log) infinity
Definition: Math.h:1541
static bool fequals(const T &a, const T &b, const float64_t eps, bool tolerant=false)
Definition: Math.h:301
static complex128_t pow(complex128_t x, complex128_t n)
Definition: Math.h:500
static float32_t invsqrt(float32_t x)
x^0.5, x being a complex128_t
Definition: Math.h:428
index_t vlen
Definition: SGVector.h:637
static uint8_t byte(T word, uint16_t p)
Definition: Math.h:989
static complex128_t pow(float64_t x, complex128_t n)
Definition: Math.h:510
#define SG_SPRINT(...)
Definition: SGIO.h:180
CMath()
Constructor - initializes log-table.
Definition: Math.cpp:66
static float64_t pow(float64_t x, int32_t n)
Definition: Math.h:475
#define ASSERT(x)
Definition: SGIO.h:201
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:112
static int32_t random(int32_t min_value, int32_t max_value)
Definition: Math.h:703
static float64_t pow(float64_t x, float64_t n)
Definition: Math.h:489
shogun vector
Definition: Parameter.h:28
static void qsort(T **vector, index_t length)
Definition: Math.h:1116
static void init_random(uint32_t initseed=0)
Definition: Math.h:673
static int32_t pow(int32_t x, int32_t n)
Definition: Math.h:465
#define radix_pop(a, n, i)
Definition: Math.h:80
double float64_t
Definition: common.h:50
static void parallel_qsort_index(complex128_t *output, T *index, uint32_t size, int32_t n_threads, int32_t limit=0)
parallel_qsort_index not implemented for complex128_t
Definition: Math.h:1242
static void min(float64_t *output, complex128_t *index, int32_t size)
complex128_t cannot be used as index
Definition: Math.h:1259
long double floatmax_t
Definition: common.h:51
#define COMPLEX128_ERROR_ONEARG(function)
Definition: Math.h:115
static float64_t cos(float64_t x)
sinh(x), x being a complex128_t
Definition: Math.h:617
static void qsort_backword_index(complex128_t *output, T *index, uint32_t size)
qsort_backword_index not implemented for complex128_t
Definition: Math.h:1212
static int32_t get_num_nonzero(T *vec, int32_t len)
Definition: Math.h:790
static T log_sum_exp(SGVector< T > values)
Definition: Math.h:862
static T max(T a, T b)
return the maximum of two integers
Definition: Math.h:163
static float64_t area_under_curve(float64_t *xy, int32_t len, bool reversed)
cosh(x), x being a complex128_t
Definition: Math.h:641
static void display_bits(T word, int32_t width=8 *sizeof(T))
display bits (useful for debugging)
Definition: Math.h:1161
static uint64_t random(uint64_t min_value, uint64_t max_value)
Definition: Math.h:688
static int64_t factorial(int32_t n)
Definition: Math.h:665
static void qsort(complex128_t **vector, index_t length)
qsort not implemented for complex128_t
Definition: Math.h:1155
static float64_t atan2(float64_t x, float64_t y)
atan(x), x being a complex128_t not implemented
Definition: Math.h:542
static int32_t arg_min(T *vec, int32_t inc, int32_t len, T *minv_ptr=NULL)
return arg_min(vec)
Definition: Math.h:255
static float64_t tan(float64_t x)
exp(x), x being a complex128_t
Definition: Math.h:524
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:1390
float float32_t
Definition: common.h:49
static float64_t sinh(float64_t x)
asin(x), x being a complex128_t not implemented
Definition: Math.h:609
: Pseudo random number geneartor
Definition: Random.h:34
static int64_t nchoosek(int32_t n, int32_t k)
Definition: Math.h:807
static T min(T *vec, int32_t len)
return the smallest element in a vector.
Definition: Math.h:171
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:1339
static int is_infinity(double f)
checks whether a float is infinity
Definition: Math.cpp:230
static float64_t acos(float64_t x)
cos(x), x being a complex128_t
Definition: Math.h:625
static float64_t abs(complex128_t a)
return the absolute value of a complex number
Definition: Math.h:223
static float64_t log2(float64_t v)
log10(x), x being a complex128_t
Definition: Math.h:567
static void display_bits(complex128_t word, int32_t width=8 *sizeof(complex128_t))
disply_bits not implemented for complex128_t
Definition: Math.h:1180
#define COMPLEX128_STDMATH(function)
Definition: Math.h:123
static void radix_sort_helper(T *array, int32_t size, uint16_t i)
Definition: Math.h:1002
static T sign(T a)
signum of type T variable a
Definition: Math.h:377
static float64_t tanh(float64_t x)
atan2(x), x being a complex128_t not implemented
Definition: Math.h:551
static int is_nan(double f)
checks whether a float is nan
Definition: Math.cpp:217
virtual const char * get_name() const
Definition: Math.h:1530
static float64_t asin(float64_t x)
sin(x), x being a complex128_t
Definition: Math.h:601
static T min(T a, T b)
return the minimum of two integers
Definition: Math.h:156
#define SG_SERROR(...)
Definition: SGIO.h:179
static float64_t exp(float64_t x)
Definition: Math.h:515
static const float64_t F_MIN_VAL64
Definition: Math.h:1561
static const float64_t F_MIN_NORM_VAL64
Definition: Math.h:1557
static float64_t log(float64_t v)
Definition: Math.h:576
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:1101
static const float64_t ALMOST_INFTY
Definition: Math.h:1538
Class which collects generic mathematical functions.
Definition: Math.h:136
static T log_mean_exp(SGVector< T > values)
Definition: Math.h:905
static void insertion_sort(T *output, int32_t size)
Definition: Math.h:961
static void swap(T &a, T &b)
swap e.g. floats a and b
Definition: Math.h:386
static complex128_t pow(complex128_t x, float64_t n)
Definition: Math.h:505
static void sort(int32_t *a, int32_t cols, int32_t sort_col=0)
Definition: Math.cpp:122
static T max(T *vec, int32_t len)
return the greatest element in a vector.
Definition: Math.h:185
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:1306
static float64_t round(float64_t d)
Definition: Math.h:360
static float32_t sqrt(float32_t x)
x^0.5
Definition: Math.h:401
static uint32_t generate_seed()
Definition: Random.cpp:315
static index_t floor_log(index_t n)
log(x), x being a complex128_t
Definition: Math.h:584
static void radix_sort(T *array, int32_t size)
Definition: Math.h:978
static uint32_t random(uint32_t min_value, uint32_t max_value)
Definition: Math.h:698
static floatmax_t powl(floatmax_t x, floatmax_t n)
x^n
Definition: Math.h:449
static float64_t logarithmic_sum(float64_t p, float64_t q)
Definition: Math.h:1485
#define radix_push(a, n, i)
Definition: Math.h:79
static T clamp(T value, T lb, T ub)
return the value clamped to interval [lb,ub]
Definition: Math.h:198
static int32_t binary_search(complex128_t **vector, index_t length, complex128_t *elem)
binary_search not implemented for complex128_t
Definition: Math.h:1370
#define COMPLEX128_ERROR_SORT_T(function)
Definition: Math.h:1843
static bool strtold(const char *str, floatmax_t *long_double_result)
Definition: Math.cpp:326
static int32_t pow(bool x, int32_t n)
Definition: Math.h:460
static const float64_t NOT_A_NUMBER
not a number
Definition: Math.h:1535
static void qsort_backward_index(T1 *output, T2 *index, int32_t size)
Definition: Math.h:1735
static const float64_t MAX_REAL_NUMBER
Definition: Math.h:1550
static T abs(T a)
return the absolute value of a number
Definition: Math.h:210
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:1271
static int64_t random(int64_t min_value, int64_t max_value)
Definition: Math.h:693
static const float64_t PI
Definition: Math.h:1544

SHOGUN 机器学习工具包 - 项目文档