IT++ Logo Newcom Logo

random.h

Go to the documentation of this file.
00001 
00033 #ifndef RANDOM_H
00034 #define RANDOM_H
00035 
00036 #include <itpp/base/operators.h>
00037 #include <ctime>
00038 
00039 
00040 namespace itpp {
00041 
00043 
00115   class Random_Generator {
00116   public:
00118     Random_Generator();
00120     void randomize();
00122     void reset() { initialize(lastSeed); reload(); }
00124     void reset(unsigned int seed) { lastSeed = seed; initialize(lastSeed); reload(); }
00126     unsigned int random_int()
00127     {
00128       if( left == 0 ) reload();
00129       --left;
00130                 
00131       register unsigned int s1;
00132       s1 = *pNext++;
00133       s1 ^= (s1 >> 11);
00134       s1 ^= (s1 <<  7) & 0x9d2c5680U;
00135       s1 ^= (s1 << 15) & 0xefc60000U;
00136       return ( s1 ^ (s1 >> 18) );
00137     }
00138 
00140     double random_01() { return (double(random_int()) + 0.5) * (1.0/4294967296.0); }
00142     double random_01_closed() { return double(random_int()) * (1.0/4294967295.0); }
00144     void get_state(ivec &out_state);
00146     void set_state(ivec &new_state);
00147   protected:
00148   private:
00150     unsigned int lastSeed;
00152     static unsigned int state[624];
00154     static unsigned int *pNext;
00156     static int left;
00157 
00159     void initialize( unsigned int seed )
00160     {
00161       // Initialize generator state with seed
00162       // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
00163       // In previous versions, most significant bits (MSBs) of the seed affect
00164       // only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
00165       register unsigned int *s = state;
00166       register unsigned int *r = state;
00167       register int i = 1;
00168       *s++ = seed & 0xffffffffU;
00169       for( ; i < 624; ++i ) 
00170         {
00171           *s++ = ( 1812433253U * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffU;
00172           r++;
00173         }
00174     }
00176     void reload()
00177     {
00178       // Generate N new values in state
00179       // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
00180       register unsigned int *p = state;
00181       register int i;
00182       for( i = 624 - 397; i--; ++p )
00183         *p = twist( p[397], p[0], p[1] );
00184       for( i = 397; --i; ++p )
00185         *p = twist( p[397-624], p[0], p[1] );
00186       *p = twist( p[397-624], p[0], state[0] );
00187       
00188       left = 624, pNext = state;
00189     }
00191     unsigned int hiBit( const unsigned int& u ) const { return u & 0x80000000U; }
00193     unsigned int loBit( const unsigned int& u ) const { return u & 0x00000001U; }
00195     unsigned int loBits( const unsigned int& u ) const { return u & 0x7fffffffU; }
00197     unsigned int mixBits( const unsigned int& u, const unsigned int& v ) const
00198     { return hiBit(u) | loBits(v); }
00199 
00200     /*
00201      * ----------------------------------------------------------------------
00202      * --- ediap - 2007/01/17 ---
00203      * ----------------------------------------------------------------------
00204      * Richard's implementation of the twist() function is was follows:
00205      *  { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfU); }
00206      * However, this code caused a warning/error under MSVC++, because
00207      * unsigned value loBit(s1) is being negated with `-' (c.f. bug report
00208      * [1635988]). I changed this to the same implementation as is used in
00209      * original C sources of Mersenne Twister RNG:
00210      *  #define MATRIX_A 0x9908b0dfUL
00211      *  #define UMASK 0x80000000UL
00212      *  #define LMASK 0x7fffffffUL
00213      *  #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
00214      *  #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) 
00215      * ----------------------------------------------------------------------
00216      */
00218     unsigned int twist( const unsigned int& m, const unsigned int& s0,
00219                         const unsigned int& s1 ) const
00220     { return m ^ (mixBits(s0,s1)>>1) ^ (loBit(s1) ? 0x9908b0dfU : 0U); }
00222     unsigned int hash( time_t t, clock_t c );
00223   };
00224 
00225 
00228 
00230   void RNG_reset(unsigned long seed);
00232   void RNG_reset();
00234   void RNG_randomize();
00236   void RNG_get_state(ivec &state);
00238   void RNG_set_state(ivec &state);
00240 
00245   class Bernoulli_RNG {
00246   public:
00248     Bernoulli_RNG(double prob) { setup(prob); }
00250     Bernoulli_RNG() { p=0.5; }
00252     void setup(double prob)
00253     {
00254       it_assert(prob>=0.0 && prob<=1.0, "The bernoulli source probability must be between 0 and 1");
00255       p = prob;
00256     }
00258     double get_setup() const { return p; }
00260     bin operator()() { return sample(); }
00262     bvec operator()(int n) { bvec temp(n); sample_vector(n, temp); return temp; }
00264     bmat operator()(int h, int w) { bmat temp(h, w); sample_matrix(h, w, temp); return temp; }
00266     bin sample() { return bin( RNG.random_01() < p ? 1 : 0 ); }
00268     void sample_vector(int size, bvec &out)
00269     {
00270       out.set_size(size, false);
00271       for (int i=0; i<size; i++) out(i) = sample(); 
00272     }
00274     void sample_matrix(int rows, int cols, bmat &out)
00275     {
00276       out.set_size(rows, cols, false);
00277       for (int i=0; i<rows*cols; i++) out(i) = sample();
00278     }
00279   protected:
00280   private:
00282     double p;
00284     Random_Generator RNG;
00285   };
00286 
00304   class I_Uniform_RNG {
00305   public:
00307     I_Uniform_RNG(int min=0, int max=1);
00309     void setup(int min, int max);
00311     void get_setup(int &min, int &max) const;
00313     int operator()() { return sample(); }
00315     ivec operator()(int n);
00317     imat operator()(int h, int w);
00319     int sample() { return ( floor_i(RNG.random_01() * (hi - lo + 1)) + lo ); }
00320   protected:
00321   private:
00323     int lo;
00325     int hi;
00327     Random_Generator RNG;
00328   };
00329 
00334   class Uniform_RNG {
00335   public:
00337     Uniform_RNG(double min=0, double max=1.0);
00339     void setup(double min, double max);
00341     void get_setup(double &min, double &max) const;
00343     double operator()() { return ( sample()* (hi_bound - lo_bound) + lo_bound ); }
00345     vec operator()(int n)
00346     { vec temp(n); sample_vector(n, temp); return (temp *(hi_bound - lo_bound) + lo_bound); }
00348     mat operator()(int h, int w)
00349     { mat temp(h,w); sample_matrix(h,w, temp); return (temp *(hi_bound - lo_bound) + lo_bound); }
00351     double sample() {  return RNG.random_01(); }
00353     void sample_vector(int size, vec &out)
00354     {
00355       out.set_size(size, false);
00356       for (int i=0; i<size; i++) out(i) = sample(); 
00357     }
00359     void sample_matrix(int rows, int cols, mat &out)
00360     {
00361       out.set_size(rows, cols, false);
00362       for (int i=0; i<rows*cols; i++) out(i) = sample(); 
00363     }
00364   protected:
00365   private:
00367     double lo_bound, hi_bound;
00369     Random_Generator RNG;
00370   };
00371 
00376   class Exponential_RNG {
00377   public:
00379     Exponential_RNG(double lambda=1.0);
00381     void setup(double lambda) { l=lambda; }
00383     double get_setup() const;
00385     double operator()() { return sample(); }
00387     vec operator()(int n);
00389     mat operator()(int h, int w);
00390   protected:
00391   private:
00393     double sample() {  return ( -std::log(RNG.random_01()) / l ); }
00395     double l;
00397     Random_Generator RNG;
00398   };
00399 
00404   class Normal_RNG {
00405   public:
00407     Normal_RNG(double meanval, double variance) { setup(meanval, variance); };
00409     Normal_RNG() { mean = 0.0; sigma = 1.0; odd = false; }
00411     void setup(double meanval, double variance)
00412     { mean = meanval; sigma = std::sqrt(variance); odd = false; }
00414     void get_setup(double &meanval, double &variance) const;
00416     double operator()() { return (sigma*sample()+mean); }
00418     vec operator()(int n) { vec temp(n); sample_vector(n, temp); return (sigma*temp+mean); }
00420     mat operator()(int h, int w) { mat temp(h,w); sample_matrix(h,w, temp); return (sigma*temp+mean); }
00422     double sample() // Box-Mueller, Polar version, See Knuth v2, 3rd ed, p122 
00423     {
00424       double x, r2;
00425       if (odd) {
00426         odd = 0;
00427         return y;
00428       }
00429       do { x = 2. * RNG.random_01() - 1.; y = 2. * RNG.random_01() - 1.; r2 = x*x + y*y; } while (r2 >= 1. || r2 < 1E-30);
00430       r2 = std::sqrt(std::log(r2)*(-2./r2)); x *= r2;  y *= r2; odd = 1;
00431       return x;
00432     }
00433  
00435     void sample_vector(int size, vec &out)
00436     {
00437       out.set_size(size, false);
00438       for (int i=0; i<size; i++) out(i) = sample(); 
00439     }
00440 
00442     void sample_matrix(int rows, int cols, mat &out)
00443     {
00444       out.set_size(rows, cols, false);
00445       for (int i=0; i<rows*cols; i++) out(i) = sample(); 
00446     }
00447   protected:
00448   private:
00450     double mean, sigma, y;
00452     int odd;
00454     Random_Generator RNG;
00455   };
00456 
00461   class Laplace_RNG {
00462   public:
00464     Laplace_RNG(double meanval=0.0, double variance=1.0);
00466     void setup(double meanval, double variance);
00468     void get_setup(double &meanval, double &variance) const;
00470     double operator()() { return sample(); }
00472     vec operator()(int n);
00474     mat operator()(int h, int w);
00476     double sample()
00477     {
00478       u=RNG.random_01();
00479       if(u<0.5)
00480         l=std::log(2.0*u);
00481       else
00482         l=-std::log(2.0*(1-u));
00483       
00484       l *= std::sqrt(var/2.0);
00485       l += mean;
00486       
00487       return l;
00488     }
00489   protected:
00490   private:
00492     double mean, var, u, l;
00494     Random_Generator RNG;
00495   };
00496 
00501   class Complex_Normal_RNG {
00502   public:
00504     Complex_Normal_RNG(std::complex<double> mean, double variance) { setup(mean, variance); }
00506     Complex_Normal_RNG() { m = 0.0; sigma=1.0; }
00508     void setup(std::complex<double> mean, double variance) { m = mean; sigma = std::sqrt(variance); }
00510     void get_setup(std::complex<double> &mean, double &variance) { mean = m; variance = sigma*sigma; }
00512     std::complex<double> operator()() { return sigma*sample()+m; }
00514     cvec operator()(int n) { cvec temp(n); sample_vector(n, temp); return (sigma*temp+m); }
00516     cmat operator()(int h, int w) { cmat temp(h,w); sample_matrix(h,w, temp); return (sigma*temp+m); }
00518     std::complex<double> sample()
00519     {
00520       double a = m_2pi * RNG.random_01();
00521       double b = std::sqrt(-std::log(RNG.random_01()));
00522       return std::complex<double>(b * std::cos(a), b * std::sin(a));
00523     }
00524 
00526     void sample_vector(int size, cvec &out)
00527     {
00528       out.set_size(size, false);
00529       for (int i=0; i<size; i++) out(i) = sample(); 
00530     }
00531 
00533     void sample_matrix(int rows, int cols, cmat &out)
00534     {
00535       out.set_size(rows, cols, false);
00536       for (int i=0; i<rows*cols; i++) out(i) = sample(); 
00537     }
00538   protected:
00539   private:
00541     std::complex<double> m;
00543     double sigma;
00545     Random_Generator RNG;
00546   };
00547 
00552   class AR1_Normal_RNG {
00553   public:
00555     AR1_Normal_RNG(double meanval=0.0, double variance=1.0, double rho=0.0);
00557     void setup(double meanval, double variance, double rho);
00559     void get_setup(double &meanval, double &variance, double &rho) const;
00561     void reset();
00563     double operator()() { return sample(); }
00565     vec operator()(int n);
00567     mat operator()(int h, int w);
00568 
00569   protected:
00570   private: 
00572     double sample();
00574     double my_mean, mem, r, factr, mean, var, r1, r2;
00576     bool odd;
00578     Random_Generator RNG;
00579   };
00580 
00585   typedef Normal_RNG Gauss_RNG;
00586 
00591   typedef AR1_Normal_RNG AR1_Gauss_RNG;
00592 
00597   class Weibull_RNG {
00598   public:
00600     Weibull_RNG(double lambda=1.0, double beta=1.0);
00601 
00603     void setup(double lambda, double beta);
00605     void get_setup(double &lambda, double &beta) { lambda=l; beta=b; }
00607     double operator()() { return sample(); }
00609     vec operator()(int n);
00611     mat operator()(int h, int w);
00612   protected:
00613   private:
00615     double sample();
00617     double l, b;
00619     double mean, var;
00621     Random_Generator RNG;
00622   };
00623 
00628   class Rayleigh_RNG {
00629   public:
00631     Rayleigh_RNG(double sigma=1.0);
00632 
00634     void setup(double sigma) { sig=sigma; }
00636     double get_setup() { return sig; }
00638     double operator()() { return sample(); }
00640     vec operator()(int n);
00642     mat operator()(int h, int w);
00643   protected:
00644   private:
00646     double sample();
00648     double sig;
00650     Random_Generator RNG;
00651   };
00652 
00657   class Rice_RNG {
00658   public:
00660     Rice_RNG(double sigma=1.0, double _s=1.0);
00662     void setup(double sigma, double _s) { sig=sigma; s=_s; }
00664     void get_setup(double &sigma, double &_s) { sigma=sig; _s=s; }
00666     double operator()() { return sample(); }
00668     vec operator()(int n);
00670     mat operator()(int h, int w);
00671   protected:
00672   private:
00674     double sample();
00676     double sig, s;
00678     Random_Generator RNG;
00679   };
00680 
00683 
00685   inline bin randb(void) { Bernoulli_RNG src; return src.sample(); }
00687   inline void randb(int size, bvec &out) { Bernoulli_RNG src; src.sample_vector(size, out); }
00689   inline bvec randb(int size) { bvec temp; randb(size, temp); return temp; }
00691   inline void randb(int rows, int cols, bmat &out) { Bernoulli_RNG src; src.sample_matrix(rows, cols, out); }
00693   inline bmat randb(int rows, int cols){ bmat temp; randb(rows, cols, temp); return temp; }
00694 
00696   inline double randu(void) { Uniform_RNG src; return src.sample(); }
00698   inline void randu(int size, vec &out) { Uniform_RNG src; src.sample_vector(size, out); }
00700   inline vec randu(int size){ vec temp; randu(size, temp); return temp; }
00702   inline void randu(int rows, int cols, mat &out) { Uniform_RNG src; src.sample_matrix(rows, cols, out); }
00704   inline mat randu(int rows, int cols) { mat temp; randu(rows, cols, temp); return temp; }
00705 
00707   inline int randi(int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(); }
00709   inline ivec randi(int size, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(size); }
00711   inline imat randi(int rows, int cols, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(rows, cols); }
00712 
00714   inline vec randray(int size, double sigma = 1.0) { Rayleigh_RNG src; src.setup(sigma); return src(size); }
00715 
00717   inline vec randrice(int size, double sigma = 1.0, double s = 1.0) { Rice_RNG src; src.setup(sigma, s); return src(size); }
00718 
00720   inline vec randexp(int size, double lambda = 1.0) { Exponential_RNG src; src.setup(lambda); return src(size); }
00721 
00723   inline double randn(void) { Normal_RNG src; return src.sample(); }
00725   inline void randn(int size, vec &out) { Normal_RNG src; src.sample_vector(size, out); }
00727   inline vec randn(int size) { vec temp; randn(size, temp); return temp; }
00729   inline void randn(int rows, int cols, mat &out)  { Normal_RNG src; src.sample_matrix(rows, cols, out); }
00731   inline mat randn(int rows, int cols){ mat temp; randn(rows, cols, temp); return temp; }
00732 
00737   inline std::complex<double> randn_c(void) { Complex_Normal_RNG src; return src.sample(); }
00739   inline void randn_c(int size, cvec &out)  { Complex_Normal_RNG src; src.sample_vector(size, out); }
00741   inline cvec randn_c(int size){ cvec temp; randn_c(size, temp); return temp; }
00743   inline void randn_c(int rows, int cols, cmat &out) { Complex_Normal_RNG src; src.sample_matrix(rows, cols, out); }
00745   inline cmat randn_c(int rows, int cols) { cmat temp; randn_c(rows, cols, temp); return temp; }
00746 
00748 
00749   // -------------------- INLINES ----------------------------------------------
00750 
00751   inline double AR1_Normal_RNG::sample()
00752   {
00753     double s;
00754 
00755     if (odd) {
00756       r1 = RNG.random_01();
00757       r2 = RNG.random_01();
00758       s = std::sqrt(factr * std::log(r2)) * std::cos(m_2pi * r1);
00759     } else
00760       s = std::sqrt(factr * std::log(r2)) * std::sin(m_2pi * r1);
00761 
00762     odd = !odd;
00763   
00764     mem = s + r * mem;
00765     s = mem + mean;
00766   
00767     return s;
00768   }
00769 
00770   inline double Weibull_RNG::sample()
00771   {
00772     return ( std::pow(-std::log(RNG.random_01()), 1.0/b) / l );
00773   }
00774 
00775   inline double Rayleigh_RNG::sample()
00776   {
00777     double r1, r2;
00778     double s1, s2, samp;
00779 
00780     r1 = RNG.random_01();
00781     r2 = RNG.random_01();
00782     s1 = std::sqrt(-2.0 * std::log(r2)) * std::cos(m_2pi * r1);
00783     s2 = std::sqrt(-2.0 * std::log(r2)) * std::sin(m_2pi * r1);
00784     // s1 and s2 are N(0,1) and independent
00785     samp = sig * std::sqrt(s1*s1 + s2*s2);
00786      
00787     return samp;
00788   }
00789 
00790   inline double Rice_RNG::sample()
00791   {
00792     double r1, r2;
00793     double s1, s2, samp;
00794     double m1 = 0.0;
00795     double m2 = std::sqrt(s*s - m1*m1);
00796 
00797     r1 = RNG.random_01();
00798     r2 = RNG.random_01();
00799     s1 = std::sqrt(-2.0 * std::log(r2)) * std::cos(m_2pi * r1);
00800     s2 = std::sqrt(-2.0 * std::log(r2)) * std::sin(m_2pi * r1);
00801     // s1 and s2 are N(0,1) and independent
00802     samp = std::sqrt((sig*s1+m1) * (sig*s1+m1) + (sig*s2+m2) * (sig*s2+m2));
00803       
00804     return samp;
00805   }
00806 
00807 } // namespace itpp
00808 
00809 #endif // #ifndef RANDOM_H
SourceForge Logo

Generated on Sat Aug 25 23:40:53 2007 for IT++ by Doxygen 1.5.2