IT++ Logo

random.cpp

Go to the documentation of this file.
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
SourceForge Logo

Generated on Thu Apr 23 20:06:40 2009 for IT++ by Doxygen 1.5.8