Quasirandom

Example of generating quasi-random numbers

Author: Lorenzo Moneta
This notebook tutorial was automatically generated with ROOTBOOK-izer (Beta) from the macro found in the ROOT repository on Thursday, January 19, 2017 at 04:34 PM.



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>



In [3]:

    
%%cpp -d
// This is a workaround to make sure the namespace is used inside functions
using namespace ROOT::Math;

Arguments are defined.



In [4]:

    
int n = 10000;
int skip = 0;



In [5]:

    
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 [6]:

    
QuasiRandomSobol r1(2);
QuasiRandomNiederreiter r2(2);

Generate n random points



In [7]:

    
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);

Check uniformity



In [8]:

    
h1->Draw("COLZ");
c1->cd(3);
h2->Draw("COLZ");
gPad->Update();

Test number of empty bins



In [9]:

    
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 [10]:

    
%jsroot on
gROOT->GetListOfCanvases()->Draw()