Logo ROOT   6.10/00
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
unuranSimple.cxx
Go to the documentation of this file.
1 // test unuran using the string interface to generate numbers according to the normal distributions
2 // compare CPU performancecwith TRandom::Gaus and opitonally GSL (using MathMore ) and CLHEP for
3 // generating normal distributed random numbers
4 //
5 // run within ROOT (.x unuranSimple.cxx+) or pass any extra parameter in the command line to get
6 // a graphics output
7 //
8 #include "TStopwatch.h"
9 #include "TUnuran.h"
10 
11 #include "TH1.h"
12 #include "TF1.h"
13 
14 #include "TRandom.h"
15 #include "TRandom1.h"
16 #include "TRandom2.h"
17 #include "TSystem.h"
18 #include "TApplication.h"
19 #include "TCanvas.h"
20 //#include "TRint.h"
21 #include "TVirtualFitter.h"
22 #include "TFitter.h"
23 #include "Math/DistFunc.h"
24 
25 #include <iostream>
26 
27 #ifdef HAVE_MATHMORE
28 #include "Math/Random.h"
29 #include "Math/GSLRndmEngines.h"
30 #endif
31 
32 
33 
34 #ifdef HAVE_CLHEP
35 #include "CLHEP/Random/RandFlat.h"
36 #include "CLHEP/Random/RandGauss.h"
37 #include "CLHEP/Random/MTwistEngine.h"
38 #include "CLHEP/Random/JamesRandom.h"
39 #include "CLHEP/Random/RanluxEngine.h"
40 #include "CLHEP/Random/Ranlux64Engine.h"
41 #include "CLHEP/Random/RanecuEngine.h"
42 #include "CLHEP/Random/Hurd160Engine.h"
43 #include "CLHEP/Random/Hurd288Engine.h"
44 #include "CLHEP/Random/RanshiEngine.h"
45 #include "CLHEP/Random/DualRand.h"
46 #include "CLHEP/Random/TripleRand.h"
47 
48 using namespace CLHEP;
49 #endif
50 
51 
52 using std::cout;
53 using std::endl;
54 
55 int unuranSimple( ) {
56 
57  // simple test of unuran
58 
59 
60  std::cout << "Test Generation of Gaussian Numbers \n\n";
61 
62  TUnuran unr;
63  if (! unr.Init( "normal()", "method=arou") ) {
64  std::cout << "Error initializing unuran" << std::endl;
65  return -1;
66  }
67 
68  // default is 10**7 but one should use 10**8 for serious timing
69  int n = 10000000;
70 
71  TStopwatch w;
72  double time;
73  w.Start();
74 
75  for (int i = 0; i < n; ++i)
76  unr.Sample();
77 
78  w.Stop();
79  time = w.CpuTime()*1.E9/n;
80  std::cout << "Time using Unuran method arou =\t " << time << "\tns/call" << std::endl;
81 
82 
83  if (! unr.Init( "normal()", "method=tdr") ) {
84  std::cout << "Error initializing unuran" << std::endl;
85  return -1;
86  }
87  w.Start();
88  for (int i = 0; i < n; ++i)
89  unr.Sample();
90 
91  w.Stop();
92  time = w.CpuTime()*1.E9/n;
93  std::cout << "Time using Unuran method tdr =\t " << time << "\tns/call" << std::endl;
94 
95  if (! unr.Init( "normal()", "method=hinv") ) {
96  std::cout << "Error initializing unuran" << std::endl;
97  return -1;
98  }
99  w.Start();
100  for (int i = 0; i < n; ++i)
101  unr.Sample();
102 
103  w.Stop();
104  time = w.CpuTime()*1.E9/n;
105  std::cout << "Time using Unuran method hinv =\t " << time << "\tns/call" << std::endl;
106 
107  w.Start();
108  for (int i = 0; i < n; ++i)
109  gRandom->Gaus(0,1);
110  w.Stop();
111  time = w.CpuTime()*1.E9/n;
112  std::cout << "Time using TRandom::Gaus =\t " << time << "\tns/call" << std::endl;
113 
114  // using Rannor
115  w.Start();
116  double x1,x2;
117  for (int i = 0; i < n/2; ++i)
118  gRandom->Rannor(x1,x2);
119  w.Stop();
120  time = w.CpuTime()*1.E9/n;
121  std::cout << "Time using TRandom::Rannor =\t " << time << "\tns/call" << std::endl;
122 
123 #ifdef HAVE_MATHMORE
124  // using GSL Ziggurat
125  ROOT::Math::Random<ROOT::Math::GSLRngMT> rgsl;
126  w.Start();
127  for (int i = 0; i < n; ++i)
128  rgsl.Gaus(0,1);
129  w.Stop();
130  time = w.CpuTime()*1.E9/n;
131  std::cout << "Time using GSL::Gaus =\t\t " << time << "\tns/call" << std::endl;
132 
133  // using GSL BoxMuller method
134  w.Start();
135  for (int i = 0; i < n; ++i)
136  rgsl.GausBM(0,1);
137  w.Stop();
138  time = w.CpuTime()*1.E9/n;
139  std::cout << "Time using GSL::GausBM = \t " << time << "\tns/call" << std::endl;
140 
141  // using GSL Ratio method
142  w.Start();
143  for (int i = 0; i < n; ++i)
144  rgsl.GausR(0,1);
145  w.Stop();
146  time = w.CpuTime()*1.E9/n;
147  std::cout << "Time using GSL::GausR =\t " << time << "\tns/call" << std::endl;
148 #endif
149 
150 
151  // run unuran standard Normal methods 2 (Polar from Marsaglia)
152 
153  if (! unr.Init( "normal()", "method=cstd;variant=2") ) {
154  std::cout << "Error initializing unuran" << std::endl;
155  return -1;
156  }
157  w.Start();
158  for (int i = 0; i < n; ++i)
159  unr.Sample();
160 
161  w.Stop();
162  time = w.CpuTime()*1.E9/n;
163  std::cout << "Time using Unuran GausPolarR =\t " << time << "\tns/call" << std::endl;
164 
165  // run unuran standard Normal method 3 (KM)
166 
167  if (! unr.Init( "normal()", "method=cstd;variant=3") ) {
168  std::cout << "Error initializing unuran" << std::endl;
169  return -1;
170  }
171  w.Start();
172  for (int i = 0; i < n; ++i)
173  unr.Sample();
174 
175  w.Stop();
176  time = w.CpuTime()*1.E9/n;
177  std::cout << "Time using Unuran Gaus K-R =\t " << time << "\tns/call" << std::endl;
178 
179  if (! unr.Init( "normal()", "method=cstd;variant=6") ) {
180  std::cout << "Error initializing unuran" << std::endl;
181  return -1;
182  }
183  w.Start();
184  for (int i = 0; i < n; ++i)
185  unr.Sample();
186 
187  w.Stop();
188  time = w.CpuTime()*1.E9/n;
189  std::cout << "Time using Unuran Gaus exp6 =\t " << time << "\tns/call" << std::endl;
190 
191 
192 // w.Start();
193 // for (int i = 0; i < n; ++i)
194 // rgsl.GausRatio(0,1);
195 // w.Stop();
196 // time = w.CpuTime()*1.E9/n;
197 // std::cout << "Time using GSL::GausRatio =\t\t " << time << "\tns/call" << std::endl;
198 
199 #ifdef HAVE_CLHEP
200  MTwistEngine eng(111);
201  RandGauss rg(eng);
202  w.Start();
203  for (int i = 0; i < n; ++i)
204  rg(0,1);
205  w.Stop();
206  time = w.CpuTime()*1.E9/n;
207  std::cout << "Time using CLHEP::Gaus =\t " << time << "\tns/call" << std::endl;
208 #endif
209 
210  std::cout << "\nTest uniform generator\n" << std::endl;
211  w.Start();
212  for (int i = 0; i < n; ++i)
213  gRandom->Rndm();
214  w.Stop();
215  time = w.CpuTime()*1.E9/n;
216  std::cout << "Time using gRandom::Rndm =\t " << time << "\tns/call" << std::endl;
217 
218  TRandom1 r1;
219  w.Start();
220  for (int i = 0; i < n; ++i)
221  r1.Rndm();
222  w.Stop();
223  time = w.CpuTime()*1.E9/n;
224  std::cout << "Time using TRandom1::Rndm =\t " << time << "\tns/call" << std::endl;
225 
226  TRandom2 r2;
227  w.Start();
228  for (int i = 0; i < n; ++i)
229  r2.Rndm();
230  w.Stop();
231  time = w.CpuTime()*1.E9/n;
232  std::cout << "Time using TRandom2::Rndm =\t " << time << "\tns/call" << std::endl;
233 
234 #ifdef HAVE_CLHEP
235  RandFlat rf(eng);
236  w.Start();
237  for (int i = 0; i < n; ++i)
238  eng.flat(); // use directly the engine (faster!)
239  w.Stop();
240  time = w.CpuTime()*1.E9/n;
241  std::cout << "Time using CLHEP::MT =\t\t " << time << "\tns/call" << std::endl;
242 
243  {
244  RanecuEngine e;
245  RandFlat rf2(e);
246  w.Start();
247  for (int i = 0; i < n; ++i)
248  e.flat();
249  w.Stop();
250  time = w.CpuTime()*1.E9/n;
251  std::cout << "Time using CLHEP::Ranecu =\t " << time << "\tns/call" << std::endl;
252  }
253  {
254  Hurd160Engine e;
255  RandFlat rf2(e);
256  w.Start();
257  for (int i = 0; i < n; ++i)
258  e.flat();
259  w.Stop();
260  time = w.CpuTime()*1.E9/n;
261  std::cout << "Time using CLHEP::Hard160 =\t " << time << "\tns/call" << std::endl;
262  }
263  {
264  Hurd288Engine e;
265  RandFlat rf2(e);
266  w.Start();
267  for (int i = 0; i < n; ++i)
268  e.flat();
269  //rf2();
270  w.Stop();
271  time = w.CpuTime()*1.E9/n;
272  std::cout << "Time using CLHEP::Hard288 =\t " << time << "\tns/call" << std::endl;
273  }
274  {
275  DualRand e;
276  RandFlat rf2(e);
277  w.Start();
278  for (int i = 0; i < n; ++i)
279  e.flat();
280  //rf2();
281  w.Stop();
282  time = w.CpuTime()*1.E9/n;
283  std::cout << "Time using CLHEP::DualRand =\t " << time << "\tns/call" << std::endl;
284  }
285  {
286  TripleRand e;
287  RandFlat rf2(e);
288  w.Start();
289  for (int i = 0; i < n; ++i)
290  e.flat();
291  //rf2();
292  w.Stop();
293  time = w.CpuTime()*1.E9/n;
294  std::cout << "Time using CLHEP::TripleRand =\t " << time << "\tns/call" << std::endl;
295  }
296  {
297  RanshiEngine e;
298  RandFlat rf2(e);
299  w.Start();
300  for (int i = 0; i < n; ++i)
301  //rf2();
302  e.flat();
303  w.Stop();
304  time = w.CpuTime()*1.E9/n;
305  std::cout << "Time using CLHEP::Runshi =\t " << time << "\tns/call" << std::endl;
306  }
307  {
308  RanluxEngine e;
309  RandFlat rf2(e);
310  w.Start();
311  for (int i = 0; i < n; ++i)
312  //rf2();
313  e.flat();
314  w.Stop();
315  time = w.CpuTime()*1.E9/n;
316  std::cout << "Time using CLHEP::RunLux =\t " << time << "\tns/call" << std::endl;
317  }
318  {
319  Ranlux64Engine e;
320  RandFlat rf2(e);
321  w.Start();
322  for (int i = 0; i < n; ++i)
323  e.flat();
324  w.Stop();
325  time = w.CpuTime()*1.E9/n;
326  std::cout << "Time using CLHEP::RunLux64 =\t " << time << "\tns/call" << std::endl;
327  }
328  {
329  HepJamesRandom e;
330  RandFlat rf2(e);
331  w.Start();
332  for (int i = 0; i < n; ++i)
333  //rf2();
334  e.flat();
335  w.Stop();
336  time = w.CpuTime()*1.E9/n;
337  std::cout << "Time using CLHEP::HepJames =\t " << time << "\tns/call" << std::endl;
338  }
339 
340 #endif
341 
342 
343  // test the quality by looking at the cdf
344  std::cout <<"\n\nTest quality of Unuran arou" << std::endl;
345  if (! unr.Init( "normal()", "method=arou") ) {
346  std::cout << "Error initializing unuran" << std::endl;
347  return -1;
348  }
349 
350  TH1D * h1 = new TH1D("h1","cdf on the data ",1000,0,1);
351  for (int i = 0; i < 1000000; ++i) {
352  double x = unr.Sample();
353  h1->Fill( ROOT::Math::normal_cdf( x , 1.0) );
354  }
355 
356  new TCanvas("c1_unuranGaus","unuran Gaus CDF");
357 
358  h1->Fit("pol0","Q");
359  TF1 * f = h1->GetFunction("pol0");
360 
361  std::cout << "CDF Uniform Fit: chi2 = " << f->GetChisquare() << " ndf = " << f->GetNDF() << std::endl;
362  std::cout << "Fit Prob = " << f->GetProb() << std::endl;
363  h1->Draw("E");
364 
365  if (f->GetProb() < 1.E-4) {
366  std::cerr << "\nERROR: UnuranSimple Test:\t Failed !!!!";
367  return -1;
368  }
369  std::cerr << "\nUnuranSimple Test:\t OK !" << std::endl;
370 
371  return 0;
372 
373 }
374 
375 #ifndef __CINT__
376 int main(int argc, char **argv)
377 {
378  int iret = 0;
379  if (argc > 1) {
380  TApplication theApp("App",&argc,argv);
381  iret = unuranSimple();
382  theApp.Run();
383  }
384  else
385  iret = unuranSimple();
386 
387  return iret;
388 }
389 #endif
390 
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3126
TRandom::Rannor
virtual void Rannor(Float_t &a, Float_t &b)
Return 2 numbers distributed following a gaussian with mean=0 and sigma=1.
Definition: TRandom.cxx:460
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
TUnuran::Sample
double Sample()
Sample 1D distribution User is responsible for having previously correctly initialized with TUnuran::...
Definition: TUnuran.cxx:389
TRandom::Gaus
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:235
TRandom2
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:27
TRandom1::Rndm
virtual Double_t Rndm()
return a random number in ]0,1]
Definition: TRandom1.cxx:360
TStopwatch::CpuTime
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
TF1::GetNDF
virtual Int_t GetNDF() const
Return the number of degrees of freedom in the fit the fNDF parameter has been previously computed du...
Definition: TF1.cxx:1616
x2
static const double x2[5]
Definition: RooGaussKronrodIntegrator1D.cxx:345
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
x
Double_t x[n]
Definition: legend1.C:17
TRandom1
The Ranlux Random number generator class.
Definition: TRandom1.h:27
TF1::GetChisquare
Double_t GetChisquare() const
Definition: TF1.h:398
TApplication::Run
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
Definition: TApplication.cxx:1139
ROOT::Math::normal_cdf
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
Definition: ProbFuncMathCore.cxx:234
ROOT::Math::Random::Gaus
double Gaus(double mean=0, double sigma=1)
Definition: Random.h:107
h1
TH1F * h1
Definition: legend1.C:5
TRandom::Rndm
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:512
unuranSimple
int unuranSimple()
Definition: unuranSimple.cxx:55
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2851
r1
unsigned int r1[N_CITIES]
Definition: simanTSP.cxx:321
TF1::GetProb
virtual Double_t GetProb() const
Return the fit probability.
Definition: TF1.cxx:1681
gRandom
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
TH1D
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:594
TCanvas
The Canvas class.
Definition: TCanvas.h:31
TUnuran
TUnuran class.
Definition: TUnuran.h:79
TRandom2::Rndm
virtual Double_t Rndm()
TausWorth generator from L&#39;Ecuyer, uses as seed 3x32bits integers Use a mask of 0xffffffffUL to make ...
Definition: TRandom2.cxx:58
x1
static const double x1[5]
Definition: RooGaussKronrodIntegrator1D.cxx:327
f
double f(double x)
Definition: testIntegration.cxx:12
e
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
TUnuran::Init
bool Init(const std::string &distr, const std::string &method)
initialize with Unuran string interface
Definition: TUnuran.cxx:79
TF1
1-Dim function class
Definition: TF1.h:150
ROOT::Experimental::TH1D
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:310
TApplication
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:39
TH1::Fit
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3564
TH1::GetFunction
virtual TF1 * GetFunction(const char *name) const
Return pointer to function with name.
Definition: TH1.cxx:8163
n
const Int_t n
Definition: legend1.C:16
main
int main(int argc, char **argv)
Definition: UnitTesting.cxx:426
r2
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
ROOT::Math::Random
Documentation for the Random class.
Definition: Random.h:39
TStopwatch
Stopwatch class.
Definition: TStopwatch.h:28