00001 00033 #ifndef RANDOM_H 00034 #define RANDOM_H 00035 00036 #include <itpp/base/operators.h> 00037 #include <ctime> 00038 00039 00040 namespace itpp { 00041 00043 00115 class Random_Generator { 00116 public: 00118 Random_Generator(); 00120 void randomize(); 00122 void reset() { initialize(lastSeed); reload(); } 00124 void reset(unsigned int seed) { lastSeed = seed; initialize(lastSeed); reload(); } 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 (double(random_int()) + 0.5) * (1.0/4294967296.0); } 00142 double random_01_closed() { return double(random_int()) * (1.0/4294967295.0); } 00144 void get_state(ivec &out_state); 00146 void set_state(ivec &new_state); 00147 protected: 00148 private: 00150 unsigned int lastSeed; 00152 static unsigned int state[624]; 00154 static unsigned int *pNext; 00156 static int left; 00157 00159 void initialize( unsigned int seed ) 00160 { 00161 // Initialize generator state with seed 00162 // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier. 00163 // In previous versions, most significant bits (MSBs) of the seed affect 00164 // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto. 00165 register unsigned int *s = state; 00166 register unsigned int *r = state; 00167 register int i = 1; 00168 *s++ = seed & 0xffffffffU; 00169 for( ; i < 624; ++i ) 00170 { 00171 *s++ = ( 1812433253U * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffU; 00172 r++; 00173 } 00174 } 00176 void reload() 00177 { 00178 // Generate N new values in state 00179 // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com) 00180 register unsigned int *p = state; 00181 register int i; 00182 for( i = 624 - 397; i--; ++p ) 00183 *p = twist( p[397], p[0], p[1] ); 00184 for( i = 397; --i; ++p ) 00185 *p = twist( p[397-624], p[0], p[1] ); 00186 *p = twist( p[397-624], p[0], state[0] ); 00187 00188 left = 624, pNext = state; 00189 } 00191 unsigned int hiBit( const unsigned int& u ) const { return u & 0x80000000U; } 00193 unsigned int loBit( const unsigned int& u ) const { return u & 0x00000001U; } 00195 unsigned int loBits( const unsigned int& u ) const { return u & 0x7fffffffU; } 00197 unsigned int mixBits( const unsigned int& u, const unsigned int& v ) const 00198 { return hiBit(u) | loBits(v); } 00199 00200 /* 00201 * ---------------------------------------------------------------------- 00202 * --- ediap - 2007/01/17 --- 00203 * ---------------------------------------------------------------------- 00204 * Richard's implementation of the twist() function is was follows: 00205 * { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfU); } 00206 * However, this code caused a warning/error under MSVC++, because 00207 * unsigned value loBit(s1) is being negated with `-' (c.f. bug report 00208 * [1635988]). I changed this to the same implementation as is used in 00209 * original C sources of Mersenne Twister RNG: 00210 * #define MATRIX_A 0x9908b0dfUL 00211 * #define UMASK 0x80000000UL 00212 * #define LMASK 0x7fffffffUL 00213 * #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) 00214 * #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) 00215 * ---------------------------------------------------------------------- 00216 */ 00218 unsigned int twist( const unsigned int& m, const unsigned int& s0, 00219 const unsigned int& s1 ) const 00220 { return m ^ (mixBits(s0,s1)>>1) ^ (loBit(s1) ? 0x9908b0dfU : 0U); } 00222 unsigned int hash( time_t t, clock_t c ); 00223 }; 00224 00225 00228 00230 void RNG_reset(unsigned long seed); 00232 void RNG_reset(); 00234 void RNG_randomize(); 00236 void RNG_get_state(ivec &state); 00238 void RNG_set_state(ivec &state); 00240 00245 class Bernoulli_RNG { 00246 public: 00248 Bernoulli_RNG(double prob) { setup(prob); } 00250 Bernoulli_RNG() { p=0.5; } 00252 void setup(double prob) 00253 { 00254 it_assert(prob>=0.0 && prob<=1.0, "The bernoulli source probability must be between 0 and 1"); 00255 p = prob; 00256 } 00258 double get_setup() const { return p; } 00260 bin operator()() { return sample(); } 00262 bvec operator()(int n) { bvec temp(n); sample_vector(n, temp); return temp; } 00264 bmat operator()(int h, int w) { bmat temp(h, w); sample_matrix(h, w, temp); return temp; } 00266 bin sample() { return bin( RNG.random_01() < p ? 1 : 0 ); } 00268 void sample_vector(int size, bvec &out) 00269 { 00270 out.set_size(size, false); 00271 for (int i=0; i<size; i++) out(i) = sample(); 00272 } 00274 void sample_matrix(int rows, int cols, bmat &out) 00275 { 00276 out.set_size(rows, cols, false); 00277 for (int i=0; i<rows*cols; i++) out(i) = sample(); 00278 } 00279 protected: 00280 private: 00282 double p; 00284 Random_Generator RNG; 00285 }; 00286 00304 class I_Uniform_RNG { 00305 public: 00307 I_Uniform_RNG(int min=0, int max=1); 00309 void setup(int min, int max); 00311 void get_setup(int &min, int &max) const; 00313 int operator()() { return sample(); } 00315 ivec operator()(int n); 00317 imat operator()(int h, int w); 00319 int sample() { return ( floor_i(RNG.random_01() * (hi - lo + 1)) + lo ); } 00320 protected: 00321 private: 00323 int lo; 00325 int hi; 00327 Random_Generator RNG; 00328 }; 00329 00334 class Uniform_RNG { 00335 public: 00337 Uniform_RNG(double min=0, double max=1.0); 00339 void setup(double min, double max); 00341 void get_setup(double &min, double &max) const; 00343 double operator()() { return ( sample()* (hi_bound - lo_bound) + lo_bound ); } 00345 vec operator()(int n) 00346 { vec temp(n); sample_vector(n, temp); return (temp *(hi_bound - lo_bound) + lo_bound); } 00348 mat operator()(int h, int w) 00349 { mat temp(h,w); sample_matrix(h,w, temp); return (temp *(hi_bound - lo_bound) + lo_bound); } 00351 double sample() { return RNG.random_01(); } 00353 void sample_vector(int size, vec &out) 00354 { 00355 out.set_size(size, false); 00356 for (int i=0; i<size; i++) out(i) = sample(); 00357 } 00359 void sample_matrix(int rows, int cols, mat &out) 00360 { 00361 out.set_size(rows, cols, false); 00362 for (int i=0; i<rows*cols; i++) out(i) = sample(); 00363 } 00364 protected: 00365 private: 00367 double lo_bound, hi_bound; 00369 Random_Generator RNG; 00370 }; 00371 00376 class Exponential_RNG { 00377 public: 00379 Exponential_RNG(double lambda=1.0); 00381 void setup(double lambda) { l=lambda; } 00383 double get_setup() const; 00385 double operator()() { return sample(); } 00387 vec operator()(int n); 00389 mat operator()(int h, int w); 00390 protected: 00391 private: 00393 double sample() { return ( -std::log(RNG.random_01()) / l ); } 00395 double l; 00397 Random_Generator RNG; 00398 }; 00399 00404 class Normal_RNG { 00405 public: 00407 Normal_RNG(double meanval, double variance) { setup(meanval, variance); }; 00409 Normal_RNG() { mean = 0.0; sigma = 1.0; odd = false; } 00411 void setup(double meanval, double variance) 00412 { mean = meanval; sigma = std::sqrt(variance); odd = false; } 00414 void get_setup(double &meanval, double &variance) const; 00416 double operator()() { return (sigma*sample()+mean); } 00418 vec operator()(int n) { vec temp(n); sample_vector(n, temp); return (sigma*temp+mean); } 00420 mat operator()(int h, int w) { mat temp(h,w); sample_matrix(h,w, temp); return (sigma*temp+mean); } 00422 double sample() // Box-Mueller, Polar version, See Knuth v2, 3rd ed, p122 00423 { 00424 double x, r2; 00425 if (odd) { 00426 odd = 0; 00427 return y; 00428 } 00429 do { x = 2. * RNG.random_01() - 1.; y = 2. * RNG.random_01() - 1.; r2 = x*x + y*y; } while (r2 >= 1. || r2 < 1E-30); 00430 r2 = std::sqrt(std::log(r2)*(-2./r2)); x *= r2; y *= r2; odd = 1; 00431 return x; 00432 } 00433 00435 void sample_vector(int size, vec &out) 00436 { 00437 out.set_size(size, false); 00438 for (int i=0; i<size; i++) out(i) = sample(); 00439 } 00440 00442 void sample_matrix(int rows, int cols, mat &out) 00443 { 00444 out.set_size(rows, cols, false); 00445 for (int i=0; i<rows*cols; i++) out(i) = sample(); 00446 } 00447 protected: 00448 private: 00450 double mean, sigma, y; 00452 int odd; 00454 Random_Generator RNG; 00455 }; 00456 00461 class Laplace_RNG { 00462 public: 00464 Laplace_RNG(double meanval=0.0, double variance=1.0); 00466 void setup(double meanval, double variance); 00468 void get_setup(double &meanval, double &variance) const; 00470 double operator()() { return sample(); } 00472 vec operator()(int n); 00474 mat operator()(int h, int w); 00476 double sample() 00477 { 00478 u=RNG.random_01(); 00479 if(u<0.5) 00480 l=std::log(2.0*u); 00481 else 00482 l=-std::log(2.0*(1-u)); 00483 00484 l *= std::sqrt(var/2.0); 00485 l += mean; 00486 00487 return l; 00488 } 00489 protected: 00490 private: 00492 double mean, var, u, l; 00494 Random_Generator RNG; 00495 }; 00496 00501 class Complex_Normal_RNG { 00502 public: 00504 Complex_Normal_RNG(std::complex<double> mean, double variance) { setup(mean, variance); } 00506 Complex_Normal_RNG() { m = 0.0; sigma=1.0; } 00508 void setup(std::complex<double> mean, double variance) { m = mean; sigma = std::sqrt(variance); } 00510 void get_setup(std::complex<double> &mean, double &variance) { mean = m; variance = sigma*sigma; } 00512 std::complex<double> operator()() { return sigma*sample()+m; } 00514 cvec operator()(int n) { cvec temp(n); sample_vector(n, temp); return (sigma*temp+m); } 00516 cmat operator()(int h, int w) { cmat temp(h,w); sample_matrix(h,w, temp); return (sigma*temp+m); } 00518 std::complex<double> sample() 00519 { 00520 double a = m_2pi * RNG.random_01(); 00521 double b = std::sqrt(-std::log(RNG.random_01())); 00522 return std::complex<double>(b * std::cos(a), b * std::sin(a)); 00523 } 00524 00526 void sample_vector(int size, cvec &out) 00527 { 00528 out.set_size(size, false); 00529 for (int i=0; i<size; i++) out(i) = sample(); 00530 } 00531 00533 void sample_matrix(int rows, int cols, cmat &out) 00534 { 00535 out.set_size(rows, cols, false); 00536 for (int i=0; i<rows*cols; i++) out(i) = sample(); 00537 } 00538 protected: 00539 private: 00541 std::complex<double> m; 00543 double sigma; 00545 Random_Generator RNG; 00546 }; 00547 00552 class AR1_Normal_RNG { 00553 public: 00555 AR1_Normal_RNG(double meanval=0.0, double variance=1.0, double rho=0.0); 00557 void setup(double meanval, double variance, double rho); 00559 void get_setup(double &meanval, double &variance, double &rho) const; 00561 void reset(); 00563 double operator()() { return sample(); } 00565 vec operator()(int n); 00567 mat operator()(int h, int w); 00568 00569 protected: 00570 private: 00572 double sample(); 00574 double my_mean, mem, r, factr, mean, var, r1, r2; 00576 bool odd; 00578 Random_Generator RNG; 00579 }; 00580 00585 typedef Normal_RNG Gauss_RNG; 00586 00591 typedef AR1_Normal_RNG AR1_Gauss_RNG; 00592 00597 class Weibull_RNG { 00598 public: 00600 Weibull_RNG(double lambda=1.0, double beta=1.0); 00601 00603 void setup(double lambda, double beta); 00605 void get_setup(double &lambda, double &beta) { lambda=l; beta=b; } 00607 double operator()() { return sample(); } 00609 vec operator()(int n); 00611 mat operator()(int h, int w); 00612 protected: 00613 private: 00615 double sample(); 00617 double l, b; 00619 double mean, var; 00621 Random_Generator RNG; 00622 }; 00623 00628 class Rayleigh_RNG { 00629 public: 00631 Rayleigh_RNG(double sigma=1.0); 00632 00634 void setup(double sigma) { sig=sigma; } 00636 double get_setup() { return sig; } 00638 double operator()() { return sample(); } 00640 vec operator()(int n); 00642 mat operator()(int h, int w); 00643 protected: 00644 private: 00646 double sample(); 00648 double sig; 00650 Random_Generator RNG; 00651 }; 00652 00657 class Rice_RNG { 00658 public: 00660 Rice_RNG(double sigma=1.0, double _s=1.0); 00662 void setup(double sigma, double _s) { sig=sigma; s=_s; } 00664 void get_setup(double &sigma, double &_s) { sigma=sig; _s=s; } 00666 double operator()() { return sample(); } 00668 vec operator()(int n); 00670 mat operator()(int h, int w); 00671 protected: 00672 private: 00674 double sample(); 00676 double sig, s; 00678 Random_Generator RNG; 00679 }; 00680 00683 00685 inline bin randb(void) { Bernoulli_RNG src; return src.sample(); } 00687 inline void randb(int size, bvec &out) { Bernoulli_RNG src; src.sample_vector(size, out); } 00689 inline bvec randb(int size) { bvec temp; randb(size, temp); return temp; } 00691 inline void randb(int rows, int cols, bmat &out) { Bernoulli_RNG src; src.sample_matrix(rows, cols, out); } 00693 inline bmat randb(int rows, int cols){ bmat temp; randb(rows, cols, temp); return temp; } 00694 00696 inline double randu(void) { Uniform_RNG src; return src.sample(); } 00698 inline void randu(int size, vec &out) { Uniform_RNG src; src.sample_vector(size, out); } 00700 inline vec randu(int size){ vec temp; randu(size, temp); return temp; } 00702 inline void randu(int rows, int cols, mat &out) { Uniform_RNG src; src.sample_matrix(rows, cols, out); } 00704 inline mat randu(int rows, int cols) { mat temp; randu(rows, cols, temp); return temp; } 00705 00707 inline int randi(int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(); } 00709 inline ivec randi(int size, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(size); } 00711 inline imat randi(int rows, int cols, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(rows, cols); } 00712 00714 inline vec randray(int size, double sigma = 1.0) { Rayleigh_RNG src; src.setup(sigma); return src(size); } 00715 00717 inline vec randrice(int size, double sigma = 1.0, double s = 1.0) { Rice_RNG src; src.setup(sigma, s); return src(size); } 00718 00720 inline vec randexp(int size, double lambda = 1.0) { Exponential_RNG src; src.setup(lambda); return src(size); } 00721 00723 inline double randn(void) { Normal_RNG src; return src.sample(); } 00725 inline void randn(int size, vec &out) { Normal_RNG src; src.sample_vector(size, out); } 00727 inline vec randn(int size) { vec temp; randn(size, temp); return temp; } 00729 inline void randn(int rows, int cols, mat &out) { Normal_RNG src; src.sample_matrix(rows, cols, out); } 00731 inline mat randn(int rows, int cols){ mat temp; randn(rows, cols, temp); return temp; } 00732 00737 inline std::complex<double> randn_c(void) { Complex_Normal_RNG src; return src.sample(); } 00739 inline void randn_c(int size, cvec &out) { Complex_Normal_RNG src; src.sample_vector(size, out); } 00741 inline cvec randn_c(int size){ cvec temp; randn_c(size, temp); return temp; } 00743 inline void randn_c(int rows, int cols, cmat &out) { Complex_Normal_RNG src; src.sample_matrix(rows, cols, out); } 00745 inline cmat randn_c(int rows, int cols) { cmat temp; randn_c(rows, cols, temp); return temp; } 00746 00748 00749 // -------------------- INLINES ---------------------------------------------- 00750 00751 inline double AR1_Normal_RNG::sample() 00752 { 00753 double s; 00754 00755 if (odd) { 00756 r1 = RNG.random_01(); 00757 r2 = RNG.random_01(); 00758 s = std::sqrt(factr * std::log(r2)) * std::cos(m_2pi * r1); 00759 } else 00760 s = std::sqrt(factr * std::log(r2)) * std::sin(m_2pi * r1); 00761 00762 odd = !odd; 00763 00764 mem = s + r * mem; 00765 s = mem + mean; 00766 00767 return s; 00768 } 00769 00770 inline double Weibull_RNG::sample() 00771 { 00772 return ( std::pow(-std::log(RNG.random_01()), 1.0/b) / l ); 00773 } 00774 00775 inline double Rayleigh_RNG::sample() 00776 { 00777 double r1, r2; 00778 double s1, s2, samp; 00779 00780 r1 = RNG.random_01(); 00781 r2 = RNG.random_01(); 00782 s1 = std::sqrt(-2.0 * std::log(r2)) * std::cos(m_2pi * r1); 00783 s2 = std::sqrt(-2.0 * std::log(r2)) * std::sin(m_2pi * r1); 00784 // s1 and s2 are N(0,1) and independent 00785 samp = sig * std::sqrt(s1*s1 + s2*s2); 00786 00787 return samp; 00788 } 00789 00790 inline double Rice_RNG::sample() 00791 { 00792 double r1, r2; 00793 double s1, s2, samp; 00794 double m1 = 0.0; 00795 double m2 = std::sqrt(s*s - m1*m1); 00796 00797 r1 = RNG.random_01(); 00798 r2 = RNG.random_01(); 00799 s1 = std::sqrt(-2.0 * std::log(r2)) * std::cos(m_2pi * r1); 00800 s2 = std::sqrt(-2.0 * std::log(r2)) * std::sin(m_2pi * r1); 00801 // s1 and s2 are N(0,1) and independent 00802 samp = std::sqrt((sig*s1+m1) * (sig*s1+m1) + (sig*s2+m2) * (sig*s2+m2)); 00803 00804 return samp; 00805 } 00806 00807 } // namespace itpp 00808 00809 #endif // #ifndef RANDOM_H
Generated on Wed Apr 18 11:45:34 2007 for IT++ by Doxygen 1.5.2