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  /* register for fast swaps */
318  register T c=a;
319  a=b;
320  b=c;
321  }
322 
324  template <class T>
325  static inline T sq(T x)
326  {
327  return x*x;
328  }
329 
331  static inline float32_t sqrt(float32_t x)
332  {
333  return ::sqrtf(x);
334  }
335 
337  static inline float64_t sqrt(float64_t x)
338  {
339  return ::sqrt(x);
340  }
341 
343  static inline floatmax_t sqrt(floatmax_t x)
344  {
345  //fall back to double precision sqrt if sqrtl is not
346  //available
347 #ifdef HAVE_SQRTL
348  return ::sqrtl(x);
349 #else
350  return ::sqrt(x);
351 #endif
352  }
353 
356 
357 
358  static inline float32_t invsqrt(float32_t x)
359  {
360  union float_to_int
361  {
362  float32_t f;
363  int32_t i;
364  };
365 
366  float_to_int tmp;
367  tmp.f=x;
368 
369  float32_t xhalf = 0.5f * x;
370  // store floating-point bits in integer tmp.i
371  // and do initial guess for Newton's method
372  tmp.i = 0x5f3759d5 - (tmp.i >> 1);
373  x = tmp.f; // convert new bits into float
374  x = x*(1.5f - xhalf*x*x); // One round of Newton's method
375  return x;
376  }
377 
379  static inline floatmax_t powl(floatmax_t x, floatmax_t n)
380  {
381  //fall back to double precision pow if powl is not
382  //available
383 #ifdef HAVE_POWL
384  return ::powl((long double) x, (long double) n);
385 #else
386  return ::pow((double) x, (double) n);
387 #endif
388  }
389 
390  static inline int32_t pow(bool x, int32_t n)
391  {
392  return 0;
393  }
394 
395  static inline int32_t pow(int32_t x, int32_t n)
396  {
397  ASSERT(n>=0)
398  int32_t result=1;
399  while (n--)
400  result*=x;
401 
402  return result;
403  }
404 
405  static inline float64_t pow(float64_t x, int32_t n)
406  {
407  if (n>=0)
408  {
409  float64_t result=1;
410  while (n--)
411  result*=x;
412 
413  return result;
414  }
415  else
416  return ::pow((double)x, (double)n);
417  }
418 
419  static inline float64_t pow(float64_t x, float64_t n)
420  {
421  return ::pow((double) x, (double) n);
422  }
423 
425  static inline complex128_t pow(complex128_t x, int32_t n)
426  {
427  return std::pow(x, n);
428  }
429 
431  {
432  return std::pow(x, n);
433  }
434 
436  {
437  return std::pow(x, n);
438  }
439 
441  {
442  return std::pow(x, n);
443  }
444 
445  static inline float64_t exp(float64_t x)
446  {
447  return ::exp((double) x);
448  }
449 
452 
453 
454  static inline float64_t tan(float64_t x)
455  {
456  return ::tan((double) x);
457  }
458 
461 
462 
463  static inline float64_t atan(float64_t x)
464  {
465  return ::atan((double) x);
466  }
467 
470 
471 
472  static inline float64_t atan2(float64_t x, float64_t y)
473  {
474  return ::atan2((double) x, (double) y);
475  }
476 
479 
480 
481  static inline float64_t tanh(float64_t x)
482  {
483  return ::tanh((double) x);
484  }
485 
488 
489  static inline float64_t log10(float64_t v)
490  {
491  return ::log(v)/::log(10.0);
492  }
493 
496 
497  static inline float64_t log2(float64_t v)
498  {
499 #ifdef HAVE_LOG2
500  return ::log2(v);
501 #else
502  return ::log(v)/::log(2.0);
503 #endif //HAVE_LOG2
504  }
505 
506  static inline float64_t log(float64_t v)
507  {
508  return ::log(v);
509  }
510 
513 
514  static inline index_t floor_log(index_t n)
515  {
516  index_t i;
517  for (i = 0; n != 0; i++)
518  n >>= 1;
519 
520  return i;
521  }
522 
523  static inline float64_t sin(float64_t x)
524  {
525  return ::sin(x);
526  }
527 
530 
531  static inline float64_t asin(float64_t x)
532  {
533  return ::asin(x);
534  }
535 
538 
539  static inline float64_t sinh(float64_t x)
540  {
541  return ::sinh(x);
542  }
543 
546 
547  static inline float64_t cos(float64_t x)
548  {
549  return ::cos(x);
550  }
551 
554 
555  static inline float64_t acos(float64_t x)
556  {
557  return ::acos(x);
558  }
559 
562 
563  static inline float64_t cosh(float64_t x)
564  {
565  return ::cosh(x);
566  }
567 
570 
571  static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
572  {
573  ASSERT(len>0 && xy)
574 
575  float64_t area = 0.0;
576 
577  if (!reversed)
578  {
579  for (int i=1; i<len; i++)
580  area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
581  }
582  else
583  {
584  for (int i=1; i<len; i++)
585  area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
586  }
587 
588  return area;
589  }
590 
591  static bool strtof(const char* str, float32_t* float_result);
592  static bool strtod(const char* str, float64_t* double_result);
593  static bool strtold(const char* str, floatmax_t* long_double_result);
594 
595  static inline int64_t factorial(int32_t n)
596  {
597  int64_t res=1;
598  for (int i=2; i<=n; i++)
599  res*=i ;
600  return res ;
601  }
602 
603  static void init_random(uint32_t initseed=0)
604  {
605  if (initseed==0)
607  else
608  seed=initseed;
609 
611  }
612 
613  static inline uint64_t random()
614  {
615  return sg_rand->random_64();
616  }
617 
618  static inline uint64_t random(uint64_t min_value, uint64_t max_value)
619  {
620  return sg_rand->random(min_value, max_value);
621  }
622 
623  static inline int64_t random(int64_t min_value, int64_t max_value)
624  {
625  return sg_rand->random(min_value, max_value);
626  }
627 
628  static inline uint32_t random(uint32_t min_value, uint32_t max_value)
629  {
630  return sg_rand->random(min_value, max_value);
631  }
632 
633  static inline int32_t random(int32_t min_value, int32_t max_value)
634  {
635  return sg_rand->random(min_value, max_value);
636  }
637 
638  static inline float32_t random(float32_t min_value, float32_t max_value)
639  {
640  return sg_rand->random(min_value, max_value);
641  }
642 
643  static inline float64_t random(float64_t min_value, float64_t max_value)
644  {
645  return sg_rand->random(min_value, max_value);
646  }
647 
648  static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
649  {
650  return sg_rand->random(min_value, max_value);
651  }
652 
656  static inline float32_t normal_random(float32_t mean, float32_t std_dev)
657  {
658  // sets up variables & makes sure rand_s.range == (0,1)
659  float32_t ret;
660  float32_t rand_u;
661  float32_t rand_v;
662  float32_t rand_s;
663  do
664  {
665  rand_u = CMath::random(-1.0, 1.0);
666  rand_v = CMath::random(-1.0, 1.0);
667  rand_s = rand_u*rand_u + rand_v*rand_v;
668  } while ((rand_s == 0) || (rand_s >= 1));
669 
670  // the meat & potatos, and then the mean & standard deviation shifting...
671  ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
672  ret = std_dev*ret + mean;
673  return ret;
674  }
675 
679  static inline float64_t normal_random(float64_t mean, float64_t std_dev)
680  {
681  return sg_rand->normal_distrib(mean, std_dev);
682  }
683 
686  static inline float32_t randn_float()
687  {
688  return normal_random(0.0, 1.0);
689  }
690 
693  static inline float64_t randn_double()
694  {
695  return sg_rand->std_normal_distrib();
696  }
697 
704  template <class T>
705  static void permute(SGVector<T> v, CRandom* rand=NULL)
706  {
707  if (rand)
708  {
709  for (index_t i=0; i<v.vlen; ++i)
710  swap(v[i], v[rand->random(i, v.vlen-1)]);
711  }
712  else
713  {
714  for (index_t i=0; i<v.vlen; ++i)
715  swap(v[i], v[random(i, v.vlen-1)]);
716  }
717  }
718 
719  template <class T>
720  static int32_t get_num_nonzero(T* vec, int32_t len)
721  {
722  int32_t nnz = 0;
723  for (index_t i=0; i<len; ++i)
724  nnz += vec[i] != 0;
725 
726  return nnz;
727  }
728 
729  static int32_t get_num_nonzero(complex128_t* vec, int32_t len)
730  {
731  int32_t nnz=0;
732  for (index_t i=0; i<len; ++i)
733  nnz+=vec[i]!=0.0;
734  return nnz;
735  }
736 
737  static inline int64_t nchoosek(int32_t n, int32_t k)
738  {
739  int64_t res=1;
740 
741  for (int32_t i=n-k+1; i<=n; i++)
742  res*=i;
743 
744  return res/factorial(k);
745  }
746 
754  static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100);
755 
762  template <class T>
763  static T log_sum_exp(SGVector<T> values)
764  {
765  REQUIRE(values.vector, "Values are empty");
766 
767  /* find minimum element index */
768  index_t min_index=0;
769  T X0=values[0];
770  if (values.vlen>1)
771  {
772  for (index_t i=1; i<values.vlen; ++i)
773  {
774  if (values[i]<X0)
775  {
776  X0=values[i];
777  min_index=i;
778  }
779  }
780  }
781 
782  /* remove element from vector copy and compute log sum exp */
783  SGVector<T> values_without_X0(values.vlen-1);
784  index_t from_idx=0;
785  index_t to_idx=0;
786  for (from_idx=0; from_idx<values.vlen; ++from_idx)
787  {
788  if (from_idx!=min_index)
789  {
790  values_without_X0[to_idx]=exp(values[from_idx]-X0);
791  to_idx++;
792  }
793  }
794 
795  return X0+log(SGVector<T>::sum(values_without_X0)+1);
796  }
797 
804  template <class T>
805  static T log_mean_exp(SGVector<T> values)
806  {
807  return log_sum_exp(values) - log(values.vlen);
808  }
809 
813  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
814  static void sort(float64_t *a, int32_t*idx, int32_t N);
815 
818  template <class T>
819  static void qsort(T* output, int32_t size)
820  {
821  if (size<=1)
822  return;
823 
824  if (size==2)
825  {
826  if (output[0] > output [1])
827  CMath::swap(output[0],output[1]);
828  return;
829  }
830  //T split=output[random(0,size-1)];
831  T split=output[size/2];
832 
833  int32_t left=0;
834  int32_t right=size-1;
835 
836  while (left<=right)
837  {
838  while (output[left] < split)
839  left++;
840  while (output[right] > split)
841  right--;
842 
843  if (left<=right)
844  {
845  CMath::swap(output[left],output[right]);
846  left++;
847  right--;
848  }
849  }
850 
851  if (right+1> 1)
852  qsort(output,right+1);
853 
854  if (size-left> 1)
855  qsort(&output[left],size-left);
856  }
857 
860  template <class T>
861  static void insertion_sort(T* output, int32_t size)
862  {
863  for (int32_t i=0; i<size-1; i++)
864  {
865  int32_t j=i-1;
866  T value=output[i];
867  while (j >= 0 && output[j] > value)
868  {
869  output[j+1] = output[j];
870  j--;
871  }
872  output[j+1]=value;
873  }
874  }
875 
877  template <class T>
878  inline static void radix_sort(T* array, int32_t size)
879  {
880  radix_sort_helper(array,size,0);
881  }
882 
883  /*
884  * Inline function to extract the byte at position p (from left)
885  * of an 64 bit integer. The function is somewhat identical to
886  * accessing an array of characters via [].
887  */
888  template <class T>
889  static inline uint8_t byte(T word, uint16_t p)
890  {
891  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
892  }
893 
895  static inline uint8_t byte(complex128_t word, uint16_t p)
896  {
897  SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
898  return uint8_t(0);
899  }
900 
901  template <class T>
902  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
903  {
904  static size_t count[256], nc, cmin;
905  T *ak;
906  uint8_t c=0;
907  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
908  T *an, *aj, *pile[256];
909  size_t *cp, cmax;
910 
911  /* Push initial array to stack */
912  sp = s;
913  radix_push(array, size, i);
914 
915  /* Loop until all digits have been sorted */
916  while (sp>s) {
917  radix_pop(array, size, i);
918  an = array + size;
919 
920  /* Make character histogram */
921  if (nc == 0) {
922  cmin = 0xff;
923  for (ak = array; ak < an; ak++) {
924  c = byte(*ak, i);
925  count[c]++;
926  if (count[c] == 1) {
927  /* Determine smallest character */
928  if (c < cmin)
929  cmin = c;
930  nc++;
931  }
932  }
933 
934  /* Sort recursively if stack size too small */
935  if (sp + nc > s + RADIX_STACK_SIZE) {
936  radix_sort_helper(array, size, i);
937  continue;
938  }
939  }
940 
941  /*
942  * Set pile[]; push incompletely sorted bins onto stack.
943  * pile[] = pointers to last out-of-place element in bins.
944  * Before permuting: pile[c-1] + count[c] = pile[c];
945  * during deal: pile[c] counts down to pile[c-1].
946  */
947  olds = bigs = sp;
948  cmax = 2;
949  ak = array;
950  pile[0xff] = an;
951  for (cp = count + cmin; nc > 0; cp++) {
952  /* Find next non-empty pile */
953  while (*cp == 0)
954  cp++;
955  /* Pile with several entries */
956  if (*cp > 1) {
957  /* Determine biggest pile */
958  if (*cp > cmax) {
959  cmax = *cp;
960  bigs = sp;
961  }
962 
963  if (i < sizeof(T)-1)
964  radix_push(ak, *cp, (uint16_t) (i + 1));
965  }
966  pile[cp - count] = ak += *cp;
967  nc--;
968  }
969 
970  /* Play it safe -- biggest bin last. */
971  swap(*olds, *bigs);
972 
973  /*
974  * Permute misplacements home. Already home: everything
975  * before aj, and in pile[c], items from pile[c] on.
976  * Inner loop:
977  * r = next element to put in place;
978  * ak = pile[r[i]] = location to put the next element.
979  * aj = bottom of 1st disordered bin.
980  * Outer loop:
981  * Once the 1st disordered bin is done, ie. aj >= ak,
982  * aj<-aj + count[c] connects the bins in array linked list;
983  * reset count[c].
984  */
985  aj = array;
986  while (aj<an)
987  {
988  T r;
989 
990  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
991  swap(*ak, r);
992 
993  *aj = r;
994  aj += count[c];
995  count[c] = 0;
996  }
997  }
998  }
999 
1001  static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
1002  {
1003  SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
1004  }
1005 
1015  template <class T>
1016  static void qsort(T** vector, index_t length)
1017  {
1018  if (length<=1)
1019  return;
1020 
1021  if (length==2)
1022  {
1023  if (*vector[0]>*vector[1])
1024  swap(vector[0],vector[1]);
1025  return;
1026  }
1027  T* split=vector[length/2];
1028 
1029  int32_t left=0;
1030  int32_t right=length-1;
1031 
1032  while (left<=right)
1033  {
1034  while (*vector[left]<*split)
1035  ++left;
1036  while (*vector[right]>*split)
1037  --right;
1038 
1039  if (left<=right)
1040  {
1041  swap(vector[left],vector[right]);
1042  ++left;
1043  --right;
1044  }
1045  }
1046 
1047  if (right+1>1)
1048  qsort(vector,right+1);
1049 
1050  if (length-left>1)
1051  qsort(&vector[left],length-left);
1052  }
1053 
1055  static void qsort(complex128_t** vector, index_t length)
1056  {
1057  SG_SERROR("CMath::qsort():: Not supported for complex128_t\n");
1058  }
1059 
1061  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1062  {
1063  ASSERT(width>=0)
1064  for (int i=0; i<width; i++)
1065  {
1066  T mask = ((T) 1)<<(sizeof(T)*8-1);
1067  while (mask)
1068  {
1069  if (mask & word)
1070  SG_SPRINT("1")
1071  else
1072  SG_SPRINT("0")
1073 
1074  mask>>=1;
1075  }
1076  }
1077  }
1078 
1080  static void display_bits(complex128_t word,
1081  int32_t width=8*sizeof(complex128_t))
1082  {
1083  SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
1084  }
1085 
1091  template <class T1,class T2>
1092  static void qsort_index(T1* output, T2* index, uint32_t size);
1093 
1095  template <class T>
1096  static void qsort_index(complex128_t* output, T* index, uint32_t size)
1097  {
1098  SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
1099  }
1100 
1106  template <class T1,class T2>
1107  static void qsort_backward_index(
1108  T1* output, T2* index, int32_t size);
1109 
1111  template <class T>
1113  complex128_t* output, T* index, uint32_t size)
1114  {
1115  SG_SERROR("CMath::qsort_backword_index():: \
1116  Not supported for complex128_t\n");
1117  }
1118 
1126  template <class T1,class T2>
1127  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1128  {
1129  int32_t n=0;
1130  thread_qsort<T1,T2> t;
1131  t.output=output;
1132  t.index=index;
1133  t.size=size;
1134  t.qsort_threads=&n;
1135  t.sort_limit=limit;
1136  t.num_threads=n_threads;
1137  parallel_qsort_index<T1,T2>(&t);
1138  }
1139 
1141  template <class T>
1142  inline static void parallel_qsort_index(complex128_t* output, T* index,
1143  uint32_t size, int32_t n_threads, int32_t limit=0)
1144  {
1145  SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1146  }
1147 
1148 
1149  template <class T1,class T2>
1150  static void* parallel_qsort_index(void* p);
1151 
1152 
1153  /* finds the smallest element in output and puts that element as the
1154  first element */
1155  template <class T>
1156  static void min(float64_t* output, T* index, int32_t size);
1157 
1159  static void min(float64_t* output, complex128_t* index, int32_t size)
1160  {
1161  SG_SERROR("CMath::min():: Not supported for complex128_t\n");
1162  }
1163 
1164  /* finds the n smallest elements in output and puts these elements as the
1165  first n elements */
1166  template <class T>
1167  static void nmin(
1168  float64_t* output, T* index, int32_t size, int32_t n);
1169 
1171  static void nmin(float64_t* output, complex128_t* index,
1172  int32_t size, int32_t n)
1173  {
1174  SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
1175  }
1176 
1177 
1178 
1179  /* finds an element in a sorted array via binary search
1180  * returns -1 if not found */
1181  template <class T>
1182  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1183  {
1184  int32_t start=0;
1185  int32_t end=size-1;
1186 
1187  if (size<1)
1188  return 0;
1189 
1190  while (start<end)
1191  {
1192  int32_t middle=(start+end)/2;
1193 
1194  if (output[middle]>elem)
1195  end=middle-1;
1196  else if (output[middle]<elem)
1197  start=middle+1;
1198  else
1199  return middle;
1200  }
1201 
1202  return start;
1203  }
1204 
1206  static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
1207  {
1208  SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
1209  return int32_t(0);
1210  }
1211 
1212  /* finds an element in a sorted array via binary search
1213  * * returns -1 if not found */
1214  template <class T>
1215  static inline int32_t binary_search(T* output, int32_t size, T elem)
1216  {
1217  int32_t ind = binary_search_helper(output, size, elem);
1218  if (ind >= 0 && output[ind] == elem)
1219  return ind;
1220  return -1;
1221  }
1222 
1224  static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
1225  {
1226  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1227  return int32_t(-1);
1228  }
1229 
1230 
1231  /* Finds an element in a sorted array of pointers via binary search
1232  * Every element is dereferenced once before being compared
1233  *
1234  * @param array array of pointers to search in (assumed being sorted)
1235  * @param length length of array
1236  * @param elem pointer to element to search for
1237  * @return index of elem, -1 if not found */
1238  template<class T>
1239  static inline int32_t binary_search(T** vector, index_t length,
1240  T* elem)
1241  {
1242  int32_t start=0;
1243  int32_t end=length-1;
1244 
1245  if (length<1)
1246  return -1;
1247 
1248  while (start<end)
1249  {
1250  int32_t middle=(start+end)/2;
1251 
1252  if (*vector[middle]>*elem)
1253  end=middle-1;
1254  else if (*vector[middle]<*elem)
1255  start=middle+1;
1256  else
1257  {
1258  start=middle;
1259  break;
1260  }
1261  }
1262 
1263  if (start>=0&&*vector[start]==*elem)
1264  return start;
1265 
1266  return -1;
1267  }
1268 
1270  static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
1271  {
1272  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1273  return int32_t(-1);
1274  }
1275 
1276  template <class T>
1278  T* output, int32_t size, T elem)
1279  {
1280  int32_t ind = binary_search_helper(output, size, elem);
1281 
1282  if (output[ind]<=elem)
1283  return ind;
1284  if (ind>0 && output[ind-1] <= elem)
1285  return ind-1;
1286  return -1;
1287  }
1288 
1291  int32_t size, complex128_t elem)
1292  {
1293  SG_SERROR("CMath::binary_search_max_lower_equal():: \
1294  Not supported for complex128_t\n");
1295  return int32_t(-1);
1296  }
1297 
1300  static float64_t Align(
1301  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1302 
1303 
1305 
1308  {
1309  return c.real();
1310  }
1311 
1314  {
1315  return c.imag();
1316  }
1317 
1319  inline static uint32_t get_seed()
1320  {
1321  return CMath::seed;
1322  }
1323 
1325  inline static uint32_t get_log_range()
1326  {
1327  return CMath::LOGRANGE;
1328  }
1329 
1330 #ifdef USE_LOGCACHE
1331 
1332  inline static uint32_t get_log_accuracy()
1333  {
1334  return CMath::LOGACCURACY;
1335  }
1336 #endif
1337 
1339  static int is_finite(double f);
1340 
1342  static int is_infinity(double f);
1343 
1345  static int is_nan(double f);
1346 
1357 #ifdef USE_LOGCACHE
1358  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1359  {
1360  float64_t diff;
1361 
1362  if (!CMath::is_finite(p))
1363  return q;
1364 
1365  if (!CMath::is_finite(q))
1366  {
1367  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
1368  return NOT_A_NUMBER;
1369  }
1370  diff = p - q;
1371  if (diff > 0)
1372  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1373  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1374  }
1375 
1377  static void init_log_table();
1378 
1380  static int32_t determine_logrange();
1381 
1383  static int32_t determine_logaccuracy(int32_t range);
1384 #else
1386  float64_t p, float64_t q)
1387  {
1388  float64_t diff;
1389 
1390  if (!CMath::is_finite(p))
1391  return q;
1392  if (!CMath::is_finite(q))
1393  return p;
1394  diff = p - q;
1395  if (diff > 0)
1396  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
1397  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
1398  }
1399 #endif
1400 #ifdef USE_LOGSUMARRAY
1401 
1406  static inline float64_t logarithmic_sum_array(
1407  float64_t *p, int32_t len)
1408  {
1409  if (len<=2)
1410  {
1411  if (len==2)
1412  return logarithmic_sum(p[0],p[1]) ;
1413  if (len==1)
1414  return p[0];
1415  return -INFTY ;
1416  }
1417  else
1418  {
1419  register float64_t *pp=p ;
1420  if (len%2==1) pp++ ;
1421  for (register int32_t j=0; j < len>>1; j++)
1422  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
1423  }
1424  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1425  }
1426 #endif //USE_LOGSUMARRAY
1427 
1428 
1430  virtual const char* get_name() const { return "Math"; }
1431  public:
1434 
1435  static const float64_t NOT_A_NUMBER;
1437  static const float64_t INFTY;
1438  static const float64_t ALMOST_INFTY;
1439 
1442 
1444  static const float64_t PI;
1445 
1448 
1449  /* largest and smallest possible float64_t */
1452 
1453  /* Floating point Limits, Normalized */
1454  static const float32_t F_MAX_VAL32;
1456  static const float64_t F_MAX_VAL64;
1458 
1459  /* Floating point limits, Denormalized */
1460  static const float32_t F_MIN_VAL32;
1461  static const float64_t F_MIN_VAL64;
1462 
1463  protected:
1465  static int32_t LOGRANGE;
1466 
1468  static uint32_t seed;
1469 
1470 #ifdef USE_LOGCACHE
1471 
1473  static int32_t LOGACCURACY;
1475 
1476  static float64_t* logtable;
1477 #endif
1478 };
1479 
1480 //implementations of template functions
1481 template <class T1,class T2>
1483  {
1484  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1485  T1* output=ps->output;
1486  T2* index=ps->index;
1487  int32_t size=ps->size;
1488  int32_t* qsort_threads=ps->qsort_threads;
1489  int32_t sort_limit=ps->sort_limit;
1490  int32_t num_threads=ps->num_threads;
1491 
1492  if (size<2)
1493  {
1494  return NULL;
1495  }
1496 
1497  if (size==2)
1498  {
1499  if (output[0] > output [1])
1500  {
1501  swap(output[0], output[1]);
1502  swap(index[0], index[1]);
1503  }
1504  return NULL;
1505  }
1506  //T1 split=output[random(0,size-1)];
1507  T1 split=output[size/2];
1508 
1509  int32_t left=0;
1510  int32_t right=size-1;
1511 
1512  while (left<=right)
1513  {
1514  while (output[left] < split)
1515  left++;
1516  while (output[right] > split)
1517  right--;
1518 
1519  if (left<=right)
1520  {
1521  swap(output[left], output[right]);
1522  swap(index[left], index[right]);
1523  left++;
1524  right--;
1525  }
1526  }
1527  bool lthread_start=false;
1528  bool rthread_start=false;
1529  pthread_t lthread;
1530  pthread_t rthread;
1531  struct thread_qsort<T1,T2> t1;
1532  struct thread_qsort<T1,T2> t2;
1533 
1534  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1535  qsort_index(output,index,right+1);
1536  else if (right+1> 1)
1537  {
1538  (*qsort_threads)++;
1539  lthread_start=true;
1540  t1.output=output;
1541  t1.index=index;
1542  t1.size=right+1;
1543  t1.qsort_threads=qsort_threads;
1544  t1.sort_limit=sort_limit;
1545  t1.num_threads=num_threads;
1546  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1547  {
1548  lthread_start=false;
1549  (*qsort_threads)--;
1550  qsort_index(output,index,right+1);
1551  }
1552  }
1553 
1554 
1555  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1556  qsort_index(&output[left],&index[left], size-left);
1557  else if (size-left> 1)
1558  {
1559  (*qsort_threads)++;
1560  rthread_start=true;
1561  t2.output=&output[left];
1562  t2.index=&index[left];
1563  t2.size=size-left;
1564  t2.qsort_threads=qsort_threads;
1565  t2.sort_limit=sort_limit;
1566  t2.num_threads=num_threads;
1567  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1568  {
1569  rthread_start=false;
1570  (*qsort_threads)--;
1571  qsort_index(&output[left],&index[left], size-left);
1572  }
1573  }
1574 
1575  if (lthread_start)
1576  {
1577  pthread_join(lthread, NULL);
1578  (*qsort_threads)--;
1579  }
1580 
1581  if (rthread_start)
1582  {
1583  pthread_join(rthread, NULL);
1584  (*qsort_threads)--;
1585  }
1586 
1587  return NULL;
1588  }
1589 
1590  template <class T1,class T2>
1591 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
1592 {
1593  if (size<=1)
1594  return;
1595 
1596  if (size==2)
1597  {
1598  if (output[0] > output [1])
1599  {
1600  swap(output[0],output[1]);
1601  swap(index[0],index[1]);
1602  }
1603  return;
1604  }
1605  //T1 split=output[random(0,size-1)];
1606  T1 split=output[size/2];
1607 
1608  int32_t left=0;
1609  int32_t right=size-1;
1610 
1611  while (left<=right)
1612  {
1613  while (output[left] < split)
1614  left++;
1615  while (output[right] > split)
1616  right--;
1617 
1618  if (left<=right)
1619  {
1620  swap(output[left],output[right]);
1621  swap(index[left],index[right]);
1622  left++;
1623  right--;
1624  }
1625  }
1626 
1627  if (right+1> 1)
1628  qsort_index(output,index,right+1);
1629 
1630  if (size-left> 1)
1631  qsort_index(&output[left],&index[left], size-left);
1632 }
1633 
1634  template <class T1,class T2>
1635 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
1636 {
1637  if (size<=1)
1638  return;
1639 
1640  if (size==2)
1641  {
1642  if (output[0] < output [1])
1643  {
1644  swap(output[0],output[1]);
1645  swap(index[0],index[1]);
1646  }
1647  return;
1648  }
1649 
1650  //T1 split=output[random(0,size-1)];
1651  T1 split=output[size/2];
1652 
1653  int32_t left=0;
1654  int32_t right=size-1;
1655 
1656  while (left<=right)
1657  {
1658  while (output[left] > split)
1659  left++;
1660  while (output[right] < split)
1661  right--;
1662 
1663  if (left<=right)
1664  {
1665  swap(output[left],output[right]);
1666  swap(index[left],index[right]);
1667  left++;
1668  right--;
1669  }
1670  }
1671 
1672  if (right+1> 1)
1673  qsort_backward_index(output,index,right+1);
1674 
1675  if (size-left> 1)
1676  qsort_backward_index(&output[left],&index[left], size-left);
1677 }
1678 
1679  template <class T>
1680 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
1681 {
1682  if (6*n*size<13*size*CMath::log(size))
1683  for (int32_t i=0; i<n; i++)
1684  min(&output[i], &index[i], size-i);
1685  else
1686  qsort_index(output, index, size);
1687 }
1688 
1689 /* move the smallest entry in the array to the beginning */
1690  template <class T>
1691 void CMath::min(float64_t* output, T* index, int32_t size)
1692 {
1693  if (size<=1)
1694  return;
1695  float64_t min_elem=output[0];
1696  int32_t min_index=0;
1697  for (int32_t i=1; i<size; i++)
1698  {
1699  if (output[i]<min_elem)
1700  {
1701  min_index=i;
1702  min_elem=output[i];
1703  }
1704  }
1705  swap(output[0], output[min_index]);
1706  swap(index[0], index[min_index]);
1707 }
1708 
1709 #define COMPLEX128_ERROR_ONEARG_T(function) \
1710 template <> \
1711 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
1712 { \
1713  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1714  #function);\
1715  return complex128_t(0.0, 0.0); \
1716 }
1717 
1718 #define COMPLEX128_ERROR_TWOARGS_T(function) \
1719 template <> \
1720 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \
1721 { \
1722  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1723  #function);\
1724  return complex128_t(0.0, 0.0); \
1725 }
1726 
1727 #define COMPLEX128_ERROR_THREEARGS_T(function) \
1728 template <> \
1729 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \
1730 { \
1731  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1732  #function);\
1733  return complex128_t(0.0, 0.0); \
1734 }
1735 
1736 #define COMPLEX128_ERROR_SORT_T(function) \
1737 template <> \
1738 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \
1739 { \
1740  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1741  #function);\
1742 }
1743 
1746 
1747 
1749 
1752 
1754 // COMPLEX128_ERROR_ONEARG_T(sign)
1755 
1758 
1761 
1764 
1765 }
1766 #undef COMPLEX128_ERROR_ONEARG
1767 #undef COMPLEX128_ERROR_ONEARG_T
1768 #undef COMPLEX128_ERROR_TWOARGS_T
1769 #undef COMPLEX128_ERROR_THREEARGS_T
1770 #undef COMPLEX128_STDMATH
1771 #undef COMPLEX128_ERROR_SORT_T
1772 #endif

SHOGUN Machine Learning Toolbox - Documentation