PLearn 0.1
|
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 :