SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros 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) 2011 Siddharth Kherada
8  * Written (W) 2011 Justin Patera
9  * Written (W) 2011 Alesis Novik
10  * Written (W) 1999-2009 Soeren Sonnenburg
11  * Written (W) 1999-2008 Gunnar Raetsch
12  * Written (W) 2007 Konrad Rieck
13  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
14  */
15 
16 #ifndef __MATHEMATICS_H_
17 #define __MATHEMATICS_H_
18 
19 #include <shogun/lib/common.h>
20 #include <shogun/io/SGIO.h>
22 #include <shogun/base/SGObject.h>
23 #include <shogun/base/Parallel.h>
24 #include <shogun/lib/DataType.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 
46 
48 #ifdef _GLIBCXX_CMATH
49 #if _GLIBCXX_USE_C99_MATH
50 #if !_GLIBCXX_USE_C99_FP_MACROS_DYNAMIC
51 
53  using std::signbit;
54 
55  using std::fpclassify;
56 
57  using std::isfinite;
58  using std::isinf;
59  using std::isnan;
60  using std::isnormal;
61 
62  using std::isgreater;
63  using std::isgreaterequal;
64  using std::isless;
65  using std::islessequal;
66  using std::islessgreater;
67  using std::isunordered;
68 #endif
69 #endif
70 #endif
71 
72 
73 #ifdef _WIN32
74 #ifndef isnan
75 #define isnan _isnan
76 #endif
77 
78 #ifndef isfinite
79 #define isfinite _isfinite
80 #endif
81 #endif //_WIN32
82 
83 #ifndef NAN
84 #include <stdlib.h>
85 #define NAN (strtod("NAN",NULL))
86 #endif
87 
88 /* Size of RNG seed */
89 #define RNG_SEED_SIZE 256
90 
91 /* Maximum stack size */
92 #define RADIX_STACK_SIZE 512
93 
94 /* Stack macros */
95 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
96 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
97 
98 #ifndef DOXYGEN_SHOULD_SKIP_THIS
99 
100 template <class T> struct radix_stack_t
101 {
103  T *sa;
105  size_t sn;
107  uint16_t si;
108 };
109 
111 
113 template <class T1, class T2> struct thread_qsort
114 {
116  T1* output;
118  T2* index;
120  uint32_t size;
121 
123  int32_t* qsort_threads;
125  int32_t sort_limit;
127  int32_t num_threads;
128 };
129 #endif // DOXYGEN_SHOULD_SKIP_THIS
130 
131 namespace shogun
132 {
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  T c=a;
226  a=b;
227  b=c;
228  }
229 
233  template <class T>
234  static inline void resize(T* &data, int64_t old_size, int64_t new_size)
235  {
236  if (old_size==new_size)
237  return;
238  T* new_data = SG_MALLOC(T, new_size);
239  for (int64_t i=0; i<old_size && i<new_size; i++)
240  new_data[i]=data[i];
241  SG_FREE(data);
242  data=new_data;
243  }
244 
246  template <class T>
247  static inline T twonorm(T* x, int32_t len)
248  {
249  float64_t result=0;
250  for (int32_t i=0; i<len; i++)
251  result+=x[i]*x[i];
252 
253  return CMath::sqrt(result);
254  }
255 
257  template <class T>
258  static inline T qsq(T* x, int32_t len, float64_t q)
259  {
260  float64_t result=0;
261  for (int32_t i=0; i<len; i++)
262  result+=CMath::pow(x[i], q);
263 
264  return result;
265  }
266 
268  template <class T>
269  static inline T qnorm(T* x, int32_t len, float64_t q)
270  {
271  ASSERT(q!=0);
272  return CMath::pow(qsq(x, len, q), 1/q);
273  }
274 
276  template <class T>
277  static inline T sq(T x)
278  {
279  return x*x;
280  }
281 
283  static inline float32_t sqrt(float32_t x)
284  {
285  return ::sqrtf(x);
286  }
287 
289  static inline float64_t sqrt(float64_t x)
290  {
291  return ::sqrt(x);
292  }
293 
295  static inline floatmax_t sqrt(floatmax_t x)
296  {
297  //fall back to double precision sqrt if sqrtl is not
298  //available
299 #ifdef HAVE_SQRTL
300  return ::sqrtl(x);
301 #else
302  return ::sqrt(x);
303 #endif
304  }
305 
307  static inline float32_t invsqrt(float32_t x)
308  {
309  union float_to_int
310  {
311  float32_t f;
312  int32_t i;
313  };
314 
315  float_to_int tmp;
316  tmp.f=x;
317 
318  float32_t xhalf = 0.5f * x;
319  // store floating-point bits in integer tmp.i
320  // and do initial guess for Newton's method
321  tmp.i = 0x5f3759d5 - (tmp.i >> 1);
322  x = tmp.f; // convert new bits into float
323  x = x*(1.5f - xhalf*x*x); // One round of Newton's method
324  return x;
325  }
326 
328  static inline floatmax_t powl(floatmax_t x, floatmax_t n)
329  {
330  //fall back to double precision pow if powl is not
331  //available
332 #ifdef HAVE_POWL
333  return ::powl((long double) x, (long double) n);
334 #else
335  return ::pow((double) x, (double) n);
336 #endif
337  }
338 
339  static inline int32_t pow(int32_t x, int32_t n)
340  {
341  ASSERT(n>=0);
342  int32_t result=1;
343  while (n--)
344  result*=x;
345 
346  return result;
347  }
348 
349  static inline float64_t pow(float64_t x, int32_t n)
350  {
351  if (n>=0)
352  {
353  float64_t result=1;
354  while (n--)
355  result*=x;
356 
357  return result;
358  }
359  else
360  return ::pow((double)x, (double)n);
361  }
362 
363  static inline float64_t pow(float64_t x, float64_t n)
364  {
365  return ::pow((double) x, (double) n);
366  }
367 
368  static inline float64_t exp(float64_t x)
369  {
370  return ::exp((double) x);
371  }
372 
374  static inline float64_t lgamma(float64_t x)
375  {
376  return ::lgamma((double) x);
377  }
378 
380  static inline float64_t tgamma(float64_t x)
381  {
382  return ::tgamma((double) x);
383  }
384 
386  static inline float64_t atan(float64_t x)
387  {
388  return ::atan((double) x);
389  }
390 
391  static inline floatmax_t lgammal(floatmax_t x)
392  {
393 #ifdef HAVE_LGAMMAL
394  return ::lgammal((long double) x);
395 #else
396  return ::lgamma((double) x);
397 #endif // HAVE_LGAMMAL
398  }
399 
400  static inline float64_t log10(float64_t v)
401  {
402  return ::log(v)/::log(10.0);
403  }
404 
405  static inline float64_t log2(float64_t v)
406  {
407 #ifdef HAVE_LOG2
408  return ::log2(v);
409 #else
410  return ::log(v)/::log(2.0);
411 #endif //HAVE_LOG2
412  }
413 
414  static inline float64_t log(float64_t v)
415  {
416  return ::log(v);
417  }
418 
419  static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
420  {
421  ASSERT(len>0 && xy);
422 
423  float64_t area = 0.0;
424 
425  if (!reversed)
426  {
427  for (int i=1; i<len; i++)
428  area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
429  }
430  else
431  {
432  for (int i=1; i<len; i++)
433  area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
434  }
435 
436  return area;
437  }
438 
439  template <class T>
440  static void transpose_matrix(
441  T*& matrix, int32_t& num_feat, int32_t& num_vec)
442  {
443  T* transposed=SG_MALLOC(T, num_vec*num_feat);
444  for (int32_t i=0; i<num_vec; i++)
445  {
446  for (int32_t j=0; j<num_feat; j++)
447  transposed[i+j*num_vec]=matrix[i*num_feat+j];
448  }
449 
450  SG_FREE(matrix);
451  matrix=transposed;
452 
453  CMath::swap(num_feat, num_vec);
454  }
455 
456 #ifdef HAVE_LAPACK
457 
458 
459  static float64_t* pinv(
460  float64_t* matrix, int32_t rows, int32_t cols,
461  float64_t* target=NULL);
462 
463 
464  //C := alpha*op( A )*op( B ) + beta*C
465  //op( X ) = X or op( X ) = X',
466  static inline void dgemm(
467  double alpha, const double* A, int rows, int cols,
468  CBLAS_TRANSPOSE transposeA, double *B, int cols_B,
469  CBLAS_TRANSPOSE transposeB, double beta, double *C)
470  {
471  cblas_dgemm(CblasColMajor, transposeA, transposeB, rows, cols, cols_B,
472  alpha, A, cols, B, cols_B, beta, C, cols);
473  }
474 
475  //y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
476  static inline void dgemv(
477  double alpha, const double *A, int rows, int cols,
478  const CBLAS_TRANSPOSE transposeA, const double* X, double beta,
479  double* Y)
480  {
481  cblas_dgemv(CblasColMajor, transposeA,
482  rows, cols, alpha, A, cols,
483  X, 1, beta, Y, 1);
484  }
485 #endif
486 
487  static inline int64_t factorial(int32_t n)
488  {
489  int64_t res=1;
490  for (int i=2; i<=n; i++)
491  res*=i ;
492  return res ;
493  }
494 
495  static void init_random(uint32_t initseed=0)
496  {
497  if (initseed==0)
498  {
499  struct timeval tv;
500  gettimeofday(&tv, NULL);
501  seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
502  }
503  else
504  seed=initseed;
505 #if !defined(CYGWIN) && !defined(__INTERIX)
506  //seed=42
507  //SG_SPRINT("initializing random number generator with %d (seed size %d)\n", seed, RNG_SEED_SIZE);
508  initstate(seed, CMath::rand_state, RNG_SEED_SIZE);
509 #endif
510  }
511 
512  static inline int64_t random()
513  {
514 #if defined(CYGWIN) || defined(__INTERIX)
515  return rand();
516 #else
517  return ::random();
518 #endif
519  }
520 
521  static inline int32_t random(int32_t min_value, int32_t max_value)
522  {
523  int32_t ret = min_value + (int32_t) ((max_value-min_value+1) * (random() / (RAND_MAX+1.0)));
524  ASSERT(ret>=min_value && ret<=max_value);
525  return ret ;
526  }
527 
528  static inline float32_t random(float32_t min_value, float32_t max_value)
529  {
530  float32_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
531 
532  if (ret<min_value || ret>max_value)
533  SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
534  ASSERT(ret>=min_value && ret<=max_value);
535  return ret;
536  }
537 
538  static inline float64_t random(float64_t min_value, float64_t max_value)
539  {
540  float64_t ret = min_value + ((max_value-min_value) * (random() / (1.0*RAND_MAX)));
541 
542  if (ret<min_value || ret>max_value)
543  SG_SPRINT("min_value:%10.10f value: %10.10f max_value:%10.10f", min_value, ret, max_value);
544  ASSERT(ret>=min_value && ret<=max_value);
545  return ret;
546  }
547 
552  {
553  // sets up variables & makes sure rand_s.range == (0,1)
554  float32_t ret;
555  float32_t rand_u;
556  float32_t rand_v;
557  float32_t rand_s;
558  do
559  {
560  rand_u = random(-1.0, 1.0);
561  rand_v = random(-1.0, 1.0);
562  rand_s = rand_u*rand_u + rand_v*rand_v;
563  } while ((rand_s == 0) || (rand_s >= 1));
564 
565  // the meat & potatos, and then the mean & standard deviation shifting...
566  ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
567  ret = std_dev*ret + mean;
568  return ret;
569  }
570 
575  {
576  float64_t ret;
577  float64_t rand_u;
578  float64_t rand_v;
579  float64_t rand_s;
580  do
581  {
582  rand_u = random(-1.0, 1.0);
583  rand_v = random(-1.0, 1.0);
584  rand_s = rand_u*rand_u + rand_v*rand_v;
585  } while ((rand_s == 0) || (rand_s >= 1));
586 
587  ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
588  ret = std_dev*ret + mean;
589  return ret;
590  }
591 
594  static inline float32_t randn_float()
595  {
596  return normal_random(0.0, 1.0);
597  }
598 
601  static inline float64_t randn_double()
602  {
603  return normal_random(0.0, 1.0);
604  }
605 
606  template <class T>
607  static T* clone_vector(const T* vec, int32_t len)
608  {
609  T* result = SG_MALLOC(T, len);
610  for (int32_t i=0; i<len; i++)
611  result[i]=vec[i];
612 
613  return result;
614  }
615  template <class T>
616  static void fill_vector(T* vec, int32_t len, T value)
617  {
618  for (int32_t i=0; i<len; i++)
619  vec[i]=value;
620  }
621  template <class T>
622  static void range_fill_vector(T* vec, int32_t len, T start=0)
623  {
624  for (int32_t i=0; i<len; i++)
625  vec[i]=i+start;
626  }
627 
628  template <class T>
629  static void random_vector(T* vec, int32_t len, T min_value, T max_value)
630  {
631  for (int32_t i=0; i<len; i++)
632  vec[i]=CMath::random(min_value, max_value);
633  }
634 
635  static inline int32_t* randperm(int32_t n)
636  {
637  int32_t* perm = SG_MALLOC(int32_t, n);
638 
639  if (!perm)
640  return NULL;
641  for (int32_t i = 0; i < n; i++)
642  perm[i] = i;
643  for (int32_t i = 0; i < n; i++)
644  swap(perm[random(0, n - 1)], perm[i]);
645  return perm;
646  }
647 
648  static inline int64_t nchoosek(int32_t n, int32_t k)
649  {
650  int64_t res=1;
651 
652  for (int32_t i=n-k+1; i<=n; i++)
653  res*=i;
654 
655  return res/factorial(k);
656  }
657 
659  template <class T>
660  static inline void vec1_plus_scalar_times_vec2(T* vec1,
661  T scalar, const T* vec2, int32_t n)
662  {
663  for (int32_t i=0; i<n; i++)
664  vec1[i]+=scalar*vec2[i];
665  }
666 
668  static inline float64_t dot(const bool* v1, const bool* v2, int32_t n)
669  {
670  float64_t r=0;
671  for (int32_t i=0; i<n; i++)
672  r+=((v1[i]) ? 1 : 0) * ((v2[i]) ? 1 : 0);
673  return r;
674  }
675 
677  static inline floatmax_t dot(const floatmax_t* v1, const floatmax_t* v2, int32_t n)
678  {
679  floatmax_t r=0;
680  for (int32_t i=0; i<n; i++)
681  r+=v1[i]*v2[i];
682  return r;
683  }
684 
685 
687  static inline float64_t dot(const float64_t* v1, const float64_t* v2, int32_t n)
688  {
689  float64_t r=0;
690 #ifdef HAVE_LAPACK
691  int32_t skip=1;
692  r = cblas_ddot(n, v1, skip, v2, skip);
693 #else
694  for (int32_t i=0; i<n; i++)
695  r+=v1[i]*v2[i];
696 #endif
697  return r;
698  }
699 
701  static inline float32_t dot(
702  const float32_t* v1, const float32_t* v2, int32_t n)
703  {
704  float64_t r=0;
705 #ifdef HAVE_LAPACK
706  int32_t skip=1;
707  r = cblas_sdot(n, v1, skip, v2, skip);
708 #else
709  for (int32_t i=0; i<n; i++)
710  r+=v1[i]*v2[i];
711 #endif
712  return r;
713  }
714 
716  static inline float64_t dot(
717  const uint64_t* v1, const uint64_t* v2, int32_t n)
718  {
719  float64_t r=0;
720  for (int32_t i=0; i<n; i++)
721  r+=((float64_t) v1[i])*v2[i];
722 
723  return r;
724  }
726  static inline float64_t dot(
727  const int64_t* v1, const int64_t* v2, int32_t n)
728  {
729  float64_t r=0;
730  for (int32_t i=0; i<n; i++)
731  r+=((float64_t) v1[i])*v2[i];
732 
733  return r;
734  }
735 
737  static inline float64_t dot(
738  const int32_t* v1, const int32_t* v2, int32_t n)
739  {
740  float64_t r=0;
741  for (int32_t i=0; i<n; i++)
742  r+=((float64_t) v1[i])*v2[i];
743 
744  return r;
745  }
746 
748  static inline float64_t dot(
749  const uint32_t* v1, const uint32_t* v2, int32_t n)
750  {
751  float64_t r=0;
752  for (int32_t i=0; i<n; i++)
753  r+=((float64_t) v1[i])*v2[i];
754 
755  return r;
756  }
757 
759  static inline float64_t dot(
760  const uint16_t* v1, const uint16_t* v2, int32_t n)
761  {
762  float64_t r=0;
763  for (int32_t i=0; i<n; i++)
764  r+=((float64_t) v1[i])*v2[i];
765 
766  return r;
767  }
768 
770  static inline float64_t dot(
771  const int16_t* v1, const int16_t* v2, int32_t n)
772  {
773  float64_t r=0;
774  for (int32_t i=0; i<n; i++)
775  r+=((float64_t) v1[i])*v2[i];
776 
777  return r;
778  }
779 
781  static inline float64_t dot(
782  const char* v1, const char* v2, int32_t n)
783  {
784  float64_t r=0;
785  for (int32_t i=0; i<n; i++)
786  r+=((float64_t) v1[i])*v2[i];
787 
788  return r;
789  }
790 
792  static inline float64_t dot(
793  const uint8_t* v1, const uint8_t* v2, int32_t n)
794  {
795  float64_t r=0;
796  for (int32_t i=0; i<n; i++)
797  r+=((float64_t) v1[i])*v2[i];
798 
799  return r;
800  }
801 
803  static inline float64_t dot(
804  const int8_t* v1, const int8_t* v2, int32_t n)
805  {
806  float64_t r=0;
807  for (int32_t i=0; i<n; i++)
808  r+=((float64_t) v1[i])*v2[i];
809 
810  return r;
811  }
812 
814  static inline float64_t dot(
815  const float64_t* v1, const char* v2, int32_t n)
816  {
817  float64_t r=0;
818  for (int32_t i=0; i<n; i++)
819  r+=((float64_t) v1[i])*v2[i];
820 
821  return r;
822  }
823 
825  template <class T>
826  static inline void vector_multiply(
827  T* target, const T* v1, const T* v2,int32_t len)
828  {
829  for (int32_t i=0; i<len; i++)
830  target[i]=v1[i]*v2[i];
831  }
832 
833 
835  template <class T>
836  static inline void add(
837  T* target, T alpha, const T* v1, T beta, const T* v2,
838  int32_t len)
839  {
840  for (int32_t i=0; i<len; i++)
841  target[i]=alpha*v1[i]+beta*v2[i];
842  }
843 
845  template <class T>
846  static inline void add_scalar(T alpha, T* vec, int32_t len)
847  {
848  for (int32_t i=0; i<len; i++)
849  vec[i]+=alpha;
850  }
851 
853  template <class T>
854  static inline void scale_vector(T alpha, T* vec, int32_t len)
855  {
856  for (int32_t i=0; i<len; i++)
857  vec[i]*=alpha;
858  }
859 
861  template <class T>
862  static inline T sum(T* vec, int32_t len)
863  {
864  T result=0;
865  for (int32_t i=0; i<len; i++)
866  result+=vec[i];
867 
868  return result;
869  }
870 
872  template <class T>
873  static inline T max(T* vec, int32_t len)
874  {
875  ASSERT(len>0);
876  T maxv=vec[0];
877 
878  for (int32_t i=1; i<len; i++)
879  maxv=CMath::max(vec[i], maxv);
880 
881  return maxv;
882  }
883 
885  template <class T>
886  static inline T sum_abs(T* vec, int32_t len)
887  {
888  T result=0;
889  for (int32_t i=0; i<len; i++)
890  result+=CMath::abs(vec[i]);
891 
892  return result;
893  }
894 
896  template <class T>
897  static inline bool fequal(T x, T y, float64_t precision=1e-6)
898  {
899  return CMath::abs(x-y)<precision;
900  }
901 
903  static inline float64_t mean(float64_t* vec, int32_t len)
904  {
906  ASSERT(vec);
907  ASSERT(len>0);
908 
909  float64_t mean=0;
910  for (int32_t i=0; i<len; i++)
911  mean+=vec[i];
912  return mean/len;
913  }
914 
915  static inline float64_t trace(
916  float64_t* mat, int32_t cols, int32_t rows)
917  {
918  float64_t trace=0;
919  for (int32_t i=0; i<rows; i++)
920  trace+=mat[i*cols+i];
921  return trace;
922  }
923 
927  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
928  static void sort(float64_t *a, int32_t*idx, int32_t N);
929 
930  /*
931  * Inline function to extract the byte at position p (from left)
932  * of an 64 bit integer. The function is somewhat identical to
933  * accessing an array of characters via [].
934  */
935 
937  template <class T>
938  inline static void radix_sort(T* array, int32_t size)
939  {
940  radix_sort_helper(array,size,0);
941  }
942 
943  template <class T>
944  static inline uint8_t byte(T word, uint16_t p)
945  {
946  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
947  }
948 
949  template <class T>
950  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
951  {
952  static size_t count[256], nc, cmin;
953  T *ak;
954  uint8_t c=0;
955  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
956  T *an, *aj, *pile[256];
957  size_t *cp, cmax;
958 
959  /* Push initial array to stack */
960  sp = s;
961  radix_push(array, size, i);
962 
963  /* Loop until all digits have been sorted */
964  while (sp>s) {
965  radix_pop(array, size, i);
966  an = array + size;
967 
968  /* Make character histogram */
969  if (nc == 0) {
970  cmin = 0xff;
971  for (ak = array; ak < an; ak++) {
972  c = byte(*ak, i);
973  count[c]++;
974  if (count[c] == 1) {
975  /* Determine smallest character */
976  if (c < cmin)
977  cmin = c;
978  nc++;
979  }
980  }
981 
982  /* Sort recursively if stack size too small */
983  if (sp + nc > s + RADIX_STACK_SIZE) {
984  radix_sort_helper(array, size, i);
985  continue;
986  }
987  }
988 
989  /*
990  * Set pile[]; push incompletely sorted bins onto stack.
991  * pile[] = pointers to last out-of-place element in bins.
992  * Before permuting: pile[c-1] + count[c] = pile[c];
993  * during deal: pile[c] counts down to pile[c-1].
994  */
995  olds = bigs = sp;
996  cmax = 2;
997  ak = array;
998  pile[0xff] = an;
999  for (cp = count + cmin; nc > 0; cp++) {
1000  /* Find next non-empty pile */
1001  while (*cp == 0)
1002  cp++;
1003  /* Pile with several entries */
1004  if (*cp > 1) {
1005  /* Determine biggest pile */
1006  if (*cp > cmax) {
1007  cmax = *cp;
1008  bigs = sp;
1009  }
1010 
1011  if (i < sizeof(T)-1)
1012  radix_push(ak, *cp, (uint16_t) (i + 1));
1013  }
1014  pile[cp - count] = ak += *cp;
1015  nc--;
1016  }
1017 
1018  /* Play it safe -- biggest bin last. */
1019  swap(*olds, *bigs);
1020 
1021  /*
1022  * Permute misplacements home. Already home: everything
1023  * before aj, and in pile[c], items from pile[c] on.
1024  * Inner loop:
1025  * r = next element to put in place;
1026  * ak = pile[r[i]] = location to put the next element.
1027  * aj = bottom of 1st disordered bin.
1028  * Outer loop:
1029  * Once the 1st disordered bin is done, ie. aj >= ak,
1030  * aj<-aj + count[c] connects the bins in array linked list;
1031  * reset count[c].
1032  */
1033  aj = array;
1034  while (aj<an)
1035  {
1036  T r;
1037 
1038  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
1039  swap(*ak, r);
1040 
1041  *aj = r;
1042  aj += count[c];
1043  count[c] = 0;
1044  }
1045  }
1046  }
1047 
1050  template <class T>
1051  static void insertion_sort(T* output, int32_t size)
1052  {
1053  for (int32_t i=0; i<size-1; i++)
1054  {
1055  int32_t j=i-1;
1056  T value=output[i];
1057  while (j >= 0 && output[j] > value)
1058  {
1059  output[j+1] = output[j];
1060  j--;
1061  }
1062  output[j+1]=value;
1063  }
1064  }
1065 
1066 
1069  template <class T>
1070  static void qsort(T* output, int32_t size)
1071  {
1072  if (size==1)
1073  return;
1074 
1075  if (size==2)
1076  {
1077  if (output[0] > output [1])
1078  swap(output[0],output[1]);
1079  return;
1080  }
1081  //T split=output[random(0,size-1)];
1082  T split=output[size/2];
1083 
1084  int32_t left=0;
1085  int32_t right=size-1;
1086 
1087  while (left<=right)
1088  {
1089  while (output[left] < split)
1090  left++;
1091  while (output[right] > split)
1092  right--;
1093 
1094  if (left<=right)
1095  {
1096  swap(output[left],output[right]);
1097  left++;
1098  right--;
1099  }
1100  }
1101 
1102  if (right+1> 1)
1103  qsort(output,right+1);
1104 
1105  if (size-left> 1)
1106  qsort(&output[left],size-left);
1107  }
1108 
1117  template <class T>
1118  static void qsort(SGVector<T*> array)
1119  {
1120  if (array.vlen==1)
1121  return;
1122 
1123  if (array.vlen==2)
1124  {
1125  if (*array.vector[0]>*array.vector[1])
1126  swap(array.vector[0],array.vector[1]);
1127  return;
1128  }
1129  T* split=array.vector[array.vlen/2];
1130 
1131  int32_t left=0;
1132  int32_t right=array.vlen-1;
1133 
1134  while (left<=right)
1135  {
1136  while (*array.vector[left]<*split)
1137  ++left;
1138  while (*array.vector[right]>*split)
1139  --right;
1140 
1141  if (left<=right)
1142  {
1143  swap(array.vector[left],array.vector[right]);
1144  ++left;
1145  --right;
1146  }
1147  }
1148 
1149  if (right+1>1)
1150  qsort(SGVector<T*>(array.vector,right+1));
1151 
1152  if (array.vlen-left>1)
1153  qsort(SGVector<T*>(&array.vector[left],array.vlen-left));
1154  }
1155 
1157  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1158  {
1159  ASSERT(width>=0);
1160  for (int i=0; i<width; i++)
1161  {
1162  T mask = ((T) 1)<<(sizeof(T)*8-1);
1163  while (mask)
1164  {
1165  if (mask & word)
1166  SG_SPRINT("1");
1167  else
1168  SG_SPRINT("0");
1169 
1170  mask>>=1;
1171  }
1172  }
1173  }
1174 
1176  template <class T> static void display_vector(
1177  const T* vector, int32_t n, const char* name="vector");
1178 
1180  template <class T> static void display_matrix(
1181  const T* matrix, int32_t rows, int32_t cols, const char* name="matrix");
1182 
1188  template <class T1,class T2>
1189  static void qsort_index(T1* output, T2* index, uint32_t size);
1190 
1196  template <class T1,class T2>
1197  static void qsort_backward_index(
1198  T1* output, T2* index, int32_t size);
1199 
1207  template <class T1,class T2>
1208  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1209  {
1210  int32_t n=0;
1211  thread_qsort<T1,T2> t;
1212  t.output=output;
1213  t.index=index;
1214  t.size=size;
1215  t.qsort_threads=&n;
1216  t.sort_limit=limit;
1217  t.num_threads=n_threads;
1218  parallel_qsort_index<T1,T2>(&t);
1219  }
1220 
1221 
1222  template <class T1,class T2>
1223  static void* parallel_qsort_index(void* p);
1224 
1225 
1226  /* finds the smallest element in output and puts that element as the
1227  first element */
1228  template <class T>
1229  static void min(float64_t* output, T* index, int32_t size);
1230 
1231  /* finds the n smallest elements in output and puts these elements as the
1232  first n elements */
1233  template <class T>
1234  static void nmin(
1235  float64_t* output, T* index, int32_t size, int32_t n);
1236 
1237  /* performs a inplace unique of a vector of type T using quicksort
1238  * returns the new number of elements */
1239  template <class T>
1240  static int32_t unique(T* output, int32_t size)
1241  {
1242  qsort(output, size);
1243  int32_t i,j=0 ;
1244  for (i=0; i<size; i++)
1245  if (i==0 || output[i]!=output[i-1])
1246  output[j++]=output[i];
1247  return j ;
1248  }
1249 
1257  static double* compute_eigenvectors(double* matrix, int n, int m)
1258  {
1259 #ifdef HAVE_LAPACK
1260  ASSERT(n == m);
1261 
1262  char V='V';
1263  char U='U';
1264  int info;
1265  int ord=n;
1266  int lda=n;
1267  double* eigenvalues=SG_CALLOC(float64_t, n+1);
1268 
1269  // lapack sym matrix eigenvalues+vectors
1270  wrap_dsyev(V, U, ord, matrix, lda,
1271  eigenvalues, &info);
1272 
1273  if (info!=0)
1274  SG_SERROR("DSYEV failed with code %d\n", info);
1275 
1276  return eigenvalues;
1277 #else
1278  SG_SERROR("Function not available - Lapack/Atlas not enabled at compile time!\n");
1279  return NULL;
1280 #endif
1281  }
1282 
1283  /* Sums up all rows of a matrix and returns the resulting rowvector */
1284  template <class T>
1285  static T* get_row_sum(T* matrix, int32_t m, int32_t n)
1286  {
1287  T* rowsums=SG_CALLOC(T, n);
1288 
1289  for (int32_t i=0; i<n; i++)
1290  {
1291  for (int32_t j=0; j<m; j++)
1292  rowsums[i]+=matrix[j+int64_t(i)*m];
1293  }
1294  return rowsums;
1295  }
1296 
1297  /* Sums up all columns of a matrix and returns the resulting columnvector */
1298  template <class T>
1299  static T* get_column_sum(T* matrix, int32_t m, int32_t n)
1300  {
1301  T* colsums=SG_CALLOC(T, m);
1302 
1303  for (int32_t i=0; i<n; i++)
1304  {
1305  for (int32_t j=0; j<m; j++)
1306  colsums[j]+=matrix[j+int64_t(i)*m];
1307  }
1308  return colsums;
1309  }
1310 
1311  /* Centers matrix (e.g. kernel matrix in feature space INPLACE */
1312  template <class T>
1313  static void center_matrix(T* matrix, int32_t m, int32_t n)
1314  {
1315  float64_t num_data=n;
1316 
1317  T* colsums=get_column_sum(matrix, m,n);
1318  T* rowsums=get_row_sum(matrix, m,n);
1319 
1320  for (int32_t i=0; i<m; i++)
1321  colsums[i]/=num_data;
1322  for (int32_t j=0; j<n; j++)
1323  rowsums[j]/=num_data;
1324 
1325  T s=sum(rowsums, n)/num_data;
1326 
1327  for (int32_t i=0; i<n; i++)
1328  {
1329  for (int32_t j=0; j<m; j++)
1330  matrix[int64_t(i)*m+j]+=s-colsums[j]-rowsums[i];
1331  }
1332 
1333  SG_FREE(rowsums);
1334  SG_FREE(colsums);
1335  }
1336 
1337  /* finds an element in a sorted array via binary search
1338  * returns -1 if not found */
1339  template <class T>
1340  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1341  {
1342  int32_t start=0;
1343  int32_t end=size-1;
1344 
1345  if (size<1)
1346  return 0;
1347 
1348  while (start<end)
1349  {
1350  int32_t middle=(start+end)/2;
1351 
1352  if (output[middle]>elem)
1353  end=middle-1;
1354  else if (output[middle]<elem)
1355  start=middle+1;
1356  else
1357  return middle;
1358  }
1359 
1360  return start;
1361  }
1362 
1363  /* finds an element in a sorted array via binary search
1364  * * returns -1 if not found */
1365  template <class T>
1366  static inline int32_t binary_search(T* output, int32_t size, T elem)
1367  {
1368  int32_t ind = binary_search_helper(output, size, elem);
1369  if (ind >= 0 && output[ind] == elem)
1370  return ind;
1371  return -1;
1372  }
1373 
1374  /* Finds an element in a sorted array of pointers via binary search
1375  * Every element is dereferenced once before being compared
1376  *
1377  * @param array array of pointers to search in (assumed being sorted)
1378  * @param elem pointer to element to search for
1379  * @return index of elem, -1 if not found */
1380  template<class T>
1381  static inline int32_t binary_search(SGVector<T*> array, T* elem)
1382  {
1383  int32_t start=0;
1384  int32_t end=array.vlen-1;
1385 
1386  if (array.vlen<1)
1387  return -1;
1388 
1389  while (start<end)
1390  {
1391  int32_t middle=(start+end)/2;
1392 
1393  if (*array.vector[middle]>*elem)
1394  end=middle-1;
1395  else if (*array.vector[middle]<*elem)
1396  start=middle+1;
1397  else
1398  {
1399  start=middle;
1400  break;
1401  }
1402  }
1403 
1404  if (start>=0&&*array.vector[start]==*elem)
1405  return start;
1406 
1407  return -1;
1408  }
1409 
1410  template <class T>
1412  T* output, int32_t size, T elem)
1413  {
1414  int32_t ind = binary_search_helper(output, size, elem);
1415 
1416  if (output[ind]<=elem)
1417  return ind;
1418  if (ind>0 && output[ind-1] <= elem)
1419  return ind-1;
1420  return -1;
1421  }
1422 
1425  static float64_t Align(
1426  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1427 
1428 
1430 
1433  static float64_t mutual_info(float64_t* p1, float64_t* p2, int32_t len);
1434 
1437  static float64_t relative_entropy(
1438  float64_t* p, float64_t* q, int32_t len);
1439 
1441  static float64_t entropy(float64_t* p, int32_t len);
1442 
1444  inline static uint32_t get_seed()
1445  {
1446  return CMath::seed;
1447  }
1448 
1450  inline static uint32_t get_log_range()
1451  {
1452  return CMath::LOGRANGE;
1453  }
1454 
1455 #ifdef USE_LOGCACHE
1456 
1457  inline static uint32_t get_log_accuracy()
1458  {
1459  return CMath::LOGACCURACY;
1460  }
1461 #endif
1462 
1464  inline static int is_finite(double f)
1465  {
1466 #if defined(isfinite) && !defined(SUNOS)
1467  return isfinite(f);
1468 #else
1469  return finite(f);
1470 #endif
1471  }
1472 
1474  inline static int is_infinity(double f)
1475  {
1476 #ifdef SUNOS
1477  if (fpclass(f) == FP_NINF || fpclass(f) == FP_PINF)
1478  return 1;
1479  else
1480  return 0;
1481 #else
1482  return isinf(f);
1483 #endif
1484  }
1485 
1487  inline static int is_nan(double f)
1488  {
1489 #ifdef SUNOS
1490  return isnand(f);
1491 #else
1492  return isnan(f);
1493 #endif
1494  }
1495 
1500 
1505 
1516 #ifdef USE_LOGCACHE
1517  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1518  {
1519  float64_t diff;
1520 
1521  if (!CMath::is_finite(p))
1522  return q;
1523 
1524  if (!CMath::is_finite(q))
1525  {
1526  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q);
1527  return NAN;
1528  }
1529  diff = p - q;
1530  if (diff > 0)
1531  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1532  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1533  }
1534 
1536  static void init_log_table();
1537 
1539  static int32_t determine_logrange();
1540 
1542  static int32_t determine_logaccuracy(int32_t range);
1543 #else
1545  float64_t p, float64_t q)
1546  {
1547  float64_t diff;
1548 
1549  if (!CMath::is_finite(p))
1550  return q;
1551  if (!CMath::is_finite(q))
1552  return p;
1553  diff = p - q;
1554  if (diff > 0)
1555  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
1556  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
1557  }
1558 #endif
1559 #ifdef USE_LOGSUMARRAY
1560 
1565  static inline float64_t logarithmic_sum_array(
1566  float64_t *p, int32_t len)
1567  {
1568  if (len<=2)
1569  {
1570  if (len==2)
1571  return logarithmic_sum(p[0],p[1]) ;
1572  if (len==1)
1573  return p[0];
1574  return -INFTY ;
1575  }
1576  else
1577  {
1578  register float64_t *pp=p ;
1579  if (len%2==1) pp++ ;
1580  for (register int32_t j=0; j < len>>1; j++)
1581  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
1582  }
1583  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1584  }
1585 #endif //USE_LOGSUMARRAY
1586 
1587 
1589  inline virtual const char* get_name() const { return "Mathematics"; }
1590  public:
1593 
1594  static const float64_t INFTY;
1595  static const float64_t ALMOST_INFTY;
1596 
1599 
1601  static const float64_t PI;
1602 
1605 
1606  /* largest and smallest possible float64_t */
1609 
1610  protected:
1612  static int32_t LOGRANGE;
1613 
1615  static uint32_t seed;
1616  static char* rand_state;
1617 
1618 #ifdef USE_LOGCACHE
1619 
1621  static int32_t LOGACCURACY;
1623 
1624  static float64_t* logtable;
1625 #endif
1626 };
1627 
1628 //implementations of template functions
1629 template <class T1,class T2>
1631  {
1632  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1633  T1* output=ps->output;
1634  T2* index=ps->index;
1635  uint32_t size=ps->size;
1636  int32_t* qsort_threads=ps->qsort_threads;
1637  int32_t sort_limit=ps->sort_limit;
1638  int32_t num_threads=ps->num_threads;
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 NULL;
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  bool lthread_start=false;
1671  bool rthread_start=false;
1672  pthread_t lthread;
1673  pthread_t rthread;
1674  struct thread_qsort<T1,T2> t1;
1675  struct thread_qsort<T1,T2> t2;
1676 
1677  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1678  qsort_index(output,index,right+1);
1679  else if (right+1> 1)
1680  {
1681  (*qsort_threads)++;
1682  lthread_start=true;
1683  t1.output=output;
1684  t1.index=index;
1685  t1.size=right+1;
1686  t1.qsort_threads=qsort_threads;
1687  t1.sort_limit=sort_limit;
1688  t1.num_threads=num_threads;
1689  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1690  {
1691  lthread_start=false;
1692  (*qsort_threads)--;
1693  qsort_index(output,index,right+1);
1694  }
1695  }
1696 
1697 
1698  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1699  qsort_index(&output[left],&index[left], size-left);
1700  else if (size-left> 1)
1701  {
1702  (*qsort_threads)++;
1703  rthread_start=true;
1704  t2.output=&output[left];
1705  t2.index=&index[left];
1706  t2.size=size-left;
1707  t2.qsort_threads=qsort_threads;
1708  t2.sort_limit=sort_limit;
1709  t2.num_threads=num_threads;
1710  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1711  {
1712  rthread_start=false;
1713  (*qsort_threads)--;
1714  qsort_index(&output[left],&index[left], size-left);
1715  }
1716  }
1717 
1718  if (lthread_start)
1719  {
1720  pthread_join(lthread, NULL);
1721  (*qsort_threads)--;
1722  }
1723 
1724  if (rthread_start)
1725  {
1726  pthread_join(rthread, NULL);
1727  (*qsort_threads)--;
1728  }
1729 
1730  return NULL;
1731  }
1732 
1733  template <class T1,class T2>
1734 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
1735 {
1736  if (size==1)
1737  return;
1738 
1739  if (size==2)
1740  {
1741  if (output[0] > output [1])
1742  {
1743  swap(output[0],output[1]);
1744  swap(index[0],index[1]);
1745  }
1746  return;
1747  }
1748  //T1 split=output[random(0,size-1)];
1749  T1 split=output[size/2];
1750 
1751  int32_t left=0;
1752  int32_t right=size-1;
1753 
1754  while (left<=right)
1755  {
1756  while (output[left] < split)
1757  left++;
1758  while (output[right] > split)
1759  right--;
1760 
1761  if (left<=right)
1762  {
1763  swap(output[left],output[right]);
1764  swap(index[left],index[right]);
1765  left++;
1766  right--;
1767  }
1768  }
1769 
1770  if (right+1> 1)
1771  qsort_index(output,index,right+1);
1772 
1773  if (size-left> 1)
1774  qsort_index(&output[left],&index[left], size-left);
1775 }
1776 
1777  template <class T1,class T2>
1778 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
1779 {
1780  if (size==2)
1781  {
1782  if (output[0] < output [1])
1783  {
1784  swap(output[0],output[1]);
1785  swap(index[0],index[1]);
1786  }
1787  return;
1788  }
1789 
1790  //T1 split=output[random(0,size-1)];
1791  T1 split=output[size/2];
1792 
1793  int32_t left=0;
1794  int32_t right=size-1;
1795 
1796  while (left<=right)
1797  {
1798  while (output[left] > split)
1799  left++;
1800  while (output[right] < split)
1801  right--;
1802 
1803  if (left<=right)
1804  {
1805  swap(output[left],output[right]);
1806  swap(index[left],index[right]);
1807  left++;
1808  right--;
1809  }
1810  }
1811 
1812  if (right+1> 1)
1813  qsort_backward_index(output,index,right+1);
1814 
1815  if (size-left> 1)
1816  qsort_backward_index(&output[left],&index[left], size-left);
1817 }
1818 
1819  template <class T>
1820 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
1821 {
1822  if (6*n*size<13*size*CMath::log(size))
1823  for (int32_t i=0; i<n; i++)
1824  min(&output[i], &index[i], size-i) ;
1825  else
1826  qsort_index(output, index, size) ;
1827 }
1828 
1829 /* move the smallest entry in the array to the beginning */
1830  template <class T>
1831 void CMath::min(float64_t* output, T* index, int32_t size)
1832 {
1833  if (size<=1)
1834  return;
1835  float64_t min_elem=output[0];
1836  int32_t min_index=0;
1837  for (int32_t i=1; i<size; i++)
1838  {
1839  if (output[i]<min_elem)
1840  {
1841  min_index=i;
1842  min_elem=output[i];
1843  }
1844  }
1845  swap(output[0], output[min_index]);
1846  swap(index[0], index[min_index]);
1847 }
1848 }
1849 #endif

SHOGUN Machine Learning Toolbox - Documentation