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"
16 #include "TUnuranMultiContDist.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 
32 TUnuran::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 
75 bool 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 
89 bool 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 
105 bool 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 
121 bool 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 
135 bool 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 
392 bool TUnuran::SampleMulti(double * x)
393 {
394  // sample multidimensional distribution
395  if (fGen == 0) return false;
396  unur_sample_vec(fGen,x);
397  return true;
398 }
399 
400 void TUnuran::SetSeed(unsigned int seed) {
401  return fRng->SetSeed(seed);
402 }
403 
404 bool 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 
419 bool 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 
434 bool 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 
449 bool 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 
MultiDist::Dpdf
static int Dpdf(double *grad, const double *x, UNUR_DISTR *dist)
Definition: UnuranDistrAdapter.h:66
TUnuran::TUnuran
TUnuran(TRandom *r=0, unsigned int log=0)
Constructor with a generator instance and given level of log output.
Definition: TUnuran.cxx:32
n
const Int_t n
Definition: legend1.C:16
TUnuran.h
TUnuran::fGen
UNUR_GEN * fGen
Definition: TUnuran.h:261
TUnuranEmpDist::Clone
TUnuranEmpDist * Clone() const
Clone (required by base class)
Definition: TUnuranEmpDist.h:108
Warning
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition: TError.cxx:231
UnuranRng
UnuranRng class for interface ROOT random generators to Unuran.
Definition: UnuranRng.h:22
TUnuranMultiContDist.h
TUnuran::SetMultiDistribution
bool SetMultiDistribution(const TUnuranMultiContDist &dist)
Definition: TUnuran.cxx:215
TUnuranEmpDist.h
r
ROOT::R::TRInterface & r
Definition: Object.C:4
xmax
float xmax
Definition: THbookFile.cxx:95
TUnuran::SetSeed
void SetSeed(unsigned int seed)
set the seed for the random number generator
Definition: TUnuran.cxx:400
ContDist::Pdf
static double Pdf(double x, const UNUR_DISTR *dist)
evaluate the probality density function
Definition: UnuranDistrAdapter.h:34
TUnuran::SampleDiscr
int SampleDiscr()
Sample discrete distributions User is responsible for having previously correctly initialized with TU...
Definition: TUnuran.cxx:378
TUnuran::~TUnuran
~TUnuran()
Destructor.
Definition: TUnuran.cxx:53
TUnuran::SetDiscreteDistribution
bool SetDiscreteDistribution(const TUnuranDiscrDist &dist)
Definition: TUnuran.cxx:300
TRandom.h
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
x
Double_t x[n]
Definition: legend1.C:17
ContDist::Dpdf
static double Dpdf(double x, const UNUR_DISTR *dist)
evaluate the derivative of the pdf
Definition: UnuranDistrAdapter.h:39
TUnuranContDist.h
TUnuranDiscrDist.h
TUnuran::SetLogLevel
bool SetLogLevel(unsigned int iflag=1)
set log level
Definition: TUnuran.cxx:404
TUnuran::fUdistr
UNUR_DISTR * fUdistr
Definition: TUnuran.h:262
TUnuran::InitPoisson
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
TUnuranMultiContDist::Clone
virtual TUnuranMultiContDist * Clone() const
Clone (required by base class)
Definition: TUnuranMultiContDist.h:86
ROOT::Math::gv_detail::dist
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
MultiDist::Pdf
static double Pdf(const double *x, UNUR_DISTR *dist)
evaluate the probality density function
Definition: UnuranDistrAdapter.h:60
TUnuran::SetRandomGenerator
bool SetRandomGenerator()
Definition: TUnuran.cxx:152
DiscrDist::Cdf
static double Cdf(int x, const UNUR_DISTR *dist)
evaluate the cumulative function
Definition: UnuranDistrAdapter.h:96
TUnuran::fRng
TRandom * fRng
Definition: TUnuran.h:265
TRandom
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
TUnuranContDist
TUnuranContDist class describing one dimensional continuous distribution.
Definition: TUnuranContDist.h:48
TUnuranEmpDist::IsBinned
bool IsBinned() const
Flag to control if data are binned.
Definition: TUnuranEmpDist.h:119
xmin
float xmin
Definition: THbookFile.cxx:95
TUnuranContDist::Clone
virtual TUnuranContDist * Clone() const
Clone (required by base class)
Definition: TUnuranContDist.h:91
TUnuran::Init
bool Init(const std::string &distr, const std::string &method)
initialize with Unuran string interface
Definition: TUnuran.cxx:75
gRandom
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
TUnuran::fUrng
UNUR_URNG * fUrng
Definition: TUnuran.h:263
TRandom::SetSeed
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:597
TUnuranMultiContDist
TUnuranMultiContDist class describing multi dimensional continuous distributions.
Definition: TUnuranMultiContDist.h:47
TUnuranDiscrDist::Clone
virtual TUnuranDiscrDist * Clone() const
Clone (required by base class)
Definition: TUnuranDiscrDist.h:101
TUnuranEmpDist::NDim
unsigned int NDim() const
Number of data dimensions.
Definition: TUnuranEmpDist.h:136
TUnuran::operator=
TUnuran & operator=(const TUnuran &rhs)
Assignment operator.
Definition: TUnuran.cxx:68
DiscrDist::Pmf
static double Pmf(int x, const UNUR_DISTR *dist)
evaluate the probality mesh function
Definition: UnuranDistrAdapter.h:90
TUnuran::InitBinomial
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
MultiDist::Pdpdf
static double Pdpdf(const double *x, int coord, UNUR_DISTR *dist)
Definition: UnuranDistrAdapter.h:73
UnuranRng.h
TUnuran::ReInitDiscrDist
bool ReInitDiscrDist(unsigned int npar, double *params)
Reinitialize UNURAN by changing the distribution parameters but mantaining same distribution and meth...
Definition: TUnuran.cxx:449
UnuranDistrAdapter.h
TUnuran::SampleMulti
bool SampleMulti(double *x)
Sample multidimensional distributions User is responsible for having previously correctly initialized...
Definition: TUnuran.cxx:392
TUnuran::SetMethodAndInit
bool SetMethodAndInit()
change the method and initialize Unuran with the previously given distribution
Definition: TUnuran.cxx:347
TUnuran::SetContDistribution
bool SetContDistribution(const TUnuranContDist &dist)
Definition: TUnuran.cxx:169
TUnuranDiscrDist
TUnuranDiscrDist class for one dimensional discrete distribution.
Definition: TUnuranDiscrDist.h:51
ContDist::Cdf
static double Cdf(double x, const UNUR_DISTR *dist)
evaluate the Cumulative distribution function, integral of the pdf
Definition: UnuranDistrAdapter.h:45
TUnuran::SetEmpiricalDistribution
bool SetEmpiricalDistribution(const TUnuranEmpDist &dist)
Definition: TUnuran.cxx:259
TUnuran::fMethod
std::string fMethod
Definition: TUnuran.h:266
TUnuran
TUnuran class.
Definition: TUnuran.h:79
TUnuranEmpDist
TUnuranEmpDist class for describing empiral distributions.
Definition: TUnuranEmpDist.h:49
TUnuran::Sample
double Sample()
Sample 1D distribution User is responsible for having previously correctly initialized with TUnuran::...
Definition: TUnuran.cxx:385
Error
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
TError.h
TUnuran::fDist
std::unique_ptr< TUnuranBaseDist > fDist
Definition: TUnuran.h:264