Logo ROOT  
Reference Guide
TUnuran.cxx
Go to the documentation of this file.
1// @(#)root/unuran:$Id$
2// Authors: L. Moneta, J. Leydold Tue Sep 26 16:25:09 2006
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10
11// Implementation file for class TUnuran
12
13#include "TUnuran.h"
14
15#include "TUnuranContDist.h"
17#include "TUnuranDiscrDist.h"
18#include "TUnuranEmpDist.h"
19
20#include "UnuranRng.h"
21#include "UnuranDistrAdapter.h"
22
23#include "TRandom.h"
24#include "TSystem.h"
25
26#include "TH1.h"
27
28#include <cassert>
29
30
31#include <unuran.h>
32
33#include "TError.h"
34
35
36TUnuran::TUnuran(TRandom * r, unsigned int debugLevel) :
37 fGen(0),
38 fUdistr(0),
39 fUrng(0),
40 fRng(r)
41{
42 // constructor implementation with a ROOT random generator
43 // if no generator is given the ROOT default is used
44 if (fRng == 0) fRng = gRandom;
45 // set debug level at global level
46 // (should be in a static initialization function of the library ? )
47 if ( debugLevel > 1)
48 unur_set_default_debug(UNUR_DEBUG_ALL);
49 else if (debugLevel == 1)
50 unur_set_default_debug(UNUR_DEBUG_INIT);
51 else
52 unur_set_default_debug(UNUR_DEBUG_OFF);
53
54}
55
56
58{
59 // Destructor implementation
60 if (fGen != 0) unur_free(fGen);
61 if (fUrng != 0) unur_urng_free(fUrng);
62 // we can delete now the distribution object
63 if (fUdistr != 0) unur_distr_free(fUdistr);
64}
65
66//private (no impl.)
68{
69 // Implementation of copy constructor.
70}
71
73{
74 // Implementation of assignment operator.
75 if (this == &rhs) return *this; // time saving self-test
76 return *this;
77}
78
79bool TUnuran::Init(const std::string & dist, const std::string & method)
80{
81 // initialize with a string
82 std::string s = dist + " & " + method;
83 fGen = unur_str2gen(s.c_str() );
84 if (fGen == 0) {
85 Error("Init","Cannot create generator object");
86 return false;
87 }
88 if (! SetRandomGenerator() ) return false;
89
90 return true;
91}
92
93bool TUnuran::Init(const TUnuranContDist & distr, const std::string & method)
94{
95 // initialization with a distribution and and generator
96 // the distribution object is copied in and managed by this class
97 // use std::unique_ptr to manage previously existing distribution objects
98 TUnuranContDist * distNew = distr.Clone();
99 fDist.reset(distNew);
100
101 fMethod = method;
102 if (! SetContDistribution(*distNew) ) return false;
103 if (! SetMethodAndInit() ) return false;
104 if (! SetRandomGenerator() ) return false;
105 return true;
106}
107
108
109bool TUnuran::Init(const TUnuranMultiContDist & distr, const std::string & method)
110{
111 // initialization with a distribution and method
112 // the distribution object is copied in and managed by this class
113 // use std::unique_ptr to manage previously existing distribution objects
114 TUnuranMultiContDist * distNew = distr.Clone();
115 fDist.reset(distNew);
116
117 fMethod = method;
118 if (! SetMultiDistribution(*distNew) ) return false;
119 if (! SetMethodAndInit() ) return false;
120 if (! SetRandomGenerator() ) return false;
121 return true;
122}
123
124
125bool TUnuran::Init(const TUnuranDiscrDist & distr, const std::string & method ) {
126 // initialization with a distribution and and generator
127 // the distribution object is copied in and managed by this class
128 // use std::unique_ptr to manage previously existing distribution objects
129 TUnuranDiscrDist * distNew = distr.Clone();
130 fDist.reset(distNew);
131
132 fMethod = method;
133 if (! SetDiscreteDistribution(*distNew) ) return false;
134 if (! SetMethodAndInit() ) return false;
135 if (! SetRandomGenerator() ) return false;
136 return true;
137}
138
139bool TUnuran::Init(const TUnuranEmpDist & distr, const std::string & method ) {
140 // initialization with a distribution and and generator
141 // the distribution object is copied in and managed by this class
142 // use std::unique_ptr to manage previously existing distribution objects
143 TUnuranEmpDist * distNew = distr.Clone();
144 fDist.reset(distNew);
145
146 fMethod = method;
147 if (distr.IsBinned()) fMethod = "hist";
148 else if (distr.NDim() > 1) fMethod = "vempk";
149 if (! SetEmpiricalDistribution(*distNew) ) return false;
150 if (! SetMethodAndInit() ) return false;
151 if (! SetRandomGenerator() ) return false;
152 return true;
153}
154
155
157{
158 // set an external random generator
159 if (fRng == 0) return false;
160 if (fGen == 0) return false;
161
162 fUrng = unur_urng_new(&UnuranRng<TRandom>::Rndm, fRng );
163 if (fUrng == 0) return false;
164 unsigned int ret = 0;
165 ret |= unur_urng_set_delete(fUrng, &UnuranRng<TRandom>::Delete);
166 ret |= unur_urng_set_seed(fUrng, &UnuranRng<TRandom>::Seed);
167 if (ret != 0) return false;
168
169 unur_chg_urng( fGen, fUrng);
170 return true;
171}
172
174{
175 // internal method to set in unuran the function pointer for a continuous univariate distribution
176 if (fUdistr != 0) unur_distr_free(fUdistr);
177 fUdistr = unur_distr_cont_new();
178 if (fUdistr == 0) return false;
179 unsigned int ret = 0;
180 ret = unur_distr_set_extobj(fUdistr, &dist);
181 if ( ! dist.IsLogPdf() ) {
182 ret |= unur_distr_cont_set_pdf(fUdistr, &ContDist::Pdf);
183 ret |= unur_distr_cont_set_dpdf(fUdistr, &ContDist::Dpdf);
184 if (dist.HasCdf() ) ret |= unur_distr_cont_set_cdf(fUdistr, &ContDist::Cdf);
185 }
186 else {
187 // case user provides log of pdf
188 ret |= unur_distr_cont_set_logpdf(fUdistr, &ContDist::Pdf);
189 ret |= unur_distr_cont_set_dlogpdf(fUdistr, &ContDist::Dpdf);
190 }
191
192 double xmin, xmax = 0;
193 if (dist.GetDomain(xmin,xmax) ) {
194 ret = unur_distr_cont_set_domain(fUdistr,xmin,xmax);
195 if (ret != 0) {
196 Error("SetContDistribution","invalid domain xmin = %g xmax = %g ",xmin,xmax);
197 return false;
198 }
199 }
200 if (dist.HasMode() ) {
201 ret = unur_distr_cont_set_mode(fUdistr, dist.Mode());
202 if (ret != 0) {
203 Error("SetContDistribution","invalid mode given, mode = %g ",dist.Mode());
204 return false;
205 }
206 }
207 if (dist.HasPdfArea() ) {
208 ret = unur_distr_cont_set_pdfarea(fUdistr, dist.PdfArea());
209 if (ret != 0) {
210 Error("SetContDistribution","invalid area given, area = %g ",dist.PdfArea());
211 return false;
212 }
213 }
214
215 return (ret ==0) ? true : false;
216}
217
218
220{
221 // internal method to set in unuran the function pointer for a multivariate distribution
222 if (fUdistr != 0) unur_distr_free(fUdistr);
223 fUdistr = unur_distr_cvec_new(dist.NDim() );
224 if (fUdistr == 0) return false;
225 unsigned int ret = 0;
226 ret |= unur_distr_set_extobj(fUdistr, &dist );
227 if ( ! dist.IsLogPdf() ) {
228 ret |= unur_distr_cvec_set_pdf(fUdistr, &MultiDist::Pdf);
229 ret |= unur_distr_cvec_set_dpdf(fUdistr, &MultiDist::Dpdf);
230 ret |= unur_distr_cvec_set_pdpdf(fUdistr, &MultiDist::Pdpdf);
231 }
232 else {
233 ret |= unur_distr_cvec_set_logpdf(fUdistr, &MultiDist::Pdf);
234 ret |= unur_distr_cvec_set_dlogpdf(fUdistr, &MultiDist::Dpdf);
235 ret |= unur_distr_cvec_set_pdlogpdf(fUdistr, &MultiDist::Pdpdf);
236 }
237
238 const double * xmin = dist.GetLowerDomain();
239 const double * xmax = dist.GetUpperDomain();
240 if ( xmin != 0 || xmax != 0 ) {
241 ret = unur_distr_cvec_set_domain_rect(fUdistr,xmin,xmax);
242 if (ret != 0) {
243 Error("SetMultiDistribution","invalid domain");
244 return false;
245 }
246#ifdef OLDVERS
247 Error("SetMultiDistribution","domain setting not available in UNURAN 0.8.1");
248#endif
249
250 }
251
252 const double * xmode = dist.GetMode();
253 if (xmode != 0) {
254 ret = unur_distr_cvec_set_mode(fUdistr, xmode);
255 if (ret != 0) {
256 Error("SetMultiDistribution","invalid mode");
257 return false;
258 }
259 }
260 return (ret ==0) ? true : false;
261}
262
264
265 // internal method to set in unuran the function pointer for am empiral distribution (from histogram)
266 if (fUdistr != 0) unur_distr_free(fUdistr);
267 if (dist.NDim() == 1)
268 fUdistr = unur_distr_cemp_new();
269 else
270 fUdistr = unur_distr_cvemp_new(dist.NDim() );
271
272 if (fUdistr == 0) return false;
273 unsigned int ret = 0;
274
275
276 // get info from histogram
277 if (dist.IsBinned() ) {
278 int nbins = dist.Data().size();
279 double min = dist.LowerBin();
280 double max = dist.UpperBin();
281 const double * pv = &(dist.Data().front());
282 ret |= unur_distr_cemp_set_hist(fUdistr, pv, nbins, min, max);
283#ifdef OLDVERS
284 Error("SetEmpiricalDistribution","hist method not available in UNURAN 0.8.1");
285#endif
286 }
287 else {
288 const double * pv = &dist.Data().front();
289 // n is number of points (size/ndim)
290 int n = dist.Data().size()/dist.NDim();
291 if (dist.NDim() == 1)
292 ret |= unur_distr_cemp_set_data(fUdistr, pv, n);
293 else
294 ret |= unur_distr_cvemp_set_data(fUdistr, pv, n);
295 }
296 if (ret != 0) {
297 Error("SetEmpiricalDistribution","invalid distribution object");
298 return false;
299 }
300 return true;
301}
302
303
305{
306 // internal method to set in unuran the function pointer for a discrete univariate distribution
307 if (fUdistr != 0) unur_distr_free(fUdistr);
308 fUdistr = unur_distr_discr_new();
309 if (fUdistr == 0) return false;
310 unsigned int ret = 0;
311 // if a probability mesh function is provided
312 if (dist.ProbVec().size() == 0) {
313 ret = unur_distr_set_extobj(fUdistr, &dist );
314 ret |= unur_distr_discr_set_pmf(fUdistr, &DiscrDist::Pmf);
315 if (dist.HasCdf() ) ret |= unur_distr_discr_set_cdf(fUdistr, &DiscrDist::Cdf);
316
317 }
318 else {
319 // case user provides vector of probabilities
320 ret |= unur_distr_discr_set_pv(fUdistr, &dist.ProbVec().front(), dist.ProbVec().size() );
321 }
322
323 int xmin, xmax = 0;
324 if (dist.GetDomain(xmin,xmax) ) {
325 ret = unur_distr_discr_set_domain(fUdistr,xmin,xmax);
326 if (ret != 0) {
327 Error("SetDiscrDistribution","invalid domain xmin = %d xmax = %d ",xmin,xmax);
328 return false;
329 }
330 }
331 if (dist.HasMode() ) {
332 ret = unur_distr_discr_set_mode(fUdistr, dist.Mode());
333 if (ret != 0) {
334 Error("SetContDistribution","invalid mode given, mode = %d ",dist.Mode());
335 return false;
336 }
337 }
338 if (dist.HasProbSum() ) {
339 ret = unur_distr_discr_set_pmfsum(fUdistr, dist.ProbSum());
340 if (ret != 0) {
341 Error("SetContDistribution","invalid sum given, mode = %g ",dist.ProbSum());
342 return false;
343 }
344 }
345
346 return (ret ==0) ? true : false;
347}
348
349
350//bool TUnuran::SetMethodAndInit(const std::string & s) {
352
353 // internal function to set a method from a distribution and
354 // initialize unuran with the given distribution.
355 if (fUdistr == 0) return false;
356
357 struct unur_slist *mlist = NULL;
358
359 UNUR_PAR * par = _unur_str2par(fUdistr, fMethod.c_str(), &mlist);
360 if (par == 0) {
361 Error("SetMethod","missing distribution information or syntax error");
362 if (mlist != 0) _unur_slist_free(mlist);
363 return false;
364 }
365
366
367 // set unuran to not use a private copy of the distribution object
368 unur_set_use_distr_privatecopy (par, false);
369
370 // need to free fGen if already existing ?
371 if (fGen != 0 ) unur_free(fGen);
372 fGen = unur_init(par);
373 _unur_slist_free(mlist);
374 if (fGen == 0) {
375 Error("SetMethod","initializing Unuran: condition for method violated");
376 return false;
377 }
378 return true;
379 }
380
381
383{
384 // sample one-dimensional distribution
385 assert(fGen != 0);
386 return unur_sample_discr(fGen);
387}
388
390{
391 // sample one-dimensional distribution
392 assert(fGen != 0);
393 return unur_sample_cont(fGen);
394}
395
397{
398 // sample multidimensional distribution
399 if (fGen == 0) return false;
400 unur_sample_vec(fGen,x);
401 return true;
402}
403
404void TUnuran::SetSeed(unsigned int seed) {
405 return fRng->SetSeed(seed);
406}
407
408bool TUnuran::SetLogLevel(unsigned int debugLevel)
409{
410 if (fGen == 0) return false;
411 int ret = 0;
412 if ( debugLevel > 1)
413 ret |= unur_chg_debug(fGen, UNUR_DEBUG_ALL);
414 else if (debugLevel == 1)
415 ret |= unur_chg_debug(fGen, UNUR_DEBUG_ALL);
416 else
417 ret |= unur_chg_debug(fGen, UNUR_DEBUG_OFF);
418
419 return (ret ==0) ? true : false;
420
421}
422
423bool TUnuran::InitPoisson(double mu, const std::string & method) {
424 // initializaton for a Poisson
425 double p[1];
426 p[0] = mu;
427
428 fUdistr = unur_distr_poisson(p,1);
429
430 fMethod = method;
431 if (fUdistr == 0) return false;
432 if (! SetMethodAndInit() ) return false;
433 if (! SetRandomGenerator() ) return false;
434 return true;
435}
436
437
438bool TUnuran::InitBinomial(unsigned int ntot, double prob, const std::string & method ) {
439 // initializaton for a Binomial
440 double par[2];
441 par[0] = ntot;
442 par[1] = prob;
443 fUdistr = unur_distr_binomial(par,2);
444
445 fMethod = method;
446 if (fUdistr == 0) return false;
447 if (! SetMethodAndInit() ) return false;
448 if (! SetRandomGenerator() ) return false;
449 return true;
450}
451
452
453bool TUnuran::ReInitDiscrDist(unsigned int npar, double * par) {
454 // re-initialization of UNURAN without freeing and creating a new fGen object
455 // works only for pre-defined distribution by changing their parameters
456 if (!fGen ) return false;
457 if (!fUdistr) return false;
458 unur_distr_discr_set_pmfparams(fUdistr,par,npar);
459 int iret = unur_reinit(fGen);
460 if (iret) Warning("ReInitDiscrDist","re-init failed - a full initizialization must be performed");
461 return (!iret);
462}
463
ROOT::R::TRInterface & r
Definition: Object.C:4
void Error(const char *location, const char *msgfmt,...)
void Warning(const char *location, const char *msgfmt,...)
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:597
TUnuranContDist class describing one dimensional continuous distribution.
virtual TUnuranContDist * Clone() const
Clone (required by base class)
TUnuranDiscrDist class for one dimensional discrete distribution.
virtual TUnuranDiscrDist * Clone() const
Clone (required by base class)
TUnuranEmpDist class for describing empiral distributions.
unsigned int NDim() const
Number of data dimensions.
TUnuranEmpDist * Clone() const
Clone (required by base class)
bool IsBinned() const
Flag to control if data are binned.
TUnuranMultiContDist class describing multi dimensional continuous distributions.
virtual TUnuranMultiContDist * Clone() const
Clone (required by base class)
TUnuran class.
Definition: TUnuran.h:79
TUnuran(TRandom *r=0, unsigned int log=0)
Constructor with a generator instance and given level of log output.
Definition: TUnuran.cxx:36
bool SetMethodAndInit()
change the method and initialize Unuran with the previously given distribution
Definition: TUnuran.cxx:351
int SampleDiscr()
Sample discrete distributions User is responsible for having previously correctly initialized with TU...
Definition: TUnuran.cxx:382
std::unique_ptr< TUnuranBaseDist > fDist
Definition: TUnuran.h:264
bool SetDiscreteDistribution(const TUnuranDiscrDist &dist)
Definition: TUnuran.cxx:304
bool InitBinomial(unsigned int ntot, double prob, const std::string &method="dstd")
Initialize method for the Binomial distribution Used to generate poisson numbers for a constant param...
Definition: TUnuran.cxx:438
bool SetContDistribution(const TUnuranContDist &dist)
Definition: TUnuran.cxx:173
bool SampleMulti(double *x)
Sample multidimensional distributions User is responsible for having previously correctly initialized...
Definition: TUnuran.cxx:396
bool ReInitDiscrDist(unsigned int npar, double *params)
Reinitialize UNURAN by changing the distribution parameters but mantaining same distribution and meth...
Definition: TUnuran.cxx:453
UNUR_DISTR * fUdistr
Definition: TUnuran.h:262
bool Init(const std::string &distr, const std::string &method)
initialize with Unuran string interface
Definition: TUnuran.cxx:79
UNUR_GEN * fGen
Definition: TUnuran.h:261
bool SetRandomGenerator()
Definition: TUnuran.cxx:156
bool SetLogLevel(unsigned int iflag=1)
set log level
Definition: TUnuran.cxx:408
UNUR_URNG * fUrng
Definition: TUnuran.h:263
TRandom * fRng
Definition: TUnuran.h:265
bool SetMultiDistribution(const TUnuranMultiContDist &dist)
Definition: TUnuran.cxx:219
double Sample()
Sample 1D distribution User is responsible for having previously correctly initialized with TUnuran::...
Definition: TUnuran.cxx:389
bool SetEmpiricalDistribution(const TUnuranEmpDist &dist)
Definition: TUnuran.cxx:263
~TUnuran()
Destructor.
Definition: TUnuran.cxx:57
bool InitPoisson(double mu, const std::string &method="dstd")
Initialize method for the Poisson distribution Used to generate poisson numbers for a constant parame...
Definition: TUnuran.cxx:423
std::string fMethod
Definition: TUnuran.h:266
TUnuran & operator=(const TUnuran &rhs)
Assignment operator.
Definition: TUnuran.cxx:72
void SetSeed(unsigned int seed)
set the seed for the random number generator
Definition: TUnuran.cxx:404
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
static constexpr double s
static double Cdf(double x, const UNUR_DISTR *dist)
evaluate the Cumulative distribution function, integral of the pdf
static double Pdf(double x, const UNUR_DISTR *dist)
evaluate the probality density function
static double Dpdf(double x, const UNUR_DISTR *dist)
evaluate the derivative of the pdf
static double Pmf(int x, const UNUR_DISTR *dist)
evaluate the probality mesh function
static double Cdf(int x, const UNUR_DISTR *dist)
evaluate the cumulative function
static double Pdf(const double *x, UNUR_DISTR *dist)
evaluate the probality density function
static int Dpdf(double *grad, const double *x, UNUR_DISTR *dist)
static double Pdpdf(const double *x, int coord, UNUR_DISTR *dist)
UnuranRng class for interface ROOT random generators to Unuran.
Definition: UnuranRng.h:22