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
Generated on Sat Feb 26 2011 16:06:31 for IT++ by Doxygen 1.7.3