Merge pull request cyclus#1827 from gonuke/binary_dists

Add binary distributions for doubles and integers

CHANGELOG.rst

@@ -16,7 +16,7 @@ Since last release
 * New packaging strategies uniform and normal (#1813)
 * Add CI support for MacOS 13, and 14 (#1814, #1819)
 * Methods to reset behavaiors of MaterialBuyPolicy (#1822)
-* More random distributions for the random number generator (#1821)
+* More random distributions for the random number generator (#1821, #1827)
 
 **Changed:**
 

src/random_number_generator.h

@@ -26,14 +26,16 @@ class RandomNumberGenerator {
     friend class UniformDoubleDist;
     friend class NormalDoubleDist;
     friend class PoissonDoubleDist;
-    friend class ExpontentialDoubleDist;
+    friend class ExponentialDoubleDist;
+    friend class BinaryDoubleDist;
     friend class FixedIntDist;
     friend class UniformIntDist;
     friend class NormalIntDist;
     friend class BinomialIntDist;
     friend class NegativeBinomialIntDist;
     friend class PoissonIntDist;
     friend class ExponentialIntDist;
+    friend class BinaryIntDist;
 
   private:
     /// Returns a random number for use in a distribution
@@ -160,14 +162,14 @@ class PoissonDoubleDist : public DoubleDistribution {
 };
 
 /// Exponential distribution requires lambda
-class ExpontentialDoubleDist : public DoubleDistribution {
+class ExponentialDoubleDist : public DoubleDistribution {
   private:
     boost::random::exponential_distribution<> dist;
     double lambda_;
   public:
-    typedef boost::shared_ptr<ExpontentialDoubleDist> Ptr;
+    typedef boost::shared_ptr<ExponentialDoubleDist> Ptr;
 
-    ExpontentialDoubleDist(double lambda) : dist(lambda), lambda_(lambda) {
+    ExponentialDoubleDist(double lambda) : dist(lambda), lambda_(lambda) {
       if (lambda_ < 0) {
         throw ValueError("Lambda must be positive");
       }
@@ -176,6 +178,26 @@ class ExpontentialDoubleDist : public DoubleDistribution {
     virtual double lambda() { return lambda_; }
 };
 
+
+/// Binary distribution requires twp options and a probability
+class BinaryDoubleDist : public DoubleDistribution {
+  private:
+    boost::random::binomial_distribution<int> dist;
+    double p_success_;
+    double choice_a_, choice_b_;
+  public:
+
+    BinaryDoubleDist(double p_success, double choice_a, double choice_b) : 
+      dist(1, p_success),choice_a_(choice_a), choice_b_(choice_b) {
+      if (p_success < 0) {
+        throw ValueError("Probability of choice A must be positive");
+      }
+    };
+    virtual double sample() { 
+      return dist(RandomNumberGenerator::gen_) ? choice_a_ : choice_b_; }
+    virtual double p() { return p_success_; }
+};
+
 class IntDistribution {
   public:
     typedef boost::shared_ptr<IntDistribution> Ptr;
@@ -306,6 +328,26 @@ class ExponentialIntDist : public IntDistribution {
     virtual double lambda() { return lambda_; }
 };
 
+/// Binary distribution requires twp options and a probability
+class BinaryIntDist : public IntDistribution {
+  private:
+    boost::random::binomial_distribution<int> dist;
+    double p_success_;
+    int choice_a_, choice_b_;
+  public:
+
+    BinaryIntDist(double p_success, int choice_a, int choice_b) : 
+      dist(1, p_success),choice_a_(choice_a), choice_b_(choice_b) {
+      if (p_success < 0) {
+        throw ValueError("Probability of choice A must be positive");
+      }
+    };
+    virtual int sample() { 
+      return dist(RandomNumberGenerator::gen_) ? choice_a_ : choice_b_; }
+    virtual double p() { return p_success_; }
+};
+
+
 }
 
 #endif // CYCLUS_SRC_RNG_H

tests/toolkit/matl_buy_policy_tests.cc

@@ -553,6 +553,50 @@ TEST_F(MatlBuyPolicyTests, BinomialActiveDormant) {
   delete a;
 }
 
+TEST_F(MatlBuyPolicyTests, BinaryActiveDormant) {
+  using cyclus::QueryResult;
+
+  boost::shared_ptr<BinaryIntDist> a_dist = boost::shared_ptr<BinaryIntDist>(new BinaryIntDist(0.5, 4, 1));
+  boost::shared_ptr<BinaryIntDist> d_dist = boost::shared_ptr<BinaryIntDist>(new BinaryIntDist(0.5, 1, 6));
+
+  int dur = 100;
+  double throughput = 1;
+
+  cyclus::MockSim sim(dur);
+  cyclus::Agent* a = new TestFacility(sim.context());
+  sim.context()->AddPrototype(a->prototype(), a);
+  sim.agent = sim.context()->CreateAgent<cyclus::Agent>(a->prototype());
+  sim.AddSource("commod1").Finalize();
+
+  TestFacility* fac = dynamic_cast<TestFacility*>(sim.agent);
+
+  cyclus::toolkit::ResBuf<cyclus::Material> inbuf;
+  TotalInvTracker buf_tracker({&inbuf});
+  cyclus::toolkit::MatlBuyPolicy policy;
+  policy.Init(fac, &inbuf, "inbuf", &buf_tracker, throughput, a_dist, d_dist, NULL)
+        .Set("commod1").Start();
+
+  EXPECT_NO_THROW(sim.Run());
+
+  QueryResult qr = sim.db().Query("Transactions", NULL);
+  // sampled cycles: 
+  // * active for 1 step, dormant for 6
+  // * active for 1 step, dormant for 1
+  // * active for 4 steps, dormant for 1 step
+  // therefore first 7 trades occur on time steps:
+  //  0, 7, 9, 10, 11, 12, 14
+  // first trade at time 0, next trade at time 7
+  EXPECT_EQ(0, qr.GetVal<int>("Time", 0));
+  EXPECT_EQ(7, qr.GetVal<int>("Time", 1));
+  EXPECT_EQ(9, qr.GetVal<int>("Time", 2));
+  EXPECT_EQ(10, qr.GetVal<int>("Time", 3));
+  EXPECT_EQ(11, qr.GetVal<int>("Time", 4));
+  EXPECT_EQ(12, qr.GetVal<int>("Time", 5));
+  EXPECT_EQ(14, qr.GetVal<int>("Time", 6));
+
+  delete a;
+}
+
 TEST_F(MatlBuyPolicyTests, MixedActiveDormant) {
   using cyclus::QueryResult;