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 {
00039 
00041 
00113 class Random_Generator
00114 {
00115 public:
00117   Random_Generator() { if (!initialized) reset(4357U); }
00119   Random_Generator(unsigned int seed) { reset(seed); }
00121   void randomize() { reset(hash(time(0), clock())); }
00123   void reset() { initialize(last_seed); reload(); initialized = true; }
00125   void reset(unsigned int seed) { last_seed = seed; reset(); }
00126 
00128   unsigned int random_int() {
00129     if (left == 0) reload();
00130     --left;
00131 
00132     register unsigned int s1;
00133     s1 = *pNext++;
00134     s1 ^= (s1 >> 11);
00135     s1 ^= (s1 <<  7) & 0x9d2c5680U;
00136     s1 ^= (s1 << 15) & 0xefc60000U;
00137     return (s1 ^(s1 >> 18));
00138   }
00139 
00141   double random_01() { return (random_int() + 0.5) * (1.0 / 4294967296.0); }
00143   double random_01_lclosed() { return random_int() * (1.0 / 4294967296.0); }
00145   double random_01_closed() { return random_int() * (1.0 / 4294967295.0); }
00147   double random53_01_lclosed() {
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     register unsigned int *s = state;
00178     register unsigned int *r = state;
00179     register int i = 1;
00180     *s++ = seed & 0xffffffffU;
00181     for (; i < 624; ++i) {
00182       *s++ = (1812433253U * (*r ^(*r >> 30)) + i) & 0xffffffffU;
00183       r++;
00184     }
00185   }
00186 
00191   void reload() {
00192     register unsigned int *p = state;
00193     register int i;
00194     for (i = 624 - 397; i--; ++p)
00195       *p = twist(p[397], p[0], p[1]);
00196     for (i = 397; --i; ++p)
00197       *p = twist(p[397-624], p[0], p[1]);
00198     *p = twist(p[397-624], p[0], state[0]);
00199 
00200     left = 624, pNext = state;
00201   }
00203   unsigned int hiBit(const unsigned int& u) const { return u & 0x80000000U; }
00205   unsigned int loBit(const unsigned int& u) const { return u & 0x00000001U; }
00207   unsigned int loBits(const unsigned int& u) const { return u & 0x7fffffffU; }
00209   unsigned int mixBits(const unsigned int& u, const unsigned int& v) const
00210   { return hiBit(u) | loBits(v); }
00211 
00212   /*
00213    * ----------------------------------------------------------------------
00214    * --- ediap - 2007/01/17 ---
00215    * ----------------------------------------------------------------------
00216    * Wagners's implementation of the twist() function was as follows:
00217    *  { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfU); }
00218    * However, this code caused a warning/error under MSVC++, because
00219    * unsigned value loBit(s1) is being negated with `-' (c.f. bug report
00220    * [1635988]). I changed this to the same implementation as is used in
00221    * original C sources of Mersenne Twister RNG:
00222    *  #define MATRIX_A 0x9908b0dfUL
00223    *  #define UMASK 0x80000000UL
00224    *  #define LMASK 0x7fffffffUL
00225    *  #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
00226    *  #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
00227    * ----------------------------------------------------------------------
00228    */
00230   unsigned int twist(const unsigned int& m, const unsigned int& s0,
00231                      const unsigned int& s1) const
00232   { return m ^(mixBits(s0, s1) >> 1) ^(loBit(s1) ? 0x9908b0dfU : 0U); }
00234   unsigned int hash(time_t t, clock_t c);
00235 };
00236 
00237 
00240 
00242 void RNG_reset(unsigned int seed);
00244 void RNG_reset();
00246 void RNG_randomize();
00248 void RNG_get_state(ivec &state);
00250 void RNG_set_state(ivec &state);
00252 
00257 class Bernoulli_RNG
00258 {
00259 public:
00261   Bernoulli_RNG(double prob) { setup(prob); }
00263   Bernoulli_RNG() { p = 0.5; }
00265   void setup(double prob) {
00266     it_assert(prob >= 0.0 && prob <= 1.0, "The bernoulli source probability must be between 0 and 1");
00267     p = prob;
00268   }
00270   double get_setup() const { return p; }
00272   bin operator()() { return sample(); }
00274   bvec operator()(int n) { bvec temp(n); sample_vector(n, temp); return temp; }
00276   bmat operator()(int h, int w) { bmat temp(h, w); sample_matrix(h, w, temp); return temp; }
00278   bin sample() { return bin(RNG.random_01() < p ? 1 : 0); }
00280   void sample_vector(int size, bvec &out) {
00281     out.set_size(size, false);
00282     for (int i = 0; i < size; i++) out(i) = sample();
00283   }
00285   void sample_matrix(int rows, int cols, bmat &out) {
00286     out.set_size(rows, cols, false);
00287     for (int i = 0; i < rows*cols; i++) out(i) = sample();
00288   }
00289 protected:
00290 private:
00292   double p;
00294   Random_Generator RNG;
00295 };
00296 
00314 class I_Uniform_RNG
00315 {
00316 public:
00318   I_Uniform_RNG(int min = 0, int max = 1);
00320   void setup(int min, int max);
00322   void get_setup(int &min, int &max) const;
00324   int operator()() { return sample(); }
00326   ivec operator()(int n);
00328   imat operator()(int h, int w);
00330   int sample() { return (floor_i(RNG.random_01() * (hi - lo + 1)) + lo); }
00331 protected:
00332 private:
00334   int lo;
00336   int hi;
00338   Random_Generator RNG;
00339 };
00340 
00345 class Uniform_RNG
00346 {
00347 public:
00349   Uniform_RNG(double min = 0, double max = 1.0);
00351   void setup(double min, double max);
00353   void get_setup(double &min, double &max) const;
00355   double operator()() { return (sample() * (hi_bound - lo_bound) + lo_bound); }
00357   vec operator()(int n) {
00358     vec temp(n);
00359     sample_vector(n, temp);
00360     temp *= hi_bound - lo_bound;
00361     temp += lo_bound;
00362     return temp;
00363   }
00365   mat operator()(int h, int w) {
00366     mat temp(h, w);
00367     sample_matrix(h, w, temp);
00368     temp *= hi_bound - lo_bound;
00369     temp += lo_bound;
00370     return temp;
00371   }
00373   double sample() {  return RNG.random_01(); }
00375   void sample_vector(int size, vec &out) {
00376     out.set_size(size, false);
00377     for (int i = 0; i < size; i++) out(i) = sample();
00378   }
00380   void sample_matrix(int rows, int cols, mat &out) {
00381     out.set_size(rows, cols, false);
00382     for (int i = 0; i < rows*cols; i++) out(i) = sample();
00383   }
00384 protected:
00385 private:
00387   double lo_bound, hi_bound;
00389   Random_Generator RNG;
00390 };
00391 
00396 class Exponential_RNG
00397 {
00398 public:
00400   Exponential_RNG(double lambda = 1.0);
00402   void setup(double lambda) { l = lambda; }
00404   double get_setup() const;
00406   double operator()() { return sample(); }
00408   vec operator()(int n);
00410   mat operator()(int h, int w);
00411 protected:
00412 private:
00414   double sample() {  return (-std::log(RNG.random_01()) / l); }
00416   double l;
00418   Random_Generator RNG;
00419 };
00420 
00435 class Normal_RNG
00436 {
00437 public:
00439   Normal_RNG(double meanval, double variance):
00440       mean(meanval), sigma(std::sqrt(variance)) {}
00442   Normal_RNG(): mean(0.0), sigma(1.0) {}
00444   void setup(double meanval, double variance)
00445   { mean = meanval; sigma = std::sqrt(variance); }
00447   void get_setup(double &meanval, double &variance) const;
00449   double operator()() { return (sigma*sample() + mean); }
00451   vec operator()(int n) {
00452     vec temp(n);
00453     sample_vector(n, temp);
00454     temp *= sigma;
00455     temp += mean;
00456     return temp;
00457   }
00459   mat operator()(int h, int w) {
00460     mat temp(h, w);
00461     sample_matrix(h, w, temp);
00462     temp *= sigma;
00463     temp += mean;
00464     return temp;
00465   }
00467   double sample();
00468 
00470   void sample_vector(int size, vec &out) {
00471     out.set_size(size, false);
00472     for (int i = 0; i < size; i++) out(i) = sample();
00473   }
00474 
00476   void sample_matrix(int rows, int cols, mat &out) {
00477     out.set_size(rows, cols, false);
00478     for (int i = 0; i < rows*cols; i++) out(i) = sample();
00479   }
00480 private:
00481   double mean, sigma;
00482   static const double ytab[128];
00483   static const unsigned int ktab[128];
00484   static const double wtab[128];
00485   static const double PARAM_R;
00486   Random_Generator RNG;
00487 };
00488 
00493 class Laplace_RNG
00494 {
00495 public:
00497   Laplace_RNG(double meanval = 0.0, double variance = 1.0);
00499   void setup(double meanval, double variance);
00501   void get_setup(double &meanval, double &variance) const;
00503   double operator()() { return sample(); }
00505   vec operator()(int n);
00507   mat operator()(int h, int w);
00509   double sample() {
00510     double u = RNG.random_01();
00511     double l = sqrt_12var;
00512     if (u < 0.5)
00513       l *= std::log(2.0 * u);
00514     else
00515       l *= -std::log(2.0 * (1 - u));
00516     return (l + mean);
00517   }
00518 private:
00519   double mean, var, sqrt_12var;
00520   Random_Generator RNG;
00521 };
00522 
00527 class Complex_Normal_RNG
00528 {
00529 public:
00531   Complex_Normal_RNG(std::complex<double> mean, double variance):
00532       norm_factor(1.0 / std::sqrt(2.0)) {
00533     setup(mean, variance);
00534   }
00536   Complex_Normal_RNG(): m(0.0), sigma(1.0), norm_factor(1.0 / std::sqrt(2.0)) {}
00538   void setup(std::complex<double> mean, double variance) {
00539     m = mean;
00540     sigma = std::sqrt(variance);
00541   }
00543   void get_setup(std::complex<double> &mean, double &variance) {
00544     mean = m;
00545     variance = sigma * sigma;
00546   }
00548   std::complex<double> operator()() { return sigma*sample() + m; }
00550   cvec operator()(int n) {
00551     cvec temp(n);
00552     sample_vector(n, temp);
00553     temp *= sigma;
00554     temp += m;
00555     return temp;
00556   }
00558   cmat operator()(int h, int w) {
00559     cmat temp(h, w);
00560     sample_matrix(h, w, temp);
00561     temp *= sigma;
00562     temp += m;
00563     return temp;
00564   }
00566   std::complex<double> sample() {
00567     double a = nRNG.sample() * norm_factor;
00568     double b = nRNG.sample() * norm_factor;
00569     return std::complex<double>(a, b);
00570   }
00571 
00573   void sample_vector(int size, cvec &out) {
00574     out.set_size(size, false);
00575     for (int i = 0; i < size; i++) out(i) = sample();
00576   }
00577 
00579   void sample_matrix(int rows, int cols, cmat &out) {
00580     out.set_size(rows, cols, false);
00581     for (int i = 0; i < rows*cols; i++) out(i) = sample();
00582   }
00583 
00585   Complex_Normal_RNG & operator=(const Complex_Normal_RNG&) { return *this; }
00586 
00587 private:
00588   std::complex<double> m;
00589   double sigma;
00590   const double norm_factor;
00591   Normal_RNG nRNG;
00592 };
00593 
00598 class AR1_Normal_RNG
00599 {
00600 public:
00602   AR1_Normal_RNG(double meanval = 0.0, double variance = 1.0,
00603                  double rho = 0.0);
00605   void setup(double meanval, double variance, double rho);
00607   void get_setup(double &meanval, double &variance, double &rho) const;
00609   void reset();
00611   double operator()() { return sample(); }
00613   vec operator()(int n);
00615   mat operator()(int h, int w);
00616 private:
00617   double sample() {
00618     mem *= r;
00619     if (odd) {
00620       r1 = m_2pi * RNG.random_01();
00621       r2 = std::sqrt(factr * std::log(RNG.random_01()));
00622       mem += r2 * std::cos(r1);
00623     }
00624     else {
00625       mem += r2 * std::sin(r1);
00626     }
00627     odd = !odd;
00628     return (mem + mean);
00629   }
00630   double mem, r, factr, mean, var, r1, r2;
00631   bool odd;
00632   Random_Generator RNG;
00633 };
00634 
00639 typedef Normal_RNG Gauss_RNG;
00640 
00645 typedef AR1_Normal_RNG AR1_Gauss_RNG;
00646 
00651 class Weibull_RNG
00652 {
00653 public:
00655   Weibull_RNG(double lambda = 1.0, double beta = 1.0);
00657   void setup(double lambda, double beta);
00659   void get_setup(double &lambda, double &beta) { lambda = l; beta = b; }
00661   double operator()() { return sample(); }
00663   vec operator()(int n);
00665   mat operator()(int h, int w);
00666 private:
00667   double sample() {
00668     return (std::pow(-std::log(RNG.random_01()), 1.0 / b) / l);
00669   }
00670   double l, b;
00671   double mean, var;
00672   Random_Generator RNG;
00673 };
00674 
00679 class Rayleigh_RNG
00680 {
00681 public:
00683   Rayleigh_RNG(double sigma = 1.0);
00685   void setup(double sigma) { sig = sigma; }
00687   double get_setup() { return sig; }
00689   double operator()() { return sample(); }
00691   vec operator()(int n);
00693   mat operator()(int h, int w);
00694 private:
00695   double sample() {
00696     double s1 = nRNG.sample();
00697     double s2 = nRNG.sample();
00698     // s1 and s2 are N(0,1) and independent
00699     return (sig * std::sqrt(s1*s1 + s2*s2));
00700   }
00701   double sig;
00702   Normal_RNG nRNG;
00703 };
00704 
00709 class Rice_RNG
00710 {
00711 public:
00713   Rice_RNG(double sigma = 1.0, double v = 1.0);
00715   void setup(double sigma, double v) { sig = sigma; s = v; }
00717   void get_setup(double &sigma, double &v) { sigma = sig; v = s; }
00719   double operator()() { return sample(); }
00721   vec operator()(int n);
00723   mat operator()(int h, int w);
00724 private:
00725   double sample() {
00726     double s1 = nRNG.sample() + s;
00727     double s2 = nRNG.sample();
00728     // s1 and s2 are N(0,1) and independent
00729     return (sig * std::sqrt(s1*s1 + s2*s2));
00730   }
00731   double sig, s;
00732   Normal_RNG nRNG;
00733 };
00734 
00737 
00739 inline bin randb(void) { Bernoulli_RNG src; return src.sample(); }
00741 inline void randb(int size, bvec &out) { Bernoulli_RNG src; src.sample_vector(size, out); }
00743 inline bvec randb(int size) { bvec temp; randb(size, temp); return temp; }
00745 inline void randb(int rows, int cols, bmat &out) { Bernoulli_RNG src; src.sample_matrix(rows, cols, out); }
00747 inline bmat randb(int rows, int cols) { bmat temp; randb(rows, cols, temp); return temp; }
00748 
00750 inline double randu(void) { Uniform_RNG src; return src.sample(); }
00752 inline void randu(int size, vec &out) { Uniform_RNG src; src.sample_vector(size, out); }
00754 inline vec randu(int size) { vec temp; randu(size, temp); return temp; }
00756 inline void randu(int rows, int cols, mat &out) { Uniform_RNG src; src.sample_matrix(rows, cols, out); }
00758 inline mat randu(int rows, int cols) { mat temp; randu(rows, cols, temp); return temp; }
00759 
00761 inline int randi(int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(); }
00763 inline ivec randi(int size, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(size); }
00765 inline imat randi(int rows, int cols, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(rows, cols); }
00766 
00768 inline vec randray(int size, double sigma = 1.0) { Rayleigh_RNG src; src.setup(sigma); return src(size); }
00769 
00771 inline vec randrice(int size, double sigma = 1.0, double s = 1.0) { Rice_RNG src; src.setup(sigma, s); return src(size); }
00772 
00774 inline vec randexp(int size, double lambda = 1.0) { Exponential_RNG src; src.setup(lambda); return src(size); }
00775 
00777 inline double randn(void) { Normal_RNG src; return src.sample(); }
00779 inline void randn(int size, vec &out) { Normal_RNG src; src.sample_vector(size, out); }
00781 inline vec randn(int size) { vec temp; randn(size, temp); return temp; }
00783 inline void randn(int rows, int cols, mat &out)  { Normal_RNG src; src.sample_matrix(rows, cols, out); }
00785 inline mat randn(int rows, int cols) { mat temp; randn(rows, cols, temp); return temp; }
00786 
00791 inline std::complex<double> randn_c(void) { Complex_Normal_RNG src; return src.sample(); }
00793 inline void randn_c(int size, cvec &out)  { Complex_Normal_RNG src; src.sample_vector(size, out); }
00795 inline cvec randn_c(int size) { cvec temp; randn_c(size, temp); return temp; }
00797 inline void randn_c(int rows, int cols, cmat &out) { Complex_Normal_RNG src; src.sample_matrix(rows, cols, out); }
00799 inline cmat randn_c(int rows, int cols) { cmat temp; randn_c(rows, cols, temp); return temp; }
00800 
00802 
00803 } // namespace itpp
00804 
00805 #endif // #ifndef RANDOM_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Tue Feb 2 09:33:29 2010 for IT++ by Doxygen 1.6.2