SHOGUN  5.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
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>
20 
21 #include <stdlib.h>
22 #include <float.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+300; //a large number
48 const float64_t CMath::ALMOST_NEG_INFTY = -1e+300;
50 const float64_t CMath::MACHINE_EPSILON=DBL_EPSILON;
51 const float64_t CMath::MAX_REAL_NUMBER=DBL_MAX;
52 const float64_t CMath::MIN_REAL_NUMBER=DBL_MIN;
53 const float32_t CMath::F_MAX_VAL32=FLT_MAX;
55 const float64_t CMath::F_MAX_VAL64=DBL_MAX;
57 const float32_t CMath::F_MIN_VAL32=(FLT_MIN * FLT_EPSILON);
58 const float64_t CMath::F_MIN_VAL64=(DBL_MIN * DBL_EPSILON);
59 
60 #ifdef USE_LOGCACHE
61 float64_t* CMath::logtable = NULL;
62 #endif
63 uint32_t CMath::seed = 0;
64 
66 : CSGObject()
67 {
68 #ifdef USE_LOGCACHE
69  LOGRANGE=CMath::determine_logrange();
70  LOGACCURACY=CMath::determine_logaccuracy(LOGRANGE);
71  CMath::logtable=SG_MALLOC(float64_t, LOGRANGE*LOGACCURACY);
72  init_log_table();
73 #else
74  int32_t i=0;
75  while ((float64_t)log(1+((float64_t)exp(-float64_t(i)))))
76  i++;
77 
78  LOGRANGE=i;
79 #endif
80 }
81 
83 {
84 #ifdef USE_LOGCACHE
85  SG_FREE(CMath::logtable);
86  CMath::logtable=NULL;
87 #endif
88 }
89 
90 float64_t CMath::dot(const float64_t* v1, const float64_t* v2, int32_t n)
91 {
92  float64_t r=0;
95  r = ev1.dot(ev2);
96  return r;
97 }
98 
99 float32_t CMath::dot(const float32_t* v1, const float32_t* v2, int32_t n)
100 {
101  float32_t r=0;
104  r = ev1.dot(ev2);
105  return r;
106 }
107 
108 #ifdef USE_LOGCACHE
109 int32_t CMath::determine_logrange()
110 {
111  int32_t i;
112  float64_t acc=0;
113  for (i=0; i<50; i++)
114  {
115  acc=((float64_t)log(1+((float64_t)exp(-float64_t(i)))));
116  if (acc<=(float64_t)LOG_TABLE_PRECISION)
117  break;
118  }
119 
120  SG_SINFO("determined range for x in table log(1+exp(-x)) is:%d (error:%G)\n",i,acc)
121  return i;
122 }
123 
124 int32_t CMath::determine_logaccuracy(int32_t range)
125 {
126  range=MAX_LOG_TABLE_SIZE/range/((int)sizeof(float64_t));
127  SG_SINFO("determined accuracy for x in table log(1+exp(-x)) is:%d (error:%G)\n",range,1.0/(double) range)
128  return range;
129 }
130 
131 //init log table of form log(1+exp(x))
132 void CMath::init_log_table()
133 {
134  for (int32_t i=0; i< LOGACCURACY*LOGRANGE; i++)
135  logtable[i]=log(1+exp(float64_t(-i)/float64_t(LOGACCURACY)));
136 }
137 #endif
138 
139 void CMath::sort(int32_t *a, int32_t cols, int32_t sort_col)
140 {
141  int32_t changed=1;
142  if (a[0]==-1) return;
143  while (changed)
144  {
145  changed=0; int32_t i=0;
146  while ((a[(i+1)*cols]!=-1) && (a[(i+1)*cols+1]!=-1)) // to be sure
147  {
148  if (a[i*cols+sort_col]>a[(i+1)*cols+sort_col])
149  {
150  for (int32_t j=0; j<cols; j++)
151  CMath::swap(a[i*cols+j],a[(i+1)*cols+j]);
152  changed=1;
153  };
154  i++;
155  };
156  };
157 }
158 
159 void CMath::sort(float64_t *a, int32_t* idx, int32_t N)
160 {
161  int32_t changed=1;
162  while (changed)
163  {
164  changed=0;
165  for (int32_t i=0; i<N-1; i++)
166  {
167  if (a[i]>a[i+1])
168  {
169  swap(a[i],a[i+1]) ;
170  swap(idx[i],idx[i+1]) ;
171  changed=1 ;
172  } ;
173  } ;
174  } ;
175 
176 }
177 
179  char* seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost)
180 {
181  float64_t actCost=0 ;
182  int32_t i1, i2 ;
183  float64_t* const gapCosts1 = SG_MALLOC(float64_t, l1 );
184  float64_t* const gapCosts2 = SG_MALLOC(float64_t, l2 );
185  float64_t* costs2_0 = SG_MALLOC(float64_t, l2 + 1 );
186  float64_t* costs2_1 = SG_MALLOC(float64_t, l2 + 1 );
187 
188  // initialize borders
189  for( i1 = 0; i1 < l1; ++i1 ) {
190  gapCosts1[ i1 ] = gapCost * i1;
191  }
192  costs2_1[ 0 ] = 0;
193  for( i2 = 0; i2 < l2; ++i2 ) {
194  gapCosts2[ i2 ] = gapCost * i2;
195  costs2_1[ i2+1 ] = costs2_1[ i2 ] + gapCosts2[ i2 ];
196  }
197  // compute alignment
198  for( i1 = 0; i1 < l1; ++i1 ) {
199  swap( costs2_0, costs2_1 );
200  actCost = costs2_0[ 0 ] + gapCosts1[ i1 ];
201  costs2_1[ 0 ] = actCost;
202  for( i2 = 0; i2 < l2; ++i2 ) {
203  const float64_t actMatch = costs2_0[ i2 ] + ( seq1[i1] == seq2[i2] );
204  const float64_t actGap1 = costs2_0[ i2+1 ] + gapCosts1[ i1 ];
205  const float64_t actGap2 = actCost + gapCosts2[ i2 ];
206  const float64_t actGap = min( actGap1, actGap2 );
207  actCost = min( actMatch, actGap );
208  costs2_1[ i2+1 ] = actCost;
209  }
210  }
211 
212  SG_FREE(gapCosts1);
213  SG_FREE(gapCosts2);
214  SG_FREE(costs2_0);
215  SG_FREE(costs2_1);
216 
217  // return the final cost
218  return actCost;
219 }
220 
221 void CMath::linspace(float64_t* output, float64_t start, float64_t end, int32_t n)
222 {
223  float64_t delta = (end-start) / (n-1);
224  float64_t v = start;
225  index_t i = 0;
226  while ( v <= end )
227  {
228  output[i++] = v;
229  v += delta;
230  }
231  output[n-1] = end;
232 }
233 
234 int CMath::is_nan(double f)
235 {
236 #ifndef HAVE_STD_ISNAN
237 #if (HAVE_DECL_ISNAN == 1) || defined(HAVE_ISNAN)
238  return ::isnan(f);
239 #else
240  return ((f != f) ? 1 : 0);
241 #endif // #if (HAVE_DECL_ISNAN == 1) || defined(HAVE_ISNAN)
242 #endif // #ifndef HAVE_STD_ISNAN
243 
244  return std::isnan(f);
245 }
246 
247 int CMath::is_infinity(double f)
248 {
249 #ifndef HAVE_STD_ISINF
250 #if (HAVE_DECL_ISINF == 1) || defined(HAVE_ISINF)
251  return ::isinf(f);
252 #elif defined(FPCLASS)
253  if (::fpclass(f) == FP_NINF) return -1;
254  else if (::fpclass(f) == FP_PINF) return 1;
255  else return 0;
256 #else
257  if ((f == f) && ((f - f) != 0.0)) return (f < 0.0 ? -1 : 1);
258  else return 0;
259 }
260 #endif // #if (HAVE_DECL_ISINF == 1) || defined(HAVE_ISINF)
261 #endif // #ifndef HAVE_STD_ISINF
262 
263  return std::isinf(f);
264 }
265 
266 int CMath::is_finite(double f)
267 {
268 #ifndef HAVE_STD_ISFINITE
269 #if (HAVE_DECL_ISFINITE == 1) || defined(HAVE_ISFINITE)
270  return ::isfinite(f);
271 #elif defined(HAVE_FINITE)
272  return ::finite(f);
273 #else
274  return ((!std::isnan(f) && !std::isinf(f)) ? 1 : 0);
275 #endif // #if (HAVE_DECL_ISFINITE == 1) || defined(HAVE_ISFINITE)
276 #endif // #ifndef HAVE_STD_ISFINITE
277 
278  return std::isfinite(f);
279 }
280 
281 bool CMath::strtof(const char* str, float32_t* float_result)
282 {
283  ASSERT(str);
284  ASSERT(float_result);
285 
286  SGVector<char> buf(strlen(str)+1);
287 
288  for (index_t i=0; i<buf.vlen-1; i++)
289  buf[i]=tolower(str[i]);
290  buf[buf.vlen-1]='\0';
291 
292  if (strstr(buf, "inf") != NULL)
293  {
294  *float_result = CMath::INFTY;
295 
296  if (strchr(buf,'-') != NULL)
297  *float_result *= -1;
298  return true;
299  }
300 
301  if (strstr(buf, "nan") != NULL)
302  {
303  *float_result = CMath::NOT_A_NUMBER;
304  return true;
305  }
306 
307  char* endptr = buf.vector;
308  *float_result=::strtof(str, &endptr);
309  return endptr != buf.vector;
310 }
311 
312 bool CMath::strtod(const char* str, float64_t* double_result)
313 {
314  ASSERT(str);
315  ASSERT(double_result);
316 
317  SGVector<char> buf(strlen(str)+1);
318 
319  for (index_t i=0; i<buf.vlen-1; i++)
320  buf[i]=tolower(str[i]);
321  buf[buf.vlen-1]='\0';
322 
323  if (strstr(buf, "inf") != NULL)
324  {
325  *double_result = CMath::INFTY;
326 
327  if (strchr(buf,'-') != NULL)
328  *double_result *= -1;
329  return true;
330  }
331 
332  if (strstr(buf, "nan") != NULL)
333  {
334  *double_result = CMath::NOT_A_NUMBER;
335  return true;
336  }
337 
338  char* endptr = buf.vector;
339  *double_result=::strtod(str, &endptr);
340  return endptr != buf.vector;
341 }
342 
343 bool CMath::strtold(const char* str, floatmax_t* long_double_result)
344 {
345  ASSERT(str);
346  ASSERT(long_double_result);
347 
348  SGVector<char> buf(strlen(str)+1);
349 
350  for (index_t i=0; i<buf.vlen-1; i++)
351  buf[i]=tolower(str[i]);
352  buf[buf.vlen-1]='\0';
353 
354  if (strstr(buf, "inf") != NULL)
355  {
356  *long_double_result = CMath::INFTY;
357 
358  if (strchr(buf,'-') != NULL)
359  *long_double_result *= -1;
360  return true;
361  }
362 
363  if (strstr(buf, "nan") != NULL)
364  {
365  *long_double_result = CMath::NOT_A_NUMBER;
366  return true;
367  }
368 
369  char* endptr = buf.vector;
370 
371 // fall back to double on win32 / cygwin since strtold is undefined there
372 #if defined(WIN32) || defined(__CYGWIN__)
373  *long_double_result=::strtod(str, &endptr);
374 #else
375  *long_double_result=::strtold(str, &endptr);
376 #endif
377 
378  return endptr != buf.vector;
379 }
380 
382 {
383  REQUIRE(rel_tolerance > 0 && rel_tolerance < 1.0,
384  "Relative tolerance (%f) should be less than 1.0 and positive\n", rel_tolerance);
385  REQUIRE(is_finite(true_value),
386  "The true_value should be finite\n");
387  float64_t abs_tolerance = rel_tolerance;
388  if (abs(true_value)>0.0)
389  {
390  if (log(abs(true_value)) + log(rel_tolerance) < log(F_MIN_VAL64))
391  abs_tolerance = F_MIN_VAL64;
392  else
393  abs_tolerance = abs(true_value * rel_tolerance);
394  }
395  return abs_tolerance;
396 }
static const float32_t F_MAX_VAL32
Definition: Math.h:2065
static const float64_t MACHINE_EPSILON
Definition: Math.h:2058
static bool strtof(const char *str, float32_t *float_result)
Definition: Math.cpp:281
static uint32_t seed
random generator seed
Definition: Math.h:2079
static int is_finite(double f)
checks whether a float is finite
Definition: Math.cpp:266
static float64_t Align(char *seq1, char *seq2, int32_t l1, int32_t l2, float64_t gapCost)
Definition: Math.cpp:178
static void linspace(float64_t *output, float64_t start, float64_t end, int32_t n=100)
Definition: Math.cpp:221
int32_t index_t
Definition: common.h:62
static bool strtod(const char *str, float64_t *double_result)
Definition: Math.cpp:312
virtual ~CMath()
Destructor - frees logtable.
Definition: Math.cpp:82
static const float64_t INFTY
infinity
Definition: Math.h:2048
static const float64_t MIN_REAL_NUMBER
Definition: Math.h:2062
#define REQUIRE(x,...)
Definition: SGIO.h:206
static const float64_t F_MAX_VAL64
Definition: Math.h:2067
static const float32_t F_MIN_VAL32
Definition: Math.h:2071
static const float32_t F_MIN_NORM_VAL32
Definition: Math.h:2066
static int32_t LOGRANGE
range for logtable: log(1+exp(x)) -LOGRANGE <= x <= 0
Definition: Math.h:2076
static const float64_t ALMOST_NEG_INFTY
almost neg (log) infinity
Definition: Math.h:2052
index_t vlen
Definition: SGVector.h:494
CMath()
Constructor - initializes log-table.
Definition: Math.cpp:65
#define ASSERT(x)
Definition: SGIO.h:201
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:115
double float64_t
Definition: common.h:50
long double floatmax_t
Definition: common.h:51
#define M_PI
workaround for log2 being a define on cygwin
Definition: Math.h:59
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
Definition: Math.h:627
float float32_t
Definition: common.h:49
static float64_t get_abs_tolerance(float64_t true_value, float64_t rel_tolerance)
Definition: Math.cpp:381
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
static int is_infinity(double f)
checks whether a float is infinity
Definition: Math.cpp:247
static int is_nan(double f)
checks whether a float is nan
Definition: Math.cpp:234
static T min(T a, T b)
Definition: Math.h:157
static float64_t exp(float64_t x)
Definition: Math.h:621
#define SG_SINFO(...)
Definition: SGIO.h:173
static const float64_t F_MIN_VAL64
Definition: Math.h:2072
static const float64_t F_MIN_NORM_VAL64
Definition: Math.h:2068
static float64_t log(float64_t v)
Definition: Math.h:922
static const float64_t ALMOST_INFTY
Definition: Math.h:2049
Class which collects generic mathematical functions.
Definition: Math.h:134
static void swap(T &a, T &b)
Definition: Math.h:438
static void sort(int32_t *a, int32_t cols, int32_t sort_col=0)
Definition: Math.cpp:139
#define NAN
Definition: Math.cpp:26
static bool strtold(const char *str, floatmax_t *long_double_result)
Definition: Math.cpp:343
static const float64_t NOT_A_NUMBER
not a number
Definition: Math.h:2046
static const float64_t MAX_REAL_NUMBER
Definition: Math.h:2061
static T abs(T a)
Definition: Math.h:179
static const float64_t PI
Definition: Math.h:2055

SHOGUN Machine Learning Toolbox - Documentation