Changeset 26514

Timestamp:

10/27/21 21:43:46 (3 years ago)

Author:

bulthuis

Message:

CHG: Implement a pseudo-random number generator for Gaussian variables in solutionsequence_sampling.cpp

File:

• issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp (modified) (3 diffs)

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

@@ -10 +10 @@
 #include <random>
 
+#undef M_PI
+#define M_PI 3.141592653589793238462643
+
+double lcg(long long* seed){/*{{{*/
+
+
+        int N = 5;
+        long long M = 1LL<<32;
+        long long a = 1812433253;
+        long long c = 1;
+
+        for(int i=0;i<N;i++){
+                *seed = a * *seed + c;
+                *seed = *seed%M;
+        }
+
+  return ((double) *seed) / ((double) M);
+
+}
+
+double BoxMuller(long long* pseed){
+
+        double u1 = lcg(pseed);
+        double u2 = lcg(pseed);
+
+        double R = sqrt(-2*log(u1));
+        double theta = 2*M_PI*u2;
+
+        return R*cos(theta);
+
+}
+
 void GaussianVector(Vector<IssmDouble>* ppf,int seed){/*{{{*/
+
 
         /*Intermediaries*/
         double      rdnumber;
 
-        /*Define random number generator*/
-        std::normal_distribution<double> distribution(0.0,1.0);  // Define probability distribution as the standard normal distribution
+        /*Define seed*/
         if(seed<0){
                 std::random_device rd;
@@ -27 +59 @@
                 seed = seed + 783728*my_rank; // change default seed for parallel simulations (by considering an arbitrary shif based on the rank number)
         }
-        //std::default_random_engine generator(seed);
-        std::mt19937 generator(seed);
+
+        long long seedbis;
+        seedbis  = (long long) seed;
 
         int        *local_indices = NULL;
@@ -37 +70 @@
         ppf->GetLocalSize(&M);
         for(int i=0;i<M;i++){
-                rdnumber = distribution(generator);
+                rdnumber = BoxMuller(&seedbis);
                 ppf->SetValue(local_indices[i],rdnumber,INS_VAL);
         }