SHOGUN  v3.2.0
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/lib/config.h>
12 #include <shogun/base/SGObject.h>
13 #include <shogun/lib/common.h>
14 #include <cmath>
17 #include <shogun/io/SGIO.h>
18 #include <shogun/lib/SGVector.h>
19 
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <math.h>
23 
24 #ifndef NAN
25 #include <stdlib.h>
26 #define NAN (strtod("NAN",NULL))
27 #endif
28 
29 
30 using namespace shogun;
31 
32 #ifdef USE_LOGCACHE
33 #ifdef USE_HMMDEBUG
34 #define MAX_LOG_TABLE_SIZE 10*1024*1024
35 #define LOG_TABLE_PRECISION 1e-6
36 #else //USE_HMMDEBUG
37 #define MAX_LOG_TABLE_SIZE 123*1024*1024
38 #define LOG_TABLE_PRECISION 1e-15
39 #endif //USE_HMMDEBUG
40 int32_t CMath::LOGACCURACY = 0; // 100000 steps per integer
41 #endif // USE_LOGCACHE
42 
43 int32_t CMath::LOGRANGE = 0; // range for logtable: log(1+exp(x)) -25 <= x <= 0
44 
46 const float64_t CMath::INFTY = INFINITY; // infinity
47 const float64_t CMath::ALMOST_INFTY = +1e+20; //a large number
53 
54 #ifdef USE_LOGCACHE
55 float64_t* CMath::logtable = NULL;
56 #endif
57 uint32_t CMath::seed = 0;
58 
60 : CSGObject()
61 {
62 #ifdef USE_LOGCACHE
63  LOGRANGE=CMath::determine_logrange();
64  LOGACCURACY=CMath::determine_logaccuracy(LOGRANGE);
65  CMath::logtable=SG_MALLOC(float64_t, LOGRANGE*LOGACCURACY);
66  init_log_table();
67 #else
68  int32_t i=0;
69  while ((float64_t)log(1+((float64_t)exp(-float64_t(i)))))
70  i++;
71 
72  LOGRANGE=i;
73 #endif
74 }
75 
77 {
78 #ifdef USE_LOGCACHE
79  SG_FREE(CMath::logtable);
80  CMath::logtable=NULL;
81 #endif
82 }
83 
84 #ifdef USE_LOGCACHE
85 int32_t CMath::determine_logrange()
86 {
87  int32_t i;
88  float64_t acc=0;
89  for (i=0; i<50; i++)
90  {
91  acc=((float64_t)log(1+((float64_t)exp(-float64_t(i)))));
92  if (acc<=(float64_t)LOG_TABLE_PRECISION)
93  break;
94  }
95 
96  SG_SINFO("determined range for x in table log(1+exp(-x)) is:%d (error:%G)\n",i,acc)
97  return i;
98 }
99 
100 int32_t CMath::determine_logaccuracy(int32_t range)
101 {
102  range=MAX_LOG_TABLE_SIZE/range/((int)sizeof(float64_t));
103  SG_SINFO("determined accuracy for x in table log(1+exp(-x)) is:%d (error:%G)\n",range,1.0/(double) range)
104  return range;
105 }
106 
107 //init log table of form log(1+exp(x))
108 void CMath::init_log_table()
109 {
110  for (int32_t i=0; i< LOGACCURACY*LOGRANGE; i++)
111  logtable[i]=log(1+exp(float64_t(-i)/float64_t(LOGACCURACY)));
112 }
113 #endif
114 
115 void CMath::sort(int32_t *a, int32_t cols, int32_t sort_col)
116 {
117  int32_t changed=1;
118  if (a[0]==-1) return ;
119  while (changed)
120  {
121  changed=0; int32_t i=0 ;
122  while ((a[(i+1)*cols]!=-1) && (a[(i+1)*cols+1]!=-1)) // to be sure
123  {
124  if (a[i*cols+sort_col]>a[(i+1)*cols+sort_col])
125  {
126  for (int32_t j=0; j<cols; j++)
127  CMath::swap(a[i*cols+j],a[(i+1)*cols+j]) ;
128  changed=1 ;
129  } ;
130  i++ ;
131  } ;
132  } ;
133 }
134 
135 void CMath::sort(float64_t *a, int32_t* idx, int32_t N)
136 {
137  int32_t changed=1;
138  while (changed)
139  {
140  changed=0;
141  for (int32_t i=0; i<N-1; i++)
142  {
143  if (a[i]>a[i+1])
144  {
145  swap(a[i],a[i+1]) ;
146  swap(idx[i],idx[i+1]) ;
147  changed=1 ;
148  } ;
149  } ;
150  } ;
151 
152 }
153 
155  char* seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost)
156 {
157  float64_t actCost=0 ;
158  int32_t i1, i2 ;
159  float64_t* const gapCosts1 = SG_MALLOC(float64_t, l1 );
160  float64_t* const gapCosts2 = SG_MALLOC(float64_t, l2 );
161  float64_t* costs2_0 = SG_MALLOC(float64_t, l2 + 1 );
162  float64_t* costs2_1 = SG_MALLOC(float64_t, l2 + 1 );
163 
164  // initialize borders
165  for( i1 = 0; i1 < l1; ++i1 ) {
166  gapCosts1[ i1 ] = gapCost * i1;
167  }
168  costs2_1[ 0 ] = 0;
169  for( i2 = 0; i2 < l2; ++i2 ) {
170  gapCosts2[ i2 ] = gapCost * i2;
171  costs2_1[ i2+1 ] = costs2_1[ i2 ] + gapCosts2[ i2 ];
172  }
173  // compute alignment
174  for( i1 = 0; i1 < l1; ++i1 ) {
175  swap( costs2_0, costs2_1 );
176  actCost = costs2_0[ 0 ] + gapCosts1[ i1 ];
177  costs2_1[ 0 ] = actCost;
178  for( i2 = 0; i2 < l2; ++i2 ) {
179  const float64_t actMatch = costs2_0[ i2 ] + ( seq1[i1] == seq2[i2] );
180  const float64_t actGap1 = costs2_0[ i2+1 ] + gapCosts1[ i1 ];
181  const float64_t actGap2 = actCost + gapCosts2[ i2 ];
182  const float64_t actGap = min( actGap1, actGap2 );
183  actCost = min( actMatch, actGap );
184  costs2_1[ i2+1 ] = actCost;
185  }
186  }
187 
188  SG_FREE(gapCosts1);
189  SG_FREE(gapCosts2);
190  SG_FREE(costs2_0);
191  SG_FREE(costs2_1);
192 
193  // return the final cost
194  return actCost;
195 }
196 
197 void CMath::linspace(float64_t* output, float64_t start, float64_t end, int32_t n)
198 {
199  float64_t delta = (end-start) / (n-1);
200  float64_t v = start;
201  index_t i = 0;
202  while ( v <= end )
203  {
204  output[i++] = v;
205  v += delta;
206  }
207  output[n-1] = end;
208 }
209 
210 
211 int CMath::is_nan(double f)
212 {
213 #ifndef HAVE_STD_ISNAN
214 #if (HAVE_DECL_ISNAN == 1) || defined(HAVE_ISNAN)
215  return ::isnan(f);
216 #else
217  return ((f != f) ? 1 : 0);
218 #endif // #if (HAVE_DECL_ISNAN == 1) || defined(HAVE_ISNAN)
219 #endif // #ifndef HAVE_STD_ISNAN
220 
221  return std::isnan(f);
222 }
223 
224 int CMath::is_infinity(double f)
225 {
226 #ifndef HAVE_STD_ISINF
227 #if (HAVE_DECL_ISINF == 1) || defined(HAVE_ISINF)
228  return ::isinf(f);
229 #elif defined(FPCLASS)
230  if (::fpclass(f) == FP_NINF) return -1;
231  else if (::fpclass(f) == FP_PINF) return 1;
232  else return 0;
233 #else
234  if ((f == f) && ((f - f) != 0.0)) return (f < 0.0 ? -1 : 1);
235  else return 0;
236 }
237 #endif // #if (HAVE_DECL_ISINF == 1) || defined(HAVE_ISINF)
238 #endif // #ifndef HAVE_STD_ISINF
239 
240  return std::isinf(f);
241 }
242 
243 int CMath::is_finite(double f)
244 {
245 #ifndef HAVE_STD_ISFINITE
246 #if (HAVE_DECL_ISFINITE == 1) || defined(HAVE_ISFINITE)
247  return ::isfinite(f);
248 #elif defined(HAVE_FINITE)
249  return ::finite(f);
250 #else
251  return ((!std::isnan(f) && !std::isinf(f)) ? 1 : 0);
252 #endif // #if (HAVE_DECL_ISFINITE == 1) || defined(HAVE_ISFINITE)
253 #endif // #ifndef HAVE_STD_ISFINITE
254 
255  return std::isfinite(f);
256 }
257 
258 bool CMath::strtof(const char* str, float32_t* float_result)
259 {
260  ASSERT(str);
261  ASSERT(float_result);
262 
263  SGVector<char> buf(strlen(str)+1);
264 
265  for (index_t i=0; i<buf.vlen-1; i++)
266  buf[i]=tolower(str[i]);
267  buf[buf.vlen-1]='\0';
268 
269  if (strstr(buf, "inf") != NULL)
270  {
271  *float_result = CMath::INFTY;
272 
273  if (strchr(buf,'-') != NULL)
274  *float_result *= -1;
275  return true;
276  }
277 
278  if (strstr(buf, "nan") != NULL)
279  {
280  *float_result = CMath::NOT_A_NUMBER;
281  return true;
282  }
283 
284  char* endptr = buf.vector;
285  *float_result=::strtof(str, &endptr);
286  return endptr != buf.vector;
287 }
288 
289 bool CMath::strtod(const char* str, float64_t* double_result)
290 {
291  ASSERT(str);
292  ASSERT(double_result);
293 
294  SGVector<char> buf(strlen(str)+1);
295 
296  for (index_t i=0; i<buf.vlen-1; i++)
297  buf[i]=tolower(str[i]);
298  buf[buf.vlen-1]='\0';
299 
300  if (strstr(buf, "inf") != NULL)
301  {
302  *double_result = CMath::INFTY;
303 
304  if (strchr(buf,'-') != NULL)
305  *double_result *= -1;
306  return true;
307  }
308 
309  if (strstr(buf, "nan") != NULL)
310  {
311  *double_result = CMath::NOT_A_NUMBER;
312  return true;
313  }
314 
315  char* endptr = buf.vector;
316  *double_result=::strtod(str, &endptr);
317  return endptr != buf.vector;
318 }
319 
320 bool CMath::strtold(const char* str, floatmax_t* long_double_result)
321 {
322  ASSERT(str);
323  ASSERT(long_double_result);
324 
325  SGVector<char> buf(strlen(str)+1);
326 
327  for (index_t i=0; i<buf.vlen-1; i++)
328  buf[i]=tolower(str[i]);
329  buf[buf.vlen-1]='\0';
330 
331  if (strstr(buf, "inf") != NULL)
332  {
333  *long_double_result = CMath::INFTY;
334 
335  if (strchr(buf,'-') != NULL)
336  *long_double_result *= -1;
337  return true;
338  }
339 
340  if (strstr(buf, "nan") != NULL)
341  {
342  *long_double_result = CMath::NOT_A_NUMBER;
343  return true;
344  }
345 
346  char* endptr = buf.vector;
347 
348 // fall back to double on win32 / cygwin since strtold is undefined there
349 #if defined(WIN32) || defined(__CYGWIN__)
350  *long_double_result=::strtod(str, &endptr);
351 #else
352  *long_double_result=::strtold(str, &endptr);
353 #endif
354 
355  return endptr != buf.vector;
356 }
357 
static const float64_t MACHINE_EPSILON
Definition: Math.h:1340
static bool strtof(const char *str, float32_t *float_result)
Definition: Math.cpp:258
static uint32_t seed
random generator seed
Definition: Math.h:1351
static int is_finite(double f)
checks whether a float is finite
Definition: Math.cpp:243
static float64_t Align(char *seq1, char *seq2, int32_t l1, int32_t l2, float64_t gapCost)
Definition: Math.cpp:154
static void linspace(float64_t *output, float64_t start, float64_t end, int32_t n=100)
Definition: Math.cpp:197
int32_t index_t
Definition: common.h:60
static bool strtod(const char *str, float64_t *double_result)
Definition: Math.cpp:289
virtual ~CMath()
Destructor - frees logtable.
Definition: Math.cpp:76
static const float64_t INFTY
infinity
Definition: Math.h:1330
static const float64_t MIN_REAL_NUMBER
Definition: Math.h:1344
static int32_t LOGRANGE
range for logtable: log(1+exp(x)) -LOGRANGE <= x <= 0
Definition: Math.h:1348
static const float64_t ALMOST_NEG_INFTY
almost neg (log) infinity
Definition: Math.h:1334
CMath()
Constructor - initializes log-table.
Definition: Math.cpp:59
#define ASSERT(x)
Definition: SGIO.h:203
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:102
double float64_t
Definition: common.h:48
long double floatmax_t
Definition: common.h:49
#define M_PI
workaround for log2 being a define on cygwin
Definition: Math.h:56
float float32_t
Definition: common.h:47
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:16
static int is_infinity(double f)
checks whether a float is infinity
Definition: Math.cpp:224
static int is_nan(double f)
checks whether a float is nan
Definition: Math.cpp:211
static T min(T a, T b)
return the minimum of two integers
Definition: Math.h:153
static float64_t exp(float64_t x)
Definition: Math.h:359
#define SG_SINFO(...)
Definition: SGIO.h:175
static float64_t log(float64_t v)
Definition: Math.h:420
static const float64_t ALMOST_INFTY
Definition: Math.h:1331
Class which collects generic mathematical functions.
Definition: Math.h:133
static void swap(T &a, T &b)
swap e.g. floats a and b
Definition: Math.h:229
static void sort(int32_t *a, int32_t cols, int32_t sort_col=0)
Definition: Math.cpp:115
#define NAN
Definition: Math.cpp:26
#define delta
Definition: sfa.cpp:23
static bool strtold(const char *str, floatmax_t *long_double_result)
Definition: Math.cpp:320
static const float64_t NOT_A_NUMBER
not a number
Definition: Math.h:1328
static const float64_t MAX_REAL_NUMBER
Definition: Math.h:1343
index_t vlen
Definition: SGVector.h:706
static const float64_t PI
Definition: Math.h:1337

SHOGUN Machine Learning Toolbox - Documentation