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