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