Logo ROOT  
Reference Guide
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
22RooMCIntegrator implements 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 "RooFit.h"
29#include "Riostream.h"
30
31#include "TMath.h"
32#include "TClass.h"
33#include "RooMCIntegrator.h"
34#include "RooArgSet.h"
35#include "RooNumber.h"
36#include "RooAbsArg.h"
37#include "RooNumIntFactory.h"
38#include "RooRealVar.h"
39#include "RooCategory.h"
40#include "RooMsgService.h"
41
42#include <math.h>
43
44
45
46using namespace std;
47
49;
50
51// Register this class with RooNumIntFactory
52
53
54////////////////////////////////////////////////////////////////////////////////
55/// This function registers class RooMCIntegrator, its configuration options
56/// and its capabilities with RooNumIntFactory
57
59{
60 RooCategory samplingMode("samplingMode","Sampling Mode") ;
61 samplingMode.defineType("Importance",RooMCIntegrator::Importance) ;
62 samplingMode.defineType("ImportanceOnly",RooMCIntegrator::ImportanceOnly) ;
63 samplingMode.defineType("Stratified",RooMCIntegrator::Stratified) ;
65
66 RooCategory genType("genType","Generator Type") ;
67 genType.defineType("QuasiRandom",RooMCIntegrator::QuasiRandom) ;
68 genType.defineType("PseudoRandom",RooMCIntegrator::PseudoRandom) ;
70
71 RooCategory verbose("verbose","Verbose flag") ;
72 verbose.defineType("true",1) ;
73 verbose.defineType("false",0) ;
74 verbose.setIndex(0) ;
75
76 RooRealVar alpha("alpha","Grid structure constant",1.5) ;
77 RooRealVar nRefineIter("nRefineIter","Number of refining iterations",5) ;
78 RooRealVar nRefinePerDim("nRefinePerDim","Number of refining samples (per dimension)",1000) ;
79 RooRealVar nIntPerDim("nIntPerDim","Number of integration samples (per dimension)",5000) ;
80
81 // Create prototype integrator
83
84 // Register prototype and default config with factory
85 fact.storeProtoIntegrator(proto,RooArgSet(samplingMode,genType,verbose,alpha,nRefineIter,nRefinePerDim,nIntPerDim)) ;
86
87 // Make this method the default for all N>2-dim integrals
89}
90
91
92
93////////////////////////////////////////////////////////////////////////////////
94/// Default constructor
95///
96/// coverity[UNINIT_CTOR]
97
99{
100}
101
102
103
104////////////////////////////////////////////////////////////////////////////////
105/// Construct an integrator over 'function' with given sampling mode
106/// and generator type. The sampling mode can be 'Importance'
107/// (default), 'ImportanceOnly' and 'Stratified'. The generator type
108/// can be 'QuasiRandom' (default) and 'PseudoRandom'. Consult the original
109/// VEGAS documentation on details of the mode and type parameters.
110
112 GeneratorType genType, Bool_t verbose) :
113 RooAbsIntegrator(function), _grid(function), _verbose(verbose),
114 _alpha(1.5), _mode(mode), _genType(genType),
115 _nRefineIter(5),_nRefinePerDim(1000),_nIntegratePerDim(5000)
116{
117 // coverity[UNINIT_CTOR]
118 if(!(_valid= _grid.isValid())) return;
119 if(_verbose) _grid.Print();
120}
121
122
123
124////////////////////////////////////////////////////////////////////////////////
125/// Construct an integrator over 'function' where the configuration details
126/// are taken from 'config'
127
130{
131 const RooArgSet& configSet = config.getConfigSection(IsA()->GetName()) ;
132 _verbose = (Bool_t) configSet.getCatIndex("verbose",0) ;
133 _alpha = configSet.getRealValue("alpha",1.5) ;
134 _mode = (SamplingMode) configSet.getCatIndex("samplingMode",Importance) ;
135 _genType = (GeneratorType) configSet.getCatIndex("genType",QuasiRandom) ;
136 _nRefineIter = (Int_t) configSet.getRealValue("nRefineIter",5) ;
137 _nRefinePerDim = (Int_t) configSet.getRealValue("nRefinePerDim",1000) ;
138 _nIntegratePerDim = (Int_t) configSet.getRealValue("nIntPerDim",5000) ;
139
140 // check that our grid initialized without errors
141 if(!(_valid= _grid.isValid())) return;
142 if(_verbose) _grid.Print();
143}
144
145
146
147////////////////////////////////////////////////////////////////////////////////
148/// Return clone of this generator operating on given function with given configuration
149/// Needed to support RooNumIntFactory
150
152{
153 return new RooMCIntegrator(function,config) ;
154}
155
156
157
158////////////////////////////////////////////////////////////////////////////////
159/// Destructor
160
162{
163}
164
165
166
167////////////////////////////////////////////////////////////////////////////////
168/// Check if we can integrate over the current domain. If return value
169/// is kTRUE we cannot handle the current limits (e.g. where the domain
170/// of one or more observables is open ended.
171
173{
174 return _grid.initialize(*integrand());
175}
176
177
178
179////////////////////////////////////////////////////////////////////////////////
180/// Evaluate the integral using a fixed number of calls to evaluate the integrand
181/// equal to about 10k per dimension. Use the first 5k calls to refine the grid
182/// over 5 iterations of 1k calls each, and the remaining 5k calls for a single
183/// high statistics integration.
184
186{
190 return ret ;
191}
192
193
194
195////////////////////////////////////////////////////////////////////////////////
196/// Perform one step of Monte Carlo integration using the specified number of iterations
197/// with (approximately) the specified number of integrand evaluation calls per iteration.
198/// Use the VEGAS algorithm, starting from the specified stage. Returns the best estimate
199/// of the integral. Also sets *absError to the estimated absolute error of the integral
200/// estimate if absError is non-zero.
201
202Double_t RooMCIntegrator::vegas(Stage stage, UInt_t calls, UInt_t iterations, Double_t *absError)
203{
204 //cout << "VEGAS stage = " << stage << " calls = " << calls << " iterations = " << iterations << endl ;
205
206 // reset the grid to its initial state if we are starting from scratch
207 if(stage == AllStages) _grid.initialize(*_function);
208
209 // reset the results of previous calculations on this grid, but reuse the grid itself.
210 if(stage <= ReuseGrid) {
211 _wtd_int_sum = 0;
212 _sum_wgts = 0;
213 _chi_sum = 0;
214 _it_num = 1;
215 _samples = 0;
216 }
217
218 // refine the results of previous calculations on the current grid.
219 if(stage <= RefineGrid) {
221 UInt_t boxes = 1;
223
224 // select the sampling mode for the next step
225 if(_mode != ImportanceOnly) {
226 // calculate the largest number of equal subdivisions ("boxes") along each
227 // axis that results in an average of no more than 2 integrand calls per cell
228 boxes = (UInt_t)floor(TMath::Power(calls/2.0,1.0/dim));
229 // use stratified sampling if we are allowed enough calls (or equivalently,
230 // if the dimension is low enough)
232 if (2*boxes >= RooGrid::maxBins) {
234 // adjust the number of bins and boxes to give an integral number >= 1 of boxes per bin
235 Int_t box_per_bin= (boxes > RooGrid::maxBins) ? boxes/RooGrid::maxBins : 1;
236 bins= boxes/box_per_bin;
237 if(bins > RooGrid::maxBins) bins= RooGrid::maxBins;
238 boxes = box_per_bin * bins;
239 oocxcoutD((TObject*)0,Integration) << "RooMCIntegrator: using stratified sampling with " << bins << " bins and "
240 << box_per_bin << " boxes/bin" << endl;
241 }
242 else {
243 oocxcoutD((TObject*)0,Integration) << "RooMCIntegrator: using importance sampling with " << bins << " bins and "
244 << boxes << " boxes" << endl;
245 }
246 }
247
248 // calculate the total number of n-dim boxes for this step
249 Double_t tot_boxes = TMath::Power((Double_t)boxes,(Double_t)dim);
250
251 // increase the total number of calls to get at least 2 calls per box, if necessary
252 _calls_per_box = (UInt_t)(calls/tot_boxes);
254 calls= (UInt_t)(_calls_per_box*tot_boxes);
255
256 // calculate the Jacobean factor: volume/(avg # of calls/bin)
257 _jac = _grid.getVolume()*TMath::Power((Double_t)bins,(Double_t)dim)/calls;
258
259 // setup our grid to use the calculated number of boxes and bins
260 _grid.setNBoxes(boxes);
261 if(bins != _grid.getNBins()) _grid.resize(bins);
262 }
263
264 // allocate memory for some book-keeping arrays
268
269 // loop over iterations for this step
270 Double_t cum_int(0),cum_sig(0);
272 _chisq = 0.0;
273 for (UInt_t it = 0; it < iterations; it++) {
274 Double_t intgrl(0),intgrl_sq(0),sig(0),jacbin(_jac);
275
276 _it_num = _it_start + it;
277
278 // reset the values associated with each grid cell
280
281 // loop over grid boxes
283 do {
284 Double_t m(0),q(0);
285 // loop over integrand evaluations within this grid box
286 for(UInt_t k = 0; k < _calls_per_box; k++) {
287 // generate a random point in this box
288 Double_t bin_vol(0);
289 _grid.generatePoint(box, x, bin, bin_vol, _genType == QuasiRandom ? kTRUE : kFALSE);
290 // evaluate the integrand at the generated point
291 Double_t fval= jacbin*bin_vol*integrand(x);
292 // update mean and variance calculations
293 Double_t d = fval - m;
294 m+= d / (k + 1.0);
295 q+= d * d * (k / (k + 1.0));
296 // accumulate the results of this evaluation (importance sampling only)
297 if (_mode != Stratified) _grid.accumulate(bin, fval*fval);
298 }
299 intgrl += m * _calls_per_box;
300 Double_t f_sq_sum = q * _calls_per_box ;
301 sig += f_sq_sum ;
302
303 // accumulate the results for this grid box (stratified sampling only)
304 if (_mode == Stratified) _grid.accumulate(bin, f_sq_sum);
305
306 // print occasional progress messages
307 if(_timer.RealTime() > 30) {
308 std::size_t index = 0;
309 std::size_t sizeOfDim = 1;
310
311 for (unsigned int i=0; i < _grid.getDimension(); ++i) {
312 index += box[i] * sizeOfDim;
313 sizeOfDim *= _grid.getNBoxes();
314 }
315 oocoutP(this,Integration) << "RooMCIntegrator: still working ... iteration "
316 << it << '/' << iterations << " box " << index << "/"<< std::pow(_grid.getNBoxes(), _grid.getDimension()) << endl;
318 }
319 else {
321 }
322
323 } while(_grid.nextBox(box));
324
325 // compute final results for this iteration
326 Double_t wgt;
327 sig = sig / (_calls_per_box - 1.0) ;
328 if (sig > 0) {
329 wgt = 1.0 / sig;
330 }
331 else if (_sum_wgts > 0) {
332 wgt = _sum_wgts / _samples;
333 }
334 else {
335 wgt = 0.0;
336 }
337 intgrl_sq = intgrl * intgrl;
338 _result = intgrl;
339 _sigma = sqrt(sig);
340
341 if (wgt > 0.0) {
342 _samples++ ;
343 _sum_wgts += wgt;
344 _wtd_int_sum += intgrl * wgt;
345 _chi_sum += intgrl_sq * wgt;
346
347 cum_int = _wtd_int_sum / _sum_wgts;
348 cum_sig = sqrt (1 / _sum_wgts);
349
350 if (_samples > 1) {
351 _chisq = (_chi_sum - _wtd_int_sum * cum_int)/(_samples - 1.0);
352 }
353 }
354 else {
355 cum_int += (intgrl - cum_int) / (it + 1.0);
356 cum_sig = 0.0;
357 }
358 oocxcoutD((TObject*)0,Integration) << "=== Iteration " << _it_num << " : I = " << intgrl << " +/- " << sqrt(sig) << endl
359 << " Cumulative : I = " << cum_int << " +/- " << cum_sig << "( chi2 = " << _chisq
360 << ")" << endl;
361 // print the grid after the final iteration
362 if (oodologD((TObject*)0,Integration)) {
363 if(it + 1 == iterations) _grid.Print("V");
364 }
366 }
367
368 // cleanup
369 delete[] bin;
370 delete[] box;
371 delete[] x;
372
373 if(absError) *absError = cum_sig;
374 return cum_int;
375}
#define d(i)
Definition: RSha256.hxx:102
#define oocxcoutD(o, a)
Definition: RooMsgService.h:83
#define oodologD(o, a)
Definition: RooMsgService.h:72
#define oocoutP(o, a)
Definition: RooMsgService.h:46
int Int_t
Definition: RtypesCore.h:43
unsigned int UInt_t
Definition: RtypesCore.h:44
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
float * q
Definition: THbookFile.cxx:87
double pow(double, double)
double floor(double)
double sqrt(double)
const char * proto
Definition: civetweb.c:16604
Int_t getCatIndex(const char *name, Int_t defVal=0, Bool_t verbose=kFALSE) const
Get index value of a RooAbsCategory stored in set with given name.
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
const RooAbsFunc * _function
const RooAbsFunc * integrand() const
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:23
bool defineType(const std::string &label)
Define a state with given name.
virtual Bool_t setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
virtual Bool_t setIndex(Int_t index, bool printError=true) override
Set value by specifying the index code of the desired state.
UInt_t getDimension() const
Definition: RooGrid.h:41
Bool_t initialize(const RooAbsFunc &function)
Calculate and store the grid dimensions and volume using the specified function, and initialize the g...
Definition: RooGrid.cxx:108
@ maxBins
Definition: RooGrid.h:61
UInt_t * createIndexVector() const
Definition: RooGrid.h:48
Bool_t isValid() const
Definition: RooGrid.h:40
UInt_t getNBins() const
Definition: RooGrid.h:43
void refine(Double_t alpha=1.5)
Refine the grid using the values that have been accumulated so far.
Definition: RooGrid.cxx:338
UInt_t getNBoxes() const
Definition: RooGrid.h:44
void resetValues()
Reset the values associated with each grid cell.
Definition: RooGrid.cxx:181
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
Definition: RooGrid.h:36
Double_t * createPoint() const
Definition: RooGrid.h:47
Double_t getVolume() const
Definition: RooGrid.h:42
void resize(UInt_t bins)
Adjust the subdivision of each axis to give the specified number of bins, using an algorithm that pre...
Definition: RooGrid.cxx:144
void generatePoint(const UInt_t box[], Double_t x[], UInt_t bin[], Double_t &vol, Bool_t useQuasiRandom=kTRUE) const
Generate a random vector in the specified box and and store its coordinates in the x[] array provided...
Definition: RooGrid.cxx:198
void setNBoxes(UInt_t boxes)
Definition: RooGrid.h:45
void accumulate(const UInt_t bin[], Double_t amount)
Add the specified amount to bin[j] of the 1D histograms associated with each axis j.
Definition: RooGrid.cxx:326
void firstBox(UInt_t box[]) const
Reset the specified array of box indices to refer to the first box in the standard traversal order.
Definition: RooGrid.cxx:246
Bool_t nextBox(UInt_t box[]) const
Update the specified array of box indices to refer to the next box in the standard traversal order an...
Definition: RooGrid.cxx:258
RooMCIntegrator implements an adaptive multi-dimensional Monte Carlo numerical integration,...
GeneratorType _genType
RooMCIntegrator()
Default constructor.
static void registerIntegrator(RooNumIntFactory &fact)
This function registers class RooMCIntegrator, its configuration options and its capabilities with Ro...
TStopwatch _timer
virtual Double_t integral(const Double_t *yvec=0)
Evaluate the integral using a fixed number of calls to evaluate the integrand equal to about 10k per ...
virtual ~RooMCIntegrator()
Destructor.
Double_t vegas(Stage stage, UInt_t calls, UInt_t iterations, Double_t *absError=0)
Perform one step of Monte Carlo integration using the specified number of iterations with (approximat...
virtual RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const
Return clone of this generator operating on given function with given configuration Needed to support...
virtual Bool_t checkLimits() const
Check if we can integrate over the current domain.
Double_t _wtd_int_sum
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
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.
RooNumIntFactory is a factory to instantiate numeric integrators from a given function binding and a ...
Bool_t storeProtoIntegrator(RooAbsIntegrator *proto, const RooArgSet &defConfig, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
Mother of all ROOT objects.
Definition: TObject.h:37
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:357
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
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
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:151
@ Integration
Definition: RooGlobalFunc.h:67
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
auto * m
Definition: textangle.C:8