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

#ifndef _theplu_yat_random_
#define _theplu_yat_random_

// $Id: random.h 3579 2017-01-16 03:54:43Z peter $

/*
  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017 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"

// always include this before <boost/excpetion_ptr.hpp> indirectly below
#include "yat/utility/boost_exception_ptr.h"

#include <boost/concept_check.hpp>
#include <boost/iterator/iterator_concepts.hpp>
#include <boost/thread.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_;
    // guard needs to be mutable because major mission for it is to protect seed_ against multi-access, and seed_ is mutable...
    mutable boost::mutex mutex_;
    static boost::mutex init_mutex_;
  };


  ///
  /// @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:
    /**
       type returned by operator()

       \since New in yat 0.10
     */
    typedef unsigned long int result_type;

    ///
    /// @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 result_type operator()(void) const = 0;

  protected:
    /// GSL random gererator
    RNG* rng_;
  };

  /**
     \brief Binomial distribution

     \see gsl_ran_binomial

     \since New in yat 0.10
   */
  class Binomial : public Discrete
  {
  public:
    /**
       \brief Constructor

       Create an object that generates random numbers from binomial
       distribution, the number of successes in \n trials each with
       success probability \a p.
     */
    Binomial(double p, unsigned int n);

    /**
       \return number from binomial distrubtion as parametrized in constructor.
     */
    unsigned long operator()(void) const;
  private:
    double p_;
    unsigned int n_;
  };

  ///
  /// @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 Geomtric Distribution

     The number of independent trials with probability \em p until the
     first success.

     Probability function \f$ p(k) = p (1-p)^(k-1) \f$

     \since New in yat 0.14
  */
  class Geometric : public Discrete
  {
  public:
    /**
       \brief Constructor

       \param p is probability for success in one trial
    */
    Geometric(double p);

    /*
      \return a number from Geomtric distribution
    */
    unsigned long operator()(void) const;

    /**
       \return a number from Geomtric distribution with success rate \a p

       \note this operator ignores parameters set in Constructor
    */
    unsigned long operator()(double p) const;

  private:
    double p_;
  };


  /**
     If we have \a n1 samples of type 1 and \a n2 samples of type 2
     and draw \a t samples with replacement, number of drawn samples
     of type 1 will follow the hyper geometric distribution.

     \since New in yat 0.14
   */
  class HyperGeometric : public Discrete
  {
  public:
    /**
       \brief Defaul constructor
     */
    HyperGeometric(void);

    /**
       \brief Constructor
       \param n1 number of samples of type 1
       \param n2 number of samples of type 2
       \param t number of samples to draw
     */
    HyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);

    /**
       \return random number from hypergeometric distribution using
       parameters set in constructor.
     */
    unsigned long int operator()(void) const;

    /**
       \return random number from hypergeometric distribution using
       parameters passed.
     */
    unsigned long int operator()(unsigned int n1, unsigned int n2,
                                 unsigned int t) const;
  private:
    unsigned int n1_;
    unsigned int n2_;
    unsigned int t_;
  };

  /**
     We have \a n1 samples of type 1 and \a n2 samples of type 2.
     Samples are drawn with replacement until \a t samles of type 2
     are drawn. Then \a k, number of drawn samples of type 1, follows
     negative hypergeometric distribution.

     \since New in yat 0.14

     \see HyperGeometric
   */
  class NegativeHyperGeometric : public Discrete
  {
  public:
    /**
       \brief Default constructor
     */
    NegativeHyperGeometric(void);

    /**
       \brief Constructor
       \param n1 number of samples of type 1
       \param n2 number of samples of type 2
       \param t number of samples of type 2 to draw
     */
    NegativeHyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);

    /**
       \return random number from negative hypergeometric distribution
       using parameters set in constructor.
     */
    unsigned long int operator()(void) const;

    /**
       \return random number from negative hypergeometric distribution
       using parameters passed.
     */
    unsigned long int operator()(unsigned int n1, unsigned int n2,
                                 unsigned int t) const;
  private:
    unsigned int n1_;
    unsigned int n2_;
    unsigned int t_;
  };


  /**
     @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:
    /**
       type returned by operator()

       \since New in yat 0.10
     */
    typedef double result_type;

    ///
    /// @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 result_type 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.

     Type Requirements:
     - RandomAccessIterator is \random_access_iterator
   */
  template<typename RandomAccessIterator>
  void random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
  {
    DiscreteUniform rnd;
    std::random_shuffle(first, last, rnd);
  }

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

#endif