Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rf903_numintcache.py File Reference

Detailed Description

View in nbviewer Open in SWAN
Numeric algorithm tuning: caching of slow numeric integrals and parameterizations of slow numeric integrals

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)
# Create, save or load workspace with pdf
# -----------------------------------------------------------------------------------
# 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()
# Use pdf from workspace for generation and fitting
# -----------------------------------------------------------------------------------
# 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)
# ROOT.This is slow in mode 0, fast in mode 1
model.fitTo(d, Verbose=True, Timer=True, PrintLevel=-1)
# 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)
# Draw x projection on canvas
c = ROOT.TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600)
framex.Draw()
c.SaveAs("rf903_numintcache.png")
# Make workspace available on command line after macro finishes
ROOT.gDirectory.Add(w)
th1 Draw()
[#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(x,y,z)
[#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:02, CP time 2.810
[#1] INFO:Minimization -- Session timer: Real time 0:00:02, CP time 2.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:00, CP time 0.990
[#1] INFO:Minimization -- Session timer: Real time 0:00:03, CP time 3.800, 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
[#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
Date
February 2018
Authors
Clemens Lange, Wouter Verkerke (C++ version)

Definition in file rf903_numintcache.py.