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

SHOGUN Machine Learning Toolbox - Documentation