Logo ROOT  
Reference Guide
RooIntegrator1D.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 RooIntegrator1D.cxx
19\class RooIntegrator1D
20\ingroup Roofitcore
21
22RooIntegrator1D implements an adaptive one-dimensional
23numerical integration algorithm.
24**/
25
26
27#include "RooFit.h"
28#include "Riostream.h"
29
30#include "TClass.h"
31#include "RooIntegrator1D.h"
32#include "RooArgSet.h"
33#include "RooRealVar.h"
34#include "RooNumber.h"
36#include "RooNumIntConfig.h"
37#include "RooNumIntFactory.h"
38#include "RooMsgService.h"
39
40#include <assert.h>
41
42
43
44using namespace std;
45
47;
48
49// Register this class with RooNumIntConfig
50
51////////////////////////////////////////////////////////////////////////////////
52/// Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory
53
55{
56 RooCategory sumRule("sumRule","Summation Rule") ;
57 sumRule.defineType("Trapezoid",RooIntegrator1D::Trapezoid) ;
58 sumRule.defineType("Midpoint",RooIntegrator1D::Midpoint) ;
59 sumRule.setLabel("Trapezoid") ;
60 RooCategory extrap("extrapolation","Extrapolation procedure") ;
61 extrap.defineType("None",0) ;
62 extrap.defineType("Wynn-Epsilon",1) ;
63 extrap.setLabel("Wynn-Epsilon") ;
64 RooRealVar maxSteps("maxSteps","Maximum number of steps",20) ;
65 RooRealVar minSteps("minSteps","Minimum number of steps",999) ;
66 RooRealVar fixSteps("fixSteps","Fixed number of steps",0) ;
67
69 fact.storeProtoIntegrator(proto,RooArgSet(sumRule,extrap,maxSteps,minSteps,fixSteps)) ;
71}
72
73
74
75////////////////////////////////////////////////////////////////////////////////
76/// coverity[UNINIT_CTOR]
77/// Default constructor
78
80 _h(0), _s(0), _c(0), _d(0), _x(0)
81{
82}
83
84
85////////////////////////////////////////////////////////////////////////////////
86/// Construct integrator on given function binding, using specified summation
87/// rule, maximum number of steps and conversion tolerance. The integration
88/// limits are taken from the function binding
89
91 Int_t maxSteps, Double_t eps) :
92 RooAbsIntegrator(function), _rule(rule), _maxSteps(maxSteps), _minStepsZero(999), _fixSteps(0), _epsAbs(eps), _epsRel(eps), _doExtrap(kTRUE)
93{
96}
97
98
99////////////////////////////////////////////////////////////////////////////////
100/// Construct integrator on given function binding for given range,
101/// using specified summation rule, maximum number of steps and
102/// conversion tolerance. The integration limits are taken from the
103/// function binding
104
106 SummationRule rule, Int_t maxSteps, Double_t eps) :
108 _rule(rule),
109 _maxSteps(maxSteps),
110 _minStepsZero(999),
111 _fixSteps(0),
112 _epsAbs(eps),
113 _epsRel(eps),
114 _doExtrap(kTRUE)
115{
117 _xmin= xmin;
118 _xmax= xmax;
120}
121
122
123////////////////////////////////////////////////////////////////////////////////
124/// Construct integrator on given function binding, using specified
125/// configuration object. The integration limits are taken from the
126/// function binding
127
129 RooAbsIntegrator(function,config.printEvalCounter()),
130 _epsAbs(config.epsAbs()),
131 _epsRel(config.epsRel())
132{
133 // Extract parameters from config object
134 const RooArgSet& configSet = config.getConfigSection(IsA()->GetName()) ;
135 _rule = (SummationRule) configSet.getCatIndex("sumRule",Trapezoid) ;
136 _maxSteps = (Int_t) configSet.getRealValue("maxSteps",20) ;
137 _minStepsZero = (Int_t) configSet.getRealValue("minSteps",999) ;
138 _fixSteps = (Int_t) configSet.getRealValue("fixSteps",0) ;
139 _doExtrap = (Bool_t) configSet.getCatIndex("extrapolation",1) ;
140
141 if (_fixSteps>_maxSteps) {
142 oocoutE((TObject*)0,Integration) << "RooIntegrator1D::ctor() ERROR: fixSteps>maxSteps, fixSteps set to maxSteps" << endl ;
144 }
145
148}
149
150
151
152////////////////////////////////////////////////////////////////////////////////
153/// Construct integrator on given function binding, using specified
154/// configuration object and integration range
155
157 const RooNumIntConfig& config) :
158 RooAbsIntegrator(function,config.printEvalCounter()),
159 _epsAbs(config.epsAbs()),
160 _epsRel(config.epsRel())
161{
162 // Extract parameters from config object
163 const RooArgSet& configSet = config.getConfigSection(IsA()->GetName()) ;
164 _rule = (SummationRule) configSet.getCatIndex("sumRule",Trapezoid) ;
165 _maxSteps = (Int_t) configSet.getRealValue("maxSteps",20) ;
166 _minStepsZero = (Int_t) configSet.getRealValue("minSteps",999) ;
167 _fixSteps = (Int_t) configSet.getRealValue("fixSteps",0) ;
168 _doExtrap = (Bool_t) configSet.getCatIndex("extrapolation",1) ;
169
171 _xmin= xmin;
172 _xmax= xmax;
174}
175
176
177
178////////////////////////////////////////////////////////////////////////////////
179/// Clone integrator with new function binding and configuration. Needed by RooNumIntFactory
180
182{
183 return new RooIntegrator1D(function,config) ;
184}
185
186
187
188////////////////////////////////////////////////////////////////////////////////
189/// Initialize the integrator
190
192{
193 // apply defaults if necessary
194 if(_maxSteps <= 0) {
195 _maxSteps= (_rule == Trapezoid) ? 20 : 14;
196 }
197
198 if(_epsRel <= 0) _epsRel= 1e-6;
199 if(_epsAbs <= 0) _epsAbs= 1e-6;
200
201 // check that the integrand is a valid function
202 if(!isValid()) {
203 oocoutE((TObject*)0,Integration) << "RooIntegrator1D::initialize: cannot integrate invalid function" << endl;
204 return kFALSE;
205 }
206
207 // Allocate coordinate buffer size after number of function dimensions
209
210
211 // Allocate workspace for numerical integration engine
212 _h= new Double_t[_maxSteps + 2];
213 _s= new Double_t[_maxSteps + 2];
214 _c= new Double_t[_nPoints + 1];
215 _d= new Double_t[_nPoints + 1];
216
217 return checkLimits();
218}
219
220
221////////////////////////////////////////////////////////////////////////////////
222/// Destructor
223
225{
226 // Release integrator workspace
227 if(_h) delete[] _h;
228 if(_s) delete[] _s;
229 if(_c) delete[] _c;
230 if(_d) delete[] _d;
231 if(_x) delete[] _x;
232}
233
234
235////////////////////////////////////////////////////////////////////////////////
236/// Change our integration limits. Return kTRUE if the new limits are
237/// ok, or otherwise kFALSE. Always returns kFALSE and does nothing
238/// if this object was constructed to always use our integrand's limits.
239
241{
243 oocoutE((TObject*)0,Integration) << "RooIntegrator1D::setLimits: cannot override integrand's limits" << endl;
244 return kFALSE;
245 }
246 _xmin= *xmin;
247 _xmax= *xmax;
248 return checkLimits();
249}
250
251
252////////////////////////////////////////////////////////////////////////////////
253/// Check that our integration range is finite and otherwise return kFALSE.
254/// Update the limits from the integrand if requested.
255
257{
259 assert(0 != integrand() && integrand()->isValid());
260 const_cast<double&>(_xmin) = integrand()->getMinLimit(0);
261 const_cast<double&>(_xmax) = integrand()->getMaxLimit(0);
262 }
263 const_cast<double&>(_range) = _xmax - _xmin;
264 if (_range < 0.) {
265 oocoutE((TObject*)0,Integration) << "RooIntegrator1D::checkLimits: bad range with min > max (_xmin = " << _xmin << " _xmax = " << _xmax << ")" << endl;
266 return kFALSE;
267 }
269}
270
271
272////////////////////////////////////////////////////////////////////////////////
273/// Calculate numeric integral at given set of function binding parameters
274
276{
277 assert(isValid());
278
279 if (_range == 0.)
280 return 0.;
281
282 // Copy yvec to xvec if provided
283 if (yvec) {
284 for (UInt_t i = 0 ; i<_function->getDimension()-1 ; i++) {
285 _x[i+1] = yvec[i] ;
286 }
287 }
288
289
290 _h[1]=1.0;
291 Double_t zeroThresh = _epsAbs/_range ;
292 for(Int_t j = 1; j <= _maxSteps; ++j) {
293 // refine our estimate using the appropriate summation rule
294 _s[j]= (_rule == Trapezoid) ? addTrapezoids(j) : addMidpoints(j);
295
296 if (j >= _minStepsZero) {
297 Bool_t allZero(kTRUE) ;
298 for (int jj=0 ; jj<=j ; jj++) {
299 if (_s[j]>=zeroThresh) {
300 allZero=kFALSE ;
301 }
302 }
303 if (allZero) {
304 //cout << "Roo1DIntegrator(" << this << "): zero convergence at step " << j << ", value = " << 0 << endl ;
305 return 0;
306 }
307 }
308
309 if (_fixSteps>0) {
310
311 // Fixed step mode, return result after fixed number of steps
312 if (j==_fixSteps) {
313 //cout << "returning result at fixed step " << j << endl ;
314 return _s[j];
315 }
316
317 } else if(j >= _nPoints) {
318
319 // extrapolate the results of recent refinements and check for a stable result
320 if (_doExtrap) {
321 extrapolate(j);
322 } else {
323 _extrapValue = _s[j] ;
324 _extrapError = _s[j]-_s[j-1] ;
325 }
326
328 return _extrapValue;
329 }
330 if(fabs(_extrapError) <= _epsAbs) {
331 return _extrapValue ;
332 }
333
334 }
335 // update the step size for the next refinement of the summation
336 _h[j+1]= (_rule == Trapezoid) ? _h[j]/4. : _h[j]/9.;
337 }
338
339 oocoutW((TObject*)0,Integration) << "RooIntegrator1D::integral: integral of " << _function->getName() << " over range (" << _xmin << "," << _xmax << ") did not converge after "
340 << _maxSteps << " steps" << endl;
341 for(Int_t j = 1; j <= _maxSteps; ++j) {
342 ooccoutW((TObject*)0,Integration) << " [" << j << "] h = " << _h[j] << " , s = " << _s[j] << endl;
343 }
344
345 return _s[_maxSteps] ;
346}
347
348
349////////////////////////////////////////////////////////////////////////////////
350/// Calculate the n-th stage of refinement of the Second Euler-Maclaurin
351/// summation rule which has the useful property of not evaluating the
352/// integrand at either of its endpoints but requires more function
353/// evaluations than the trapezoidal rule. This rule can be used with
354/// a suitable change of variables to estimate improper integrals.
355
357{
358 Double_t x,tnm,sum,del,ddel;
359 Int_t it,j;
360
361 if(n == 1) {
362 Double_t xmid= 0.5*(_xmin + _xmax);
363 return (_savedResult= _range*integrand(xvec(xmid)));
364 }
365 else {
366 for(it=1, j=1; j < n-1; j++) it*= 3;
367 tnm= it;
368 del= _range/(3.*tnm);
369 ddel= del+del;
370 x= _xmin + 0.5*del;
371 for(sum= 0, j= 1; j <= it; j++) {
372 sum+= integrand(xvec(x));
373 x+= ddel;
374 sum+= integrand(xvec(x));
375 x+= del;
376 }
377 return (_savedResult= (_savedResult + _range*sum/tnm)/3.);
378 }
379}
380
381
382////////////////////////////////////////////////////////////////////////////////
383/// Calculate the n-th stage of refinement of the extended trapezoidal
384/// summation rule. This is the most efficient rule for a well behaved
385/// integrands that can be evaluated over its entire range, including the
386/// endpoints.
387
389{
390 if (n == 1) {
391 // use a single trapezoid to cover the full range
393 }
394 else {
395 // break the range down into several trapezoids using 2**(n-2)
396 // equally-spaced interior points
397 const int nInt = 1 << (n-2);
398 const double del = _range/nInt;
399 const double xmin = _xmin;
400
401 double sum = 0.;
402 // TODO Replace by batch computation
403 for (int j=0; j<nInt; ++j) {
404 double x = xmin + (0.5+j)*del;
405 sum += integrand(xvec(x));
406 }
407
408 return (_savedResult= 0.5*(_savedResult + _range*sum/nInt));
409 }
410}
411
412
413
414////////////////////////////////////////////////////////////////////////////////
415/// Extrapolate result to final value
416
418{
419 Double_t *xa= &_h[n-_nPoints];
420 Double_t *ya= &_s[n-_nPoints];
421 Int_t i,m,ns=1;
422 Double_t den,dif,dift,ho,hp,w;
423
424 dif=fabs(xa[1]);
425 for (i= 1; i <= _nPoints; i++) {
426 if ((dift=fabs(xa[i])) < dif) {
427 ns=i;
428 dif=dift;
429 }
430 _c[i]=ya[i];
431 _d[i]=ya[i];
432 }
433 _extrapValue= ya[ns--];
434 for(m= 1; m < _nPoints; m++) {
435 for(i= 1; i <= _nPoints-m; i++) {
436 ho=xa[i];
437 hp=xa[i+m];
438 w=_c[i+1]-_d[i];
439 if((den=ho-hp) == 0.0) {
440 oocoutE((TObject*)0,Integration) << "RooIntegrator1D::extrapolate: internal error" << endl;
441 }
442 den=w/den;
443 _d[i]=hp*den;
444 _c[i]=ho*den;
445 }
446 _extrapValue += (_extrapError=(2*ns < (_nPoints-m) ? _c[ns+1] : _d[ns--]));
447 }
448}
#define e(i)
Definition: RSha256.hxx:103
#define oocoutW(o, a)
Definition: RooMsgService.h:47
#define oocoutE(o, a)
Definition: RooMsgService.h:48
#define ooccoutW(o, a)
Definition: RooMsgService.h:55
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 xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
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
virtual Double_t getMinLimit(UInt_t dimension) const =0
virtual Double_t getMaxLimit(UInt_t dimension) const =0
UInt_t getDimension() const
Definition: RooAbsFunc.h:29
virtual const char * getName() const
Definition: RooAbsFunc.h:59
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
const RooAbsFunc * _function
const RooAbsFunc * integrand() const
Bool_t isValid() 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.
RooIntegrator1D implements an adaptive one-dimensional numerical integration algorithm.
Double_t addTrapezoids(Int_t n)
Calculate the n-th stage of refinement of the extended trapezoidal summation rule.
Double_t _range
Upper integration bound.
virtual Bool_t checkLimits() const
Check that our integration range is finite and otherwise return kFALSE.
virtual RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const
Clone integrator with new function binding and configuration. Needed by RooNumIntFactory.
void extrapolate(Int_t n)
Extrapolate result to final value.
Double_t * xvec(Double_t &xx)
Integrator workspace.
Bool_t setLimits(Double_t *xmin, Double_t *xmax)
Change our integration limits.
SummationRule _rule
Double_t _extrapValue
Size of integration range.
virtual ~RooIntegrator1D()
Destructor.
Bool_t initialize()
Initialize the integrator.
virtual Double_t integral(const Double_t *yvec=0)
Calculate numeric integral at given set of function binding parameters.
Double_t addMidpoints(Int_t n)
Calculate the n-th stage of refinement of the Second Euler-Maclaurin summation rule which has the use...
Double_t * _d
Integrator workspace.
Bool_t _useIntegrandLimits
static void registerIntegrator(RooNumIntFactory &fact)
Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory.
Double_t * _h
Error on extrapolated value.
Double_t _savedResult
Integrator workspace.
Double_t _extrapError
Extrapolated value.
RooIntegrator1D()
coverity[UNINIT_CTOR] Default constructor
Double_t * _c
Integrator workspace.
Double_t * _s
Integrator workspace.
Double_t _xmax
Lower integration bound.
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
RooCategory & method1D()
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...
static Int_t isInfinite(Double_t x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:58
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 x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:151
@ Integration
Definition: RooGlobalFunc.h:67
static constexpr double ns
auto * m
Definition: textangle.C:8
static long int sum(long int i)
Definition: Factory.cxx:2275