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  template <class T>
170  static inline T clamp(T value, T lb, T ub)
171  {
172  if (value<=lb)
173  return lb;
174  else if (value>=ub)
175  return ub;
176  else
177  return value;
178  }
179 
181  template <class T>
182  static inline T abs(T a)
183  {
184  // can't be a>=0?(a):(-a), because compiler complains about
185  // 'comparison always true' when T is unsigned
186  if (a==0)
187  return 0;
188  else if (a>0)
189  return a;
190  else
191  return -a;
192  }
193 
195  static inline float64_t abs(complex128_t a)
196  {
197  float64_t a_real=a.real();
198  float64_t a_imag=a.imag();
199  return (CMath::sqrt(a_real*a_real+a_imag*a_imag));
200  }
202 
205 
212  template <class T>
213  static inline bool fequals_abs(const T& a, const T& b,
214  const float64_t eps)
215  {
216  const T diff = CMath::abs<T>((a-b));
217  return (diff < eps);
218  }
219 
229  template <class T>
230  static inline bool fequals(const T& a, const T& b,
231  const float64_t eps, bool tolerant=false)
232  {
233  const T absA = CMath::abs<T>(a);
234  const T absB = CMath::abs<T>(b);
235  const T diff = CMath::abs<T>((a-b));
236  T comp;
237 
238  // Handle this separately since NAN is unordered
240  return true;
241 
242  // Required for JSON Serialization Tests
243  if (tolerant)
244  return CMath::fequals_abs<T>(a, b, eps);
245 
246  // handles float32_t and float64_t separately
247  if (sizeof(T) == 4)
249 
250  else
252 
253  if (a==b)
254  return true;
255 
256  // both a and b are 0 and relative error is less meaningful
257  else if ( (a==0) || (b==0) || (diff < comp) )
258  return (diff<(eps * comp));
259 
260  // use max(relative error, diff) to handle large eps
261  else
262  {
263  T check = ((diff/(absA + absB)) > diff)?
264  (diff/(absA + absB)):diff;
265  return (check < eps);
266  }
267  }
268 
269  /* Get the corresponding absolute tolerance for unit test given a relative tolerance
270  *
271  * Note that a unit test will be passed only when
272  * \f[
273  * |v_\text{true} - v_\text{predict}| \leq tol_\text{relative} * |v_\text{true}|
274  * \f] which is equivalent to
275  * \f[
276  * |v_\text{true} - v_\text{predict}| \leq tol_\text{absolute}
277  * \f] where
278  * \f[
279  * tol_\text{absolute} = tol_\text{relative} * |v_\text{true}|
280  * \f]
281  *
282  * @param true_value true value should be finite (neither NAN nor INF)
283  * @param rel_tolorance relative tolerance should be positive and less than 1.0
284  *
285  * @return the corresponding absolute tolerance
286  */
287  static float64_t get_abs_tolorance(float64_t true_value, float64_t rel_tolorance);
288 
289  static inline float64_t round(float64_t d)
290  {
291  return ::floor(d+0.5);
292  }
293 
294  static inline float64_t floor(float64_t d)
295  {
296  return ::floor(d);
297  }
298 
299  static inline float64_t ceil(float64_t d)
300  {
301  return ::ceil(d);
302  }
303 
305  template <class T>
306  static inline T sign(T a)
307  {
308  if (a==0)
309  return 0;
310  else return (a<0) ? (-1) : (+1);
311  }
312 
314  template <class T>
315  static inline void swap(T &a,T &b)
316  {
317  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 
703  template <class T>
704  static void permute(SGVector<T> v, CRandom* rand=NULL)
705  {
706  if (rand)
707  {
708  for (index_t i=0; i<v.vlen; ++i)
709  swap(v[i], v[rand->random(i, v.vlen-1)]);
710  }
711  else
712  {
713  for (index_t i=0; i<v.vlen; ++i)
714  swap(v[i], v[random(i, v.vlen-1)]);
715  }
716  }
717 
718  template <class T>
719  static int32_t get_num_nonzero(T* vec, int32_t len)
720  {
721  int32_t nnz = 0;
722  for (index_t i=0; i<len; ++i)
723  nnz += vec[i] != 0;
724 
725  return nnz;
726  }
727 
728  static int32_t get_num_nonzero(complex128_t* vec, int32_t len)
729  {
730  int32_t nnz=0;
731  for (index_t i=0; i<len; ++i)
732  nnz+=vec[i]!=0.0;
733  return nnz;
734  }
735 
736  static inline int64_t nchoosek(int32_t n, int32_t k)
737  {
738  int64_t res=1;
739 
740  for (int32_t i=n-k+1; i<=n; i++)
741  res*=i;
742 
743  return res/factorial(k);
744  }
745 
753  static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100);
754 
761  template <class T>
762  static T log_sum_exp(SGVector<T> values)
763  {
764  REQUIRE(values.vector, "Values are empty");
765 
766  /* find minimum element index */
767  index_t min_index=0;
768  T X0=values[0];
769  if (values.vlen>1)
770  {
771  for (index_t i=1; i<values.vlen; ++i)
772  {
773  if (values[i]<X0)
774  {
775  X0=values[i];
776  min_index=i;
777  }
778  }
779  }
780 
781  /* remove element from vector copy and compute log sum exp */
782  SGVector<T> values_without_X0(values.vlen-1);
783  index_t from_idx=0;
784  index_t to_idx=0;
785  for (from_idx=0; from_idx<values.vlen; ++from_idx)
786  {
787  if (from_idx!=min_index)
788  {
789  values_without_X0[to_idx]=exp(values[from_idx]-X0);
790  to_idx++;
791  }
792  }
793 
794  return X0+log(SGVector<T>::sum(values_without_X0)+1);
795  }
796 
803  template <class T>
804  static T log_mean_exp(SGVector<T> values)
805  {
806  return log_sum_exp(values) - log(values.vlen);
807  }
808 
812  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
813  static void sort(float64_t *a, int32_t*idx, int32_t N);
814 
817  template <class T>
818  static void qsort(T* output, int32_t size)
819  {
820  if (size<=1)
821  return;
822 
823  if (size==2)
824  {
825  if (output[0] > output [1])
826  CMath::swap(output[0],output[1]);
827  return;
828  }
829  //T split=output[random(0,size-1)];
830  T split=output[size/2];
831 
832  int32_t left=0;
833  int32_t right=size-1;
834 
835  while (left<=right)
836  {
837  while (output[left] < split)
838  left++;
839  while (output[right] > split)
840  right--;
841 
842  if (left<=right)
843  {
844  CMath::swap(output[left],output[right]);
845  left++;
846  right--;
847  }
848  }
849 
850  if (right+1> 1)
851  qsort(output,right+1);
852 
853  if (size-left> 1)
854  qsort(&output[left],size-left);
855  }
856 
859  template <class T>
860  static void insertion_sort(T* output, int32_t size)
861  {
862  for (int32_t i=0; i<size-1; i++)
863  {
864  int32_t j=i-1;
865  T value=output[i];
866  while (j >= 0 && output[j] > value)
867  {
868  output[j+1] = output[j];
869  j--;
870  }
871  output[j+1]=value;
872  }
873  }
874 
876  template <class T>
877  inline static void radix_sort(T* array, int32_t size)
878  {
879  radix_sort_helper(array,size,0);
880  }
881 
882  /*
883  * Inline function to extract the byte at position p (from left)
884  * of an 64 bit integer. The function is somewhat identical to
885  * accessing an array of characters via [].
886  */
887  template <class T>
888  static inline uint8_t byte(T word, uint16_t p)
889  {
890  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
891  }
892 
894  static inline uint8_t byte(complex128_t word, uint16_t p)
895  {
896  SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
897  return uint8_t(0);
898  }
899 
900  template <class T>
901  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
902  {
903  static size_t count[256], nc, cmin;
904  T *ak;
905  uint8_t c=0;
906  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
907  T *an, *aj, *pile[256];
908  size_t *cp, cmax;
909 
910  /* Push initial array to stack */
911  sp = s;
912  radix_push(array, size, i);
913 
914  /* Loop until all digits have been sorted */
915  while (sp>s) {
916  radix_pop(array, size, i);
917  an = array + size;
918 
919  /* Make character histogram */
920  if (nc == 0) {
921  cmin = 0xff;
922  for (ak = array; ak < an; ak++) {
923  c = byte(*ak, i);
924  count[c]++;
925  if (count[c] == 1) {
926  /* Determine smallest character */
927  if (c < cmin)
928  cmin = c;
929  nc++;
930  }
931  }
932 
933  /* Sort recursively if stack size too small */
934  if (sp + nc > s + RADIX_STACK_SIZE) {
935  radix_sort_helper(array, size, i);
936  continue;
937  }
938  }
939 
940  /*
941  * Set pile[]; push incompletely sorted bins onto stack.
942  * pile[] = pointers to last out-of-place element in bins.
943  * Before permuting: pile[c-1] + count[c] = pile[c];
944  * during deal: pile[c] counts down to pile[c-1].
945  */
946  olds = bigs = sp;
947  cmax = 2;
948  ak = array;
949  pile[0xff] = an;
950  for (cp = count + cmin; nc > 0; cp++) {
951  /* Find next non-empty pile */
952  while (*cp == 0)
953  cp++;
954  /* Pile with several entries */
955  if (*cp > 1) {
956  /* Determine biggest pile */
957  if (*cp > cmax) {
958  cmax = *cp;
959  bigs = sp;
960  }
961 
962  if (i < sizeof(T)-1)
963  radix_push(ak, *cp, (uint16_t) (i + 1));
964  }
965  pile[cp - count] = ak += *cp;
966  nc--;
967  }
968 
969  /* Play it safe -- biggest bin last. */
970  swap(*olds, *bigs);
971 
972  /*
973  * Permute misplacements home. Already home: everything
974  * before aj, and in pile[c], items from pile[c] on.
975  * Inner loop:
976  * r = next element to put in place;
977  * ak = pile[r[i]] = location to put the next element.
978  * aj = bottom of 1st disordered bin.
979  * Outer loop:
980  * Once the 1st disordered bin is done, ie. aj >= ak,
981  * aj<-aj + count[c] connects the bins in array linked list;
982  * reset count[c].
983  */
984  aj = array;
985  while (aj<an)
986  {
987  T r;
988 
989  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
990  swap(*ak, r);
991 
992  *aj = r;
993  aj += count[c];
994  count[c] = 0;
995  }
996  }
997  }
998 
1000  static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
1001  {
1002  SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
1003  }
1004 
1014  template <class T>
1015  static void qsort(T** vector, index_t length)
1016  {
1017  if (length<=1)
1018  return;
1019 
1020  if (length==2)
1021  {
1022  if (*vector[0]>*vector[1])
1023  swap(vector[0],vector[1]);
1024  return;
1025  }
1026  T* split=vector[length/2];
1027 
1028  int32_t left=0;
1029  int32_t right=length-1;
1030 
1031  while (left<=right)
1032  {
1033  while (*vector[left]<*split)
1034  ++left;
1035  while (*vector[right]>*split)
1036  --right;
1037 
1038  if (left<=right)
1039  {
1040  swap(vector[left],vector[right]);
1041  ++left;
1042  --right;
1043  }
1044  }
1045 
1046  if (right+1>1)
1047  qsort(vector,right+1);
1048 
1049  if (length-left>1)
1050  qsort(&vector[left],length-left);
1051  }
1052 
1054  static void qsort(complex128_t** vector, index_t length)
1055  {
1056  SG_SERROR("CMath::qsort():: Not supported for complex128_t\n");
1057  }
1058 
1060  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1061  {
1062  ASSERT(width>=0)
1063  for (int i=0; i<width; i++)
1064  {
1065  T mask = ((T) 1)<<(sizeof(T)*8-1);
1066  while (mask)
1067  {
1068  if (mask & word)
1069  SG_SPRINT("1")
1070  else
1071  SG_SPRINT("0")
1072 
1073  mask>>=1;
1074  }
1075  }
1076  }
1077 
1079  static void display_bits(complex128_t word,
1080  int32_t width=8*sizeof(complex128_t))
1081  {
1082  SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
1083  }
1084 
1090  template <class T1,class T2>
1091  static void qsort_index(T1* output, T2* index, uint32_t size);
1092 
1094  template <class T>
1095  static void qsort_index(complex128_t* output, T* index, uint32_t size)
1096  {
1097  SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
1098  }
1099 
1105  template <class T1,class T2>
1106  static void qsort_backward_index(
1107  T1* output, T2* index, int32_t size);
1108 
1110  template <class T>
1112  complex128_t* output, T* index, uint32_t size)
1113  {
1114  SG_SERROR("CMath::qsort_backword_index():: \
1115  Not supported for complex128_t\n");
1116  }
1117 
1125  template <class T1,class T2>
1126  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1127  {
1128  int32_t n=0;
1129  thread_qsort<T1,T2> t;
1130  t.output=output;
1131  t.index=index;
1132  t.size=size;
1133  t.qsort_threads=&n;
1134  t.sort_limit=limit;
1135  t.num_threads=n_threads;
1136  parallel_qsort_index<T1,T2>(&t);
1137  }
1138 
1140  template <class T>
1141  inline static void parallel_qsort_index(complex128_t* output, T* index,
1142  uint32_t size, int32_t n_threads, int32_t limit=0)
1143  {
1144  SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1145  }
1146 
1147 
1148  template <class T1,class T2>
1149  static void* parallel_qsort_index(void* p);
1150 
1151 
1152  /* finds the smallest element in output and puts that element as the
1153  first element */
1154  template <class T>
1155  static void min(float64_t* output, T* index, int32_t size);
1156 
1158  static void min(float64_t* output, complex128_t* index, int32_t size)
1159  {
1160  SG_SERROR("CMath::min():: Not supported for complex128_t\n");
1161  }
1162 
1163  /* finds the n smallest elements in output and puts these elements as the
1164  first n elements */
1165  template <class T>
1166  static void nmin(
1167  float64_t* output, T* index, int32_t size, int32_t n);
1168 
1170  static void nmin(float64_t* output, complex128_t* index,
1171  int32_t size, int32_t n)
1172  {
1173  SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
1174  }
1175 
1176 
1177 
1178  /* finds an element in a sorted array via binary search
1179  * returns -1 if not found */
1180  template <class T>
1181  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1182  {
1183  int32_t start=0;
1184  int32_t end=size-1;
1185 
1186  if (size<1)
1187  return 0;
1188 
1189  while (start<end)
1190  {
1191  int32_t middle=(start+end)/2;
1192 
1193  if (output[middle]>elem)
1194  end=middle-1;
1195  else if (output[middle]<elem)
1196  start=middle+1;
1197  else
1198  return middle;
1199  }
1200 
1201  return start;
1202  }
1203 
1205  static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
1206  {
1207  SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
1208  return int32_t(0);
1209  }
1210 
1211  /* finds an element in a sorted array via binary search
1212  * * returns -1 if not found */
1213  template <class T>
1214  static inline int32_t binary_search(T* output, int32_t size, T elem)
1215  {
1216  int32_t ind = binary_search_helper(output, size, elem);
1217  if (ind >= 0 && output[ind] == elem)
1218  return ind;
1219  return -1;
1220  }
1221 
1223  static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
1224  {
1225  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1226  return int32_t(-1);
1227  }
1228 
1229 
1230  /* Finds an element in a sorted array of pointers via binary search
1231  * Every element is dereferenced once before being compared
1232  *
1233  * @param array array of pointers to search in (assumed being sorted)
1234  * @param length length of array
1235  * @param elem pointer to element to search for
1236  * @return index of elem, -1 if not found */
1237  template<class T>
1238  static inline int32_t binary_search(T** vector, index_t length,
1239  T* elem)
1240  {
1241  int32_t start=0;
1242  int32_t end=length-1;
1243 
1244  if (length<1)
1245  return -1;
1246 
1247  while (start<end)
1248  {
1249  int32_t middle=(start+end)/2;
1250 
1251  if (*vector[middle]>*elem)
1252  end=middle-1;
1253  else if (*vector[middle]<*elem)
1254  start=middle+1;
1255  else
1256  {
1257  start=middle;
1258  break;
1259  }
1260  }
1261 
1262  if (start>=0&&*vector[start]==*elem)
1263  return start;
1264 
1265  return -1;
1266  }
1267 
1269  static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
1270  {
1271  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1272  return int32_t(-1);
1273  }
1274 
1275  template <class T>
1277  T* output, int32_t size, T elem)
1278  {
1279  int32_t ind = binary_search_helper(output, size, elem);
1280 
1281  if (output[ind]<=elem)
1282  return ind;
1283  if (ind>0 && output[ind-1] <= elem)
1284  return ind-1;
1285  return -1;
1286  }
1287 
1290  int32_t size, complex128_t elem)
1291  {
1292  SG_SERROR("CMath::binary_search_max_lower_equal():: \
1293  Not supported for complex128_t\n");
1294  return int32_t(-1);
1295  }
1296 
1299  static float64_t Align(
1300  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1301 
1302 
1304 
1307  {
1308  return c.real();
1309  }
1310 
1313  {
1314  return c.imag();
1315  }
1316 
1318  inline static uint32_t get_seed()
1319  {
1320  return CMath::seed;
1321  }
1322 
1324  inline static uint32_t get_log_range()
1325  {
1326  return CMath::LOGRANGE;
1327  }
1328 
1329 #ifdef USE_LOGCACHE
1330 
1331  inline static uint32_t get_log_accuracy()
1332  {
1333  return CMath::LOGACCURACY;
1334  }
1335 #endif
1336 
1338  static int is_finite(double f);
1339 
1341  static int is_infinity(double f);
1342 
1344  static int is_nan(double f);
1345 
1356 #ifdef USE_LOGCACHE
1357  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1358  {
1359  float64_t diff;
1360 
1361  if (!CMath::is_finite(p))
1362  return q;
1363 
1364  if (!CMath::is_finite(q))
1365  {
1366  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
1367  return NOT_A_NUMBER;
1368  }
1369  diff = p - q;
1370  if (diff > 0)
1371  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1372  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1373  }
1374 
1376  static void init_log_table();
1377 
1379  static int32_t determine_logrange();
1380 
1382  static int32_t determine_logaccuracy(int32_t range);
1383 #else
1385  float64_t p, float64_t q)
1386  {
1387  float64_t diff;
1388 
1389  if (!CMath::is_finite(p))
1390  return q;
1391  if (!CMath::is_finite(q))
1392  return p;
1393  diff = p - q;
1394  if (diff > 0)
1395  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
1396  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
1397  }
1398 #endif
1399 #ifdef USE_LOGSUMARRAY
1400 
1405  static inline float64_t logarithmic_sum_array(
1406  float64_t *p, int32_t len)
1407  {
1408  if (len<=2)
1409  {
1410  if (len==2)
1411  return logarithmic_sum(p[0],p[1]) ;
1412  if (len==1)
1413  return p[0];
1414  return -INFTY ;
1415  }
1416  else
1417  {
1418  float64_t *pp=p ;
1419  if (len%2==1) pp++ ;
1420  for (int32_t j=0; j < len>>1; j++)
1421  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
1422  }
1423  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1424  }
1425 #endif //USE_LOGSUMARRAY
1426 
1427 
1429  virtual const char* get_name() const { return "Math"; }
1430  public:
1433 
1434  static const float64_t NOT_A_NUMBER;
1436  static const float64_t INFTY;
1437  static const float64_t ALMOST_INFTY;
1438 
1441 
1443  static const float64_t PI;
1444 
1447 
1448  /* largest and smallest possible float64_t */
1451 
1452  /* Floating point Limits, Normalized */
1453  static const float32_t F_MAX_VAL32;
1455  static const float64_t F_MAX_VAL64;
1457 
1458  /* Floating point limits, Denormalized */
1459  static const float32_t F_MIN_VAL32;
1460  static const float64_t F_MIN_VAL64;
1461 
1462  protected:
1464  static int32_t LOGRANGE;
1465 
1467  static uint32_t seed;
1468 
1469 #ifdef USE_LOGCACHE
1470 
1472  static int32_t LOGACCURACY;
1474 
1475  static float64_t* logtable;
1476 #endif
1477 };
1478 
1479 //implementations of template functions
1480 template <class T1,class T2>
1482  {
1483  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1484  T1* output=ps->output;
1485  T2* index=ps->index;
1486  int32_t size=ps->size;
1487  int32_t* qsort_threads=ps->qsort_threads;
1488  int32_t sort_limit=ps->sort_limit;
1489  int32_t num_threads=ps->num_threads;
1490 
1491  if (size<2)
1492  {
1493  return NULL;
1494  }
1495 
1496  if (size==2)
1497  {
1498  if (output[0] > output [1])
1499  {
1500  swap(output[0], output[1]);
1501  swap(index[0], index[1]);
1502  }
1503  return NULL;
1504  }
1505  //T1 split=output[random(0,size-1)];
1506  T1 split=output[size/2];
1507 
1508  int32_t left=0;
1509  int32_t right=size-1;
1510 
1511  while (left<=right)
1512  {
1513  while (output[left] < split)
1514  left++;
1515  while (output[right] > split)
1516  right--;
1517 
1518  if (left<=right)
1519  {
1520  swap(output[left], output[right]);
1521  swap(index[left], index[right]);
1522  left++;
1523  right--;
1524  }
1525  }
1526  bool lthread_start=false;
1527  bool rthread_start=false;
1528  pthread_t lthread;
1529  pthread_t rthread;
1530  struct thread_qsort<T1,T2> t1;
1531  struct thread_qsort<T1,T2> t2;
1532 
1533  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1534  qsort_index(output,index,right+1);
1535  else if (right+1> 1)
1536  {
1537  (*qsort_threads)++;
1538  lthread_start=true;
1539  t1.output=output;
1540  t1.index=index;
1541  t1.size=right+1;
1542  t1.qsort_threads=qsort_threads;
1543  t1.sort_limit=sort_limit;
1544  t1.num_threads=num_threads;
1545  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1546  {
1547  lthread_start=false;
1548  (*qsort_threads)--;
1549  qsort_index(output,index,right+1);
1550  }
1551  }
1552 
1553 
1554  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1555  qsort_index(&output[left],&index[left], size-left);
1556  else if (size-left> 1)
1557  {
1558  (*qsort_threads)++;
1559  rthread_start=true;
1560  t2.output=&output[left];
1561  t2.index=&index[left];
1562  t2.size=size-left;
1563  t2.qsort_threads=qsort_threads;
1564  t2.sort_limit=sort_limit;
1565  t2.num_threads=num_threads;
1566  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1567  {
1568  rthread_start=false;
1569  (*qsort_threads)--;
1570  qsort_index(&output[left],&index[left], size-left);
1571  }
1572  }
1573 
1574  if (lthread_start)
1575  {
1576  pthread_join(lthread, NULL);
1577  (*qsort_threads)--;
1578  }
1579 
1580  if (rthread_start)
1581  {
1582  pthread_join(rthread, NULL);
1583  (*qsort_threads)--;
1584  }
1585 
1586  return NULL;
1587  }
1588 
1589  template <class T1,class T2>
1590 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
1591 {
1592  if (size<=1)
1593  return;
1594 
1595  if (size==2)
1596  {
1597  if (output[0] > output [1])
1598  {
1599  swap(output[0],output[1]);
1600  swap(index[0],index[1]);
1601  }
1602  return;
1603  }
1604  //T1 split=output[random(0,size-1)];
1605  T1 split=output[size/2];
1606 
1607  int32_t left=0;
1608  int32_t right=size-1;
1609 
1610  while (left<=right)
1611  {
1612  while (output[left] < split)
1613  left++;
1614  while (output[right] > split)
1615  right--;
1616 
1617  if (left<=right)
1618  {
1619  swap(output[left],output[right]);
1620  swap(index[left],index[right]);
1621  left++;
1622  right--;
1623  }
1624  }
1625 
1626  if (right+1> 1)
1627  qsort_index(output,index,right+1);
1628 
1629  if (size-left> 1)
1630  qsort_index(&output[left],&index[left], size-left);
1631 }
1632 
1633  template <class T1,class T2>
1634 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
1635 {
1636  if (size<=1)
1637  return;
1638 
1639  if (size==2)
1640  {
1641  if (output[0] < output [1])
1642  {
1643  swap(output[0],output[1]);
1644  swap(index[0],index[1]);
1645  }
1646  return;
1647  }
1648 
1649  //T1 split=output[random(0,size-1)];
1650  T1 split=output[size/2];
1651 
1652  int32_t left=0;
1653  int32_t right=size-1;
1654 
1655  while (left<=right)
1656  {
1657  while (output[left] > split)
1658  left++;
1659  while (output[right] < split)
1660  right--;
1661 
1662  if (left<=right)
1663  {
1664  swap(output[left],output[right]);
1665  swap(index[left],index[right]);
1666  left++;
1667  right--;
1668  }
1669  }
1670 
1671  if (right+1> 1)
1672  qsort_backward_index(output,index,right+1);
1673 
1674  if (size-left> 1)
1675  qsort_backward_index(&output[left],&index[left], size-left);
1676 }
1677 
1678  template <class T>
1679 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
1680 {
1681  if (6*n*size<13*size*CMath::log(size))
1682  for (int32_t i=0; i<n; i++)
1683  min(&output[i], &index[i], size-i);
1684  else
1685  qsort_index(output, index, size);
1686 }
1687 
1688 /* move the smallest entry in the array to the beginning */
1689  template <class T>
1690 void CMath::min(float64_t* output, T* index, int32_t size)
1691 {
1692  if (size<=1)
1693  return;
1694  float64_t min_elem=output[0];
1695  int32_t min_index=0;
1696  for (int32_t i=1; i<size; i++)
1697  {
1698  if (output[i]<min_elem)
1699  {
1700  min_index=i;
1701  min_elem=output[i];
1702  }
1703  }
1704  swap(output[0], output[min_index]);
1705  swap(index[0], index[min_index]);
1706 }
1707 
1708 #define COMPLEX128_ERROR_ONEARG_T(function) \
1709 template <> \
1710 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
1711 { \
1712  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1713  #function);\
1714  return complex128_t(0.0, 0.0); \
1715 }
1716 
1717 #define COMPLEX128_ERROR_TWOARGS_T(function) \
1718 template <> \
1719 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \
1720 { \
1721  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1722  #function);\
1723  return complex128_t(0.0, 0.0); \
1724 }
1725 
1726 #define COMPLEX128_ERROR_THREEARGS_T(function) \
1727 template <> \
1728 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \
1729 { \
1730  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1731  #function);\
1732  return complex128_t(0.0, 0.0); \
1733 }
1734 
1735 #define COMPLEX128_ERROR_SORT_T(function) \
1736 template <> \
1737 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \
1738 { \
1739  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1740  #function);\
1741 }
1742 
1745 
1746 
1748 
1751 
1753 // COMPLEX128_ERROR_ONEARG_T(sign)
1754 
1757 
1760 
1763 
1764 }
1765 #undef COMPLEX128_ERROR_ONEARG
1766 #undef COMPLEX128_ERROR_ONEARG_T
1767 #undef COMPLEX128_ERROR_TWOARGS_T
1768 #undef COMPLEX128_ERROR_THREEARGS_T
1769 #undef COMPLEX128_STDMATH
1770 #undef COMPLEX128_ERROR_SORT_T
1771 #endif

SHOGUN Machine Learning Toolbox - Documentation