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
Generated on Sun Dec 9 17:26:17 2007 for IT++ by Doxygen 1.5.4