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

Last change on this file since 1184 was 1004, checked in by Peter, 14 years ago

adding a wrapper around std::random_shuffle that uses the RNG class. Changed all calls to random_shuffle to use the added function. Acctually all calls already used the RNG class so essentially this change is only cosmetic, but by providing a function I think it is easier to avoid using multiple random generators.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.2 KB
Line 
1#ifndef _theplu_yat_random_
2#define _theplu_yat_random_
3
4// $Id: random.h 1004 2008-01-23 18:24:53Z peter $
5
6/*
7  Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://trac.thep.lu.se/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 2 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include "yat/statistics/Histogram.h"
28
29#include <gsl/gsl_rng.h>
30#include <gsl/gsl_randist.h>
31
32#include <algorithm>
33#include <string>
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    u_long max(void) const;
101
102    ///
103    /// @brief Returns the smallest number that the random number
104    /// generator can return.
105    ///
106    u_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(u_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    u_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    virtual ~RNG(void);
158
159    static RNG* instance_;
160    gsl_rng* rng_;
161  };
162
163
164  ///
165  /// @brief Class holding state of a random generator
166  ///
167  class RNG_state
168  {
169  public:
170    ///
171    /// @brief Constructor
172    ///
173    RNG_state(const RNG*);
174
175    ///
176    /// @brief Destructor
177    ///
178    ~RNG_state(void);
179
180    ///
181    /// @return const pointer to underlying GSL random generator.
182    ///
183    const gsl_rng* rng(void) const;
184
185  private:
186    gsl_rng* rng_;
187
188  };
189   
190
191  // --------------------- Discrete distribtuions ---------------------
192
193  ///
194  /// @brief Discrete random number distributions.
195  ///
196  /// Abstract base class for discrete random number
197  /// distributions. Given K discrete events with different
198  /// probabilities \f$ P[k] \f$, produce a random value k consistent
199  /// with its probability.
200  ///
201  class Discrete
202  {
203  public:
204    ///
205    /// @brief Constructor
206    ///
207    Discrete(void);
208
209    ///
210    /// @brief The destructor
211    ///
212    virtual ~Discrete(void);
213
214    ///
215    /// @brief Set the seed to \a s.
216    ///
217    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
218    /// default value from the rng's original implementation is used
219    /// (cf. GSL documentation).
220    ///
221    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
222    ///
223    void seed(u_long s) const;
224
225    ///
226    /// @brief Set the seed using the /dev/urandom device.
227    ///
228    /// @return The seed acquired from /dev/urandom.
229    ///
230    /// @see seed, RNG::seed_from_devurandom, RNG::seed
231    ///
232    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
233
234    ///
235    /// @return A random number.
236    ///
237    virtual u_long operator()(void) const = 0;
238   
239  protected:
240    /// GSL random gererator
241    RNG* rng_;
242  };
243
244  ///
245  /// @brief General
246  ///
247  class DiscreteGeneral : public Discrete
248  {
249  public:
250    ///
251    /// @brief Constructor
252    ///
253    /// @param hist is a Histogram defining the probability distribution
254    ///
255    DiscreteGeneral(const statistics::Histogram& hist);
256   
257    ///
258    /// @brief Destructor
259    ///
260    ~DiscreteGeneral(void);
261
262    ///
263    /// The generated number is an integer and proportinal to the
264    /// frequency in the corresponding histogram bin. In other words,
265    /// the probability that 0 is returned is proportinal to the size
266    /// of the first bin.
267    ///
268    /// @return A random number.
269    ///
270    u_long operator()(void) const;
271
272  private:
273     gsl_ran_discrete_t* gen_;
274  };
275
276  ///
277  /// @brief Discrete uniform distribution
278  ///
279  /// Discrete uniform distribution also known as the "equally likely
280  /// outcomes" distribution. Each outcome, in this case an integer
281  /// from [0,n-1] , have equal probability to occur.
282  ///
283  /// Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
284  /// k < n \f$ \n
285  /// Expectation value: \f$ \frac{n-1}{2} \f$ \n
286  /// Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
287  ///
288  class DiscreteUniform : public Discrete
289  {
290  public:
291    /**
292       \brief Constructor.
293
294       The generator will generate integers within the range \f$
295       [0,n-1] \f$. If \a n is zero, then the whole range of the
296       underlying RNG will be used \f$ [min,max] \f$. Setting \a n to
297       zero is the preferred way to sample the whole range of the
298       underlying RNG, i.e. not setting \n to RNG.max.
299
300       \throw If \a n is larger than the maximum number the underlying
301       random number generator can return, then a GSL_error exception
302       is thrown.
303    */
304    DiscreteUniform(u_long n=0);
305
306    /**
307       \brief Get a random number
308
309       The returned integer is either in the range [RNG.min,RNG.max]
310       or [0,n-1] depending on how the random number generator was
311       created.
312
313       \see DiscreteUniform(const u_long n=0)
314    */
315    u_long operator()(void) const;
316
317    /**
318       \brief Get a random integer in the range \f$ [0,n-1] \f$.
319
320       All integers in the range [0,n-1] are equally likely. This
321       function should be avoided for sampling the whole range of the
322       underlying RNG.
323
324       \throw GSL_error if \a n is larger than the range of the
325       underlying generator.
326    */
327    u_long operator()(u_long n) const;
328
329  private:
330    u_long range_;
331  };
332
333  ///
334  /// @brief Poisson Distribution
335  ///
336  /// Having a Poisson process (i.e. no memory), number of occurences
337  /// within a given time window is Poisson distributed. This
338  /// distribution is the limit of a Binomial distribution when number
339  /// of attempts is large, and the probability for one attempt to be
340  /// succesful is small (in such a way that the expected number of
341  /// succesful attempts is \f$ m \f$.
342  ///
343  /// Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
344  /// k  \f$ \n
345  /// Expectation value: \f$ m \f$ \n
346  /// Variance: \f$ m \f$
347  ///
348  class Poisson : public Discrete
349  {
350  public:
351    ///
352    /// @brief Constructor
353    ///
354    /// @param m is expectation value
355    ///
356    Poisson(const double m=1);
357
358    ///
359    /// @return A Poisson distributed number.
360    ///
361    u_long operator()(void) const;
362
363    ///
364    /// @return A Poisson distributed number with expectation value
365    /// \a m
366    ///
367    /// @note this operator ignores parameters set in Constructor
368    ///
369    u_long operator()(const double m) const;
370
371  private:
372    double m_;
373  };
374
375  // --------------------- Continuous distribtuions ---------------------
376
377  ///
378  /// @brief Continuous random number distributions.
379  ///
380  /// Abstract base class for continuous random number distributions.
381  ///
382  class Continuous
383  {
384  public:
385
386    ///
387    /// @brief Constructor
388    ///
389    Continuous(void);
390
391    ///
392    /// @brief The destructor
393    ///
394    virtual ~Continuous(void);
395
396    ///
397    /// @brief Set the seed to \a s.
398    ///
399    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
400    /// default value from the rng's original implementation is used
401    /// (cf. GSL documentation).
402    ///
403    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
404    ///
405    void seed(u_long s) const;
406
407    ///
408    /// @brief Set the seed using the /dev/urandom device.
409    ///
410    /// @return The seed acquired from /dev/urandom.
411    ///
412    /// @see seed, RNG::seed_from_devurandom, RNG::seed
413    ///
414    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
415
416    ///
417    /// @return A random number
418    ///
419    virtual double operator()(void) const = 0;
420
421  protected:
422    /// pointer to GSL random generator
423    RNG* rng_;
424  };
425
426  // ContinuousUniform is declared before ContinuousGeneral to avoid
427  // forward declaration
428  ///
429  /// @brief Uniform distribution
430  ///
431  /// Class for generating a random number from a uniform distribution
432  /// in the range [0,1), i.e. zero is included but not 1.
433  ///
434  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
435  /// Expectation value: 0.5 \n
436  /// Variance: \f$ \frac{1}{12} \f$
437  ///
438  class ContinuousUniform : public Continuous
439  {
440  public:
441    double operator()(void) const;
442  };
443
444  ///
445  /// @brief Generates numbers from a histogram in a continuous manner.
446  ///
447  class ContinuousGeneral : public Continuous
448  {
449  public:
450    ///
451    /// @brief Constructor
452    ///
453    /// @param hist is a Histogram defining the probability distribution
454    ///
455    ContinuousGeneral(const statistics::Histogram& hist);
456
457    ///
458    /// The number is generated in a two step process. First the bin
459    /// in the histogram is randomly selected (see
460    /// DiscreteGeneral). Then a number is generated uniformly from
461    /// the interval defined by the bin.
462    ///
463    /// @return A random number.
464    ///
465    double operator()(void) const;
466
467  private:
468    const DiscreteGeneral discrete_;
469    const statistics::Histogram hist_;
470    ContinuousUniform u_;
471  };
472
473  ///
474  /// @brief Generator of random numbers from an exponential
475  /// distribution.
476  ///
477  /// The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
478  /// \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
479  /// variance \f$ m^2 \f$
480  ///
481  class Exponential : public Continuous
482  {
483  public:
484    ///
485    /// @brief Constructor
486    ///
487    /// @param m is the expectation value of the distribution.
488    ///
489    Exponential(const double m=1);
490
491    ///
492    /// @return A random number from exponential distribution.
493    ///
494    double operator()(void) const;
495
496    ///
497    /// @return A random number from exponential distribution, with
498    /// expectation value \a m
499    ///
500    /// @note This operator ignores parameters given in constructor.
501    ///
502    double operator()(const double m) const;
503
504  private:
505    double m_;
506  };
507
508  ///
509  /// @brief Gaussian distribution
510  ///
511  /// Class for generating a random number from a Gaussian
512  /// distribution between zero and unity. Utilizes the Box-Muller
513  /// algorithm, which needs two calls to random generator.
514  ///
515  /// Distribution function \f$ f(x) =
516  /// \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
517  /// \f$ \n
518  /// Expectation value: \f$ \mu \f$ \n
519  /// Variance: \f$ \sigma^2 \f$
520  ///
521  class Gaussian : public Continuous
522  {
523  public:
524    ///
525    /// @brief Constructor
526    ///
527    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
528    /// @param m is the expectation value \f$ \mu \f$ of the distribution
529    ///
530    Gaussian(const double s=1, const double m=0);
531
532    ///
533    /// @return A random Gaussian number
534    ///
535    double operator()(void) const;
536
537    ///
538    /// @return A random Gaussian number with standard deviation \a s
539    /// and expectation value 0.
540    ///
541    /// @note this operator ignores parameters given in Constructor
542    ///
543    double operator()(const double s) const;
544
545    ///
546    /// @return A random Gaussian number with standard deviation \a s
547    /// and expectation value \a m.
548    ///
549    /// @note this operator ignores parameters given in Constructor
550    ///
551    double operator()(const double s, const double m) const;
552
553  private:
554    double m_;
555    double s_;
556  };
557
558  /**
559     \brief Convenience function to shuffle a range with singleton RNG.
560
561     Wrapper around std::random_shuffle using DiscreteUniform as
562     random generator and thereby using the underlying RNG class,
563     which is singleton.
564   */
565  template<typename T>
566  void random_shuffle(T first, T last)
567  {
568    DiscreteUniform rnd;
569    std::random_shuffle(first, last, rnd);
570  }
571
572}}} // of namespace random, yat, and theplu
573
574#endif
Note: See TracBrowser for help on using the repository browser.