#ifndef _theplu_yat_random_
#define _theplu_yat_random_
// $Id: random.h 4010 2020-10-19 02:33:47Z peter $
/*
Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2020 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 .
*/
#include "yat/utility/config_public.h"
#include "yat/statistics/Histogram.h"
#include "yat/utility/deprecate.h"
#include
#include
#include
#include
#include
#include
#include
#include
#include
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 GSL
/// Manual). 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 GSL
/// Manual. 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). GSL documentation.
///
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 std::atomic instance_;
// holds one gsl_rng per thread. Access through rng(void) so a
// gsl_rng is allocated if necessary.
static thread_local std::unique_ptr 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 std::mutex mutex_;
static std::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:
/**
Default constructor leaves the object uninitialized and
behaviour of operator() is undefined.
\since New in yat 0.19
*/
DiscreteGeneral(void);
///
/// @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 Move constructor
\since New in yat 0.19
*/
DiscreteGeneral(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&);
/**
\brief Move assignment operator
\since New in yat 0.19
*/
DiscreteGeneral& operator=(DiscreteGeneral&&);
private:
void free(void);
void preproc(void);
gsl_ran_discrete_t* gen_;
std::vector p_;
friend void swap(DiscreteGeneral& lhs, DiscreteGeneral& rhs);
};
/**
\relates DiscreteGeneral
\since Specialisation new in yat 0.19
*/
void swap(DiscreteGeneral& lhs, DiscreteGeneral& rhs);
/**
@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:
/// argument type is unsigned long int
typedef unsigned long argument_type;
/**
\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;
/**
For n>0 the min value is 0. For n=0 (default) the min value
is RNG::min(), typically zero but depending on generator in use.
\return smallest possible value
\since New in yat 0.18
*/
unsigned long int min(void) const;
/**
For n>0 the max value is n-1. For n=0 (default) the max value
is RNG::max().
\return maximal value that class returns
\since New in yat 0.18
*/
unsigned long int max(void) 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
void random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
{
DiscreteUniform rnd;
std::random_shuffle(first, last, rnd);
}
}}} // of namespace random, yat, and theplu
#endif