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

Last change on this file since 2747 was 2747, checked in by Peter, 9 years ago

doc typo

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 15.5 KB
Line 
1#ifndef _theplu_yat_random_
2#define _theplu_yat_random_
3
4// $Id: random.h 2747 2012-06-18 05:32:07Z peter $
5
6/*
7  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
8  Copyright (C) 2009, 2010, 2011 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/statistics/Histogram.h"
27#include "yat/utility/deprecate.h"
28
29#include <boost/concept_check.hpp>
30
31#include <gsl/gsl_rng.h>
32#include <gsl/gsl_randist.h>
33
34#include <algorithm>
35#include <string>
36#include <vector>
37
38namespace theplu {
39namespace yat {
40namespace random {
41
42  //forward declaration
43  class RNG_state;
44
45  ///
46  /// @brief Random Number Generator
47  ///
48  /// The RNG class is wrapper to the GSL random number generator
49  /// (rng). This class provides a single global instance of the rng,
50  /// and makes sure there is only one point of access to the
51  /// generator.
52  ///
53  /// There is information about how to change seeding and generators
54  /// at run time without recompilation using environment variables in
55  /// the GSL manual (Chapter on random number generators). RNG of
56  /// course support seeding at compile time if you don't want to
57  /// bother about environment variables and GSL.
58  ///
59  /// There are many different rng's available in GSL. Currently only
60  /// the default generator is implemented and no other one is
61  /// choosable through the class interface. This means that you have
62  /// to fall back to the use of environment variables as described in
63  /// the GSL documentation, or be bold and request support for other
64  /// rng's through the class interface.
65  ///
66  /// Not all GSL functionality is implemented, we'll add
67  /// functionality when needed and may do it when requested. Better
68  /// yet, supply us with code and we will probably add it to the code
69  /// (BUT remember to implement reasonable tests for your code and
70  /// follow the coding style.)
71  ///
72  /// The current implementation is NOT thread safe since the RNG is
73  /// implemented as a singleton. However, the underlying GSL rng's
74  /// support thread safety since each instance of GSL rng's keep
75  /// track of their own state accordning to GSL documentation.
76  ///
77  /// @see Design Patterns (the singleton and adapter pattern). GSL
78  /// documentation.
79  ///
80  class RNG
81  {
82  public:
83
84    ///
85    /// @brief Get an instance of the random number generator.
86    ///
87    /// Get an instance of the random number generator. If the random
88    /// number generator is not already created, the call will create
89    /// a new generator and use the default seed. The seed must be
90    /// changed with the seed or seed_from_devurandom member
91    /// functions.
92    ///
93    /// @return A pointer to the random number generator.
94    ///
95    /// @see seed and seed_from_devurandom
96    ///
97    static RNG* instance(void);
98
99    ///
100    /// @brief Returns the largest number that the random number
101    /// generator can return.
102    ///
103    unsigned long max(void) const;
104
105    ///
106    /// @brief Returns the smallest number that the random number
107    /// generator can return.
108    ///
109    unsigned long min(void) const;
110
111    ///
112    /// @brief Returns the name of the random number generator
113    ///
114    std::string name(void) const;
115
116    ///
117    /// @return const pointer to underlying GSL random generator.
118    ///
119    const gsl_rng* rng(void) const;
120
121    ///
122    /// @brief Set the seed \a s for the rng.
123    ///
124    /// Set the seed \a s for the rng. If \a s is zero, a default
125    /// value from the rng's original implementation is used (cf. GSL
126    /// documentation).
127    ///
128    /// @see seed_from_devurandom
129    ///
130    void seed(unsigned long s) const;
131
132    ///
133    /// @brief Seed the rng using the /dev/urandom device.
134    ///
135    /// @return The seed acquired from /dev/urandom.
136    ///
137    unsigned long seed_from_devurandom(void);
138
139    /**
140       \brief Set the state to \a state.
141
142       \return 0 on success, non-zero otherwise.
143
144       \see gsl_rng_memcpy
145    */
146    int set_state(const RNG_state&);
147
148  private:
149    RNG(void);
150
151    /**
152       \brief Not implemented.
153
154       This copy contructor is not implemented. The constructor is
155       declared in order to avoid compiler generated default copy
156       constructor.
157     */
158    RNG(const RNG&);
159
160    /**
161       There can be only one RNG so assignment is always
162       self-assignment and we do not allow it
163    */
164    RNG& operator=(const RNG&);
165
166    virtual ~RNG(void);
167
168    static RNG* instance_;
169    gsl_rng* rng_;
170  };
171
172
173  ///
174  /// @brief Class holding state of a random generator
175  ///
176  class RNG_state
177  {
178  public:
179    ///
180    /// @brief Constructor
181    ///
182    explicit RNG_state(const RNG*);
183
184    /**
185       Copy Constructor
186
187       \since Explicitely declared since yat 0.5
188     */
189    RNG_state(const RNG_state&);
190
191    ///
192    /// @brief Destructor
193    ///
194    ~RNG_state(void);
195
196    ///
197    /// @return const pointer to underlying GSL random generator.
198    ///
199    const gsl_rng* rng(void) const;
200
201    /**
202       Assignment operator
203
204       \since Explicitely declared since yat 0.5
205     */
206    RNG_state& operator=(const RNG_state&);
207
208  private:
209    gsl_rng* rng_;
210
211    void clone(const gsl_rng&);
212  };
213   
214
215  // --------------------- Discrete distribtuions ---------------------
216
217  ///
218  /// @brief Discrete random number distributions.
219  ///
220  /// Abstract base class for discrete random number
221  /// distributions. Given K discrete events with different
222  /// probabilities \f$ P[k] \f$, produce a random value k consistent
223  /// with its probability.
224  ///
225  class Discrete
226  {
227  public:
228    ///
229    /// @brief Constructor
230    ///
231    Discrete(void);
232
233    ///
234    /// @brief The destructor
235    ///
236    virtual ~Discrete(void);
237
238    ///
239    /// @brief Set the seed to \a s.
240    ///
241    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
242    /// default value from the rng's original implementation is used
243    /// (cf. GSL documentation).
244    ///
245    /// \deprecated Provided for backward compatibility with the 0.7
246    /// API. Use RNG::instance()->seed(s) instead.
247    ///
248    void seed(unsigned long s) const YAT_DEPRECATE;
249
250    ///
251    /// @brief Set the seed using the /dev/urandom device.
252    ///
253    /// @return The seed acquired from /dev/urandom.
254    ///
255    /// \deprecated Provided for backward compatibility with the 0.7
256    /// API. Use RNG::instance()->seed_from_devurandom() instead.
257    ///
258    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;
259
260    ///
261    /// @return A random number.
262    ///
263    virtual unsigned long operator()(void) const = 0;
264   
265  protected:
266    /// GSL random gererator
267    RNG* rng_;
268  };
269
270  ///
271  /// @brief General
272  ///
273  class DiscreteGeneral : public Discrete
274  {
275  public:
276    ///
277    /// @brief Constructor
278    ///
279    /// @param hist is a Histogram defining the probability distribution
280    ///
281    explicit DiscreteGeneral(const statistics::Histogram& hist);
282   
283    /**
284       \brief Copy constructor
285
286       \since Explicitely implemented in yat 0.5
287     */
288    DiscreteGeneral(const DiscreteGeneral&);
289
290    ///
291    /// @brief Destructor
292    ///
293    ~DiscreteGeneral(void);
294
295    ///
296    /// The generated number is an integer and proportional to the
297    /// frequency in the corresponding histogram bin. In other words,
298    /// the probability that 0 is returned is proportinal to the size
299    /// of the first bin.
300    ///
301    /// @return A random number.
302    ///
303    unsigned long operator()(void) const;
304
305    /**
306       \brief Assignment operator
307
308       \since Explicitely implemented in yat 0.5
309     */
310    DiscreteGeneral& operator=(const DiscreteGeneral&);
311
312  private:
313    void free(void);
314    void preproc(void);
315
316    gsl_ran_discrete_t* gen_;
317    std::vector<double> p_;
318  };
319
320  /**
321     @brief Discrete uniform distribution
322 
323     Discrete uniform distribution also known as the "equally likely
324     outcomes" distribution. Each outcome, in this case an integer
325     from [0,n-1] , have equal probability to occur.
326     
327     Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
328     k < n \f$ \n
329     Expectation value: \f$ \frac{n-1}{2} \f$ \n
330     Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
331  */
332  class DiscreteUniform 
333    : public Discrete, 
334      public std::unary_function<unsigned long, unsigned long>
335  {
336  public:
337    /**
338       \brief Constructor.
339
340       The generator will generate integers within the range \f$
341       [0,n-1] \f$. If \a n is zero, then the whole range of the
342       underlying RNG will be used \f$ [min,max] \f$. Setting \a n to
343       zero is the preferred way to sample the whole range of the
344       underlying RNG, i.e. not setting \n to RNG.max.
345
346       \throw If \a n is larger than the maximum number the underlying
347       random number generator can return, then a GSL_error exception
348       is thrown.
349    */
350    explicit DiscreteUniform(unsigned long n=0);
351
352    /**
353       \brief Get a random number
354
355       The returned integer is either in the range [RNG.min,RNG.max]
356       or [0,n-1] depending on how the random number generator was
357       created.
358
359       \see DiscreteUniform(const unsigned long n=0)
360    */
361    unsigned long operator()(void) const;
362
363    /**
364       \brief Get a random integer in the range \f$ [0,n-1] \f$.
365
366       All integers in the range [0,n-1] are equally likely. This
367       function should be avoided for sampling the whole range of the
368       underlying RNG.
369
370       \throw GSL_error if \a n is larger than the range of the
371       underlying generator.
372    */
373    unsigned long operator()(unsigned long n) const;
374
375  private:
376    unsigned long range_;
377  };
378
379  /**
380     @brief Poisson Distribution
381 
382     Having a Poisson process (i.e. no memory), number of occurences
383     within a given time window is Poisson distributed. This
384     distribution is the limit of a Binomial distribution when number
385     of attempts is large, and the probability for one attempt to be
386     succesful is small (in such a way that the expected number of
387     succesful attempts is \f$ m \f$.
388     
389     Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
390     k  \f$ \n
391     Expectation value: \f$ m \f$ \n
392     Variance: \f$ m \f$
393  */
394  class Poisson : public Discrete
395  {
396  public:
397    ///
398    /// @brief Constructor
399    ///
400    /// @param m is expectation value
401    ///
402    explicit Poisson(const double m=1);
403
404    ///
405    /// @return A Poisson distributed number.
406    ///
407    unsigned long operator()(void) const;
408
409    ///
410    /// @return A Poisson distributed number with expectation value
411    /// \a m
412    ///
413    /// @note this operator ignores parameters set in Constructor
414    ///
415    unsigned long operator()(const double m) const;
416
417  private:
418    double m_;
419  };
420
421  // --------------------- Continuous distribtuions ---------------------
422
423  ///
424  /// @brief Continuous random number distributions.
425  ///
426  /// Abstract base class for continuous random number distributions.
427  ///
428  class Continuous
429  {
430  public:
431
432    ///
433    /// @brief Constructor
434    ///
435    Continuous(void);
436
437    ///
438    /// @brief The destructor
439    ///
440    virtual ~Continuous(void);
441
442    ///
443    /// @brief Set the seed to \a s.
444    ///
445    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
446    /// default value from the rng's original implementation is used
447    /// (cf. GSL documentation).
448    ///
449    /// \deprecated Provided for backward compatibility with the 0.7
450    /// API. Use RNG::instance()->seed(s) instead.
451    ///
452    void seed(unsigned long s) const YAT_DEPRECATE;
453
454    ///
455    /// @brief Set the seed using the /dev/urandom device.
456    ///
457    /// @return The seed acquired from /dev/urandom.
458    ///
459    /// \deprecated Provided for backward compatibility with the 0.7
460    /// API. Use RNG::instance()->seed_from_devurandom() instead.
461    ///
462    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;
463
464    ///
465    /// @return A random number
466    ///
467    virtual double operator()(void) const = 0;
468
469  protected:
470    /// pointer to GSL random generator
471    RNG* rng_;
472  };
473
474  // ContinuousUniform is declared before ContinuousGeneral to avoid
475  // forward declaration
476  ///
477  /// @brief Uniform distribution
478  ///
479  /// Class for generating a random number from a uniform distribution
480  /// in the range [0,1), i.e. zero is included but not 1.
481  ///
482  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
483  /// Expectation value: 0.5 \n
484  /// Variance: \f$ \frac{1}{12} \f$
485  ///
486  class ContinuousUniform : public Continuous
487  {
488  public:
489    double operator()(void) const;
490  };
491
492  ///
493  /// @brief Generates numbers from a histogram in a continuous manner.
494  ///
495  class ContinuousGeneral : public Continuous
496  {
497  public:
498    ///
499    /// @brief Constructor
500    ///
501    /// @param hist is a Histogram defining the probability distribution
502    ///
503    explicit ContinuousGeneral(const statistics::Histogram& hist);
504
505    ///
506    /// The number is generated in a two step process. First the bin
507    /// in the histogram is randomly selected (see
508    /// DiscreteGeneral). Then a number is generated uniformly from
509    /// the interval defined by the bin.
510    ///
511    /// @return A random number.
512    ///
513    double operator()(void) const;
514
515  private:
516    const DiscreteGeneral discrete_;
517    const statistics::Histogram hist_;
518    ContinuousUniform u_;
519  };
520
521  /**
522     \brief Generator of random numbers from an exponential
523     distribution.
524     
525     The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
526     \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
527     variance \f$ m^2 \f$
528  */
529  class Exponential : public Continuous
530  {
531  public:
532    ///
533    /// @brief Constructor
534    ///
535    /// @param m is the expectation value of the distribution.
536    ///
537    explicit Exponential(const double m=1);
538
539    ///
540    /// @return A random number from exponential distribution.
541    ///
542    double operator()(void) const;
543
544    ///
545    /// @return A random number from exponential distribution, with
546    /// expectation value \a m
547    ///
548    /// @note This operator ignores parameters given in constructor.
549    ///
550    double operator()(const double m) const;
551
552  private:
553    double m_;
554  };
555
556  /**
557     @brief Gaussian distribution
558     
559     Class for generating a random number from a Gaussian distribution
560     between zero and unity. Utilizes the Box-Muller algorithm, which
561     needs two calls to random generator.
562     
563     Distribution function \f$ f(x) =
564     \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
565     \f$ \n
566     Expectation value: \f$ \mu \f$ \n
567     Variance: \f$ \sigma^2 \f$
568  */
569  class Gaussian : public Continuous
570  {
571  public:
572    ///
573    /// @brief Constructor
574    ///
575    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
576    /// @param m is the expectation value \f$ \mu \f$ of the distribution
577    ///
578    explicit Gaussian(const double s=1, const double m=0);
579
580    ///
581    /// @return A random Gaussian number
582    ///
583    double operator()(void) const;
584
585    ///
586    /// @return A random Gaussian number with standard deviation \a s
587    /// and expectation value 0.
588    ///
589    /// @note this operator ignores parameters given in Constructor
590    ///
591    double operator()(const double s) const;
592
593    ///
594    /// @return A random Gaussian number with standard deviation \a s
595    /// and expectation value \a m.
596    ///
597    /// @note this operator ignores parameters given in Constructor
598    ///
599    double operator()(const double s, const double m) const;
600
601  private:
602    double m_;
603    double s_;
604  };
605
606  /**
607     \brief Convenience function to shuffle a range with singleton RNG.
608
609     Wrapper around std::random_shuffle using DiscreteUniform as
610     random generator and thereby using the underlying RNG class,
611     which is singleton.
612
613     RandomAccessIterator must be mutable
614   */
615  template<typename RandomAccessIterator>
616  void random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
617  {
618    typedef RandomAccessIterator rai;
619    BOOST_CONCEPT_ASSERT((boost::Mutable_RandomAccessIterator<rai>));
620    DiscreteUniform rnd;
621    std::random_shuffle(first, last, rnd);
622  }
623
624}}} // of namespace random, yat, and theplu
625
626#endif
Note: See TracBrowser for help on using the repository browser.