1 /*
2  * scrm is an implementation of the Sequential-Coalescent-with-Recombination Model.
3  *
4  * Copyright (C) 2013, 2014 Paul R. Staab, Sha (Joe) Zhu, Dirk Metzler and Gerton Lunter
5  *
6  * This file is part of scrm.
7  *
8  * scrm is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifndef scrm_src_random_random_generator
24 #define scrm_src_random_random_generator
25 
26 #include "../macros.h" // Needs to be before cassert
27 
28 #include <cassert>
29 #include <cmath>
30 #include <memory>
31 
32 #include "fastfunc.h"
33 
34 
35 class RandomGenerator
36 {
37  friend class MersenneTwister;
38  public:
RandomGenerator()39   RandomGenerator() : ff_(std::make_shared<FastFunc>()) { };
RandomGenerator(std::shared_ptr<FastFunc> ff)40   RandomGenerator(std::shared_ptr<FastFunc> ff) : ff_(ff) { };
41 
~RandomGenerator()42   virtual ~RandomGenerator() {}
43 
44   //Getters & Setters
seed()45   size_t seed() const { return seed_; }
46 
set_seed(const size_t seed)47   virtual void set_seed(const size_t seed) {
48     this->seed_ = seed;
49   }
50 
51   virtual double sample() =0;
52 
53   // Base class methods
54   // Initialize unit_exponential; must be called when the random generator is up and running
initializeUnitExponential()55   void initializeUnitExponential() {
56     this->unit_exponential_ = sampleUnitExponential();
57   }
58 
59   // Uniformly samples a number out of 0, ..., range-1
60   // Unit tested
sampleInt(const int max_value)61   int sampleInt(const int max_value) {
62     return(static_cast<int>(this->sample()*max_value));
63   }
64 
65   // Samples from an exponential distribution
66   // Unit tested
sampleExpo(const double lambda)67   double sampleExpo(const double lambda) {
68     return sampleUnitExponential() / lambda;
69   }
70 
71   // Samples from an exponential distribution; return -1 if beyond limit
72   // If a limit is known, this version is faster than the standard one
73   // Unit tested
sampleExpoLimit(const double lambda,const double limit)74   double sampleExpoLimit(const double lambda, const double limit) {
75     return sampleExpoExpoLimit(lambda, 0, limit);
76   }
77 
78   double sampleExpoExpoLimit(const double b, const double c, const double limit);
79 
80 #ifdef UNITTEST
81   friend class TestRandomGenerator;
82 #endif
83 
84   // fast functions
ff()85   std::shared_ptr<FastFunc> ff() { return this->ff_; }
86 
87 
88  protected:
89   // Sample from a unit exponential distribution
90   // Unit tested
91   // fastlog can return 0, which causes a bug in scrm.
92   // log or fastlog does seems to have an influence on the runtime.
sampleUnitExponential()93   virtual double sampleUnitExponential() {
94     return -(this->ff()->fastlog(sample()));
95     //return -std::log(sample());
96   }
97 
98  protected:
99   // seed
100   size_t seed_;
101   // cache for a unit-exponentially distributed variable
102   double unit_exponential_;
103   std::shared_ptr<FastFunc> ff_;
104 };
105 
106 #endif
107