00001 00029 #ifndef RANDOM_H 00030 #define RANDOM_H 00031 00032 #include <itpp/base/random_dsfmt.h> 00033 #include <itpp/base/operators.h> 00034 00035 00036 namespace itpp 00037 { 00038 00040 00050 typedef DSFMT_19937_RNG Random_Generator; 00051 00054 00056 void RNG_reset(unsigned int seed); 00058 void RNG_reset(); 00060 void RNG_randomize(); 00062 void RNG_get_state(ivec &state); 00064 void RNG_set_state(const ivec &state); 00066 00071 class Bernoulli_RNG 00072 { 00073 public: 00075 Bernoulli_RNG(double prob) { setup(prob); } 00077 Bernoulli_RNG() { p = 0.5; } 00079 void setup(double prob) { 00080 it_assert(prob >= 0.0 && prob <= 1.0, "The Bernoulli source probability " 00081 "must be between 0 and 1"); 00082 p = prob; 00083 } 00085 double get_setup() const { return p; } 00087 bin operator()() { return sample(); } 00089 bvec operator()(int n) { bvec temp(n); sample_vector(n, temp); return temp; } 00091 bmat operator()(int h, int w) { bmat temp(h, w); sample_matrix(h, w, temp); return temp; } 00093 bin sample() { return RNG.genrand_close_open() < p ? bin(1) : bin(0); } 00095 void sample_vector(int size, bvec &out) { 00096 out.set_size(size, false); 00097 for (int i = 0; i < size; i++) out(i) = sample(); 00098 } 00100 void sample_matrix(int rows, int cols, bmat &out) { 00101 out.set_size(rows, cols, false); 00102 for (int i = 0; i < rows*cols; i++) out(i) = sample(); 00103 } 00104 protected: 00105 private: 00107 double p; 00109 Random_Generator RNG; 00110 }; 00111 00129 class I_Uniform_RNG 00130 { 00131 public: 00133 I_Uniform_RNG(int min = 0, int max = 1); 00135 void setup(int min, int max); 00137 void get_setup(int &min, int &max) const; 00139 int operator()() { return sample(); } 00141 ivec operator()(int n); 00143 imat operator()(int h, int w); 00145 int sample() { 00146 return floor_i(RNG.genrand_close_open() * (hi - lo + 1)) + lo; 00147 } 00148 private: 00150 int lo; 00152 int hi; 00154 Random_Generator RNG; 00155 }; 00156 00161 class Uniform_RNG 00162 { 00163 public: 00165 Uniform_RNG(double min = 0, double max = 1.0); 00167 void setup(double min, double max); 00169 void get_setup(double &min, double &max) const; 00171 double operator()() { return (sample() * (hi_bound - lo_bound) + lo_bound); } 00173 vec operator()(int n) { 00174 vec temp(n); 00175 sample_vector(n, temp); 00176 temp *= hi_bound - lo_bound; 00177 temp += lo_bound; 00178 return temp; 00179 } 00181 mat operator()(int h, int w) { 00182 mat temp(h, w); 00183 sample_matrix(h, w, temp); 00184 temp *= hi_bound - lo_bound; 00185 temp += lo_bound; 00186 return temp; 00187 } 00189 double sample() { return RNG.genrand_close_open(); } 00191 void sample_vector(int size, vec &out) { 00192 out.set_size(size, false); 00193 for (int i = 0; i < size; i++) out(i) = sample(); 00194 } 00196 void sample_matrix(int rows, int cols, mat &out) { 00197 out.set_size(rows, cols, false); 00198 for (int i = 0; i < rows*cols; i++) out(i) = sample(); 00199 } 00200 protected: 00201 private: 00203 double lo_bound, hi_bound; 00205 Random_Generator RNG; 00206 }; 00207 00212 class Exponential_RNG 00213 { 00214 public: 00216 Exponential_RNG(double lambda = 1.0); 00218 void setup(double lambda) { l = lambda; } 00220 double get_setup() const; 00222 double operator()() { return sample(); } 00224 vec operator()(int n); 00226 mat operator()(int h, int w); 00227 private: 00228 double sample() { return (-std::log(RNG.genrand_open_close()) / l); } 00229 double l; 00230 Random_Generator RNG; 00231 }; 00232 00247 class Normal_RNG 00248 { 00249 public: 00251 Normal_RNG(double meanval, double variance): 00252 mean(meanval), sigma(std::sqrt(variance)) {} 00254 Normal_RNG(): mean(0.0), sigma(1.0) {} 00256 void setup(double meanval, double variance) 00257 { mean = meanval; sigma = std::sqrt(variance); } 00259 void get_setup(double &meanval, double &variance) const; 00261 double operator()() { return (sigma*sample() + mean); } 00263 vec operator()(int n) { 00264 vec temp(n); 00265 sample_vector(n, temp); 00266 temp *= sigma; 00267 temp += mean; 00268 return temp; 00269 } 00271 mat operator()(int h, int w) { 00272 mat temp(h, w); 00273 sample_matrix(h, w, temp); 00274 temp *= sigma; 00275 temp += mean; 00276 return temp; 00277 } 00279 double sample(); 00280 00282 void sample_vector(int size, vec &out) { 00283 out.set_size(size, false); 00284 for (int i = 0; i < size; i++) out(i) = sample(); 00285 } 00286 00288 void sample_matrix(int rows, int cols, mat &out) { 00289 out.set_size(rows, cols, false); 00290 for (int i = 0; i < rows*cols; i++) out(i) = sample(); 00291 } 00292 private: 00293 double mean, sigma; 00294 static const double ytab[128]; 00295 static const unsigned int ktab[128]; 00296 static const double wtab[128]; 00297 static const double PARAM_R; 00298 Random_Generator RNG; 00299 }; 00300 00317 class Gamma_RNG 00318 { 00319 public: 00321 Gamma_RNG(double a = 1.0, double b = 1.0): alpha(a), beta(b) {} 00323 void setup(double a, double b) { alpha = a; beta = b; } 00325 double operator()() { return sample(); } 00327 vec operator()(int n); 00329 mat operator()(int r, int c); 00331 double sample(); 00332 private: 00333 double alpha; 00334 double beta; 00335 Random_Generator RNG; 00336 Normal_RNG NRNG; 00337 }; 00338 00343 class Laplace_RNG 00344 { 00345 public: 00347 Laplace_RNG(double meanval = 0.0, double variance = 1.0); 00349 void setup(double meanval, double variance); 00351 void get_setup(double &meanval, double &variance) const; 00353 double operator()() { return sample(); } 00355 vec operator()(int n); 00357 mat operator()(int h, int w); 00359 double sample() { 00360 double u = RNG.genrand_open_open(); 00361 double l = sqrt_12var; 00362 if (u < 0.5) 00363 l *= std::log(2.0 * u); 00364 else 00365 l *= -std::log(2.0 * (1 - u)); 00366 return (l + mean); 00367 } 00368 private: 00369 double mean, var, sqrt_12var; 00370 Random_Generator RNG; 00371 }; 00372 00377 class Complex_Normal_RNG 00378 { 00379 public: 00381 Complex_Normal_RNG(std::complex<double> mean, double variance): 00382 norm_factor(1.0 / std::sqrt(2.0)) { 00383 setup(mean, variance); 00384 } 00386 Complex_Normal_RNG(): m(0.0), sigma(1.0), norm_factor(1.0 / std::sqrt(2.0)) {} 00388 void setup(std::complex<double> mean, double variance) { 00389 m = mean; 00390 sigma = std::sqrt(variance); 00391 } 00393 void get_setup(std::complex<double> &mean, double &variance) { 00394 mean = m; 00395 variance = sigma * sigma; 00396 } 00398 std::complex<double> operator()() { return sigma*sample() + m; } 00400 cvec operator()(int n) { 00401 cvec temp(n); 00402 sample_vector(n, temp); 00403 temp *= sigma; 00404 temp += m; 00405 return temp; 00406 } 00408 cmat operator()(int h, int w) { 00409 cmat temp(h, w); 00410 sample_matrix(h, w, temp); 00411 temp *= sigma; 00412 temp += m; 00413 return temp; 00414 } 00416 std::complex<double> sample() { 00417 double a = nRNG.sample() * norm_factor; 00418 double b = nRNG.sample() * norm_factor; 00419 return std::complex<double>(a, b); 00420 } 00421 00423 void sample_vector(int size, cvec &out) { 00424 out.set_size(size, false); 00425 for (int i = 0; i < size; i++) out(i) = sample(); 00426 } 00427 00429 void sample_matrix(int rows, int cols, cmat &out) { 00430 out.set_size(rows, cols, false); 00431 for (int i = 0; i < rows*cols; i++) out(i) = sample(); 00432 } 00433 00435 Complex_Normal_RNG & operator=(const Complex_Normal_RNG&) { return *this; } 00436 00437 private: 00438 std::complex<double> m; 00439 double sigma; 00440 const double norm_factor; 00441 Normal_RNG nRNG; 00442 }; 00443 00448 class AR1_Normal_RNG 00449 { 00450 public: 00452 AR1_Normal_RNG(double meanval = 0.0, double variance = 1.0, 00453 double rho = 0.0); 00455 void setup(double meanval, double variance, double rho); 00457 void get_setup(double &meanval, double &variance, double &rho) const; 00459 void reset(); 00461 double operator()() { return sample(); } 00463 vec operator()(int n); 00465 mat operator()(int h, int w); 00466 private: 00467 double sample() { 00468 mem *= r; 00469 if (odd) { 00470 r1 = m_2pi * RNG.genrand_open_close(); 00471 r2 = std::sqrt(factr * std::log(RNG.genrand_open_close())); 00472 mem += r2 * std::cos(r1); 00473 } 00474 else { 00475 mem += r2 * std::sin(r1); 00476 } 00477 odd = !odd; 00478 return (mem + mean); 00479 } 00480 double mem, r, factr, mean, var, r1, r2; 00481 bool odd; 00482 Random_Generator RNG; 00483 }; 00484 00489 typedef Normal_RNG Gauss_RNG; 00490 00495 typedef AR1_Normal_RNG AR1_Gauss_RNG; 00496 00501 class Weibull_RNG 00502 { 00503 public: 00505 Weibull_RNG(double lambda = 1.0, double beta = 1.0); 00507 void setup(double lambda, double beta); 00509 void get_setup(double &lambda, double &beta) { lambda = l; beta = b; } 00511 double operator()() { return sample(); } 00513 vec operator()(int n); 00515 mat operator()(int h, int w); 00516 private: 00517 double sample() { 00518 return (std::pow(-std::log(RNG.genrand_open_close()), 1.0 / b) / l); 00519 } 00520 double l, b; 00521 double mean, var; 00522 Random_Generator RNG; 00523 }; 00524 00529 class Rayleigh_RNG 00530 { 00531 public: 00533 Rayleigh_RNG(double sigma = 1.0); 00535 void setup(double sigma) { sig = sigma; } 00537 double get_setup() { return sig; } 00539 double operator()() { return sample(); } 00541 vec operator()(int n); 00543 mat operator()(int h, int w); 00544 private: 00545 double sample() { 00546 double s1 = nRNG.sample(); 00547 double s2 = nRNG.sample(); 00548 // s1 and s2 are N(0,1) and independent 00549 return (sig * std::sqrt(s1*s1 + s2*s2)); 00550 } 00551 double sig; 00552 Normal_RNG nRNG; 00553 }; 00554 00559 class Rice_RNG 00560 { 00561 public: 00563 Rice_RNG(double sigma = 1.0, double v = 1.0); 00565 void setup(double sigma, double v) { sig = sigma; s = v; } 00567 void get_setup(double &sigma, double &v) { sigma = sig; v = s; } 00569 double operator()() { return sample(); } 00571 vec operator()(int n); 00573 mat operator()(int h, int w); 00574 private: 00575 double sample() { 00576 double s1 = nRNG.sample() + s; 00577 double s2 = nRNG.sample(); 00578 // s1 and s2 are N(0,1) and independent 00579 return (sig * std::sqrt(s1*s1 + s2*s2)); 00580 } 00581 double sig, s; 00582 Normal_RNG nRNG; 00583 }; 00584 00587 00589 inline bin randb(void) { Bernoulli_RNG src; return src.sample(); } 00591 inline void randb(int size, bvec &out) { Bernoulli_RNG src; src.sample_vector(size, out); } 00593 inline bvec randb(int size) { bvec temp; randb(size, temp); return temp; } 00595 inline void randb(int rows, int cols, bmat &out) { Bernoulli_RNG src; src.sample_matrix(rows, cols, out); } 00597 inline bmat randb(int rows, int cols) { bmat temp; randb(rows, cols, temp); return temp; } 00598 00600 inline double randu(void) { Uniform_RNG src; return src.sample(); } 00602 inline void randu(int size, vec &out) { Uniform_RNG src; src.sample_vector(size, out); } 00604 inline vec randu(int size) { vec temp; randu(size, temp); return temp; } 00606 inline void randu(int rows, int cols, mat &out) { Uniform_RNG src; src.sample_matrix(rows, cols, out); } 00608 inline mat randu(int rows, int cols) { mat temp; randu(rows, cols, temp); return temp; } 00609 00611 inline int randi(int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(); } 00613 inline ivec randi(int size, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(size); } 00615 inline imat randi(int rows, int cols, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(rows, cols); } 00616 00618 inline vec randray(int size, double sigma = 1.0) { Rayleigh_RNG src; src.setup(sigma); return src(size); } 00619 00621 inline vec randrice(int size, double sigma = 1.0, double s = 1.0) { Rice_RNG src; src.setup(sigma, s); return src(size); } 00622 00624 inline vec randexp(int size, double lambda = 1.0) { Exponential_RNG src; src.setup(lambda); return src(size); } 00625 00627 inline double randn(void) { Normal_RNG src; return src.sample(); } 00629 inline void randn(int size, vec &out) { Normal_RNG src; src.sample_vector(size, out); } 00631 inline vec randn(int size) { vec temp; randn(size, temp); return temp; } 00633 inline void randn(int rows, int cols, mat &out) { Normal_RNG src; src.sample_matrix(rows, cols, out); } 00635 inline mat randn(int rows, int cols) { mat temp; randn(rows, cols, temp); return temp; } 00636 00641 inline std::complex<double> randn_c(void) { Complex_Normal_RNG src; return src.sample(); } 00643 inline void randn_c(int size, cvec &out) { Complex_Normal_RNG src; src.sample_vector(size, out); } 00645 inline cvec randn_c(int size) { cvec temp; randn_c(size, temp); return temp; } 00647 inline void randn_c(int rows, int cols, cmat &out) { Complex_Normal_RNG src; src.sample_matrix(rows, cols, out); } 00649 inline cmat randn_c(int rows, int cols) { cmat temp; randn_c(rows, cols, temp); return temp; } 00650 00652 00653 } // namespace itpp 00654 00655 #endif // #ifndef RANDOM_H
Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4