Changeset 26583


Ignore:
Timestamp:
11/09/21 16:07:14 (3 years ago)
Author:
bulthuis
Message:

BUG: Try to fix bugs on tests

Location:
issm/trunk-jpl
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r26579 r26583  
    55#include <iostream>
    66#include "./randomgenerator.h"
    7 
    8 #undef M_PI
    9 #define M_PI 3.141592653589793238462643
    10 
    11 
    12         uniform_distribution_rnd::uniform_distribution_rnd(){/*{{{*/
    13 
    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){/*{{{*/
    24 
    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         }
    41 
    42 
    43         normal_distribution_rnd::normal_distribution_rnd(){/*{{{*/
    44 
    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){/*{{{*/
    52 
    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(){/*{{{*/
    62 
    63                         uniform_distribution_rnd unifdistri;
    64                         unifdistri.seed(_seed);
    65 
    66                         double u1 = unifdistri.generator();
    67                         double u2 = unifdistri.generator();
    68 
    69                         double R = sqrt(-2*log(u1));
    70                         double theta = 2*M_PI*u2;
    71 
    72                         seed(unifdistri.get_seed());
    73 
    74                         return mean + sdev * (R*cos(theta));
    75 
    76                 }
    77                 /*}}}*/
    78                        

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

    r26579 r26583  
    1010
    1111class uniform_distribution_rnd
    12   {
     12{
    1313
    14       private:
     14  private:
     15    int a;
     16    int c;
     17    unsigned int m;
     18    unsigned _seed;
     19    double a1;
     20    double a2;
    1521
    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
     22    int drnd() { return( _seed = ( a * _seed + c ) % m ); }
    2223
    23       public:
     24  public:
     25    uniform_distribution_rnd() : _seed( 0 ), a( 1103515245 ), c( 12345 ), m( 2147483648 ), a1(0.0), a2(1.0) {}
     26    uniform_distribution_rnd(double a_1,double a_2) : _seed( 0 ), a( 1103515245 ), c( 12345 ), m( 2147483648 ), a1(a_1), a2(a_2) {}
     27    void seed( unsigned int s ) { _seed = s; }
     28    unsigned int get_seed() { return _seed; }
     29    double rnd() { return (a2-a1)*(double) drnd()/ m + a1; }
    2430
    25         /*constructors, destructors: */
    26         uniform_distribution_rnd();
    27         uniform_distribution_rnd(double a_1, double a_2);
    28         ~uniform_distribution_rnd();
     31};
    2932
    30         void seed( unsigned int s );
    31         unsigned int get_seed();
    32         double generator();
     33class normal_distribution_rnd
     34{
    3335
    34   };
     36  private:
     37    unsigned _seed;
     38    double mean;
     39    double sdev;
    3540
    36   class normal_distribution_rnd
    37   {
     41  public:
     42    normal_distribution_rnd() : _seed( 0 ), mean( 0), sdev(1.0) {}
     43    normal_distribution_rnd(double m,double s) : _seed( 0 ), mean( m ), sdev(s) {}
     44    void seed( unsigned int s ) { _seed = s; }
     45    double rnd()
     46    {
     47      uniform_distribution_rnd unifdistri;
     48      unifdistri.seed(_seed);
    3849
    39       private:
    40         unsigned int _seed; // seed value
    41         double mean;        // mean value
    42         double sdev;        // standard deviation
     50      double u1 = unifdistri.rnd();
     51      double u2 = unifdistri.rnd();
    4352
    44       public:
     53      double R = sqrt(-2*log(u1));
     54      double theta = 2*M_PI*u2;
    4555
    46         /*constructors, destructors: */
    47         normal_distribution_rnd();
    48         normal_distribution_rnd(double m,double s);
    49         ~normal_distribution_rnd();
     56      seed(unifdistri.get_seed());
    5057
    51         void seed( unsigned int s );
    52         double generator();
     58      return mean + sdev * (R*cos(theta));
    5359
    54   };
     60    }
     61
     62};
    5563
    5664#endif //ifndef _RANDOMGENERATOR_H_

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

    r26582 r26583  
    1818
    1919        /*Define seed*/
    20         //normal_distribution_rnd distribution;
     20        normal_distribution_rnd 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 = 1.0;//distribution.generator();
     40                rdnumber = distribution.rnd();
    4141                ppf->SetValue(local_indices[i],rdnumber,INS_VAL);
    4242        }
Note: See TracChangeset for help on using the changeset viewer.