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
377
379{
380 // sample one-dimensional distribution
381 assert(fGen != 0);
382 return unur_sample_discr(fGen);
383}
384
386{
387 // sample one-dimensional distribution
388 assert(fGen != 0);
389 return unur_sample_cont(fGen);
390}
391
393{
394 // sample multidimensional distribution
395 if (fGen == 0) return false;
396 unur_sample_vec(fGen,x);
397 return true;
398}
399
400void TUnuran::SetSeed(unsigned int seed) {
401 return fRng->SetSeed(seed);
402}
403
404bool TUnuran::SetLogLevel(unsigned int debugLevel)
405{
406 if (fGen == 0) return false;
407 int ret = 0;
408 if ( debugLevel > 1)
409 ret |= unur_chg_debug(fGen, UNUR_DEBUG_ALL);
410 else if (debugLevel == 1)
411 ret |= unur_chg_debug(fGen, UNUR_DEBUG_ALL);
412 else
413 ret |= unur_chg_debug(fGen, UNUR_DEBUG_OFF);
414
415 return (ret ==0) ? true : false;
416
417}
418
419bool TUnuran::InitPoisson(double mu, const std::string & method) {
420 // initializaton for a Poisson
421 double p[1];
422 p[0] = mu;
423
424 fUdistr = unur_distr_poisson(p,1);
425
426 fMethod = method;
427 if (fUdistr == 0) return false;
428 if (! SetMethodAndInit() ) return false;
429 if (! SetRandomGenerator() ) return false;
430 return true;
431}
432
433
434bool TUnuran::InitBinomial(unsigned int ntot, double prob, const std::string & method ) {
435 // initializaton for a Binomial
436 double par[2];
437 par[0] = ntot;
438 par[1] = prob;
439 fUdistr = unur_distr_binomial(par,2);
440
441 fMethod = method;
442 if (fUdistr == 0) return false;
443 if (! SetMethodAndInit() ) return false;
444 if (! SetRandomGenerator() ) return false;
445 return true;
446}
447
448
449bool TUnuran::ReInitDiscrDist(unsigned int npar, double * par) {
450 // re-initialization of UNURAN without freeing and creating a new fGen object
451 // works only for pre-defined distribution by changing their parameters
452 if (!fGen ) return false;
453 if (!fUdistr) return false;
454 unur_distr_discr_set_pmfparams(fUdistr,par,npar);
455 int iret = unur_reinit(fGen);
456 if (iret) Warning("ReInitDiscrDist","re-init failed - a full initizialization must be performed");
457 return (!iret);
458}
459
ROOT::R::TRInterface & r
Definition: Object.C:4
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.
bool IsBinned() const
Flag to control if data are binned.
TUnuranEmpDist * Clone() const
Clone (required by base class)
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:32
bool SetMethodAndInit()
change the method and initialize Unuran with the previously given distribution
Definition: TUnuran.cxx:347
int SampleDiscr()
Sample discrete distributions User is responsible for having previously correctly initialized with TU...
Definition: TUnuran.cxx:378
std::unique_ptr< TUnuranBaseDist > fDist
Definition: TUnuran.h:264
bool SetDiscreteDistribution(const TUnuranDiscrDist &dist)
Definition: TUnuran.cxx:300
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:434
bool SetContDistribution(const TUnuranContDist &dist)
Definition: TUnuran.cxx:169
bool SampleMulti(double *x)
Sample multidimensional distributions User is responsible for having previously correctly initialized...
Definition: TUnuran.cxx:392
bool ReInitDiscrDist(unsigned int npar, double *params)
Reinitialize UNURAN by changing the distribution parameters but mantaining same distribution and meth...
Definition: TUnuran.cxx:449
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:75
UNUR_GEN * fGen
Definition: TUnuran.h:261
bool SetRandomGenerator()
Definition: TUnuran.cxx:152
bool SetLogLevel(unsigned int iflag=1)
set log level
Definition: TUnuran.cxx:404
UNUR_URNG * fUrng
Definition: TUnuran.h:263
TRandom * fRng
Definition: TUnuran.h:265
bool SetMultiDistribution(const TUnuranMultiContDist &dist)
Definition: TUnuran.cxx:215
double Sample()
Sample 1D distribution User is responsible for having previously correctly initialized with TUnuran::...
Definition: TUnuran.cxx:385
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 Used to generate poisson numbers for a constant parame...
Definition: TUnuran.cxx:419
std::string fMethod
Definition: TUnuran.h:266
TUnuran & operator=(const TUnuran &rhs)
Assignment operator.
Definition: TUnuran.cxx:68
void SetSeed(unsigned int seed)
set the seed for the random number generator
Definition: TUnuran.cxx:400
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
void Error(const char *location, const char *va_(fmt),...)
Definition: TClingUtils.h:789
void Warning(const char *location, const char *va_(fmt),...)
Definition: TClingUtils.h:819
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