00001 00030 #include <itpp/base/random.h> 00031 #include <itpp/base/math/elem_math.h> 00032 #include <limits> 00033 00034 00035 namespace itpp 00036 { 00037 00039 // Random_Generator 00041 00042 bool Random_Generator::initialized = false; 00043 int Random_Generator::left = 0; 00044 unsigned int Random_Generator::state[624]; 00045 unsigned int *Random_Generator::pNext; 00046 00047 unsigned int Random_Generator::hash(time_t t, clock_t c) 00048 { 00049 // Get a unsigned int from t and c 00050 // Better than uint(x) in case x is floating point in [0,1] 00051 // Based on code by Lawrence Kirby (fred@genesis.demon.co.uk) 00052 static unsigned int differ = 0; // guarantee time-based seeds will change 00053 00054 unsigned int h1 = 0; 00055 unsigned char *p = (unsigned char *) & t; 00056 for (size_t i = 0; i < sizeof(t); ++i) { 00057 h1 *= std::numeric_limits<unsigned char>::max() + 2U; 00058 h1 += p[i]; 00059 } 00060 unsigned int h2 = 0; 00061 p = (unsigned char *) & c; 00062 for (size_t j = 0; j < sizeof(c); ++j) { 00063 h2 *= std::numeric_limits<unsigned char>::max() + 2U; 00064 h2 += p[j]; 00065 } 00066 return (h1 + differ++) ^ h2; 00067 } 00068 00069 void Random_Generator::get_state(ivec &out_state) 00070 { 00071 out_state.set_size(625, false); 00072 for (int i = 0; i < 624; i++) 00073 out_state(i) = state[i]; 00074 00075 out_state(624) = left; // the number of elements left in state before reload 00076 } 00077 00078 void Random_Generator::set_state(ivec &new_state) 00079 { 00080 it_assert(new_state.size() == 625, "Random_Generator::set_state(): Not a valid state vector"); 00081 00082 for (int i = 0; i < 624; i++) 00083 state[i] = new_state(i); 00084 00085 left = new_state(624); 00086 pNext = &state[624-left]; 00087 } 00088 00089 00090 // Set the seed of the Global Random Number Generator 00091 void RNG_reset(unsigned int seed) 00092 { 00093 Random_Generator RNG; 00094 RNG.reset(seed); 00095 } 00096 00097 // Set the seed of the Global Random Number Generator to the same as last time 00098 void RNG_reset() 00099 { 00100 Random_Generator RNG; 00101 RNG.reset(); 00102 } 00103 00104 // Set a random seed for the Global Random Number Generator 00105 void RNG_randomize() 00106 { 00107 Random_Generator RNG; 00108 RNG.randomize(); 00109 } 00110 00111 // Save current full state of generator in memory 00112 void RNG_get_state(ivec &state) 00113 { 00114 Random_Generator RNG; 00115 RNG.get_state(state); 00116 } 00117 00118 // Resume the state saved in memory 00119 void RNG_set_state(ivec &state) 00120 { 00121 Random_Generator RNG; 00122 RNG.set_state(state); 00123 } 00124 00126 // I_Uniform_RNG 00128 00129 I_Uniform_RNG::I_Uniform_RNG(int min, int max) 00130 { 00131 setup(min, max); 00132 } 00133 00134 void I_Uniform_RNG::setup(int min, int max) 00135 { 00136 if (min <= max) { 00137 lo = min; 00138 hi = max; 00139 } 00140 else { 00141 lo = max; 00142 hi = min; 00143 } 00144 } 00145 00146 void I_Uniform_RNG::get_setup(int &min, int &max) const 00147 { 00148 min = lo; 00149 max = hi; 00150 } 00151 00152 ivec I_Uniform_RNG::operator()(int n) 00153 { 00154 ivec vv(n); 00155 00156 for (int i=0; i < n; i++) 00157 vv(i) = sample(); 00158 00159 return vv; 00160 } 00161 00162 imat I_Uniform_RNG::operator()(int h, int w) 00163 { 00164 imat mm(h, w); 00165 int i, j; 00166 00167 for (i = 0; i < h; i++) 00168 for (j = 0; j < w; j++) 00169 mm(i, j) = sample(); 00170 00171 return mm; 00172 } 00173 00175 // Uniform_RNG 00177 00178 Uniform_RNG::Uniform_RNG(double min, double max) 00179 { 00180 setup(min, max); 00181 } 00182 00183 void Uniform_RNG::setup(double min, double max) 00184 { 00185 if (min <= max) { 00186 lo_bound = min; 00187 hi_bound = max; 00188 } 00189 else { 00190 lo_bound = max; 00191 hi_bound = min; 00192 } 00193 } 00194 00195 void Uniform_RNG::get_setup(double &min, double &max) const 00196 { 00197 min = lo_bound; 00198 max = hi_bound; 00199 } 00200 00202 // Exp_RNG 00204 00205 Exponential_RNG::Exponential_RNG(double lambda) 00206 { 00207 setup(lambda); 00208 } 00209 00210 vec Exponential_RNG::operator()(int n) 00211 { 00212 vec vv(n); 00213 00214 for (int i=0; i < n; i++) 00215 vv(i) = sample(); 00216 00217 return vv; 00218 } 00219 00220 mat Exponential_RNG::operator()(int h, int w) 00221 { 00222 mat mm(h, w); 00223 int i, j; 00224 00225 for (i = 0; i < h; i++) 00226 for (j = 0; j < w; j++) 00227 mm(i, j) = sample(); 00228 00229 return mm; 00230 } 00231 00233 // Normal_RNG 00235 00236 void Normal_RNG::get_setup(double &meanval, double &variance) const 00237 { 00238 meanval = mean; 00239 variance = sigma * sigma; 00240 } 00241 00242 // (Ziggurat) tabulated values for the heigt of the Ziggurat levels 00243 const double Normal_RNG::ytab[128] = { 00244 1, 0.963598623011, 0.936280813353, 0.913041104253, 00245 0.892278506696, 0.873239356919, 0.855496407634, 0.838778928349, 00246 0.822902083699, 0.807732738234, 0.793171045519, 0.779139726505, 00247 0.765577436082, 0.752434456248, 0.739669787677, 0.727249120285, 00248 0.715143377413, 0.703327646455, 0.691780377035, 0.68048276891, 00249 0.669418297233, 0.65857233912, 0.647931876189, 0.637485254896, 00250 0.62722199145, 0.617132611532, 0.607208517467, 0.597441877296, 00251 0.587825531465, 0.578352913803, 0.569017984198, 0.559815170911, 00252 0.550739320877, 0.541785656682, 0.532949739145, 0.524227434628, 00253 0.515614886373, 0.507108489253, 0.498704867478, 0.490400854812, 00254 0.482193476986, 0.47407993601, 0.466057596125, 0.458123971214, 00255 0.450276713467, 0.442513603171, 0.434832539473, 0.427231532022, 00256 0.419708693379, 0.41226223212, 0.404890446548, 0.397591718955, 00257 0.390364510382, 0.383207355816, 0.376118859788, 0.369097692334, 00258 0.362142585282, 0.355252328834, 0.348425768415, 0.341661801776, 00259 0.334959376311, 0.328317486588, 0.321735172063, 0.31521151497, 00260 0.308745638367, 0.302336704338, 0.29598391232, 0.289686497571, 00261 0.283443729739, 0.27725491156, 0.271119377649, 0.265036493387, 00262 0.259005653912, 0.253026283183, 0.247097833139, 0.241219782932, 00263 0.235391638239, 0.229612930649, 0.223883217122, 0.218202079518, 00264 0.212569124201, 0.206983981709, 0.201446306496, 0.195955776745, 00265 0.190512094256, 0.185114984406, 0.179764196185, 0.174459502324, 00266 0.169200699492, 0.1639876086, 0.158820075195, 0.153697969964, 00267 0.148621189348, 0.143589656295, 0.138603321143, 0.133662162669, 00268 0.128766189309, 0.123915440582, 0.119109988745, 0.114349940703, 00269 0.10963544023, 0.104966670533, 0.100343857232, 0.0957672718266, 00270 0.0912372357329, 0.0867541250127, 0.082318375932, 0.0779304915295, 00271 0.0735910494266, 0.0693007111742, 0.065060233529, 0.0608704821745, 00272 0.056732448584, 0.05264727098, 0.0486162607163, 0.0446409359769, 00273 0.0407230655415, 0.0368647267386, 0.0330683839378, 0.0293369977411, 00274 0.0256741818288, 0.0220844372634, 0.0185735200577, 0.0151490552854, 00275 0.0118216532614, 0.00860719483079, 0.00553245272614, 0.00265435214565 00276 }; 00277 00278 /* 00279 * (Ziggurat) tabulated values for 2^24 times x[i]/x[i+1], used to accept 00280 * for U*x[i+1]<=x[i] without any floating point operations 00281 */ 00282 const unsigned int Normal_RNG::ktab[128] = { 00283 0, 12590644, 14272653, 14988939, 00284 15384584, 15635009, 15807561, 15933577, 00285 16029594, 16105155, 16166147, 16216399, 00286 16258508, 16294295, 16325078, 16351831, 00287 16375291, 16396026, 16414479, 16431002, 00288 16445880, 16459343, 16471578, 16482744, 00289 16492970, 16502368, 16511031, 16519039, 00290 16526459, 16533352, 16539769, 16545755, 00291 16551348, 16556584, 16561493, 16566101, 00292 16570433, 16574511, 16578353, 16581977, 00293 16585398, 16588629, 16591685, 16594575, 00294 16597311, 16599901, 16602354, 16604679, 00295 16606881, 16608968, 16610945, 16612818, 00296 16614592, 16616272, 16617861, 16619363, 00297 16620782, 16622121, 16623383, 16624570, 00298 16625685, 16626730, 16627708, 16628619, 00299 16629465, 16630248, 16630969, 16631628, 00300 16632228, 16632768, 16633248, 16633671, 00301 16634034, 16634340, 16634586, 16634774, 00302 16634903, 16634972, 16634980, 16634926, 00303 16634810, 16634628, 16634381, 16634066, 00304 16633680, 16633222, 16632688, 16632075, 00305 16631380, 16630598, 16629726, 16628757, 00306 16627686, 16626507, 16625212, 16623794, 00307 16622243, 16620548, 16618698, 16616679, 00308 16614476, 16612071, 16609444, 16606571, 00309 16603425, 16599973, 16596178, 16591995, 00310 16587369, 16582237, 16576520, 16570120, 00311 16562917, 16554758, 16545450, 16534739, 00312 16522287, 16507638, 16490152, 16468907, 00313 16442518, 16408804, 16364095, 16301683, 00314 16207738, 16047994, 15704248, 15472926 00315 }; 00316 00317 // (Ziggurat) tabulated values of 2^{-24}*x[i] 00318 const double Normal_RNG::wtab[128] = { 00319 1.62318314817e-08, 2.16291505214e-08, 2.54246305087e-08, 2.84579525938e-08, 00320 3.10340022482e-08, 3.33011726243e-08, 3.53439060345e-08, 3.72152672658e-08, 00321 3.8950989572e-08, 4.05763964764e-08, 4.21101548915e-08, 4.35664624904e-08, 00322 4.49563968336e-08, 4.62887864029e-08, 4.75707945735e-08, 4.88083237257e-08, 00323 5.00063025384e-08, 5.11688950428e-08, 5.22996558616e-08, 5.34016475624e-08, 00324 5.44775307871e-08, 5.55296344581e-08, 5.65600111659e-08, 5.75704813695e-08, 00325 5.85626690412e-08, 5.95380306862e-08, 6.04978791776e-08, 6.14434034901e-08, 00326 6.23756851626e-08, 6.32957121259e-08, 6.42043903937e-08, 6.51025540077e-08, 00327 6.59909735447e-08, 6.68703634341e-08, 6.77413882848e-08, 6.8604668381e-08, 00328 6.94607844804e-08, 7.03102820203e-08, 7.11536748229e-08, 7.1991448372e-08, 00329 7.2824062723e-08, 7.36519550992e-08, 7.44755422158e-08, 7.52952223703e-08, 00330 7.61113773308e-08, 7.69243740467e-08, 7.77345662086e-08, 7.85422956743e-08, 00331 7.93478937793e-08, 8.01516825471e-08, 8.09539758128e-08, 8.17550802699e-08, 00332 8.25552964535e-08, 8.33549196661e-08, 8.41542408569e-08, 8.49535474601e-08, 00333 8.57531242006e-08, 8.65532538723e-08, 8.73542180955e-08, 8.8156298059e-08, 00334 8.89597752521e-08, 8.97649321908e-08, 9.05720531451e-08, 9.138142487e-08, 00335 9.21933373471e-08, 9.30080845407e-08, 9.38259651738e-08, 9.46472835298e-08, 00336 9.54723502847e-08, 9.63014833769e-08, 9.71350089201e-08, 9.79732621669e-08, 00337 9.88165885297e-08, 9.96653446693e-08, 1.00519899658e-07, 1.0138063623e-07, 00338 1.02247952126e-07, 1.03122261554e-07, 1.04003996769e-07, 1.04893609795e-07, 00339 1.05791574313e-07, 1.06698387725e-07, 1.07614573423e-07, 1.08540683296e-07, 00340 1.09477300508e-07, 1.1042504257e-07, 1.11384564771e-07, 1.12356564007e-07, 00341 1.13341783071e-07, 1.14341015475e-07, 1.15355110887e-07, 1.16384981291e-07, 00342 1.17431607977e-07, 1.18496049514e-07, 1.19579450872e-07, 1.20683053909e-07, 00343 1.21808209468e-07, 1.2295639141e-07, 1.24129212952e-07, 1.25328445797e-07, 00344 1.26556042658e-07, 1.27814163916e-07, 1.29105209375e-07, 1.30431856341e-07, 00345 1.31797105598e-07, 1.3320433736e-07, 1.34657379914e-07, 1.36160594606e-07, 00346 1.37718982103e-07, 1.39338316679e-07, 1.41025317971e-07, 1.42787873535e-07, 00347 1.44635331499e-07, 1.4657889173e-07, 1.48632138436e-07, 1.50811780719e-07, 00348 1.53138707402e-07, 1.55639532047e-07, 1.58348931426e-07, 1.61313325908e-07, 00349 1.64596952856e-07, 1.68292495203e-07, 1.72541128694e-07, 1.77574279496e-07, 00350 1.83813550477e-07, 1.92166040885e-07, 2.05295471952e-07, 2.22600839893e-07 00351 }; 00352 00353 // (Ziggurat) position of right-most step 00354 const double Normal_RNG::PARAM_R = 3.44428647676; 00355 00356 // Get a Normal distributed (0,1) sample 00357 double Normal_RNG::sample() 00358 { 00359 unsigned int u, sign, i, j; 00360 double x, y; 00361 00362 while (true) { 00363 u = RNG.random_int(); 00364 sign = u & 0x80; // 1 bit for the sign 00365 i = u & 0x7f; // 7 bits to choose the step 00366 j = u >> 8; // 24 bits for the x-value 00367 00368 x = j * wtab[i]; 00369 00370 if (j < ktab[i]) 00371 break; 00372 00373 if (i < 127) { 00374 y = ytab[i+1] + (ytab[i] - ytab[i+1]) * RNG.random_01(); 00375 } 00376 else { 00377 x = PARAM_R - std::log(1.0 - RNG.random_01()) / PARAM_R; 00378 y = std::exp(-PARAM_R * (x - 0.5 * PARAM_R)) * RNG.random_01(); 00379 } 00380 00381 if (y < std::exp(-0.5 * x * x)) 00382 break; 00383 } 00384 return sign ? x : -x; 00385 } 00386 00387 00389 // Laplace_RNG 00391 00392 Laplace_RNG::Laplace_RNG(double meanval, double variance) 00393 { 00394 setup(meanval, variance); 00395 } 00396 00397 void Laplace_RNG::setup(double meanval, double variance) 00398 { 00399 mean = meanval; 00400 var = variance; 00401 sqrt_12var = std::sqrt(variance / 2.0); 00402 } 00403 00404 void Laplace_RNG::get_setup(double &meanval, double &variance) const 00405 { 00406 meanval = mean; 00407 variance = var; 00408 } 00409 00410 00411 00412 vec Laplace_RNG::operator()(int n) 00413 { 00414 vec vv(n); 00415 00416 for (int i=0; i < n; i++) 00417 vv(i) = sample(); 00418 00419 return vv; 00420 } 00421 00422 mat Laplace_RNG::operator()(int h, int w) 00423 { 00424 mat mm(h, w); 00425 int i, j; 00426 00427 for (i = 0; i < h; i++) 00428 for (j = 0; j < w; j++) 00429 mm(i, j) = sample(); 00430 00431 return mm; 00432 } 00433 00435 // AR1_Normal_RNG 00437 00438 AR1_Normal_RNG::AR1_Normal_RNG(double meanval, double variance, double rho) 00439 { 00440 setup(meanval, variance, rho); 00441 } 00442 00443 void AR1_Normal_RNG::setup(double meanval, double variance, double rho) 00444 { 00445 mean = meanval; 00446 var = variance; 00447 r = rho; 00448 factr = -2.0 * var * (1.0 - rho * rho); 00449 mem = 0.0; 00450 odd = true; 00451 } 00452 00453 void AR1_Normal_RNG::get_setup(double &meanval, double &variance, 00454 double &rho) const 00455 { 00456 meanval = mean; 00457 variance = var; 00458 rho = r; 00459 } 00460 00461 vec AR1_Normal_RNG::operator()(int n) 00462 { 00463 vec vv(n); 00464 00465 for (int i=0; i < n; i++) 00466 vv(i) = sample(); 00467 00468 return vv; 00469 } 00470 00471 mat AR1_Normal_RNG::operator()(int h, int w) 00472 { 00473 mat mm(h, w); 00474 int i, j; 00475 00476 for (i = 0; i < h; i++) 00477 for (j = 0; j < w; j++) 00478 mm(i, j) = sample(); 00479 00480 return mm; 00481 } 00482 00483 void AR1_Normal_RNG::reset() 00484 { 00485 mem = 0.0; 00486 } 00487 00489 // Weibull_RNG 00491 00492 Weibull_RNG::Weibull_RNG(double lambda, double beta) 00493 { 00494 setup(lambda, beta); 00495 } 00496 00497 void Weibull_RNG::setup(double lambda, double beta) 00498 { 00499 l = lambda; 00500 b = beta; 00501 mean = tgamma(1.0 + 1.0 / b) / l; 00502 var = tgamma(1.0 + 2.0 / b) / (l * l) - mean; 00503 } 00504 00505 00506 vec Weibull_RNG::operator()(int n) 00507 { 00508 vec vv(n); 00509 00510 for (int i=0; i < n; i++) 00511 vv(i) = sample(); 00512 00513 return vv; 00514 } 00515 00516 mat Weibull_RNG::operator()(int h, int w) 00517 { 00518 mat mm(h, w); 00519 int i, j; 00520 00521 for (i = 0; i < h; i++) 00522 for (j = 0; j < w; j++) 00523 mm(i, j) = sample(); 00524 00525 return mm; 00526 } 00527 00529 // Rayleigh_RNG 00531 00532 Rayleigh_RNG::Rayleigh_RNG(double sigma) 00533 { 00534 setup(sigma); 00535 } 00536 00537 vec Rayleigh_RNG::operator()(int n) 00538 { 00539 vec vv(n); 00540 00541 for (int i=0; i < n; i++) 00542 vv(i) = sample(); 00543 00544 return vv; 00545 } 00546 00547 mat Rayleigh_RNG::operator()(int h, int w) 00548 { 00549 mat mm(h, w); 00550 int i, j; 00551 00552 for (i = 0; i < h; i++) 00553 for (j = 0; j < w; j++) 00554 mm(i, j) = sample(); 00555 00556 return mm; 00557 } 00558 00560 // Rice_RNG 00562 00563 Rice_RNG::Rice_RNG(double lambda, double beta) 00564 { 00565 setup(lambda, beta); 00566 } 00567 00568 vec Rice_RNG::operator()(int n) 00569 { 00570 vec vv(n); 00571 00572 for (int i=0; i < n; i++) 00573 vv(i) = sample(); 00574 00575 return vv; 00576 } 00577 00578 mat Rice_RNG::operator()(int h, int w) 00579 { 00580 mat mm(h, w); 00581 int i, j; 00582 00583 for (i = 0; i < h; i++) 00584 for (j = 0; j < w; j++) 00585 mm(i, j) = sample(); 00586 00587 return mm; 00588 } 00589 00590 } // namespace itpp
Generated on Tue Dec 6 2011 16:51:48 for IT++ by Doxygen 1.7.4