PLearn 0.1
PRandom.cc
Go to the documentation of this file.
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(&ltime);
00371     today = localtime(&ltime);
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 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines