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

SHOGUN Machine Learning Toolbox - Documentation