Logo ROOT   6.19/01
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 
UNUR_URNG * fUrng
Definition: TUnuran.h:263
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 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
static double Pdf(const double *x, UNUR_DISTR *dist)
evaluate the probality density function
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
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
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:597
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
static constexpr double s
TUnuranEmpDist class for describing empiral distributions.
void Error(const char *location, const char *msgfmt,...)
unsigned int NDim() const
Number of data dimensions.
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
ROOT::R::TRInterface & r
Definition: Object.C:4
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
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
TUnuran class.
Definition: TUnuran.h:79
bool SetDiscreteDistribution(const TUnuranDiscrDist &dist)
Definition: TUnuran.cxx:304
UNUR_GEN * fGen
Definition: TUnuran.h:261
bool SetRandomGenerator()
Definition: TUnuran.cxx:156
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
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)
const Int_t n
Definition: legend1.C:16