Logo ROOT  
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 <ctime>
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 namespace ROOT {
50 namespace Math {
51 
52 
53 
54 
55 
56  // default constructor (need to call set type later)
58  fRng(nullptr),
59  fCurTime(0)
60  { }
61 
62  // constructor from external rng
63  // internal generator will be managed or not depending on
64  // how the GSLRngWrapper is created
66  fRng(new GSLRngWrapper(*rng) ),
67  fCurTime(0)
68  {}
69 
70  // copy constructor
72  fRng(new GSLRngWrapper(*eng.fRng) ),
73  fCurTime(0)
74  {}
75 
77  // destructor : call terminate if not yet called
78  if (fRng) Terminate();
79  }
80 
81  // assignment operator
83  if (this == &eng) return *this;
84  if (fRng)
85  *fRng = *eng.fRng;
86  else
87  fRng = new GSLRngWrapper(*eng.fRng);
88  fCurTime = eng.fCurTime;
89  return *this;
90  }
91 
92 
94  // initialize the generator by allocating the GSL object
95  // if type was not passed create with default generator
96  if (!fRng) fRng = new GSLRngWrapper();
97  fRng->Allocate();
98  }
99 
101  // terminate the generator by freeing the GSL object
102  if (!fRng) return;
103  fRng->Free();
104  delete fRng;
105  fRng = 0;
106  }
107 
108 
110  // generate random between 0 and 1.
111  // 0 is excluded
112  return gsl_rng_uniform_pos(fRng->Rng() );
113  }
114 
115 
116  unsigned long GSLRandomEngine::RndmInt(unsigned long max) const {
117  // generate a random integer number between 0 and MAX
118  return gsl_rng_uniform_int( fRng->Rng(), max );
119  }
120 
121  unsigned long GSLRandomEngine::MinInt() const {
122  // return minimum integer value used in RndmInt
123  return gsl_rng_min( fRng->Rng() );
124  }
125 
126  unsigned long GSLRandomEngine::MaxInt() const {
127  // return maximum integr value used in RndmInt
128  return gsl_rng_max( fRng->Rng() );
129  }
130 
131  void GSLRandomEngine::RandomArray(double * begin, double * end ) const {
132  // generate array of randoms betweeen 0 and 1. 0 is excluded
133  // specialization for double * (to be faster)
134  for ( double * itr = begin; itr != end; ++itr ) {
135  *itr = gsl_rng_uniform_pos(fRng->Rng() );
136  }
137  }
138 
139  void GSLRandomEngine::SetSeed(unsigned int seed) const {
140  // set the seed, if = 0then the seed is set randomly using an std::rand()
141  // seeded with the current time. Be carefuk in case the current time is
142  // the same in consecutive calls
143  if (seed == 0) {
144  // use like in root (use time)
145  time_t curtime;
146  time(&curtime);
147  unsigned int ct = static_cast<unsigned int>(curtime);
148  if (ct != fCurTime) {
149  fCurTime = ct;
150  // set the seed for rand
151  srand(ct);
152  }
153  seed = rand();
154  }
155 
156  assert(fRng);
157  gsl_rng_set(fRng->Rng(), seed );
158  }
159 
160  std::string GSLRandomEngine::Name() const {
161  //////////////////////////////////////////////////////////////////////////
162 
163  assert ( fRng != 0);
164  assert ( fRng->Rng() != 0 );
165  return std::string( gsl_rng_name( fRng->Rng() ) );
166  }
167 
168  unsigned int GSLRandomEngine::Size() const {
169  //////////////////////////////////////////////////////////////////////////
170 
171  assert (fRng != 0);
172  return gsl_rng_size( fRng->Rng() );
173  }
174 
175 
176  // Random distributions
177 
178  double GSLRandomEngine::GaussianZig(double sigma) const
179  {
180  // Gaussian distribution. Use fast ziggurat algorithm implemented since GSL 1.8
181  return gsl_ran_gaussian_ziggurat( fRng->Rng(), sigma);
182  }
183 
184  double GSLRandomEngine::Gaussian(double sigma) const
185  {
186  // Gaussian distribution. Use default Box-Muller method
187  return gsl_ran_gaussian( fRng->Rng(), sigma);
188  }
189 
191  {
192  // Gaussian distribution. Use ratio method
193  return gsl_ran_gaussian_ratio_method( fRng->Rng(), sigma);
194  }
195 
196 
197  double GSLRandomEngine::GaussianTail(double a , double sigma) const
198  {
199  // Gaussian Tail distribution: eeturn values larger than a distributed
200  // according to the gaussian
201  return gsl_ran_gaussian_tail( fRng->Rng(), a, sigma);
202  }
203 
204  void GSLRandomEngine::Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const
205  {
206  // Gaussian Bivariate distribution, with correlation coefficient rho
207  gsl_ran_bivariate_gaussian( fRng->Rng(), sigmaX, sigmaY, rho, &x, &y);
208  }
209 
210  double GSLRandomEngine::Exponential(double mu) const
211  {
212  // Exponential distribution
213  return gsl_ran_exponential( fRng->Rng(), mu);
214  }
215 
216  double GSLRandomEngine::Cauchy(double a) const
217  {
218  // Cauchy distribution
219  return gsl_ran_cauchy( fRng->Rng(), a);
220  }
221 
222  double GSLRandomEngine::Landau() const
223  {
224  // Landau distribution
225  return gsl_ran_landau( fRng->Rng());
226  }
227 
228  double GSLRandomEngine::Beta(double a, double b) const
229  {
230  // Beta distribution
231  return gsl_ran_beta( fRng->Rng(), a, b);
232  }
233 
234  double GSLRandomEngine::Gamma(double a, double b) const
235  {
236  // Gamma distribution
237  return gsl_ran_gamma( fRng->Rng(), a, b);
238  }
239 
240  double GSLRandomEngine::LogNormal(double zeta, double sigma) const
241  {
242  // Log normal distribution
243  return gsl_ran_lognormal( fRng->Rng(), zeta, sigma);
244  }
245 
246  double GSLRandomEngine::ChiSquare(double nu) const
247  {
248  // Chi square distribution
249  return gsl_ran_chisq( fRng->Rng(), nu);
250  }
251 
252 
253  double GSLRandomEngine::FDist(double nu1, double nu2) const
254  {
255  // F distribution
256  return gsl_ran_fdist( fRng->Rng(), nu1, nu2);
257  }
258 
259  double GSLRandomEngine::tDist(double nu) const
260  {
261  // t distribution
262  return gsl_ran_tdist( fRng->Rng(), nu);
263  }
264 
265  double GSLRandomEngine::Rayleigh(double sigma) const
266  {
267  // Rayleigh distribution
268  return gsl_ran_rayleigh( fRng->Rng(), sigma);
269  }
270 
271  double GSLRandomEngine::Logistic(double a) const
272  {
273  // Logistic distribution
274  return gsl_ran_logistic( fRng->Rng(), a);
275  }
276 
277  double GSLRandomEngine::Pareto(double a, double b) const
278  {
279  // Pareto distribution
280  return gsl_ran_pareto( fRng->Rng(), a, b);
281  }
282 
283  void GSLRandomEngine::Dir2D(double &x, double &y) const
284  {
285  // generate random numbers in a 2D circle of radious 1
286  gsl_ran_dir_2d( fRng->Rng(), &x, &y);
287  }
288 
289  void GSLRandomEngine::Dir3D(double &x, double &y, double &z) const
290  {
291  // generate random numbers in a 3D sphere of radious 1
292  gsl_ran_dir_3d( fRng->Rng(), &x, &y, &z);
293  }
294 
295  unsigned int GSLRandomEngine::Poisson(double mu) const
296  {
297  // Poisson distribution
298  return gsl_ran_poisson( fRng->Rng(), mu);
299  }
300 
301  unsigned int GSLRandomEngine::Binomial(double p, unsigned int n) const
302  {
303  // Binomial distribution
304  return gsl_ran_binomial( fRng->Rng(), p, n);
305  }
306 
307  unsigned int GSLRandomEngine::NegativeBinomial(double p, double n) const
308  {
309  // Negative Binomial distribution
310  return gsl_ran_negative_binomial( fRng->Rng(), p, n);
311  }
312 
313 
314  std::vector<unsigned int> GSLRandomEngine::Multinomial( unsigned int ntot, const std::vector<double> & p ) const
315  {
316  // Multinomial distribution return vector of integers which sum is ntot
317  std::vector<unsigned int> ival( p.size());
318  gsl_ran_multinomial( fRng->Rng(), p.size(), ntot, &p.front(), &ival[0]);
319  return ival;
320  }
321 
322 
323 
324  //----------------------------------------------------
325  // generators
326  //----------------------------------------------------
327 
328  /////////////////////////////////////////////////////////////////////////////
329 
331  {
332  SetType(new GSLRngWrapper(gsl_rng_mt19937));
333  Initialize();
334  }
335 
336 
337  // old ranlux - equivalent to TRandom1
339  {
340  SetType(new GSLRngWrapper(gsl_rng_ranlux) );
341  Initialize();
342  }
343 
344  // second generation of Ranlux (single precision version - luxury 1)
346  {
347  SetType(new GSLRngWrapper(gsl_rng_ranlxs1) );
348  Initialize();
349  }
350 
351  // second generation of Ranlux (single precision version - luxury 2)
353  {
354  SetType(new GSLRngWrapper(gsl_rng_ranlxs2) );
355  Initialize();
356  }
357 
358  // double precision version - luxury 1
360  {
361  SetType(new GSLRngWrapper(gsl_rng_ranlxd1) );
362  Initialize();
363  }
364 
365  // double precision version - luxury 2
367  {
368  SetType(new GSLRngWrapper(gsl_rng_ranlxd2) );
369  Initialize();
370  }
371 
372  /////////////////////////////////////////////////////////////////////////////
373 
375  {
376  SetType(new GSLRngWrapper(gsl_rng_taus2) );
377  Initialize();
378  }
379 
380  /////////////////////////////////////////////////////////////////////////////
381 
383  {
384  SetType(new GSLRngWrapper(gsl_rng_gfsr4) );
385  Initialize();
386  }
387 
388  /////////////////////////////////////////////////////////////////////////////
389 
391  {
392  SetType(new GSLRngWrapper(gsl_rng_cmrg) );
393  Initialize();
394  }
395 
396  /////////////////////////////////////////////////////////////////////////////
397 
399  {
400  SetType(new GSLRngWrapper(gsl_rng_mrg) );
401  Initialize();
402  }
403 
404 
405  /////////////////////////////////////////////////////////////////////////////
406 
408  {
409  SetType(new GSLRngWrapper(gsl_rng_rand) );
410  Initialize();
411  }
412 
413  /////////////////////////////////////////////////////////////////////////////
414 
416  {
417  SetType(new GSLRngWrapper(gsl_rng_ranmar) );
418  Initialize();
419  }
420 
421  /////////////////////////////////////////////////////////////////////////////
422 
424  {
425  SetType(new GSLRngWrapper(gsl_rng_minstd) );
426  Initialize();
427  }
428 
429 
430  // for extra engines based on ROOT
432  {
434  Initialize(); // this creates the gsl_rng structure
435  // no real need to call CreateEngine since the underlined MIXMAX engine is created
436  // by calling GSLMixMaxWrapper::Seed(gsl_default_seed) that is called
437  // when gsl_rng is allocated (in Initialize)
439  }
441  // we need to explicitly delete the ROOT wrapper class
443  }
444 
445 } // namespace Math
446 } // namespace ROOT
447 
448 
449 
ROOT::Math::GSLRandomEngine::GSLRandomEngine
GSLRandomEngine()
default constructor.
Definition: GSLRndmEngines.cxx:57
ROOT::Math::GSLRandomEngine::fCurTime
unsigned int fCurTime
Definition: GSLRndmEngines.h:322
ROOT::Math::GSLRandomEngine::RndmInt
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 ...
Definition: GSLRndmEngines.cxx:116
ROOT::Math::GSLRngRanLuxS1::GSLRngRanLuxS1
GSLRngRanLuxS1()
Definition: GSLRndmEngines.cxx:345
n
const Int_t n
Definition: legend1.C:16
ROOT::Math::GSLRngRanLux::GSLRngRanLux
GSLRngRanLux()
Definition: GSLRndmEngines.cxx:338
ROOT::Math::GSLRngMRG::GSLRngMRG
GSLRngMRG()
Definition: GSLRndmEngines.cxx:398
ROOT::Math::GSLRandomEngine::Name
std::string Name() const
return name of generator
Definition: GSLRndmEngines.cxx:160
ROOT::Math::GSLRandomEngine::SetSeed
void SetSeed(unsigned int seed) const
set the random generator seed
Definition: GSLRndmEngines.cxx:139
ROOT::Math::GSLRngMixMax::GSLRngMixMax
GSLRngMixMax()
Definition: GSLRndmEngines.cxx:431
ROOT::Math::GSLRngGFSR4::GSLRngGFSR4
GSLRngGFSR4()
Definition: GSLRndmEngines.cxx:382
ROOT::Math::GSLRngROOTWrapper::FreeEngine
static void FreeEngine(gsl_rng *r)
Definition: GSLRngROOTWrapper.h:71
ROOT::Math::GSLRandomEngine::Binomial
unsigned int Binomial(double p, unsigned int n) const
Binomial distribution.
Definition: GSLRndmEngines.cxx:301
ROOT::Math::GSLRandomEngine::Pareto
double Pareto(double a, double b) const
Pareto distribution.
Definition: GSLRndmEngines.cxx:277
ROOT::Math::GSLRandomEngine::Gamma
double Gamma(double a, double b) const
Gamma distribution.
Definition: GSLRndmEngines.cxx:234
r
ROOT::R::TRInterface & r
Definition: Object.C:4
ROOT::Math::GSLRngWrapper
GSLRngWrapper class to wrap gsl_rng structure.
Definition: GSLRngWrapper.h:25
ROOT::Math::GSLRandomEngine::NegativeBinomial
unsigned int NegativeBinomial(double p, double n) const
Negative Binomial distribution.
Definition: GSLRndmEngines.cxx:307
ROOT::Math::GSLRandomEngine::Multinomial
std::vector< unsigned int > Multinomial(unsigned int ntot, const std::vector< double > &p) const
Multinomial distribution.
Definition: GSLRndmEngines.cxx:314
ROOT::Math::GSLRngWrapper::Allocate
void Allocate()
Definition: GSLRngWrapper.h:96
ROOT::Math::GSLRandomEngine::Poisson
unsigned int Poisson(double mu) const
Poisson distribution.
Definition: GSLRndmEngines.cxx:295
ROOT::Math::GSLRngMinStd::GSLRngMinStd
GSLRngMinStd()
Definition: GSLRndmEngines.cxx:423
gsl_rng_mixmax
const gsl_rng_type * gsl_rng_mixmax
Definition: GSLRngROOTWrapper.h:106
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Math::GSLRandomEngine::tDist
double tDist(double nu) const
t student distribution
Definition: GSLRndmEngines.cxx:259
ROOT::Math::GSLRngRand::GSLRngRand
GSLRngRand()
Definition: GSLRndmEngines.cxx:407
ROOT::Math::GSLRandomEngine::RandomArray
void RandomArray(Iterator begin, Iterator end) const
Generate an array of random numbers.
Definition: GSLRndmEngines.h:142
ROOT::Math::GSLRngWrapper::Rng
gsl_rng * Rng()
Definition: GSLRngWrapper.h:125
ROOT::Math::GSLRandomEngine::Dir2D
void Dir2D(double &x, double &y) const
generate random numbers in a 2D circle of radious 1
Definition: GSLRndmEngines.cxx:283
GSLRngROOTWrapper.h
ROOT::Math::GSLRandomEngine::Logistic
double Logistic(double a) const
Logistic distribution.
Definition: GSLRndmEngines.cxx:271
ROOT::Math::GSLRandomEngine::Gaussian2D
void Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const
Bivariate Gaussian distribution with correlation.
Definition: GSLRndmEngines.cxx:204
ROOT::Math::GSLRngRanMar::GSLRngRanMar
GSLRngRanMar()
Definition: GSLRndmEngines.cxx:415
b
#define b(i)
Definition: RSha256.hxx:100
ROOT::Math::GSLRandomEngine::Size
unsigned int Size() const
return the state size of generator
Definition: GSLRndmEngines.cxx:168
ROOT::Math::GSLRandomEngine::Beta
double Beta(double a, double b) const
Beta distribution.
Definition: GSLRndmEngines.cxx:228
gsl_ran_gaussian_acr
double gsl_ran_gaussian_acr(const gsl_rng *r, const double sigma)
ROOT::Math::GSLRandomEngine::FDist
double FDist(double nu1, double nu2) const
F distrbution.
Definition: GSLRndmEngines.cxx:253
ROOT::Math::GSLRngRanLuxS2::GSLRngRanLuxS2
GSLRngRanLuxS2()
Definition: GSLRndmEngines.cxx:352
ROOT::Math::GSLRngRanLuxD2::GSLRngRanLuxD2
GSLRngRanLuxD2()
Definition: GSLRndmEngines.cxx:366
ROOT::Math::GSLRandomEngine::SetType
void SetType(GSLRngWrapper *r)
internal method used by the derived class to set the type of generators
Definition: GSLRndmEngines.h:309
ROOT::Math::GSLRandomEngine::Terminate
void Terminate()
delete pointer to contained rng
Definition: GSLRndmEngines.cxx:100
ROOT::Math::GSLRandomEngine::Dir3D
void Dir3D(double &x, double &y, double &z) const
generate random numbers in a 3D sphere of radious 1
Definition: GSLRndmEngines.cxx:289
ROOT::Math::GSLRngRanLuxD1::GSLRngRanLuxD1
GSLRngRanLuxD1()
Definition: GSLRndmEngines.cxx:359
ROOT::Math::GSLRandomEngine::Landau
double Landau() const
Landau distribution.
Definition: GSLRndmEngines.cxx:222
ROOT::Math::GSLRandomEngine::Gaussian
double Gaussian(double sigma) const
Gaussian distribution - default method is Box-Muller (polar method)
Definition: GSLRndmEngines.cxx:184
ROOT::Math::GSLRandomEngine::Cauchy
double Cauchy(double a) const
Cauchy distribution.
Definition: GSLRndmEngines.cxx:216
a
auto * a
Definition: textangle.C:12
ROOT::Math::GSLRandomEngine::operator()
double operator()() const
Generate a random number between ]0,1] 0 is excluded and 1 is included.
Definition: GSLRndmEngines.cxx:109
y
Double_t y[n]
Definition: legend1.C:17
ROOT::Math::GSLRngWrapper::Free
void Free()
Definition: GSLRngWrapper.h:103
ROOT::Math::GSLRandomEngine::ChiSquare
double ChiSquare(double nu) const
Chi square distribution.
Definition: GSLRndmEngines.cxx:246
ROOT::Math::GSLRandomEngine::Exponential
double Exponential(double mu) const
Exponential distribution.
Definition: GSLRndmEngines.cxx:210
sigma
const Double_t sigma
Definition: h1analysisProxy.h:11
ROOT::Math::GSLRandomEngine::operator=
GSLRandomEngine & operator=(const GSLRandomEngine &eng)
Assignment operator : make a deep copy of the contained GSL generator.
Definition: GSLRndmEngines.cxx:82
ROOT::Math::GSLRandomEngine::GaussianZig
double GaussianZig(double sigma) const
Gaussian distribution - Ziggurat method.
Definition: GSLRndmEngines.cxx:178
ROOT::Math::GSLRandomEngine::Engine
GSLRngWrapper * Engine()
internal method to return the engine Used by class like GSLMCIntegrator to set the engine
Definition: GSLRndmEngines.h:315
GSLRndmEngines.h
ROOT::Math::GSLRandomEngine::Initialize
void Initialize()
initialize the generator If no rng is present the default one based on Mersenne and Twister is create...
Definition: GSLRndmEngines.cxx:93
ROOT::Math::GSLRandomEngine::MaxInt
unsigned long MaxInt() const
return the maximum integer +1 a generator can handle
Definition: GSLRndmEngines.cxx:126
ROOT::Math::GSLRandomEngine::MinInt
unsigned long MinInt() const
return the minimum integer a generator can handle typically this value is 0
Definition: GSLRndmEngines.cxx:121
ROOT::Math::GSLRandomEngine::GaussianTail
double GaussianTail(double a, double sigma) const
Gaussian Tail distribution.
Definition: GSLRndmEngines.cxx:197
ROOT::Math::GSLRngCMRG::GSLRngCMRG
GSLRngCMRG()
Definition: GSLRndmEngines.cxx:390
ROOT::Math::GSLRngROOTWrapper::CreateEngine
static void CreateEngine(gsl_rng *r)
Definition: GSLRngROOTWrapper.h:46
ROOT::Math::GSLRandomEngine::fRng
GSLRngWrapper * fRng
Definition: GSLRndmEngines.h:321
ROOT::Math::GSLRandomEngine::LogNormal
double LogNormal(double zeta, double sigma) const
Log Normal distribution.
Definition: GSLRndmEngines.cxx:240
ROOT::Math::GSLRandomEngine::Rayleigh
double Rayleigh(double sigma) const
Rayleigh distribution.
Definition: GSLRndmEngines.cxx:265
ROOT::Math::GSLRngMixMax::~GSLRngMixMax
virtual ~GSLRngMixMax()
Definition: GSLRndmEngines.cxx:440
ROOT::Math::GSLRandomEngine
GSLRandomEngine Base class for all GSL random engines, normally user instantiate the derived classes ...
Definition: GSLRndmEngines.h:65
ROOT::Math::GSLRandomEngine::~GSLRandomEngine
virtual ~GSLRandomEngine()
call Terminate()
Definition: GSLRndmEngines.cxx:76
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
GSLRngWrapper.h
ROOT::Math::GSLRngMT::GSLRngMT
GSLRngMT()
Definition: GSLRndmEngines.cxx:330
Math
Namespace for new Math classes and functions.
ROOT::Math::GSLRngTaus::GSLRngTaus
GSLRngTaus()
Definition: GSLRndmEngines.cxx:374
ROOT::Math::GSLRandomEngine::GaussianRatio
double GaussianRatio(double sigma) const
Gaussian distribution - Ratio method.
Definition: GSLRndmEngines.cxx:190