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.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 :