PLearn 0.1
|
00001 // -*- C++ -*- 00002 00003 // PLearn (A C++ Machine Learning Library) 00004 // Copyright (C) 1998 Pascal Vincent 00005 // Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal 00006 // 00007 00008 // Redistribution and use in source and binary forms, with or without 00009 // modification, are permitted provided that the following conditions are met: 00010 // 00011 // 1. Redistributions of source code must retain the above copyright 00012 // notice, this list of conditions and the following disclaimer. 00013 // 00014 // 2. Redistributions in binary form must reproduce the above copyright 00015 // notice, this list of conditions and the following disclaimer in the 00016 // documentation and/or other materials provided with the distribution. 00017 // 00018 // 3. The name of the authors may not be used to endorse or promote 00019 // products derived from this software without specific prior written 00020 // permission. 00021 // 00022 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR 00023 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 00024 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN 00025 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 00027 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00028 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00029 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00030 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00031 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00032 // 00033 // This file is part of the PLearn library. For more information on the PLearn 00034 // library, go to the PLearn Web site at www.plearn.org 00035 00036 00037 /* ******************************************************* 00038 * $Id: pl_math.h 9593 2008-10-21 16:54:57Z nouiz $ 00039 * This file is part of the PLearn library. 00040 ******************************************************* */ 00041 00042 00045 #ifndef pl_math_INC 00046 #define pl_math_INC 00047 00048 #include <cmath> 00049 #include <cfloat> 00050 #include <climits> 00051 #include <plearn/base/plerror.h> 00052 00053 #ifdef WIN32 00054 #include <limits> 00055 #endif 00056 00057 namespace PLearn { 00058 using namespace std; 00059 00060 // Define 'real' constants. 00061 // NB: the max is not really the max, but it's large enough for most applications. 00062 #if defined(USEDOUBLE) 00063 //#define real double 00064 typedef double real; 00065 #define REAL_MAX DBL_MAX 00066 #define REAL_EPSILON DBL_EPSILON 00067 #define ABSOLUTE_TOLERANCE 1e-8 00068 #define RELATIVE_TOLERANCE ABSOLUTE_TOLERANCE 00069 #elif defined(USEFLOAT) 00070 //#define real float 00071 typedef float real; 00072 #define REAL_MAX FLT_MAX 00073 #define REAL_EPSILON FLT_EPSILON 00074 #define ABSOLUTE_TOLERANCE 1e-6 00075 #define RELATIVE_TOLERANCE ABSOLUTE_TOLERANCE 00076 #else 00077 #error You must define either USEDOUBLE or USEFLOAT 00078 #endif 00079 00080 union _plearn_nan_type { unsigned char c[4]; float d; }; 00081 extern _plearn_nan_type plearn_nan; 00082 00085 #if defined(NAN) && !defined(WIN32) && !defined(__INTEL_COMPILER) 00086 #define MISSING_VALUE NAN 00087 #else 00088 #define MISSING_VALUE (plearn_nan.d) 00089 #endif 00090 00093 #ifdef WIN32 00094 #ifdef INFINITY 00095 #undef INFINITY 00096 #endif 00097 #define INFINITY numeric_limits<real>::infinity() 00098 #endif 00099 00101 #ifndef M_SQRT2 00102 #define M_SQRT2 1.41421356237309504880 00103 #endif 00104 00106 #ifndef M_PI 00107 #define M_PI 3.14159265358979323846 00108 #endif 00109 00113 inline real log_force_nan_if_negative(real a) 00114 { return a < 0 ? MISSING_VALUE : std::log(a); } 00115 #if defined(__CYGWIN__) && defined(__GNUC__) 00116 #define pl_log log_force_nan_if_negative 00117 #else 00118 #define pl_log std::log 00119 #endif 00120 // Every now and then, uncomment the line below in order to spot places where 00121 // 'log' is used instead of 'pl_log'. 00122 // #define log USE_pl_log_INSTEAD_OF_log 00123 00124 using namespace std; 00125 00126 using std::sqrt; 00127 using std::pow; 00128 using std::exp; 00129 using std::tanh; 00130 using std::abs; 00131 00132 inline real sign(real a) { 00133 if (a>0) return 1; 00134 if (a<0) return -1; 00135 return 0; 00136 } 00137 inline real positive(real a) { if (a>0) return a; return 0; } 00138 inline real negative(real a) { if (a<0) return a; return 0; } 00139 00141 #if !defined(MIN) 00142 # define MIN(a,b) ((a)<(b)?(a):(b)) 00143 #endif 00144 00145 #if !defined(MAX) 00146 # define MAX(a,b) ((a)>(b)?(a):(b)) 00147 #endif 00148 00149 #if !defined(SIGN) 00150 # define SIGN(a) ((a)>=0?1:-1) 00151 #endif 00152 00153 #if !defined(Pi) 00154 # define Pi 3.141592653589793 00155 #endif 00156 00157 #if !defined(LogPi) 00158 # define LogPi 1.14472988585 00159 #endif 00160 00161 #if !defined(Log2Pi) 00162 # define Log2Pi 1.837877066409 00163 #endif 00164 00165 #if !defined(LOG_2) 00166 # define LOG_2 0.693147180559945 00167 #endif 00168 00169 #if !defined(LOG_INIT) 00170 # define LOG_INIT -REAL_MAX 00171 #endif 00172 00173 #if !defined(MINUS_LOG_THRESHOLD) 00174 # define MINUS_LOG_THRESHOLD -18.42 00175 #endif 00176 00177 #if !defined(DEG2RAD) 00178 # define DEG2RAD Pi/180.0 00179 #endif 00180 00181 #if !defined(RAD2DEG) 00182 # define RAD2DEG 57.29578 00183 #endif 00184 00185 // Specific definitions for some non-Linux OS. 00186 #if defined(_MSC_VER) // Microsoft Visual Studio. 00187 #define isnan(x) (_isnan(x) != 0) 00188 00189 // 'isinf' cannot be just redefined, as there is apparently no existing 00190 // function exhibiting the same behavior. 00191 template<class T> 00192 inline int isinf(T x) 00193 { 00194 return x == numeric_limits<T>::infinity() ? 1 00195 : -x == numeric_limits<T>::infinity() ? -1 00196 : 0; 00197 } 00198 00199 #define finite _finite 00200 #define round(x) int(x + 0.5) 00201 #define rint round 00202 #define log1p(x) pl_log(1.0 + x) 00203 #endif 00204 00206 bool is_equal(real a, real b, real absolute_tolerance_threshold = 1.0, 00207 real absolute_tolerance = ABSOLUTE_TOLERANCE, 00208 real relative_tolerance = RELATIVE_TOLERANCE); 00209 00212 inline bool fast_is_equal(real a, real b, real absolute_tolerance_threshold = 1.0, 00213 real absolute_tolerance = ABSOLUTE_TOLERANCE, 00214 real relative_tolerance = RELATIVE_TOLERANCE) 00215 { 00216 // Check for 'nan' or 'inf' values in debug mode. 00217 #ifdef BOUNDCHECK 00218 if (isnan(a) || isinf(a) || isnan(b) || isinf(b)) 00219 PLERROR("In fast_is_equal - Either 'a' or 'b' is 'nan' or 'inf'"); 00220 #endif 00221 real a_abs = fabs(a); 00222 real b_abs = fabs(b); 00223 if (a_abs < absolute_tolerance_threshold && b_abs < absolute_tolerance_threshold) 00224 return (fabs(a-b) <= absolute_tolerance); 00225 real diff_abs = fabs(a - b); 00226 return diff_abs <= relative_tolerance*a_abs && diff_abs <= relative_tolerance*b_abs; 00227 } 00228 00234 inline bool fast_exact_is_equal(real a, real b) 00235 { 00236 // TODO This could (and should) be optimized. 00237 return (a <= b && b <= a); 00238 } 00239 00241 inline bool fast_is_more(real a, real b, real absolute_tolerance_threshold = 1.0, 00242 real absolute_tolerance = ABSOLUTE_TOLERANCE, 00243 real relative_tolerance = RELATIVE_TOLERANCE) 00244 { 00245 return a>b && !fast_is_equal(a,b,absolute_tolerance_threshold,absolute_tolerance,relative_tolerance); 00246 } 00247 00249 inline bool fast_is_less(real a, real b, real absolute_tolerance_threshold = 1.0, 00250 real absolute_tolerance = ABSOLUTE_TOLERANCE, 00251 real relative_tolerance = RELATIVE_TOLERANCE) 00252 { 00253 return a<b && !fast_is_equal(a,b,absolute_tolerance_threshold,absolute_tolerance,relative_tolerance); 00254 } 00255 00257 inline bool is_more(real a, real b, real absolute_tolerance_threshold = 1.0, 00258 real absolute_tolerance = ABSOLUTE_TOLERANCE, 00259 real relative_tolerance = RELATIVE_TOLERANCE) 00260 { 00261 return a>b && !is_equal(a,b,absolute_tolerance_threshold,absolute_tolerance,relative_tolerance); 00262 00263 } 00264 00265 template<class T> 00266 inline T square(const T& x) 00267 { return x*x; } 00268 00269 // Wish I could get rid of this, but it's used a s some function pointer 00270 // in calls to some apply function in RandomVar. So I'd need to change 00271 // those calls first to sth more stl like. 00272 real square_f(real x); 00273 00274 template<class T> 00275 inline T two(const T& x) 00276 { return x+x; } 00277 00278 #define TANHTABLESIZE 5000 00279 #define MAXTANHX 10. 00280 00281 class PLMathInitializer 00282 { 00283 public: 00284 PLMathInitializer(); 00285 ~PLMathInitializer(); 00286 }; 00287 00288 /* 00289 #define DOUBLE2INT(i,d) {maniax = ((d)+6755399441055744.0); i=*((int *)(&maniax));} 00290 00291 inline int double2int(double d) 00292 { 00293 double maniax = d+6755399441055744.0; 00294 return *((int *)&maniax); 00295 } 00296 */ 00297 00298 #if defined(LINUX) && !defined(__INTEL_COMPILER) // note: intel compiler on SGI does not like that 00299 #define DOUBLE_TO_INT(in,out) __asm__ __volatile__ ("fistpl %0" : "=m" (out) : "t" (in) : "st") 00300 #else 00301 #define DOUBLE_TO_INT(in, out) out = int(round(in)) 00302 #endif 00303 00304 extern float tanhtable[TANHTABLESIZE]; 00305 extern PLMathInitializer pl_math_initializer; 00306 00307 inline real fasttanh(const real& x) 00308 { 00309 if (isnan(x)) return x; // tanh(nan)=nan 00310 int is_inf=isinf(x); 00311 if (is_inf>0) return 1; // tanh(inf)=1 00312 if (is_inf<0) return -1; // tanh(-inf)=-1 00313 if(x>0) 00314 { 00315 if(x>MAXTANHX) 00316 return real(tanhtable[TANHTABLESIZE-1]); 00317 else 00318 { 00319 int i; 00320 DOUBLE_TO_INT( double(x*((TANHTABLESIZE-1)/MAXTANHX)), i); 00321 return real(tanhtable[i]); 00322 } 00323 } 00324 else 00325 { 00326 real nx = -x; 00327 if(nx>MAXTANHX) 00328 return real(-tanhtable[TANHTABLESIZE-1]); 00329 else 00330 { 00331 int i; 00332 DOUBLE_TO_INT( double(nx*((TANHTABLESIZE-1)/MAXTANHX)), i); 00333 return real(-tanhtable[i]); 00334 } 00335 } 00336 } 00337 00338 inline real fastsigmoid(const real& x) 00339 { return (real)0.5*(fasttanh(0.5*x)+1.); } 00340 00341 00342 // These are quadratic approximations to tanh and sigmoid 00343 00344 /* 00345 inline real ultrafasttanh(const real& x) 00346 { 00347 if(x>1.92033) return 0.96016; 00348 else if (x>0) return 0.96016 - 0.26037 * square(x - 1.92033); 00349 else if (x<=-1.92033) return -0.96016; 00350 else return 0.26037 * square(x + 1.92033) - 0.96016; 00351 } 00352 00353 inline real ultrafastsigmoid(const real& x) 00354 { return (real)0.5*(ultrafasttanh(0.5*x)+1.); } 00355 */ 00356 00357 inline real ultrafasttanh(const real& x) 00358 { 00359 if(x>=0) 00360 return (x<1.7 ? (1.5*x/(1+x)) : ( x<3 ? (0.935409070603099 + 0.0458812946797165*(x-1.7)) :0.99505475368673)); 00361 else 00362 { 00363 real xx = -x; 00364 return -(xx<1.7 ? (1.5*xx/(1+xx)) : ( xx<3 ? (0.935409070603099 + 0.0458812946797165*(xx-1.7)) :0.99505475368673)); 00365 } 00366 } 00367 00368 /* 00369 inline real ultrafasttanh(const real& x) 00370 { 00371 return x/(1+fabs(x)); 00372 } 00373 */ 00374 00375 inline real ultrafastsigmoid(const real& x) 00376 { 00377 // return 0.5*x / (1. + fabs(x)) + 0.5; 00378 return (real)0.5*(ultrafasttanh(0.5*x)+1.); 00379 // return fastsigmoid(x); 00380 } 00381 00382 // target is -1 or 1 00383 inline real hinge_loss(const real& output, int target) 00384 { 00385 real off_margin = 1-target*output; 00386 return off_margin>0?off_margin:0; 00387 } 00388 00389 // return d(hinge_loss(output,target))/doutput 00390 // target is -1 or 1 00391 inline real d_hinge_loss(const real& output, int target) 00392 { 00393 real margin = target*output; 00394 if (margin<1) return -target; 00395 return 0; 00396 } 00397 00400 template<class T> 00401 inline bool is_missing(const T& x) { return false; } 00402 00404 inline bool is_missing(double x) { return isnan(x); } 00405 00407 inline bool is_missing(float x) { return isnan(x); } 00408 00409 inline bool is_integer(real x) { return fast_exact_is_equal(real(int(x)), x); } 00410 00411 inline real FABS(real x) 00412 { return x>=0. ?x :-x; } 00413 00414 #define FSWAP(a,b) do {real _c; _c = *(a); *(a) = *(b); *(b) = _c;} while(0) 00415 00416 inline real mypow(real x, real p) 00417 { return fast_exact_is_equal(x, 0) ? 0 :pow(x,p); } 00418 00419 inline real ipow(real x, int p) 00420 { 00421 PLASSERT( p >= 0 ); 00422 real result = 1.0; 00423 while(p--) 00424 result *= x; 00425 return result; 00426 } 00427 00428 inline int ipow(int x, int p) 00429 { 00430 PLASSERT( p >= 0 ); 00431 int result = 1; 00432 while(p--) 00433 result *= x; 00434 return result; 00435 } 00436 00438 inline real sigmoid(real x) 00439 { return (real)0.5*(tanh(0.5*x)+1.); } 00440 00443 inline real is_positive(real x) { return x>0? 1 : 0; } 00444 00446 inline real inverse_sigmoid(real x) 00447 { 00448 #ifdef BOUNDCHECK 00449 if (x < 0. || x > 1. || is_missing(x)) 00450 PLERROR("In inv_sigmoid_value: a should be in [0,1]"); 00451 #endif 00452 // We specify an absolute 1e-5 threshold to have the same behavior as with 00453 // the old FEQUAL macro. 00454 if (fast_is_equal(x,0.,REAL_MAX,1e-5)) 00455 return -88.; 00456 else if (fast_is_equal(x,1.,REAL_MAX,1e-5)) 00457 return 14.5; 00458 else 00459 return real(-pl_log(1./x - 1.)); 00460 } 00461 00463 inline real softplus(real x) 00464 { 00465 if(x<=-30.) 00466 return 0.0; 00467 else if(x>=30.) 00468 return x; 00469 else 00470 return log1p(exp(x)); 00472 } 00473 00474 inline real tabulated_softplus(real x) 00475 { 00476 static const int n_softplus_values = 1000000; 00477 static const real min_softplus_arg = -10; 00478 static const real max_softplus_arg = 10; 00479 static const real softplus_delta = (n_softplus_values-1)/(max_softplus_arg-min_softplus_arg); 00480 static real softplus_values[n_softplus_values]; 00481 static bool computed_softplus_table = false; 00482 if (isnan(x)) return x; // softplus(nan)=nan 00483 int is_inf=isinf(x); 00484 if (is_inf>0) return x; // softplus(inf)=inf 00485 if (is_inf<0) return 0; // softplus(-inf)=0 00486 if (!computed_softplus_table) 00487 { 00488 real y=min_softplus_arg; 00489 real dy=1.0/softplus_delta; 00490 for (int i=0;i<n_softplus_values;i++,y+=dy) 00491 softplus_values[i] = softplus(y); 00492 computed_softplus_table=true; 00493 } 00494 if (x<min_softplus_arg) return 0; 00495 if (x>max_softplus_arg) return x; 00496 int bin = int(rint((x-min_softplus_arg)*softplus_delta)); 00497 return softplus_values[bin]; 00498 } 00499 00501 inline real inverse_softplus(real y) 00502 { 00503 if (y<0) 00504 return MISSING_VALUE; 00505 if (y>=30) 00506 return y; 00507 if (fast_exact_is_equal(y, 0)) 00508 return -30; 00509 return pl_log(exp(y)-1); 00510 } 00511 00512 inline real hard_slope(real x, real left=0, real right=1) 00513 { 00514 if (x<left) return 0; 00515 if (x>right) return 1; 00516 return (x-left)/(right-left); 00517 } 00518 00520 inline real log_sigmoid(real x) 00521 { 00522 return -softplus(-x); 00523 } 00524 00525 // as smoothness-->infty this becomes the linear by part function that 00526 // is 0 in [-infty,left], linear in [left,right], and 1 in [right,infty]. 00527 // For finite smoothness, it is a smoother function, always with value in the interval [0,1]. 00528 // It is always monotonically increasing wrt x (positive derivative in x). 00529 inline real soft_slope(real x, real smoothness=1, real left=0, real right=1) 00530 { 00531 if (fast_exact_is_equal(smoothness, 0)) 00532 return 0.5; 00533 if (smoothness>1000) 00534 return hard_slope(x,left,right); 00535 return 1 + (softplus(-smoothness*(x-left))-softplus(-smoothness*(x-right)))/(smoothness*(right-left)); 00536 } 00537 00538 inline real tabulated_soft_slope(real x, real smoothness=1, real left=0, real right=1) 00539 { 00540 if (fast_exact_is_equal(smoothness, 0)) 00541 return 0.5; 00542 if (smoothness>1000) 00543 return hard_slope(x,left,right); 00544 return 1 + (tabulated_softplus(-smoothness*(x-left))-tabulated_softplus(-smoothness*(x-right)))/(smoothness*(right-left)); 00545 } 00546 00547 // This is the derivative of soft_slope with respect to x. 00548 inline real d_soft_slope(real x, real smoothness=1, real left=0, real right=1) 00549 { 00550 // note that d(softplus(z))/dz = sigmoid(z) 00551 return (-sigmoid(-smoothness*(x-left))+sigmoid(-smoothness*(x-right)))/(right-left); 00552 } 00553 00555 inline int n_choose(int M,int N) 00556 { 00557 int k=M-N; 00558 float res=1; 00559 int i; 00560 for (i=1;i<=k;i++) { 00561 res *= (i+N)/(float)i; 00562 } 00563 return (int)(res+0.499); 00564 } 00565 00566 real safeflog(real a); 00567 inline real safelog(real a) { return safeflog(a); } 00568 real safeexp(real a); 00569 00570 real log(real base, real a); 00571 real logtwo(real a); 00572 real safeflog(real base, real a); 00573 real safeflog2(real a); 00574 00575 typedef real (*tRealFunc)(real); 00576 typedef real (*tRealReadFunc)(real,real); 00577 00580 real logadd(double log_a, double log_b); 00581 00582 #ifdef USEFLOAT 00583 00584 real logadd(real log_a, real log_b); 00585 #endif 00586 00588 real logsub(real log_a, real log_b); 00589 00590 <<<<<<< .mine 00593 ======= 00596 >>>>>>> .r9426 00599 real dilogarithm(real x); 00600 00601 inline real softplus_primitive(real x) { 00602 return -dilogarithm(-exp(x)); 00603 } 00604 00605 real tabulated_softplus_primitive(real x); 00606 00607 real hard_slope_integral(real left=0, real right=1, real a=0, real b=1); 00608 00609 00610 // integral of the soft_slope function between a and b 00611 real soft_slope_integral(real smoothness=1, real left=0, real right=1, real a=0, real b=1); 00612 real tabulated_soft_slope_integral(real smoothness=1, real left=0, real right=1, real a=0, real b=1); 00613 00614 } // end of namespace PLearn 00615 00616 00617 #endif 00618 00619 00620 /* 00621 Local Variables: 00622 mode:c++ 00623 c-basic-offset:4 00624 c-file-style:"stroustrup" 00625 c-file-offsets:((innamespace . 0)(inline-open . 0)) 00626 indent-tabs-mode:nil 00627 fill-column:79 00628 End: 00629 */ 00630 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :