Changeset 26597


Ignore:
Timestamp:
11/10/21 15:55:39 (3 years ago)
Author:
bulthuis
Message:

BUG: Try to fix bugs pseudo-random number generators

Location:
issm/trunk-jpl
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r26596 r26597  
    22 * \brief random number generating functions
    33 */
    4 
    5 #include <iostream>
    6 #include "./randomgenerator.h"
    7 
    8 #undef M_PI
    9 #define M_PI 3.141592653589793238462643
    10 
    11 //rnd_uniform_distribution::rnd_uniform_distribution(){/*{{{*/
    12 
    13 //              a   = 1103515245;       // BSD Formula
    14 //              c  = 12345;                                     // BSD Formula
    15 //              m = 2147483648;                 // BSD Formula
    16 //              _seed = 0;
    17 //              lbound = 0.0;
    18 //              ubound = 1.0;
    19 //              return;
    20 //}
    21 /*}}}*/

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

    r26596 r26597  
    22 * \brief prototypes for randomgenerator.h
    33 */
    4 
    5 //#ifndef _RANDOMGENERATOR_H_
    6 //#define _RANDOMGENERATOR_H_
    7 
    8 #undef M_PI
    9 #define M_PI 3.141592653589793238462643
    10 
    11 class rnd_uniform_distribution
    12 {
    13 
    14   private:
    15     int a;
    16     int c;
    17     unsigned int m;
    18     unsigned _seed;
    19     double lbound;
    20     double ubound;
    21 
    22   public:
    23 
    24     /*constructors, destructors: */
    25     rnd_uniform_distribution() : _seed( 0 ), a( 1103515245 ), c( 12345 ), m( 2147483648 ), lbound(0.0), ubound(1.0) {}
    26     rnd_uniform_distribution(double a_1,double a_2) : _seed( 0 ), a( 1103515245 ), c( 12345 ), m( 2147483648 ), lbound(a_1), ubound(a_2) {}
    27     ~rnd_uniform_distribution(){}
    28 
    29     void seed( unsigned int s ) { _seed = s; }
    30     unsigned int get_seed() { return _seed; }
    31     double generator(){
    32       _seed = ( a * _seed + c ) % m;
    33       return (ubound-lbound)*(double) _seed/ m + lbound;
    34     }
    35 
    36 };
    37 
    38 class rnd_normal_distribution
    39 {
    40 
    41   private:
    42     unsigned _seed;
    43     double mean;
    44     double sdev;
    45 
    46   public:
    47 
    48     /*constructors, destructors: */
    49     rnd_normal_distribution() : _seed( 0 ), mean( 0), sdev(1.0) {}
    50     rnd_normal_distribution(double m,double s) : _seed( 0 ), mean( m ), sdev(s) {}
    51     ~rnd_normal_distribution(){}
    52 
    53     void seed( unsigned int s ) { _seed = s; }
    54     double generator()
    55     {
    56       rnd_uniform_distribution unifdistri;
    57       unifdistri.seed(_seed);
    58 
    59       double u1 = unifdistri.generator();
    60       double u2 = unifdistri.generator();
    61 
    62       double R = sqrt(-2*log(u1));
    63       double theta = 2*M_PI*u2;
    64 
    65       seed(unifdistri.get_seed());
    66 
    67       return mean + sdev * (R*cos(theta));
    68 
    69     }
    70 
    71 };
    72 
    73 //#endif //ifndef _RANDOMGENERATOR_H_

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

    r26592 r26597  
    1818
    1919        /*Define seed*/
    20         rnd_normal_distribution distribution;
     20        //rnd_normal_distribution distribution;
    2121        if(seed<0){
    2222                std::random_device rd;
    … …  
    2929                seed = seed + 783728*my_rank; // change default seed for parallel simulations (by considering an arbitrary shif based on the rank number)
    3030        }
    31         distribution.seed(seed);
     31        //distribution.seed(seed);
    3232
    3333        int        *local_indices = NULL;
    … …  
    3838        ppf->GetLocalSize(&M);
    3939        for(int i=0;i<M;i++){
    40                 rdnumber = distribution.generator();
     40                rdnumber = 1.0;// distribution.generator();
    4141                ppf->SetValue(local_indices[i],rdnumber,INS_VAL);
    4242        }
Note: See TracChangeset for help on using the changeset viewer.