Example macro to show unuran capabilities The results are compared with what is obtained using TRandom or TF1::GetRandom The macro is divided in 3 parts:
#include <iostream>
#include <cassert>
using std::cout;
using std::endl;
#define NGEN 1000000
int izone = 0;
void testStringAPI() {
TH1D * h1 =
new TH1D(
"h1G",
"gaussian distribution from Unuran",100,-10,10);
TH1D * h2 =
new TH1D(
"h2G",
"gaussian distribution from TRandom",100,-10,10);
cout << "\nTest using UNURAN string API \n\n";
if (! unr.
Init(
"normal()",
"method=arou") ) {
cout << "Error initializing unuran" << endl;
return;
}
int n = NGEN;
for (
int i = 0; i <
n; ++i) {
}
cout <<
"Time using Unuran method " << unr.
MethodName() <<
"\t=\t " << w.
CpuTime() << endl;
for (
int i = 0; i <
n; ++i) {
double x = gRandom->Gaus(0,1);
}
cout <<
"Time using TRandom::Gaus \t=\t " << w.
CpuTime() << endl;
assert(c1 != 0);
}
double distr(double *x, double *p) {
}
double cdf(double *x, double *p) {
}
void testDistr1D() {
cout << "\nTest 1D Continous distributions\n\n";
TH1D * h1 =
new TH1D(
"h1BW",
"Breit-Wigner distribution from Unuran",100,-10,10);
TH1D * h2 =
new TH1D(
"h2BW",
"Breit-Wigner distribution from GetRandom",100,-10,10);
TF1 * f =
new TF1(
"distrFunc",distr,-10,10,2);
double par[2] = {1,0};
TF1 * fc =
new TF1(
"cdfFunc",cdf,-10,10,2);
int logLevel = 2;
std::string method = "tdr";
cout << "Error initializing unuran" << endl;
return;
}
int n = NGEN;
for (
int i = 0; i <
n; ++i) {
}
cout <<
"Time using Unuran method " << unr.
MethodName() <<
"\t=\t " << w.
CpuTime() << endl;
for (
int i = 0; i <
n; ++i) {
}
cout <<
"Time using TF1::GetRandom() \t=\t " << w.
CpuTime() << endl;
std::cout << " chi2 test of UNURAN vs GetRandom generated histograms: " << std::endl;
}
double gaus3d(double *x, double *p) {
double sigma_x = p[0];
double sigma_y = p[1];
double sigma_z = p[2];
double rho = p[2];
double u = x[0] / sigma_x ;
double v = x[1] / sigma_y ;
double w = x[2] / sigma_z ;
double c = 1 - rho*rho ;
double result = (1 / (2 *
TMath::Pi() * sigma_x * sigma_y * sigma_z *
sqrt(c)))
*
exp (-(u * u - 2 * rho * u * v + v * v + w*w) / (2 * c));
return result;
}
void testDistrMultiDim() {
cout << "\nTest Multidimensional distributions\n\n";
TH3D * h1 =
new TH3D(
"h13D",
"gaussian 3D distribution from Unuran",50,-10,10,50,-10,10,50,-10,10);
TH3D * h2 =
new TH3D(
"h23D",
"gaussian 3D distribution from GetRandom",50,-10,10,50,-10,10,50,-10,10);
TF3 * f =
new TF3(
"g3d",gaus3d,-10,10,-10,10,-10,10,3);
double par[3] = {2,2,0.5};
std::string method = "hitro";
cout << "Error initializing unuran" << endl;
return;
}
double x[3];
for (int i = 0; i < NGEN; ++i) {
h1->
Fill(x[0],x[1],x[2]);
}
cout <<
"Time using Unuran method " << unr.
MethodName() <<
"\t=\t\t " << w.
CpuTime() << endl;
assert(c1 != 0);
int np = 200;
for (int i = 0; i < NGEN; ++i) {
h2->
Fill(x[0],x[1],x[2]);
}
cout <<
"Time using TF1::GetRandom \t\t=\t " << w.
CpuTime() << endl;
std::cout << " chi2 test of UNURAN vs GetRandom generated histograms: " << std::endl;
}
double poisson(double * x, double * p) {
}
void testDiscDistr() {
cout << "\nTest Discrete distributions\n\n";
TH1D * h1 =
new TH1D(
"h1PS",
"Unuran Poisson prob",20,0,20);
TH1D * h2 =
new TH1D(
"h2PS",
"Poisson dist from TRandom",20,0,20);
double mu = 5;
TF1 * f =
new TF1(
"fps",poisson,1,0,1);
bool ret = unr.
Init(dist2,
"dari");
if (!ret) return;
int n = NGEN;
for (
int i = 0; i <
n; ++i) {
}
cout <<
"Time using Unuran method " << unr.
MethodName() <<
"\t=\t\t " << w.
CpuTime() << endl;
for (
int i = 0; i <
n; ++i) {
h2->
Fill( gRandom->Poisson(mu) );
}
cout <<
"Time using TRandom::Poisson " <<
"\t=\t\t " << w.
CpuTime() << endl;
std::cout << " chi2 test of UNURAN vs TRandom generated histograms: " << std::endl;
}
void testEmpDistr() {
cout << "\nTest Empirical distributions using smoothing\n\n";
const int Ndata = 1000;
double x[Ndata];
for (int i = 0; i < Ndata; ++i) {
if (i < 0.5*Ndata )
x[i] = gRandom->Gaus(-1.,1.);
else
x[i] = gRandom->Gaus(1.,3.);
}
TH1D * h0 =
new TH1D(
"h0Ref",
"Starting data",100,-10,10);
TH1D * h1 =
new TH1D(
"h1Unr",
"Unuran unbin Generated data",100,-10,10);
TH1D * h1b =
new TH1D(
"h1bUnr",
"Unuran bin Generated data",100,-10,10);
TH1D * h2 =
new TH1D(
"h2GR",
"Data from TH1::GetRandom",100,-10,10);
int n = NGEN;
for (
int i = 0; i <
n; ++i) {
}
cout <<
"Time using Unuran unbin " << unr.
MethodName() <<
"\t=\t\t " << w.
CpuTime() << endl;
if (!unr.
Init(binDist))
return;
for (
int i = 0; i <
n; ++i) {
}
cout <<
"Time using Unuran bin " << unr.
MethodName() <<
"\t=\t\t " << w.
CpuTime() << endl;
for (
int i = 0; i <
n; ++i) {
}
cout <<
"Time using TH1::GetRandom " <<
"\t=\t\t " << w.
CpuTime() << endl;
}
void unuranDemo() {
c1 =
new TCanvas(
"c1_unuranMulti",
"Multidimensional distribution",10,10,1000,1000);
testStringAPI();
testDistr1D();
testDistrMultiDim();
testDiscDistr();
testEmpDistr();
}