random_number_generator.h

Go to the documentation of this file.

#ifndef CYCLUS_SRC_RNG_H
#define CYCLUS_SRC_RNG_H
 
#include <boost/random.hpp>
#include <cstdint>
 
#include "error.h"
 
class SimInitTest;
class RandomTest;
 
namespace cyclus {
 
class SimInfo;
 
typedef boost::random::mt19937 Generator;
typedef boost::random::variate_generator<Generator&, boost::random::normal_distribution<> > NormalDist;
 
/// A random number generator.
class RandomNumberGenerator {
    friend class ::SimInitTest;
    friend class ::RandomTest;
 
    // all distributions are fiends
    friend class FixedDoubleDist;
    friend class UniformDoubleDist;
    friend class NormalDoubleDist;
    friend class PoissonDoubleDist;
    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
    static Generator gen_;
 
  public:
    /// Initialize from seed
    void Initialize(SimInfo si);
 
    std::uint32_t random();
 
    // in archetype code
    //
    // #include "random_number_generator.h"
    //
    // class Foo:
    //    private:
    //       NormalDist myNormalDist;
    //
    // Foo::Foo() {
    //    myNormalDist = RandomNumberGenerator::makeNormalDist(mean, std_dev, min, max);
    // }
    //
    // void Foo::Tick() {
    //     bar = myNormalDist();
    // }
    //
    
    static NormalDist makeNormalDist(double mean, double std_dev, double min, double max) {
        boost::random::normal_distribution<> dist(mean, std_dev);
        NormalDist rn(gen_, dist);
 
        return rn;
    }
 
    /// wrappers for boost::random distributions
 
    /// generate a random number between [0, 1)
    double random_01();
    
    /// 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());
 
};
  
class DoubleDistribution {
  public:
    typedef boost::shared_ptr<DoubleDistribution> Ptr;
 
    virtual double sample() = 0;
    virtual double max() = 0;
};
 
class FixedDoubleDist : public DoubleDistribution {
  private:
    double value;
  public:
    typedef boost::shared_ptr<FixedDoubleDist> Ptr;
 
    FixedDoubleDist(double value_) : value(value_) {};
    virtual double sample() { return value; };
    virtual double max() { return value; };
};
 
class UniformDoubleDist : public DoubleDistribution {
  private:
    boost::random::uniform_real_distribution<> dist;
  public:
    typedef boost::shared_ptr<UniformDoubleDist> Ptr;
 
    UniformDoubleDist(double min = 0, double max=1) : dist(min, max) {};
    virtual double sample() { return dist(RandomNumberGenerator::gen_); }
    virtual double max() { return dist.max(); }
};
 
class NormalDoubleDist : public DoubleDistribution {
  private:
    boost::random::normal_distribution<> dist;
    double min_;
    double max_;
  public:
    typedef boost::shared_ptr<NormalDoubleDist> Ptr;
 
    NormalDoubleDist(double mean, double std_dev, double min=0, double max=1) : dist(mean, std_dev), min_(min), max_(max) {
      if (min_ == max_) {
        throw ValueError("Min and max cannot be equal for a normal distribution. Either use FixedDoubleDist or change the min/max.");
      }
      if (max_ < (mean - 3*std_dev) || min_ > (mean + 3*std_dev)) {
        Warn<VALUE_WARNING>("Dist is sampling from a tail of a truncated normal more than 3 standard deviations from the mean. Drawing sampling may be inefficient");
      }
    };
    virtual double sample();
    virtual double max() { return max_; }
};
 
/// Poisson distribution requires a mean
class PoissonDoubleDist : public DoubleDistribution {
  private:
    boost::random::poisson_distribution<> dist;
    double mean_;
  public:
    typedef boost::shared_ptr<PoissonDoubleDist> Ptr;
 
    PoissonDoubleDist(double mean) : dist(mean_) {
      if (mean_ < 0) {
        throw ValueError("Mean must be positive");
      }
    };
    virtual double sample() { return dist(RandomNumberGenerator::gen_); }
    virtual double mean() { return dist.mean(); }
};
 
/// Exponential distribution requires lambda
class ExponentialDoubleDist : public DoubleDistribution {
  private:
    boost::random::exponential_distribution<> dist;
    double lambda_;
  public:
    typedef boost::shared_ptr<ExponentialDoubleDist> Ptr;
 
    ExponentialDoubleDist(double lambda) : dist(lambda), lambda_(lambda) {
      if (lambda_ < 0) {
        throw ValueError("Lambda must be positive");
      }
    };
    virtual double sample() { return dist(RandomNumberGenerator::gen_); }
    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;
    virtual int sample() = 0;
};
 
class FixedIntDist : public IntDistribution {
  private:
    int value;
  public:
    typedef boost::shared_ptr<FixedIntDist> Ptr;
 
    FixedIntDist(int value_) : value(value_) {};
    virtual int sample() { return value; };
};
 
class UniformIntDist : public IntDistribution {
  private:
    boost::random::uniform_int_distribution<> dist;
  public:
    typedef boost::shared_ptr<UniformIntDist> Ptr;
 
    UniformIntDist(int min = 0, int max=1) : dist(min, max) {};
    virtual int sample() { return dist(RandomNumberGenerator::gen_); }
    virtual int max() { return dist.max(); }
};
 
class NormalIntDist : public IntDistribution {
  private:
    boost::random::normal_distribution<> dist;
    int min_;
    int max_;
  public:
    typedef boost::shared_ptr<NormalIntDist> Ptr;
 
    NormalIntDist(double mean, double std_dev, int min=0, int max=1) : dist(mean, std_dev), min_(min), max_(max) {
      if (min_ == max_) {
        throw ValueError("Min and max cannot be equal for a normal distribution. Either use FixedIntDist or change the min/max.");
      }
      if (max_ < (mean - 3*std_dev) || min_ > (mean + 3*std_dev)) {
        Warn<VALUE_WARNING>("Dist is sampling from a tail of a truncated normal more than 3 standard deviations from the mean. Drawing sampling may be inefficient");
      }
    };
    virtual int sample();
    virtual int max() { return max_; }
};
 
/// Binomial distribution requries an integer number of trials and a
/// probability of success for each trial. Sampling returns the number of 
/// successes. When trials is one, this is equivalent to a Bernoulli dist
class BinomialIntDist : public IntDistribution {
  private:
    boost::random::binomial_distribution<int> dist;
    int trials_;
    double p_success_;
  public:
    typedef boost::shared_ptr<BinomialIntDist> Ptr;
 
    BinomialIntDist(int trials, double p_success) : dist(trials, p_success), trials_(trials), p_success_(p_success) {
      if (trials_ < 1) {
        throw ValueError("Max must be positive and greater than zero");
      }
      if (p_success_ > 1 || p_success_ < 0) {
        throw ValueError("Probability must be between zero and one");
      }
    };
    virtual int sample() { return dist(RandomNumberGenerator::gen_); }
    virtual int trials() { return trials_; }
    virtual int p() { return p_success_; }
};
 
/// NegativeBinomialIntDist takes the number of successes desired and a 
/// probability of success on a single trial and returns the number of trials 
/// needed to reach the desired successes. When successes is one, this is
/// equivalent to a Geometric distribution.
class NegativeBinomialIntDist : public IntDistribution {
  private:
    boost::random::negative_binomial_distribution<> dist;
    int successes_;
    double p_success_;
  public:
    typedef boost::shared_ptr<NegativeBinomialIntDist> Ptr;
 
    NegativeBinomialIntDist(int successes, double p_success) : dist(successes, p_success), successes_(successes), p_success_(p_success) {
      if (successes < 1) {
        throw ValueError("Successes must be positive and greater than zero");
      }
      if (p_success > 1 || p_success < 0) {
        throw ValueError("Probability must be between zero and one");
      }
    };
    virtual int sample() { return dist(RandomNumberGenerator::gen_); }
    virtual int successes() { return successes_; }
    virtual int p() { return p_success_; }
};
 
/// Poisson distribution requires a mean 
class PoissonIntDist : public IntDistribution {
  private:
    boost::random::poisson_distribution<> dist;
    double mean_;
  public:
    typedef boost::shared_ptr<PoissonIntDist> Ptr;
 
    PoissonIntDist(double mean) : dist(mean_) {
      if (mean_ < 0) {
        throw ValueError("Mean must be positive");
      }
    };
    virtual int sample() { return dist(RandomNumberGenerator::gen_); }
    virtual double mean() { return dist.mean(); }
};
 
/// Exponential distribution requires lambda
class ExponentialIntDist : public IntDistribution {
  private:
    boost::random::exponential_distribution<> dist;
    double lambda_;
  public:
    typedef boost::shared_ptr<ExponentialIntDist> Ptr;
 
    ExponentialIntDist(double lambda) : dist(lambda), lambda_(lambda) {
      if (lambda_ < 0) {
        throw ValueError("Lambda must be positive");
      }
    };
    virtual int sample() { return dist(RandomNumberGenerator::gen_); }
    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