rf903_numintcache

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, November 30, 2022 at 11:24 AM.

In [1]:
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
Welcome to JupyROOT 6.27/01

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

In [2]:
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)

In [3]:
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

In [4]:
model.fitTo(d, Verbose=True, Timer=True)
Out[4]:
<cppyy.gbl.RooFitResult object at 0x(nil)>
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#0] WARNING:Minimization -- RooAbsMinimizerFcn::synchronize: WARNING: no initial error estimate available for a: using 1
 **********
 **    1 **SET PRINT           1
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 a            0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00
 **********
 **    3 **SET ERR         0.5
 **********
 **********
 **    4 **SET PRINT           1
 **********
 **********
 **    5 **SET STR           1
 **********
 NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY
 **********
 **    6 **MIGRAD         500           1
 **********
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.

prevFCN = 1659.930708   START MIGRAD MINIMIZATION.  STRATEGY  1.  CONVERGENCE WHEN EDM .LT. 1.00e-03
a=0.02003, [#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(x,y,z)

prevFCN = 1668.090831  a=-0.02003, 
prevFCN = 1666.341378  a=0.002003, 
prevFCN = 1660.094465  a=-0.002003, 
prevFCN = 1659.913535   FCN=1659.93 FROM MIGRAD    STATUS=INITIATE        4 CALLS           5 TOTAL
                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX       
  EXT PARAMETER               CURRENT GUESS       STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a            0.00000e+00   1.00000e+00   2.01358e-01   2.25805e+02
                               ERR DEF= 0.5
a=-0.001236, 
prevFCN = 1659.902781  a=-0.001036, 
prevFCN = 1659.903514  a=-0.001437, 
prevFCN = 1659.903515  a=-0.001089, 
prevFCN = 1659.903177  a=-0.001383, 
prevFCN = 1659.903178   MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
a=-0.001236, 
prevFCN = 1659.902781  a=-0.001089, 
prevFCN = 1659.903177  a=-0.001383, 
prevFCN = 1659.903178  a=-0.001207, 
prevFCN = 1659.902797  a=-0.001266, 
prevFCN = 1659.902797   COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=1659.9 FROM MIGRAD    STATUS=CONVERGED      14 CALLS          15 TOTAL
                     EDM=2.26676e-10    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a           -1.23627e-03   5.23021e-03   2.94345e-05  -1.43931e-02
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  1    ERR DEF=0.5
  2.736e-05 
[#1] INFO:Minimization -- Command timer: Real time 0:00:02, CP time 2.750
[#1] INFO:Minimization -- Session timer: Real time 0:00:02, CP time 2.750
a=-0.001236,  **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE         500
 **********

prevFCN = 1659.902781  a=-0.001207, 
prevFCN = 1659.902797  a=-0.001266, 
prevFCN = 1659.902797  a=-0.00123, 
prevFCN = 1659.902782  a=-0.001242, 
prevFCN = 1659.902782   COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=1659.9 FROM HESSE     STATUS=OK              5 CALLS          20 TOTAL
                     EDM=2.26062e-10    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  a           -1.23627e-03   5.23023e-03   5.88690e-06  -2.47254e-04
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  1    ERR DEF=0.5
  2.736e-05 
[#1] INFO:Minimization -- Command timer: Real time 0:00:00, CP time 0.960
[#1] INFO:Minimization -- Session timer: Real time 0:00:03, CP time 3.710, 2 slices
a=-0.001236, [#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 20

Projection on x (always slow as 2D integral over Y, at fitted value of a is not cached)

In [5]:
framex = w["x"].frame(Title="Projection of 3D model on X")
d.plotOn(framex)
model.plotOn(framex)
Out[5]:
<cppyy.gbl.RooPlot object at 0xa0fb5a0>
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) plot on x integrates over variables (z,y)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[y,z]_Norm[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(y,z)

Draw x projection on canvas

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

In [7]:
ROOT.gDirectory.Add(w)

Draw all canvases

In [8]:
from ROOT import gROOT 
gROOT.GetListOfCanvases().Draw()