IT++ Logo

random.h

Go to the documentation of this file.
00001 
00030 #ifndef RANDOM_H
00031 #define RANDOM_H
00032 
00033 #include <itpp/base/operators.h>
00034 #include <ctime>
00035 
00036 
00037 namespace itpp {
00038 
00040 
00112   class Random_Generator {
00113   public:
00115     Random_Generator() { if (!initialized) reset(4357U); }
00117     Random_Generator(unsigned int seed) { reset(seed); }
00119     void randomize() { reset(hash(time(0), clock())); }
00121     void reset() { initialize(last_seed); reload(); initialized = true; }
00123     void reset(unsigned int seed) { last_seed = seed; reset(); }
00124 
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 (random_int() + 0.5) * (1.0/4294967296.0); }
00142     double random_01_lclosed() { return random_int() * (1.0/4294967296.0); }
00144     double random_01_closed() { return random_int() * (1.0/4294967295.0); }
00146     double random53_01_lclosed()
00147     {
00148       return ((random_int() >> 5) * 67108864.0 + (random_int() >> 6))
00149         * (1.0/9007199254740992.0); // by Isaku Wada
00150     }
00151 
00153     void get_state(ivec &out_state);
00155     void set_state(ivec &new_state);
00156 
00157   private:
00159     static bool initialized;
00161     unsigned int last_seed;
00163     static unsigned int state[624];
00165     static unsigned int *pNext;
00167     static int left;
00168 
00176     void initialize( unsigned int seed )
00177     {
00178       register unsigned int *s = state;
00179       register unsigned int *r = state;
00180       register int i = 1;
00181       *s++ = seed & 0xffffffffU;
00182       for( ; i < 624; ++i )
00183         {
00184           *s++ = ( 1812433253U * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffU;
00185           r++;
00186         }
00187     }
00188 
00193     void reload()
00194     {
00195       register unsigned int *p = state;
00196       register int i;
00197       for( i = 624 - 397; i--; ++p )
00198         *p = twist( p[397], p[0], p[1] );
00199       for( i = 397; --i; ++p )
00200         *p = twist( p[397-624], p[0], p[1] );
00201       *p = twist( p[397-624], p[0], state[0] );
00202 
00203       left = 624, pNext = state;
00204     }
00206     unsigned int hiBit( const unsigned int& u ) const { return u & 0x80000000U; }
00208     unsigned int loBit( const unsigned int& u ) const { return u & 0x00000001U; }
00210     unsigned int loBits( const unsigned int& u ) const { return u & 0x7fffffffU; }
00212     unsigned int mixBits( const unsigned int& u, const unsigned int& v ) const
00213     { return hiBit(u) | loBits(v); }
00214 
00215     /*
00216      * ----------------------------------------------------------------------
00217      * --- ediap - 2007/01/17 ---
00218      * ----------------------------------------------------------------------
00219      * Wagners's implementation of the twist() function was as follows:
00220      *  { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfU); }
00221      * However, this code caused a warning/error under MSVC++, because
00222      * unsigned value loBit(s1) is being negated with `-' (c.f. bug report
00223      * [1635988]). I changed this to the same implementation as is used in
00224      * original C sources of Mersenne Twister RNG:
00225      *  #define MATRIX_A 0x9908b0dfUL
00226      *  #define UMASK 0x80000000UL
00227      *  #define LMASK 0x7fffffffUL
00228      *  #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
00229      *  #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
00230      * ----------------------------------------------------------------------
00231      */
00233     unsigned int twist( const unsigned int& m, const unsigned int& s0,
00234                         const unsigned int& s1 ) const
00235     { return m ^ (mixBits(s0,s1)>>1) ^ (loBit(s1) ? 0x9908b0dfU : 0U); }
00237     unsigned int hash( time_t t, clock_t c );
00238   };
00239 
00240 
00243 
00245   void RNG_reset(unsigned int seed);
00247   void RNG_reset();
00249   void RNG_randomize();
00251   void RNG_get_state(ivec &state);
00253   void RNG_set_state(ivec &state);
00255 
00260   class Bernoulli_RNG {
00261   public:
00263     Bernoulli_RNG(double prob) { setup(prob); }
00265     Bernoulli_RNG() { p=0.5; }
00267     void setup(double prob)
00268     {
00269       it_assert(prob>=0.0 && prob<=1.0, "The bernoulli source probability must be between 0 and 1");
00270       p = prob;
00271     }
00273     double get_setup() const { return p; }
00275     bin operator()() { return sample(); }
00277     bvec operator()(int n) { bvec temp(n); sample_vector(n, temp); return temp; }
00279     bmat operator()(int h, int w) { bmat temp(h, w); sample_matrix(h, w, temp); return temp; }
00281     bin sample() { return bin( RNG.random_01() < p ? 1 : 0 ); }
00283     void sample_vector(int size, bvec &out)
00284     {
00285       out.set_size(size, false);
00286       for (int i=0; i<size; i++) out(i) = sample();
00287     }
00289     void sample_matrix(int rows, int cols, bmat &out)
00290     {
00291       out.set_size(rows, cols, false);
00292       for (int i=0; i<rows*cols; i++) out(i) = sample();
00293     }
00294   protected:
00295   private:
00297     double p;
00299     Random_Generator RNG;
00300   };
00301 
00319   class I_Uniform_RNG {
00320   public:
00322     I_Uniform_RNG(int min=0, int max=1);
00324     void setup(int min, int max);
00326     void get_setup(int &min, int &max) const;
00328     int operator()() { return sample(); }
00330     ivec operator()(int n);
00332     imat operator()(int h, int w);
00334     int sample() { return ( floor_i(RNG.random_01() * (hi - lo + 1)) + lo ); }
00335   protected:
00336   private:
00338     int lo;
00340     int hi;
00342     Random_Generator RNG;
00343   };
00344 
00349   class Uniform_RNG {
00350   public:
00352     Uniform_RNG(double min=0, double max=1.0);
00354     void setup(double min, double max);
00356     void get_setup(double &min, double &max) const;
00358     double operator()() { return (sample() * (hi_bound - lo_bound) + lo_bound); }
00360     vec operator()(int n)
00361     {
00362       vec temp(n);
00363       sample_vector(n, temp);
00364       temp *= hi_bound - lo_bound;
00365       temp += lo_bound;
00366       return temp;
00367     }
00369     mat operator()(int h, int w)
00370     {
00371       mat temp(h, w);
00372       sample_matrix(h, w, temp);
00373       temp *= hi_bound - lo_bound;
00374       temp += lo_bound;
00375       return temp;
00376     }
00378     double sample() {  return RNG.random_01(); }
00380     void sample_vector(int size, vec &out)
00381     {
00382       out.set_size(size, false);
00383       for (int i=0; i<size; i++) out(i) = sample();
00384     }
00386     void sample_matrix(int rows, int cols, mat &out)
00387     {
00388       out.set_size(rows, cols, false);
00389       for (int i=0; i<rows*cols; i++) out(i) = sample();
00390     }
00391   protected:
00392   private:
00394     double lo_bound, hi_bound;
00396     Random_Generator RNG;
00397   };
00398 
00403   class Exponential_RNG {
00404   public:
00406     Exponential_RNG(double lambda=1.0);
00408     void setup(double lambda) { l=lambda; }
00410     double get_setup() const;
00412     double operator()() { return sample(); }
00414     vec operator()(int n);
00416     mat operator()(int h, int w);
00417   protected:
00418   private:
00420     double sample() {  return ( -std::log(RNG.random_01()) / l ); }
00422     double l;
00424     Random_Generator RNG;
00425   };
00426 
00441   class Normal_RNG {
00442   public:
00444     Normal_RNG(double meanval, double variance):
00445       mean(meanval), sigma(std::sqrt(variance)) {}
00447     Normal_RNG(): mean(0.0), sigma(1.0) {}
00449     void setup(double meanval, double variance)
00450     { mean = meanval; sigma = std::sqrt(variance); }
00452     void get_setup(double &meanval, double &variance) const;
00454     double operator()() { return (sigma*sample()+mean); }
00456     vec operator()(int n)
00457     {
00458       vec temp(n);
00459       sample_vector(n, temp);
00460       temp *= sigma;
00461       temp += mean;
00462       return temp;
00463     }
00465     mat operator()(int h, int w)
00466     {
00467       mat temp(h, w);
00468       sample_matrix(h, w, temp);
00469       temp *= sigma;
00470       temp += mean;
00471       return temp;
00472     }
00474     double sample();
00475 
00477     void sample_vector(int size, vec &out)
00478     {
00479       out.set_size(size, false);
00480       for (int i=0; i<size; i++) out(i) = sample();
00481     }
00482 
00484     void sample_matrix(int rows, int cols, mat &out)
00485     {
00486       out.set_size(rows, cols, false);
00487       for (int i=0; i<rows*cols; i++) out(i) = sample();
00488     }
00489   private:
00490     double mean, sigma;
00491     static const double ytab[128];
00492     static const unsigned int ktab[128];
00493     static const double wtab[128];
00494     static const double PARAM_R;
00495     Random_Generator RNG;
00496   };
00497 
00502   class Laplace_RNG {
00503   public:
00505     Laplace_RNG(double meanval = 0.0, double variance = 1.0);
00507     void setup(double meanval, double variance);
00509     void get_setup(double &meanval, double &variance) const;
00511     double operator()() { return sample(); }
00513     vec operator()(int n);
00515     mat operator()(int h, int w);
00517     double sample()
00518     {
00519       double u = RNG.random_01();
00520       double l = sqrt_12var;
00521       if (u < 0.5)
00522         l *= std::log(2.0 * u);
00523       else
00524         l *= -std::log(2.0 * (1-u));
00525       return (l + mean);
00526     }
00527   private:
00528     double mean, var, sqrt_12var;
00529     Random_Generator RNG;
00530   };
00531 
00536   class Complex_Normal_RNG {
00537   public:
00539     Complex_Normal_RNG(std::complex<double> mean, double variance):
00540       norm_factor(1.0/std::sqrt(2.0))
00541     {
00542       setup(mean, variance);
00543     }
00545     Complex_Normal_RNG(): m(0.0), sigma(1.0), norm_factor(1.0/std::sqrt(2.0)) {}
00547     void setup(std::complex<double> mean, double variance)
00548     {
00549       m = mean; sigma = std::sqrt(variance);
00550     }
00552     void get_setup(std::complex<double> &mean, double &variance)
00553     {
00554       mean = m; variance = sigma*sigma;
00555     }
00557     std::complex<double> operator()() { return sigma*sample()+m; }
00559     cvec operator()(int n)
00560     {
00561       cvec temp(n);
00562       sample_vector(n, temp);
00563       temp *= sigma;
00564       temp += m;
00565       return temp;
00566     }
00568     cmat operator()(int h, int w)
00569     {
00570       cmat temp(h, w);
00571       sample_matrix(h, w, temp);
00572       temp *= sigma;
00573       temp += m;
00574       return temp;
00575     }
00577     std::complex<double> sample()
00578     {
00579       double a = nRNG.sample() * norm_factor;
00580       double b = nRNG.sample() * norm_factor;
00581       return std::complex<double>(a, b);
00582     }
00583 
00585     void sample_vector(int size, cvec &out)
00586     {
00587       out.set_size(size, false);
00588       for (int i=0; i<size; i++) out(i) = sample();
00589     }
00590 
00592     void sample_matrix(int rows, int cols, cmat &out)
00593     {
00594       out.set_size(rows, cols, false);
00595       for (int i=0; i<rows*cols; i++) out(i) = sample();
00596     }
00597   private:
00598     std::complex<double> m;
00599     double sigma;
00600     const double norm_factor;
00601     Normal_RNG nRNG;
00602   };
00603 
00608   class AR1_Normal_RNG {
00609   public:
00611     AR1_Normal_RNG(double meanval = 0.0, double variance = 1.0,
00612                    double rho = 0.0);
00614     void setup(double meanval, double variance, double rho);
00616     void get_setup(double &meanval, double &variance, double &rho) const;
00618     void reset();
00620     double operator()() { return sample(); }
00622     vec operator()(int n);
00624     mat operator()(int h, int w);
00625   private:
00626     double sample()
00627     {
00628       mem *= r;
00629       if (odd) {
00630         r1 = m_2pi * RNG.random_01();
00631         r2 = std::sqrt(factr * std::log(RNG.random_01()));
00632         mem += r2 * std::cos(r1);
00633       }
00634       else {
00635         mem += r2 * std::sin(r1);
00636       }
00637       odd = !odd;
00638       return (mem + mean);
00639     }
00640     double mem, r, factr, mean, var, r1, r2;
00641     bool odd;
00642     Random_Generator RNG;
00643   };
00644 
00649   typedef Normal_RNG Gauss_RNG;
00650 
00655   typedef AR1_Normal_RNG AR1_Gauss_RNG;
00656 
00661   class Weibull_RNG {
00662   public:
00664     Weibull_RNG(double lambda = 1.0, double beta = 1.0);
00666     void setup(double lambda, double beta);
00668     void get_setup(double &lambda, double &beta) { lambda = l; beta = b; }
00670     double operator()() { return sample(); }
00672     vec operator()(int n);
00674     mat operator()(int h, int w);
00675   private:
00676     double sample()
00677     {
00678       return (std::pow(-std::log(RNG.random_01()), 1.0/b) / l);
00679     }
00680     double l, b;
00681     double mean, var;
00682     Random_Generator RNG;
00683   };
00684 
00689   class Rayleigh_RNG {
00690   public:
00692     Rayleigh_RNG(double sigma = 1.0);
00694     void setup(double sigma) { sig = sigma; }
00696     double get_setup() { return sig; }
00698     double operator()() { return sample(); }
00700     vec operator()(int n);
00702     mat operator()(int h, int w);
00703   private:
00704     double sample()
00705     {
00706       double s1 = nRNG.sample();
00707       double s2 = nRNG.sample();
00708       // s1 and s2 are N(0,1) and independent
00709       return (sig * std::sqrt(s1*s1 + s2*s2));
00710     }
00711     double sig;
00712     Normal_RNG nRNG;
00713   };
00714 
00719   class Rice_RNG {
00720   public:
00722     Rice_RNG(double sigma = 1.0, double v = 1.0);
00724     void setup(double sigma, double v) { sig = sigma; s = v; }
00726     void get_setup(double &sigma, double &v) { sigma = sig; v = s; }
00728     double operator()() { return sample(); }
00730     vec operator()(int n);
00732     mat operator()(int h, int w);
00733   private:
00734     double sample()
00735     {
00736       double s1 = nRNG.sample() + s;
00737       double s2 = nRNG.sample();
00738       // s1 and s2 are N(0,1) and independent
00739       return (sig * std::sqrt(s1*s1 + s2*s2));
00740     }
00741     double sig, s;
00742     Normal_RNG nRNG;
00743   };
00744 
00747 
00749   inline bin randb(void) { Bernoulli_RNG src; return src.sample(); }
00751   inline void randb(int size, bvec &out) { Bernoulli_RNG src; src.sample_vector(size, out); }
00753   inline bvec randb(int size) { bvec temp; randb(size, temp); return temp; }
00755   inline void randb(int rows, int cols, bmat &out) { Bernoulli_RNG src; src.sample_matrix(rows, cols, out); }
00757   inline bmat randb(int rows, int cols){ bmat temp; randb(rows, cols, temp); return temp; }
00758 
00760   inline double randu(void) { Uniform_RNG src; return src.sample(); }
00762   inline void randu(int size, vec &out) { Uniform_RNG src; src.sample_vector(size, out); }
00764   inline vec randu(int size){ vec temp; randu(size, temp); return temp; }
00766   inline void randu(int rows, int cols, mat &out) { Uniform_RNG src; src.sample_matrix(rows, cols, out); }
00768   inline mat randu(int rows, int cols) { mat temp; randu(rows, cols, temp); return temp; }
00769 
00771   inline int randi(int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(); }
00773   inline ivec randi(int size, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(size); }
00775   inline imat randi(int rows, int cols, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(rows, cols); }
00776 
00778   inline vec randray(int size, double sigma = 1.0) { Rayleigh_RNG src; src.setup(sigma); return src(size); }
00779 
00781   inline vec randrice(int size, double sigma = 1.0, double s = 1.0) { Rice_RNG src; src.setup(sigma, s); return src(size); }
00782 
00784   inline vec randexp(int size, double lambda = 1.0) { Exponential_RNG src; src.setup(lambda); return src(size); }
00785 
00787   inline double randn(void) { Normal_RNG src; return src.sample(); }
00789   inline void randn(int size, vec &out) { Normal_RNG src; src.sample_vector(size, out); }
00791   inline vec randn(int size) { vec temp; randn(size, temp); return temp; }
00793   inline void randn(int rows, int cols, mat &out)  { Normal_RNG src; src.sample_matrix(rows, cols, out); }
00795   inline mat randn(int rows, int cols){ mat temp; randn(rows, cols, temp); return temp; }
00796 
00801   inline std::complex<double> randn_c(void) { Complex_Normal_RNG src; return src.sample(); }
00803   inline void randn_c(int size, cvec &out)  { Complex_Normal_RNG src; src.sample_vector(size, out); }
00805   inline cvec randn_c(int size){ cvec temp; randn_c(size, temp); return temp; }
00807   inline void randn_c(int rows, int cols, cmat &out) { Complex_Normal_RNG src; src.sample_matrix(rows, cols, out); }
00809   inline cmat randn_c(int rows, int cols) { cmat temp; randn_c(rows, cols, temp); return temp; }
00810 
00812 
00813 } // namespace itpp
00814 
00815 #endif // #ifndef RANDOM_H
SourceForge Logo

Generated on Sat Apr 19 10:41:12 2008 for IT++ by Doxygen 1.5.5