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