ROOT  6.06/09
Reference Guide
unuranMultiDim.cxx
Go to the documentation of this file.
1 // test using Multi-dim Distribution object interface
2 // and compare results and CPU performances using TF3::GetRandom in case of 3D
3 // and test also case of dim = 10 and 100
4 //
5 // run within ROOT (.x unuranMultiDim.cxx+) or pass any extra parameter in the command line to get
6 // a graphics output
7 
8 
9 #include "TStopwatch.h"
10 #include "TUnuran.h"
11 #include "TUnuranMultiContDist.h"
12 
13 #include "TH3.h"
14 #include "TF3.h"
15 #include "TCanvas.h"
16 #include "TMath.h"
17 
18 #include "TRandom3.h"
19 #include "TSystem.h"
20 #include "TApplication.h"
21 #include "TError.h"
22 
23 // #include "Math/ProbFunc.h"
24 // #include "Math/DistFunc.h"
25 
26 
27 #include <iostream>
28 #include <cmath>
29 #include <cassert>
30 
31 #include <vector>
32 
33 using std::cout;
34 using std::endl;
35 
36 int n;
37 
38 bool useRandomSeed = false; // to use a random seed different every time
39 
40 double gaus3d(double *x, double *p) {
41 
42  double sigma_x = p[0];
43  double sigma_y = p[1];
44  double sigma_z = p[2];
45  double rho = p[3];
46  double u = x[0] / sigma_x ;
47  double v = x[1] / sigma_y ;
48  double w = x[2] / sigma_z ;
49  double c = 1 - rho*rho ;
50  double result = (1 / (2 * TMath::Pi() * sigma_x * sigma_y * sigma_z * sqrt(c)))
51  * exp (-(u * u - 2 * rho * u * v + v * v + w*w) / (2 * c));
52  return result;
53 }
54 double log_gaus3d(double *x, double *p) {
55  return std::log(gaus3d(x,p) );
56 }
57 
58 double gaus10d(double * x, double *) {
59  int i;
60  double tmp = 0.;
61  for (i=0; i<10; i++)
62  tmp -= x[i] * x[i];
63  return exp(tmp/2.);
64 }
65 
66 double gaus100d(double * x, double *) {
67  int i;
68  double tmp = 0.;
69  for (i=0; i<100; i++)
70  tmp -= x[i] * x[i];
71  return exp(tmp/2.);
72 }
73 
74 
75 int testUnuran(TUnuran & unr, const std::string & method, const TUnuranMultiContDist & dist, TH3 * h1, const TH3 * href = 0, const int dim = 3 ) {
76 
77 
78  assert (dim >=3);
79  // init unuran
80  bool ret = unr.Init(dist,method);
81  if (!ret) {
82  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
83  return -1;
84  }
85 
86  h1->Reset();
87 
88 
89  // test first the time
90  TStopwatch w;
91 
92  w.Start();
93  std::vector<double> x(dim);
94  for (int i = 0; i < n; ++i) {
95  unr.SampleMulti( &x[0]);
96  h1->Fill(x[0],x[1],x[2]);
97  }
98 
99  w.Stop();
100  double time = w.CpuTime()*1.E9/n;
101 
102  std::cout << "Time using Unuran " << unr.MethodName() << " \t=\t " << time << "\tns/call\t";
103  if (href != 0) {
104  double prob = href->Chi2Test(h1,"UU");
105  double ksprob = href->KolmogorovTest(h1);
106  std::cout << "\tChi2 Prob = "<< prob << "\tKS Prob = " << ksprob << std::endl;
107  // use lower value since hitro is not very precise
108  // use ks for hitro since chi2 gives too big error
109  if (unr.MethodName() == "hitro") prob = ksprob;
110  if (prob < 1.E-6 ) {
111  std::cout << "\nChi2 Test failed ! " << std::endl;
112  href->Chi2Test(h1,"UUP"); // print all chi2 test info
113  return 1;
114  }
115  }
116  else
117  std::cout << std::endl;
118 
119  return 0;
120 }
121 
122 int testGetRandom(TF3 * f, TH3 * h1, const TH3 * href = 0) {
123 
124 
125  // test first the time
126  h1->Reset();
127 
128  TStopwatch w;
129  w.Start();
130  double x[3] = {0,0,0};
131  for (int i = 0; i < n; ++i) {
132  f->GetRandom3(x[0],x[1],x[2]);
133  h1->Fill(x[0],x[1],x[2]);
134  }
135 
136  w.Stop();
137  double time = w.CpuTime()*1.E9/n;
138 
139 
140  if (href != 0) {
141  double prob = href->Chi2Test(h1,"UU");
142  double ksprob = href->KolmogorovTest(h1);
143  std::cout << "Time using TF1::GetRandom() \t=\t " << time << "\tns/call \t\tChi2 Prob = "<< prob
144  << "\tKS Prob = " << ksprob << std::endl;
145  if (prob < 1E-06) {
146  std::cout << "\tChi2 Test failed ! " << std::endl;
147  href->Chi2Test(h1,"UUP"); // print all chi2 test info
148  return 2;
149  }
150  }
151  else
152  std::cout << "Time using TF1::GetRandom() \t=\t " << time << "\tns/call\n";
153  return 0;
154 }
155 
156 
158 
159  // switch off printing of info messages from chi2 test
160  gErrorIgnoreLevel = 1001;
161 
162  // check if using a random seed
163  if (useRandomSeed) gRandom->SetSeed(0);
164 
165  gSystem->Load("libMathCore");
166  gSystem->Load("libUnuran");
167 
168  // simple test of unuran
169  n = 100000;
170 
171 
172 
173  TH3D * h1 = new TH3D("h1","UNURAN gaussian 3D distribution",50,-10,10,50,-10,10,50,-10,10);
174  TH3D * h2 = new TH3D("h2","TF1::GetRandom gaussian 3D distribution",50,-10,10,50,-10,10,50,-10,10);
175 
176 
177  TH3D * h3 = new TH3D("h3","UNURAN truncated gaussian 3D distribution",50,-1,1,50,-1,1,50,-1,1);
178  TH3D * h4 = new TH3D("h4","TF1::GetRandom truncated gaussian 3D distribution",50,-1,1,50,-1,1,50,-1,1);
179 
180 
181  TF3 * f = new TF3("g3d",gaus3d,-10,10,-10,10,-10,10,4);
182  double par[4] = {2,2,2,0.5};
183  f->SetParameters(par);
184  TF3 * flog = new TF3("logg3d",log_gaus3d,-10,10,-10,10,-10,10,4);
185  flog->SetParameters(par);
186 
187 
188 
189  std::cout << "Test using an undefined domain :\n\n";
190 
191  std::cout << "Function Points used in GetRandom: [ " << f->GetNpx() << " , "
192  << f->GetNpy() << " , " << f->GetNpz() << " ]" << std::endl;
193  testGetRandom(f,h1);
194 
195  // test first with getrandom
196  // need to have a larger value to get good quality
197  int np = 100;
198  f->SetNpx(np); f->SetNpy(np); f->SetNpz(np);
199  std::cout << "Function Points used in GetRandom: [ " << f->GetNpx() << " , "
200  << f->GetNpy() << " , " << f->GetNpz() << " ]" << std::endl;
201 
202  testGetRandom(f,h2,h1);
203 
204  *h1 = *h2;
205  np = 200;
206  f->SetNpx(np); f->SetNpy(np); f->SetNpz(np);
207  std::cout << "Function Points used in GetRandom: [ " << f->GetNpx() << " , "
208  << f->GetNpy() << " , " << f->GetNpz() << " ]" << std::endl;
209 
210  testGetRandom(f,h2,h1);
211 
212 
213  // create multi-dim distribution
215 
216 
217  TUnuran unr(gRandom,2); // 2 is debug level
218 
219  int iret = 0;
220  TH3 * href = new TH3D(*h2);
221 
222  //vnrou method (requires only pdf)
223  std::string method = "vnrou";
224  iret |= testUnuran(unr, method, dist, h1, href);
225 
226 
227  //hitro method (requires only pdf)
228  method = "hitro";
229  iret |= testUnuran(unr, method, dist, h1, href);
230 
231  //gibbs requires log of pdf and derivative
232 #ifdef USE_GIBBS
233  method = "gibbs";
234  // need to create a new multi-dim distribution with log of pdf
235  TUnuranMultiContDist logdist(flog,0,true);
236  iret |= testUnuran(unr, method, logdist, h1, href);
237 #endif
238 
239  // test setting the mode
240  std::cout << "\nTest setting the mode in Unuran distribution:\n\n";
241  double m[3] = {0,0,0};
242  dist.SetMode(m);
243 
244  method = "vnrou";
245  iret |= testUnuran(unr, method, dist, h1, href);
246 
247  method = "hitro";
248  iret |= testUnuran(unr, method, dist, h1, href);
249 
250 #ifdef USE_GIBBS
251  method = "gibbs";
252  logdist.SetMode(m);
253  iret |= testUnuran(unr, method, logdist, h1, href);
254 #endif
255 
256  std::cout << "\nTest truncated distribution:\n\n";
257 
258  double xmin[3] = { -1, -1, -1 };
259  double xmax[3] = { 1, 1, 1 };
260 
261  f->SetRange(xmin[0],xmin[1],xmin[2],xmax[0],xmax[1],xmax[2]);
262  // change function domain (not yet implemented in unuran for multidim functions)
263  dist.SetDomain(xmin,xmax);
264 
265  np = 30;
266  f->SetNpx(np); f->SetNpy(np); f->SetNpz(np);
267  std::cout << "Function Points used in GetRandom: [ " << f->GetNpx() << " , "
268  << f->GetNpy() << " , " << f->GetNpz() << " ]" << std::endl;
269  testGetRandom(f,h3);
270 
271  href = h3;
272  np = 100;
273  f->SetNpx(np); f->SetNpy(np); f->SetNpz(np);
274  std::cout << "Function Points used in GetRandom: [ " << f->GetNpx() << " , "
275  << f->GetNpy() << " , " << f->GetNpz() << " ]" << std::endl;
276  testGetRandom(f,h4,href);
277  href = h4;
278 
279  method = "vnrou";
280  iret |= testUnuran(unr, method, dist, h3, href);
281  method = "hitro";
282  iret |= testUnuran(unr, method, dist, h3, href);
283 
284 
285 #ifdef USE_GIBBS
286  method = "gibbs";
287  logdist.SetDomain(xmin,xmax);
288  iret |= testUnuran(unr, method, logdist, h3, href);
289 #endif
290 
291 
292  TCanvas * c1 = new TCanvas("c1_unuranMulti","Multidimensional distribution",10,10,900,900);
293  c1->Divide(2,2);
294 
295  c1->cd(1); h1->Draw();
296  c1->cd(2); h2->Draw();
297  c1->cd(3); h3->Draw();
298  c1->cd(4); h4->Draw();
299 
300  // make a ref histo for 10 dim using first 3 dim
301  c1->Update();
302 
303 
304 
305  TH3D * hrefN = new TH3D("hrefN","UNURAN gaussian ref N-Dim distribution (first 3 dim)",30,-3,3,30,-3,3,30,-3,3);
306  TH3D * h10v = new TH3D("h10v","UNURAN gaussian N-Dim distribution (first 3 dim)",30,-3,3,30,-3,3,30,-3,3);
307  TH3D * h10h = new TH3D("h10h","UNURAN gaussian N-Dim distribution (first 3 dim)",30,-3,3,30,-3,3,30,-3,3);
308  TH3D * h100 = new TH3D("h100","UNURAN gaussian N-Dim distribution (first 3 dim)",30,-3,3,30,-3,3,30,-3,3);
309 
310  int scale = 5;
311 
312  double par2[4] = {1,1,1,0.};
313  f->SetParameters(par2);
314  TUnuranMultiContDist dist3(f);
315  method = "vnrou";
316  n/= scale;
317  iret |= testUnuran(unr, method, dist3, hrefN);
318 
319  std::cout << "\nTest 10 dimension: (be patient , it takes time....)\n\n";
320 
321  TF1 * f10 = new TF1("g10d",gaus10d,-10,10,0);
322  TUnuranMultiContDist dist10(f10,10);
323 
324 
325  TCanvas * c2 = new TCanvas("c2_unuranMulti","Multidimensional distribution",100,10,900,900);
326  c2->Divide(2,2);
327 
328  c2->cd(1); hrefN->Draw();
329 
330  //n/= scale;
331  method = "vnrou";
332  iret |= testUnuran(unr, method, dist10, h10v, hrefN,10);
333  c2->cd(2); h10v->Draw();
334  c2->Update();
335 
336  //n*=scale;
337  method = "hitro;thinning=20";
338  iret |= testUnuran(unr, method, dist10, h10h, hrefN,10);
339  c2->cd(3); h10h->Draw();
340  c2->Update();
341 
342 
343  // 100 dim
344  std::cout << "\nTest 100 dimension: ( be patient , it takes time....)\n\n";
345  TF1 * f100 = new TF1("g100d",gaus100d,-10,10,0);
346  TUnuranMultiContDist dist100(f100,100);
347 
348  // scale = 5;
349  // n/=scale;
350  std::cout << "number of events to be generated = " << n << std::endl;
351  method = "hitro;thinning=200";
352  iret |= testUnuran(unr, method, dist100, h100, hrefN,100);
353 // n/= 100;
354 // method = "vnrou";
355 // iret |= testUnuran(unr, method, dist100, hN, hrefN,100);
356 
357 
358  c2->cd(4); h100->Draw();
359  c2->Update();
360 
361 
362  if (iret != 0)
363  std::cerr <<"\n\nUnuRan MultiDim Continous Distribution Test:\t Failed !!!!!!!\n" << std::endl;
364  else
365  std::cerr << "\n\nUnuRan MultiDim Continous Distribution Test:\t OK\n" << std::endl;
366  return iret;
367 
368 }
369 
370 #ifndef __CINT__
371 int main(int argc, char **argv)
372 {
373  int iret = 0;
374  if (argc > 1) {
375  TApplication theApp("App",&argc,argv);
376  iret = unuranMultiDim();
377  theApp.Run();
378  }
379  else
380  iret = unuranMultiDim();
381 
382  return iret;
383 }
384 #endif
385 
double par[1]
Definition: unuranDistr.cxx:38
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:439
float xmin
Definition: THbookFile.cxx:93
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
Definition: TF1.cxx:3139
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
R__EXTERN Int_t gErrorIgnoreLevel
Definition: TError.h:107
TCanvas * c1
Definition: legend1.C:2
#define assert(cond)
Definition: unittest.h:542
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:659
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition: TSystem.cxx:1818
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:123
int testGetRandom(TF3 *f, TH3 *h1, const TH3 *href=0)
bool SampleMulti(double *x)
Sample multidimensional distributions User is responsible for having previously correctly initialized...
Definition: TUnuran.cxx:396
virtual void SetSeed(UInt_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:568
double gaus100d(double *x, double *)
Int_t GetNpz() const
Definition: TF3.h:95
double sqrt(double)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
Double_t x[n]
Definition: legend1.C:17
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
TString flog
Definition: pq2main.cxx:37
TH1F * h1
Definition: legend1.C:5
The 3-D histogram classes derived from the 1-D histogram classes.
Definition: TH3.h:35
double gaus3d(double *x, double *p)
R__EXTERN TSystem * gSystem
Definition: TSystem.h:549
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2878
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:309
SVector< double, 2 > v
Definition: Dict.h:5
double log_gaus3d(double *x, double *p)
A 3-Dim function with parameters.
Definition: TF3.h:30
TMarker * m
Definition: textangle.C:8
virtual void SetNpy(Int_t npy=100)
Set the number of points used to draw the function.
Definition: TF2.cxx:907
virtual void GetRandom3(Double_t &xrandom, Double_t &yrandom, Double_t &zrandom)
Return 3 random numbers following this function shape.
Definition: TF3.cxx:322
Double_t E()
Definition: TMath.h:54
const std::string & MethodName() const
used Unuran method
Definition: TUnuran.h:239
void SetDomain(const double *xmin, const double *xmax)
set the domain of the distribution giving an array of minimum and maximum values By default otherwise...
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH3.cxx:275
float xmax
Definition: THbookFile.cxx:93
TUnuranMultiContDist class describing multi dimensional continuous distributions. ...
int unuranMultiDim()
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
int testUnuran(TUnuran &unr, const std::string &method, const TUnuranMultiContDist &dist, TH3 *h1, const TH3 *href=0, const int dim=3)
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH3.cxx:3476
Double_t Pi()
Definition: TMath.h:44
THist< 3, double > TH3D
Definition: THist.h:326
virtual Int_t GetNpx() const
Definition: TF1.h:352
The Canvas class.
Definition: TCanvas.h:48
return c2
Definition: legend2.C:14
TUnuran class.
Definition: TUnuran.h:81
double f(double x)
bool useRandomSeed
double gaus10d(double *x, double *)
int main(int argc, char **argv)
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1077
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition: TF3.h:145
bool Init(const std::string &distr, const std::string &method)
initialize with Unuran string interface
Definition: TUnuran.cxx:79
1-Dim function class
Definition: TF1.h:149
void SetMode(const double *x)
set the mode of the distribution (coordinates of the distribution maximum values) ...
Int_t GetNpy() const
Definition: TF2.h:104
int n
double result[121]
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:45
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2179
double exp(double)
double log(double)
virtual void SetNpz(Int_t npz=30)
Set the number of points used to draw the function.
Definition: TF3.cxx:664
Stopwatch class.
Definition: TStopwatch.h:30