source: issm/oecreview/Archive/25834-26739/ISSM-26611-26612.diff@ 26740

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 <random>
+//#include <random>
 
 void GaussianVector(Vector<IssmDouble>* 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;i<M;i++){
-               rdnumber = distribution.generator();
+               rdnumber = distribution.generator(random_engine);
                ppf->SetValue(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_ */