source: issm/oecreview/Archive/25834-26739/ISSM-26577-26578.diff@ 27230

Last change on this file since 27230 was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 7.6 KB

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

    Cannot display: file marked as a binary type.
svn:mime-type = application/octet-stream
     
    88#undef M_PI
    99#define M_PI 3.141592653589793238462643
    1010
     11namespace rdngen{
     12        uniform_distribution_rnd::uniform_distribution_rnd(){/*{{{*/
    1113
    12 uniform_distribution_rnd::uniform_distribution_rnd(){/*{{{*/
     14                        a   = 1103515245;       // BSD Formula
     15                        c  = 12345;                                     // BSD Formula
     16                        m = 2147483648;                 // BSD Formula
     17                        _seed = 0;
     18                        lbound = 0.0;
     19                        ubound = 1.0;
     20                        return;
     21        }
     22                /*}}}*/
     23        uniform_distribution_rnd::uniform_distribution_rnd(double lower,double upper){/*{{{*/
    1324
    14                 a   = 1103515245;       // BSD Formula
    15                 c  = 12345;                                     // BSD Formula
    16                 m = 2147483648;                 // BSD Formula
    17                 _seed = 0;
    18                 lbound = 0.0;
    19                 ubound = 1.0;
    20                 return;
    21 }
    22         /*}}}*/
    23 uniform_distribution_rnd::uniform_distribution_rnd(double lower,double upper){/*{{{*/
     25                        a   = 1103515245;               // BSD Formula
     26                        c  = 12345;                                     // BSD Formula
     27                        m = 2147483648;                 // BSD Formula
     28                        _seed = 0;
     29                        lbound = lower;
     30                        ubound = upper;
     31                        return;
     32        }
     33                /*}}}*/
     34        uniform_distribution_rnd::~uniform_distribution_rnd(){}
     35        void uniform_distribution_rnd::seed( unsigned int s ) { _seed = s; }
     36        unsigned int uniform_distribution_rnd::get_seed() { return _seed; }
     37        double uniform_distribution_rnd::generator() {
     38                        _seed = ( a * _seed + c ) % m ;
     39                        return (ubound-lbound)*(double) _seed/ m + lbound;
     40        }
    2441
    25                 a   = 1103515245;               // BSD Formula
    26                 c  = 12345;                                     // BSD Formula
    27                 m = 2147483648;                 // BSD Formula
    28                 _seed = 0;
    29                 lbound = lower;
    30                 ubound = upper;
    31                 return;
    32 }
    33         /*}}}*/
    34 uniform_distribution_rnd::~uniform_distribution_rnd(){}
    35 void uniform_distribution_rnd::seed( unsigned int s ) { _seed = s; }
    36 unsigned int uniform_distribution_rnd::get_seed() { return _seed; }
    37 double uniform_distribution_rnd::generator() {
    38                 _seed = ( a * _seed + c ) % m ;
    39                 return (ubound-lbound)*(double) _seed/ m + lbound;
    40 }
    4142
     43        normal_distribution_rnd::normal_distribution_rnd(){/*{{{*/
    4244
    43 normal_distribution_rnd::normal_distribution_rnd(){/*{{{*/
     45                        _seed = 0;
     46                        mean   = 0;
     47                        sdev  = 1.0;
     48                        return;
     49        }
     50        /*}}}*/
     51        normal_distribution_rnd::normal_distribution_rnd(double m,double s){/*{{{*/
    4452
    45                 _seed = 0;
    46                 mean   = 0;
    47                 sdev  = 1.0;
    48                 return;
    49 }
    50 /*}}}*/
    51 normal_distribution_rnd::normal_distribution_rnd(double m,double s){/*{{{*/
     53                        _seed = 0;
     54                        mean   = m;
     55                        sdev  = s;
     56                        return;
     57        }
     58                /*}}}*/
     59        normal_distribution_rnd::~normal_distribution_rnd(){}
     60        void normal_distribution_rnd::seed( unsigned int s ) { _seed = s; }
     61        double normal_distribution_rnd::generator(){/*{{{*/
    5262
    53                 _seed = 0;
    54                 mean   = m;
    55                 sdev  = s;
    56                 return;
    57 }
    58         /*}}}*/
    59 normal_distribution_rnd::~normal_distribution_rnd(){}
    60 void normal_distribution_rnd::seed( unsigned int s ) { _seed = s; }
    61 double normal_distribution_rnd::generator(){/*{{{*/
     63                        uniform_distribution_rnd unifdistri;
     64                        unifdistri.seed(_seed);
    6265
    63                 uniform_distribution_rnd unifdistri;
    64                 unifdistri.seed(_seed);
     66                        double u1 = unifdistri.generator();
     67                        double u2 = unifdistri.generator();
    6568
    66                 double u1 = unifdistri.generator();
    67                 double u2 = unifdistri.generator();
     69                        double R = sqrt(-2*log(u1));
     70                        double theta = 2*M_PI*u2;
    6871
    69                 double R = sqrt(-2*log(u1));
    70                 double theta = 2*M_PI*u2;
     72                        seed(unifdistri.get_seed());
    7173
    72                 seed(unifdistri.get_seed());
     74                        return mean + sdev * (R*cos(theta));
    7375
    74                 return mean + sdev * (R*cos(theta));
    75 
    76         }
    77         /*}}}*/
     76                }
     77                /*}}}*/
     78}                               

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

     
    88#undef M_PI
    99#define M_PI 3.141592653589793238462643
    1010
    11 class uniform_distribution_rnd
    12 {
     11namespace rdngen{
     12  class uniform_distribution_rnd
     13  {
    1314
    14     private:
     15      private:
    1516
    16       unsigned int a;       //multiplier of the linear congruential generator
    17       unsigned int c;       //increment of the linear congruential generator
    18       unsigned int m;       // modulo of the linear congruential generator
    19       unsigned int _seed;   // seed value
    20       double lbound;        // lower bound of uniform distribution
    21       double ubound;        // upper bound of uniform distribution
     17        unsigned int a;       //multiplier of the linear congruential generator
     18        unsigned int c;       //increment of the linear congruential generator
     19        unsigned int m;       // modulo of the linear congruential generator
     20        unsigned int _seed;   // seed value
     21        double lbound;        // lower bound of uniform distribution
     22        double ubound;        // upper bound of uniform distribution
    2223
    23     public:
     24      public:
    2425
    25       /*constructors, destructors: */
    26       uniform_distribution_rnd();
    27       uniform_distribution_rnd(double a_1, double a_2);
    28       ~uniform_distribution_rnd();
     26        /*constructors, destructors: */
     27        uniform_distribution_rnd();
     28        uniform_distribution_rnd(double a_1, double a_2);
     29        ~uniform_distribution_rnd();
    2930
    30       void seed( unsigned int s );
    31       unsigned int get_seed();
    32       double generator();
     31        void seed( unsigned int s );
     32        unsigned int get_seed();
     33        double generator();
    3334
    34 };
     35  };
    3536
    36 class normal_distribution_rnd
    37 {
     37  class normal_distribution_rnd
     38  {
    3839
    39     private:
    40       unsigned int _seed; // seed value
    41       double mean;        // mean value
    42       double sdev;        // standard deviation
     40      private:
     41        unsigned int _seed; // seed value
     42        double mean;        // mean value
     43        double sdev;        // standard deviation
    4344
    44     public:
     45      public:
    4546
    46       /*constructors, destructors: */
    47       normal_distribution_rnd();
    48       normal_distribution_rnd(double m,double s);
    49       ~normal_distribution_rnd();
     47        /*constructors, destructors: */
     48        normal_distribution_rnd();
     49        normal_distribution_rnd(double m,double s);
     50        ~normal_distribution_rnd();
    5051
    51       void seed( unsigned int s );
    52       double generator();
     52        void seed( unsigned int s );
     53        double generator();
    5354
    54 };
     55  };
    5556
     57}
    5658#endif //ifndef _RANDOMGENERATOR_H_

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

     
    1717        double      rdnumber;
    1818
    1919        /*Define seed*/
    20         //normal_distribution_rnd distribution;
     20        rdngen::normal_distribution_rnd distribution;
    2121        if(seed<0){
    2222                std::random_device rd;
    2323                seed = rd();
    … …  
    2828                ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&my_rank);
    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;
    3434        IssmDouble *local_vector = NULL;
    … …  
    3737  int M;
    3838        ppf->GetLocalSize(&M);
    3939        for(int i=0;i<M;i++){
    40                 rdnumber = 1.0;//distribution.rnd();
     40                rdnumber = distribution.generator();
    4141                ppf->SetValue(local_indices[i],rdnumber,INS_VAL);
    4242        }
    4343        ppf->Assemble();
Note: See TracBrowser for help on using the repository browser.