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

Last change on this file since 2142 was 2142, checked in by Peter, 12 years ago

clarify iterator type in random_shuffle. refs #489

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