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

../trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp

@@ -8 +8 @@
 #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 +17 @@
         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 +37 @@
   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();

../trunk-jpl/src/c/shared/Random/randomgenerator.cpp

@@ -12 +12 @@
 
 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));
 
         }

../trunk-jpl/src/c/shared/Random/randomgenerator.h

@@ -10 +10 @@
 
 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 +53 @@
   {
 
     private:
-      unsigned _seed;
       double mean;
       double sdev;
 
@@ -49 +63 @@
       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_ */