PLearn 0.1
|
00001 // -*- C++ -*- 00002 00003 // PRandom.cc 00004 // 00005 // Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio, University of Montreal 00006 // Copyright (C) 2005 Olivier Delalleau 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 * $Id: .pyskeleton_header 544 2003-09-01 00:05:31Z plearner $ 00038 ******************************************************* */ 00039 00040 // Authors: Olivier Delalleau 00041 00044 // Constants used for random numbers generation. 00045 #define RAND_EPS 1.2e-7 00046 #define RAND_RNMX (1.0 - RAND_EPS) 00047 00048 #include "PRandom.h" 00049 00050 namespace PLearn { 00051 using namespace std; 00052 00054 // PRandom // 00056 PRandom::PRandom(int32_t seed): 00057 #ifdef BOUNDCHECK 00058 samples_count(0), 00059 #endif 00060 exponential_distribution(0), 00061 normal_distribution(0), 00062 uniform_01(0), 00063 the_seed(0), 00064 fixed_seed(0), 00065 seed_(seed) 00066 { 00067 // For convenience, we systematically call build() in the constructor. 00068 build(); 00069 } 00070 00071 PRandom::PRandom(const PRandom& rhs): 00072 #ifdef BOUNDCHECK 00073 samples_count (rhs.get_samples_count()), 00074 #endif 00075 rgen (*(rhs.get_rgen())), 00076 the_seed (rhs.get_the_seed()), 00077 fixed_seed (rhs.get_fixed_seed()), 00078 seed_ (rhs.get_seed()) 00079 { 00080 // Note: the extra parentheses are here to tell the compiler that the 00081 // assignments are meant to be used as truth values. 00082 if ((exponential_distribution = rhs.get_exponential_distribution())) 00083 exponential_distribution = new boost::exponential_distribution<> 00084 (*exponential_distribution); 00085 if ((normal_distribution = rhs.get_normal_distribution())) 00086 normal_distribution = new boost::normal_distribution<> 00087 (*normal_distribution); 00088 if ((uniform_01 = rhs.get_uniform_01())) 00089 uniform_01 = new boost::uniform_01<boost::mt19937> 00090 (*uniform_01); 00091 } 00092 00093 PRandom PRandom::operator=(const PRandom& rhs) 00094 { 00095 #ifdef BOUNDCHECK 00096 samples_count = rhs.get_samples_count(); 00097 #endif 00098 rgen = *(rhs.get_rgen()); 00099 the_seed = rhs.get_the_seed(); 00100 fixed_seed = rhs.get_fixed_seed(); 00101 seed_ = rhs.get_seed(); 00102 00103 if ((exponential_distribution = rhs.get_exponential_distribution())) 00104 exponential_distribution = new boost::exponential_distribution<> 00105 (*exponential_distribution); 00106 if ((normal_distribution = rhs.get_normal_distribution())) 00107 normal_distribution = new boost::normal_distribution<> 00108 (*normal_distribution); 00109 if ((uniform_01 = rhs.get_uniform_01())) 00110 uniform_01 = new boost::uniform_01<boost::mt19937> 00111 (*uniform_01); 00112 00113 return (*this); 00114 } 00115 00116 00117 PLEARN_IMPLEMENT_OBJECT(PRandom, 00118 "Perform a number of random operations, including generating random numbers", 00119 "" 00120 ); 00121 00123 // declareOptions // 00125 void PRandom::declareOptions(OptionList& ol) 00126 { 00127 declareOption(ol, "seed", &PRandom::seed_, OptionBase::buildoption, 00128 "Seed for the random number generator, set at build time:\n" 00129 " - -1 : initialized with the current CPU time\n" 00130 " - 0 : the current seed is left intact\n" 00131 " - x > 0 : the seed is changed to 'x'"); 00132 00133 // Declared as a learnt option to hide some complexity to the novice. 00134 declareOption(ol, "fixed_seed", &PRandom::fixed_seed, OptionBase::learntoption, 00135 "If set to 0, will be ignored. If set to -2, its value will be copied from\n" 00136 "the 'seed' option. If set to any other value, it must always be equal to\n" 00137 "'seed' when build() is called. This allows one to prevent the seed from\n" 00138 "being accidentally modified by setting 'fixed_seed' to -2. Someone modifying\n" 00139 "the seed afterwards will then get an error.\n"); 00140 00141 // Now call the parent class' declareOptions 00142 inherited::declareOptions(ol); 00143 } 00144 00146 // bounded_uniform // 00148 real PRandom::bounded_uniform(real a, real b) { 00149 real res = uniform_sample()*(b-a) + a; 00150 if (res >= b) 00151 return b*RAND_RNMX; 00152 else 00153 return res; 00154 } 00155 00157 // build // 00159 void PRandom::build() 00160 { 00161 inherited::build(); 00162 build_(); 00163 } 00164 00166 // build_ // 00168 void PRandom::build_() 00169 { 00170 if (fixed_seed) { 00171 if (fixed_seed == -2) 00172 fixed_seed = seed_; 00173 else { 00174 if (seed_ != fixed_seed) 00175 PLERROR("In PRandom::build_ - You are not allowed to modify the seed of " 00176 "a PRandom object whose seed has been fixed"); 00177 } 00178 } 00179 if (seed_ == -1) 00180 this->time_seed_(); 00181 else if (seed_ == 0) {} 00182 else if (seed_ > 0) 00183 this->manual_seed_(seed_); 00184 else 00185 PLERROR("In PRandom::build_ - The only value allowed for the seed are " 00186 "-1, 0 or a strictly positive int32_t integer"); 00187 } 00188 00190 // common // 00192 PP<PRandom> PRandom::common(bool random_seed) 00193 { 00194 static PP<PRandom> gen_random = 0; 00195 static PP<PRandom> gen_const = 0; 00196 if (random_seed) { 00197 if (!gen_random) { 00198 gen_random = new PRandom(); 00199 gen_random->fixed_seed = -2; 00200 gen_random->build(); 00201 } 00202 return gen_random; 00203 } else { 00204 if (!gen_const) { 00205 gen_const = new PRandom(12345678); 00206 gen_const->fixed_seed = -2; 00207 gen_const->build(); 00208 } 00209 return gen_const; 00210 } 00211 } 00212 00214 // exp_sample // 00216 real PRandom::exp_sample() { 00217 ensure_exponential_distribution(); 00218 #ifdef BOUNDCHECK 00219 samples_count++; 00220 #endif 00221 return real((*exponential_distribution)(*uniform_01)); 00222 } 00223 00225 // fill_random_discrete // 00227 void PRandom::fill_random_discrete(const Vec& dest, const Vec& set) 00228 { 00229 PLASSERT( dest.isEmpty() || !set.isEmpty() ); 00230 Vec::iterator it = dest.begin(); 00231 Vec::iterator itend = dest.end(); 00232 int n = set.length(); 00233 for(; it != itend; ++it) 00234 *it = set[this->uniform_multinomial_sample(n)]; 00235 } 00236 00238 // fill_random_normal // 00240 void PRandom::fill_random_normal(const Vec& dest, real mean, real stddev) { 00241 for (int i = 0; i < dest.length(); i++) 00242 dest[i] = gaussian_mu_sigma(mean, stddev); 00243 } 00244 00245 void PRandom::fill_random_normal(const Mat& dest, real mean, real stddev) { 00246 for (int i = 0; i < dest.length(); i++) 00247 fill_random_normal(dest(i), mean, stddev); 00248 } 00249 00251 // fill_random_uniform // 00253 void PRandom::fill_random_uniform(const Vec& dest, real min, real max) 00254 { 00255 for (int i = 0; i < dest.length(); i++) 00256 dest[i] = bounded_uniform(min, max); 00257 } 00258 00259 void PRandom::fill_random_uniform(const Mat& dest, real min, real max) 00260 { 00261 for (int i = 0; i < dest.length(); i++) 00262 fill_random_uniform(dest(i), min, max); 00263 } 00264 00266 // gaussian_01 // 00268 real PRandom::gaussian_01() { 00269 ensure_normal_distribution(); 00270 #ifdef BOUNDCHECK 00271 // Drawing one Gaussian sample has the same effect on the underlying 00272 // generator as drawing two uniform or exponential samples 00273 samples_count += 2; 00274 #endif 00275 return real((*normal_distribution)(*uniform_01)); 00276 } 00277 00279 // gaussian_mu_sigma // 00281 real PRandom::gaussian_mu_sigma(real mu, real sigma) { 00282 return gaussian_01() * sigma + mu; 00283 } 00284 00286 // makeDeepCopyFromShallowCopy // 00288 void PRandom::makeDeepCopyFromShallowCopy(CopiesMap& copies) 00289 { 00290 inherited::makeDeepCopyFromShallowCopy(copies); 00291 // Nothing more should be added here: this object is meant to be properly 00292 // copied directly within the copy constructor (this is where you should add 00293 // any statement needed for a proper copy). 00294 } 00295 00297 // manual_seed // 00299 void PRandom::manual_seed(int32_t x) 00300 { 00301 if (fixed_seed) 00302 PLERROR("In PRandom::manual_seed - You are not allowed to change the seed " 00303 "of a PRandom object whose seed is fixed"); 00304 seed_ = x; 00305 build(); 00306 } 00307 00309 // manual_seed_ // 00311 void PRandom::manual_seed_(int32_t x) 00312 { 00313 the_seed = uint32_t(x); 00314 rgen.seed(the_seed); 00315 if (uniform_01) { 00316 // The boost::uniform_01 object must be re-constructed from the updated 00317 // random number generator. 00318 delete uniform_01; 00319 uniform_01 = 0; 00320 } 00321 // Systematically construct the uniform_01 member, which is the basis for most 00322 // of the random operations. 00323 ensure_uniform_01(); 00324 #ifdef BOUNDCHECK 00325 samples_count = 0; 00326 #endif 00327 } 00328 00330 // multinomial_sample // 00332 int PRandom::multinomial_sample(const Vec& distribution) { 00333 real u = this->uniform_sample(); 00334 real* pi = distribution.data(); 00335 real s = *pi; 00336 int n = distribution.length(); 00337 int i = 0; 00338 while ((i<n) && (s<u)) { 00339 i++; 00340 pi++; 00341 s += *pi; 00342 } 00343 if (i == n) 00344 i = n - 1; // Improbable, but... 00345 return i; 00346 } 00347 00349 // binomial_sample // 00351 int PRandom::binomial_sample(real pp) { 00352 if( pp < 0 || pp > 1 ) 00353 PLERROR("In PRandom::binomial_sample, pp should be between 0 and 1, " 00354 "but is %f.", pp); 00355 00356 real u = this->uniform_sample(); 00357 if( pp < u ) 00358 return 0; 00359 else 00360 return 1; 00361 } 00362 00364 // time_seed_ // 00366 void PRandom::time_seed_() 00367 { 00368 time_t ltime; 00369 struct tm *today; 00370 time(<ime); 00371 today = localtime(<ime); 00372 manual_seed_((int32_t)today->tm_sec+ 00373 60*today->tm_min+ 00374 60*60*today->tm_hour+ 00375 60*60*24*today->tm_mday); 00376 } 00377 00379 // uniform_sample // 00381 real PRandom::uniform_sample() { 00382 #ifdef BOUNDCHECK 00383 samples_count++; 00384 #endif 00385 return real((*uniform_01)()); 00386 } 00387 00388 00389 00391 // newRGen // 00393 PP<PRandom> PRandom::split() 00394 { 00395 // Draw twice from underlying rgen, 00396 // just in case it outputs its state 00397 // directly... (any better ideas?) 00398 int s; 00399 do 00400 { 00401 //no need for a specific distribution 00402 s= rgen() + rgen(); 00403 if(s < 0) s= -s; 00404 } while(s == 0); 00405 00406 return new PRandom(s); 00407 } 00408 00409 00411 // ~ // 00413 PRandom::~PRandom() { 00414 if (uniform_01) { 00415 delete uniform_01; 00416 uniform_01 = 0; 00417 } 00418 if (normal_distribution) { 00419 delete normal_distribution; 00420 normal_distribution = 0; 00421 } 00422 if (exponential_distribution) { 00423 delete exponential_distribution; 00424 exponential_distribution = 0; 00425 } 00426 } 00427 00428 } // end of namespace PLearn 00429 00430 00431 /* 00432 Local Variables: 00433 mode:c++ 00434 c-basic-offset:4 00435 c-file-style:"stroustrup" 00436 c-file-offsets:((innamespace . 0)(inline-open . 0)) 00437 indent-tabs-mode:nil 00438 fill-column:79 00439 End: 00440 */ 00441 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :