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 
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
SourceForge Logo

Generated on Sat Apr 19 10:42:05 2008 for IT++ by Doxygen 1.5.5