SHOGUN  4.1.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/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;
93 #ifdef HAVE_EIGEN3
94  Eigen::Map<const Eigen::VectorXd> ev1(v1,n);
95  Eigen::Map<const Eigen::VectorXd> ev2(v2,n);
96  r = ev1.dot(ev2);
97 #elif HAVE_LAPACK
98  int32_t skip=1;
99  r = cblas_ddot(n, v1, skip, v2, skip);
100 #else
101  for (int32_t i=0; i<n; i++)
102  r+=v1[i]*v2[i];
103 #endif
104  return r;
105 }
106 
107 float32_t CMath::dot(const float32_t* v1, const float32_t* v2, int32_t n)
108 {
109  float32_t r=0;
110 #ifdef HAVE_EIGEN3
111  Eigen::Map<const Eigen::VectorXf> ev1(v1,n);
112  Eigen::Map<const Eigen::VectorXf> ev2(v2,n);
113  r = ev1.dot(ev2);
114 #elif HAVE_LAPACK
115  int32_t skip=1;
116  r = cblas_sdot(n, v1, skip, v2, skip);
117 #else
118  for (int32_t i=0; i<n; i++)
119  r+=v1[i]*v2[i];
120 #endif
121  return r;
122 }
123 
124 #ifdef USE_LOGCACHE
125 int32_t CMath::determine_logrange()
126 {
127  int32_t i;
128  float64_t acc=0;
129  for (i=0; i<50; i++)
130  {
131  acc=((float64_t)log(1+((float64_t)exp(-float64_t(i)))));
132  if (acc<=(float64_t)LOG_TABLE_PRECISION)
133  break;
134  }
135 
136  SG_SINFO("determined range for x in table log(1+exp(-x)) is:%d (error:%G)\n",i,acc)
137  return i;
138 }
139 
140 int32_t CMath::determine_logaccuracy(int32_t range)
141 {
142  range=MAX_LOG_TABLE_SIZE/range/((int)sizeof(float64_t));
143  SG_SINFO("determined accuracy for x in table log(1+exp(-x)) is:%d (error:%G)\n",range,1.0/(double) range)
144  return range;
145 }
146 
147 //init log table of form log(1+exp(x))
148 void CMath::init_log_table()
149 {
150  for (int32_t i=0; i< LOGACCURACY*LOGRANGE; i++)
151  logtable[i]=log(1+exp(float64_t(-i)/float64_t(LOGACCURACY)));
152 }
153 #endif
154 
155 void CMath::sort(int32_t *a, int32_t cols, int32_t sort_col)
156 {
157  int32_t changed=1;
158  if (a[0]==-1) return;
159  while (changed)
160  {
161  changed=0; int32_t i=0;
162  while ((a[(i+1)*cols]!=-1) && (a[(i+1)*cols+1]!=-1)) // to be sure
163  {
164  if (a[i*cols+sort_col]>a[(i+1)*cols+sort_col])
165  {
166  for (int32_t j=0; j<cols; j++)
167  CMath::swap(a[i*cols+j],a[(i+1)*cols+j]);
168  changed=1;
169  };
170  i++;
171  };
172  };
173 }
174 
175 void CMath::sort(float64_t *a, int32_t* idx, int32_t N)
176 {
177  int32_t changed=1;
178  while (changed)
179  {
180  changed=0;
181  for (int32_t i=0; i<N-1; i++)
182  {
183  if (a[i]>a[i+1])
184  {
185  swap(a[i],a[i+1]) ;
186  swap(idx[i],idx[i+1]) ;
187  changed=1 ;
188  } ;
189  } ;
190  } ;
191 
192 }
193 
195  char* seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost)
196 {
197  float64_t actCost=0 ;
198  int32_t i1, i2 ;
199  float64_t* const gapCosts1 = SG_MALLOC(float64_t, l1 );
200  float64_t* const gapCosts2 = SG_MALLOC(float64_t, l2 );
201  float64_t* costs2_0 = SG_MALLOC(float64_t, l2 + 1 );
202  float64_t* costs2_1 = SG_MALLOC(float64_t, l2 + 1 );
203 
204  // initialize borders
205  for( i1 = 0; i1 < l1; ++i1 ) {
206  gapCosts1[ i1 ] = gapCost * i1;
207  }
208  costs2_1[ 0 ] = 0;
209  for( i2 = 0; i2 < l2; ++i2 ) {
210  gapCosts2[ i2 ] = gapCost * i2;
211  costs2_1[ i2+1 ] = costs2_1[ i2 ] + gapCosts2[ i2 ];
212  }
213  // compute alignment
214  for( i1 = 0; i1 < l1; ++i1 ) {
215  swap( costs2_0, costs2_1 );
216  actCost = costs2_0[ 0 ] + gapCosts1[ i1 ];
217  costs2_1[ 0 ] = actCost;
218  for( i2 = 0; i2 < l2; ++i2 ) {
219  const float64_t actMatch = costs2_0[ i2 ] + ( seq1[i1] == seq2[i2] );
220  const float64_t actGap1 = costs2_0[ i2+1 ] + gapCosts1[ i1 ];
221  const float64_t actGap2 = actCost + gapCosts2[ i2 ];
222  const float64_t actGap = min( actGap1, actGap2 );
223  actCost = min( actMatch, actGap );
224  costs2_1[ i2+1 ] = actCost;
225  }
226  }
227 
228  SG_FREE(gapCosts1);
229  SG_FREE(gapCosts2);
230  SG_FREE(costs2_0);
231  SG_FREE(costs2_1);
232 
233  // return the final cost
234  return actCost;
235 }
236 
237 void CMath::linspace(float64_t* output, float64_t start, float64_t end, int32_t n)
238 {
239  float64_t delta = (end-start) / (n-1);
240  float64_t v = start;
241  index_t i = 0;
242  while ( v <= end )
243  {
244  output[i++] = v;
245  v += delta;
246  }
247  output[n-1] = end;
248 }
249 
250 int CMath::is_nan(double f)
251 {
252 #ifndef HAVE_STD_ISNAN
253 #if (HAVE_DECL_ISNAN == 1) || defined(HAVE_ISNAN)
254  return ::isnan(f);
255 #else
256  return ((f != f) ? 1 : 0);
257 #endif // #if (HAVE_DECL_ISNAN == 1) || defined(HAVE_ISNAN)
258 #endif // #ifndef HAVE_STD_ISNAN
259 
260  return std::isnan(f);
261 }
262 
263 int CMath::is_infinity(double f)
264 {
265 #ifndef HAVE_STD_ISINF
266 #if (HAVE_DECL_ISINF == 1) || defined(HAVE_ISINF)
267  return ::isinf(f);
268 #elif defined(FPCLASS)
269  if (::fpclass(f) == FP_NINF) return -1;
270  else if (::fpclass(f) == FP_PINF) return 1;
271  else return 0;
272 #else
273  if ((f == f) && ((f - f) != 0.0)) return (f < 0.0 ? -1 : 1);
274  else return 0;
275 }
276 #endif // #if (HAVE_DECL_ISINF == 1) || defined(HAVE_ISINF)
277 #endif // #ifndef HAVE_STD_ISINF
278 
279  return std::isinf(f);
280 }
281 
282 int CMath::is_finite(double f)
283 {
284 #ifndef HAVE_STD_ISFINITE
285 #if (HAVE_DECL_ISFINITE == 1) || defined(HAVE_ISFINITE)
286  return ::isfinite(f);
287 #elif defined(HAVE_FINITE)
288  return ::finite(f);
289 #else
290  return ((!std::isnan(f) && !std::isinf(f)) ? 1 : 0);
291 #endif // #if (HAVE_DECL_ISFINITE == 1) || defined(HAVE_ISFINITE)
292 #endif // #ifndef HAVE_STD_ISFINITE
293 
294  return std::isfinite(f);
295 }
296 
297 bool CMath::strtof(const char* str, float32_t* float_result)
298 {
299  ASSERT(str);
300  ASSERT(float_result);
301 
302  SGVector<char> buf(strlen(str)+1);
303 
304  for (index_t i=0; i<buf.vlen-1; i++)
305  buf[i]=tolower(str[i]);
306  buf[buf.vlen-1]='\0';
307 
308  if (strstr(buf, "inf") != NULL)
309  {
310  *float_result = CMath::INFTY;
311 
312  if (strchr(buf,'-') != NULL)
313  *float_result *= -1;
314  return true;
315  }
316 
317  if (strstr(buf, "nan") != NULL)
318  {
319  *float_result = CMath::NOT_A_NUMBER;
320  return true;
321  }
322 
323  char* endptr = buf.vector;
324  *float_result=::strtof(str, &endptr);
325  return endptr != buf.vector;
326 }
327 
328 bool CMath::strtod(const char* str, float64_t* double_result)
329 {
330  ASSERT(str);
331  ASSERT(double_result);
332 
333  SGVector<char> buf(strlen(str)+1);
334 
335  for (index_t i=0; i<buf.vlen-1; i++)
336  buf[i]=tolower(str[i]);
337  buf[buf.vlen-1]='\0';
338 
339  if (strstr(buf, "inf") != NULL)
340  {
341  *double_result = CMath::INFTY;
342 
343  if (strchr(buf,'-') != NULL)
344  *double_result *= -1;
345  return true;
346  }
347 
348  if (strstr(buf, "nan") != NULL)
349  {
350  *double_result = CMath::NOT_A_NUMBER;
351  return true;
352  }
353 
354  char* endptr = buf.vector;
355  *double_result=::strtod(str, &endptr);
356  return endptr != buf.vector;
357 }
358 
359 bool CMath::strtold(const char* str, floatmax_t* long_double_result)
360 {
361  ASSERT(str);
362  ASSERT(long_double_result);
363 
364  SGVector<char> buf(strlen(str)+1);
365 
366  for (index_t i=0; i<buf.vlen-1; i++)
367  buf[i]=tolower(str[i]);
368  buf[buf.vlen-1]='\0';
369 
370  if (strstr(buf, "inf") != NULL)
371  {
372  *long_double_result = CMath::INFTY;
373 
374  if (strchr(buf,'-') != NULL)
375  *long_double_result *= -1;
376  return true;
377  }
378 
379  if (strstr(buf, "nan") != NULL)
380  {
381  *long_double_result = CMath::NOT_A_NUMBER;
382  return true;
383  }
384 
385  char* endptr = buf.vector;
386 
387 // fall back to double on win32 / cygwin since strtold is undefined there
388 #if defined(WIN32) || defined(__CYGWIN__)
389  *long_double_result=::strtod(str, &endptr);
390 #else
391  *long_double_result=::strtold(str, &endptr);
392 #endif
393 
394  return endptr != buf.vector;
395 }
396 
398 {
399  REQUIRE(rel_tolerance > 0 && rel_tolerance < 1.0,
400  "Relative tolerance (%f) should be less than 1.0 and positive\n", rel_tolerance);
401  REQUIRE(is_finite(true_value),
402  "The true_value should be finite\n");
403  float64_t abs_tolerance = rel_tolerance;
404  if (abs(true_value)>0.0)
405  {
406  if (log(abs(true_value)) + log(rel_tolerance) < log(F_MIN_VAL64))
407  abs_tolerance = F_MIN_VAL64;
408  else
409  abs_tolerance = abs(true_value * rel_tolerance);
410  }
411  return abs_tolerance;
412 }

SHOGUN Machine Learning Toolbox - Documentation