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
Generated on Sat Aug 25 23:40:03 2007 for IT++ by Doxygen 1.5.2