PLearn 0.1
pl_math.cc
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.cc 9203 2008-07-03 16:39:04Z nouiz $
00039  * This file is part of the PLearn library.
00040  ******************************************************* */
00041 
00042 
00045 #include "pl_math.h"
00046 
00047 namespace PLearn {
00048 using namespace std;
00049 
00050 
00051 # ifdef BIGENDIAN
00052 _plearn_nan_type plearn_nan = { {0x7f, 0xc0, 0, 0} };
00053 # endif
00054 # ifdef LITTLEENDIAN
00055 _plearn_nan_type plearn_nan = { {0, 0, 0xc0, 0x7f} };
00056 # endif
00057 
00058 float tanhtable[TANHTABLESIZE];
00059 
00060 PLMathInitializer::PLMathInitializer()
00061 {
00063     real scaling = MAXTANHX/(TANHTABLESIZE-1);
00064     for(int i=0; i<TANHTABLESIZE; i++)
00065         tanhtable[i] = (float) tanh(i*scaling);
00066 }
00067   
00068 PLMathInitializer::~PLMathInitializer() 
00069 {}
00070   
00071 PLMathInitializer pl_math_initializer;
00072 
00074 // is_equal //
00076 bool is_equal(real a, real b, real absolute_tolerance_threshold, 
00077               real absolute_tolerance,
00078               real relative_tolerance)
00079 {
00080     if (isnan(a)){
00081         if (isnan(b))
00082             return true;
00083         else
00084             return false;
00085     }
00086     if (isnan(b))
00087         return false;
00088     if (int inf_a = isinf(a))
00089         return inf_a == isinf(b);
00090     if (isinf(b))
00091         return false;
00092     return fast_is_equal(a, b, absolute_tolerance_threshold, absolute_tolerance, relative_tolerance);
00093 }
00094 
00095 real safeflog(real a)
00096 {
00097     if (a < 0.0)
00098         PLERROR("safeflog: negative argument (%f)", a);
00099     if (a < 1e-25)
00100         return -57.5;
00101     else return (real)pl_log((double)a);
00102 }
00103 
00104 real safeexp(real a)
00105 {
00106 #ifdef USEDOUBLE
00107     if (a < -300) return 0;
00108     if (a > 300) return 1e38;
00109 #else
00110     if (a < -87) return 0;
00111     if (a > 43) return 5e18;
00112 #endif
00113     return exp(a);
00114 }
00115 
00116 real log(real base, real a)
00117 {
00118     return pl_log(a) / pl_log(base);
00119 }
00120 
00121 real logtwo(real a)
00122 {
00123     return pl_log(a) / LOG_2;
00124 }
00125 
00126 real safeflog(real base, real a)
00127 {
00128     return safeflog(a) / safeflog(base);
00129 }
00130 
00131 real safeflog2(real a)
00132 {
00133     return safeflog(a) / LOG_2;
00134 }
00135 
00136 real tabulated_softplus_primitive(real x) {
00137     static const int n_softplus_primitive_values = 10000;
00138     static const real min_softplus_primitive_arg = -20;
00139     static const real max_softplus_primitive_arg = 10;
00140     static const real max_offset = max_softplus_primitive_arg*max_softplus_primitive_arg*0.5;
00141     static const real softplus_primitive_delta = (n_softplus_primitive_values-1)/(max_softplus_primitive_arg-min_softplus_primitive_arg);
00142     static real softplus_primitive_values[n_softplus_primitive_values];
00143     static bool computed_softplus_primitive_table = false;
00144     if (!computed_softplus_primitive_table)
00145     {
00146         real y=min_softplus_primitive_arg;
00147         real dy=1.0/softplus_primitive_delta;
00148         for (int i=0;i<n_softplus_primitive_values;i++,y+=dy)
00149             softplus_primitive_values[i] = softplus_primitive(y);
00150         computed_softplus_primitive_table=true;
00151     }
00152     if (x<min_softplus_primitive_arg) return 0;
00153     if (x>max_softplus_primitive_arg) return softplus_primitive_values[n_softplus_primitive_values-1]+x*x*0.5 - max_offset;
00154     int bin = int(rint((x-min_softplus_primitive_arg)*softplus_primitive_delta));
00155     return softplus_primitive_values[bin];
00156 }
00157 
00158 // compute log(exp(log_a)+exp(log_b)) without losing too much precision
00159 real logadd(double log_a, double log_b)
00160 {
00161     if (log_a < log_b)
00162     { // swap them
00163         double tmp = log_a;
00164         log_a = log_b;
00165         log_b = tmp;
00166     } else if (fast_exact_is_equal(log_a, log_b)) {
00167         // Special case when log_a == log_b. In particular this works when both
00168         // log_a and log_b are (+-) INFINITY: it will return (+-) INFINITY
00169         // instead of NaN.
00170         return LOG_2 + log_a;
00171     }
00172     double negative_absolute_difference = log_b - log_a;
00173     if (negative_absolute_difference < MINUS_LOG_THRESHOLD)
00174         return real(log_a);
00175     return (real)(log_a + log1p(exp(negative_absolute_difference)));
00176 }
00177 
00178 #ifdef USEFLOAT
00179 real logadd(real log_a, real log_b)
00180 {
00181     return logadd(double(log_a), double(log_b));
00182 }
00183 #endif
00184 
00185 real square_f(real x)
00186 { return x*x; }
00187 
00188 // compute log(exp(log_a)-exp(log_b)) without losing too much precision
00189 real logsub(real log_a, real log_b)
00190 {
00191     if (log_a < log_b)
00192         PLERROR("log_sub: log_a (%f) should be greater than log_b (%f)", log_a, log_b);
00193  
00194     real negative_absolute_difference = log_b - log_a;
00195  
00196     // We specify an absolute 1e-5 threshold to have the same behavior as with
00197     // the old FEQUAL macro.
00198     if (fast_is_equal(log_a, log_b, REAL_MAX, 1e-5))
00199         return -REAL_MAX;
00200     else if (negative_absolute_difference < MINUS_LOG_THRESHOLD)
00201         return log_a;
00202     else
00203         return log_a + log1p(-exp(negative_absolute_difference));
00204 }
00205 
00206 real small_dilogarithm(real x)
00207 {
00208     // TODO Deal with x == 0.
00209     real somme = x;
00210     real prod = x;
00211     int i=2;
00212     for (;i<=999;i++)
00213     {
00214         real coef = (i-1.0)/i;
00215         prod *= x*coef*coef;
00216         somme += prod;
00217         if (fabs(prod/somme)<1e-16) break; // tolerance
00218     }
00219     static bool warning_was_raised=false;
00220     if (i==1000 && !warning_was_raised) 
00221     {
00222         warning_was_raised=true;
00223         PLWARNING("dilogarithm (%f): insufficient precision", x);
00224     }
00225     return somme;
00226 }
00227 
00228 real positive_dilogarithm(real x)
00229 {
00230     if (x<0.5)
00231         return small_dilogarithm(x);
00232     else if (x<1.0)
00233         return Pi*Pi/6.0 - small_dilogarithm(1.0-x) - pl_log(x)*pl_log(1-x);
00234     else if (fast_exact_is_equal(x, 1.0))
00235         return Pi*Pi/6.0;
00236     else if (x<=1.01)
00237     {
00238         real delta=x-1.0;
00239         real log_delta=pl_log(delta);
00240         return Pi*Pi/6.0 + delta*(1-log_delta+delta*
00241                                   ((2*log_delta-1)/4 + delta*
00242                                    ((1-3*log_delta)/9 + delta*
00243                                     ((4*log_delta-1)/16 + delta*
00244                                      ((1-5*log_delta)/25 + delta*
00245                                       ((6*log_delta-1)/36 + delta*
00246                                        ((1-7*log_delta)/49 + delta*
00247                                         (8*log_delta-1)/64)))))));
00248     }
00249     else if (x<=2.0)
00250     {
00251         real logx = pl_log(x);
00252         return Pi*Pi/6.0 + small_dilogarithm(1.0-1.0/x) - logx*(0.5*logx+pl_log(1-1/x));
00253     } else 
00254     {
00255         real logx = pl_log(x);
00256         return Pi*Pi/3.0 - small_dilogarithm(1.0/x) - 0.5*logx*logx;
00257     }
00258 }
00259 
00260 real dilogarithm(real x)
00261 {
00262     if (is_missing(x))
00263     {
00264 #ifdef BOUNDCHECK
00265         PLWARNING("Dilogarithm taking NaN as input");
00266 #endif
00267         return MISSING_VALUE;
00268     }
00269     if (x<0)
00270         return -positive_dilogarithm(-x) + 0.5*positive_dilogarithm(x*x);
00271     else 
00272         if (fast_exact_is_equal(x, 0)) return 0;
00273         else
00274             return positive_dilogarithm(x);
00275 }
00276 
00277 real hard_slope_integral(real l, real r, real a, real b)
00278 {
00279     if (b<l) return 0;
00280     if (b<r)
00281     {
00282         if (a<l) 
00283             return 0.5*(b-l)*(b-l)/(r-l);
00284         else // a>=l
00285             return 0.5*((b-l)*(b-l)-(a-l)*(a-l))/(r-l);
00286     }
00287     else // b>=r
00288     {
00289         if (a<l)
00290             return 0.5*(r-l)+(b-r);
00291         else if (a<r) // l<a<r
00292             return 0.5*((r-l) - (a-l)*(a-l)/(r-l)) + (b-r);
00293         else // a>r
00294             return b-a;
00295     }
00296 }
00297 
00298 real soft_slope_integral(real smoothness, real left, real right, real a, real b)
00299 {
00300     if (fast_exact_is_equal(smoothness, 0))
00301         return 0.5*(b-a);
00302     if (smoothness<100)
00303         return 
00304             (b - a) + (softplus_primitive(-smoothness*(b-right)) - softplus_primitive(-smoothness*(b-left))
00305                        -softplus_primitive(-smoothness*(a-right)) + softplus_primitive(-smoothness*(a-left)))/
00306             (smoothness*smoothness*(right-left));
00307     // else do the integral of the hard slope function
00308     return hard_slope_integral(left,right,a,b);
00309 }
00310 
00311 real tabulated_soft_slope_integral(real smoothness, real left, real right, real a, real b)
00312 {
00313     if (fast_exact_is_equal(smoothness, 0))
00314         return 0.5*(b-a);
00315     if (smoothness<100)
00316         return 
00317             (b - a) + (tabulated_softplus_primitive(-smoothness*(b-right)) - tabulated_softplus_primitive(-smoothness*(b-left))
00318                        -tabulated_softplus_primitive(-smoothness*(a-right)) + tabulated_softplus_primitive(-smoothness*(a-left)))/
00319             (smoothness*smoothness*(right-left));
00320     // else do the integral of the hard slope function
00321     return hard_slope_integral(left,right,a,b);
00322 }
00323 
00324 } // end of namespace PLearn
00325 
00326 
00327 /*
00328   Local Variables:
00329   mode:c++
00330   c-basic-offset:4
00331   c-file-style:"stroustrup"
00332   c-file-offsets:((innamespace . 0)(inline-open . 0))
00333   indent-tabs-mode:nil
00334   fill-column:79
00335   End:
00336 */
00337 // 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