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