SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Math.cpp
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) 1999-2009 Soeren Sonnenburg
8  * Written (W) 1999-2008 Gunnar Raetsch
9  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
10  */
11 #include <shogun/base/SGObject.h>
12 #include <shogun/lib/common.h>
15 #include <shogun/io/SGIO.h>
16 
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <math.h>
20 
21 using namespace shogun;
22 
23 #ifdef USE_LOGCACHE
24 #ifdef USE_HMMDEBUG
25 #define MAX_LOG_TABLE_SIZE 10*1024*1024
26 #define LOG_TABLE_PRECISION 1e-6
27 #else //USE_HMMDEBUG
28 #define MAX_LOG_TABLE_SIZE 123*1024*1024
29 #define LOG_TABLE_PRECISION 1e-15
30 #endif //USE_HMMDEBUG
31 int32_t CMath::LOGACCURACY = 0; // 100000 steps per integer
32 #endif // USE_LOGCACHE
33 
34 int32_t CMath::LOGRANGE = 0; // range for logtable: log(1+exp(x)) -25 <= x <= 0
35 
36 const float64_t CMath::INFTY = -log(0.0); // infinity
37 const float64_t CMath::ALMOST_INFTY = +1e+20; //a large number
39 const float64_t CMath::PI=M_PI;
43 
44 #ifdef USE_LOGCACHE
45 float64_t* CMath::logtable = NULL;
46 #endif
47 char* CMath::rand_state = NULL;
48 uint32_t CMath::seed = 0;
49 
51 : CSGObject()
52 {
54  init_random();
55 
56 #ifdef USE_LOGCACHE
57  LOGRANGE=CMath::determine_logrange();
58  LOGACCURACY=CMath::determine_logaccuracy(LOGRANGE);
59  CMath::logtable=SG_MALLOC(float64_t, LOGRANGE*LOGACCURACY);
60  init_log_table();
61 #else
62  int32_t i=0;
63  while ((float64_t)log(1+((float64_t)exp(-float64_t(i)))))
64  i++;
65 
66  LOGRANGE=i;
67 #endif
68 }
69 
71 {
73  CMath::rand_state=NULL;
74 #ifdef USE_LOGCACHE
75  SG_FREE(CMath::logtable);
76  CMath::logtable=NULL;
77 #endif
78 }
79 
80 #ifdef USE_LOGCACHE
81 int32_t CMath::determine_logrange()
82 {
83  int32_t i;
84  float64_t acc=0;
85  for (i=0; i<50; i++)
86  {
87  acc=((float64_t)log(1+((float64_t)exp(-float64_t(i)))));
88  if (acc<=(float64_t)LOG_TABLE_PRECISION)
89  break;
90  }
91 
92  SG_SINFO( "determined range for x in table log(1+exp(-x)) is:%d (error:%G)\n",i,acc);
93  return i;
94 }
95 
96 int32_t CMath::determine_logaccuracy(int32_t range)
97 {
98  range=MAX_LOG_TABLE_SIZE/range/((int)sizeof(float64_t));
99  SG_SINFO( "determined accuracy for x in table log(1+exp(-x)) is:%d (error:%G)\n",range,1.0/(double) range);
100  return range;
101 }
102 
103 //init log table of form log(1+exp(x))
104 void CMath::init_log_table()
105 {
106  for (int32_t i=0; i< LOGACCURACY*LOGRANGE; i++)
107  logtable[i]=log(1+exp(float64_t(-i)/float64_t(LOGACCURACY)));
108 }
109 #endif
110 
111 void CMath::sort(int32_t *a, int32_t cols, int32_t sort_col)
112 {
113  int32_t changed=1;
114  if (a[0]==-1) return ;
115  while (changed)
116  {
117  changed=0; int32_t i=0 ;
118  while ((a[(i+1)*cols]!=-1) && (a[(i+1)*cols+1]!=-1)) // to be sure
119  {
120  if (a[i*cols+sort_col]>a[(i+1)*cols+sort_col])
121  {
122  for (int32_t j=0; j<cols; j++)
123  CMath::swap(a[i*cols+j],a[(i+1)*cols+j]) ;
124  changed=1 ;
125  } ;
126  i++ ;
127  } ;
128  } ;
129 }
130 
131 void CMath::sort(float64_t *a, int32_t* idx, int32_t N)
132 {
133  int32_t changed=1;
134  while (changed)
135  {
136  changed=0;
137  for (int32_t i=0; i<N-1; i++)
138  {
139  if (a[i]>a[i+1])
140  {
141  swap(a[i],a[i+1]) ;
142  swap(idx[i],idx[i+1]) ;
143  changed=1 ;
144  } ;
145  } ;
146  } ;
147 
148 }
149 
151  char* seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost)
152 {
153  float64_t actCost=0 ;
154  int32_t i1, i2 ;
155  float64_t* const gapCosts1 = SG_MALLOC(float64_t, l1 );
156  float64_t* const gapCosts2 = SG_MALLOC(float64_t, l2 );
157  float64_t* costs2_0 = SG_MALLOC(float64_t, l2 + 1 );
158  float64_t* costs2_1 = SG_MALLOC(float64_t, l2 + 1 );
159 
160  // initialize borders
161  for( i1 = 0; i1 < l1; ++i1 ) {
162  gapCosts1[ i1 ] = gapCost * i1;
163  }
164  costs2_1[ 0 ] = 0;
165  for( i2 = 0; i2 < l2; ++i2 ) {
166  gapCosts2[ i2 ] = gapCost * i2;
167  costs2_1[ i2+1 ] = costs2_1[ i2 ] + gapCosts2[ i2 ];
168  }
169  // compute alignment
170  for( i1 = 0; i1 < l1; ++i1 ) {
171  swap( costs2_0, costs2_1 );
172  actCost = costs2_0[ 0 ] + gapCosts1[ i1 ];
173  costs2_1[ 0 ] = actCost;
174  for( i2 = 0; i2 < l2; ++i2 ) {
175  const float64_t actMatch = costs2_0[ i2 ] + ( seq1[i1] == seq2[i2] );
176  const float64_t actGap1 = costs2_0[ i2+1 ] + gapCosts1[ i1 ];
177  const float64_t actGap2 = actCost + gapCosts2[ i2 ];
178  const float64_t actGap = min( actGap1, actGap2 );
179  actCost = min( actMatch, actGap );
180  costs2_1[ i2+1 ] = actCost;
181  }
182  }
183 
184  SG_FREE(gapCosts1);
185  SG_FREE(gapCosts2);
186  SG_FREE(costs2_0);
187  SG_FREE(costs2_1);
188 
189  // return the final cost
190  return actCost;
191 }

SHOGUN Machine Learning Toolbox - Documentation