Changeset 26574


Ignore:
Timestamp:
11/09/21 10:40:43 (3 years ago)
Author:
bulthuis
Message:

CHG: write pseudrandom number generators as a namespace

Location:
issm/trunk-jpl/src/c
Files:
3 edited

Legend:

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

    r26554 r26574  
    55#include <iostream>
    66#include "./randomgenerator.h"
     7
     8#undef M_PI
     9#define M_PI 3.141592653589793238462643
     10
     11namespace rdn
     12{
     13        uniform_distribution::uniform_distribution(){/*{{{*/
     14
     15                a   = 1103515245;       // BSD Formula
     16                c  = 12345;                                     // BSD Formula
     17                m = 2147483648;                 // BSD Formula
     18                _seed = 0;
     19                lbound = 0.0;
     20                ubound = 1.0;
     21                return;
     22        }
     23        /*}}}*/
     24        uniform_distribution::uniform_distribution(double lower,double upper){/*{{{*/
     25
     26                a   = 1103515245;               // BSD Formula
     27                c  = 12345;                                     // BSD Formula
     28                m = 2147483648;                 // BSD Formula
     29                _seed = 0;
     30                lbound = lower;
     31                ubound = upper;
     32                return;
     33        }
     34        /*}}}*/
     35        uniform_distribution::~uniform_distribution(){}
     36        void uniform_distribution::seed( unsigned int s ) { _seed = s; }
     37        unsigned int uniform_distribution::get_seed() { return _seed; }
     38        double uniform_distribution::generator() {
     39                _seed = ( a * _seed + c ) % m ;
     40                return (ubound-lbound)*(double) _seed/ m + lbound;
     41        }
     42
     43
     44        normal_distribution::normal_distribution(){/*{{{*/
     45
     46                _seed = 0;
     47                mean   = 0;
     48                sdev  = 1.0;
     49                return;
     50        }
     51        /*}}}*/
     52        normal_distribution::normal_distribution(double m,double s){/*{{{*/
     53
     54                _seed = 0;
     55                mean   = m;
     56                sdev  = s;
     57                return;
     58        }
     59        /*}}}*/
     60        normal_distribution::~normal_distribution(){}
     61        void normal_distribution::seed( unsigned int s ) { _seed = s; }
     62        double normal_distribution::generator(){/*{{{*/
     63
     64                uniform_distribution unifdistri;
     65                unifdistri.seed(_seed);
     66
     67                double u1 = unifdistri.generator();
     68                double u2 = unifdistri.generator();
     69
     70                double R = sqrt(-2*log(u1));
     71                double theta = 2*M_PI*u2;
     72
     73                seed(unifdistri.get_seed());
     74
     75                return mean + sdev * (R*cos(theta));
     76
     77        }
     78        /*}}}*/
     79}
  • issm/trunk-jpl/src/c/shared/Random/randomgenerator.h

    r26554 r26574  
    66#define _RANDOMGENERATOR_H_
    77
    8 #undef M_PI
    9 #define M_PI 3.141592653589793238462643
     8namespace rdn
     9{
     10  class uniform_distribution
     11  {
    1012
    11 class uniform_distribution_rnd
    12 {
     13    private:
    1314
    14   private:
    15     int a;
    16     int c;
    17     unsigned int m;
    18     unsigned _seed;
    19     double a1;
    20     double a2;
     15      unsigned int a;       //multiplier of the linear congruential generator
     16      unsigned int c;       //increment of the linear congruential generator
     17      unsigned int m;       // modulo of the linear congruential generator
     18      unsigned int _seed;   // seed value
     19      double lbound;        // lower bound of uniform distribution
     20      double ubound;        // upper bound of uniform distribution
    2121
    22     int rnd_int() { return( _seed = ( a * _seed + c ) % m ); }
     22    public:
    2323
    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) rnd_int()/ m + a1; }
     24      /*constructors, destructors: */
     25      uniform_distribution();
     26      uniform_distribution(double a_1, double a_2);
     27      ~uniform_distribution();
    3028
    31 };
     29      void seed( unsigned int s );
     30      unsigned int get_seed();
     31      double generator();
    3232
    33 class normal_distribution_rnd
    34 {
     33  };
    3534
    36   private:
    37     unsigned _seed;
    38     double mean;
    39     double sdev;
     35  class normal_distribution
     36  {
    4037
    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);
     38    private:
     39      unsigned int _seed; // seed value
     40      double mean;        // mean value
     41      double sdev;        // standard deviation
    4942
    50       double u1 = unifdistri.rnd();
    51       double u2 = unifdistri.rnd();
     43    public:
    5244
    53       double R = sqrt(-2*log(u1));
    54       double theta = 2*M_PI*u2;
     45      /*constructors, destructors: */
     46      normal_distribution();
     47      normal_distribution(double m,double s);
     48      ~normal_distribution();
    5549
    56       seed(unifdistri.get_seed());
     50      void seed( unsigned int s );
     51      double generator();
    5752
    58       return mean + sdev * (R*cos(theta));
     53  };
    5954
    60     }
    6155
    62 };
    63 
     56}
    6457#endif //ifndef _RANDOMGENERATOR_H_
  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp

    r26554 r26574  
    1818
    1919        /*Define seed*/
    20         normal_distribution_rnd distribution;
     20        rdn::normal_distribution distribution;
    2121        if(seed<0){
    2222                std::random_device rd;
     
    3838        ppf->GetLocalSize(&M);
    3939        for(int i=0;i<M;i++){
    40                 rdnumber = distribution.rnd();
     40                rdnumber = distribution.generator();
    4141                ppf->SetValue(local_indices[i],rdnumber,INS_VAL);
    4242        }
Note: See TracChangeset for help on using the changeset viewer.