Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooMCIntegrator.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooMCIntegrator.cxx
19\class RooMCIntegrator
20\ingroup Roofitcore
21
22Implements an adaptive multi-dimensional Monte Carlo
23numerical integration, following the VEGAS algorithm originally described
24in G. P. Lepage, J. Comp. Phys. 27, 192(1978). This implementation is
25based on a C version from the 0.9 beta release of the GNU scientific library.
26**/
27
28#include "Riostream.h"
29
30#include "TMath.h"
31#include "RooMCIntegrator.h"
32#include "RooArgSet.h"
33#include "RooNumber.h"
34#include "RooNumIntFactory.h"
35#include "RooRealVar.h"
36#include "RooCategory.h"
37#include "RooMsgService.h"
38
39#include <cmath>
40
41
42using namespace std;
43
44
45// Register this class with RooNumIntFactory
46
47
48////////////////////////////////////////////////////////////////////////////////
49/// This function registers class RooMCIntegrator, its configuration options
50/// and its capabilities with RooNumIntFactory
51
53{
54 RooCategory samplingMode("samplingMode","Sampling Mode") ;
55 samplingMode.defineType("Importance",RooMCIntegrator::Importance) ;
56 samplingMode.defineType("ImportanceOnly",RooMCIntegrator::ImportanceOnly) ;
57 samplingMode.defineType("Stratified",RooMCIntegrator::Stratified) ;
59
60 RooCategory genType("genType","Generator Type") ;
61 genType.defineType("QuasiRandom",RooMCIntegrator::QuasiRandom) ;
62 genType.defineType("PseudoRandom",RooMCIntegrator::PseudoRandom) ;
64
65 RooCategory verbose("verbose","Verbose flag") ;
66 verbose.defineType("true",1) ;
67 verbose.defineType("false",0) ;
68 verbose.setIndex(0) ;
69
70 RooRealVar alpha("alpha","Grid structure constant",1.5) ;
71 RooRealVar nRefineIter("nRefineIter","Number of refining iterations",5) ;
72 RooRealVar nRefinePerDim("nRefinePerDim","Number of refining samples (per dimension)",1000) ;
73 RooRealVar nIntPerDim("nIntPerDim","Number of integration samples (per dimension)",5000) ;
74
75 // Create prototype integrator
76 auto creator = [](const RooAbsFunc& function, const RooNumIntConfig& config) {
77 return std::make_unique<RooMCIntegrator>(function,config);
78 };
79
80 // Register prototype and default config with factory
81 std::string name = "RooMCIntegrator";
82 fact.registerPlugin(name, creator, {samplingMode,genType,verbose,alpha,nRefineIter,nRefinePerDim,nIntPerDim},
83 /*canIntegrate1D=*/true,
84 /*canIntegrate2D=*/true,
85 /*canIntegrateND=*/true,
86 /*canIntegrateOpenEnded=*/false);
87
88 // Make this method the default for all N>2-dim integrals
90}
91
92
93////////////////////////////////////////////////////////////////////////////////
94/// Construct an integrator over 'function' with given sampling mode
95/// and generator type. The sampling mode can be 'Importance'
96/// (default), 'ImportanceOnly' and 'Stratified'. The generator type
97/// can be 'QuasiRandom' (default) and 'PseudoRandom'. Consult the original
98/// VEGAS documentation on details of the mode and type parameters.
99
101 GeneratorType genType, bool verbose) :
102 RooAbsIntegrator(function), _grid(function), _verbose(verbose),
103 _alpha(1.5), _mode(mode), _genType(genType),
104 _nRefineIter(5),_nRefinePerDim(1000),_nIntegratePerDim(5000)
105{
106 // coverity[UNINIT_CTOR]
107 if(!(_valid= _grid.isValid())) return;
108 if(_verbose) _grid.print(std::cout);
109}
110
111
112
113////////////////////////////////////////////////////////////////////////////////
114/// Construct an integrator over 'function' where the configuration details
115/// are taken from 'config'
116
119{
120 const RooArgSet& configSet = config.getConfigSection("RooMCIntegrator") ;
121 _verbose = (bool) configSet.getCatIndex("verbose",0) ;
122 _alpha = configSet.getRealValue("alpha",1.5) ;
123 _mode = (SamplingMode) configSet.getCatIndex("samplingMode",Importance) ;
124 _genType = (GeneratorType) configSet.getCatIndex("genType",QuasiRandom) ;
125 _nRefineIter = (Int_t) configSet.getRealValue("nRefineIter",5) ;
126 _nRefinePerDim = (Int_t) configSet.getRealValue("nRefinePerDim",1000) ;
127 _nIntegratePerDim = (Int_t) configSet.getRealValue("nIntPerDim",5000) ;
128
129 // check that our grid initialized without errors
130 if(!(_valid= _grid.isValid())) return;
131 if(_verbose) _grid.print(std::cout);
132}
133
134////////////////////////////////////////////////////////////////////////////////
135/// Check if we can integrate over the current domain. If return value
136/// is true we cannot handle the current limits (e.g. where the domain
137/// of one or more observables is open ended.
138
140{
141 return _grid.initialize(*integrand());
142}
143
144
145
146////////////////////////////////////////////////////////////////////////////////
147/// Evaluate the integral using a fixed number of calls to evaluate the integrand
148/// equal to about 10k per dimension. Use the first 5k calls to refine the grid
149/// over 5 iterations of 1k calls each, and the remaining 5k calls for a single
150/// high statistics integration.
151
152double RooMCIntegrator::integral(const double* /*yvec*/)
153{
154 _timer.Start(true);
156 double ret = vegas(ReuseGrid,_nIntegratePerDim*_grid.getDimension(),1);
157 return ret ;
158}
159
160
161
162////////////////////////////////////////////////////////////////////////////////
163/// Perform one step of Monte Carlo integration using the specified number of iterations
164/// with (approximately) the specified number of integrand evaluation calls per iteration.
165/// Use the VEGAS algorithm, starting from the specified stage. Returns the best estimate
166/// of the integral. Also sets *absError to the estimated absolute error of the integral
167/// estimate if absError is non-zero.
168
169double RooMCIntegrator::vegas(Stage stage, UInt_t calls, UInt_t iterations, double *absError)
170{
171 //cout << "VEGAS stage = " << stage << " calls = " << calls << " iterations = " << iterations << endl ;
172
173 // reset the grid to its initial state if we are starting from scratch
174 if(stage == AllStages) _grid.initialize(*_function);
175
176 // reset the results of previous calculations on this grid, but reuse the grid itself.
177 if(stage <= ReuseGrid) {
178 _wtd_int_sum = 0;
179 _sum_wgts = 0;
180 _chi_sum = 0;
181 _it_num = 1;
182 _samples = 0;
183 }
184
185 // refine the results of previous calculations on the current grid.
186 if(stage <= RefineGrid) {
187 UInt_t bins = RooGrid::maxBins;
188 UInt_t boxes = 1;
189 UInt_t dim(_grid.getDimension());
190
191 // select the sampling mode for the next step
192 if(_mode != ImportanceOnly) {
193 // calculate the largest number of equal subdivisions ("boxes") along each
194 // axis that results in an average of no more than 2 integrand calls per cell
195 boxes = (UInt_t)floor(TMath::Power(calls/2.0,1.0/dim));
196 // use stratified sampling if we are allowed enough calls (or equivalently,
197 // if the dimension is low enough)
199 if (2*boxes >= RooGrid::maxBins) {
201 // adjust the number of bins and boxes to give an integral number >= 1 of boxes per bin
202 Int_t box_per_bin= (boxes > RooGrid::maxBins) ? boxes/RooGrid::maxBins : 1;
203 bins= boxes/box_per_bin;
204 if(bins > RooGrid::maxBins) bins= RooGrid::maxBins;
205 boxes = box_per_bin * bins;
206 oocxcoutD((TObject*)nullptr,Integration) << "RooMCIntegrator: using stratified sampling with " << bins << " bins and "
207 << box_per_bin << " boxes/bin" << endl;
208 }
209 else {
210 oocxcoutD((TObject*)nullptr,Integration) << "RooMCIntegrator: using importance sampling with " << bins << " bins and "
211 << boxes << " boxes" << endl;
212 }
213 }
214
215 // calculate the total number of n-dim boxes for this step
216 double tot_boxes = TMath::Power((double)boxes,(double)dim);
217
218 // increase the total number of calls to get at least 2 calls per box, if necessary
219 _calls_per_box = (UInt_t)(calls/tot_boxes);
221 calls= (UInt_t)(_calls_per_box*tot_boxes);
222
223 // calculate the Jacobean factor: volume/(avg # of calls/bin)
224 _jac = _grid.getVolume()*TMath::Power((double)bins,(double)dim)/calls;
225
226 // setup our grid to use the calculated number of boxes and bins
227 _grid.setNBoxes(boxes);
228 if(bins != _grid.getNBins()) _grid.resize(bins);
229 }
230
231 // allocate memory for some book-keeping arrays
232 std::vector<UInt_t> box(_grid.getDimension());
233 std::vector<UInt_t> bin(_grid.getDimension());
234 std::vector<double> x(_grid.getDimension());
235
236
237 // loop over iterations for this step
238 double cum_int(0);
239 double cum_sig(0);
241 _chisq = 0.0;
242 for (UInt_t it = 0; it < iterations; it++) {
243 double intgrl(0);
244 double intgrl_sq(0);
245 double sig(0);
246 double jacbin(_jac);
247
248 _it_num = _it_start + it;
249
250 // reset the values associated with each grid cell
251 _grid.resetValues();
252
253 // loop over grid boxes
254 _grid.firstBox(box.data());
255 do {
256 double m(0);
257 double q(0);
258 // loop over integrand evaluations within this grid box
259 for(UInt_t k = 0; k < _calls_per_box; k++) {
260 // generate a random point in this box
261 double bin_vol(0);
262 _grid.generatePoint(box.data(), x.data(), bin.data(), bin_vol, _genType == QuasiRandom ? true : false);
263 // evaluate the integrand at the generated point
264 double fval= jacbin*bin_vol*integrand(x.data());
265 // update mean and variance calculations
266 double d = fval - m;
267 m+= d / (k + 1.0);
268 q+= d * d * (k / (k + 1.0));
269 // accumulate the results of this evaluation (importance sampling only)
270 if (_mode != Stratified) _grid.accumulate(bin.data(), fval*fval);
271 }
272 intgrl += m * _calls_per_box;
273 double f_sq_sum = q * _calls_per_box ;
274 sig += f_sq_sum ;
275
276 // accumulate the results for this grid box (stratified sampling only)
277 if (_mode == Stratified) _grid.accumulate(bin.data(), f_sq_sum);
278
279 // print occasional progress messages
280 if(_timer.RealTime() > 30) {
281 std::size_t index = 0;
282 std::size_t sizeOfDim = 1;
283
284 for (unsigned int i=0; i < _grid.getDimension(); ++i) {
285 index += box[i] * sizeOfDim;
286 sizeOfDim *= _grid.getNBoxes();
287 }
288 oocoutP(nullptr, Integration) << "RooMCIntegrator: still working ... iteration "
289 << it << '/' << iterations << " box " << index << "/"<< std::pow(_grid.getNBoxes(), _grid.getDimension()) << endl;
290 _timer.Start(true);
291 }
292 else {
293 _timer.Start(false);
294 }
295
296 } while(_grid.nextBox(box.data()));
297
298 // compute final results for this iteration
299 double wgt;
300 sig = sig / (_calls_per_box - 1.0) ;
301 if (sig > 0) {
302 wgt = 1.0 / sig;
303 }
304 else if (_sum_wgts > 0) {
305 wgt = _sum_wgts / _samples;
306 }
307 else {
308 wgt = 0.0;
309 }
310 intgrl_sq = intgrl * intgrl;
311 _result = intgrl;
312 _sigma = sqrt(sig);
313
314 if (wgt > 0.0) {
315 _samples++ ;
316 _sum_wgts += wgt;
317 _wtd_int_sum += intgrl * wgt;
318 _chi_sum += intgrl_sq * wgt;
319
320 cum_int = _wtd_int_sum / _sum_wgts;
321 cum_sig = sqrt (1 / _sum_wgts);
322
323 if (_samples > 1) {
324 _chisq = (_chi_sum - _wtd_int_sum * cum_int)/(_samples - 1.0);
325 }
326 }
327 else {
328 cum_int += (intgrl - cum_int) / (it + 1.0);
329 cum_sig = 0.0;
330 }
331 oocxcoutD((TObject*)nullptr,Integration) << "=== Iteration " << _it_num << " : I = " << intgrl << " +/- " << sqrt(sig) << endl
332 << " Cumulative : I = " << cum_int << " +/- " << cum_sig << "( chi2 = " << _chisq
333 << ")" << endl;
334 // print the grid after the final iteration
335 if (oodologD((TObject*)nullptr,Integration)) {
336 if(it + 1 == iterations) _grid.print(std::cout, true);
337 }
338 _grid.refine(_alpha);
339 }
340
341 if(absError) *absError = cum_sig;
342 return cum_int;
343}
#define d(i)
Definition RSha256.hxx:102
RooAbsReal & function()
bool _verbose
Verbose messaging if true.
#define oocxcoutD(o, a)
#define oodologD(o, a)
#define oocoutP(o, a)
int Int_t
Definition RtypesCore.h:45
unsigned int UInt_t
Definition RtypesCore.h:46
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char mode
char name[80]
Definition TGX11.cxx:110
float * q
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
Int_t getCatIndex(const char *name, Int_t defVal=0, bool verbose=false) const
Get index value of a RooAbsCategory stored in set with given name.
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition RooAbsFunc.h:27
Abstract interface for integrators of real-valued functions that implement the RooAbsFunc interface.
const RooAbsFunc * _function
Pointer to function binding of integrand.
const RooAbsFunc * integrand() const
Return integrand function binding.
bool _valid
Is integrator in valid state?
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Object to represent discrete states.
Definition RooCategory.h:28
bool setIndex(Int_t index, bool printError=true) override
Set value by specifying the index code of the desired state.
bool defineType(const std::string &label)
Define a state with given name.
bool setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
GeneratorType _genType
Generator type.
static void registerIntegrator(RooNumIntFactory &fact)
This function registers class RooMCIntegrator, its configuration options and its capabilities with Ro...
TStopwatch _timer
Timer.
Int_t _nIntegratePerDim
Number of integration samplings (per dim)
RooMCIntegrator(const RooAbsFunc &function, SamplingMode mode=Importance, GeneratorType genType=QuasiRandom, bool verbose=false)
Construct an integrator over 'function' with given sampling mode and generator type.
bool checkLimits() const override
Check if we can integrate over the current domain.
double _sigma
Scratch variables preserved between calls to vegas1/2/2.
double vegas(Stage stage, UInt_t calls, UInt_t iterations, double *absError=nullptr)
Perform one step of Monte Carlo integration using the specified number of iterations with (approximat...
bool _verbose
Verbosity control.
Int_t _mode
Sampling mode.
Int_t _nRefinePerDim
Number of refinement samplings (per dim)
double _alpha
Grid stiffness parameter.
UInt_t _calls_per_box
Scratch variables preserved between calls to vegas1/2/2.
Int_t _nRefineIter
Number of refinement iterations.
double integral(const double *yvec=nullptr) override
Evaluate the integral using a fixed number of calls to evaluate the integrand equal to about 10k per ...
Holds the configuration parameters of the various numeric integrators used by RooRealIntegral.
RooCategory & methodND()
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
static RooNumIntConfig & defaultConfig()
Return reference to instance of default numeric integrator configuration object.
Factory to instantiate numeric integrators from a given function binding and a given configuration.
bool registerPlugin(std::string const &name, Creator const &creator, const RooArgSet &defConfig, bool canIntegrate1D, bool canIntegrate2D, bool canIntegrateND, bool canIntegrateOpenEnded, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
Variable that can be changed from the outside.
Definition RooRealVar.h:37
Mother of all ROOT objects.
Definition TObject.h:41
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition fillpatterns.C:1
Double_t x[n]
Definition legend1.C:17
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721
TMarker m
Definition textangle.C:8