add distributions (not yet working)

src/context.cc

@@ -232,16 +232,38 @@ void Context::InitSim(SimInfo si) {
 
   si_ = si;
   ti_->Initialize(this, si);
-  rng_->Initialize(this,si);
+  rng_->Initialize(si);
 
 }
 
 int Context::time() {
   return ti_->time();
 }
 
-int Context::randomnumber() {
-  return rng_->randomnumber();
+int Context::random() {
+  return rng_->random();
+}
+
+double Context::random_01() {
+  return rng_->random_01();
+}
+
+int Context::random_uniform_int(int low, int high) {
+  return rng_->random_uniform_int(low, high);
+}
+
+double Context::random_uniform_real(double low, double high) {
+  return rng_->random_uniform_real(low, high);
+}
+
+double Context::random_normal_real(double mean, double std_dev, double low,
+                                     double high) {
+  return rng_->random_normal_real(mean, std_dev, low, high);
+}
+
+int Context::random_normal_int(double mean, double std_dev, int low,
+                                 int high) {
+  return rng_->random_normal_int(mean, std_dev, low, high);
 }
 
 void Context::RegisterTimeListener(TimeListener* tl) {

src/context.h

@@ -249,8 +249,24 @@ class Context {
   /// Returns the current simulation timestep.
   virtual int time();
 
-  /// Returns the next random number.
-  virtual int randomnumber();
+  int random();
+
+  /// Generates a random number on the range [0,1)]
+  double random_01();
+
+  /// Returns a random number from a uniform integer distribution.
+  int random_uniform_int(int low, int high);
+
+  /// Returns a random number from a uniform real distribution.
+  double random_uniform_real(double low, double high);
+
+  /// Returns a random number from a normal distribution.
+  double random_normal_real(double mean, double std_dev, double low=0,
+                            double high=std::numeric_limits<double>::max());
+
+  /// Returns a random number from a lognormal distribution.
+  int random_normal_int(double mean, double std_dev, int low=0,
+                        int high=std::numeric_limits<int>::max());
 
   /// Returns the duration of a single time step in seconds.
   inline uint64_t dt() {return si_.dt;};

src/random_number_generator.cc

@@ -2,24 +2,60 @@
 
 #include <iostream>
 #include <string>
-#include <boost/random/mersenne_twister.hpp>
 
 #include "logger.h"
 #include "sim_init.h"
 #include "context.h"
 
 namespace cyclus{
 
-    void RandomNumberGenerator::Initialize(Context* ctx, SimInfo si){
+    void RandomNumberGenerator::Initialize(SimInfo si){
 
-        boost::mt19937 gen_(si.seed);
+        Generator gen_(si.seed);
 
         CLOG(LEV_INFO1) << "Pseudo random number generator initialized with seed: " << si.seed;
     }
 
-    int RandomNumberGenerator::randomnumber(){
+    std::uint32_t RandomNumberGenerator::random(){
         return gen_();
+    }
 
-}
+    double RandomNumberGenerator::random_01(){
+        boost::uniform_01<> dist;
+        return dist(gen_);
+    }
+
+    int RandomNumberGenerator::random_uniform_int(int low, int high){
+        boost::uniform_int<> dist(low, high);
+        boost::variate_generator<boost::mt19937&, boost::uniform_int<> > rn(gen_, dist);
+        return rn();
+    }
+
+    double RandomNumberGenerator::random_uniform_real(double low, double high){
+        boost::uniform_real<> dist(low, high);
+        boost::variate_generator<boost::mt19937&, boost::uniform_real<> > rn(gen_, dist);
+        return rn();
+    }
+
+    double RandomNumberGenerator::random_normal_real(double mean, double std_dev, double low, double high){
+        boost::normal_distribution<> dist(mean, std_dev);
+        boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > rn(gen_, dist);
+        double val = rn();
+        while (val < low || val > high){
+            val = rn();
+        }
+        return val;
+    }
+
+    int RandomNumberGenerator::random_normal_int(double mean, double std_dev, int low, int high){
+        boost::normal_distribution<> dist(mean, std_dev);
+        boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > rn(gen_, dist);
+        double val = rn();
+        while (val < low || val > high){
+            val = rn();
+        }
+        int rounded_val = std::lrint(val);
+        return val;
+    }
 
 }

src/random_number_generator.h

@@ -10,21 +10,44 @@ namespace cyclus {
 class Context;
 class SimInfo;
 
+typedef boost::mt19937 Generator;
+
 /// A random number generator.
 class RandomNumberGenerator {
     friend class ::SimInitTest;
+  private:
+    /// Returns a random number for use in a distribution
+    Generator& gen_;
   public:
     RandomNumberGenerator();
 
     ~RandomNumberGenerator();
 
     /// Initialize from seed
-    void Initialize(Context* ctx, SimInfo si);
+    void Initialize(SimInfo si);
+
+    std::uint32_t random();
+
+    /// wrappers for boost::random distributions
 
-    boost::mt19937 gen_;
+    /// geenerate a random number between [0, 1)
+    double random_01();
 
-    /// Returns a random number for use in a distribution
-    int randomnumber();
+    /// generate a random integer between [low, high)
+    int random_uniform_int(int low, int high);
+
+    /// generate a random real number between [low, high)
+    double random_uniform_real(double low, double high);
+
+    /// generate a double from a normal distribution, with truncation
+    /// at low and high
+    double random_normal_real(double mean, double std_dev, double low=0,
+                              double high=std::numeric_limits<double>::max());
+
+    /// generates an integer from a normal distribution, with truncation
+    /// uses rounding to convert double to int
+    int random_normal_int(double mean, double std_dev, int low=0,
+                          int high=std::numeric_limits<int>::max());
 
 };
 

tests/sim_init_tests.cc

@@ -440,3 +440,9 @@ TEST_F(SimInitTest, RestartSimInfo) {
   EXPECT_EQ("restart", info.parent_type);
   EXPECT_EQ(2, info.branch_time);
 }
+
+TEST_F(SimInitTest, GetRandom01) {
+  std::uint32_t r = ctx->random_01();
+  EXPECT_GE(r, 0);
+  EXPECT_LE(r, 1);
+}