Logo ROOT   6.10/09
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 #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 
36 TUnuran::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 
79 bool 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 
93 bool 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 
109 bool 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 
125 bool 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 
139 bool 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 
396 bool TUnuran::SampleMulti(double * x)
397 {
398  // sample multidimensional distribution
399  if (fGen == 0) return false;
400  unur_sample_vec(fGen,x);
401  return true;
402 }
403 
404 void TUnuran::SetSeed(unsigned int seed) {
405  return fRng->SetSeed(seed);
406 }
407 
408 bool 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 
423 bool 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 
438 bool 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 
453 bool 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 
bool HasPdfArea() const
check if distribution has a pre-computed area below the Pdf
UNUR_URNG * fUrng
Definition: TUnuran.h:263
double par[1]
Definition: unuranDistr.cxx:38
static double Pdpdf(const double *x, int coord, UNUR_DISTR *dist)
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
float xmin
Definition: THbookFile.cxx:93
bool IsLogPdf() const
flag to control if given function represent the log of a pdf
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
double Sample()
Sample 1D distribution User is responsible for having previously correctly initialized with TUnuran::...
Definition: TUnuran.cxx:389
bool ReInitDiscrDist(unsigned int npar, double *params)
Reinitialize UNURAN by changing the distribution parameters but mantaining same distribution and meth...
Definition: TUnuran.cxx:453
static double Cdf(int x, const UNUR_DISTR *dist)
evaluate the cumulative function
TUnuran & operator=(const TUnuran &rhs)
Assignment operator.
Definition: TUnuran.cxx:72
const double * GetMode() const
get the mode (vector of coordinate positions of the maxima of the distribution) If a mode has not def...
static double Pdf(const double *x, UNUR_DISTR *dist)
evaluate the probality density function
bool GetDomain(double &xmin, double &xmax) const
check if distribution has a defined domain and return in case its domain
const double * GetLowerDomain() const
get the distribution lower domain values.
unsigned int NDim() const
get number of dimension of the distribution
bool GetDomain(int &xmin, int &xmax) const
check if distribution has domain and return in case its domain
double UpperBin() const
upper value of binned data (return 0 for unbinned data)
int nbins[3]
int SampleDiscr()
Sample discrete distributions User is responsible for having previously correctly initialized with TU...
Definition: TUnuran.cxx:382
UnuranRng class for interface ROOT random generators to Unuran.
Definition: UnuranRng.h:22
#define NULL
Definition: RtypesCore.h:88
bool SampleMulti(double *x)
Sample multidimensional distributions User is responsible for having previously correctly initialized...
Definition: TUnuran.cxx:396
static double Dpdf(double x, const UNUR_DISTR *dist)
evaluate the derivative of the pdf
bool SetLogLevel(unsigned int iflag=1)
set log level
Definition: TUnuran.cxx:408
static int Dpdf(double *grad, const double *x, UNUR_DISTR *dist)
Double_t x[n]
Definition: legend1.C:17
const std::vector< double > & Data() const
Return reference to data vector (unbinned or binned data)
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
virtual TUnuranDiscrDist * Clone() const
Clone (required by base class)
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:568
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
void SetSeed(unsigned int seed)
set the seed for the random number generator
Definition: TUnuran.cxx:404
virtual TUnuranMultiContDist * Clone() const
Clone (required by base class)
TRandom * fRng
Definition: TUnuran.h:265
TUnuranDiscrDist class for one dimensional discrete distribution.
bool SetMultiDistribution(const TUnuranMultiContDist &dist)
Definition: TUnuran.cxx:219
bool HasCdf() const
check if a cdf function is provided for the distribution
TUnuranEmpDist class for describing empiral distributions.
bool IsLogPdf() const
flag to control if given function represent the log of a pdf
void Error(const char *location, const char *msgfmt,...)
unsigned int NDim() const
Number of data dimensions.
double LowerBin() const
Min value of binned data (return 0 for unbinned data)
bool SetContDistribution(const TUnuranContDist &dist)
Definition: TUnuran.cxx:173
~TUnuran()
Destructor.
Definition: TUnuran.cxx:57
static double Pmf(int x, const UNUR_DISTR *dist)
evaluate the probality mesh function
bool HasProbSum() const
flag to control if distribution provides the total area of the probability function ...
double ProbSum() const
return area of the pdf
const std::vector< double > & ProbVec() const
retrieve a reference to the vector of the probabilities : Prob(i) If the distribution is defined from...
TRandom2 r(17)
bool HasMode() const
flag to control if distribution provides the mode
static double Pdf(double x, const UNUR_DISTR *dist)
evaluate the probality density function
float xmax
Definition: THbookFile.cxx:93
TUnuranMultiContDist class describing multi dimensional continuous distributions. ...
static double Cdf(double x, const UNUR_DISTR *dist)
evaluate the Cumulative distribution function, integral of the pdf
std::string fMethod
Definition: TUnuran.h:266
void Warning(const char *location, const char *msgfmt,...)
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
int Mode() const
get the mode (x location of function maximum)
TUnuran(TRandom *r=0, unsigned int log=0)
Constructor with a generator instance and given level of log output.
Definition: TUnuran.cxx:36
bool HasCdf() const
flag to control if distribution provides also a Cdf
bool SetMethodAndInit()
change the method and initialize Unuran with the previously given distribution
Definition: TUnuran.cxx:351
TUnuran class.
Definition: TUnuran.h:79
bool SetDiscreteDistribution(const TUnuranDiscrDist &dist)
Definition: TUnuran.cxx:304
double PdfArea() const
return area below the pdf
UNUR_GEN * fGen
Definition: TUnuran.h:261
bool SetRandomGenerator()
Definition: TUnuran.cxx:156
const double * GetUpperDomain() const
get the distribution upper domain values.
UNUR_DISTR * fUdistr
Definition: TUnuran.h:262
TUnuranContDist class describing one dimensional continuous distribution.
bool IsBinned() const
Flag to control if data are binned.
bool Init(const std::string &distr, const std::string &method)
initialize with Unuran string interface
Definition: TUnuran.cxx:79
bool HasMode() const
check if distribution has a pre-computed mode
std::unique_ptr< TUnuranBaseDist > fDist
Definition: TUnuran.h:264
bool SetEmpiricalDistribution(const TUnuranEmpDist &dist)
Definition: TUnuran.cxx:263
virtual TUnuranContDist * Clone() const
Clone (required by base class)
TUnuranEmpDist * Clone() const
Clone (required by base class)
double Mode() const
return the mode (x location of maximum of the pdf)
const Int_t n
Definition: legend1.C:16