Numeric algorithm tuning: caching of slow numeric integrals and parameterizations of slow numeric integrals
Author: Clemens Lange, Wouter Verkerke (C++ version)
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, April 17, 2024 at 11:20 AM.
import sys
import ROOT
def getWorkspace(mode):
# Create, save or load workspace with pdf
# -----------------------------------------------------------------------------------
#
# Mode = 0 : Create workspace for plain running (no integral caching)
# Mode = 1 : Generate workspace with precalculated integral and store it on file
# Mode = 2 : Load previously stored workspace from file
w = ROOT.RooWorkspace()
if mode != 2:
# Create empty workspace workspace
w = ROOT.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 precalculate normalization integral that integrate at least
# two dimensions numerically. In self 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)
model = w["model"]
model.setStringAttribute("CACHEPARMINT", "x:y:z")
# Evaluate pdf once to trigger filling of cache
normSet = {w["x"], w["y"], w["z"]}
model.getVal(normSet)
w.writeToFile("rf903_numintcache.root")
if mode == 2:
# Load preexisting workspace from file in mode==2
f = ROOT.TFile("rf903_numintcache.root")
w = f.Get("w")
# Return created or loaded workspace
return w
mode = 0
Mode = 0 : Run plain fit (slow) Mode = 1 : Generate workspace with precalculated integral and store it on file (prepare for accelerated running) Mode = 2 : Run fit from previously stored workspace including cached integrals (fast, run in mode=1 first)
Make/load workspace, here in mode 1
w = getWorkspace(mode)
if mode == 1:
# Show workspace that was created
w.Print()
# Show plot of cached integral values
hhcache = w.expensiveObjectCache().getObj(1)
if hhcache:
ROOT.TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600)
hhcache.createHistogram("a").Draw()
else:
ROOT.RooFit.Error("rf903_numintcache", "Cached histogram is not existing in workspace")
sys.exit()
ROOT.This is always slow (need to find maximum function value empirically in 3D space)
model = w["model"]
d = model.generate({w["x"], w["y"], w["z"]}, 1000)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(x,y,z)
ROOT.This is slow in mode 0, fast in mode 1
model.fitTo(d, Verbose=True, Timer=True, PrintLevel=-1)
<cppyy.gbl.RooFitResult object at 0x(nil)>
[#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.810 [#1] INFO:Minimization -- Session timer: Real time 0:00:03, CP time 3.810 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.050 [#1] INFO:Minimization -- Session timer: Real time 0:00:04, CP time 4.860, 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, at fitted value of a is not cached)
framex = w["x"].frame(Title="Projection of 3D model on X")
d.plotOn(framex)
model.plotOn(framex)
<cppyy.gbl.RooPlot object at 0x998cb40>
[#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
c = ROOT.TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600)
framex.Draw()
c.SaveAs("rf903_numintcache.png")
Info in <TCanvas::Print>: png file rf903_numintcache.png has been created
Make workspace available on command line after macro finishes
ROOT.gDirectory.Add(w)
Draw all canvases
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()