SHOGUN  v2.0.0
 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) 2012 Fernando José Iglesias García
8  * Written (W) 2011 Siddharth Kherada
9  * Written (W) 2011 Justin Patera
10  * Written (W) 2011 Alesis Novik
11  * Written (W) 2011-2012 Heiko Strathmann
12  * Written (W) 1999-2009 Soeren Sonnenburg
13  * Written (W) 1999-2008 Gunnar Raetsch
14  * Written (W) 2007 Konrad Rieck
15  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
16  */
17 
18 #ifndef __MATHEMATICS_H_
19 #define __MATHEMATICS_H_
20 
21 #include <shogun/base/SGObject.h>
22 #include <shogun/lib/common.h>
23 #include <shogun/io/SGIO.h>
24 #include <shogun/base/Parallel.h>
25 
26 #include <math.h>
27 #include <stdio.h>
28 #include <float.h>
29 #include <pthread.h>
30 #include <unistd.h>
31 #include <sys/types.h>
32 #include <sys/time.h>
33 #include <time.h>
34 
35 #ifdef SUNOS
36 #include <ieeefp.h>
37 #endif
38 
40 #ifdef log2
41 #define cygwin_log2 log2
42 #undef log2
43 #endif
44 
45 
47 #ifdef _GLIBCXX_CMATH
48 #if _GLIBCXX_USE_C99_MATH
49 #if !_GLIBCXX_USE_C99_FP_MACROS_DYNAMIC
50 
52  using std::signbit;
53 
54  using std::fpclassify;
55 
56  using std::isfinite;
57  using std::isinf;
58  using std::isnan;
59  using std::isnormal;
60 
61  using std::isgreater;
62  using std::isgreaterequal;
63  using std::isless;
64  using std::islessequal;
65  using std::islessgreater;
66  using std::isunordered;
67 #endif
68 #endif
69 #endif
70 
71 
72 #ifdef _WIN32
73 #ifndef isnan
74 #define isnan _isnan
75 #endif
76 
77 #ifndef isfinite
78 #define isfinite _isfinite
79 #endif
80 #endif //_WIN32
81 
82 #ifndef NAN
83 #include <stdlib.h>
84 #define NAN (strtod("NAN",NULL))
85 #endif
86 
87 /* Size of RNG seed */
88 #define RNG_SEED_SIZE 256
89 
90 /* Maximum stack size */
91 #define RADIX_STACK_SIZE 512
92 
93 /* Stack macros */
94 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
95 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
96 
97 #ifndef DOXYGEN_SHOULD_SKIP_THIS
98 
99 template <class T> struct radix_stack_t
100 {
102  T *sa;
104  size_t sn;
106  uint16_t si;
107 };
108 
110 
112 template <class T1, class T2> struct thread_qsort
113 {
115  T1* output;
117  T2* index;
119  uint32_t size;
120 
122  int32_t* qsort_threads;
124  int32_t sort_limit;
126  int32_t num_threads;
127 };
128 #endif // DOXYGEN_SHOULD_SKIP_THIS
129 
130 namespace shogun
131 {
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  }
193 
196 
197  static inline float64_t round(float64_t d)
198  {
199  return ::floor(d+0.5);
200  }
201 
202  static inline float64_t floor(float64_t d)
203  {
204  return ::floor(d);
205  }
206 
207  static inline float64_t ceil(float64_t d)
208  {
209  return ::ceil(d);
210  }
211 
213  template <class T>
214  static inline T sign(T a)
215  {
216  if (a==0)
217  return 0;
218  else return (a<0) ? (-1) : (+1);
219  }
220 
222  template <class T>
223  static inline void swap(T &a,T &b)
224  {
225  /* register for fast swaps */
226  register T c=a;
227  a=b;
228  b=c;
229  }
230 
232  template <class T>
233  static inline T sq(T x)
234  {
235  return x*x;
236  }
237 
239  static inline float32_t sqrt(float32_t x)
240  {
241  return ::sqrtf(x);
242  }
243 
245  static inline float64_t sqrt(float64_t x)
246  {
247  return ::sqrt(x);
248  }
249 
251  static inline floatmax_t sqrt(floatmax_t x)
252  {
253  //fall back to double precision sqrt if sqrtl is not
254  //available
255 #ifdef HAVE_SQRTL
256  return ::sqrtl(x);
257 #else
258  return ::sqrt(x);
259 #endif
260  }
261 
263  static inline float32_t invsqrt(float32_t x)
264  {
265  union float_to_int
266  {
267  float32_t f;
268  int32_t i;
269  };
270 
271  float_to_int tmp;
272  tmp.f=x;
273 
274  float32_t xhalf = 0.5f * x;
275  // store floating-point bits in integer tmp.i
276  // and do initial guess for Newton's method
277  tmp.i = 0x5f3759d5 - (tmp.i >> 1);
278  x = tmp.f; // convert new bits into float
279  x = x*(1.5f - xhalf*x*x); // One round of Newton's method
280  return x;
281  }
282 
284  static inline floatmax_t powl(floatmax_t x, floatmax_t n)
285  {
286  //fall back to double precision pow if powl is not
287  //available
288 #ifdef HAVE_POWL
289  return ::powl((long double) x, (long double) n);
290 #else
291  return ::pow((double) x, (double) n);
292 #endif
293  }
294 
295  static inline int32_t pow(int32_t x, int32_t n)
296  {
297  ASSERT(n>=0);
298  int32_t result=1;
299  while (n--)
300  result*=x;
301 
302  return result;
303  }
304 
305  static inline float64_t pow(float64_t x, int32_t n)
306  {
307  if (n>=0)
308  {
309  float64_t result=1;
310  while (n--)
311  result*=x;
312 
313  return result;
314  }
315  else
316  return ::pow((double)x, (double)n);
317  }
318 
319  static inline float64_t pow(float64_t x, float64_t n)
320  {
321  return ::pow((double) x, (double) n);
322  }
323 
324  static inline float64_t exp(float64_t x)
325  {
326  return ::exp((double) x);
327  }
328 
330  static inline float64_t atan(float64_t x)
331  {
332  return ::atan((double) x);
333  }
334 
335  static inline float64_t log10(float64_t v)
336  {
337  return ::log(v)/::log(10.0);
338  }
339 
340  static inline float64_t log2(float64_t v)
341  {
342 #ifdef HAVE_LOG2
343  return ::log2(v);
344 #else
345  return ::log(v)/::log(2.0);
346 #endif //HAVE_LOG2
347  }
348 
349  static inline float64_t log(float64_t v)
350  {
351  return ::log(v);
352  }
353 
354  static inline float64_t sin(float64_t x)
355  {
356  return ::sin(x);
357  }
358 
359  static inline float64_t cos(float64_t x)
360  {
361  return ::cos(x);
362  }
363 
364  static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
365  {
366  ASSERT(len>0 && xy);
367 
368  float64_t area = 0.0;
369 
370  if (!reversed)
371  {
372  for (int i=1; i<len; i++)
373  area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
374  }
375  else
376  {
377  for (int i=1; i<len; i++)
378  area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
379  }
380 
381  return area;
382  }
383 
384  static inline int64_t factorial(int32_t n)
385  {
386  int64_t res=1;
387  for (int i=2; i<=n; i++)
388  res*=i ;
389  return res ;
390  }
391 
392  static void init_random(uint32_t initseed=0)
393  {
394  if (initseed==0)
395  {
396  struct timeval tv;
397  gettimeofday(&tv, NULL);
398  seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
399  }
400  else
401  seed=initseed;
402 #if !defined(CYGWIN) && !defined(__INTERIX)
403  //seed=42
404  //SG_SPRINT("initializing random number generator with %d (seed size %d)\n", seed, RNG_SEED_SIZE);
405  initstate(seed, CMath::rand_state, RNG_SEED_SIZE);
406 #endif
407  }
408 
409  static inline int64_t random()
410  {
411 #if defined(CYGWIN) || defined(__INTERIX)
412  return rand();
413 #else
414  return ::random();
415 #endif
416  }
417 
418  static inline uint64_t random(uint64_t min_value, uint64_t max_value)
419  {
420  uint64_t ret = min_value + (uint64_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
421  ASSERT(ret>=min_value && ret<=max_value);
422  return ret ;
423  }
424 
425  static inline int64_t random(int64_t min_value, int64_t max_value)
426  {
427  int64_t ret = min_value + (int64_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
428  ASSERT(ret>=min_value && ret<=max_value);
429  return ret ;
430  }
431 
432  static inline uint32_t random(uint32_t min_value, uint32_t max_value)
433  {
434  uint32_t ret = min_value + (uint32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
435  ASSERT(ret>=min_value && ret<=max_value);
436  return ret ;
437  }
438 
439  static inline int32_t random(int32_t min_value, int32_t max_value)
440  {
441  int32_t ret = min_value + (int32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
442  ASSERT(ret>=min_value && ret<=max_value);
443  return ret ;
444  }
445 
446  static inline float32_t random(float32_t min_value, float32_t max_value)
447  {
448  float32_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
449 
450  if (ret<min_value || ret>max_value)
451  SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
452  ASSERT(ret>=min_value && ret<=max_value);
453  return ret;
454  }
455 
456  static inline float64_t random(float64_t min_value, float64_t max_value)
457  {
458  float64_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
459 
460  if (ret<min_value || ret>max_value)
461  SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
462  ASSERT(ret>=min_value && ret<=max_value);
463  return ret;
464  }
465 
466  static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
467  {
468  floatmax_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
469 
470  if (ret<min_value || ret>max_value)
471  SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
472  ASSERT(ret>=min_value && ret<=max_value);
473  return ret;
474  }
475 
479  static inline float32_t normal_random(float32_t mean, float32_t std_dev)
480  {
481  // sets up variables & makes sure rand_s.range == (0,1)
482  float32_t ret;
483  float32_t rand_u;
484  float32_t rand_v;
485  float32_t rand_s;
486  do
487  {
488  rand_u = random(-1.0, 1.0);
489  rand_v = random(-1.0, 1.0);
490  rand_s = rand_u*rand_u + rand_v*rand_v;
491  } while ((rand_s == 0) || (rand_s >= 1));
492 
493  // the meat & potatos, and then the mean & standard deviation shifting...
494  ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
495  ret = std_dev*ret + mean;
496  return ret;
497  }
498 
502  static inline float64_t normal_random(float64_t mean, float64_t std_dev)
503  {
504  float64_t ret;
505  float64_t rand_u;
506  float64_t rand_v;
507  float64_t rand_s;
508  do
509  {
510  rand_u = random(-1.0, 1.0);
511  rand_v = random(-1.0, 1.0);
512  rand_s = rand_u*rand_u + rand_v*rand_v;
513  } while ((rand_s == 0) || (rand_s >= 1));
514 
515  ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
516  ret = std_dev*ret + mean;
517  return ret;
518  }
519 
522  static inline float32_t randn_float()
523  {
524  return normal_random(0.0, 1.0);
525  }
526 
529  static inline float64_t randn_double()
530  {
531  return normal_random(0.0, 1.0);
532  }
533 
534  template <class T>
535  static int32_t get_num_nonzero(T* vec, int32_t len)
536  {
537  int32_t nnz = 0;
538  for (index_t i=0; i<len; ++i)
539  nnz += vec[i] != 0;
540 
541  return nnz;
542  }
543 
548  static inline SGVector<int32_t> randperm_vec(int32_t n)
549  {
550  SGVector<int32_t> result(randperm(n), n);
551  return result;
552  }
553 
560  static inline int32_t* randperm(int32_t n)
561  {
562  int32_t* perm = SG_MALLOC(int32_t, n);
563 
564  if (!perm)
565  return NULL;
566  for (int32_t i = 0; i < n; i++)
567  perm[i] = i;
568  for (int32_t i = 0; i < n; i++)
569  swap(perm[random(0, n - 1)], perm[i]);
570  return perm;
571  }
572 
573  static inline int64_t nchoosek(int32_t n, int32_t k)
574  {
575  int64_t res=1;
576 
577  for (int32_t i=n-k+1; i<=n; i++)
578  res*=i;
579 
580  return res/factorial(k);
581  }
582 
586  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
587  static void sort(float64_t *a, int32_t*idx, int32_t N);
588 
589  /*
590  * Inline function to extract the byte at position p (from left)
591  * of an 64 bit integer. The function is somewhat identical to
592  * accessing an array of characters via [].
593  */
594 
596  template <class T>
597  inline static void radix_sort(T* array, int32_t size)
598  {
599  radix_sort_helper(array,size,0);
600  }
601 
602  template <class T>
603  static inline uint8_t byte(T word, uint16_t p)
604  {
605  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
606  }
607 
608  template <class T>
609  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
610  {
611  static size_t count[256], nc, cmin;
612  T *ak;
613  uint8_t c=0;
614  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
615  T *an, *aj, *pile[256];
616  size_t *cp, cmax;
617 
618  /* Push initial array to stack */
619  sp = s;
620  radix_push(array, size, i);
621 
622  /* Loop until all digits have been sorted */
623  while (sp>s) {
624  radix_pop(array, size, i);
625  an = array + size;
626 
627  /* Make character histogram */
628  if (nc == 0) {
629  cmin = 0xff;
630  for (ak = array; ak < an; ak++) {
631  c = byte(*ak, i);
632  count[c]++;
633  if (count[c] == 1) {
634  /* Determine smallest character */
635  if (c < cmin)
636  cmin = c;
637  nc++;
638  }
639  }
640 
641  /* Sort recursively if stack size too small */
642  if (sp + nc > s + RADIX_STACK_SIZE) {
643  radix_sort_helper(array, size, i);
644  continue;
645  }
646  }
647 
648  /*
649  * Set pile[]; push incompletely sorted bins onto stack.
650  * pile[] = pointers to last out-of-place element in bins.
651  * Before permuting: pile[c-1] + count[c] = pile[c];
652  * during deal: pile[c] counts down to pile[c-1].
653  */
654  olds = bigs = sp;
655  cmax = 2;
656  ak = array;
657  pile[0xff] = an;
658  for (cp = count + cmin; nc > 0; cp++) {
659  /* Find next non-empty pile */
660  while (*cp == 0)
661  cp++;
662  /* Pile with several entries */
663  if (*cp > 1) {
664  /* Determine biggest pile */
665  if (*cp > cmax) {
666  cmax = *cp;
667  bigs = sp;
668  }
669 
670  if (i < sizeof(T)-1)
671  radix_push(ak, *cp, (uint16_t) (i + 1));
672  }
673  pile[cp - count] = ak += *cp;
674  nc--;
675  }
676 
677  /* Play it safe -- biggest bin last. */
678  swap(*olds, *bigs);
679 
680  /*
681  * Permute misplacements home. Already home: everything
682  * before aj, and in pile[c], items from pile[c] on.
683  * Inner loop:
684  * r = next element to put in place;
685  * ak = pile[r[i]] = location to put the next element.
686  * aj = bottom of 1st disordered bin.
687  * Outer loop:
688  * Once the 1st disordered bin is done, ie. aj >= ak,
689  * aj<-aj + count[c] connects the bins in array linked list;
690  * reset count[c].
691  */
692  aj = array;
693  while (aj<an)
694  {
695  T r;
696 
697  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
698  swap(*ak, r);
699 
700  *aj = r;
701  aj += count[c];
702  count[c] = 0;
703  }
704  }
705  }
706 
709  template <class T>
710  static void insertion_sort(T* output, int32_t size)
711  {
712  for (int32_t i=0; i<size-1; i++)
713  {
714  int32_t j=i-1;
715  T value=output[i];
716  while (j >= 0 && output[j] > value)
717  {
718  output[j+1] = output[j];
719  j--;
720  }
721  output[j+1]=value;
722  }
723  }
724 
732  template <class T>
733  static index_t find_position_to_insert(SGVector<T> vector, T element)
734  {
735  index_t i;
736  for (i=0; i<vector.vlen; ++i)
737  {
738  if (vector[i]>element)
739  break;
740  }
741  return i;
742  }
743 
749  template <class T>
750  static void qsort(SGVector<T> vector)
751  {
752  qsort(vector.vector, vector.vlen);
753  }
754 
757  template <class T>
758  static void qsort(T* output, int32_t size)
759  {
760  if (size==1)
761  return;
762 
763  if (size==2)
764  {
765  if (output[0] > output [1])
766  swap(output[0],output[1]);
767  return;
768  }
769  //T split=output[random(0,size-1)];
770  T split=output[size/2];
771 
772  int32_t left=0;
773  int32_t right=size-1;
774 
775  while (left<=right)
776  {
777  while (output[left] < split)
778  left++;
779  while (output[right] > split)
780  right--;
781 
782  if (left<=right)
783  {
784  swap(output[left],output[right]);
785  left++;
786  right--;
787  }
788  }
789 
790  if (right+1> 1)
791  qsort(output,right+1);
792 
793  if (size-left> 1)
794  qsort(&output[left],size-left);
795  }
796 
806  template <class T>
807  static void qsort(T** vector, index_t length)
808  {
809  if (length==1)
810  return;
811 
812  if (length==2)
813  {
814  if (*vector[0]>*vector[1])
815  swap(vector[0],vector[1]);
816  return;
817  }
818  T* split=vector[length/2];
819 
820  int32_t left=0;
821  int32_t right=length-1;
822 
823  while (left<=right)
824  {
825  while (*vector[left]<*split)
826  ++left;
827  while (*vector[right]>*split)
828  --right;
829 
830  if (left<=right)
831  {
832  swap(vector[left],vector[right]);
833  ++left;
834  --right;
835  }
836  }
837 
838  if (right+1>1)
839  qsort(vector,right+1);
840 
841  if (length-left>1)
842  qsort(&vector[left],length-left);
843  }
844 
846  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
847  {
848  ASSERT(width>=0);
849  for (int i=0; i<width; i++)
850  {
851  T mask = ((T) 1)<<(sizeof(T)*8-1);
852  while (mask)
853  {
854  if (mask & word)
855  SG_SPRINT("1");
856  else
857  SG_SPRINT("0");
858 
859  mask>>=1;
860  }
861  }
862  }
863 
869  template <class T1,class T2>
870  static void qsort_index(T1* output, T2* index, uint32_t size);
871 
877  template <class T1,class T2>
878  static void qsort_backward_index(
879  T1* output, T2* index, int32_t size);
880 
888  template <class T1,class T2>
889  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
890  {
891  int32_t n=0;
892  thread_qsort<T1,T2> t;
893  t.output=output;
894  t.index=index;
895  t.size=size;
896  t.qsort_threads=&n;
897  t.sort_limit=limit;
898  t.num_threads=n_threads;
899  parallel_qsort_index<T1,T2>(&t);
900  }
901 
902 
903  template <class T1,class T2>
904  static void* parallel_qsort_index(void* p);
905 
906 
907  /* finds the smallest element in output and puts that element as the
908  first element */
909  template <class T>
910  static void min(float64_t* output, T* index, int32_t size);
911 
912  /* finds the n smallest elements in output and puts these elements as the
913  first n elements */
914  template <class T>
915  static void nmin(
916  float64_t* output, T* index, int32_t size, int32_t n);
917 
918 
919 
920 
921  /* finds an element in a sorted array via binary search
922  * returns -1 if not found */
923  template <class T>
924  static int32_t binary_search_helper(T* output, int32_t size, T elem)
925  {
926  int32_t start=0;
927  int32_t end=size-1;
928 
929  if (size<1)
930  return 0;
931 
932  while (start<end)
933  {
934  int32_t middle=(start+end)/2;
935 
936  if (output[middle]>elem)
937  end=middle-1;
938  else if (output[middle]<elem)
939  start=middle+1;
940  else
941  return middle;
942  }
943 
944  return start;
945  }
946 
947  /* finds an element in a sorted array via binary search
948  * * returns -1 if not found */
949  template <class T>
950  static inline int32_t binary_search(T* output, int32_t size, T elem)
951  {
952  int32_t ind = binary_search_helper(output, size, elem);
953  if (ind >= 0 && output[ind] == elem)
954  return ind;
955  return -1;
956  }
957 
958  /* Finds an element in a sorted array of pointers via binary search
959  * Every element is dereferenced once before being compared
960  *
961  * @param array array of pointers to search in (assumed being sorted)
962  * @param length length of array
963  * @param elem pointer to element to search for
964  * @return index of elem, -1 if not found */
965  template<class T>
966  static inline int32_t binary_search(T** vector, index_t length,
967  T* elem)
968  {
969  int32_t start=0;
970  int32_t end=length-1;
971 
972  if (length<1)
973  return -1;
974 
975  while (start<end)
976  {
977  int32_t middle=(start+end)/2;
978 
979  if (*vector[middle]>*elem)
980  end=middle-1;
981  else if (*vector[middle]<*elem)
982  start=middle+1;
983  else
984  {
985  start=middle;
986  break;
987  }
988  }
989 
990  if (start>=0&&*vector[start]==*elem)
991  return start;
992 
993  return -1;
994  }
995 
996  template <class T>
997  static int32_t binary_search_max_lower_equal(
998  T* output, int32_t size, T elem)
999  {
1000  int32_t ind = binary_search_helper(output, size, elem);
1001 
1002  if (output[ind]<=elem)
1003  return ind;
1004  if (ind>0 && output[ind-1] <= elem)
1005  return ind-1;
1006  return -1;
1007  }
1008 
1011  static float64_t Align(
1012  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1013 
1014 
1016 
1018  inline static uint32_t get_seed()
1019  {
1020  return CMath::seed;
1021  }
1022 
1024  inline static uint32_t get_log_range()
1025  {
1026  return CMath::LOGRANGE;
1027  }
1028 
1029 #ifdef USE_LOGCACHE
1030 
1031  inline static uint32_t get_log_accuracy()
1032  {
1033  return CMath::LOGACCURACY;
1034  }
1035 #endif
1036 
1038  inline static int is_finite(double f)
1039  {
1040 #if defined(isfinite) && !defined(SUNOS)
1041  return isfinite(f);
1042 #else
1043  return finite(f);
1044 #endif
1045  }
1046 
1048  inline static int is_infinity(double f)
1049  {
1050 #ifdef SUNOS
1051  if (fpclass(f) == FP_NINF || fpclass(f) == FP_PINF)
1052  return 1;
1053  else
1054  return 0;
1055 #else
1056  return isinf(f);
1057 #endif
1058  }
1059 
1061  inline static int is_nan(double f)
1062  {
1063 #ifdef SUNOS
1064  return isnand(f);
1065 #else
1066  return isnan(f);
1067 #endif
1068  }
1069 
1080 #ifdef USE_LOGCACHE
1081  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1082  {
1083  float64_t diff;
1084 
1085  if (!CMath::is_finite(p))
1086  return q;
1087 
1088  if (!CMath::is_finite(q))
1089  {
1090  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
1091  return NAN;
1092  }
1093  diff = p - q;
1094  if (diff > 0)
1095  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1096  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1097  }
1098 
1100  static void init_log_table();
1101 
1103  static int32_t determine_logrange();
1104 
1106  static int32_t determine_logaccuracy(int32_t range);
1107 #else
1108  static inline float64_t logarithmic_sum(
1109  float64_t p, float64_t q)
1110  {
1111  float64_t diff;
1112 
1113  if (!CMath::is_finite(p))
1114  return q;
1115  if (!CMath::is_finite(q))
1116  return p;
1117  diff = p - q;
1118  if (diff > 0)
1119  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
1120  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
1121  }
1122 #endif
1123 #ifdef USE_LOGSUMARRAY
1124 
1129  static inline float64_t logarithmic_sum_array(
1130  float64_t *p, int32_t len)
1131  {
1132  if (len<=2)
1133  {
1134  if (len==2)
1135  return logarithmic_sum(p[0],p[1]) ;
1136  if (len==1)
1137  return p[0];
1138  return -INFTY ;
1139  }
1140  else
1141  {
1142  register float64_t *pp=p ;
1143  if (len%2==1) pp++ ;
1144  for (register int32_t j=0; j < len>>1; j++)
1145  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
1146  }
1147  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1148  }
1149 #endif //USE_LOGSUMARRAY
1150 
1151 
1153  inline virtual const char* get_name() const { return "Mathematics"; }
1154  public:
1157 
1158  static const float64_t INFTY;
1159  static const float64_t ALMOST_INFTY;
1160 
1162  static const float64_t ALMOST_NEG_INFTY;
1163 
1165  static const float64_t PI;
1166 
1168  static const float64_t MACHINE_EPSILON;
1169 
1170  /* largest and smallest possible float64_t */
1171  static const float64_t MAX_REAL_NUMBER;
1172  static const float64_t MIN_REAL_NUMBER;
1173 
1174  protected:
1176  static int32_t LOGRANGE;
1177 
1179  static uint32_t seed;
1180  static char* rand_state;
1181 
1182 #ifdef USE_LOGCACHE
1183 
1185  static int32_t LOGACCURACY;
1187 
1188  static float64_t* logtable;
1189 #endif
1190 };
1191 
1192 //implementations of template functions
1193 template <class T1,class T2>
1194 void* CMath::parallel_qsort_index(void* p)
1195  {
1196  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1197  T1* output=ps->output;
1198  T2* index=ps->index;
1199  uint32_t size=ps->size;
1200  int32_t* qsort_threads=ps->qsort_threads;
1201  int32_t sort_limit=ps->sort_limit;
1202  int32_t num_threads=ps->num_threads;
1203 
1204  if (size==2)
1205  {
1206  if (output[0] > output [1])
1207  {
1208  swap(output[0], output[1]);
1209  swap(index[0], index[1]);
1210  }
1211  return NULL;
1212  }
1213  //T1 split=output[random(0,size-1)];
1214  T1 split=output[size/2];
1215 
1216  int32_t left=0;
1217  int32_t right=size-1;
1218 
1219  while (left<=right)
1220  {
1221  while (output[left] < split)
1222  left++;
1223  while (output[right] > split)
1224  right--;
1225 
1226  if (left<=right)
1227  {
1228  swap(output[left], output[right]);
1229  swap(index[left], index[right]);
1230  left++;
1231  right--;
1232  }
1233  }
1234  bool lthread_start=false;
1235  bool rthread_start=false;
1236  pthread_t lthread;
1237  pthread_t rthread;
1238  struct thread_qsort<T1,T2> t1;
1239  struct thread_qsort<T1,T2> t2;
1240 
1241  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1242  qsort_index(output,index,right+1);
1243  else if (right+1> 1)
1244  {
1245  (*qsort_threads)++;
1246  lthread_start=true;
1247  t1.output=output;
1248  t1.index=index;
1249  t1.size=right+1;
1250  t1.qsort_threads=qsort_threads;
1251  t1.sort_limit=sort_limit;
1252  t1.num_threads=num_threads;
1253  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1254  {
1255  lthread_start=false;
1256  (*qsort_threads)--;
1257  qsort_index(output,index,right+1);
1258  }
1259  }
1260 
1261 
1262  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1263  qsort_index(&output[left],&index[left], size-left);
1264  else if (size-left> 1)
1265  {
1266  (*qsort_threads)++;
1267  rthread_start=true;
1268  t2.output=&output[left];
1269  t2.index=&index[left];
1270  t2.size=size-left;
1271  t2.qsort_threads=qsort_threads;
1272  t2.sort_limit=sort_limit;
1273  t2.num_threads=num_threads;
1274  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1275  {
1276  rthread_start=false;
1277  (*qsort_threads)--;
1278  qsort_index(&output[left],&index[left], size-left);
1279  }
1280  }
1281 
1282  if (lthread_start)
1283  {
1284  pthread_join(lthread, NULL);
1285  (*qsort_threads)--;
1286  }
1287 
1288  if (rthread_start)
1289  {
1290  pthread_join(rthread, NULL);
1291  (*qsort_threads)--;
1292  }
1293 
1294  return NULL;
1295  }
1296 
1297  template <class T1,class T2>
1298 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
1299 {
1300  if (size==1)
1301  return;
1302 
1303  if (size==2)
1304  {
1305  if (output[0] > output [1])
1306  {
1307  swap(output[0],output[1]);
1308  swap(index[0],index[1]);
1309  }
1310  return;
1311  }
1312  //T1 split=output[random(0,size-1)];
1313  T1 split=output[size/2];
1314 
1315  int32_t left=0;
1316  int32_t right=size-1;
1317 
1318  while (left<=right)
1319  {
1320  while (output[left] < split)
1321  left++;
1322  while (output[right] > split)
1323  right--;
1324 
1325  if (left<=right)
1326  {
1327  swap(output[left],output[right]);
1328  swap(index[left],index[right]);
1329  left++;
1330  right--;
1331  }
1332  }
1333 
1334  if (right+1> 1)
1335  qsort_index(output,index,right+1);
1336 
1337  if (size-left> 1)
1338  qsort_index(&output[left],&index[left], size-left);
1339 }
1340 
1341  template <class T1,class T2>
1342 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
1343 {
1344  if (size==2)
1345  {
1346  if (output[0] < output [1])
1347  {
1348  swap(output[0],output[1]);
1349  swap(index[0],index[1]);
1350  }
1351  return;
1352  }
1353 
1354  //T1 split=output[random(0,size-1)];
1355  T1 split=output[size/2];
1356 
1357  int32_t left=0;
1358  int32_t right=size-1;
1359 
1360  while (left<=right)
1361  {
1362  while (output[left] > split)
1363  left++;
1364  while (output[right] < split)
1365  right--;
1366 
1367  if (left<=right)
1368  {
1369  swap(output[left],output[right]);
1370  swap(index[left],index[right]);
1371  left++;
1372  right--;
1373  }
1374  }
1375 
1376  if (right+1> 1)
1377  qsort_backward_index(output,index,right+1);
1378 
1379  if (size-left> 1)
1380  qsort_backward_index(&output[left],&index[left], size-left);
1381 }
1382 
1383  template <class T>
1384 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
1385 {
1386  if (6*n*size<13*size*CMath::log(size))
1387  for (int32_t i=0; i<n; i++)
1388  min(&output[i], &index[i], size-i) ;
1389  else
1390  qsort_index(output, index, size) ;
1391 }
1392 
1393 /* move the smallest entry in the array to the beginning */
1394  template <class T>
1395 void CMath::min(float64_t* output, T* index, int32_t size)
1396 {
1397  if (size<=1)
1398  return;
1399  float64_t min_elem=output[0];
1400  int32_t min_index=0;
1401  for (int32_t i=1; i<size; i++)
1402  {
1403  if (output[i]<min_elem)
1404  {
1405  min_index=i;
1406  min_elem=output[i];
1407  }
1408  }
1409  swap(output[0], output[min_index]);
1410  swap(index[0], index[min_index]);
1411 }
1412 }
1413 #endif
SHOGUN Machine Learning Toolbox - Documentation