SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Math.h
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Thoralf Klein
8  * Written (W) 2013 Soumyajit De
9  * Written (W) 2012 Fernando Jose Iglesias Garcia
10  * Written (W) 2011 Siddharth Kherada
11  * Written (W) 2011 Justin Patera
12  * Written (W) 2011 Alesis Novik
13  * Written (W) 2011-2013 Heiko Strathmann
14  * Written (W) 1999-2009 Soeren Sonnenburg
15  * Written (W) 1999-2008 Gunnar Raetsch
16  * Written (W) 2007 Konrad Rieck
17  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
18  */
19 
20 #ifndef __MATHEMATICS_H_
21 #define __MATHEMATICS_H_
22 
23 #include <shogun/lib/config.h>
24 
25 #include <shogun/base/SGObject.h>
26 #include <shogun/lib/common.h>
27 #include <shogun/io/SGIO.h>
28 #include <shogun/base/Parallel.h>
30 #include <shogun/lib/SGVector.h>
31 
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 
143  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 
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 
832  template <class T>
833  static T log_sum_exp(SGVector<T> values)
834  {
835  REQUIRE(values.vector, "Values are empty");
836  REQUIRE(values.vlen>0,"number of values supplied is 0\n");
837 
838  /* find minimum element index */
839  index_t min_index=0;
840  T X0=values[0];
841  if (values.vlen>1)
842  {
843  for (index_t i=1; i<values.vlen; ++i)
844  {
845  if (values[i]<X0)
846  {
847  X0=values[i];
848  min_index=i;
849  }
850  }
851  }
852 
853  /* remove element from vector copy and compute log sum exp */
854  SGVector<T> values_without_X0(values.vlen-1);
855  index_t from_idx=0;
856  index_t to_idx=0;
857  for (from_idx=0; from_idx<values.vlen; ++from_idx)
858  {
859  if (from_idx!=min_index)
860  {
861  values_without_X0[to_idx]=exp(values[from_idx]-X0);
862  to_idx++;
863  }
864  }
865 
866  return X0+log(SGVector<T>::sum(values_without_X0)+1);
867  }
868 
875  template <class T>
876  static T log_mean_exp(SGVector<T> values)
877  {
878  return log_sum_exp(values) - log(values.vlen);
879  }
880 
884  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
885  static void sort(float64_t *a, int32_t*idx, int32_t N);
886 
889  template <class T>
890  static void qsort(T* output, int32_t size)
891  {
892  if (size<=1)
893  return;
894 
895  if (size==2)
896  {
897  if (output[0] > output [1])
898  CMath::swap(output[0],output[1]);
899  return;
900  }
901  //T split=output[random(0,size-1)];
902  T split=output[size/2];
903 
904  int32_t left=0;
905  int32_t right=size-1;
906 
907  while (left<=right)
908  {
909  while (output[left] < split)
910  left++;
911  while (output[right] > split)
912  right--;
913 
914  if (left<=right)
915  {
916  CMath::swap(output[left],output[right]);
917  left++;
918  right--;
919  }
920  }
921 
922  if (right+1> 1)
923  qsort(output,right+1);
924 
925  if (size-left> 1)
926  qsort(&output[left],size-left);
927  }
928 
931  template <class T>
932  static void insertion_sort(T* output, int32_t size)
933  {
934  for (int32_t i=0; i<size-1; i++)
935  {
936  int32_t j=i-1;
937  T value=output[i];
938  while (j >= 0 && output[j] > value)
939  {
940  output[j+1] = output[j];
941  j--;
942  }
943  output[j+1]=value;
944  }
945  }
946 
948  template <class T>
949  inline static void radix_sort(T* array, int32_t size)
950  {
951  radix_sort_helper(array,size,0);
952  }
953 
954  /*
955  * Inline function to extract the byte at position p (from left)
956  * of an 64 bit integer. The function is somewhat identical to
957  * accessing an array of characters via [].
958  */
959  template <class T>
960  static inline uint8_t byte(T word, uint16_t p)
961  {
962  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
963  }
964 
966  static inline uint8_t byte(complex128_t word, uint16_t p)
967  {
968  SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
969  return uint8_t(0);
970  }
971 
972  template <class T>
973  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
974  {
975  static size_t count[256], nc, cmin;
976  T *ak;
977  uint8_t c=0;
978  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
979  T *an, *aj, *pile[256];
980  size_t *cp, cmax;
981 
982  /* Push initial array to stack */
983  sp = s;
984  radix_push(array, size, i);
985 
986  /* Loop until all digits have been sorted */
987  while (sp>s) {
988  radix_pop(array, size, i);
989  an = array + size;
990 
991  /* Make character histogram */
992  if (nc == 0) {
993  cmin = 0xff;
994  for (ak = array; ak < an; ak++) {
995  c = byte(*ak, i);
996  count[c]++;
997  if (count[c] == 1) {
998  /* Determine smallest character */
999  if (c < cmin)
1000  cmin = c;
1001  nc++;
1002  }
1003  }
1004 
1005  /* Sort recursively if stack size too small */
1006  if (sp + nc > s + RADIX_STACK_SIZE) {
1007  radix_sort_helper(array, size, i);
1008  continue;
1009  }
1010  }
1011 
1012  /*
1013  * Set pile[]; push incompletely sorted bins onto stack.
1014  * pile[] = pointers to last out-of-place element in bins.
1015  * Before permuting: pile[c-1] + count[c] = pile[c];
1016  * during deal: pile[c] counts down to pile[c-1].
1017  */
1018  olds = bigs = sp;
1019  cmax = 2;
1020  ak = array;
1021  pile[0xff] = an;
1022  for (cp = count + cmin; nc > 0; cp++) {
1023  /* Find next non-empty pile */
1024  while (*cp == 0)
1025  cp++;
1026  /* Pile with several entries */
1027  if (*cp > 1) {
1028  /* Determine biggest pile */
1029  if (*cp > cmax) {
1030  cmax = *cp;
1031  bigs = sp;
1032  }
1033 
1034  if (i < sizeof(T)-1)
1035  radix_push(ak, *cp, (uint16_t) (i + 1));
1036  }
1037  pile[cp - count] = ak += *cp;
1038  nc--;
1039  }
1040 
1041  /* Play it safe -- biggest bin last. */
1042  swap(*olds, *bigs);
1043 
1044  /*
1045  * Permute misplacements home. Already home: everything
1046  * before aj, and in pile[c], items from pile[c] on.
1047  * Inner loop:
1048  * r = next element to put in place;
1049  * ak = pile[r[i]] = location to put the next element.
1050  * aj = bottom of 1st disordered bin.
1051  * Outer loop:
1052  * Once the 1st disordered bin is done, ie. aj >= ak,
1053  * aj<-aj + count[c] connects the bins in array linked list;
1054  * reset count[c].
1055  */
1056  aj = array;
1057  while (aj<an)
1058  {
1059  T r;
1060 
1061  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
1062  swap(*ak, r);
1063 
1064  *aj = r;
1065  aj += count[c];
1066  count[c] = 0;
1067  }
1068  }
1069  }
1070 
1072  static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
1073  {
1074  SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
1075  }
1076 
1086  template <class T>
1087  static void qsort(T** vector, index_t length)
1088  {
1089  if (length<=1)
1090  return;
1091 
1092  if (length==2)
1093  {
1094  if (*vector[0]>*vector[1])
1095  swap(vector[0],vector[1]);
1096  return;
1097  }
1098  T* split=vector[length/2];
1099 
1100  int32_t left=0;
1101  int32_t right=length-1;
1102 
1103  while (left<=right)
1104  {
1105  while (*vector[left]<*split)
1106  ++left;
1107  while (*vector[right]>*split)
1108  --right;
1109 
1110  if (left<=right)
1111  {
1112  swap(vector[left],vector[right]);
1113  ++left;
1114  --right;
1115  }
1116  }
1117 
1118  if (right+1>1)
1119  qsort(vector,right+1);
1120 
1121  if (length-left>1)
1122  qsort(&vector[left],length-left);
1123  }
1124 
1126  static void qsort(complex128_t** vector, index_t length)
1127  {
1128  SG_SERROR("CMath::qsort():: Not supported for complex128_t\n");
1129  }
1130 
1132  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1133  {
1134  ASSERT(width>=0)
1135  for (int i=0; i<width; i++)
1136  {
1137  T mask = ((T) 1)<<(sizeof(T)*8-1);
1138  while (mask)
1139  {
1140  if (mask & word)
1141  SG_SPRINT("1")
1142  else
1143  SG_SPRINT("0")
1144 
1145  mask>>=1;
1146  }
1147  }
1148  }
1149 
1151  static void display_bits(complex128_t word,
1152  int32_t width=8*sizeof(complex128_t))
1153  {
1154  SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
1155  }
1156 
1162  template <class T1,class T2>
1163  static void qsort_index(T1* output, T2* index, uint32_t size);
1164 
1166  template <class T>
1167  static void qsort_index(complex128_t* output, T* index, uint32_t size)
1168  {
1169  SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
1170  }
1171 
1177  template <class T1,class T2>
1178  static void qsort_backward_index(
1179  T1* output, T2* index, int32_t size);
1180 
1182  template <class T>
1184  complex128_t* output, T* index, uint32_t size)
1185  {
1186  SG_SERROR("CMath::qsort_backword_index():: \
1187  Not supported for complex128_t\n");
1188  }
1189 
1197  template <class T1,class T2>
1198  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1199  {
1200  int32_t n=0;
1201  thread_qsort<T1,T2> t;
1202  t.output=output;
1203  t.index=index;
1204  t.size=size;
1205  t.qsort_threads=&n;
1206  t.sort_limit=limit;
1207  t.num_threads=n_threads;
1208  parallel_qsort_index<T1,T2>(&t);
1209  }
1210 
1212  template <class T>
1213  inline static void parallel_qsort_index(complex128_t* output, T* index,
1214  uint32_t size, int32_t n_threads, int32_t limit=0)
1215  {
1216  SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1217  }
1218 
1219 
1220  template <class T1,class T2>
1221  static void* parallel_qsort_index(void* p);
1222 
1223 
1224  /* finds the smallest element in output and puts that element as the
1225  first element */
1226  template <class T>
1227  static void min(float64_t* output, T* index, int32_t size);
1228 
1230  static void min(float64_t* output, complex128_t* index, int32_t size)
1231  {
1232  SG_SERROR("CMath::min():: Not supported for complex128_t\n");
1233  }
1234 
1235  /* finds the n smallest elements in output and puts these elements as the
1236  first n elements */
1237  template <class T>
1238  static void nmin(
1239  float64_t* output, T* index, int32_t size, int32_t n);
1240 
1242  static void nmin(float64_t* output, complex128_t* index,
1243  int32_t size, int32_t n)
1244  {
1245  SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
1246  }
1247 
1248 
1249 
1250  /* finds an element in a sorted array via binary search
1251  * returns -1 if not found */
1252  template <class T>
1253  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1254  {
1255  int32_t start=0;
1256  int32_t end=size-1;
1257 
1258  if (size<1)
1259  return 0;
1260 
1261  while (start<end)
1262  {
1263  int32_t middle=(start+end)/2;
1264 
1265  if (output[middle]>elem)
1266  end=middle-1;
1267  else if (output[middle]<elem)
1268  start=middle+1;
1269  else
1270  return middle;
1271  }
1272 
1273  return start;
1274  }
1275 
1277  static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
1278  {
1279  SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
1280  return int32_t(0);
1281  }
1282 
1283  /* finds an element in a sorted array via binary search
1284  * * returns -1 if not found */
1285  template <class T>
1286  static inline int32_t binary_search(T* output, int32_t size, T elem)
1287  {
1288  int32_t ind = binary_search_helper(output, size, elem);
1289  if (ind >= 0 && output[ind] == elem)
1290  return ind;
1291  return -1;
1292  }
1293 
1295  static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
1296  {
1297  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1298  return int32_t(-1);
1299  }
1300 
1301 
1302  /* Finds an element in a sorted array of pointers via binary search
1303  * Every element is dereferenced once before being compared
1304  *
1305  * @param array array of pointers to search in (assumed being sorted)
1306  * @param length length of array
1307  * @param elem pointer to element to search for
1308  * @return index of elem, -1 if not found */
1309  template<class T>
1310  static inline int32_t binary_search(T** vector, index_t length,
1311  T* elem)
1312  {
1313  int32_t start=0;
1314  int32_t end=length-1;
1315 
1316  if (length<1)
1317  return -1;
1318 
1319  while (start<end)
1320  {
1321  int32_t middle=(start+end)/2;
1322 
1323  if (*vector[middle]>*elem)
1324  end=middle-1;
1325  else if (*vector[middle]<*elem)
1326  start=middle+1;
1327  else
1328  {
1329  start=middle;
1330  break;
1331  }
1332  }
1333 
1334  if (start>=0&&*vector[start]==*elem)
1335  return start;
1336 
1337  return -1;
1338  }
1339 
1341  static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
1342  {
1343  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1344  return int32_t(-1);
1345  }
1346 
1347  template <class T>
1349  T* output, int32_t size, T elem)
1350  {
1351  int32_t ind = binary_search_helper(output, size, elem);
1352 
1353  if (output[ind]<=elem)
1354  return ind;
1355  if (ind>0 && output[ind-1] <= elem)
1356  return ind-1;
1357  return -1;
1358  }
1359 
1362  int32_t size, complex128_t elem)
1363  {
1364  SG_SERROR("CMath::binary_search_max_lower_equal():: \
1365  Not supported for complex128_t\n");
1366  return int32_t(-1);
1367  }
1368 
1371  static float64_t Align(
1372  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1373 
1374 
1376 
1379  {
1380  return c.real();
1381  }
1382 
1385  {
1386  return c.imag();
1387  }
1388 
1390  inline static uint32_t get_seed()
1391  {
1392  return CMath::seed;
1393  }
1394 
1396  inline static uint32_t get_log_range()
1397  {
1398  return CMath::LOGRANGE;
1399  }
1400 
1401 #ifdef USE_LOGCACHE
1402 
1403  inline static uint32_t get_log_accuracy()
1404  {
1405  return CMath::LOGACCURACY;
1406  }
1407 #endif
1408 
1410  static int is_finite(double f);
1411 
1413  static int is_infinity(double f);
1414 
1416  static int is_nan(double f);
1417 
1428 #ifdef USE_LOGCACHE
1429  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1430  {
1431  float64_t diff;
1432 
1433  if (!CMath::is_finite(p))
1434  return q;
1435 
1436  if (!CMath::is_finite(q))
1437  {
1438  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
1439  return NOT_A_NUMBER;
1440  }
1441  diff = p - q;
1442  if (diff > 0)
1443  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1444  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1445  }
1446 
1448  static void init_log_table();
1449 
1451  static int32_t determine_logrange();
1452 
1454  static int32_t determine_logaccuracy(int32_t range);
1455 #else
1457  float64_t p, float64_t q)
1458  {
1459  float64_t diff;
1460 
1461  if (!CMath::is_finite(p))
1462  return q;
1463  if (!CMath::is_finite(q))
1464  return p;
1465  diff = p - q;
1466  if (diff > 0)
1467  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
1468  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
1469  }
1470 #endif
1471 #ifdef USE_LOGSUMARRAY
1472 
1477  static inline float64_t logarithmic_sum_array(
1478  float64_t *p, int32_t len)
1479  {
1480  if (len<=2)
1481  {
1482  if (len==2)
1483  return logarithmic_sum(p[0],p[1]) ;
1484  if (len==1)
1485  return p[0];
1486  return -INFTY ;
1487  }
1488  else
1489  {
1490  float64_t *pp=p ;
1491  if (len%2==1) pp++ ;
1492  for (int32_t j=0; j < len>>1; j++)
1493  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
1494  }
1495  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1496  }
1497 #endif //USE_LOGSUMARRAY
1498 
1499 
1501  virtual const char* get_name() const { return "Math"; }
1502  public:
1505 
1506  static const float64_t NOT_A_NUMBER;
1508  static const float64_t INFTY;
1509  static const float64_t ALMOST_INFTY;
1510 
1513 
1515  static const float64_t PI;
1516 
1519 
1520  /* largest and smallest possible float64_t */
1523 
1524  /* Floating point Limits, Normalized */
1525  static const float32_t F_MAX_VAL32;
1527  static const float64_t F_MAX_VAL64;
1529 
1530  /* Floating point limits, Denormalized */
1531  static const float32_t F_MIN_VAL32;
1532  static const float64_t F_MIN_VAL64;
1533 
1534  protected:
1536  static int32_t LOGRANGE;
1537 
1539  static uint32_t seed;
1540 
1541 #ifdef USE_LOGCACHE
1542 
1544  static int32_t LOGACCURACY;
1546 
1547  static float64_t* logtable;
1548 #endif
1549 };
1550 
1551 //implementations of template functions
1552 template <class T1,class T2>
1554  {
1555  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1556  T1* output=ps->output;
1557  T2* index=ps->index;
1558  int32_t size=ps->size;
1559  int32_t* qsort_threads=ps->qsort_threads;
1560  int32_t sort_limit=ps->sort_limit;
1561  int32_t num_threads=ps->num_threads;
1562 
1563  if (size<2)
1564  {
1565  return NULL;
1566  }
1567 
1568  if (size==2)
1569  {
1570  if (output[0] > output [1])
1571  {
1572  swap(output[0], output[1]);
1573  swap(index[0], index[1]);
1574  }
1575  return NULL;
1576  }
1577  //T1 split=output[random(0,size-1)];
1578  T1 split=output[size/2];
1579 
1580  int32_t left=0;
1581  int32_t right=size-1;
1582 
1583  while (left<=right)
1584  {
1585  while (output[left] < split)
1586  left++;
1587  while (output[right] > split)
1588  right--;
1589 
1590  if (left<=right)
1591  {
1592  swap(output[left], output[right]);
1593  swap(index[left], index[right]);
1594  left++;
1595  right--;
1596  }
1597  }
1598  bool lthread_start=false;
1599  bool rthread_start=false;
1600  pthread_t lthread;
1601  pthread_t rthread;
1602  struct thread_qsort<T1,T2> t1;
1603  struct thread_qsort<T1,T2> t2;
1604 
1605  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1606  qsort_index(output,index,right+1);
1607  else if (right+1> 1)
1608  {
1609  (*qsort_threads)++;
1610  lthread_start=true;
1611  t1.output=output;
1612  t1.index=index;
1613  t1.size=right+1;
1614  t1.qsort_threads=qsort_threads;
1615  t1.sort_limit=sort_limit;
1616  t1.num_threads=num_threads;
1617  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1618  {
1619  lthread_start=false;
1620  (*qsort_threads)--;
1621  qsort_index(output,index,right+1);
1622  }
1623  }
1624 
1625 
1626  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1627  qsort_index(&output[left],&index[left], size-left);
1628  else if (size-left> 1)
1629  {
1630  (*qsort_threads)++;
1631  rthread_start=true;
1632  t2.output=&output[left];
1633  t2.index=&index[left];
1634  t2.size=size-left;
1635  t2.qsort_threads=qsort_threads;
1636  t2.sort_limit=sort_limit;
1637  t2.num_threads=num_threads;
1638  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1639  {
1640  rthread_start=false;
1641  (*qsort_threads)--;
1642  qsort_index(&output[left],&index[left], size-left);
1643  }
1644  }
1645 
1646  if (lthread_start)
1647  {
1648  pthread_join(lthread, NULL);
1649  (*qsort_threads)--;
1650  }
1651 
1652  if (rthread_start)
1653  {
1654  pthread_join(rthread, NULL);
1655  (*qsort_threads)--;
1656  }
1657 
1658  return NULL;
1659  }
1660 
1661  template <class T1,class T2>
1662 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
1663 {
1664  if (size<=1)
1665  return;
1666 
1667  if (size==2)
1668  {
1669  if (output[0] > output [1])
1670  {
1671  swap(output[0],output[1]);
1672  swap(index[0],index[1]);
1673  }
1674  return;
1675  }
1676  //T1 split=output[random(0,size-1)];
1677  T1 split=output[size/2];
1678 
1679  int32_t left=0;
1680  int32_t right=size-1;
1681 
1682  while (left<=right)
1683  {
1684  while (output[left] < split)
1685  left++;
1686  while (output[right] > split)
1687  right--;
1688 
1689  if (left<=right)
1690  {
1691  swap(output[left],output[right]);
1692  swap(index[left],index[right]);
1693  left++;
1694  right--;
1695  }
1696  }
1697 
1698  if (right+1> 1)
1699  qsort_index(output,index,right+1);
1700 
1701  if (size-left> 1)
1702  qsort_index(&output[left],&index[left], size-left);
1703 }
1704 
1705  template <class T1,class T2>
1706 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
1707 {
1708  if (size<=1)
1709  return;
1710 
1711  if (size==2)
1712  {
1713  if (output[0] < output [1])
1714  {
1715  swap(output[0],output[1]);
1716  swap(index[0],index[1]);
1717  }
1718  return;
1719  }
1720 
1721  //T1 split=output[random(0,size-1)];
1722  T1 split=output[size/2];
1723 
1724  int32_t left=0;
1725  int32_t right=size-1;
1726 
1727  while (left<=right)
1728  {
1729  while (output[left] > split)
1730  left++;
1731  while (output[right] < split)
1732  right--;
1733 
1734  if (left<=right)
1735  {
1736  swap(output[left],output[right]);
1737  swap(index[left],index[right]);
1738  left++;
1739  right--;
1740  }
1741  }
1742 
1743  if (right+1> 1)
1744  qsort_backward_index(output,index,right+1);
1745 
1746  if (size-left> 1)
1747  qsort_backward_index(&output[left],&index[left], size-left);
1748 }
1749 
1750  template <class T>
1751 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
1752 {
1753  if (6*n*size<13*size*CMath::log(size))
1754  for (int32_t i=0; i<n; i++)
1755  min(&output[i], &index[i], size-i);
1756  else
1757  qsort_index(output, index, size);
1758 }
1759 
1760 /* move the smallest entry in the array to the beginning */
1761  template <class T>
1762 void CMath::min(float64_t* output, T* index, int32_t size)
1763 {
1764  if (size<=1)
1765  return;
1766  float64_t min_elem=output[0];
1767  int32_t min_index=0;
1768  for (int32_t i=1; i<size; i++)
1769  {
1770  if (output[i]<min_elem)
1771  {
1772  min_index=i;
1773  min_elem=output[i];
1774  }
1775  }
1776  swap(output[0], output[min_index]);
1777  swap(index[0], index[min_index]);
1778 }
1779 
1780 #define COMPLEX128_ERROR_ONEARG_T(function) \
1781 template <> \
1782 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
1783 { \
1784  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1785  #function);\
1786  return complex128_t(0.0, 0.0); \
1787 }
1788 
1789 #define COMPLEX128_ERROR_TWOARGS_T(function) \
1790 template <> \
1791 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \
1792 { \
1793  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1794  #function);\
1795  return complex128_t(0.0, 0.0); \
1796 }
1797 
1798 #define COMPLEX128_ERROR_THREEARGS_T(function) \
1799 template <> \
1800 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \
1801 { \
1802  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1803  #function);\
1804  return complex128_t(0.0, 0.0); \
1805 }
1806 
1807 #define COMPLEX128_ERROR_SORT_T(function) \
1808 template <> \
1809 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \
1810 { \
1811  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1812  #function);\
1813 }
1814 
1815 #define COMPLEX128_ERROR_ARG_MAX_MIN(function) \
1816 template <> \
1817 inline int32_t CMath::function<complex128_t>(complex128_t * a, int32_t b, int32_t c, complex128_t * d) \
1818 { \
1819  int32_t maxIdx=0; \
1820  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1821  #function);\
1822  return maxIdx; \
1823 }
1824 
1827 
1828 
1830 
1833 
1835 // COMPLEX128_ERROR_ONEARG_T(sign)
1836 
1839 
1842 
1845 
1848 
1851 
1852 }
1853 #undef COMPLEX128_ERROR_ONEARG
1854 #undef COMPLEX128_ERROR_ONEARG_T
1855 #undef COMPLEX128_ERROR_TWOARGS_T
1856 #undef COMPLEX128_ERROR_THREEARGS_T
1857 #undef COMPLEX128_STDMATH
1858 #undef COMPLEX128_ERROR_SORT_T
1859 #endif

SHOGUN Machine Learning Toolbox - Documentation