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