# rf903_numintcache¶

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

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

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

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

## Create, save or load workspace with p.d.f.¶

Make/load workspace, exit here in mode 1

In [4]:
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;
}

input_line_52:13:16: warning: 'createHistogram' is deprecated: will be removed in ROOT v6.30: Use the overload of RooAbsData::createHistogram that takes RooFit command arguments. [-Wdeprecated-declarations]
hhcache->createHistogram("a")->Draw();
^
/home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/include/RooAbsData.h:229:7: note: 'createHistogram' has been explicitly marked deprecated here
R__DEPRECATED(6, 30, "Use the overload of RooAbsData::createHistogram that takes RooFit command arguments.");
^
/home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/include/ROOT/RConfig.hxx:517:3: note: expanded from macro 'R__DEPRECATED'
_R__JOIN3_(_R__DEPRECATED_,MAJOR,MINOR)("will be removed in ROOT v" #MAJOR "." #MINOR ": " REASON)
^
/home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/include/ROOT/RConfig.hxx:467:30: note: expanded from macro '_R__JOIN3_'
#   define _R__JOIN3_(F,X,Y) _NAME3_(F,X,Y)
^
/home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/include/ROOT/RConfig.hxx:449:39: note: expanded from macro '_NAME3_'
#   define _NAME3_(name1,name2,name3) name1##name2##name3
^
<scratch space>:9:1: note: expanded from here
_R__DEPRECATED_630
^
/home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/include/ROOT/RConfig.hxx:501:37: note: expanded from macro '_R__DEPRECATED_630'
# define _R__DEPRECATED_630(REASON) _R__DEPRECATED_LATER(REASON)
^
/home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/include/ROOT/RConfig.hxx:483:56: note: expanded from macro '_R__DEPRECATED_LATER'
#   define _R__DEPRECATED_LATER(REASON) __attribute__((deprecated(REASON)))
^


## Use p.d.f. from workspace for generation and fitting¶

This is always slow (need to find maximum function value empirically in 3D space)

In [5]:
RooDataSet *d = w1->pdf("model")->generate(RooArgSet(*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

In [6]:
w1->pdf("model")->fitTo(*d, Verbose(true), Timer(true));

[#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
**********
**********
**********
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
**********
**********
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.830
[#1] INFO:Minimization -- Session timer: Real time 0:00:02, CP time 2.830
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:01, CP time 1.180
[#1] INFO:Minimization -- Session timer: Real time 0:00:03, CP time 4.010, 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,Z at fitted value of a is not cached)

In [7]:
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[y,z]_Norm[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(y,z)


Draw x projection on canvas

In [8]:
auto canv = new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
framex->Draw();
canv->Draw();


Make workspace available on command line after macro finishes

In [9]:
gDirectory->Add(w1);


Draw all canvases

In [10]:
%jsroot on
gROOT->GetListOfCanvases()->Draw()