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