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 00359 // Laplace_RNG 00361 00362 Laplace_RNG::Laplace_RNG(double meanval, double variance) 00363 { 00364 setup(meanval, variance); 00365 } 00366 00367 void Laplace_RNG::setup(double meanval, double variance) 00368 { 00369 mean = meanval; 00370 var = variance; 00371 sqrt_12var = std::sqrt(variance / 2.0); 00372 } 00373 00374 void Laplace_RNG::get_setup(double &meanval, double &variance) const 00375 { 00376 meanval = mean; 00377 variance = var; 00378 } 00379 00380 00381 00382 vec Laplace_RNG::operator()(int n) 00383 { 00384 vec vv(n); 00385 00386 for (int i=0; i<n; i++) 00387 vv(i) = sample(); 00388 00389 return vv; 00390 } 00391 00392 mat Laplace_RNG::operator()(int h, int w) 00393 { 00394 mat mm(h,w); 00395 int i,j; 00396 00397 for (i=0; i<h; i++) 00398 for (j=0; j<w; j++) 00399 mm(i,j) = sample(); 00400 00401 return mm; 00402 } 00403 00405 // AR1_Normal_RNG 00407 00408 AR1_Normal_RNG::AR1_Normal_RNG(double meanval, double variance, double rho) 00409 { 00410 setup(meanval, variance, rho); 00411 } 00412 00413 void AR1_Normal_RNG::setup(double meanval, double variance, double rho) 00414 { 00415 mean = meanval; 00416 var = variance; 00417 r = rho; 00418 factr = -2.0 * var * (1.0 - rho*rho); 00419 mem = 0.0; 00420 odd = true; 00421 } 00422 00423 void AR1_Normal_RNG::get_setup(double &meanval, double &variance, 00424 double &rho) const 00425 { 00426 meanval = mean; 00427 variance = var; 00428 rho = r; 00429 } 00430 00431 vec AR1_Normal_RNG::operator()(int n) 00432 { 00433 vec vv(n); 00434 00435 for (int i=0; i<n; i++) 00436 vv(i) = sample(); 00437 00438 return vv; 00439 } 00440 00441 mat AR1_Normal_RNG::operator()(int h, int w) 00442 { 00443 mat mm(h,w); 00444 int i,j; 00445 00446 for (i=0; i<h; i++) 00447 for (j=0; j<w; j++) 00448 mm(i,j) = sample(); 00449 00450 return mm; 00451 } 00452 00453 void AR1_Normal_RNG::reset() 00454 { 00455 mem = 0.0; 00456 } 00457 00459 // Weibull_RNG 00461 00462 Weibull_RNG::Weibull_RNG(double lambda, double beta) 00463 { 00464 setup(lambda, beta); 00465 } 00466 00467 void Weibull_RNG::setup(double lambda, double beta) 00468 { 00469 l=lambda; 00470 b=beta; 00471 mean = gamma(1.0 + 1.0 / b) / l; 00472 var = gamma(1.0+2.0/b)/(l*l) - mean; 00473 } 00474 00475 00476 vec Weibull_RNG::operator()(int n) 00477 { 00478 vec vv(n); 00479 00480 for (int i=0; i<n; i++) 00481 vv(i) = sample(); 00482 00483 return vv; 00484 } 00485 00486 mat Weibull_RNG::operator()(int h, int w) 00487 { 00488 mat mm(h,w); 00489 int i,j; 00490 00491 for (i=0; i<h; i++) 00492 for (j=0; j<w; j++) 00493 mm(i,j) = sample(); 00494 00495 return mm; 00496 } 00497 00499 // Rayleigh_RNG 00501 00502 Rayleigh_RNG::Rayleigh_RNG(double sigma) 00503 { 00504 setup(sigma); 00505 } 00506 00507 vec Rayleigh_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 Rayleigh_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 // Rice_RNG 00532 00533 Rice_RNG::Rice_RNG(double lambda, double beta) 00534 { 00535 setup(lambda, beta); 00536 } 00537 00538 vec Rice_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 Rice_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 00560 } // namespace itpp
Generated on Sun Dec 9 17:31:00 2007 for IT++ by Doxygen 1.5.4