IT++ Logo Newcom Logo

random.cpp

Go to the documentation of this file.
00001 
00033 #include <itpp/base/random.h>
00034 #include <limits>
00035 
00036 
00037 namespace itpp { 
00038 
00039   int Random_Generator::left = 0;
00040   unsigned int Random_Generator::state[624];
00041   unsigned int *Random_Generator::pNext;
00042 
00044   static bool __Random_Generator_seed_is_initialized = false; 
00045 
00047   // Random_Generator
00049 
00050   Random_Generator::Random_Generator() : lastSeed( 4357U )
00051   {
00052     if (!__Random_Generator_seed_is_initialized){
00053       reset();
00054       __Random_Generator_seed_is_initialized = true;
00055     }
00056   }
00057 
00058   unsigned int Random_Generator::hash( time_t t, clock_t c )
00059   {
00060     // Get a unsigned int from t and c
00061     // Better than uint(x) in case x is floating point in [0,1]
00062     // Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
00063     static unsigned int differ = 0; // guarantee time-based seeds will change
00064 
00065     unsigned int h1 = 0;
00066     unsigned char *p = (unsigned char *) &t;
00067     for( size_t i = 0; i < sizeof(t); ++i )
00068       {
00069         h1 *= std::numeric_limits<unsigned char>::max() + 2U;
00070         h1 += p[i];
00071       }
00072     unsigned int h2 = 0;
00073     p = (unsigned char *) &c;
00074     for( size_t j = 0; j < sizeof(c); ++j )
00075       {
00076         h2 *= std::numeric_limits<unsigned char>::max() + 2U;
00077         h2 += p[j];
00078       }
00079     return ( h1 + differ++ ) ^ h2;
00080   }
00081 
00082   void Random_Generator::randomize()
00083   {
00084     lastSeed = hash(time(0), clock());
00085     reset();
00086   }
00087 
00088   void Random_Generator::get_state(ivec &out_state)
00089   {
00090     out_state.set_size(625, false);
00091     for (int i=0; i<624; i++)
00092       out_state(i) = state[i];
00093 
00094     out_state(624) = left; // the number of elements left in state before reload
00095   }
00096 
00097   void Random_Generator::set_state(ivec &new_state)
00098   {
00099     it_assert(new_state.size()==625, "Random_Generator::set_state(): Not a valid state vector");
00100 
00101     for (int i=0; i<624; i++)
00102       state[i] = new_state(i);
00103 
00104     left = new_state(624);
00105     pNext = &state[624-left];
00106   }
00107 
00108 
00109   // Set the seed of the Global Random Number Generator
00110   void RNG_reset(unsigned long seed)
00111   {
00112     Random_Generator RNG;
00113     RNG.reset(seed);
00114   }
00115 
00116   // Set the seed of the Global Random Number Generator to the same as last time
00117   void RNG_reset()
00118   {
00119     Random_Generator RNG;
00120     RNG.reset();
00121   }
00122 
00123   // Set a random seed for the Global Random Number Generator
00124   void RNG_randomize()
00125   {
00126     Random_Generator RNG;
00127     RNG.randomize();
00128   }
00129 
00130   // Save current full state of generator in memory
00131   void RNG_get_state(ivec &state)
00132   {
00133     Random_Generator RNG;
00134     RNG.get_state(state);
00135   }
00136 
00137   // Resume the state saved in memory
00138   void RNG_set_state(ivec &state)
00139   {
00140     Random_Generator RNG;
00141     RNG.set_state(state);
00142   }
00143 
00145   // I_Uniform_RNG
00147 
00148   I_Uniform_RNG::I_Uniform_RNG(int min, int max)
00149   {
00150     setup(min, max);
00151   }
00152 
00153   void I_Uniform_RNG::setup(int min, int max)
00154   {
00155     if (min <= max) {
00156       lo = min;
00157       hi = max;
00158     }
00159     else {
00160       lo = max;
00161       hi = min;
00162     }
00163   }
00164 
00165   void I_Uniform_RNG::get_setup(int &min, int &max) const
00166   {
00167     min = lo;
00168     max = hi;
00169   }
00170 
00171   ivec I_Uniform_RNG::operator()(int n)
00172   {
00173     ivec vv(n);
00174 
00175     for (int i=0; i<n; i++)
00176       vv(i) = sample();
00177 
00178     return vv;
00179   }
00180 
00181   imat I_Uniform_RNG::operator()(int h, int w)
00182   {
00183     imat mm(h,w);
00184     int i,j;
00185 
00186     for (i=0; i<h; i++)
00187       for (j=0; j<w; j++)
00188         mm(i,j) = sample();
00189 
00190     return mm;
00191   }
00192 
00194   // Uniform_RNG
00196 
00197   Uniform_RNG::Uniform_RNG(double min, double max)
00198   {
00199     setup(min, max);
00200   }
00201 
00202   void Uniform_RNG::setup(double min, double max)
00203   {
00204     if (min <= max) {
00205       lo_bound = min;
00206       hi_bound = max;
00207     }
00208     else {
00209       lo_bound = max;
00210       hi_bound = min;
00211     }
00212   }
00213 
00214   void Uniform_RNG::get_setup(double &min, double &max) const
00215   {
00216     min = lo_bound;
00217     max = hi_bound;
00218   }
00219 
00221   // Exp_RNG
00223 
00224   Exponential_RNG::Exponential_RNG(double lambda)
00225   {
00226     setup(lambda);
00227   }
00228 
00229   vec Exponential_RNG::operator()(int n)
00230   {
00231     vec vv(n);
00232 
00233     for (int i=0; i<n; i++)
00234       vv(i) = sample();
00235 
00236     return vv;
00237   }
00238 
00239   mat Exponential_RNG::operator()(int h, int w)
00240   {
00241     mat mm(h,w);
00242     int i,j;
00243 
00244     for (i=0; i<h; i++)
00245       for (j=0; j<w; j++)
00246         mm(i,j) = sample();
00247 
00248     return mm;
00249   }
00250 
00252   // Normal_RNG
00254 
00255   void Normal_RNG::get_setup(double &meanval, double &variance) const
00256   {
00257     meanval = mean;
00258     variance = sigma*sigma;
00259   }
00260 
00262   // Laplace_RNG
00264 
00265   Laplace_RNG::Laplace_RNG(double meanval, double variance)
00266   {
00267     setup(meanval, variance);
00268   }
00269 
00270   void Laplace_RNG::setup(double meanval, double variance)
00271   {
00272     mean = meanval;
00273     var = variance;
00274   }
00275 
00276   void Laplace_RNG::get_setup(double &meanval, double &variance) const
00277   {
00278     meanval = mean;
00279     variance = var;
00280   }
00281 
00282 
00283 
00284   vec Laplace_RNG::operator()(int n)
00285   {
00286     vec vv(n);
00287 
00288     for (int i=0; i<n; i++)
00289       vv(i) = sample();
00290 
00291     return vv;
00292   }
00293 
00294   mat Laplace_RNG::operator()(int h, int w)
00295   {
00296     mat mm(h,w);
00297     int i,j;
00298 
00299     for (i=0; i<h; i++)
00300       for (j=0; j<w; j++)
00301         mm(i,j) = sample();
00302 
00303     return mm;
00304   }
00305 
00307   // AR1_Normal_RNG
00309 
00310   AR1_Normal_RNG::AR1_Normal_RNG(double meanval, double variance, double rho)
00311   {
00312     mean = meanval;
00313     var = variance;
00314     r = rho;
00315     mem = 0.0;
00316     factr = -2.0 * var * (1.0 - rho*rho);
00317     odd = true;
00318   }
00319 
00320   void AR1_Normal_RNG::setup(double meanval, double variance, double rho)
00321   {
00322     mean = meanval;
00323     var = variance;
00324     r = rho;
00325     factr = -2.0 * var * (1.0 - rho*rho);
00326     mem = 0.0;
00327     odd = true;
00328   }
00329 
00330   void AR1_Normal_RNG::get_setup(double &meanval, double &variance, double &rho) const
00331   {
00332     meanval = mean;
00333     variance = var;
00334     rho = r;
00335   }
00336 
00337   vec AR1_Normal_RNG::operator()(int n)
00338   {
00339     vec vv(n);
00340 
00341     for (int i=0; i<n; i++)
00342       vv(i) = sample();
00343 
00344     return vv;
00345   }
00346 
00347   mat AR1_Normal_RNG::operator()(int h, int w)
00348   {
00349     mat mm(h,w);
00350     int i,j;
00351 
00352     for (i=0; i<h; i++)
00353       for (j=0; j<w; j++)
00354         mm(i,j) = sample();
00355 
00356     return mm;
00357   }
00358 
00359   void AR1_Normal_RNG::reset()
00360   {
00361     mem = 0.0;
00362   }
00363 
00365   // Weibull_RNG
00367 
00368   Weibull_RNG::Weibull_RNG(double lambda, double beta)
00369   {
00370     setup(lambda, beta);
00371   }
00372 
00373   void Weibull_RNG::setup(double lambda, double beta)
00374   {
00375     l=lambda;
00376     b=beta;
00377     mean = gamma(1.0 + 1.0 / b) / l;
00378     var = gamma(1.0+2.0/b)/(l*l) - mean;
00379   }
00380 
00381 
00382   vec Weibull_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 Weibull_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   // Rayleigh_RNG
00407 
00408   Rayleigh_RNG::Rayleigh_RNG(double sigma)
00409   {
00410     setup(sigma);
00411   }
00412 
00413   vec Rayleigh_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 Rayleigh_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   // Rice_RNG
00438 
00439   Rice_RNG::Rice_RNG(double lambda, double beta)
00440   {
00441     setup(lambda, beta);
00442   }
00443 
00444   vec Rice_RNG::operator()(int n)
00445   {
00446     vec vv(n);
00447 
00448     for (int i=0; i<n; i++)
00449       vv(i) = sample();
00450 
00451     return vv;
00452   }
00453 
00454   mat Rice_RNG::operator()(int h, int w)
00455   {
00456     mat mm(h,w);
00457     int i,j;
00458 
00459     for (i=0; i<h; i++)
00460       for (j=0; j<w; j++)
00461         mm(i,j) = sample();
00462 
00463     return mm;
00464   }
00465 
00466 } // namespace itpp
SourceForge Logo

Generated on Sat Aug 25 23:40:53 2007 for IT++ by Doxygen 1.5.2