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

Last change on this file since 4009 was 4009, checked in by Peter, 12 months ago

closes #965; Default constructor for DiscreteGeneral?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 21.0 KB
Line 
1#ifndef _theplu_yat_random_
2#define _theplu_yat_random_
3
4// $Id: random.h 4009 2020-10-19 00:54:47Z peter $
5
6/*
7  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
8  Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2020 Peter Johansson
9
10  This file is part of the yat library, http://dev.thep.lu.se/yat
11
12  The yat library is free software; you can redistribute it and/or
13  modify it under the terms of the GNU General Public License as
14  published by the Free Software Foundation; either version 3 of the
15  License, or (at your option) any later version.
16
17  The yat library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20  General Public License for more details.
21
22  You should have received a copy of the GNU General Public License
23  along with yat. If not, see <http://www.gnu.org/licenses/>.
24*/
25
26#include "yat/utility/config_public.h"
27
28#include "yat/statistics/Histogram.h"
29#include "yat/utility/deprecate.h"
30
31#include <boost/concept_check.hpp>
32
33#include <gsl/gsl_rng.h>
34#include <gsl/gsl_randist.h>
35
36#include <algorithm>
37#include <atomic>
38#include <memory>
39#include <mutex>
40#include <string>
41#include <vector>
42
43namespace theplu {
44namespace yat {
45namespace random {
46
47  //forward declaration
48  class RNG_state;
49
50  ///
51  /// @brief Random Number Generator
52  ///
53  /// The RNG class is wrapper to the GSL random number generator
54  /// (rng). In yat 0.8 (or older) this class provided a single global
55  /// instance of the rng, and made sure there was only one point of
56  /// access to the generator. Since version 0.9 this class provides
57  /// one rng per thread in order to avoid collisions in multi-thread
58  /// applications.
59  ///
60  /// There are many different rng's available in GSL. RNG uses the
61  /// default generator, unless the global variable \c gsl_rng_default
62  /// has been modified (see <a
63  /// href=\gsl_url/Random-number-environment-variables.html>GSL
64  /// Manual</a>). Note, \c gsl_rng_default should be changed before
65  /// RNG creates its generator and safest way to achieve this is to
66  /// modify \c gsl_rng_default prior calling instance() the first
67  /// time.
68  ///
69  /// There is information about how to change seeding and generators
70  /// at run time without recompilation using environment variables in
71  /// the <a
72  /// href=\gsl_url/Random-number-environment-variables.html>GSL
73  /// Manual</a>. RNG supports seeding at compile time if you don't
74  /// want to bother about environment variables and GSL.
75  ///
76  /// The class provides one generator per thread. The first generator
77  /// created is seeded with \c gsl_rng_default_seed and subsequent
78  /// generators are seeded with \c gsl_rng_default_seed + 1, \c
79  /// gsl_rng_default_seed + 2 etc, unless the seed has been modified
80  /// with seed() or seed_from_devurandom().
81  ///
82  /// @see Design Patterns (the singleton and adapter pattern). <a
83  /// href=\gsl_url/Random-Number-Generation.html>GSL documentation</a>.
84  ///
85  class RNG
86  {
87  public:
88
89    ///
90    /// @brief Get an instance of the Random Number Generator.
91    ///
92    /// Get an instance of RNG. If a random number generator is not
93    /// already created for current thread, the call will create a new
94    /// generator of type \c gsl_rng_default. If it is the first
95    /// generator created it will be seeded with \c
96    /// gsl_rng_default_seed; otherwise created generator will be
97    /// seeded with \c seed + \c n, where \c seed is the latest seed
98    /// set (with seed() or seed_from_devurandom()) The seed may be
99    /// changed with the seed or seed_from_devurandom member
100    /// functions.
101    ///
102    /// @return A pointer to the random number generator.
103    ///
104    /// @see seed and seed_from_devurandom
105    ///
106    static RNG* instance(void);
107
108    ///
109    /// @brief Returns the largest number that the random number
110    /// generator can return.
111    ///
112    unsigned long max(void) const;
113
114    ///
115    /// @brief Returns the smallest number that the random number
116    /// generator can return.
117    ///
118    unsigned long min(void) const;
119
120    ///
121    /// @brief Returns the name of the random number generator
122    ///
123    std::string name(void) const;
124
125    ///
126    /// Access underlying GSL random number generator speicific to
127    /// current thread. Behaviour of returned generator is undefined
128    /// outside current thread.
129    ///
130    /// @return const pointer to underlying GSL random generator.
131    ///
132    const gsl_rng* rng(void) const;
133
134    ///
135    /// @brief Set the seed \a s for the rng.
136    ///
137    /// Set the seed \a s for the rng. If \a s is zero, a default
138    /// value from the rng's original implementation is used (cf. GSL
139    /// documentation).
140    ///
141    /// This function will also effect generators created subsequently
142    /// in other threads. The seed \a s is saved and subsequent
143    /// generators will be created with seed \c s + 1, \c s + 2, etc.
144    ///
145    /// @see seed_from_devurandom
146    ///
147    void seed(unsigned long s) const;
148
149    ///
150    /// @brief Seed the rng using the /dev/urandom device.
151    ///
152    /// This function will also effect generators in other threads
153    /// created subsequntly (see seed()).
154    ///
155    /// @return The seed acquired from /dev/urandom.
156    ///
157    unsigned long seed_from_devurandom(void);
158
159    /**
160       \brief Set the state to \a state.
161
162       \return 0 always.
163
164       \note this function only effects the RNG in current thread
165
166       \throw utility::GSL_error on error
167
168       \see gsl_rng_memcpy
169    */
170    // return int for backward compatibility with yat 0.8
171    int set_state(const RNG_state&);
172
173  private:
174    RNG(void);
175
176    /**
177       \brief Not implemented.
178
179       This copy contructor is not implemented. The constructor is
180       declared in order to avoid compiler generated default copy
181       constructor.
182     */
183    RNG(const RNG&);
184
185    /**
186       There can be only one RNG so assignment is always
187       self-assignment and we do not allow it
188    */
189    RNG& operator=(const RNG&);
190
191    virtual ~RNG(void);
192    void rng_alloc(void) const;
193
194    static std::atomic<RNG*> instance_;
195    // holds one gsl_rng per thread. Access through rng(void) so a
196    // gsl_rng is allocated if necessary.
197    static thread_local std::unique_ptr<gsl_rng, void(*)(gsl_rng*)> rng_;
198    mutable unsigned long seed_;
199    // guard needs to be mutable because major mission for it is to protect seed_ against multi-access, and seed_ is mutable...
200    mutable std::mutex mutex_;
201    static std::mutex init_mutex_;
202  };
203
204
205  ///
206  /// @brief Class holding state of a random generator
207  ///
208  class RNG_state
209  {
210  public:
211    ///
212    /// @brief Constructor
213    ///
214    explicit RNG_state(const RNG*);
215
216    /**
217       Copy Constructor
218
219       \since Explicitely declared since yat 0.5
220     */
221    RNG_state(const RNG_state&);
222
223    ///
224    /// @brief Destructor
225    ///
226    ~RNG_state(void);
227
228    ///
229    /// @return const pointer to underlying GSL random generator.
230    ///
231    const gsl_rng* rng(void) const;
232
233    /**
234       Assignment operator
235
236       \since Explicitely declared since yat 0.5
237     */
238    RNG_state& operator=(const RNG_state&);
239
240  private:
241    gsl_rng* rng_;
242
243    void clone(const gsl_rng&);
244  };
245
246
247  // --------------------- Discrete distribtuions ---------------------
248
249  ///
250  /// @brief Discrete random number distributions.
251  ///
252  /// Abstract base class for discrete random number
253  /// distributions. Given K discrete events with different
254  /// probabilities \f$ P[k] \f$, produce a random value k consistent
255  /// with its probability.
256  ///
257  class Discrete
258  {
259  public:
260    /**
261       type returned by operator()
262
263       \since New in yat 0.10
264     */
265    typedef unsigned long int result_type;
266
267    ///
268    /// @brief Constructor
269    ///
270    Discrete(void);
271
272    ///
273    /// @brief The destructor
274    ///
275    virtual ~Discrete(void);
276
277    ///
278    /// @brief Set the seed to \a s.
279    ///
280    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
281    /// default value from the rng's original implementation is used
282    /// (cf. GSL documentation).
283    ///
284    /// \deprecated Provided for backward compatibility with the 0.7
285    /// API. Use RNG::instance()->seed(s) instead.
286    ///
287    void seed(unsigned long s) const YAT_DEPRECATE;
288
289    ///
290    /// @brief Set the seed using the /dev/urandom device.
291    ///
292    /// @return The seed acquired from /dev/urandom.
293    ///
294    /// \deprecated Provided for backward compatibility with the 0.7
295    /// API. Use RNG::instance()->seed_from_devurandom() instead.
296    ///
297    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;
298
299    ///
300    /// @return A random number.
301    ///
302    virtual result_type operator()(void) const = 0;
303
304  protected:
305    /// GSL random gererator
306    RNG* rng_;
307  };
308
309  /**
310     \brief Binomial distribution
311
312     \see gsl_ran_binomial
313
314     \since New in yat 0.10
315   */
316  class Binomial : public Discrete
317  {
318  public:
319    /**
320       \brief Constructor
321
322       Create an object that generates random numbers from binomial
323       distribution, the number of successes in \n trials each with
324       success probability \a p.
325     */
326    Binomial(double p, unsigned int n);
327
328    /**
329       \return number from binomial distrubtion as parametrized in constructor.
330     */
331    unsigned long operator()(void) const;
332  private:
333    double p_;
334    unsigned int n_;
335  };
336
337  ///
338  /// @brief General
339  ///
340  class DiscreteGeneral : public Discrete
341  {
342  public:
343    /**
344       Default constructor leaves the object uninitialized and
345       behaviour of operator() is undefined.
346
347       \since New in yat 0.19
348     */
349    DiscreteGeneral(void);
350
351    ///
352    /// @brief Constructor
353    ///
354    /// @param hist is a Histogram defining the probability distribution
355    ///
356    explicit DiscreteGeneral(const statistics::Histogram& hist);
357
358    /**
359       \brief Copy constructor
360
361       \since Explicitely implemented in yat 0.5
362     */
363    DiscreteGeneral(const DiscreteGeneral&);
364
365    ///
366    /// @brief Destructor
367    ///
368    ~DiscreteGeneral(void);
369
370    ///
371    /// The generated number is an integer and proportional to the
372    /// frequency in the corresponding histogram bin. In other words,
373    /// the probability that 0 is returned is proportinal to the size
374    /// of the first bin.
375    ///
376    /// @return A random number.
377    ///
378    unsigned long operator()(void) const;
379
380    /**
381       \brief Assignment operator
382
383       \since Explicitely implemented in yat 0.5
384     */
385    DiscreteGeneral& operator=(const DiscreteGeneral&);
386
387  private:
388    void free(void);
389    void preproc(void);
390
391    gsl_ran_discrete_t* gen_;
392    std::vector<double> p_;
393  };
394
395  /**
396     @brief Discrete uniform distribution
397
398     Discrete uniform distribution also known as the "equally likely
399     outcomes" distribution. Each outcome, in this case an integer
400     from [0,n-1] , have equal probability to occur.
401
402     Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
403     k < n \f$ \n
404     Expectation value: \f$ \frac{n-1}{2} \f$ \n
405     Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
406  */
407  class DiscreteUniform : public Discrete
408  {
409  public:
410    /// argument type is unsigned long int
411    typedef unsigned long argument_type;
412
413    /**
414       \brief Constructor.
415
416       The generator will generate integers within the range \f$
417       [0,n-1] \f$. If \a n is zero, then the whole range of the
418       underlying RNG will be used \f$ [min,max] \f$. Setting \a n to
419       zero is the preferred way to sample the whole range of the
420       underlying RNG, i.e. not setting \n to RNG.max.
421
422       \throw If \a n is larger than the maximum number the underlying
423       random number generator can return, then a GSL_error exception
424       is thrown.
425    */
426    explicit DiscreteUniform(unsigned long n=0);
427
428    /**
429       \brief Get a random number
430
431       The returned integer is either in the range [RNG.min,RNG.max]
432       or [0,n-1] depending on how the random number generator was
433       created.
434
435       \see DiscreteUniform(const unsigned long n=0)
436    */
437    unsigned long operator()(void) const;
438
439    /**
440       \brief Get a random integer in the range \f$ [0,n-1] \f$.
441
442       All integers in the range [0,n-1] are equally likely. This
443       function should be avoided for sampling the whole range of the
444       underlying RNG.
445
446       \throw GSL_error if \a n is larger than the range of the
447       underlying generator.
448    */
449    unsigned long operator()(unsigned long n) const;
450
451    /**
452       For n>0 the min value is 0. For n=0 (default) the min value
453       is RNG::min(), typically zero but depending on generator in use.
454
455       \return smallest possible value
456
457       \since New in yat 0.18
458     */
459    unsigned long int min(void) const;
460
461    /**
462       For n>0 the max value is n-1. For n=0 (default) the max value
463       is RNG::max().
464
465       \return maximal value that class returns
466
467       \since New in yat 0.18
468     */
469    unsigned long int max(void) const;
470
471  private:
472    unsigned long range_;
473  };
474
475
476  /**
477     @brief Geomtric Distribution
478
479     The number of independent trials with probability \em p until the
480     first success.
481
482     Probability function \f$ p(k) = p (1-p)^(k-1) \f$
483
484     \since New in yat 0.14
485  */
486  class Geometric : public Discrete
487  {
488  public:
489    /**
490       \brief Constructor
491
492       \param p is probability for success in one trial
493    */
494    Geometric(double p);
495
496    /*
497      \return a number from Geomtric distribution
498    */
499    unsigned long operator()(void) const;
500
501    /**
502       \return a number from Geomtric distribution with success rate \a p
503
504       \note this operator ignores parameters set in Constructor
505    */
506    unsigned long operator()(double p) const;
507
508  private:
509    double p_;
510  };
511
512
513  /**
514     If we have \a n1 samples of type 1 and \a n2 samples of type 2
515     and draw \a t samples with replacement, number of drawn samples
516     of type 1 will follow the hyper geometric distribution.
517
518     \since New in yat 0.14
519   */
520  class HyperGeometric : public Discrete
521  {
522  public:
523    /**
524       \brief Defaul constructor
525     */
526    HyperGeometric(void);
527
528    /**
529       \brief Constructor
530       \param n1 number of samples of type 1
531       \param n2 number of samples of type 2
532       \param t number of samples to draw
533     */
534    HyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
535
536    /**
537       \return random number from hypergeometric distribution using
538       parameters set in constructor.
539     */
540    unsigned long int operator()(void) const;
541
542    /**
543       \return random number from hypergeometric distribution using
544       parameters passed.
545     */
546    unsigned long int operator()(unsigned int n1, unsigned int n2,
547                                 unsigned int t) const;
548  private:
549    unsigned int n1_;
550    unsigned int n2_;
551    unsigned int t_;
552  };
553
554  /**
555     We have \a n1 samples of type 1 and \a n2 samples of type 2.
556     Samples are drawn with replacement until \a t samles of type 2
557     are drawn. Then \a k, number of drawn samples of type 1, follows
558     negative hypergeometric distribution.
559
560     \since New in yat 0.14
561
562     \see HyperGeometric
563   */
564  class NegativeHyperGeometric : public Discrete
565  {
566  public:
567    /**
568       \brief Default constructor
569     */
570    NegativeHyperGeometric(void);
571
572    /**
573       \brief Constructor
574       \param n1 number of samples of type 1
575       \param n2 number of samples of type 2
576       \param t number of samples of type 2 to draw
577     */
578    NegativeHyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
579
580    /**
581       \return random number from negative hypergeometric distribution
582       using parameters set in constructor.
583     */
584    unsigned long int operator()(void) const;
585
586    /**
587       \return random number from negative hypergeometric distribution
588       using parameters passed.
589     */
590    unsigned long int operator()(unsigned int n1, unsigned int n2,
591                                 unsigned int t) const;
592  private:
593    unsigned int n1_;
594    unsigned int n2_;
595    unsigned int t_;
596  };
597
598
599  /**
600     @brief Poisson Distribution
601
602     Having a Poisson process (i.e. no memory), number of occurences
603     within a given time window is Poisson distributed. This
604     distribution is the limit of a Binomial distribution when number
605     of attempts is large, and the probability for one attempt to be
606     succesful is small (in such a way that the expected number of
607     succesful attempts is \f$ m \f$.
608
609     Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
610     k  \f$ \n
611     Expectation value: \f$ m \f$ \n
612     Variance: \f$ m \f$
613  */
614  class Poisson : public Discrete
615  {
616  public:
617    ///
618    /// @brief Constructor
619    ///
620    /// @param m is expectation value
621    ///
622    explicit Poisson(const double m=1);
623
624    ///
625    /// @return A Poisson distributed number.
626    ///
627    unsigned long operator()(void) const;
628
629    ///
630    /// @return A Poisson distributed number with expectation value
631    /// \a m
632    ///
633    /// @note this operator ignores parameters set in Constructor
634    ///
635    unsigned long operator()(const double m) const;
636
637  private:
638    double m_;
639  };
640
641  // --------------------- Continuous distribtuions ---------------------
642
643  ///
644  /// @brief Continuous random number distributions.
645  ///
646  /// Abstract base class for continuous random number distributions.
647  ///
648  class Continuous
649  {
650  public:
651    /**
652       type returned by operator()
653
654       \since New in yat 0.10
655     */
656    typedef double result_type;
657
658    ///
659    /// @brief Constructor
660    ///
661    Continuous(void);
662
663    ///
664    /// @brief The destructor
665    ///
666    virtual ~Continuous(void);
667
668    ///
669    /// @brief Set the seed to \a s.
670    ///
671    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
672    /// default value from the rng's original implementation is used
673    /// (cf. GSL documentation).
674    ///
675    /// \deprecated Provided for backward compatibility with the 0.7
676    /// API. Use RNG::instance()->seed(s) instead.
677    ///
678    void seed(unsigned long s) const YAT_DEPRECATE;
679
680    ///
681    /// @brief Set the seed using the /dev/urandom device.
682    ///
683    /// @return The seed acquired from /dev/urandom.
684    ///
685    /// \deprecated Provided for backward compatibility with the 0.7
686    /// API. Use RNG::instance()->seed_from_devurandom() instead.
687    ///
688    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;
689
690    ///
691    /// @return A random number
692    ///
693    virtual result_type operator()(void) const = 0;
694
695  protected:
696    /// pointer to GSL random generator
697    RNG* rng_;
698  };
699
700  // ContinuousUniform is declared before ContinuousGeneral to avoid
701  // forward declaration
702  ///
703  /// @brief Uniform distribution
704  ///
705  /// Class for generating a random number from a uniform distribution
706  /// in the range [0,1), i.e. zero is included but not 1.
707  ///
708  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
709  /// Expectation value: 0.5 \n
710  /// Variance: \f$ \frac{1}{12} \f$
711  ///
712  class ContinuousUniform : public Continuous
713  {
714  public:
715    double operator()(void) const;
716  };
717
718  ///
719  /// @brief Generates numbers from a histogram in a continuous manner.
720  ///
721  class ContinuousGeneral : public Continuous
722  {
723  public:
724    ///
725    /// @brief Constructor
726    ///
727    /// @param hist is a Histogram defining the probability distribution
728    ///
729    explicit ContinuousGeneral(const statistics::Histogram& hist);
730
731    ///
732    /// The number is generated in a two step process. First the bin
733    /// in the histogram is randomly selected (see
734    /// DiscreteGeneral). Then a number is generated uniformly from
735    /// the interval defined by the bin.
736    ///
737    /// @return A random number.
738    ///
739    double operator()(void) const;
740
741  private:
742    const DiscreteGeneral discrete_;
743    const statistics::Histogram hist_;
744    ContinuousUniform u_;
745  };
746
747  /**
748     \brief Generator of random numbers from an exponential
749     distribution.
750
751     The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
752     \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
753     variance \f$ m^2 \f$
754  */
755  class Exponential : public Continuous
756  {
757  public:
758    ///
759    /// @brief Constructor
760    ///
761    /// @param m is the expectation value of the distribution.
762    ///
763    explicit Exponential(const double m=1);
764
765    ///
766    /// @return A random number from exponential distribution.
767    ///
768    double operator()(void) const;
769
770    ///
771    /// @return A random number from exponential distribution, with
772    /// expectation value \a m
773    ///
774    /// @note This operator ignores parameters given in constructor.
775    ///
776    double operator()(const double m) const;
777
778  private:
779    double m_;
780  };
781
782  /**
783     @brief Gaussian distribution
784
785     Class for generating a random number from a Gaussian distribution
786     between zero and unity. Utilizes the Box-Muller algorithm, which
787     needs two calls to random generator.
788
789     Distribution function \f$ f(x) =
790     \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
791     \f$ \n
792     Expectation value: \f$ \mu \f$ \n
793     Variance: \f$ \sigma^2 \f$
794  */
795  class Gaussian : public Continuous
796  {
797  public:
798    ///
799    /// @brief Constructor
800    ///
801    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
802    /// @param m is the expectation value \f$ \mu \f$ of the distribution
803    ///
804    explicit Gaussian(const double s=1, const double m=0);
805
806    ///
807    /// @return A random Gaussian number
808    ///
809    double operator()(void) const;
810
811    ///
812    /// @return A random Gaussian number with standard deviation \a s
813    /// and expectation value 0.
814    ///
815    /// @note this operator ignores parameters given in Constructor
816    ///
817    double operator()(const double s) const;
818
819    ///
820    /// @return A random Gaussian number with standard deviation \a s
821    /// and expectation value \a m.
822    ///
823    /// @note this operator ignores parameters given in Constructor
824    ///
825    double operator()(const double s, const double m) const;
826
827  private:
828    double m_;
829    double s_;
830  };
831
832  /**
833     \brief Convenience function to shuffle a range with singleton RNG.
834
835     Wrapper around std::random_shuffle using DiscreteUniform as
836     random generator and thereby using the underlying RNG class,
837     which is singleton.
838
839     Type Requirements:
840     - RandomAccessIterator is \random_access_iterator
841   */
842  template<typename RandomAccessIterator>
843  void random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
844  {
845    DiscreteUniform rnd;
846    std::random_shuffle(first, last, rnd);
847  }
848
849}}} // of namespace random, yat, and theplu
850
851#endif
Note: See TracBrowser for help on using the repository browser.