RNG, add seed in simulation init and initialize PRNG

CHANGELOG.rst

@@ -19,6 +19,7 @@ Since last release
 * build and test are now fown on githubAction in place or CircleCI (#1569)
 * Have separate workflows for testing, publishing dependency images, and publishing release images (#1597, #1602, #1606, #1609)
 * Add Ubuntu 20.04 to the list of supported platforms (#1605, #1608)
+* Add random number generator (Mersenne Twister 19937, from boost) and the ability to set the seed in the simulation control block (#1599)
 
 **Changed:**
 

src/agent.h

@@ -21,6 +21,7 @@
 #endif
 
 class SimInitTest;
+class RandomTest;
 
 namespace cyclus {
 
@@ -49,6 +50,7 @@ typedef std::map<std::string, std::vector<Resource::Ptr> > Inventories;
 class Agent : public StateWrangler, virtual public Ider {
   friend class SimInit;
   friend class ::SimInitTest;
+  friend class ::RandomTest;
 
  public:
   /// Creates a new agent that is managed by the given context. Note that the

src/composition.h

@@ -6,6 +6,7 @@
 #include <boost/shared_ptr.hpp>
 
 class SimInitTest;
+class RandomTest;
 
 namespace cyclus {
 
@@ -38,6 +39,7 @@ typedef std::map<Nuc, double> CompMap;
 class Composition {
   friend class SimInit;
   friend class ::SimInitTest;
+  friend class ::RandomTest;
 
  public:
   typedef boost::shared_ptr<Composition> Ptr;

src/context.cc

@@ -9,6 +9,7 @@
 #include "pyhooks.h"
 #include "sim_init.h"
 #include "timer.h"
+#include "random_number_generator.h"
 #include "version.h"
 
 namespace cyclus {
@@ -26,7 +27,9 @@ SimInfo::SimInfo()
       explicit_inventory(false),
       explicit_inventory_compact(false),
       parent_sim(boost::uuids::nil_uuid()),
-      parent_type("init") {}
+      parent_type("init"),
+      seed(kDefaultSeed),
+      stride(kDefaultStride) {}
 
 SimInfo::SimInfo(int dur, int y0, int m0, std::string handle)
     : duration(dur),
@@ -39,7 +42,9 @@ SimInfo::SimInfo(int dur, int y0, int m0, std::string handle)
       explicit_inventory(false),
       explicit_inventory_compact(false),
       parent_sim(boost::uuids::nil_uuid()),
-      parent_type("init") {}
+      parent_type("init"),
+      seed(kDefaultSeed),
+      stride(kDefaultStride) {}
 
 SimInfo::SimInfo(int dur, int y0, int m0, std::string handle, std::string d)
     : duration(dur),
@@ -52,7 +57,9 @@ SimInfo::SimInfo(int dur, int y0, int m0, std::string handle, std::string d)
       explicit_inventory(false),
       explicit_inventory_compact(false),
       parent_sim(boost::uuids::nil_uuid()),
-      parent_type("init") {}
+      parent_type("init"),
+      seed(kDefaultSeed),
+      stride(kDefaultStride) {}
 
 SimInfo::SimInfo(int dur, boost::uuids::uuid parent_sim,
                  int branch_time, std::string parent_type,
@@ -67,20 +74,27 @@ SimInfo::SimInfo(int dur, boost::uuids::uuid parent_sim,
       branch_time(branch_time),
       explicit_inventory(false),
       explicit_inventory_compact(false),
-      handle(handle) {}
+      handle(handle),
+      seed(kDefaultSeed),
+      stride(kDefaultStride) {}
 
 Context::Context(Timer* ti, Recorder* rec)
     : ti_(ti),
       rec_(rec),
       solver_(NULL),
       trans_id_(0),
-      si_(0) {}
+      si_(0) {
+
+        rng_ = new RandomNumberGenerator();
+      }
 
 Context::~Context() {
   if (solver_ != NULL) {
     delete solver_;
   }
-
+  if (rng_ != NULL) {
+    delete rng_;
+  }
   // initiate deletion of agents that don't have parents.
   // dealloc will propagate through hierarchy as agents delete their children
   std::vector<Agent*> to_del;
@@ -184,6 +198,8 @@ void Context::InitSim(SimInfo si) {
       ->AddVal("InitialYear", si.y0)
       ->AddVal("InitialMonth", si.m0)
       ->AddVal("Duration", si.duration)
+      ->AddVal("Seed", static_cast<int>(si.seed))
+      ->AddVal("Stride", static_cast<int>(si.stride))
       ->AddVal("ParentSimId", si.parent_sim)
       ->AddVal("ParentType", si.parent_type)
       ->AddVal("BranchTime", si.branch_time)
@@ -222,12 +238,40 @@ void Context::InitSim(SimInfo si) {
 
   si_ = si;
   ti_->Initialize(this, si);
+  rng_->Initialize(si);
+
 }
 
 int Context::time() {
   return ti_->time();
 }
 
+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) {
   ti_->RegisterTimeListener(tl);
 }

src/context.h

@@ -21,7 +21,12 @@
 
 const uint64_t kDefaultTimeStepDur = 2629846;
 
+const uint64_t kDefaultSeed = 20160212;
+
+const uint64_t kDefaultStride = 10000;
+
 class SimInitTest;
+class RandomTest;
 
 namespace cyclus {
 
@@ -33,6 +38,7 @@ class Timer;
 class TimeListener;
 class SimInit;
 class DynamicModule;
+class RandomNumberGenerator;
 
 /// Container for a static simulation-global parameters that both describe
 /// the simulation and affect its behavior.
@@ -114,6 +120,15 @@ class SimInfo {
   /// every time step in a table (i.e. agent ID, Time, Quantity,
   /// Composition-object and/or reference).
   bool explicit_inventory_compact;
+
+  /// Seed for random number generator
+  uint64_t seed;
+
+  /// Stride length. Currently unused, but available for future development 
+  /// that may wish to initiate multiple random number generators from the
+  /// same seed, skipping forward in the sequence by the stride length times
+  /// some parameter, such as the agent_id. 
+  uint64_t stride;
 };
 
 /// A simulation context provides access to necessary simulation-global
@@ -130,6 +145,7 @@ class SimInfo {
 class Context {
  public:
   friend class ::SimInitTest;
+  friend class ::RandomTest;
   friend class SimInit;
   friend class Agent;
   friend class Timer;
@@ -234,9 +250,34 @@ class Context {
   /// Returns the current simulation timestep.
   virtual int time();
 
+  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;};
 
+  /// Returns the seed for the random number generator.
+  inline uint64_t seed() {return si_.seed;};
+
+  /// Returns the stride for the random number generator.
+  inline uint64_t stride() {return si_.stride;};
+
   /// Return static simulation info.
   inline SimInfo sim_info() const {
     return si_;
@@ -312,6 +353,7 @@ class Context {
   ExchangeSolver* solver_;
   Recorder* rec_;
   int trans_id_;
+  RandomNumberGenerator* rng_;
 };
 
 }  // namespace cyclus

src/dynamic_module.h

@@ -8,6 +8,7 @@
 
 // for testing
 class SimInitTest;
+class RandomTest;
 
 namespace cyclus {
 
@@ -94,9 +95,10 @@ class DynamicModule {
   /// added to this map when loaded.
   static std::map<std::string, DynamicModule*> modules_;
 
-  /// for testing - see sim_init_tests
+  /// for testing - see sim_init_tests and random_tests
   friend class ::SimInitTest;
-  /// for testing - see sim_init_tests
+  friend class ::RandomTest;
+  /// for testing - see sim_init_tests and random_tests
   static std::map<std::string, AgentCtor*> man_ctors_;
 
   /// the name of the module

src/product.h

@@ -8,6 +8,7 @@
 #include "res_tracker.h"
 
 class SimInitTest;
+class RandomTest;
 
 namespace cyclus {
 
@@ -18,6 +19,7 @@ namespace cyclus {
 class Product : public Resource {
   friend class SimInit;
   friend class ::SimInitTest;
+  friend class ::RandomTest;
 
  public:
   typedef

src/random_number_generator.cc

@@ -0,0 +1,63 @@
+#include "random_number_generator.h"
+
+#include <iostream>
+#include <string>
+
+#include "logger.h"
+#include "sim_init.h"
+#include "context.h"
+
+namespace cyclus{
+
+    Generator RandomNumberGenerator::gen_;
+
+    void RandomNumberGenerator::Initialize(SimInfo si){
+
+        gen_.seed(si.seed);
+
+        CLOG(LEV_INFO1) << "Pseudo random number generator initialized with seed: " << si.seed;
+    }
+
+    std::uint32_t RandomNumberGenerator::random(){
+        return gen_();
+    }
+
+    double RandomNumberGenerator::random_01(){
+        boost::random::uniform_01<> dist;
+        return dist(gen_);
+    }
+
+    int RandomNumberGenerator::random_uniform_int(int low, int high){
+        boost::random::uniform_int_distribution<> dist(low, high);
+        boost::random::variate_generator<Generator&, boost::random::uniform_int_distribution<> > rn(gen_, dist);
+        return rn();
+    }
+
+    double RandomNumberGenerator::random_uniform_real(double low, double high){
+        boost::random::uniform_real_distribution<> dist(low, high);
+        boost::random::variate_generator<Generator&, boost::random::uniform_real_distribution<> > rn(gen_, dist);
+        return rn();
+    }
+
+    double RandomNumberGenerator::random_normal_real(double mean, double std_dev, double low, double high){
+        boost::random::normal_distribution<> dist(mean, std_dev);
+        boost::random::variate_generator<Generator&, boost::random::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::random::normal_distribution<> dist(mean, std_dev);
+        boost::random::variate_generator<Generator&, boost::random::normal_distribution<> > rn(gen_, dist);
+        double val = rn();
+        while (val < low || val > high){
+            val = rn();
+        }
+        int rounded_val = std::lrint(val);
+        return val;
+    }
+
+}