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_linalg.h"
39#include "gsl/gsl_matrix.h"
40#include "gsl/gsl_rng.h"
41#include "gsl/gsl_randist.h"
42#include "gsl/gsl_vector.h"
43#include "gsl/gsl_version.h"
44
45#include "Math/GSLRndmEngines.h"
46#include "GSLRngWrapper.h"
47// for wrapping in GSL ROOT engines
48#include "GSLRngROOTWrapper.h"
49
50extern double gsl_ran_gaussian_acr( const gsl_rng * r, const double sigma);
51
52// gsl_multivarate_gaussian was added in GSL 2.2
53// For older GSL versions (e.g. Ubuntu 16.04 comes with GSL 2.1) we can add it here by hand
54// from: http://git.savannah.gnu.org/cgit/gsl.git/tree/randist/mvgauss.c?h=release-2-6&id=8f0165f5cb2ae02e386cd33ff10e47ffb46ea7da
55#if (GSL_MAJOR_VERSION == 1) || ((GSL_MAJOR_VERSION == 2) && (GSL_MINOR_VERSION < 2))
56#include <gsl/gsl_blas.h>
57extern int
58gsl_ran_multivariate_gaussian(const gsl_rng *r, const gsl_vector *mu, const gsl_matrix *L, gsl_vector *result)
59{
60 const size_t M = L->size1;
61 const size_t N = L->size2;
62
63 if (M != N) {
64 GSL_ERROR("requires square matrix", GSL_ENOTSQR);
65 } else if (mu->size != M) {
66 GSL_ERROR("incompatible dimension of mean vector with variance-covariance matrix", GSL_EBADLEN);
67 } else if (result->size != M) {
68 GSL_ERROR("incompatible dimension of result vector", GSL_EBADLEN);
69 } else {
70 size_t i;
71
72 for (i = 0; i < M; ++i)
73 gsl_vector_set(result, i, gsl_ran_ugaussian(r));
74
75 gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasNonUnit, L, result);
76 gsl_vector_add(result, mu);
77
78 return GSL_SUCCESS;
79 }
80}
81#endif
82
83namespace ROOT {
84namespace Math {
85
86
87
88
89
90 // default constructor (need to call set type later)
92 fRng(nullptr),
93 fCurTime(0)
94 { }
95
96 // constructor from external rng
97 // internal generator will be managed or not depending on
98 // how the GSLRngWrapper is created
100 fRng(new GSLRngWrapper(*rng) ),
101 fCurTime(0)
102 {}
103
104 // copy constructor
106 fRng(new GSLRngWrapper(*eng.fRng) ),
107 fCurTime(0)
108 {}
109
111 // destructor : call terminate if not yet called
112 if (fRng) Terminate();
113 }
114
115 // assignment operator
117 if (this == &eng) return *this;
118 if (fRng)
119 *fRng = *eng.fRng;
120 else
121 fRng = new GSLRngWrapper(*eng.fRng);
122 fCurTime = eng.fCurTime;
123 return *this;
124 }
125
126
128 // initialize the generator by allocating the GSL object
129 // if type was not passed create with default generator
130 if (!fRng) fRng = new GSLRngWrapper();
131 fRng->Allocate();
132 }
133
135 // terminate the generator by freeing the GSL object
136 if (!fRng) return;
137 fRng->Free();
138 delete fRng;
139 fRng = 0;
140 }
141
142
144 // generate random between 0 and 1.
145 // 0 is excluded
146 return gsl_rng_uniform_pos(fRng->Rng() );
147 }
148
149
150 unsigned long GSLRandomEngine::RndmInt(unsigned long max) const {
151 // generate a random integer number between 0 and MAX
152 return gsl_rng_uniform_int( fRng->Rng(), max );
153 }
154
155 unsigned long GSLRandomEngine::MinInt() const {
156 // return minimum integer value used in RndmInt
157 return gsl_rng_min( fRng->Rng() );
158 }
159
160 unsigned long GSLRandomEngine::MaxInt() const {
161 // return maximum integr value used in RndmInt
162 return gsl_rng_max( fRng->Rng() );
163 }
164
165 void GSLRandomEngine::RandomArray(double * begin, double * end ) const {
166 // generate array of randoms betweeen 0 and 1. 0 is excluded
167 // specialization for double * (to be faster)
168 for ( double * itr = begin; itr != end; ++itr ) {
169 *itr = gsl_rng_uniform_pos(fRng->Rng() );
170 }
171 }
172
173 void GSLRandomEngine::SetSeed(unsigned int seed) const {
174 // set the seed, if = 0then the seed is set randomly using an std::rand()
175 // seeded with the current time. Be carefuk in case the current time is
176 // the same in consecutive calls
177 if (seed == 0) {
178 // use like in root (use time)
179 time_t curtime;
180 time(&curtime);
181 unsigned int ct = static_cast<unsigned int>(curtime);
182 if (ct != fCurTime) {
183 fCurTime = ct;
184 // set the seed for rand
185 srand(ct);
186 }
187 seed = rand();
188 }
189
190 assert(fRng);
191 gsl_rng_set(fRng->Rng(), seed );
192 }
193
194 std::string GSLRandomEngine::Name() const {
195 //////////////////////////////////////////////////////////////////////////
196
197 assert ( fRng != 0);
198 assert ( fRng->Rng() != 0 );
199 return std::string( gsl_rng_name( fRng->Rng() ) );
200 }
201
202 unsigned int GSLRandomEngine::Size() const {
203 //////////////////////////////////////////////////////////////////////////
204
205 assert (fRng != 0);
206 return gsl_rng_size( fRng->Rng() );
207 }
208
209
210 // Random distributions
211
213 {
214 // Gaussian distribution. Use fast ziggurat algorithm implemented since GSL 1.8
215 return gsl_ran_gaussian_ziggurat( fRng->Rng(), sigma);
216 }
217
218 double GSLRandomEngine::Gaussian(double sigma) const
219 {
220 // Gaussian distribution. Use default Box-Muller method
221 return gsl_ran_gaussian( fRng->Rng(), sigma);
222 }
223
225 {
226 // Gaussian distribution. Use ratio method
227 return gsl_ran_gaussian_ratio_method( fRng->Rng(), sigma);
228 }
229
230
231 double GSLRandomEngine::GaussianTail(double a , double sigma) const
232 {
233 // Gaussian Tail distribution: eeturn values larger than a distributed
234 // according to the gaussian
235 return gsl_ran_gaussian_tail( fRng->Rng(), a, sigma);
236 }
237
238 void GSLRandomEngine::Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const
239 {
240 // Gaussian Bivariate distribution, with correlation coefficient rho
241 gsl_ran_bivariate_gaussian( fRng->Rng(), sigmaX, sigmaY, rho, &x, &y);
242 }
243
244 void GSLRandomEngine::GaussianND(const int dim, double *pars, double *covmat, double *genpars) const
245 {
246 // Gaussian Multivariate distribution
247 gsl_vector *mu = gsl_vector_alloc(dim);
248 gsl_vector *genpars_vec = gsl_vector_alloc(dim);
249 gsl_matrix *L = gsl_matrix_alloc(dim, dim);
250 for (int i = 0; i < dim; ++i) {
251 gsl_vector_set(mu, i, pars[i]);
252 for (int j = 0; j < dim; ++j) {
253 gsl_matrix_set(L, i, j, covmat[i * dim + j]);
254 }
255 }
256 gsl_linalg_cholesky_decomp(L);
257 gsl_ran_multivariate_gaussian(fRng->Rng(), mu, L, genpars_vec);
258 for (int i = 0; i < dim; ++i) {
259 genpars[i] = gsl_vector_get(genpars_vec, i);
260 }
261 }
262
263 double GSLRandomEngine::Exponential(double mu) const
264 {
265 // Exponential distribution
266 return gsl_ran_exponential( fRng->Rng(), mu);
267 }
268
269 double GSLRandomEngine::Cauchy(double a) const
270 {
271 // Cauchy distribution
272 return gsl_ran_cauchy( fRng->Rng(), a);
273 }
274
276 {
277 // Landau distribution
278 return gsl_ran_landau( fRng->Rng());
279 }
280
281 double GSLRandomEngine::Beta(double a, double b) const
282 {
283 // Beta distribution
284 return gsl_ran_beta( fRng->Rng(), a, b);
285 }
286
287 double GSLRandomEngine::Gamma(double a, double b) const
288 {
289 // Gamma distribution
290 return gsl_ran_gamma( fRng->Rng(), a, b);
291 }
292
293 double GSLRandomEngine::LogNormal(double zeta, double sigma) const
294 {
295 // Log normal distribution
296 return gsl_ran_lognormal( fRng->Rng(), zeta, sigma);
297 }
298
299 double GSLRandomEngine::ChiSquare(double nu) const
300 {
301 // Chi square distribution
302 return gsl_ran_chisq( fRng->Rng(), nu);
303 }
304
305
306 double GSLRandomEngine::FDist(double nu1, double nu2) const
307 {
308 // F distribution
309 return gsl_ran_fdist( fRng->Rng(), nu1, nu2);
310 }
311
312 double GSLRandomEngine::tDist(double nu) const
313 {
314 // t distribution
315 return gsl_ran_tdist( fRng->Rng(), nu);
316 }
317
318 double GSLRandomEngine::Rayleigh(double sigma) const
319 {
320 // Rayleigh distribution
321 return gsl_ran_rayleigh( fRng->Rng(), sigma);
322 }
323
324 double GSLRandomEngine::Logistic(double a) const
325 {
326 // Logistic distribution
327 return gsl_ran_logistic( fRng->Rng(), a);
328 }
329
330 double GSLRandomEngine::Pareto(double a, double b) const
331 {
332 // Pareto distribution
333 return gsl_ran_pareto( fRng->Rng(), a, b);
334 }
335
336 void GSLRandomEngine::Dir2D(double &x, double &y) const
337 {
338 // generate random numbers in a 2D circle of radious 1
339 gsl_ran_dir_2d( fRng->Rng(), &x, &y);
340 }
341
342 void GSLRandomEngine::Dir3D(double &x, double &y, double &z) const
343 {
344 // generate random numbers in a 3D sphere of radious 1
345 gsl_ran_dir_3d( fRng->Rng(), &x, &y, &z);
346 }
347
348 unsigned int GSLRandomEngine::Poisson(double mu) const
349 {
350 // Poisson distribution
351 return gsl_ran_poisson( fRng->Rng(), mu);
352 }
353
354 unsigned int GSLRandomEngine::Binomial(double p, unsigned int n) const
355 {
356 // Binomial distribution
357 return gsl_ran_binomial( fRng->Rng(), p, n);
358 }
359
360 unsigned int GSLRandomEngine::NegativeBinomial(double p, double n) const
361 {
362 // Negative Binomial distribution
363 return gsl_ran_negative_binomial( fRng->Rng(), p, n);
364 }
365
366
367 std::vector<unsigned int> GSLRandomEngine::Multinomial( unsigned int ntot, const std::vector<double> & p ) const
368 {
369 // Multinomial distribution return vector of integers which sum is ntot
370 std::vector<unsigned int> ival( p.size());
371 gsl_ran_multinomial( fRng->Rng(), p.size(), ntot, &p.front(), &ival[0]);
372 return ival;
373 }
374
375
376
377 //----------------------------------------------------
378 // generators
379 //----------------------------------------------------
380
381 /////////////////////////////////////////////////////////////////////////////
382
384 {
385 SetType(new GSLRngWrapper(gsl_rng_mt19937));
386 Initialize();
387 }
388
389
390 // old ranlux - equivalent to TRandom1
392 {
393 SetType(new GSLRngWrapper(gsl_rng_ranlux) );
394 Initialize();
395 }
396
397 // second generation of Ranlux (single precision version - luxury 1)
399 {
400 SetType(new GSLRngWrapper(gsl_rng_ranlxs1) );
401 Initialize();
402 }
403
404 // second generation of Ranlux (single precision version - luxury 2)
406 {
407 SetType(new GSLRngWrapper(gsl_rng_ranlxs2) );
408 Initialize();
409 }
410
411 // double precision version - luxury 1
413 {
414 SetType(new GSLRngWrapper(gsl_rng_ranlxd1) );
415 Initialize();
416 }
417
418 // double precision version - luxury 2
420 {
421 SetType(new GSLRngWrapper(gsl_rng_ranlxd2) );
422 Initialize();
423 }
424
425 /////////////////////////////////////////////////////////////////////////////
426
428 {
429 SetType(new GSLRngWrapper(gsl_rng_taus2) );
430 Initialize();
431 }
432
433 /////////////////////////////////////////////////////////////////////////////
434
436 {
437 SetType(new GSLRngWrapper(gsl_rng_gfsr4) );
438 Initialize();
439 }
440
441 /////////////////////////////////////////////////////////////////////////////
442
444 {
445 SetType(new GSLRngWrapper(gsl_rng_cmrg) );
446 Initialize();
447 }
448
449 /////////////////////////////////////////////////////////////////////////////
450
452 {
453 SetType(new GSLRngWrapper(gsl_rng_mrg) );
454 Initialize();
455 }
456
457
458 /////////////////////////////////////////////////////////////////////////////
459
461 {
462 SetType(new GSLRngWrapper(gsl_rng_rand) );
463 Initialize();
464 }
465
466 /////////////////////////////////////////////////////////////////////////////
467
469 {
470 SetType(new GSLRngWrapper(gsl_rng_ranmar) );
471 Initialize();
472 }
473
474 /////////////////////////////////////////////////////////////////////////////
475
477 {
478 SetType(new GSLRngWrapper(gsl_rng_minstd) );
479 Initialize();
480 }
481
482
483 // for extra engines based on ROOT
485 {
487 Initialize(); // this creates the gsl_rng structure
488 // no real need to call CreateEngine since the underlined MIXMAX engine is created
489 // by calling GSLMixMaxWrapper::Seed(gsl_default_seed) that is called
490 // when gsl_rng is allocated (in Initialize)
492 }
494 // we need to explicitly delete the ROOT wrapper class
496 }
497
498} // namespace Math
499} // namespace ROOT
double gsl_ran_gaussian_acr(const gsl_rng *r, const double sigma)
const gsl_rng_type * gsl_rng_mixmax
#define GSL_ERROR(a, b)
#define N
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
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 GaussianND(const int dim, double *pars, double *covmat, double *genpars) const
Multivariate Gaussian distribution.
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.
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
static constexpr double L
static void FreeEngine(gsl_rng *r)
static void CreateEngine(gsl_rng *r)
auto * a
Definition: textangle.C:12