Changeset 26621


Ignore:
Timestamp:
11/13/21 07:25:07 (3 years ago)
Author:
vverjans
Message:

CHG: making random number generation in ISSM dependent on classes of randomgenerator.h (should allow cross-platform tests of stochasticity)

Location:
issm/trunk-jpl
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/shared/Random/random.cpp

    r26483 r26621  
    99#include <math.h>
    1010#include <float.h>    /*  DBL_EPSILON  */
    11 #include <chrono>
    1211#include <cstdarg>
    1312#include <iostream>
    14 #include <random>
    1513
    1614#include "../Matrix/matrix.h"
    … …  
    1816#include "../MemOps/MemOps.h"
    1917#include "../io/io.h"
     18#include "./randomgenerator.h"
    2019/*}}}*/
    2120
    22 void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev, int seedfixed=-1) { /*{{{*/
     21void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev, int seed=-1) { /*{{{*/
    2322
    24         unsigned seed;
    25         /*Random seed using time_since_epoch*/
    26         if(seedfixed==-1) seed = std::chrono::steady_clock::now().time_since_epoch().count();
    27         /*Seed fixed by input argument*/
    28         else seed = seedfixed;
    29         std::default_random_engine generator(seed);
    30         /*Normal Probability Distribution*/
    31    std::normal_distribution<IssmPDouble> normdistri(mean,sdev);
    32         *prand = normdistri(generator);
     23        /*Seed the pseudo-random number generator*/
     24        rnd::linear_congruential_engine randomengine;
     25        randomengine.seed(seed);
     26        /*Normal distribution*/
     27        rnd::normal_distribution distriNormal(mean,sdev);
     28        *prand = distriNormal.generator(randomengine);
    3329} /*}}}*/
    34 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seedfixed=-1) { /*{{{*/
     30void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seed=-1) { /*{{{*/
    3531   
    3632        IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
    … …  
    3834   IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
    3935
    40         /*True randomness if seedfixed==-1, otherwise random seed is fixed at seedfixed*/
    41         for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0,seedfixed);
     36        /*True randomness if seed<0, otherwise random seed is fixed at seed*/
     37        /*Seed the pseudo-random number generator, repeatedly calling univariateNormal does not ensure randomness*/
     38        rnd::linear_congruential_engine randomengine;
     39        randomengine.seed(seed);
     40        /*Normal distribution*/
     41        rnd::normal_distribution distriNormal(0.0,1.0);
     42        for(int i=0;i<dim;i++){
     43                sampleStandardNormal[i] = distriNormal.generator(randomengine);
     44        }
     45
     46        /*Cholsesky decomposition of the covariance matrix*/
    4247        CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
    4348   
    … …  
    5560   xDelete<IssmDouble>(Lchol);
    5661} /*}}}*/
    57 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seedfixed=-1) { /*{{{*/
     62void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seed=-1) { /*{{{*/
    5863       
    5964        IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
    6065        IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
    6166        IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
     67       
     68        /*True randomness if seed<0, otherwise random seed is fixed at seed*/
     69        /*Seed the pseudo-random number generator, repeatedly calling univariateNormal does not ensure randomness*/
     70        rnd::linear_congruential_engine randomengine;
     71        randomengine.seed(seed);
     72        /*Normal distribution*/
     73        rnd::normal_distribution distriNormal(0.0,1.0);
     74        for(int i=0;i<dim;i++){
     75                sampleStandardNormal[i] = distriNormal.generator(randomengine);
     76        }
    6277
    63         /*True randomness if seedfixed==-1, otherwise random seed is fixed at seedfixed*/
    64         for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0,seedfixed);
     78        /*Cholsesky decomposition of the covariance matrix*/
    6579        CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
    6680

  • issm/trunk-jpl/src/c/shared/Random/random.h

    r26483 r26621  
    66#define _RANDOM_H_
    77
    8 void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev, int seedfixed=-1);
    9 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seedfixed=-1);
    10 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seedfixed=-1);
     8void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev, int seed=-1);
     9void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seed=-1);
     10void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seed=-1);
    1111
    1212#endif //ifndef _RANDOM_H_
Note: See TracChangeset for help on using the changeset viewer.