Changeset 26554

Timestamp:

11/08/21 14:41:51 (3 years ago)

Author:

bulthuis

Message:

NEW: add pseudo-random number generators + use in solutionsequence_sampling.cpp

Location:

issm/trunk-jpl

Files:

• src/c/shared/Random/randomgenerator.cpp (added)
• src/c/shared/Random/randomgenerator.h (added)
• src/c/solutionsequences/solutionsequence_sampling.cpp (modified) (4 diffs)
• test/Archives/Archive134.arch (modified) ( previous)

issm/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>
-
-#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){/*{{{*/
@@ -49 +18 @@
 
         /*Define seed*/
+        normal_distribution_rnd distribution;
         if(seed<0){
                 std::random_device rd;
@@ -59 +29 @@
                 seed = seed + 783728*my_rank; // change default seed for parallel simulations (by considering an arbitrary shif based on the rank number)
         }
-
-        long long seedbis;
-        seedbis  = (long long) seed;
+        distribution.seed(seed);
 
         int        *local_indices = NULL;
@@ -70 +38 @@
         ppf->GetLocalSize(&M);
         for(int i=0;i<M;i++){
-                rdnumber = BoxMuller(&seedbis);
+                rdnumber = distribution.rnd();
                 ppf->SetValue(local_indices[i],rdnumber,INS_VAL);
         }