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

Generated on Sun Sep 14 18:57:02 2008 for IT++ by Doxygen 1.5.6