Rf 9 0 3_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, January 19, 2022 at 10:39 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"
In [2]:
%%cpp -d
// This is a workaround to make sure the namespace is used inside functions
using namespace RooFit;
In [3]:
RooWorkspace *getWorkspace(Int_t mode);

A helper function is created:

In [4]:
%%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 [5]:
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 [6]:
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;
}
RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby 
                Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
                All rights reserved, please read http://roofit.sourceforge.net/license.txt

[#1] INFO:ObjectHandling -- RooWorkspace::exportToCint(w) INFO: references to all objects in this workspace will be created in CINT in 'namespace w'

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 [7]:
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 [8]:
w1->pdf("model")->fitTo(*d, Verbose(kTRUE), Timer(kTRUE));
[#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.740
[#1] INFO:Minimization -- Session timer: Real time 0:00:02, CP time 2.740
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.950
[#1] INFO:Minimization -- Session timer: Real time 0:00:03, CP time 3.690, 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 [9]:
RooPlot *framex = w1->var("x")->frame(Title("Projection of 3D model on X"));
d->plotOn(framex);
w1->pdf("model")->plotOn(framex);
[#0] WARNING:Plotting -- Cannot apply a bin width correction and use Poisson errors. Not correcting for bin width.
[#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 [10]:
auto canv = new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
framex->Draw();
canv->Draw();

Make workspace available on command line after macro finishes

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

Draw all canvases

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