#ifndef _theplu_yat_random_
#define _theplu_yat_random_
// $Id: random.h 818 2007-03-17 00:55:57Z peter $
/*
Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson
This file is part of the yat library, http://lev.thep.lu.se/trac/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 2 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 this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
02111-1307, USA.
*/
#include
#include
#include
namespace theplu {
namespace yat {
namespace statistics {
class Histogram;
}
namespace random {
//forward declarion
class RNG_state;
///
/// @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 in
/// the GSL manual (Chapter on random number generators). 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.)
///
/// The current implementation is NOT thread safe since the RNG is
/// implemented as a singleton. However, the underlying GSL rng's
/// support thread safety since each instance of GSL rng's keep
/// track of their own state accordning to GSL documentation.
///
/// @see Design Patterns (the singleton and adapter pattern). GSL
/// documentation.
///
class RNG
{
public:
virtual ~RNG(void);
///
/// @brief Get an instance of the random number generator.
///
/// Get an instance of the random number generator. If the random
/// number generator is not already created, the call will create
/// a new generator and use the default seed. The seed must 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.
///
u_long max(void) const;
///
/// @brief Returns the smallest number that the random number
/// generator can return.
///
u_long min(void) const;
///
/// @brief Returns the name of the random number generator
///
std::string name(void) const;
///
/// @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).
///
/// @see seed_from_devurandom
///
void seed(u_long s) const;
///
/// @brief Seed the rng using the /dev/urandom device.
///
/// @return The seed acquired from /dev/urandom.
///
u_long seed_from_devurandom(void);
/**
\brief Set the state to \a state.
\return 0 on success, non-zero otherwise.
\see gsl_rng_memcpy
*/
int set_state(const RNG_state&);
private:
RNG(void);
static RNG* instance_;
gsl_rng* rng_;
};
///
/// @brief Class holding state of a random generator
///
class RNG_state
{
public:
///
/// @brief Constructor
///
RNG_state(const RNG*);
///
/// @brief Destructor
///
~RNG_state(void);
///
/// @return const pointer to underlying GSL random generator.
///
const gsl_rng* rng(void) const;
private:
gsl_rng* 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).
///
/// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
///
void seed(u_long s) const;
///
/// @brief Set the seed using the /dev/urandom device.
///
/// @return The seed acquired from /dev/urandom.
///
/// @see seed, RNG::seed_from_devurandom, RNG::seed
///
u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
///
/// @return A random number.
///
virtual u_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
///
DiscreteGeneral(const statistics::Histogram& hist);
///
/// @brief Destructor
///
~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.
///
u_long operator()(void) const;
private:
gsl_ran_discrete_t* gen_;
};
///
/// @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.
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.
*/
DiscreteUniform(u_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 u_long n=0)
*/
u_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 \an is larger than the range of the
underlying generator.
*/
u_long operator()(u_long n) const;
private:
u_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
///
Poisson(const double m=1);
///
/// @return A Poisson distributed number.
///
u_long operator()(void) const;
///
/// @return A Poisson distributed number with expectation value
/// \a m
///
/// @note this operator ignores parameters set in Constructor
///
u_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).
///
/// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
///
void seed(u_long s) const;
///
/// @brief Set the seed using the /dev/urandom device.
///
/// @return The seed acquired from /dev/urandom.
///
/// @see seed, RNG::seed_from_devurandom, RNG::seed
///
u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
///
/// @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
///
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.
///
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
///
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_;
};
}}} // of namespace random, yat, and theplu
#endif