Logo ROOT   6.10/09
Reference Guide
GSLRndmEngines.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: L. Moneta, A. Zsenei 08/2005
3 
4  /**********************************************************************
5  * *
6  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
7  * *
8  * This library is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License *
10  * as published by the Free Software Foundation; either version 2 *
11  * of the License, or (at your option) any later version. *
12  * *
13  * This library is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
16  * General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this library (see file COPYING); if not, write *
20  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21  * 330, Boston, MA 02111-1307 USA, or contact the author. *
22  * *
23  **********************************************************************/
24 
25 // Header file for class GSLRandom
26 //
27 // Created by: moneta at Sun Nov 21 16:26:03 2004
28 //
29 // Last update: Sun Nov 21 16:26:03 2004
30 //
31 
32 
33 
34 // need to be included later
35 #include <time.h>
36 #include <cassert>
37 
38 #include "gsl/gsl_rng.h"
39 #include "gsl/gsl_randist.h"
40 
41 
42 #include "Math/GSLRndmEngines.h"
43 #include "GSLRngWrapper.h"
44 // for wrapping in GSL ROOT engines
45 #include "GSLRngROOTWrapper.h"
46 
47 extern double gsl_ran_gaussian_acr( const gsl_rng * r, const double sigma);
48 
49 //#include <iostream>
50 
51 namespace ROOT {
52 namespace Math {
53 
54 
55 
56 
57 
58  // default constructor (need to call set type later)
60  fRng(nullptr),
61  fCurTime(0)
62  { }
63 
64  // constructor from external rng
65  // internal generator will be managed or not depending on
66  // how the GSLRngWrapper is created
68  fRng(new GSLRngWrapper(*rng) ),
69  fCurTime(0)
70  {}
71 
72  // copy constructor
74  fRng(new GSLRngWrapper(*eng.fRng) ),
75  fCurTime(0)
76  {}
77 
79  // destructor : call terminate if not yet called
80  if (fRng) Terminate();
81  }
82 
83  // assignment operator
85  if (this == &eng) return *this;
86  if (fRng)
87  *fRng = *eng.fRng;
88  else
89  fRng = new GSLRngWrapper(*eng.fRng);
90  fCurTime = eng.fCurTime;
91  return *this;
92  }
93 
94 
96  // initialize the generator by allocating the GSL object
97  // if type was not passed create with default generator
98  if (!fRng) fRng = new GSLRngWrapper();
99  fRng->Allocate();
100  }
101 
103  // terminate the generator by freeing the GSL object
104  if (!fRng) return;
105  fRng->Free();
106  delete fRng;
107  fRng = 0;
108  }
109 
110 
112  // generate random between 0 and 1.
113  // 0 is excluded
114  return gsl_rng_uniform_pos(fRng->Rng() );
115  }
116 
117 
118  unsigned long GSLRandomEngine::RndmInt(unsigned long max) const {
119  // generate a random integer number between 0 and MAX
120  return gsl_rng_uniform_int( fRng->Rng(), max );
121  }
122 
123  unsigned long GSLRandomEngine::MinInt() const {
124  // return minimum integer value used in RndmInt
125  return gsl_rng_min( fRng->Rng() );
126  }
127 
128  unsigned long GSLRandomEngine::MaxInt() const {
129  // return maximum integr value used in RndmInt
130  return gsl_rng_max( fRng->Rng() );
131  }
132 
133  void GSLRandomEngine::RandomArray(double * begin, double * end ) const {
134  // generate array of randoms betweeen 0 and 1. 0 is excluded
135  // specialization for double * (to be faster)
136  for ( double * itr = begin; itr != end; ++itr ) {
137  *itr = gsl_rng_uniform_pos(fRng->Rng() );
138  }
139  }
140 
141  void GSLRandomEngine::SetSeed(unsigned int seed) const {
142  // set the seed, if = 0then the seed is set randomly using an std::rand()
143  // seeded with the current time. Be carefuk in case the current time is
144  // the same in consecutive calls
145  if (seed == 0) {
146  // use like in root (use time)
147  time_t curtime;
148  time(&curtime);
149  unsigned int ct = static_cast<unsigned int>(curtime);
150  if (ct != fCurTime) {
151  fCurTime = ct;
152  // set the seed for rand
153  srand(ct);
154  }
155  seed = rand();
156  }
157 
158  assert(fRng);
159  gsl_rng_set(fRng->Rng(), seed );
160  }
161 
162  std::string GSLRandomEngine::Name() const {
163  //////////////////////////////////////////////////////////////////////////
164 
165  assert ( fRng != 0);
166  assert ( fRng->Rng() != 0 );
167  return std::string( gsl_rng_name( fRng->Rng() ) );
168  }
169 
170  unsigned int GSLRandomEngine::Size() const {
171  //////////////////////////////////////////////////////////////////////////
172 
173  assert (fRng != 0);
174  return gsl_rng_size( fRng->Rng() );
175  }
176 
177 
178  // Random distributions
179 
180  double GSLRandomEngine::GaussianZig(double sigma) const
181  {
182  // Gaussian distribution. Use fast ziggurat algorithm implemented since GSL 1.8
183  return gsl_ran_gaussian_ziggurat( fRng->Rng(), sigma);
184  }
185 
186  double GSLRandomEngine::Gaussian(double sigma) const
187  {
188  // Gaussian distribution. Use default Box-Muller method
189  return gsl_ran_gaussian( fRng->Rng(), sigma);
190  }
191 
193  {
194  // Gaussian distribution. Use ratio method
195  return gsl_ran_gaussian_ratio_method( fRng->Rng(), sigma);
196  }
197 
198 
199  double GSLRandomEngine::GaussianTail(double a , double sigma) const
200  {
201  // Gaussian Tail distribution: eeturn values larger than a distributed
202  // according to the gaussian
203  return gsl_ran_gaussian_tail( fRng->Rng(), a, sigma);
204  }
205 
206  void GSLRandomEngine::Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const
207  {
208  // Gaussian Bivariate distribution, with correlation coefficient rho
209  gsl_ran_bivariate_gaussian( fRng->Rng(), sigmaX, sigmaY, rho, &x, &y);
210  }
211 
212  double GSLRandomEngine::Exponential(double mu) const
213  {
214  // Exponential distribution
215  return gsl_ran_exponential( fRng->Rng(), mu);
216  }
217 
218  double GSLRandomEngine::Cauchy(double a) const
219  {
220  // Cauchy distribution
221  return gsl_ran_cauchy( fRng->Rng(), a);
222  }
223 
224  double GSLRandomEngine::Landau() const
225  {
226  // Landau distribution
227  return gsl_ran_landau( fRng->Rng());
228  }
229 
230  double GSLRandomEngine::Beta(double a, double b) const
231  {
232  // Beta distribution
233  return gsl_ran_beta( fRng->Rng(), a, b);
234  }
235 
236  double GSLRandomEngine::Gamma(double a, double b) const
237  {
238  // Gamma distribution
239  return gsl_ran_gamma( fRng->Rng(), a, b);
240  }
241 
242  double GSLRandomEngine::LogNormal(double zeta, double sigma) const
243  {
244  // Log normal distribution
245  return gsl_ran_lognormal( fRng->Rng(), zeta, sigma);
246  }
247 
248  double GSLRandomEngine::ChiSquare(double nu) const
249  {
250  // Chi square distribution
251  return gsl_ran_chisq( fRng->Rng(), nu);
252  }
253 
254 
255  double GSLRandomEngine::FDist(double nu1, double nu2) const
256  {
257  // F distribution
258  return gsl_ran_fdist( fRng->Rng(), nu1, nu2);
259  }
260 
261  double GSLRandomEngine::tDist(double nu) const
262  {
263  // t distribution
264  return gsl_ran_tdist( fRng->Rng(), nu);
265  }
266 
267  double GSLRandomEngine::Rayleigh(double sigma) const
268  {
269  // Rayleigh distribution
270  return gsl_ran_rayleigh( fRng->Rng(), sigma);
271  }
272 
273  double GSLRandomEngine::Logistic(double a) const
274  {
275  // Logistic distribution
276  return gsl_ran_logistic( fRng->Rng(), a);
277  }
278 
279  double GSLRandomEngine::Pareto(double a, double b) const
280  {
281  // Pareto distribution
282  return gsl_ran_pareto( fRng->Rng(), a, b);
283  }
284 
285  void GSLRandomEngine::Dir2D(double &x, double &y) const
286  {
287  // generate random numbers in a 2D circle of radious 1
288  gsl_ran_dir_2d( fRng->Rng(), &x, &y);
289  }
290 
291  void GSLRandomEngine::Dir3D(double &x, double &y, double &z) const
292  {
293  // generate random numbers in a 3D sphere of radious 1
294  gsl_ran_dir_3d( fRng->Rng(), &x, &y, &z);
295  }
296 
297  unsigned int GSLRandomEngine::Poisson(double mu) const
298  {
299  // Poisson distribution
300  return gsl_ran_poisson( fRng->Rng(), mu);
301  }
302 
303  unsigned int GSLRandomEngine::Binomial(double p, unsigned int n) const
304  {
305  // Binomial distribution
306  return gsl_ran_binomial( fRng->Rng(), p, n);
307  }
308 
309  unsigned int GSLRandomEngine::NegativeBinomial(double p, double n) const
310  {
311  // Negative Binomial distribution
312  return gsl_ran_negative_binomial( fRng->Rng(), p, n);
313  }
314 
315 
316  std::vector<unsigned int> GSLRandomEngine::Multinomial( unsigned int ntot, const std::vector<double> & p ) const
317  {
318  // Multinomial distribution return vector of integers which sum is ntot
319  std::vector<unsigned int> ival( p.size());
320  gsl_ran_multinomial( fRng->Rng(), p.size(), ntot, &p.front(), &ival[0]);
321  return ival;
322  }
323 
324 
325 
326  //----------------------------------------------------
327  // generators
328  //----------------------------------------------------
329 
330  /////////////////////////////////////////////////////////////////////////////
331 
333  {
334  SetType(new GSLRngWrapper(gsl_rng_mt19937));
335  Initialize();
336  }
337 
338 
339  // old ranlux - equivalent to TRandom1
341  {
342  SetType(new GSLRngWrapper(gsl_rng_ranlux) );
343  Initialize();
344  }
345 
346  // second generation of Ranlux (single precision version - luxury 1)
348  {
349  SetType(new GSLRngWrapper(gsl_rng_ranlxs1) );
350  Initialize();
351  }
352 
353  // second generation of Ranlux (single precision version - luxury 2)
355  {
356  SetType(new GSLRngWrapper(gsl_rng_ranlxs2) );
357  Initialize();
358  }
359 
360  // double precision version - luxury 1
362  {
363  SetType(new GSLRngWrapper(gsl_rng_ranlxd1) );
364  Initialize();
365  }
366 
367  // double precision version - luxury 2
369  {
370  SetType(new GSLRngWrapper(gsl_rng_ranlxd2) );
371  Initialize();
372  }
373 
374  /////////////////////////////////////////////////////////////////////////////
375 
377  {
378  SetType(new GSLRngWrapper(gsl_rng_taus2) );
379  Initialize();
380  }
381 
382  /////////////////////////////////////////////////////////////////////////////
383 
385  {
386  SetType(new GSLRngWrapper(gsl_rng_gfsr4) );
387  Initialize();
388  }
389 
390  /////////////////////////////////////////////////////////////////////////////
391 
393  {
394  SetType(new GSLRngWrapper(gsl_rng_cmrg) );
395  Initialize();
396  }
397 
398  /////////////////////////////////////////////////////////////////////////////
399 
401  {
402  SetType(new GSLRngWrapper(gsl_rng_mrg) );
403  Initialize();
404  }
405 
406 
407  /////////////////////////////////////////////////////////////////////////////
408 
410  {
411  SetType(new GSLRngWrapper(gsl_rng_rand) );
412  Initialize();
413  }
414 
415  /////////////////////////////////////////////////////////////////////////////
416 
418  {
419  SetType(new GSLRngWrapper(gsl_rng_ranmar) );
420  Initialize();
421  }
422 
423  /////////////////////////////////////////////////////////////////////////////
424 
426  {
427  SetType(new GSLRngWrapper(gsl_rng_minstd) );
428  Initialize();
429  }
430 
431 
432  // for extra engines based on ROOT
434  {
436  Initialize();
437  }
439  // we need to explicitly delete the ROOT wrapper class
440  GSLMixMaxWrapper::Free(Engine()->Rng()->state);
441  }
442 
443 } // namespace Math
444 } // namespace ROOT
445 
446 
447 
double Beta(double a, double b) const
Beta distribution.
double Rayleigh(double sigma) const
Rayleigh distribution.
GSLRandomEngine & operator=(const GSLRandomEngine &eng)
Assignment operator : make a deep copy of the contained GSL generator.
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double Gamma(double a, double b) const
Gamma distribution.
double Logistic(double a) const
Logistic distribution.
void SetType(GSLRngWrapper *r)
internal method used by the derived class to set the type of generators
double tDist(double nu) const
t student distribution
TArc * a
Definition: textangle.C:12
void RandomArray(Iterator begin, Iterator end) const
Generate an array of random numbers.
double FDist(double nu1, double nu2) const
F distrbution.
unsigned int Poisson(double mu) const
Poisson distribution.
double Pareto(double a, double b) const
Pareto distribution.
Double_t x[n]
Definition: legend1.C:17
unsigned int NegativeBinomial(double p, double n) const
Negative Binomial distribution.
void Dir3D(double &x, double &y, double &z) const
generate random numbers in a 3D sphere of radious 1
void Terminate()
delete pointer to contained rng
double ChiSquare(double nu) const
Chi square distribution.
GSLRandomEngine Base class for all GSL random engines, normally user instantiate the derived classes ...
GSLRngWrapper class to wrap gsl_rng structure.
Definition: GSLRngWrapper.h:25
unsigned int Size() const
return the state size of generator
double GaussianZig(double sigma) const
Gaussian distribution - Ziggurat method.
const Double_t sigma
double Landau() const
Landau distribution.
double Exponential(double mu) const
Exponential distribution.
std::vector< unsigned int > Multinomial(unsigned int ntot, const std::vector< double > &p) const
Multinomial distribution.
unsigned long MinInt() const
return the minimum integer a generator can handle typically this value is 0
virtual ~GSLRandomEngine()
call Terminate()
TRandom2 r(17)
double LogNormal(double zeta, double sigma) const
Log Normal distribution.
unsigned long RndmInt(unsigned long max) const
Generate an integer number between [0,max-1] (including 0 and max-1) if max is larger than available ...
GSLRngWrapper * Engine()
internal method to return the engine Used by class like GSLMCIntegrator to set the engine ...
void Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const
Bivariate Gaussian distribution with correlation.
unsigned int Binomial(double p, unsigned int n) const
Binomial distribution.
double GaussianTail(double a, double sigma) const
Gaussian Tail distribution.
double operator()() const
Generate a random number between ]0,1] 0 is excluded and 1 is included.
void SetSeed(unsigned int seed) const
set the random generator seed
Double_t y[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
std::string Name() const
return name of generator
double Cauchy(double a) const
Cauchy distribution.
unsigned long MaxInt() const
return the maximum integer +1 a generator can handle
double GaussianRatio(double sigma) const
Gaussian distribution - Ratio method.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
void Initialize()
initialize the generator If no rng is present the default one based on Mersenne and Twister is create...
double Gaussian(double sigma) const
Gaussian distribution - default method is Box-Muller (polar method)
const Int_t n
Definition: legend1.C:16
const gsl_rng_type * gsl_rng_mixmax
GSLRandomEngine()
default constructor.
double gsl_ran_gaussian_acr(const gsl_rng *r, const double sigma)
void Dir2D(double &x, double &y) const
generate random numbers in a 2D circle of radious 1