PLearn 0.1
PRandom.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // PRandom.h
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 
00045 #ifndef PRandom_INC
00046 #define PRandom_INC
00047 
00048 #include <plearn/base/Object.h>
00049 #include <plearn/math/TMat_maths.h>          // sum()
00050 #include <boost/random/exponential_distribution.hpp>
00051 #include <boost/random/mersenne_twister.hpp>
00052 #include <boost/random/normal_distribution.hpp>
00053 #include <boost/random/uniform_01.hpp>
00054 
00055 namespace PLearn {
00056 
00064 class PRandom: public Object
00065 {
00066 
00067 private:
00068 
00069     typedef Object inherited;
00070 
00071 protected:
00072 
00073 #ifdef BOUNDCHECK
00074 
00075     int samples_count;
00076 #endif
00077 
00079     boost::mt19937 rgen;
00080 
00082     boost::exponential_distribution<>* exponential_distribution;
00083 
00085     boost::normal_distribution<>* normal_distribution;
00086 
00088     boost::uniform_01<boost::mt19937>* uniform_01;
00089 
00091     uint32_t the_seed;
00092 
00093     // *********************
00094     // * protected options *
00095     // *********************
00096 
00097     int32_t fixed_seed;
00098 
00099 public:
00100 
00101     // ************************
00102     // * public build options *
00103     // ************************
00104 
00105     int32_t seed_;
00106 
00107     // ****************
00108     // * Constructors *
00109     // ****************
00110 
00112     PRandom(int32_t seed = -1);
00113 
00117     PRandom(const PRandom& rhs);
00118 
00121     PRandom operator=(const PRandom& rhs);
00122 
00124     virtual ~PRandom();
00125 
00126     // ******************
00127     // * Object methods *
00128     // ******************
00129 
00131     const boost::mt19937*                     get_rgen()                     const
00132     { return &rgen; }
00133     boost::exponential_distribution<>*        get_exponential_distribution() const
00134     { return exponential_distribution; }
00135     boost::normal_distribution<>*             get_normal_distribution()      const
00136     { return normal_distribution; }
00137     boost::uniform_01<boost::mt19937>*        get_uniform_01()               const
00138     { return uniform_01; }
00139 
00140     uint32_t get_the_seed()   const { return the_seed; }
00141     int32_t  get_fixed_seed() const { return fixed_seed; }
00142     int32_t  get_seed()       const { return seed_; }
00143 #ifdef BOUNDCHECK
00144     int get_samples_count()   const { return samples_count; }
00145 #endif
00146 
00147 private:
00148 
00150     void build_();
00151 
00152 protected:
00153 
00155     static void declareOptions(OptionList& ol);
00156 
00159     void time_seed_();
00160 
00163     void manual_seed_(int32_t x);
00164 
00168     inline void ensure_uniform_01() {
00169         if (!uniform_01)
00170             uniform_01 = new boost::uniform_01<boost::mt19937>(rgen);
00171     }
00172 
00176     inline void ensure_normal_distribution() {
00177         if (!normal_distribution)
00178             normal_distribution = new boost::normal_distribution<>();
00179     }
00180 
00184     inline void ensure_exponential_distribution() {
00185         if (!exponential_distribution)
00186             exponential_distribution = new boost::exponential_distribution<>();
00187     }
00188 
00189 public:
00190 
00191     // Declares other standard object methods.
00192     PLEARN_DECLARE_OBJECT(PRandom);
00193 
00199     virtual void build();
00200 
00202     virtual void makeDeepCopyFromShallowCopy(CopiesMap& copies);
00203 
00208     void manual_seed(int32_t x);
00209 
00211     inline void time_seed() { manual_seed(-1); }
00212 
00214     real uniform_sample();
00216     real bounded_uniform(real a, real b);
00217 
00219     real gaussian_01();
00220     inline real normal_sample() { return this->gaussian_01(); }
00222     real gaussian_mu_sigma(real mu, real sigma);
00223 
00225     void fill_random_discrete(const Vec& dest, const Vec& set);
00226 
00229     void fill_random_normal(const Vec& dest, real mean = 0, real stddev = 1);
00230 
00233     void fill_random_normal(const Mat& dest, real mean = 0, real stddev = 1);
00234 
00237     void fill_random_uniform(const Vec& dest, real min = 0, real max = 1);
00238 
00241     void fill_random_uniform(const Mat& dest, real min = 0, real max = 1);
00242 
00245     real exp_sample();
00246 
00247     /* TODO Implement.
00249     real gamma_sample(int ia);
00251     real poisson_sample(real xm);
00254     real binom_sample(real pp, int n = 1);
00255     */
00256 
00262     int multinomial_sample(const Vec& distribution);
00263 
00265     inline int uniform_multinomial_sample(int n)
00266     { return int(n * this->uniform_sample()); }
00267 
00269     int binomial_sample(real pp);
00270 
00272     template<class T>
00273     void shuffleElements(const TVec<T>& vec);
00274 
00281     template <class T>
00282     TVec<T> weightedShuffleElements(const TVec<T>& vec, const Vec& weights,
00283                                     int number_to_keep = -1 /* all by default */);
00284 
00286     template<class T>
00287     void shuffleRows(const TMat<T>& mat);
00288 
00291     PP<PRandom> split();
00292 
00293     /*** Static methods ***/
00294 
00303     static PP<PRandom> common(bool random_seed = true);
00304 };
00305 
00306 // Declares a few other classes and functions related to this class
00307 DECLARE_OBJECT_PTR(PRandom);
00308 
00309 
00310 //#####  shuffleElements Implementation  ######################################
00311 
00312 template<class T>
00313 void PRandom::shuffleElements(const TVec<T>& vec)
00314 {
00315     if (vec.isEmpty())
00316         return;
00317     T* v = vec.data();
00318     T tmp;
00319     int n = vec.length();
00320     for (int i = 0; i < vec.length(); i++) {
00321         int j = i + this->uniform_multinomial_sample(n - i);
00322         tmp = v[i];
00323         v[i] = v[j];
00324         v[j] = tmp;
00325     }
00326 }
00327 
00328 
00329 //#####  weightedShuffleElements Implementation  ##############################
00330 
00331 template <class T>
00332 TVec<T> PRandom::weightedShuffleElements(const TVec<T>& vec, const Vec& weights,
00333                                          int number_to_keep)
00334 {
00335     PLASSERT( vec.size() == weights.size() );
00336     if (number_to_keep < 0)
00337         number_to_keep = vec.size();
00338 
00339     Vec w = weights.copy();
00340     TVec<T> r(number_to_keep);
00341 
00342     // Normalize the weights
00343     real s = sum(w);
00344     w /= s;
00345 
00346     // The algorithm is fairly simple-minded and is O(N^2), so be careful: we
00347     // pick the first element according to a multinomial distribution given by
00348     // the weights.  Then we set the probability of that element to 0 and
00349     // renormalize the weights.  We pick the second element according to the
00350     // new distribution, and so forth.
00351     for (int i=0 ; i<number_to_keep ; ++i) {
00352         int j = this->multinomial_sample(w);
00353         r[i] = vec[j];
00354 
00355         // Set its weight to zero and renormalize.  To keep numerical
00356         // precision, every 20 time-steps we do a full renormalization (and not
00357         // only an incremental one).
00358         if (i % 20 == 0) {
00359             w[j] = 0;
00360             real s = sum(w);
00361             w /= s;
00362         }
00363         else {
00364             real old_w = w[j];
00365             w[j] = 0;
00366             w /= (1. - old_w); // New sum of weights is 1 minus what we took out
00367         }
00368     }
00369 
00370     return r;
00371 }
00372 
00373 
00374 //#####  shuffleRows Implementation  ##########################################
00375 
00376 template<class T>
00377 void PRandom::shuffleRows(const TMat<T>& mat)
00378 {
00379     int n = mat.length();
00380     for (int i = 0; i < n; i++) {
00381         int j = i + int(this->uniform_sample() * (n - i));
00382         mat.swapRows(i,j);
00383     }
00384 }
00385 
00386 
00387 } // end of namespace PLearn
00388 
00389 #endif
00390 
00391 
00392 /*
00393   Local Variables:
00394   mode:c++
00395   c-basic-offset:4
00396   c-file-style:"stroustrup"
00397   c-file-offsets:((innamespace . 0)(inline-open . 0))
00398   indent-tabs-mode:nil
00399   fill-column:79
00400   End:
00401 */
00402 // 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