Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp (revision 26611) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp (revision 26612) @@ -8,7 +8,7 @@ #include "../modules/modules.h" #include "../analyses/analyses.h" #include "../shared/Random/randomgenerator.h" -#include +//#include void GaussianVector(Vector* ppf,int seed){/*{{{*/ @@ -17,19 +17,19 @@ double rdnumber; /*Define seed*/ - rnd::normal_distribution distribution; - if(seed<0){ - std::random_device rd; - seed = rd(); - } - else + rnd::linear_congruential_engine random_engine; + if(seed>=0) { int my_rank; ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&my_rank); seed = seed + 783728*my_rank; // change default seed for parallel simulations (by considering an arbitrary shif based on the rank number) } - distribution.seed(seed); + random_engine.seed(seed); + /* Define univariate distribution */ + + rnd::normal_distribution distribution(0.0,1.0); + int *local_indices = NULL; IssmDouble *local_vector = NULL; ppf->GetLocalVector(&local_vector,&local_indices); @@ -37,7 +37,7 @@ int M; ppf->GetLocalSize(&M); for(int i=0;iSetValue(local_indices[i],rdnumber,INS_VAL); } ppf->Assemble(); Index: ../trunk-jpl/src/c/shared/Random/randomgenerator.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Random/randomgenerator.cpp (revision 26611) +++ ../trunk-jpl/src/c/shared/Random/randomgenerator.cpp (revision 26612) @@ -12,68 +12,94 @@ namespace rnd{ - uniform_distribution::uniform_distribution(){/*{{{*/ + /* Linear congruential engine */ - a = 1103515245; // BSD Formula - c = 12345; // BSD Formula - m = 2147483648; // BSD Formula - _seed = 0; - lbound = 0.0; - ubound = 1.0; + linear_congruential_engine::linear_congruential_engine(){/*{{{*/ + + pseed = new unsigned int; + *pseed = std::chrono::steady_clock::now().time_since_epoch()/std::chrono::milliseconds(1); + a = 1103515245; // BSD Formula + c = 12345; // BSD Formula + m = 2147483648; // BSD Formula return; + }/*}}}*/ + + linear_congruential_engine::linear_congruential_engine(unsigned int _a, unsigned int _b, unsigned int _m){/*{{{*/ + + pseed = new unsigned int; + *pseed = std::chrono::steady_clock::now().time_since_epoch()/std::chrono::milliseconds(1); + a = _a; + c = _b; + m = _m; + return; + }/*}}}*/ + + linear_congruential_engine::~linear_congruential_engine(){} + + unsigned int linear_congruential_engine::get_m() { return m; } + + void linear_congruential_engine::seed( int s ) { + if(s<0) + *pseed = std::chrono::steady_clock::now().time_since_epoch()/std::chrono::milliseconds(1); + else + *pseed = (unsigned) s; + } + + unsigned int linear_congruential_engine::generator(){ + *pseed = ( a * *pseed + c ) % m ; + return *pseed; + } + + /* Uniform distribution */ + + uniform_distribution::uniform_distribution(){/*{{{*/ + a = 0.0; + b = 1.0; + return; } /*}}}*/ - uniform_distribution::uniform_distribution(double lower,double upper){/*{{{*/ - a = 1103515245; // BSD Formula - c = 12345; // BSD Formula - m = 2147483648; // BSD Formula - _seed = 0; - lbound = lower; - ubound = upper; + uniform_distribution::uniform_distribution(double _a,double _b){/*{{{*/ + a = _a; + b = _b; return; } /*}}}*/ + uniform_distribution::~uniform_distribution(){} - void uniform_distribution::seed( unsigned int s ) { _seed = s; } - unsigned int uniform_distribution::get_seed() { return _seed; } - double uniform_distribution::generator() { - _seed = ( a * _seed + c ) % m ; - return (ubound-lbound)*(double) _seed/ m + lbound; + + double uniform_distribution::generator(rnd::linear_congruential_engine random_engine) { + return (b-a)*(double) random_engine.generator()/ random_engine.get_m() + a; } + /* Normal distribution */ normal_distribution::normal_distribution(){/*{{{*/ - - _seed = 0; mean = 0; sdev = 1.0; return; } /*}}}*/ + normal_distribution::normal_distribution(double m,double s){/*{{{*/ - - _seed = 0; mean = m; sdev = s; return; } - /*}}}*/ + /*}}}*/ + normal_distribution::~normal_distribution(){} - void normal_distribution::seed( unsigned int s ) { _seed = s; } - double normal_distribution::generator(){/*{{{*/ + double normal_distribution::generator(rnd::linear_congruential_engine random_engine){/*{{{*/ + rnd::uniform_distribution unifdistri; - unifdistri.seed(_seed); - double u1 = unifdistri.generator(); - double u2 = unifdistri.generator(); + double u1 = unifdistri.generator(random_engine); + double u2 = unifdistri.generator(random_engine); double R = sqrt(-2*log(u1)); double theta = 2*M_PI*u2; - seed(unifdistri.get_seed()); - return mean + sdev * (R*cos(theta)); } Index: ../trunk-jpl/src/c/shared/Random/randomgenerator.h =================================================================== --- ../trunk-jpl/src/c/shared/Random/randomgenerator.h (revision 26611) +++ ../trunk-jpl/src/c/shared/Random/randomgenerator.h (revision 26612) @@ -10,27 +10,42 @@ namespace rnd{ + class linear_congruential_engine + { + private: + unsigned int a; + unsigned int c; + unsigned int m; + unsigned int *pseed; + + public: + + /*constructors, destructors: */ + linear_congruential_engine(); + linear_congruential_engine(unsigned int _a, unsigned int _b, unsigned int _m); + ~linear_congruential_engine(); + + unsigned int get_m(); + void seed( int s ); + unsigned int generator(); + + }; + class uniform_distribution { private: - int a; - int c; - unsigned int m; - unsigned _seed; - double lbound; - double ubound; + double a; // lower bound of range + double b; // upper bound of range public: /*constructors, destructors: */ uniform_distribution(); - uniform_distribution(double a_1,double a_2); + uniform_distribution(double _a,double _b); ~uniform_distribution(); - void seed( unsigned int s ); - unsigned int get_seed(); - double generator(); + double generator(rnd::linear_congruential_engine random_engine); }; @@ -38,7 +53,6 @@ { private: - unsigned _seed; double mean; double sdev; @@ -49,11 +63,11 @@ normal_distribution(double m,double s); ~normal_distribution(); - void seed( unsigned int s ); - double generator(); + double generator(rnd::linear_congruential_engine random_engine); }; + } -#endif //ifndef _RANDOMGENERATOR_H_ +#endif //* _RANDOMGENERATOR_H_ */