source: trunk/lib/random/random.h @ 389

// $Id: random.h 389 2005-08-15 11:37:07Z peter $

#ifndef _theplu_random_
#define _theplu_random_

#include <c++_tools/statistics/Histogram.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include <string>

namespace theplu {
namespace random {  

  ///
  /// @brief Random Number Generator
  ///
  /// The RNG class is wrapper to the GSL random number generator
  /// (rng). This class provides a single global instance of the rng,
  /// and makes sure there is only one point of access to the
  /// generator.
  ///
  /// There is information about how to change seeding and generators
  /// at run time without recompilation using environment
  /// variables. RNG of course support seeding at compile time if you
  /// don't want to bother about environment variables and GSL.
  ///
  /// There are many different rng's available in GSL. Currently only
  /// the default generator is implemented and no other one is
  /// choosable through the class interface. This means that you have
  /// to fall back to the use of environment variables as described in
  /// the GSL documentation, or be bold and request support for other
  /// rng's through the class interface.
  ///
  /// Not all GSL functionality is implemented, we'll add
  /// functionality when needed and may do it when requested. Better
  /// yet, supply us with code and we will probably add it to the code
  /// (BUT remember to implement reasonable tests for your code and
  /// follow the coding style.)
  ///
  /// This implementation may be thread safe (according to GSL
  /// documentation), but should be checked to be so before trusting
  /// thread safety.
  ///
  /// @see Design Patterns (the singleton and adapter pattern). GSL
  /// documentation.
  ///
  class RNG
  {
  public:

    virtual ~RNG(void);

    static RNG* instance(u_long seed=0);

    ///
    /// @brief Returns the largest number that the random number
    /// generator can return.
    ///
    inline u_long max(void) const { return gsl_rng_max(rng_); }

    ///
    /// @brief Returns the smallest number that the random number
    /// generator can return.
    ///
    inline u_long min(void) const { return gsl_rng_min(rng_); }

    ///
    /// @brief Returns the name of the random number generator
    ///
    inline std::string name(void) const { return gsl_rng_name(rng_); }

    inline const gsl_rng* rng(void) const { return rng_; }

    ///
    /// Set the seed \a s for the rng. If \a s is zero, a default
    /// value from the rng's original implementation is used (cf. GSL
    /// documentation).
    ///
    /// @brief Set the seed \a s for the rng.
    ///
    /// @see seed_dev_urandom
    ///
    inline void seed(u_long s) const { gsl_rng_set(rng_,s); }

    ///
    /// @brief Seed the rng using the /dev/urandom device.
    ///
    /// @return The seed acquired from /dev/urandom.
    ///
    u_long seed_from_devurandom(void);

  private:

    RNG(u_long seed);

    static RNG* instance_;
    gsl_rng* rng_;
  };



  ///
  /// @brief Continuous random number distributions.
  ///
  /// Abstract base class for continuous random number distributions.
  ///
  class Continuous
  {
  public:

    ///
    /// @brief Constructor
    ///
    inline Continuous(void) { rng_=RNG::instance(); }

    ///
    /// @return A random number
    ///
    virtual double operator()(void) const = 0;

  protected:
    RNG* rng_;
  };


  ///
  /// @brief Uniform distribution
  ///
  /// Class for generating a random number from a uniform distribution
  /// in the range [0,1), i.e. zero is included but not 1.
  ///
  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
  /// Expectation value: 0.5 \n
  /// Variance: \f$ \frac{1}{12} \f$
  ///
  class ContinuousUniform : public Continuous
  {
  public:

    inline double operator()(void) const { return gsl_rng_uniform(rng_->rng());}

  };

  ///
  /// @brief Gaussian distribution
  ///
  /// Class for generating a random number from a Gaussian
  /// distribution between zero and unity. Utilizes the Box-Muller
  /// algorithm, which needs two calls to random generator.
  ///
  /// Distribution function \f$ f(x) =
  /// \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
  /// \f$ \n
  /// Expectation value: \f$ \mu \f$ \n
  /// Variance: \f$ \sigma^2 \f$
  ///
  class Gaussian : public Continuous
  {
  public:
    ///
    /// @brief Constructor
    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
    /// m is the expectation value \f$ \mu \f$ of the distribution
    ///
    inline Gaussian(const double s=1, const double m=0) 
      : m_(m), s_(s) {}

    ///
    /// @return A random Gaussian number
    ///
    inline double operator()(void) const 
    { return gsl_ran_gaussian(rng_->rng(), s_)+m_; }

    ///
    /// @return A random Gaussian number with standard deviation \a s
    /// and expectation value 0.
    /// @note this operator ignores parameters given in Constructor
    ///
    inline double operator()(const double s) const 
    { return gsl_ran_gaussian(rng_->rng(), s); }

    ///
    /// @return A random Gaussian number with standard deviation \a s
    /// and expectation value \a m.
    /// @note this operator ignores parameters given in Constructor
    ///
    inline double operator()(const double s, const double m) const 
    { return gsl_ran_gaussian(rng_->rng(), s)+m; }

  private:
    double m_;
    double s_;
  };

  ///
  /// @brief Exponential distribution
  ///
  /// Class for generating a random number from an Exponential
  /// distribution. 
  ///
  /// Distribution function \f$ f(x) = \frac{1}{m}\exp(-x/a) \f$ for
  /// \f$ x \f$ \n 
  /// Expectation value: \f$ m \f$ \n 
  /// Variance: \f$ m^2 \f$
  ///
  class Exponential : public Continuous
  {
  public:
    ///
    /// @brief Constructor
    /// @param m is the expectation value of the distribution.
    ///
    inline Exponential(const double m=1) 
      : m_(m){}

    ///
    /// @return A random number from exponential distribution
    ///
    inline double operator()(void) const 
    { return gsl_ran_exponential(rng_->rng(), m_); }

    ///
    /// @return A random number from exponential distribution, with
    /// expectation value \a m
    /// @note this operator ignores parameters given in Constructor
    /// 
    inline double operator()(const double m) const 
    { return gsl_ran_exponential(rng_->rng(), m); }

  private:
    double m_;
  };

  ///
  /// @brief Discrete random number distributions.
  ///
  /// Abstract base class for discrete random number
  /// distributions. Given K discrete events with different
  /// probabilities \f$ P[k] \f$, produce a random value k consistent
  /// with its probability.
  ///
  class Discrete 
  {
  public:
    ///
    /// @brief Constructor
    ///
    inline Discrete(void) { rng_=RNG::instance(); }

    ///
    /// @return A random number.
    ///
    virtual u_long operator()(void) const = 0;
    
  protected:
    RNG* rng_;
  };

  ///
  /// @brief Discrete uniform distribution
  ///
  /// Discrete uniform distribution also known as the "equally likely
  /// outcomes" distribution. Each outcome, in this case an integer
  /// from [0,n-1] , have equal probability to occur.
  ///
  /// Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
  /// k < n \f$ \n 
  /// Expectation value: \f$ \frac{n-1}{2} \f$ \n 
  /// Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
  ///

  class DiscreteUniform : public Discrete 
  {
  public:
    ///
    /// @brief Constructor.
    ///
    DiscreteUniform(void);

    ///
    /// @brief Constructor.
    ///
    /// @param n sets the range. Object will generate integers from
    /// [0,n-1].
    ///
    /// @todo exception should be thrown, when n is out of range
    ///
    DiscreteUniform(const u_long n);
      
    ///
    /// This function returns a random integer from 0 to n-1
    /// inclusive. All integers in the range [0,n-1] are equally
    /// likely. n is set in constructor.
    ///
    inline u_long operator()(void) const 
    { return gsl_rng_uniform_int(rng_->rng(), range_); }

    ///
    /// This function returns a random integer from 0 to n-1
    /// inclusive. All integers in the range [0,n-1] are equally
    /// likely.
    ///
    /// @todo exception should be thrown, when n is out of range
    ///
    inline u_long operator()(const u_long n) const 
    { return gsl_rng_uniform_int(rng_->rng(), n); }


  private:
    u_long range_;
  };


  ///
  /// @brief Poisson Distribution
  ///
  /// Having a Poisson process (no memory), number of occurences
  /// within a given time window is Poisson distributed. This
  /// distribution is the limit of a Binomial distribution when number
  /// of attempts is large, and the probability for one attempt to be
  /// succesful is small (in such a way that the expected number of
  /// succesful attempts is \f$ m \f$.
  ///
  /// Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
  /// k  \f$ \n 
  /// Expectation value: \f$ m \f$ \n 
  /// Variance: \f$ m \f$
  ///
  class Poisson : public Discrete
  {
  public:
    ///
    /// @brief Constructor
    ///
    /// @param m is expectation value
    inline Poisson(const double m=1) 
      : m_(m){}

    ///
    /// @return A Poisson distributed number.
    ///
    inline u_long operator()(void) const 
    { return gsl_ran_poisson(rng_->rng(), m_); }

    ///
    /// @return A Poisson distributed number with expectation value \a
    /// m
    /// @note this operator ignores parameters set in Constructor
    inline u_long operator()(const double m) const 
    { return gsl_ran_poisson(rng_->rng(), m); }

  private:
    double m_;
  };

  ///
  /// @brief General
  ///
  class DiscreteGeneral : public Discrete 
  {
  public:
    ///
    /// @brief Constructor
    ///
    /// @param hist is a Histogram defining the probability distribution
    ///
    DiscreteGeneral(const statistics::Histogram& hist) ;
    
    ///
    /// @brief Destructor
    ///
    virtual ~DiscreteGeneral(void);

    ///
    /// The generated number is an integer and proportinal to the
    /// frequency in the corresponding histogram bin. In other words,
    /// the probability that 0 is returned is proportinal to the size
    /// of the first bin.
    ///
    /// @return A random number.
    ///
    inline u_long operator()(void) const 
    { return gsl_ran_discrete(rng_->rng(), gen_); }

  private:
     gsl_ran_discrete_t* gen_;
  };


  ///
  /// Class to generate numbers from a histogram in a continuous manner.
  ///
  class ContinuousGeneral : public Continuous
  {
  public:
    ///
    /// @brief Constructor
    ///
    /// @param hist is a Histogram defining the probability distribution
    ///
    inline ContinuousGeneral(const statistics::Histogram& hist) 
      : discrete_(DiscreteGeneral(hist)), hist_(hist) {}
    
    ///
    /// @brief Destructor
    ///
    virtual ~ContinuousGeneral(void);

    ///
    /// The number is generated in a two step process. First the bin
    /// in the histogram is randomly selected (see
    /// DiscreteGeneral). Then a number is generated uniformly from
    /// the interval defined by the bin.
    ///
    /// @return A random number.
    ///
    inline double operator()(void) const 
    { return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing(); }

  private:
    const DiscreteGeneral discrete_;
    const statistics::Histogram& hist_;
    ContinuousUniform u_;
  };


}} // of namespace random and namespace theplu

#endif