quasirandom

Example of generating quasi-random numbers

Author: Lorenzo Moneta
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, August 17, 2022 at 09:34 AM.

In [1]:
gSystem->Load("libMathMore");
In [2]:
%%cpp -d
#include "Math/QuasiRandom.h"
#include "Math/Random.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TStopwatch.h"

#include <iostream>

using namespace ROOT::Math;

Arguments are defined.

In [3]:
int n = 10000;
int skip = 0;
In [4]:
TH2D * h0 = new TH2D("h0","Pseudo-random Sequence",200,0,1,200,0,1);
TH2D * h1 = new TH2D("h1","Sobol Sequence",200,0,1,200,0,1);
TH2D * h2 = new TH2D("h2","Niederrer Sequence",200,0,1,200,0,1);

RandomMT         r0;

quasi random numbers need to be created giving the dimension of the sequence in this case we generate 2-d sequence

In [5]:
QuasiRandomSobol r1(2);
QuasiRandomNiederreiter r2(2);

generate n random points

In [6]:
double x[2];
TStopwatch w; w.Start();
for (int i = 0; i < n; ++i)  {
   r0.RndmArray(2,x);
   h0->Fill(x[0],x[1]);
}
std::cout << "Time for gRandom ";
w.Print();

w.Start();
if( skip>0) r1.Skip(skip);
for (int i = 0; i < n; ++i)  {
   r1.Next(x);
   h1->Fill(x[0],x[1]);
}
std::cout << "Time for Sobol ";
w.Print();

w.Start();
if( skip>0) r2.Skip(skip);
for (int i = 0; i < n; ++i)  {
   r2.Next(x);
   h2->Fill(x[0],x[1]);
}
std::cout << "Time for Niederreiter ";
w.Print();

TCanvas * c1 = new TCanvas("c1","Random sequence",600,1200);
c1->Divide(1,3);
c1->cd(1);
h0->Draw("COLZ");
c1->cd(2);
Time for gRandom Real time 0:00:00, CP time 0.000
Time for Sobol Real time 0:00:00, CP time 0.000
Time for Niederreiter Real time 0:00:00, CP time 0.000

check uniformity

In [7]:
h1->Draw("COLZ");
c1->cd(3);
h2->Draw("COLZ");
gPad->Update();

test number of empty bins

In [8]:
int nzerobins0 = 0;
int nzerobins1 = 0;
int nzerobins2 = 0;
for (int i = 1; i <= h1->GetNbinsX(); ++i) {
   for (int j = 1; j <= h1->GetNbinsY(); ++j) {
      if (h0->GetBinContent(i,j) == 0 ) nzerobins0++;
      if (h1->GetBinContent(i,j) == 0 ) nzerobins1++;
      if (h2->GetBinContent(i,j) == 0 ) nzerobins2++;
   }
}

std::cout << "number of empty bins for pseudo-random = " << nzerobins0 << std::endl;
std::cout << "number of empty bins for " << r1.Name() << "\t= " << nzerobins1 << std::endl;
std::cout << "number of empty bins for " << r2.Name() << "\t= " << nzerobins2 << std::endl;

int iret = 0;
if (nzerobins1 >= nzerobins0 ) iret += 1;
if (nzerobins2 >= nzerobins0 ) iret += 2;
return iret;
number of empty bins for pseudo-random = 31139
number of empty bins for sobol	= 30512
number of empty bins for niederreiter-base-2	= 30512

Draw all canvases

In [9]:
%jsroot on
gROOT->GetListOfCanvases()->Draw()