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 
31 #ifndef _USE_MATH_DEFINES
32 #define _USE_MATH_DEFINES
33 #endif
34 
35 #include <math.h>
36 #include <stdio.h>
37 #include <float.h>
38 #include <sys/types.h>
39 #ifndef _WIN32
40 #include <unistd.h>
41 #endif
42 
43 #ifdef HAVE_PTHREAD
44 #include <pthread.h>
45 #endif
46 
47 #ifdef SUNOS
48 #include <ieeefp.h>
49 #endif
50 
52 #ifdef log2
53 #define cygwin_log2 log2
54 #undef log2
55 #endif
56 
57 #ifndef M_PI
58 #define M_PI 3.14159265358979323846
59 #endif
60 
61 #ifdef _WIN32
62 #ifndef isnan
63 #define isnan _isnan
64 #endif
65 
66 #ifndef isfinite
67 #define isfinite _isfinite
68 #endif
69 #endif //_WIN32
70 
71 /* Size of RNG seed */
72 #define RNG_SEED_SIZE 256
73 
74 /* Maximum stack size */
75 #define RADIX_STACK_SIZE 512
76 
77 /* Stack macros */
78 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
79 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
80 
81 #ifndef DOXYGEN_SHOULD_SKIP_THIS
82 
83 template <class T> struct radix_stack_t
84 {
86  T *sa;
88  size_t sn;
90  uint16_t si;
91 };
92 
94 
96 template <class T1, class T2> struct thread_qsort
97 {
99  T1* output;
101  T2* index;
103  uint32_t size;
104 
106  int32_t* qsort_threads;
108  int32_t sort_limit;
110  int32_t num_threads;
111 };
112 #endif // DOXYGEN_SHOULD_SKIP_THIS
113 
114 #define COMPLEX128_ERROR_ONEARG(function) \
115 static inline complex128_t function(complex128_t a) \
116 { \
117  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
118  #function);\
119  return complex128_t(0.0, 0.0); \
120 }
121 
122 #define COMPLEX128_STDMATH(function) \
123 static inline complex128_t function(complex128_t a) \
124 { \
125  return std::function(a); \
126 }
127 
128 namespace shogun
129 {
131  extern CRandom* sg_rand;
132  class CSGObject;
135 class CMath : public CSGObject
136 {
137  public:
141 
142  CMath();
143 
145  virtual ~CMath();
147 
151 
153  //
154  template <class T>
155  static inline T min(T a, T b)
156  {
157  return (a<=b) ? a : b;
158  }
159 
161  template <class T>
162  static inline T max(T a, T b)
163  {
164  return (a>=b) ? a : b;
165  }
166 
168  template <class T>
169  static inline T clamp(T value, T lb, T ub)
170  {
171  if (value<=lb)
172  return lb;
173  else if (value>=ub)
174  return ub;
175  else
176  return value;
177  }
178 
180  template <class T>
181  static inline T abs(T a)
182  {
183  // can't be a>=0?(a):(-a), because compiler complains about
184  // 'comparison always true' when T is unsigned
185  if (a==0)
186  return 0;
187  else if (a>0)
188  return a;
189  else
190  return -a;
191  }
192 
194  static inline float64_t abs(complex128_t a)
195  {
196  float64_t a_real=a.real();
197  float64_t a_imag=a.imag();
198  return (CMath::sqrt(a_real*a_real+a_imag*a_imag));
199  }
201 
204 
211  template <class T>
212  static inline bool fequals_abs(const T& a, const T& b,
213  const float64_t eps)
214  {
215  const T diff = CMath::abs<T>((a-b));
216  return (diff < eps);
217  }
218 
228  template <class T>
229  static inline bool fequals(const T& a, const T& b,
230  const float64_t eps, bool tolerant=false)
231  {
232  const T absA = CMath::abs<T>(a);
233  const T absB = CMath::abs<T>(b);
234  const T diff = CMath::abs<T>((a-b));
235  T comp;
236 
237  // Handle this separately since NAN is unordered
239  return true;
240 
241  // Required for JSON Serialization Tests
242  if (tolerant)
243  return CMath::fequals_abs<T>(a, b, eps);
244 
245  // handles float32_t and float64_t separately
246  if (sizeof(T) == 4)
248 
249  else
251 
252  if (a==b)
253  return true;
254 
255  // both a and b are 0 and relative error is less meaningful
256  else if ( (a==0) || (b==0) || (diff < comp) )
257  return (diff<(eps * comp));
258 
259  // use max(relative error, diff) to handle large eps
260  else
261  {
262  T check = ((diff/(absA + absB)) > diff)?
263  (diff/(absA + absB)):diff;
264  return (check < eps);
265  }
266  }
267 
268  /* Get the corresponding absolute tolerance for unit test given a relative tolerance
269  *
270  * Note that a unit test will be passed only when
271  * \f[
272  * |v_\text{true} - v_\text{predict}| \leq tol_\text{relative} * |v_\text{true}|
273  * \f] which is equivalent to
274  * \f[
275  * |v_\text{true} - v_\text{predict}| \leq tol_\text{absolute}
276  * \f] where
277  * \f[
278  * tol_\text{absolute} = tol_\text{relative} * |v_\text{true}|
279  * \f]
280  *
281  * @param true_value true value should be finite (neither NAN nor INF)
282  * @param rel_tolorance relative tolerance should be positive and less than 1.0
283  *
284  * @return the corresponding absolute tolerance
285  */
286  static float64_t get_abs_tolorance(float64_t true_value, float64_t rel_tolorance);
287 
288  static inline float64_t round(float64_t d)
289  {
290  return ::floor(d+0.5);
291  }
292 
293  static inline float64_t floor(float64_t d)
294  {
295  return ::floor(d);
296  }
297 
298  static inline float64_t ceil(float64_t d)
299  {
300  return ::ceil(d);
301  }
302 
304  template <class T>
305  static inline T sign(T a)
306  {
307  if (a==0)
308  return 0;
309  else return (a<0) ? (-1) : (+1);
310  }
311 
313  template <class T>
314  static inline void swap(T &a,T &b)
315  {
316  /* register for fast swaps */
317  register T c=a;
318  a=b;
319  b=c;
320  }
321 
323  template <class T>
324  static inline T sq(T x)
325  {
326  return x*x;
327  }
328 
330  static inline float32_t sqrt(float32_t x)
331  {
332  return ::sqrtf(x);
333  }
334 
336  static inline float64_t sqrt(float64_t x)
337  {
338  return ::sqrt(x);
339  }
340 
342  static inline floatmax_t sqrt(floatmax_t x)
343  {
344  //fall back to double precision sqrt if sqrtl is not
345  //available
346 #ifdef HAVE_SQRTL
347  return ::sqrtl(x);
348 #else
349  return ::sqrt(x);
350 #endif
351  }
352 
355 
356 
357  static inline float32_t invsqrt(float32_t x)
358  {
359  union float_to_int
360  {
361  float32_t f;
362  int32_t i;
363  };
364 
365  float_to_int tmp;
366  tmp.f=x;
367 
368  float32_t xhalf = 0.5f * x;
369  // store floating-point bits in integer tmp.i
370  // and do initial guess for Newton's method
371  tmp.i = 0x5f3759d5 - (tmp.i >> 1);
372  x = tmp.f; // convert new bits into float
373  x = x*(1.5f - xhalf*x*x); // One round of Newton's method
374  return x;
375  }
376 
378  static inline floatmax_t powl(floatmax_t x, floatmax_t n)
379  {
380  //fall back to double precision pow if powl is not
381  //available
382 #ifdef HAVE_POWL
383  return ::powl((long double) x, (long double) n);
384 #else
385  return ::pow((double) x, (double) n);
386 #endif
387  }
388 
389  static inline int32_t pow(bool x, int32_t n)
390  {
391  return 0;
392  }
393 
394  static inline int32_t pow(int32_t x, int32_t n)
395  {
396  ASSERT(n>=0)
397  int32_t result=1;
398  while (n--)
399  result*=x;
400 
401  return result;
402  }
403 
404  static inline float64_t pow(float64_t x, int32_t n)
405  {
406  if (n>=0)
407  {
408  float64_t result=1;
409  while (n--)
410  result*=x;
411 
412  return result;
413  }
414  else
415  return ::pow((double)x, (double)n);
416  }
417 
418  static inline float64_t pow(float64_t x, float64_t n)
419  {
420  return ::pow((double) x, (double) n);
421  }
422 
424  static inline complex128_t pow(complex128_t x, int32_t n)
425  {
426  return std::pow(x, n);
427  }
428 
430  {
431  return std::pow(x, n);
432  }
433 
435  {
436  return std::pow(x, n);
437  }
438 
440  {
441  return std::pow(x, n);
442  }
443 
444  static inline float64_t exp(float64_t x)
445  {
446  return ::exp((double) x);
447  }
448 
451 
452 
453  static inline float64_t tan(float64_t x)
454  {
455  return ::tan((double) x);
456  }
457 
460 
461 
462  static inline float64_t atan(float64_t x)
463  {
464  return ::atan((double) x);
465  }
466 
469 
470 
471  static inline float64_t atan2(float64_t x, float64_t y)
472  {
473  return ::atan2((double) x, (double) y);
474  }
475 
478 
479 
480  static inline float64_t tanh(float64_t x)
481  {
482  return ::tanh((double) x);
483  }
484 
487 
488  static inline float64_t log10(float64_t v)
489  {
490  return ::log(v)/::log(10.0);
491  }
492 
495 
496  static inline float64_t log2(float64_t v)
497  {
498 #ifdef HAVE_LOG2
499  return ::log2(v);
500 #else
501  return ::log(v)/::log(2.0);
502 #endif //HAVE_LOG2
503  }
504 
505  static inline float64_t log(float64_t v)
506  {
507  return ::log(v);
508  }
509 
512 
513  static inline index_t floor_log(index_t n)
514  {
515  index_t i;
516  for (i = 0; n != 0; i++)
517  n >>= 1;
518 
519  return i;
520  }
521 
522  static inline float64_t sin(float64_t x)
523  {
524  return ::sin(x);
525  }
526 
529 
530  static inline float64_t asin(float64_t x)
531  {
532  return ::asin(x);
533  }
534 
537 
538  static inline float64_t sinh(float64_t x)
539  {
540  return ::sinh(x);
541  }
542 
545 
546  static inline float64_t cos(float64_t x)
547  {
548  return ::cos(x);
549  }
550 
553 
554  static inline float64_t acos(float64_t x)
555  {
556  return ::acos(x);
557  }
558 
561 
562  static inline float64_t cosh(float64_t x)
563  {
564  return ::cosh(x);
565  }
566 
569 
570  static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
571  {
572  ASSERT(len>0 && xy)
573 
574  float64_t area = 0.0;
575 
576  if (!reversed)
577  {
578  for (int i=1; i<len; i++)
579  area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
580  }
581  else
582  {
583  for (int i=1; i<len; i++)
584  area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
585  }
586 
587  return area;
588  }
589 
590  static bool strtof(const char* str, float32_t* float_result);
591  static bool strtod(const char* str, float64_t* double_result);
592  static bool strtold(const char* str, floatmax_t* long_double_result);
593 
594  static inline int64_t factorial(int32_t n)
595  {
596  int64_t res=1;
597  for (int i=2; i<=n; i++)
598  res*=i ;
599  return res ;
600  }
601 
602  static void init_random(uint32_t initseed=0)
603  {
604  if (initseed==0)
606  else
607  seed=initseed;
608 
610  }
611 
612  static inline uint64_t random()
613  {
614  return sg_rand->random_64();
615  }
616 
617  static inline uint64_t random(uint64_t min_value, uint64_t max_value)
618  {
619  return sg_rand->random(min_value, max_value);
620  }
621 
622  static inline int64_t random(int64_t min_value, int64_t max_value)
623  {
624  return sg_rand->random(min_value, max_value);
625  }
626 
627  static inline uint32_t random(uint32_t min_value, uint32_t max_value)
628  {
629  return sg_rand->random(min_value, max_value);
630  }
631 
632  static inline int32_t random(int32_t min_value, int32_t max_value)
633  {
634  return sg_rand->random(min_value, max_value);
635  }
636 
637  static inline float32_t random(float32_t min_value, float32_t max_value)
638  {
639  return sg_rand->random(min_value, max_value);
640  }
641 
642  static inline float64_t random(float64_t min_value, float64_t max_value)
643  {
644  return sg_rand->random(min_value, max_value);
645  }
646 
647  static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
648  {
649  return sg_rand->random(min_value, max_value);
650  }
651 
655  static inline float32_t normal_random(float32_t mean, float32_t std_dev)
656  {
657  // sets up variables & makes sure rand_s.range == (0,1)
658  float32_t ret;
659  float32_t rand_u;
660  float32_t rand_v;
661  float32_t rand_s;
662  do
663  {
664  rand_u = CMath::random(-1.0, 1.0);
665  rand_v = CMath::random(-1.0, 1.0);
666  rand_s = rand_u*rand_u + rand_v*rand_v;
667  } while ((rand_s == 0) || (rand_s >= 1));
668 
669  // the meat & potatos, and then the mean & standard deviation shifting...
670  ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
671  ret = std_dev*ret + mean;
672  return ret;
673  }
674 
678  static inline float64_t normal_random(float64_t mean, float64_t std_dev)
679  {
680  return sg_rand->normal_distrib(mean, std_dev);
681  }
682 
685  static inline float32_t randn_float()
686  {
687  return normal_random(0.0, 1.0);
688  }
689 
692  static inline float64_t randn_double()
693  {
694  return sg_rand->std_normal_distrib();
695  }
696 
697  template <class T>
698  static int32_t get_num_nonzero(T* vec, int32_t len)
699  {
700  int32_t nnz = 0;
701  for (index_t i=0; i<len; ++i)
702  nnz += vec[i] != 0;
703 
704  return nnz;
705  }
706 
707  static int32_t get_num_nonzero(complex128_t* vec, int32_t len)
708  {
709  int32_t nnz=0;
710  for (index_t i=0; i<len; ++i)
711  nnz+=vec[i]!=0.0;
712  return nnz;
713  }
714 
715  static inline int64_t nchoosek(int32_t n, int32_t k)
716  {
717  int64_t res=1;
718 
719  for (int32_t i=n-k+1; i<=n; i++)
720  res*=i;
721 
722  return res/factorial(k);
723  }
724 
732  static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100);
733 
740  template <class T>
741  static T log_sum_exp(SGVector<T> values)
742  {
743  REQUIRE(values.vector, "Values are empty");
744 
745  /* find minimum element index */
746  index_t min_index=0;
747  T X0=values[0];
748  if (values.vlen>1)
749  {
750  for (index_t i=1; i<values.vlen; ++i)
751  {
752  if (values[i]<X0)
753  {
754  X0=values[i];
755  min_index=i;
756  }
757  }
758  }
759 
760  /* remove element from vector copy and compute log sum exp */
761  SGVector<T> values_without_X0(values.vlen-1);
762  index_t from_idx=0;
763  index_t to_idx=0;
764  for (from_idx=0; from_idx<values.vlen; ++from_idx)
765  {
766  if (from_idx!=min_index)
767  {
768  values_without_X0[to_idx]=exp(values[from_idx]-X0);
769  to_idx++;
770  }
771  }
772 
773  return X0+log(SGVector<T>::sum(values_without_X0)+1);
774  }
775 
782  template <class T>
783  static T log_mean_exp(SGVector<T> values)
784  {
785  return log_sum_exp(values) - log(values.vlen);
786  }
787 
791  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
792  static void sort(float64_t *a, int32_t*idx, int32_t N);
793 
796  template <class T>
797  static void qsort(T* output, int32_t size)
798  {
799  if (size<=1)
800  return;
801 
802  if (size==2)
803  {
804  if (output[0] > output [1])
805  CMath::swap(output[0],output[1]);
806  return;
807  }
808  //T split=output[random(0,size-1)];
809  T split=output[size/2];
810 
811  int32_t left=0;
812  int32_t right=size-1;
813 
814  while (left<=right)
815  {
816  while (output[left] < split)
817  left++;
818  while (output[right] > split)
819  right--;
820 
821  if (left<=right)
822  {
823  CMath::swap(output[left],output[right]);
824  left++;
825  right--;
826  }
827  }
828 
829  if (right+1> 1)
830  qsort(output,right+1);
831 
832  if (size-left> 1)
833  qsort(&output[left],size-left);
834  }
835 
838  template <class T>
839  static void insertion_sort(T* output, int32_t size)
840  {
841  for (int32_t i=0; i<size-1; i++)
842  {
843  int32_t j=i-1;
844  T value=output[i];
845  while (j >= 0 && output[j] > value)
846  {
847  output[j+1] = output[j];
848  j--;
849  }
850  output[j+1]=value;
851  }
852  }
853 
855  template <class T>
856  inline static void radix_sort(T* array, int32_t size)
857  {
858  radix_sort_helper(array,size,0);
859  }
860 
861  /*
862  * Inline function to extract the byte at position p (from left)
863  * of an 64 bit integer. The function is somewhat identical to
864  * accessing an array of characters via [].
865  */
866  template <class T>
867  static inline uint8_t byte(T word, uint16_t p)
868  {
869  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
870  }
871 
873  static inline uint8_t byte(complex128_t word, uint16_t p)
874  {
875  SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
876  return uint8_t(0);
877  }
878 
879  template <class T>
880  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
881  {
882  static size_t count[256], nc, cmin;
883  T *ak;
884  uint8_t c=0;
885  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
886  T *an, *aj, *pile[256];
887  size_t *cp, cmax;
888 
889  /* Push initial array to stack */
890  sp = s;
891  radix_push(array, size, i);
892 
893  /* Loop until all digits have been sorted */
894  while (sp>s) {
895  radix_pop(array, size, i);
896  an = array + size;
897 
898  /* Make character histogram */
899  if (nc == 0) {
900  cmin = 0xff;
901  for (ak = array; ak < an; ak++) {
902  c = byte(*ak, i);
903  count[c]++;
904  if (count[c] == 1) {
905  /* Determine smallest character */
906  if (c < cmin)
907  cmin = c;
908  nc++;
909  }
910  }
911 
912  /* Sort recursively if stack size too small */
913  if (sp + nc > s + RADIX_STACK_SIZE) {
914  radix_sort_helper(array, size, i);
915  continue;
916  }
917  }
918 
919  /*
920  * Set pile[]; push incompletely sorted bins onto stack.
921  * pile[] = pointers to last out-of-place element in bins.
922  * Before permuting: pile[c-1] + count[c] = pile[c];
923  * during deal: pile[c] counts down to pile[c-1].
924  */
925  olds = bigs = sp;
926  cmax = 2;
927  ak = array;
928  pile[0xff] = an;
929  for (cp = count + cmin; nc > 0; cp++) {
930  /* Find next non-empty pile */
931  while (*cp == 0)
932  cp++;
933  /* Pile with several entries */
934  if (*cp > 1) {
935  /* Determine biggest pile */
936  if (*cp > cmax) {
937  cmax = *cp;
938  bigs = sp;
939  }
940 
941  if (i < sizeof(T)-1)
942  radix_push(ak, *cp, (uint16_t) (i + 1));
943  }
944  pile[cp - count] = ak += *cp;
945  nc--;
946  }
947 
948  /* Play it safe -- biggest bin last. */
949  swap(*olds, *bigs);
950 
951  /*
952  * Permute misplacements home. Already home: everything
953  * before aj, and in pile[c], items from pile[c] on.
954  * Inner loop:
955  * r = next element to put in place;
956  * ak = pile[r[i]] = location to put the next element.
957  * aj = bottom of 1st disordered bin.
958  * Outer loop:
959  * Once the 1st disordered bin is done, ie. aj >= ak,
960  * aj<-aj + count[c] connects the bins in array linked list;
961  * reset count[c].
962  */
963  aj = array;
964  while (aj<an)
965  {
966  T r;
967 
968  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
969  swap(*ak, r);
970 
971  *aj = r;
972  aj += count[c];
973  count[c] = 0;
974  }
975  }
976  }
977 
979  static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
980  {
981  SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
982  }
983 
993  template <class T>
994  static void qsort(T** vector, index_t length)
995  {
996  if (length<=1)
997  return;
998 
999  if (length==2)
1000  {
1001  if (*vector[0]>*vector[1])
1002  swap(vector[0],vector[1]);
1003  return;
1004  }
1005  T* split=vector[length/2];
1006 
1007  int32_t left=0;
1008  int32_t right=length-1;
1009 
1010  while (left<=right)
1011  {
1012  while (*vector[left]<*split)
1013  ++left;
1014  while (*vector[right]>*split)
1015  --right;
1016 
1017  if (left<=right)
1018  {
1019  swap(vector[left],vector[right]);
1020  ++left;
1021  --right;
1022  }
1023  }
1024 
1025  if (right+1>1)
1026  qsort(vector,right+1);
1027 
1028  if (length-left>1)
1029  qsort(&vector[left],length-left);
1030  }
1031 
1033  static void qsort(complex128_t** vector, index_t length)
1034  {
1035  SG_SERROR("CMath::qsort():: Not supported for complex128_t\n");
1036  }
1037 
1039  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1040  {
1041  ASSERT(width>=0)
1042  for (int i=0; i<width; i++)
1043  {
1044  T mask = ((T) 1)<<(sizeof(T)*8-1);
1045  while (mask)
1046  {
1047  if (mask & word)
1048  SG_SPRINT("1")
1049  else
1050  SG_SPRINT("0")
1051 
1052  mask>>=1;
1053  }
1054  }
1055  }
1056 
1058  static void display_bits(complex128_t word,
1059  int32_t width=8*sizeof(complex128_t))
1060  {
1061  SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
1062  }
1063 
1069  template <class T1,class T2>
1070  static void qsort_index(T1* output, T2* index, uint32_t size);
1071 
1073  template <class T>
1074  static void qsort_index(complex128_t* output, T* index, uint32_t size)
1075  {
1076  SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
1077  }
1078 
1084  template <class T1,class T2>
1085  static void qsort_backward_index(
1086  T1* output, T2* index, int32_t size);
1087 
1089  template <class T>
1091  complex128_t* output, T* index, uint32_t size)
1092  {
1093  SG_SERROR("CMath::qsort_backword_index():: \
1094  Not supported for complex128_t\n");
1095  }
1096 
1104  template <class T1,class T2>
1105  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1106  {
1107  int32_t n=0;
1108  thread_qsort<T1,T2> t;
1109  t.output=output;
1110  t.index=index;
1111  t.size=size;
1112  t.qsort_threads=&n;
1113  t.sort_limit=limit;
1114  t.num_threads=n_threads;
1115  parallel_qsort_index<T1,T2>(&t);
1116  }
1117 
1119  template <class T>
1120  inline static void parallel_qsort_index(complex128_t* output, T* index,
1121  uint32_t size, int32_t n_threads, int32_t limit=0)
1122  {
1123  SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1124  }
1125 
1126 
1127  template <class T1,class T2>
1128  static void* parallel_qsort_index(void* p);
1129 
1130 
1131  /* finds the smallest element in output and puts that element as the
1132  first element */
1133  template <class T>
1134  static void min(float64_t* output, T* index, int32_t size);
1135 
1137  static void min(float64_t* output, complex128_t* index, int32_t size)
1138  {
1139  SG_SERROR("CMath::min():: Not supported for complex128_t\n");
1140  }
1141 
1142  /* finds the n smallest elements in output and puts these elements as the
1143  first n elements */
1144  template <class T>
1145  static void nmin(
1146  float64_t* output, T* index, int32_t size, int32_t n);
1147 
1149  static void nmin(float64_t* output, complex128_t* index,
1150  int32_t size, int32_t n)
1151  {
1152  SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
1153  }
1154 
1155 
1156 
1157  /* finds an element in a sorted array via binary search
1158  * returns -1 if not found */
1159  template <class T>
1160  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1161  {
1162  int32_t start=0;
1163  int32_t end=size-1;
1164 
1165  if (size<1)
1166  return 0;
1167 
1168  while (start<end)
1169  {
1170  int32_t middle=(start+end)/2;
1171 
1172  if (output[middle]>elem)
1173  end=middle-1;
1174  else if (output[middle]<elem)
1175  start=middle+1;
1176  else
1177  return middle;
1178  }
1179 
1180  return start;
1181  }
1182 
1184  static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
1185  {
1186  SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
1187  return int32_t(0);
1188  }
1189 
1190  /* finds an element in a sorted array via binary search
1191  * * returns -1 if not found */
1192  template <class T>
1193  static inline int32_t binary_search(T* output, int32_t size, T elem)
1194  {
1195  int32_t ind = binary_search_helper(output, size, elem);
1196  if (ind >= 0 && output[ind] == elem)
1197  return ind;
1198  return -1;
1199  }
1200 
1202  static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
1203  {
1204  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1205  return int32_t(-1);
1206  }
1207 
1208 
1209  /* Finds an element in a sorted array of pointers via binary search
1210  * Every element is dereferenced once before being compared
1211  *
1212  * @param array array of pointers to search in (assumed being sorted)
1213  * @param length length of array
1214  * @param elem pointer to element to search for
1215  * @return index of elem, -1 if not found */
1216  template<class T>
1217  static inline int32_t binary_search(T** vector, index_t length,
1218  T* elem)
1219  {
1220  int32_t start=0;
1221  int32_t end=length-1;
1222 
1223  if (length<1)
1224  return -1;
1225 
1226  while (start<end)
1227  {
1228  int32_t middle=(start+end)/2;
1229 
1230  if (*vector[middle]>*elem)
1231  end=middle-1;
1232  else if (*vector[middle]<*elem)
1233  start=middle+1;
1234  else
1235  {
1236  start=middle;
1237  break;
1238  }
1239  }
1240 
1241  if (start>=0&&*vector[start]==*elem)
1242  return start;
1243 
1244  return -1;
1245  }
1246 
1248  static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
1249  {
1250  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1251  return int32_t(-1);
1252  }
1253 
1254  template <class T>
1256  T* output, int32_t size, T elem)
1257  {
1258  int32_t ind = binary_search_helper(output, size, elem);
1259 
1260  if (output[ind]<=elem)
1261  return ind;
1262  if (ind>0 && output[ind-1] <= elem)
1263  return ind-1;
1264  return -1;
1265  }
1266 
1269  int32_t size, complex128_t elem)
1270  {
1271  SG_SERROR("CMath::binary_search_max_lower_equal():: \
1272  Not supported for complex128_t\n");
1273  return int32_t(-1);
1274  }
1275 
1278  static float64_t Align(
1279  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1280 
1281 
1283 
1286  {
1287  return c.real();
1288  }
1289 
1292  {
1293  return c.imag();
1294  }
1295 
1297  inline static uint32_t get_seed()
1298  {
1299  return CMath::seed;
1300  }
1301 
1303  inline static uint32_t get_log_range()
1304  {
1305  return CMath::LOGRANGE;
1306  }
1307 
1308 #ifdef USE_LOGCACHE
1309 
1310  inline static uint32_t get_log_accuracy()
1311  {
1312  return CMath::LOGACCURACY;
1313  }
1314 #endif
1315 
1317  static int is_finite(double f);
1318 
1320  static int is_infinity(double f);
1321 
1323  static int is_nan(double f);
1324 
1335 #ifdef USE_LOGCACHE
1336  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1337  {
1338  float64_t diff;
1339 
1340  if (!CMath::is_finite(p))
1341  return q;
1342 
1343  if (!CMath::is_finite(q))
1344  {
1345  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
1346  return NOT_A_NUMBER;
1347  }
1348  diff = p - q;
1349  if (diff > 0)
1350  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1351  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1352  }
1353 
1355  static void init_log_table();
1356 
1358  static int32_t determine_logrange();
1359 
1361  static int32_t determine_logaccuracy(int32_t range);
1362 #else
1364  float64_t p, float64_t q)
1365  {
1366  float64_t diff;
1367 
1368  if (!CMath::is_finite(p))
1369  return q;
1370  if (!CMath::is_finite(q))
1371  return p;
1372  diff = p - q;
1373  if (diff > 0)
1374  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
1375  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
1376  }
1377 #endif
1378 #ifdef USE_LOGSUMARRAY
1379 
1384  static inline float64_t logarithmic_sum_array(
1385  float64_t *p, int32_t len)
1386  {
1387  if (len<=2)
1388  {
1389  if (len==2)
1390  return logarithmic_sum(p[0],p[1]) ;
1391  if (len==1)
1392  return p[0];
1393  return -INFTY ;
1394  }
1395  else
1396  {
1397  register float64_t *pp=p ;
1398  if (len%2==1) pp++ ;
1399  for (register int32_t j=0; j < len>>1; j++)
1400  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
1401  }
1402  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1403  }
1404 #endif //USE_LOGSUMARRAY
1405 
1406 
1408  virtual const char* get_name() const { return "Math"; }
1409  public:
1412 
1413  static const float64_t NOT_A_NUMBER;
1415  static const float64_t INFTY;
1416  static const float64_t ALMOST_INFTY;
1417 
1420 
1422  static const float64_t PI;
1423 
1426 
1427  /* largest and smallest possible float64_t */
1430 
1431  /* Floating point Limits, Normalized */
1432  static const float32_t F_MAX_VAL32;
1434  static const float64_t F_MAX_VAL64;
1436 
1437  /* Floating point limits, Denormalized */
1438  static const float32_t F_MIN_VAL32;
1439  static const float64_t F_MIN_VAL64;
1440 
1441  protected:
1443  static int32_t LOGRANGE;
1444 
1446  static uint32_t seed;
1447 
1448 #ifdef USE_LOGCACHE
1449 
1451  static int32_t LOGACCURACY;
1453 
1454  static float64_t* logtable;
1455 #endif
1456 };
1457 
1458 //implementations of template functions
1459 template <class T1,class T2>
1461  {
1462  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1463  T1* output=ps->output;
1464  T2* index=ps->index;
1465  int32_t size=ps->size;
1466  int32_t* qsort_threads=ps->qsort_threads;
1467  int32_t sort_limit=ps->sort_limit;
1468  int32_t num_threads=ps->num_threads;
1469 
1470  if (size<2)
1471  {
1472  return NULL;
1473  }
1474 
1475  if (size==2)
1476  {
1477  if (output[0] > output [1])
1478  {
1479  swap(output[0], output[1]);
1480  swap(index[0], index[1]);
1481  }
1482  return NULL;
1483  }
1484  //T1 split=output[random(0,size-1)];
1485  T1 split=output[size/2];
1486 
1487  int32_t left=0;
1488  int32_t right=size-1;
1489 
1490  while (left<=right)
1491  {
1492  while (output[left] < split)
1493  left++;
1494  while (output[right] > split)
1495  right--;
1496 
1497  if (left<=right)
1498  {
1499  swap(output[left], output[right]);
1500  swap(index[left], index[right]);
1501  left++;
1502  right--;
1503  }
1504  }
1505  bool lthread_start=false;
1506  bool rthread_start=false;
1507  pthread_t lthread;
1508  pthread_t rthread;
1509  struct thread_qsort<T1,T2> t1;
1510  struct thread_qsort<T1,T2> t2;
1511 
1512  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1513  qsort_index(output,index,right+1);
1514  else if (right+1> 1)
1515  {
1516  (*qsort_threads)++;
1517  lthread_start=true;
1518  t1.output=output;
1519  t1.index=index;
1520  t1.size=right+1;
1521  t1.qsort_threads=qsort_threads;
1522  t1.sort_limit=sort_limit;
1523  t1.num_threads=num_threads;
1524  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1525  {
1526  lthread_start=false;
1527  (*qsort_threads)--;
1528  qsort_index(output,index,right+1);
1529  }
1530  }
1531 
1532 
1533  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1534  qsort_index(&output[left],&index[left], size-left);
1535  else if (size-left> 1)
1536  {
1537  (*qsort_threads)++;
1538  rthread_start=true;
1539  t2.output=&output[left];
1540  t2.index=&index[left];
1541  t2.size=size-left;
1542  t2.qsort_threads=qsort_threads;
1543  t2.sort_limit=sort_limit;
1544  t2.num_threads=num_threads;
1545  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1546  {
1547  rthread_start=false;
1548  (*qsort_threads)--;
1549  qsort_index(&output[left],&index[left], size-left);
1550  }
1551  }
1552 
1553  if (lthread_start)
1554  {
1555  pthread_join(lthread, NULL);
1556  (*qsort_threads)--;
1557  }
1558 
1559  if (rthread_start)
1560  {
1561  pthread_join(rthread, NULL);
1562  (*qsort_threads)--;
1563  }
1564 
1565  return NULL;
1566  }
1567 
1568  template <class T1,class T2>
1569 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
1570 {
1571  if (size<=1)
1572  return;
1573 
1574  if (size==2)
1575  {
1576  if (output[0] > output [1])
1577  {
1578  swap(output[0],output[1]);
1579  swap(index[0],index[1]);
1580  }
1581  return;
1582  }
1583  //T1 split=output[random(0,size-1)];
1584  T1 split=output[size/2];
1585 
1586  int32_t left=0;
1587  int32_t right=size-1;
1588 
1589  while (left<=right)
1590  {
1591  while (output[left] < split)
1592  left++;
1593  while (output[right] > split)
1594  right--;
1595 
1596  if (left<=right)
1597  {
1598  swap(output[left],output[right]);
1599  swap(index[left],index[right]);
1600  left++;
1601  right--;
1602  }
1603  }
1604 
1605  if (right+1> 1)
1606  qsort_index(output,index,right+1);
1607 
1608  if (size-left> 1)
1609  qsort_index(&output[left],&index[left], size-left);
1610 }
1611 
1612  template <class T1,class T2>
1613 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
1614 {
1615  if (size<=1)
1616  return;
1617 
1618  if (size==2)
1619  {
1620  if (output[0] < output [1])
1621  {
1622  swap(output[0],output[1]);
1623  swap(index[0],index[1]);
1624  }
1625  return;
1626  }
1627 
1628  //T1 split=output[random(0,size-1)];
1629  T1 split=output[size/2];
1630 
1631  int32_t left=0;
1632  int32_t right=size-1;
1633 
1634  while (left<=right)
1635  {
1636  while (output[left] > split)
1637  left++;
1638  while (output[right] < split)
1639  right--;
1640 
1641  if (left<=right)
1642  {
1643  swap(output[left],output[right]);
1644  swap(index[left],index[right]);
1645  left++;
1646  right--;
1647  }
1648  }
1649 
1650  if (right+1> 1)
1651  qsort_backward_index(output,index,right+1);
1652 
1653  if (size-left> 1)
1654  qsort_backward_index(&output[left],&index[left], size-left);
1655 }
1656 
1657  template <class T>
1658 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
1659 {
1660  if (6*n*size<13*size*CMath::log(size))
1661  for (int32_t i=0; i<n; i++)
1662  min(&output[i], &index[i], size-i) ;
1663  else
1664  qsort_index(output, index, size) ;
1665 }
1666 
1667 /* move the smallest entry in the array to the beginning */
1668  template <class T>
1669 void CMath::min(float64_t* output, T* index, int32_t size)
1670 {
1671  if (size<=1)
1672  return;
1673  float64_t min_elem=output[0];
1674  int32_t min_index=0;
1675  for (int32_t i=1; i<size; i++)
1676  {
1677  if (output[i]<min_elem)
1678  {
1679  min_index=i;
1680  min_elem=output[i];
1681  }
1682  }
1683  swap(output[0], output[min_index]);
1684  swap(index[0], index[min_index]);
1685 }
1686 
1687 #define COMPLEX128_ERROR_ONEARG_T(function) \
1688 template <> \
1689 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
1690 { \
1691  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1692  #function);\
1693  return complex128_t(0.0, 0.0); \
1694 }
1695 
1696 #define COMPLEX128_ERROR_TWOARGS_T(function) \
1697 template <> \
1698 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \
1699 { \
1700  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1701  #function);\
1702  return complex128_t(0.0, 0.0); \
1703 }
1704 
1705 #define COMPLEX128_ERROR_THREEARGS_T(function) \
1706 template <> \
1707 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \
1708 { \
1709  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1710  #function);\
1711  return complex128_t(0.0, 0.0); \
1712 }
1713 
1714 #define COMPLEX128_ERROR_SORT_T(function) \
1715 template <> \
1716 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \
1717 { \
1718  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1719  #function);\
1720 }
1721 
1724 
1725 
1727 
1730 
1732 // COMPLEX128_ERROR_ONEARG_T(sign)
1733 
1736 
1739 
1742 
1743 }
1744 #undef COMPLEX128_ERROR_ONEARG
1745 #undef COMPLEX128_ERROR_ONEARG_T
1746 #undef COMPLEX128_ERROR_TWOARGS_T
1747 #undef COMPLEX128_ERROR_THREEARGS_T
1748 #undef COMPLEX128_STDMATH
1749 #undef COMPLEX128_ERROR_SORT_T
1750 #endif

SHOGUN Machine Learning Toolbox - Documentation