source: trunk/yat/random/random.h @ 2804

#ifndef _theplu_yat_random_
#define _theplu_yat_random_

// $Id: random.h 2804 2012-07-31 10:06:52Z peter $

/*
  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009, 2010, 2011, 2012 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "yat/statistics/Histogram.h"
#include "yat/utility/deprecate.h"

#include <boost/concept_check.hpp>
#include <boost/thread/tss.hpp>

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

#include <algorithm>
#include <string>
#include <vector>

namespace theplu {
namespace yat {
namespace random {

  //forward declaration
  class RNG_state;

  ///
  /// @brief Random Number Generator
  ///
  /// The RNG class is wrapper to the GSL random number generator
  /// (rng). In yat 0.8 (or older) this class provided a single global
  /// instance of the rng, and made sure there was only one point of
  /// access to the generator. Since version 0.9 this class provides
  /// one rng per thread in order to avoid collisions in multi-thread
  /// applications.
  ///
  /// There are many different rng's available in GSL. RNG uses the
  /// default generator, unless the global variable \c gsl_rng_default
  /// has been modified (see <a
  /// href=\gsl_url/Random-number-environment-variables.html>GSL
  /// Manual</a>). Note, \c gsl_rng_default should be changed before
  /// RNG creates its generator and safest way to achieve this is to
  /// modify \c gsl_rng_default prior calling instance() the first
  /// time.
  ///
  /// There is information about how to change seeding and generators
  /// at run time without recompilation using environment variables in
  /// the <a
  /// href=\gsl_url/Random-number-environment-variables.html>GSL
  /// Manual</a>. RNG supports seeding at compile time if you don't
  /// want to bother about environment variables and GSL.
  ///
  /// The class provides one generator per thread. The first generator
  /// created is seeded with \c gsl_rng_default_seed and subsequent
  /// generators are seeded with \c gsl_rng_default_seed + 1, \c
  /// gsl_rng_default_seed + 2 etc, unless the seed has been modified
  /// with seed() or seed_from_devurandom().
  ///
  /// @see Design Patterns (the singleton and adapter pattern). <a
  /// href=\gsl_url/Random-Number-Generation.html>GSL documentation</a>.
  ///
  class RNG
  {
  public:

    ///
    /// @brief Get an instance of the Random Number Generator.
    ///
    /// Get an instance of RNG. If a random number generator is not
    /// already created for current thread, the call will create a new
    /// generator of type \c gsl_rng_default. If it is the first
    /// generator created it will be seeded with \c
    /// gsl_rng_default_seed; otherwise created generator will be
    /// seeded with \c seed + \c n, where \c seed is the latest seed
    /// set (with seed() or seed_from_devurandom()) The seed may be
    /// changed with the seed or seed_from_devurandom member
    /// functions.
    ///
    /// @return A pointer to the random number generator.
    ///
    /// @see seed and seed_from_devurandom
    ///
    static RNG* instance(void);

    ///
    /// @brief Returns the largest number that the random number
    /// generator can return.
    ///
    unsigned long max(void) const;

    ///
    /// @brief Returns the smallest number that the random number
    /// generator can return.
    ///
    unsigned long min(void) const;

    ///
    /// @brief Returns the name of the random number generator
    ///
    std::string name(void) const;

    ///
    /// Access underlying GSL random number generator speicific to
    /// current thread. Behaviour of returned generator is undefined
    /// outside current thread.
    ///
    /// @return const pointer to underlying GSL random generator.
    ///
    const gsl_rng* rng(void) const;

    ///
    /// @brief Set the seed \a s for the 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).
    ///
    /// This function will also effect generators created subsequently
    /// in other threads. The seed \a s is saved and subsequent
    /// generators will be created with seed \c s + 1, \c s + 2, etc.
    ///
    /// @see seed_from_devurandom
    ///
    void seed(unsigned long s) const;

    ///
    /// @brief Seed the rng using the /dev/urandom device.
    ///
    /// This function will also effect generators in other threads
    /// created subsequntly (see seed()).
    ///
    /// @return The seed acquired from /dev/urandom.
    ///
    unsigned long seed_from_devurandom(void);

    /**
       \brief Set the state to \a state.

       \return 0 always.

       \note this function only effects the RNG in current thread

       \throw utility::GSL_error on error

       \see gsl_rng_memcpy
    */
    // return int for backward compatibility with yat 0.8
    int set_state(const RNG_state&);

  private:
    RNG(void);

    /**
       \brief Not implemented.

       This copy contructor is not implemented. The constructor is
       declared in order to avoid compiler generated default copy
       constructor.
     */
    RNG(const RNG&);

    /**
       There can be only one RNG so assignment is always
       self-assignment and we do not allow it
    */
    RNG& operator=(const RNG&);

    virtual ~RNG(void);
    void rng_alloc(void) const;

    static RNG* instance_;
    // holds one gsl_rng per thread. Access through rng(void) so a
    // gsl_rng is allocated if necessary.
    mutable boost::thread_specific_ptr<gsl_rng> rng_;
    mutable unsigned long seed_;
  };


  ///
  /// @brief Class holding state of a random generator
  ///
  class RNG_state
  {
  public:
    ///
    /// @brief Constructor
    ///
    explicit RNG_state(const RNG*);

    /**
       Copy Constructor

       \since Explicitely declared since yat 0.5
     */
    RNG_state(const RNG_state&);

    ///
    /// @brief Destructor
    ///
    ~RNG_state(void);

    ///
    /// @return const pointer to underlying GSL random generator.
    ///
    const gsl_rng* rng(void) const;

    /**
       Assignment operator

       \since Explicitely declared since yat 0.5
     */
    RNG_state& operator=(const RNG_state&);

  private:
    gsl_rng* rng_;

    void clone(const gsl_rng&);
  };
    

  // --------------------- Discrete distribtuions ---------------------

  ///
  /// @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
    ///
    Discrete(void);

    ///
    /// @brief The destructor
    ///
    virtual ~Discrete(void);

    ///
    /// @brief Set the seed to \a s.
    ///
    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
    /// default value from the rng's original implementation is used
    /// (cf. GSL documentation).
    ///
    /// \deprecated Provided for backward compatibility with the 0.7
    /// API. Use RNG::instance()->seed(s) instead.
    ///
    void seed(unsigned long s) const YAT_DEPRECATE;

    ///
    /// @brief Set the seed using the /dev/urandom device.
    ///
    /// @return The seed acquired from /dev/urandom.
    ///
    /// \deprecated Provided for backward compatibility with the 0.7
    /// API. Use RNG::instance()->seed_from_devurandom() instead.
    ///
    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;

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

  ///
  /// @brief General
  ///
  class DiscreteGeneral : public Discrete 
  {
  public:
    ///
    /// @brief Constructor
    ///
    /// @param hist is a Histogram defining the probability distribution
    ///
    explicit DiscreteGeneral(const statistics::Histogram& hist);
    
    /**
       \brief Copy constructor

       \since Explicitely implemented in yat 0.5
     */
    DiscreteGeneral(const DiscreteGeneral&);

    ///
    /// @brief Destructor
    ///
    ~DiscreteGeneral(void);

    ///
    /// The generated number is an integer and proportional 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.
    ///
    unsigned long operator()(void) const;

    /**
       \brief Assignment operator

       \since Explicitely implemented in yat 0.5
     */
    DiscreteGeneral& operator=(const DiscreteGeneral&);

  private:
    void free(void);
    void preproc(void);

    gsl_ran_discrete_t* gen_;
    std::vector<double> p_;
  };

  /**
     @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 std::unary_function<unsigned long, unsigned long>
  {
  public:
    /**
       \brief Constructor.

       The generator will generate integers within the range \f$
       [0,n-1] \f$. If \a n is zero, then the whole range of the
       underlying RNG will be used \f$ [min,max] \f$. Setting \a n to
       zero is the preferred way to sample the whole range of the
       underlying RNG, i.e. not setting \n to RNG.max.

       \throw If \a n is larger than the maximum number the underlying
       random number generator can return, then a GSL_error exception
       is thrown.
    */
    explicit DiscreteUniform(unsigned long n=0);

    /**
       \brief Get a random number

       The returned integer is either in the range [RNG.min,RNG.max]
       or [0,n-1] depending on how the random number generator was
       created.

       \see DiscreteUniform(const unsigned long n=0)
    */
    unsigned long operator()(void) const;

    /**
       \brief Get a random integer in the range \f$ [0,n-1] \f$.

       All integers in the range [0,n-1] are equally likely. This
       function should be avoided for sampling the whole range of the
       underlying RNG.

       \throw GSL_error if \a n is larger than the range of the
       underlying generator.
    */
    unsigned long operator()(unsigned long n) const;

  private:
    unsigned long range_;
  };

  /**
     @brief Poisson Distribution
  
     Having a Poisson process (i.e. 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
    ///
    explicit Poisson(const double m=1);

    ///
    /// @return A Poisson distributed number.
    ///
    unsigned long operator()(void) const;

    ///
    /// @return A Poisson distributed number with expectation value
    /// \a m
    ///
    /// @note this operator ignores parameters set in Constructor
    ///
    unsigned long operator()(const double m) const;

  private:
    double m_;
  };

  // --------------------- Continuous distribtuions ---------------------

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

    ///
    /// @brief Constructor
    ///
    Continuous(void);

    ///
    /// @brief The destructor
    ///
    virtual ~Continuous(void);

    ///
    /// @brief Set the seed to \a s.
    ///
    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
    /// default value from the rng's original implementation is used
    /// (cf. GSL documentation).
    ///
    /// \deprecated Provided for backward compatibility with the 0.7
    /// API. Use RNG::instance()->seed(s) instead.
    ///
    void seed(unsigned long s) const YAT_DEPRECATE;

    ///
    /// @brief Set the seed using the /dev/urandom device.
    ///
    /// @return The seed acquired from /dev/urandom.
    ///
    /// \deprecated Provided for backward compatibility with the 0.7
    /// API. Use RNG::instance()->seed_from_devurandom() instead.
    ///
    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;

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

  protected:
    /// pointer to GSL random generator
    RNG* rng_;
  };

  // ContinuousUniform is declared before ContinuousGeneral to avoid
  // forward declaration
  ///
  /// @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:
    double operator()(void) const;
  };

  ///
  /// @brief Generates numbers from a histogram in a continuous manner.
  ///
  class ContinuousGeneral : public Continuous
  {
  public:
    ///
    /// @brief Constructor
    ///
    /// @param hist is a Histogram defining the probability distribution
    ///
    explicit ContinuousGeneral(const statistics::Histogram& hist);

    ///
    /// 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.
    ///
    double operator()(void) const;

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

  /**
     \brief Generator of random numbers from an exponential
     distribution.
     
     The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
     \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
     variance \f$ m^2 \f$
  */
  class Exponential : public Continuous
  {
  public:
    ///
    /// @brief Constructor
    ///
    /// @param m is the expectation value of the distribution.
    ///
    explicit Exponential(const double m=1);

    ///
    /// @return A random number from exponential distribution.
    ///
    double operator()(void) const;

    ///
    /// @return A random number from exponential distribution, with
    /// expectation value \a m
    ///
    /// @note This operator ignores parameters given in constructor.
    ///
    double operator()(const double m) const;

  private:
    double m_;
  };

  /**
     @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
    /// @param m is the expectation value \f$ \mu \f$ of the distribution
    ///
    explicit Gaussian(const double s=1, const double m=0);

    ///
    /// @return A random Gaussian number
    ///
    double operator()(void) const;

    ///
    /// @return A random Gaussian number with standard deviation \a s
    /// and expectation value 0.
    ///
    /// @note this operator ignores parameters given in Constructor
    ///
    double operator()(const double s) const;

    ///
    /// @return A random Gaussian number with standard deviation \a s
    /// and expectation value \a m.
    ///
    /// @note this operator ignores parameters given in Constructor
    ///
    double operator()(const double s, const double m) const;

  private:
    double m_;
    double s_;
  };

  /**
     \brief Convenience function to shuffle a range with singleton RNG.

     Wrapper around std::random_shuffle using DiscreteUniform as
     random generator and thereby using the underlying RNG class,
     which is singleton.

     RandomAccessIterator must be mutable
   */
  template<typename RandomAccessIterator>
  void random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
  {
    typedef RandomAccessIterator rai;
    BOOST_CONCEPT_ASSERT((boost::Mutable_RandomAccessIterator<rai>));
    DiscreteUniform rnd;
    std::random_shuffle(first, last, rnd);
  }

}}} // of namespace random, yat, and theplu

#endif