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     {
00476       unsigned int u, sign, i, j;
00477       double x, y;
00478 
00479       while(true) {
00480   u = RNG.random_int();
00481   sign = u & 0x80; // 1 bit for the sign
00482   i = u & 0x7f; // 7 bits to choose the step
00483   j = u >> 8; // 24 bits for the x-value
00484 
00485   x = j * wtab[i];
00486 
00487   if (j < ktab[i])
00488     break;
00489 
00490   if (i < 127) {
00491     y = ytab[i+1] + (ytab[i] - ytab[i+1]) * RNG.random_01();
00492   }
00493   else {
00494     x = PARAM_R - std::log(1.0 - RNG.random_01()) / PARAM_R;
00495     y = std::exp(-PARAM_R * (x - 0.5 * PARAM_R)) * RNG.random_01();
00496   }
00497 
00498   if (y < std::exp(-0.5 * x * x))
00499     break;
00500       }
00501       return sign ? x : -x;
00502     }
00503 
00505     void sample_vector(int size, vec &out)
00506     {
00507       out.set_size(size, false);
00508       for (int i=0; i<size; i++) out(i) = sample();
00509     }
00510 
00512     void sample_matrix(int rows, int cols, mat &out)
00513     {
00514       out.set_size(rows, cols, false);
00515       for (int i=0; i<rows*cols; i++) out(i) = sample();
00516     }
00517   private:
00518     double mean, sigma;
00519     static const double ytab[128];
00520     static const unsigned int ktab[128];
00521     static const double wtab[128];
00522     static const double PARAM_R;
00523     Random_Generator RNG;
00524   };
00525 
00530   class Laplace_RNG {
00531   public:
00533     Laplace_RNG(double meanval = 0.0, double variance = 1.0);
00535     void setup(double meanval, double variance);
00537     void get_setup(double &meanval, double &variance) const;
00539     double operator()() { return sample(); }
00541     vec operator()(int n);
00543     mat operator()(int h, int w);
00545     double sample()
00546     {
00547       double u = RNG.random_01();
00548       double l = sqrt_12var;
00549       if (u < 0.5)
00550   l *= std::log(2.0 * u);
00551       else
00552   l *= -std::log(2.0 * (1-u));
00553       return (l + mean);
00554     }
00555   private:
00556     double mean, var, sqrt_12var;
00557     Random_Generator RNG;
00558   };
00559 
00564   class Complex_Normal_RNG {
00565   public:
00567     Complex_Normal_RNG(std::complex<double> mean, double variance):
00568       norm_factor(1.0/std::sqrt(2.0))
00569     {
00570       setup(mean, variance);
00571     }
00573     Complex_Normal_RNG(): m(0.0), sigma(1.0), norm_factor(1.0/std::sqrt(2.0)) {}
00575     void setup(std::complex<double> mean, double variance)
00576     {
00577       m = mean; sigma = std::sqrt(variance);
00578     }
00580     void get_setup(std::complex<double> &mean, double &variance)
00581     {
00582       mean = m; variance = sigma*sigma;
00583     }
00585     std::complex<double> operator()() { return sigma*sample()+m; }
00587     cvec operator()(int n)
00588     {
00589       cvec temp(n);
00590       sample_vector(n, temp);
00591       temp *= sigma;
00592       temp += m;
00593       return temp;
00594     }
00596     cmat operator()(int h, int w)
00597     {
00598       cmat temp(h, w);
00599       sample_matrix(h, w, temp);
00600       temp *= sigma;
00601       temp += m;
00602       return temp;
00603     }
00605     std::complex<double> sample()
00606     {
00607       double a = nRNG.sample() * norm_factor;
00608       double b = nRNG.sample() * norm_factor;
00609       return std::complex<double>(a, b);
00610     }
00611 
00613     void sample_vector(int size, cvec &out)
00614     {
00615       out.set_size(size, false);
00616       for (int i=0; i<size; i++) out(i) = sample();
00617     }
00618 
00620     void sample_matrix(int rows, int cols, cmat &out)
00621     {
00622       out.set_size(rows, cols, false);
00623       for (int i=0; i<rows*cols; i++) out(i) = sample();
00624     }
00625   private:
00626     std::complex<double> m;
00627     double sigma;
00628     const double norm_factor;
00629     Normal_RNG nRNG;
00630   };
00631 
00636   class AR1_Normal_RNG {
00637   public:
00639     AR1_Normal_RNG(double meanval = 0.0, double variance = 1.0,
00640        double rho = 0.0);
00642     void setup(double meanval, double variance, double rho);
00644     void get_setup(double &meanval, double &variance, double &rho) const;
00646     void reset();
00648     double operator()() { return sample(); }
00650     vec operator()(int n);
00652     mat operator()(int h, int w);
00653   private:
00654     double sample()
00655     {
00656       mem *= r;
00657       if (odd) {
00658   r1 = m_2pi * RNG.random_01();
00659   r2 = std::sqrt(factr * std::log(RNG.random_01()));
00660   mem += r2 * std::cos(r1);
00661       }
00662       else {
00663   mem += r2 * std::sin(r1);
00664       }
00665       odd = !odd;
00666       return (mem + mean);
00667     }
00668     double mem, r, factr, mean, var, r1, r2;
00669     bool odd;
00670     Random_Generator RNG;
00671   };
00672 
00677   typedef Normal_RNG Gauss_RNG;
00678 
00683   typedef AR1_Normal_RNG AR1_Gauss_RNG;
00684 
00689   class Weibull_RNG {
00690   public:
00692     Weibull_RNG(double lambda = 1.0, double beta = 1.0);
00694     void setup(double lambda, double beta);
00696     void get_setup(double &lambda, double &beta) { lambda = l; beta = b; }
00698     double operator()() { return sample(); }
00700     vec operator()(int n);
00702     mat operator()(int h, int w);
00703   private:
00704     double sample()
00705     {
00706       return (std::pow(-std::log(RNG.random_01()), 1.0/b) / l);
00707     }
00708     double l, b;
00709     double mean, var;
00710     Random_Generator RNG;
00711   };
00712 
00717   class Rayleigh_RNG {
00718   public:
00720     Rayleigh_RNG(double sigma = 1.0);
00722     void setup(double sigma) { sig = sigma; }
00724     double get_setup() { return sig; }
00726     double operator()() { return sample(); }
00728     vec operator()(int n);
00730     mat operator()(int h, int w);
00731   private:
00732     double sample()
00733     {
00734       double s1 = nRNG.sample();
00735       double s2 = nRNG.sample();
00736       // s1 and s2 are N(0,1) and independent
00737       return (sig * std::sqrt(s1*s1 + s2*s2));
00738     }
00739     double sig;
00740     Normal_RNG nRNG;
00741   };
00742 
00747   class Rice_RNG {
00748   public:
00750     Rice_RNG(double sigma = 1.0, double v = 1.0);
00752     void setup(double sigma, double v) { sig = sigma; s = v; }
00754     void get_setup(double &sigma, double &v) { sigma = sig; v = s; }
00756     double operator()() { return sample(); }
00758     vec operator()(int n);
00760     mat operator()(int h, int w);
00761   private:
00762     double sample()
00763     {
00764       double s1 = nRNG.sample() + s;
00765       double s2 = nRNG.sample();
00766       // s1 and s2 are N(0,1) and independent
00767       return (sig * std::sqrt(s1*s1 + s2*s2));
00768     }
00769     double sig, s;
00770     Normal_RNG nRNG;
00771   };
00772 
00775 
00777   inline bin randb(void) { Bernoulli_RNG src; return src.sample(); }
00779   inline void randb(int size, bvec &out) { Bernoulli_RNG src; src.sample_vector(size, out); }
00781   inline bvec randb(int size) { bvec temp; randb(size, temp); return temp; }
00783   inline void randb(int rows, int cols, bmat &out) { Bernoulli_RNG src; src.sample_matrix(rows, cols, out); }
00785   inline bmat randb(int rows, int cols){ bmat temp; randb(rows, cols, temp); return temp; }
00786 
00788   inline double randu(void) { Uniform_RNG src; return src.sample(); }
00790   inline void randu(int size, vec &out) { Uniform_RNG src; src.sample_vector(size, out); }
00792   inline vec randu(int size){ vec temp; randu(size, temp); return temp; }
00794   inline void randu(int rows, int cols, mat &out) { Uniform_RNG src; src.sample_matrix(rows, cols, out); }
00796   inline mat randu(int rows, int cols) { mat temp; randu(rows, cols, temp); return temp; }
00797 
00799   inline int randi(int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(); }
00801   inline ivec randi(int size, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(size); }
00803   inline imat randi(int rows, int cols, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(rows, cols); }
00804 
00806   inline vec randray(int size, double sigma = 1.0) { Rayleigh_RNG src; src.setup(sigma); return src(size); }
00807 
00809   inline vec randrice(int size, double sigma = 1.0, double s = 1.0) { Rice_RNG src; src.setup(sigma, s); return src(size); }
00810 
00812   inline vec randexp(int size, double lambda = 1.0) { Exponential_RNG src; src.setup(lambda); return src(size); }
00813 
00815   inline double randn(void) { Normal_RNG src; return src.sample(); }
00817   inline void randn(int size, vec &out) { Normal_RNG src; src.sample_vector(size, out); }
00819   inline vec randn(int size) { vec temp; randn(size, temp); return temp; }
00821   inline void randn(int rows, int cols, mat &out)  { Normal_RNG src; src.sample_matrix(rows, cols, out); }
00823   inline mat randn(int rows, int cols){ mat temp; randn(rows, cols, temp); return temp; }
00824 
00829   inline std::complex<double> randn_c(void) { Complex_Normal_RNG src; return src.sample(); }
00831   inline void randn_c(int size, cvec &out)  { Complex_Normal_RNG src; src.sample_vector(size, out); }
00833   inline cvec randn_c(int size){ cvec temp; randn_c(size, temp); return temp; }
00835   inline void randn_c(int rows, int cols, cmat &out) { Complex_Normal_RNG src; src.sample_matrix(rows, cols, out); }
00837   inline cmat randn_c(int rows, int cols) { cmat temp; randn_c(rows, cols, temp); return temp; }
00838 
00840 
00841 } // namespace itpp
00842 
00843 #endif // #ifndef RANDOM_H
SourceForge Logo

Generated on Sun Sep 14 18:52:26 2008 for IT++ by Doxygen 1.5.6