source: trunk/yat/random/random.h

Last change on this file was 4010, checked in by Peter, 11 months ago

Implement move constructo and move-assignment operator for
DiscreteGeneral? (closes #964). Using swap idiom also for
copy-assignment as it gives stronger exception safety, i.e., if
preproc throws, lhs is not modified (whereas before ::p_ was).

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 21.5 KB
Line 
1#ifndef _theplu_yat_random_
2#define _theplu_yat_random_
3
4// $Id: random.h 4010 2020-10-19 02:33: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 Move constructor
367
368       \since New in yat 0.19
369     */
370    DiscreteGeneral(DiscreteGeneral&&);
371
372    ///
373    /// @brief Destructor
374    ///
375    ~DiscreteGeneral(void);
376
377    ///
378    /// The generated number is an integer and proportional to the
379    /// frequency in the corresponding histogram bin. In other words,
380    /// the probability that 0 is returned is proportinal to the size
381    /// of the first bin.
382    ///
383    /// @return A random number.
384    ///
385    unsigned long operator()(void) const;
386
387    /**
388       \brief Assignment operator
389
390       \since Explicitely implemented in yat 0.5
391     */
392    DiscreteGeneral& operator=(const DiscreteGeneral&);
393
394    /**
395       \brief Move assignment operator
396
397       \since New in yat 0.19
398     */
399    DiscreteGeneral& operator=(DiscreteGeneral&&);
400
401  private:
402    void free(void);
403    void preproc(void);
404
405    gsl_ran_discrete_t* gen_;
406    std::vector<double> p_;
407    friend void swap(DiscreteGeneral& lhs, DiscreteGeneral& rhs);
408  };
409
410  /**
411     \relates DiscreteGeneral
412
413     \since Specialisation new in yat 0.19
414   */
415  void swap(DiscreteGeneral& lhs, DiscreteGeneral& rhs);
416
417  /**
418     @brief Discrete uniform distribution
419
420     Discrete uniform distribution also known as the "equally likely
421     outcomes" distribution. Each outcome, in this case an integer
422     from [0,n-1] , have equal probability to occur.
423
424     Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
425     k < n \f$ \n
426     Expectation value: \f$ \frac{n-1}{2} \f$ \n
427     Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
428  */
429  class DiscreteUniform : public Discrete
430  {
431  public:
432    /// argument type is unsigned long int
433    typedef unsigned long argument_type;
434
435    /**
436       \brief Constructor.
437
438       The generator will generate integers within the range \f$
439       [0,n-1] \f$. If \a n is zero, then the whole range of the
440       underlying RNG will be used \f$ [min,max] \f$. Setting \a n to
441       zero is the preferred way to sample the whole range of the
442       underlying RNG, i.e. not setting \n to RNG.max.
443
444       \throw If \a n is larger than the maximum number the underlying
445       random number generator can return, then a GSL_error exception
446       is thrown.
447    */
448    explicit DiscreteUniform(unsigned long n=0);
449
450    /**
451       \brief Get a random number
452
453       The returned integer is either in the range [RNG.min,RNG.max]
454       or [0,n-1] depending on how the random number generator was
455       created.
456
457       \see DiscreteUniform(const unsigned long n=0)
458    */
459    unsigned long operator()(void) const;
460
461    /**
462       \brief Get a random integer in the range \f$ [0,n-1] \f$.
463
464       All integers in the range [0,n-1] are equally likely. This
465       function should be avoided for sampling the whole range of the
466       underlying RNG.
467
468       \throw GSL_error if \a n is larger than the range of the
469       underlying generator.
470    */
471    unsigned long operator()(unsigned long n) const;
472
473    /**
474       For n>0 the min value is 0. For n=0 (default) the min value
475       is RNG::min(), typically zero but depending on generator in use.
476
477       \return smallest possible value
478
479       \since New in yat 0.18
480     */
481    unsigned long int min(void) const;
482
483    /**
484       For n>0 the max value is n-1. For n=0 (default) the max value
485       is RNG::max().
486
487       \return maximal value that class returns
488
489       \since New in yat 0.18
490     */
491    unsigned long int max(void) const;
492
493  private:
494    unsigned long range_;
495  };
496
497
498  /**
499     @brief Geomtric Distribution
500
501     The number of independent trials with probability \em p until the
502     first success.
503
504     Probability function \f$ p(k) = p (1-p)^(k-1) \f$
505
506     \since New in yat 0.14
507  */
508  class Geometric : public Discrete
509  {
510  public:
511    /**
512       \brief Constructor
513
514       \param p is probability for success in one trial
515    */
516    Geometric(double p);
517
518    /*
519      \return a number from Geomtric distribution
520    */
521    unsigned long operator()(void) const;
522
523    /**
524       \return a number from Geomtric distribution with success rate \a p
525
526       \note this operator ignores parameters set in Constructor
527    */
528    unsigned long operator()(double p) const;
529
530  private:
531    double p_;
532  };
533
534
535  /**
536     If we have \a n1 samples of type 1 and \a n2 samples of type 2
537     and draw \a t samples with replacement, number of drawn samples
538     of type 1 will follow the hyper geometric distribution.
539
540     \since New in yat 0.14
541   */
542  class HyperGeometric : public Discrete
543  {
544  public:
545    /**
546       \brief Defaul constructor
547     */
548    HyperGeometric(void);
549
550    /**
551       \brief Constructor
552       \param n1 number of samples of type 1
553       \param n2 number of samples of type 2
554       \param t number of samples to draw
555     */
556    HyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
557
558    /**
559       \return random number from hypergeometric distribution using
560       parameters set in constructor.
561     */
562    unsigned long int operator()(void) const;
563
564    /**
565       \return random number from hypergeometric distribution using
566       parameters passed.
567     */
568    unsigned long int operator()(unsigned int n1, unsigned int n2,
569                                 unsigned int t) const;
570  private:
571    unsigned int n1_;
572    unsigned int n2_;
573    unsigned int t_;
574  };
575
576  /**
577     We have \a n1 samples of type 1 and \a n2 samples of type 2.
578     Samples are drawn with replacement until \a t samles of type 2
579     are drawn. Then \a k, number of drawn samples of type 1, follows
580     negative hypergeometric distribution.
581
582     \since New in yat 0.14
583
584     \see HyperGeometric
585   */
586  class NegativeHyperGeometric : public Discrete
587  {
588  public:
589    /**
590       \brief Default constructor
591     */
592    NegativeHyperGeometric(void);
593
594    /**
595       \brief Constructor
596       \param n1 number of samples of type 1
597       \param n2 number of samples of type 2
598       \param t number of samples of type 2 to draw
599     */
600    NegativeHyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
601
602    /**
603       \return random number from negative hypergeometric distribution
604       using parameters set in constructor.
605     */
606    unsigned long int operator()(void) const;
607
608    /**
609       \return random number from negative hypergeometric distribution
610       using parameters passed.
611     */
612    unsigned long int operator()(unsigned int n1, unsigned int n2,
613                                 unsigned int t) const;
614  private:
615    unsigned int n1_;
616    unsigned int n2_;
617    unsigned int t_;
618  };
619
620
621  /**
622     @brief Poisson Distribution
623
624     Having a Poisson process (i.e. no memory), number of occurences
625     within a given time window is Poisson distributed. This
626     distribution is the limit of a Binomial distribution when number
627     of attempts is large, and the probability for one attempt to be
628     succesful is small (in such a way that the expected number of
629     succesful attempts is \f$ m \f$.
630
631     Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
632     k  \f$ \n
633     Expectation value: \f$ m \f$ \n
634     Variance: \f$ m \f$
635  */
636  class Poisson : public Discrete
637  {
638  public:
639    ///
640    /// @brief Constructor
641    ///
642    /// @param m is expectation value
643    ///
644    explicit Poisson(const double m=1);
645
646    ///
647    /// @return A Poisson distributed number.
648    ///
649    unsigned long operator()(void) const;
650
651    ///
652    /// @return A Poisson distributed number with expectation value
653    /// \a m
654    ///
655    /// @note this operator ignores parameters set in Constructor
656    ///
657    unsigned long operator()(const double m) const;
658
659  private:
660    double m_;
661  };
662
663  // --------------------- Continuous distribtuions ---------------------
664
665  ///
666  /// @brief Continuous random number distributions.
667  ///
668  /// Abstract base class for continuous random number distributions.
669  ///
670  class Continuous
671  {
672  public:
673    /**
674       type returned by operator()
675
676       \since New in yat 0.10
677     */
678    typedef double result_type;
679
680    ///
681    /// @brief Constructor
682    ///
683    Continuous(void);
684
685    ///
686    /// @brief The destructor
687    ///
688    virtual ~Continuous(void);
689
690    ///
691    /// @brief Set the seed to \a s.
692    ///
693    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
694    /// default value from the rng's original implementation is used
695    /// (cf. GSL documentation).
696    ///
697    /// \deprecated Provided for backward compatibility with the 0.7
698    /// API. Use RNG::instance()->seed(s) instead.
699    ///
700    void seed(unsigned long s) const YAT_DEPRECATE;
701
702    ///
703    /// @brief Set the seed using the /dev/urandom device.
704    ///
705    /// @return The seed acquired from /dev/urandom.
706    ///
707    /// \deprecated Provided for backward compatibility with the 0.7
708    /// API. Use RNG::instance()->seed_from_devurandom() instead.
709    ///
710    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;
711
712    ///
713    /// @return A random number
714    ///
715    virtual result_type operator()(void) const = 0;
716
717  protected:
718    /// pointer to GSL random generator
719    RNG* rng_;
720  };
721
722  // ContinuousUniform is declared before ContinuousGeneral to avoid
723  // forward declaration
724  ///
725  /// @brief Uniform distribution
726  ///
727  /// Class for generating a random number from a uniform distribution
728  /// in the range [0,1), i.e. zero is included but not 1.
729  ///
730  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
731  /// Expectation value: 0.5 \n
732  /// Variance: \f$ \frac{1}{12} \f$
733  ///
734  class ContinuousUniform : public Continuous
735  {
736  public:
737    double operator()(void) const;
738  };
739
740  ///
741  /// @brief Generates numbers from a histogram in a continuous manner.
742  ///
743  class ContinuousGeneral : public Continuous
744  {
745  public:
746    ///
747    /// @brief Constructor
748    ///
749    /// @param hist is a Histogram defining the probability distribution
750    ///
751    explicit ContinuousGeneral(const statistics::Histogram& hist);
752
753    ///
754    /// The number is generated in a two step process. First the bin
755    /// in the histogram is randomly selected (see
756    /// DiscreteGeneral). Then a number is generated uniformly from
757    /// the interval defined by the bin.
758    ///
759    /// @return A random number.
760    ///
761    double operator()(void) const;
762
763  private:
764    const DiscreteGeneral discrete_;
765    const statistics::Histogram hist_;
766    ContinuousUniform u_;
767  };
768
769  /**
770     \brief Generator of random numbers from an exponential
771     distribution.
772
773     The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
774     \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
775     variance \f$ m^2 \f$
776  */
777  class Exponential : public Continuous
778  {
779  public:
780    ///
781    /// @brief Constructor
782    ///
783    /// @param m is the expectation value of the distribution.
784    ///
785    explicit Exponential(const double m=1);
786
787    ///
788    /// @return A random number from exponential distribution.
789    ///
790    double operator()(void) const;
791
792    ///
793    /// @return A random number from exponential distribution, with
794    /// expectation value \a m
795    ///
796    /// @note This operator ignores parameters given in constructor.
797    ///
798    double operator()(const double m) const;
799
800  private:
801    double m_;
802  };
803
804  /**
805     @brief Gaussian distribution
806
807     Class for generating a random number from a Gaussian distribution
808     between zero and unity. Utilizes the Box-Muller algorithm, which
809     needs two calls to random generator.
810
811     Distribution function \f$ f(x) =
812     \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
813     \f$ \n
814     Expectation value: \f$ \mu \f$ \n
815     Variance: \f$ \sigma^2 \f$
816  */
817  class Gaussian : public Continuous
818  {
819  public:
820    ///
821    /// @brief Constructor
822    ///
823    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
824    /// @param m is the expectation value \f$ \mu \f$ of the distribution
825    ///
826    explicit Gaussian(const double s=1, const double m=0);
827
828    ///
829    /// @return A random Gaussian number
830    ///
831    double operator()(void) const;
832
833    ///
834    /// @return A random Gaussian number with standard deviation \a s
835    /// and expectation value 0.
836    ///
837    /// @note this operator ignores parameters given in Constructor
838    ///
839    double operator()(const double s) const;
840
841    ///
842    /// @return A random Gaussian number with standard deviation \a s
843    /// and expectation value \a m.
844    ///
845    /// @note this operator ignores parameters given in Constructor
846    ///
847    double operator()(const double s, const double m) const;
848
849  private:
850    double m_;
851    double s_;
852  };
853
854  /**
855     \brief Convenience function to shuffle a range with singleton RNG.
856
857     Wrapper around std::random_shuffle using DiscreteUniform as
858     random generator and thereby using the underlying RNG class,
859     which is singleton.
860
861     Type Requirements:
862     - RandomAccessIterator is \random_access_iterator
863   */
864  template<typename RandomAccessIterator>
865  void random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
866  {
867    DiscreteUniform rnd;
868    std::random_shuffle(first, last, rnd);
869  }
870
871}}} // of namespace random, yat, and theplu
872
873#endif
Note: See TracBrowser for help on using the repository browser.