Numeric algorithm tuning: caching of slow numeric integrals and parameterization of slow numeric integrals
Author: Wouter Verkerke
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:18 PM.
%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooWorkspace.h"
#include "RooExpensiveObjectCache.h"
#include "TFile.h"
#include "TH1.h"
using namespace RooFit;
RooWorkspace *getWorkspace(Int_t mode);
Definition of a helper function:
%%cpp -d
RooWorkspace *getWorkspace(Int_t mode)
{
// C r e a t e , s a v e o r l o a d w o r k s p a c e w i t h p . d . f .
// -----------------------------------------------------------------------------------
//
// Mode = 0 : Create workspace for plain running (no integral caching)
// Mode = 1 : Generate workspace with pre-calculated integral and store it on file
// Mode = 2 : Load previously stored workspace from file
RooWorkspace *w(0);
if (mode != 2) {
// Create empty workspace workspace
w = new RooWorkspace("w", 1);
// Make a difficult to normalize pdf in 3 dimensions that is integrated numerically.
w->factory("EXPR::model('1/((x-a)*(x-a)+0.01)+1/((y-a)*(y-a)+0.01)+1/"
"((z-a)*(z-a)+0.01)',x[-1,1],y[-1,1],z[-1,1],a[-5,5])");
}
if (mode == 1) {
// Instruct model to pre-calculate normalization integral that integrate at least
// two dimensions numerically. In this specific case the integral value for
// all values of parameter 'a' are stored in a histogram and available for use
// in subsequent fitting and plotting operations (interpolation is applied)
// w->pdf("model")->setNormValueCaching(3) ;
w->pdf("model")->setStringAttribute("CACHEPARMINT", "x:y:z");
// Evaluate pdf once to trigger filling of cache
RooArgSet normSet(*w->var("x"), *w->var("y"), *w->var("z"));
w->pdf("model")->getVal(&normSet);
w->writeToFile("rf903_numintcache.root");
}
if (mode == 2) {
// Load preexisting workspace from file in mode==2
TFile *f = new TFile("rf903_numintcache.root");
w = (RooWorkspace *)f->Get("w");
}
// Return created or loaded workspace
return w;
}
Arguments are defined.
Int_t mode = 0;
Mode = 0 : Run plain fit (slow) Mode = 1 : Generate workspace with pre-calculated integral and store it on file (prepare for accelerated running) Mode = 2 : Run fit from previously stored workspace including cached integrals (fast, requires run in mode=1 first)
Make/load workspace, exit here in mode 1
RooWorkspace *w1 = getWorkspace(mode);
if (mode == 1) {
// Show workspace that was created
w1->Print();
// Show plot of cached integral values
RooDataHist *hhcache = (RooDataHist *)w1->expensiveObjectCache().getObj(1);
if (hhcache) {
new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
hhcache->createHistogram("a")->Draw();
} else {
Error("rf903_numintcache", "Cached histogram is not existing in workspace");
}
return;
}
This is always slow (need to find maximum function value empirically in 3D space)
std::unique_ptr<RooDataSet> d{w1->pdf("model")->generate({*w1->var("x"), *w1->var("y"), *w1->var("z")}, 1000)};
[#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(x,y,z)
This is slow in mode 0, but fast in mode 1
w1->pdf("model")->fitTo(*d, Verbose(true), Timer(true), PrintLevel(-1));
[#1] INFO:Fitting -- RooAbsPdf::fitTo(model_over_model_Int[x,y,z]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- using CPU computation library compiled with -mavx2 [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_model_over_model_Int[x,y,z]_modelData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#0] WARNING:Minimization -- RooAbsMinimizerFcn::synchronize: WARNING: no initial error estimate available for a: using 1 prevFCN = 1659.930708 a=0.02833, [#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(x,y,z) prevFCN = 1675.611563 a=-0.02833, prevFCN = 1673.217894 a=0.002833, prevFCN = 1660.205177 a=-0.002833, prevFCN = 1659.94939 a=0.0002833, prevFCN = 1659.944972 a=-0.0002833, prevFCN = 1659.919376 a=-0.001237, prevFCN = 1659.902781 a=-0.001089, prevFCN = 1659.903175 a=-0.001384, prevFCN = 1659.90318 a=-0.001237, prevFCN = 1659.902781 a=-0.001089, prevFCN = 1659.903175 a=-0.001384, prevFCN = 1659.90318 a=-0.001207, prevFCN = 1659.902797 a=-0.001266, prevFCN = 1659.902798 [#1] INFO:Minimization -- Command timer: Real time 0:00:03, CP time 3.300 [#1] INFO:Minimization -- Session timer: Real time 0:00:03, CP time 3.300 a=-0.001237, prevFCN = 1659.902781 a=-0.001207, prevFCN = 1659.902797 a=-0.001266, prevFCN = 1659.902798 a=-0.001231, prevFCN = 1659.902782 a=-0.001243, prevFCN = 1659.902782 [#1] INFO:Minimization -- Command timer: Real time 0:00:01, CP time 1.130 [#1] INFO:Minimization -- Session timer: Real time 0:00:04, CP time 4.430, 2 slices a=-0.001237, [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#0] WARNING:NumericIntegration -- RooAdaptiveIntegratorND::dtor(model) WARNING: Number of suppressed warningings about integral evaluations where target precision was not reached is 17
Projection on x (always slow as 2D integral over Y,Z at fitted value of a is not cached)
RooPlot *framex = w1->var("x")->frame(Title("Projection of 3D model on X"));
d->plotOn(framex);
w1->pdf("model")->plotOn(framex);
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) plot on x integrates over variables (y,z) [#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(x,y,z) [#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[y,z]_Norm[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(y,z) [#0] WARNING:NumericIntegration -- RooAdaptiveIntegratorND::dtor(model) WARNING: Number of suppressed warningings about integral evaluations where target precision was not reached is 1
Draw x projection on canvas
auto canv = new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
framex->Draw();
canv->Draw();
Make workspace available on command line after macro finishes
gDirectory->Add(w1);
Draw all canvases
%jsroot on
gROOT->GetListOfCanvases()->Draw()