Logo ROOT   6.18/05
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
47extern double gsl_ran_gaussian_acr( const gsl_rng * r, const double sigma);
48
49//#include <iostream>
50
51namespace ROOT {
52namespace 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
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
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 gsl_ran_gaussian_acr(const gsl_rng *r, const double sigma)
const gsl_rng_type * gsl_rng_mixmax
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
GSLRandomEngine Base class for all GSL random engines, normally user instantiate the derived classes ...
double operator()() const
Generate a random number between ]0,1] 0 is excluded and 1 is included.
void Dir3D(double &x, double &y, double &z) const
generate random numbers in a 3D sphere of radious 1
void SetSeed(unsigned int seed) const
set the random generator seed
void Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const
Bivariate Gaussian distribution with correlation.
unsigned int Poisson(double mu) const
Poisson distribution.
double ChiSquare(double nu) const
Chi square distribution.
double FDist(double nu1, double nu2) const
F distrbution.
GSLRngWrapper * Engine()
internal method to return the engine Used by class like GSLMCIntegrator to set the engine
double Gamma(double a, double b) const
Gamma distribution.
GSLRandomEngine()
default constructor.
double Exponential(double mu) const
Exponential distribution.
double GaussianTail(double a, double sigma) const
Gaussian Tail distribution.
double Rayleigh(double sigma) const
Rayleigh distribution.
std::vector< unsigned int > Multinomial(unsigned int ntot, const std::vector< double > &p) const
Multinomial distribution.
unsigned long MaxInt() const
return the maximum integer +1 a generator can handle
double Cauchy(double a) const
Cauchy distribution.
double GaussianZig(double sigma) const
Gaussian distribution - Ziggurat method.
unsigned long MinInt() const
return the minimum integer a generator can handle typically this value is 0
double Gaussian(double sigma) const
Gaussian distribution - default method is Box-Muller (polar method)
double Beta(double a, double b) const
Beta distribution.
double Landau() const
Landau distribution.
unsigned int NegativeBinomial(double p, double n) const
Negative Binomial distribution.
void RandomArray(Iterator begin, Iterator end) const
Generate an array of random numbers.
unsigned int Size() const
return the state size of generator
double Pareto(double a, double b) const
Pareto distribution.
double LogNormal(double zeta, double sigma) const
Log Normal distribution.
unsigned int Binomial(double p, unsigned int n) const
Binomial distribution.
void Dir2D(double &x, double &y) const
generate random numbers in a 2D circle of radious 1
void Terminate()
delete pointer to contained rng
std::string Name() const
return name of generator
double GaussianRatio(double sigma) const
Gaussian distribution - Ratio method.
double Logistic(double a) const
Logistic 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 ...
void SetType(GSLRngWrapper *r)
internal method used by the derived class to set the type of generators
void Initialize()
initialize the generator If no rng is present the default one based on Mersenne and Twister is create...
double tDist(double nu) const
t student distribution
GSLRandomEngine & operator=(const GSLRandomEngine &eng)
Assignment operator : make a deep copy of the contained GSL generator.
virtual ~GSLRandomEngine()
call Terminate()
GSLRngWrapper class to wrap gsl_rng structure.
Definition: GSLRngWrapper.h:25
const Double_t sigma
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
auto * a
Definition: textangle.C:12