PLearn 0.1
pl_math.h
Go to the documentation of this file.
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 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines