ROOT   Reference Guide
RooGaussKronrodIntegrator1D.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 *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
15 *****************************************************************************/
16
17/**
18\file RooGaussKronrodIntegrator1D.cxx
19\class RooGaussKronrodIntegrator1D
20\ingroup Roofitcore
21
22RooGaussKronrodIntegrator1D implements the Gauss-Kronrod integration algorithm.
23
24An Gaussian quadrature method for numerical integration in which
25error is estimation based on evaluation at special points known as
26"Kronrod points." By suitably picking these points, abscissas from
27previous iterations can be reused as part of the new set of points,
28whereas usual Gaussian quadrature would require recomputation of
29all abscissas at each iteration.
30
31This class automatically handles (-inf,+inf) integrals by dividing
32the integration in three regions (-inf,-1), (-1,1), (1,inf) and
33calculating the 1st and 3rd term using a x -> 1/x coordinate
34transformation
35
36This class embeds the Gauss-Kronrod integrator from the GNU
37Scientific Library version 1.5 and applies the 10-, 21-, 43- and
3887-point rule in succession until the required target precision is
39reached
40**/
41
42
43
44#include "RooFit.h"
45
46#include <assert.h>
47#include <math.h>
48#include <float.h>
49#include "Riostream.h"
50#include "TMath.h"
52#include "RooArgSet.h"
53#include "RooRealVar.h"
54#include "RooNumber.h"
55#include "RooNumIntFactory.h"
57#include "RooMsgService.h"
58
59
60using namespace std;
61
63;
64
65
66// --- From GSL_MATH.h -------------------------------------------
68{
69 double (* function) (double x, void * params);
70 void * params;
71};
73#define GSL_FN_EVAL(F,x) (*((F)->function))(x,(F)->params)
74//----From GSL_INTEGRATION.h ---------------------------------------
76 double a, double b,
77 double epsabs, double epsrel,
78 double *result, double *abserr,
79 size_t * neval);
80//-------------------------------------------------------------------
81
82// register integrator class
83// create a derived class in order to call the protected method of the
84// RoodaptiveGaussKronrodIntegrator1D
85namespace RooFit_internal {
87
88 static void registerIntegrator()
89 {
90 auto &intFactory = RooNumIntFactory::instance();
92 }
93};
94// class used to register integrator at loafing time
97};
98
100}
101
102
103////////////////////////////////////////////////////////////////////////////////
104/// Register RooGaussKronrodIntegrator1D, its parameters and capabilities with RooNumIntConfig
105
107{
109 oocoutI((TObject*)nullptr,Integration) << "RooGaussKronrodIntegrator1D has been registered" << std::endl;
110}
111
112
113
114////////////////////////////////////////////////////////////////////////////////
115/// coverity[UNINIT_CTOR]
116/// Default constructor
117
119{
120}
121
122
123
124////////////////////////////////////////////////////////////////////////////////
125/// Construct integral on 'function' using given configuration object. The integration
126/// range is taken from the definition in the function binding
127
130 _epsAbs(config.epsRel()),
131 _epsRel(config.epsAbs())
132{
135}
136
137
138
139////////////////////////////////////////////////////////////////////////////////
140/// Construct integral on 'function' using given configuration object in the given range
141
143 Double_t xmin, Double_t xmax, const RooNumIntConfig& config) :
145 _epsAbs(config.epsRel()),
146 _epsRel(config.epsAbs()),
147 _xmin(xmin),
148 _xmax(xmax)
149{
152}
153
154
155
156////////////////////////////////////////////////////////////////////////////////
157/// Clone integrator with given function and configuration. Needed for RooNumIntFactory
158
160{
161 return new RooGaussKronrodIntegrator1D(function,config) ;
162}
163
164
165
166////////////////////////////////////////////////////////////////////////////////
167/// Perform one-time initialization of integrator
168
170{
171 // Allocate coordinate buffer size after number of function dimensions
173
174 return checkLimits();
175}
176
177
178
179////////////////////////////////////////////////////////////////////////////////
180/// Destructor
181
183{
184 if (_x) {
185 delete[] _x ;
186 }
187}
188
189
190
191////////////////////////////////////////////////////////////////////////////////
192/// Change our integration limits. Return kTRUE if the new limits are
193/// ok, or otherwise kFALSE. Always returns kFALSE and does nothing
194/// if this object was constructed to always use our integrand's limits.
195
197{
199 oocoutE((TObject*)0,Eval) << "RooGaussKronrodIntegrator1D::setLimits: cannot override integrand's limits" << endl;
200 return kFALSE;
201 }
202 _xmin= *xmin;
203 _xmax= *xmax;
204 return checkLimits();
205}
206
207
208
209////////////////////////////////////////////////////////////////////////////////
210/// Check that our integration range is finite and otherwise return kFALSE.
211/// Update the limits from the integrand if requested.
212
214{
216 assert(0 != integrand() && integrand()->isValid());
219 }
220 return kTRUE ;
221}
222
223
224
226{
228 return instance->integrand(instance->xvec(x)) ;
229}
230
231
232
233////////////////////////////////////////////////////////////////////////////////
234/// Calculate and return integral
235
237{
238 assert(isValid());
239
240 // Copy yvec to xvec if provided
241 if (yvec) {
242 UInt_t i ; for (i=0 ; i<_function->getDimension()-1 ; i++) {
243 _x[i+1] = yvec[i] ;
244 }
245 }
246
247 // Setup glue function
250 F.params = this ;
251
252 // Return values
253 double result, error;
254 size_t neval = 0 ;
255
256 // Call GSL implementation of integeator
257 gsl_integration_qng (&F, _xmin, _xmax, _epsAbs, _epsRel, &result, &error, &neval);
258
259 return result;
260}
261
262
263// ----------------------------------------------------------------------------
264// ---------- Code below imported from GSL version 1.5 ------------------------
265// ----------------------------------------------------------------------------
266
267/*
268 *
269 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
270 *
271 * This program is free software; you can redistribute it and/or modify
273 * the Free Software Foundation; either version 2 of the License, or (at
274 * your option) any later version.
275 *
276 * This program is distributed in the hope that it will be useful, but
277 * WITHOUT ANY WARRANTY; without even the implied warranty of
278 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
279 * General Public License for more details.
280 *
281 * You should have received a copy of the GNU General Public License
282 * along with this program; if not, write to the Free Software
283 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
284 */
285
286#define GSL_SUCCESS 0
287#define GSL_EBADTOL 13 /* user specified an invalid tolerance */
288#define GSL_ETOL 14 /* failed to reach the specified tolerance */
289#define GSL_ERROR(a,b) oocoutE((TObject*)0,Eval) << "RooGaussKronrodIntegrator1D::integral() ERROR: " << a << endl ; return b ;
290#define GSL_DBL_MIN 2.2250738585072014e-308
291#define GSL_DBL_EPSILON 2.2204460492503131e-16
292
293
294// INCLUDED BELOW #include "qng.c"
295
297 double a, double b,
298 double epsabs, double epsrel,
299 double *result, double *abserr,
300 size_t * neval);
301
302
303
304// INCLUDED BELOW #include "err.c"
305static double rescale_error (double err, const double result_abs, const double result_asc) ;
306
307static double
308rescale_error (double err, const double result_abs, const double result_asc)
309{
310 err = fabs(err) ;
311
312 if (result_asc != 0 && err != 0)
313 {
314 double scale = TMath::Power((200 * err / result_asc), 1.5) ;
315
316 if (scale < 1)
317 {
318 err = result_asc * scale ;
319 }
320 else
321 {
322 err = result_asc ;
323 }
324 }
325 if (result_abs > GSL_DBL_MIN / (50 * GSL_DBL_EPSILON))
326 {
327 double min_err = 50 * GSL_DBL_EPSILON * result_abs ;
328
329 if (min_err > err)
330 {
331 err = min_err ;
332 }
333 }
334
335 return err ;
336}
337
338
339// INCLUDED BELOW #include "qng.h"
340/* Gauss-Kronrod-Patterson quadrature coefficients for use in
341 quadpack routine qng. These coefficients were calculated with
342 101 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov
343 1981. */
344
345/* x1, abscissae common to the 10-, 21-, 43- and 87-point rule */
346static const double x1[5] = {
347 0.973906528517171720077964012084452,
348 0.865063366688984510732096688423493,
349 0.679409568299024406234327365114874,
350 0.433395394129247190799265943165784,
351 0.148874338981631210884826001129720
352} ;
353
354/* w10, weights of the 10-point formula */
355static const double w10[5] = {
356 0.066671344308688137593568809893332,
357 0.149451349150580593145776339657697,
358 0.219086362515982043995534934228163,
359 0.269266719309996355091226921569469,
360 0.295524224714752870173892994651338
361} ;
362
363/* x2, abscissae common to the 21-, 43- and 87-point rule */
364static const double x2[5] = {
365 0.995657163025808080735527280689003,
366 0.930157491355708226001207180059508,
367 0.780817726586416897063717578345042,
368 0.562757134668604683339000099272694,
369 0.294392862701460198131126603103866
370} ;
371
372/* w21a, weights of the 21-point formula for abscissae x1 */
373static const double w21a[5] = {
374 0.032558162307964727478818972459390,
375 0.075039674810919952767043140916190,
376 0.109387158802297641899210590325805,
377 0.134709217311473325928054001771707,
378 0.147739104901338491374841515972068
379} ;
380
381/* w21b, weights of the 21-point formula for abscissae x2 */
382static const double w21b[6] = {
383 0.011694638867371874278064396062192,
384 0.054755896574351996031381300244580,
385 0.093125454583697605535065465083366,
386 0.123491976262065851077958109831074,
387 0.142775938577060080797094273138717,
388 0.149445554002916905664936468389821
389} ;
390
391/* x3, abscissae common to the 43- and 87-point rule */
392static const double x3[11] = {
393 0.999333360901932081394099323919911,
394 0.987433402908088869795961478381209,
395 0.954807934814266299257919200290473,
396 0.900148695748328293625099494069092,
397 0.825198314983114150847066732588520,
398 0.732148388989304982612354848755461,
399 0.622847970537725238641159120344323,
400 0.499479574071056499952214885499755,
401 0.364901661346580768043989548502644,
402 0.222254919776601296498260928066212,
403 0.074650617461383322043914435796506
404} ;
405
406/* w43a, weights of the 43-point formula for abscissae x1, x3 */
407static const double w43a[10] = {
408 0.016296734289666564924281974617663,
409 0.037522876120869501461613795898115,
410 0.054694902058255442147212685465005,
411 0.067355414609478086075553166302174,
412 0.073870199632393953432140695251367,
413 0.005768556059769796184184327908655,
414 0.027371890593248842081276069289151,
415 0.046560826910428830743339154433824,
416 0.061744995201442564496240336030883,
417 0.071387267268693397768559114425516
418} ;
419
420/* w43b, weights of the 43-point formula for abscissae x3 */
421static const double w43b[12] = {
422 0.001844477640212414100389106552965,
423 0.010798689585891651740465406741293,
424 0.021895363867795428102523123075149,
425 0.032597463975345689443882222526137,
426 0.042163137935191811847627924327955,
427 0.050741939600184577780189020092084,
428 0.058379395542619248375475369330206,
429 0.064746404951445885544689259517511,
430 0.069566197912356484528633315038405,
431 0.072824441471833208150939535192842,
432 0.074507751014175118273571813842889,
433 0.074722147517403005594425168280423
434} ;
435
436/* x4, abscissae of the 87-point rule */
437static const double x4[22] = {
438 0.999902977262729234490529830591582,
439 0.997989895986678745427496322365960,
440 0.992175497860687222808523352251425,
441 0.981358163572712773571916941623894,
442 0.965057623858384619128284110607926,
443 0.943167613133670596816416634507426,
444 0.915806414685507209591826430720050,
445 0.883221657771316501372117548744163,
446 0.845710748462415666605902011504855,
447 0.803557658035230982788739474980964,
448 0.757005730685495558328942793432020,
449 0.706273209787321819824094274740840,
450 0.651589466501177922534422205016736,
451 0.593223374057961088875273770349144,
452 0.531493605970831932285268948562671,
453 0.466763623042022844871966781659270,
454 0.399424847859218804732101665817923,
455 0.329874877106188288265053371824597,
456 0.258503559202161551802280975429025,
457 0.185695396568346652015917141167606,
458 0.111842213179907468172398359241362,
459 0.037352123394619870814998165437704
460} ;
461
462/* w87a, weights of the 87-point formula for abscissae x1, x2, x3 */
463static const double w87a[21] = {
464 0.008148377384149172900002878448190,
465 0.018761438201562822243935059003794,
466 0.027347451050052286161582829741283,
467 0.033677707311637930046581056957588,
468 0.036935099820427907614589586742499,
469 0.002884872430211530501334156248695,
470 0.013685946022712701888950035273128,
471 0.023280413502888311123409291030404,
472 0.030872497611713358675466394126442,
473 0.035693633639418770719351355457044,
474 0.000915283345202241360843392549948,
475 0.005399280219300471367738743391053,
476 0.010947679601118931134327826856808,
477 0.016298731696787335262665703223280,
478 0.021081568889203835112433060188190,
479 0.025370969769253827243467999831710,
480 0.029189697756475752501446154084920,
481 0.032373202467202789685788194889595,
482 0.034783098950365142750781997949596,
483 0.036412220731351787562801163687577,
484 0.037253875503047708539592001191226
485} ;
486
487/* w87b, weights of the 87-point formula for abscissae x4 */
488static const double w87b[23] = {
489 0.000274145563762072350016527092881,
490 0.001807124155057942948341311753254,
491 0.004096869282759164864458070683480,
492 0.006758290051847378699816577897424,
493 0.009549957672201646536053581325377,
494 0.012329447652244853694626639963780,
495 0.015010447346388952376697286041943,
496 0.017548967986243191099665352925900,
497 0.019938037786440888202278192730714,
498 0.022194935961012286796332102959499,
499 0.024339147126000805470360647041454,
500 0.026374505414839207241503786552615,
501 0.028286910788771200659968002987960,
502 0.030052581128092695322521110347341,
503 0.031646751371439929404586051078883,
504 0.033050413419978503290785944862689,
505 0.034255099704226061787082821046821,
506 0.035262412660156681033782717998428,
507 0.036076989622888701185500318003895,
508 0.036698604498456094498018047441094,
509 0.037120549269832576114119958413599,
510 0.037334228751935040321235449094698,
511 0.037361073762679023410321241766599
512} ;
513
514
515int
517 double a, double b,
518 double epsabs, double epsrel,
519 double * result, double * abserr, size_t * neval)
520{
521 double fv1[5], fv2[5], fv3[5], fv4[5];
522 double savfun[21]; /* array of function values which have been computed */
523 double res10, res21, res43, res87; /* 10, 21, 43 and 87 point results */
524 double result_kronrod, err ;
525 double resabs; /* approximation to the integral of abs(f) */
526 double resasc; /* approximation to the integral of abs(f-i/(b-a)) */
527
528 const double half_length = 0.5 * (b - a);
529 const double abs_half_length = fabs (half_length);
530 const double center = 0.5 * (b + a);
531 const double f_center = GSL_FN_EVAL(f, center);
532
533 int k ;
534
535 if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
536 {
537 * result = 0;
538 * abserr = 0;
539 * neval = 0;
540 GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
542 };
543
544 /* Compute the integral using the 10- and 21-point formula. */
545
546 res10 = 0;
547 res21 = w21b[5] * f_center;
548 resabs = w21b[5] * fabs (f_center);
549
550 for (k = 0; k < 5; k++)
551 {
552 const double abscissa = half_length * x1[k];
553 const double fval1 = GSL_FN_EVAL(f, center + abscissa);
554 const double fval2 = GSL_FN_EVAL(f, center - abscissa);
555 const double fval = fval1 + fval2;
556 res10 += w10[k] * fval;
557 res21 += w21a[k] * fval;
558 resabs += w21a[k] * (fabs (fval1) + fabs (fval2));
559 savfun[k] = fval;
560 fv1[k] = fval1;
561 fv2[k] = fval2;
562 }
563
564 for (k = 0; k < 5; k++)
565 {
566 const double abscissa = half_length * x2[k];
567 const double fval1 = GSL_FN_EVAL(f, center + abscissa);
568 const double fval2 = GSL_FN_EVAL(f, center - abscissa);
569 const double fval = fval1 + fval2;
570 res21 += w21b[k] * fval;
571 resabs += w21b[k] * (fabs (fval1) + fabs (fval2));
572 savfun[k + 5] = fval;
573 fv3[k] = fval1;
574 fv4[k] = fval2;
575 }
576
577 resabs *= abs_half_length ;
578
579 {
580 const double mean = 0.5 * res21;
581
582 resasc = w21b[5] * fabs (f_center - mean);
583
584 for (k = 0; k < 5; k++)
585 {
586 resasc +=
587 (w21a[k] * (fabs (fv1[k] - mean) + fabs (fv2[k] - mean))
588 + w21b[k] * (fabs (fv3[k] - mean) + fabs (fv4[k] - mean)));
589 }
590 resasc *= abs_half_length ;
591 }
592
593 result_kronrod = res21 * half_length;
594
595 err = rescale_error ((res21 - res10) * half_length, resabs, resasc) ;
596
597 /* test for convergence. */
598
599 if (err < epsabs || err < epsrel * fabs (result_kronrod))
600 {
601 * result = result_kronrod ;
602 * abserr = err ;
603 * neval = 21;
604 return GSL_SUCCESS;
605 }
606
607 /* compute the integral using the 43-point formula. */
608
609 res43 = w43b[11] * f_center;
610
611 for (k = 0; k < 10; k++)
612 {
613 res43 += savfun[k] * w43a[k];
614 }
615
616 for (k = 0; k < 11; k++)
617 {
618 const double abscissa = half_length * x3[k];
619 const double fval = (GSL_FN_EVAL(f, center + abscissa)
620 + GSL_FN_EVAL(f, center - abscissa));
621 res43 += fval * w43b[k];
622 savfun[k + 10] = fval;
623 }
624
625 /* test for convergence */
626
627 result_kronrod = res43 * half_length;
628 err = rescale_error ((res43 - res21) * half_length, resabs, resasc);
629
630 if (err < epsabs || err < epsrel * fabs (result_kronrod))
631 {
632 * result = result_kronrod ;
633 * abserr = err ;
634 * neval = 43;
635 return GSL_SUCCESS;
636 }
637
638 /* compute the integral using the 87-point formula. */
639
640 res87 = w87b[22] * f_center;
641
642 for (k = 0; k < 21; k++)
643 {
644 res87 += savfun[k] * w87a[k];
645 }
646
647 for (k = 0; k < 22; k++)
648 {
649 const double abscissa = half_length * x4[k];
650 res87 += w87b[k] * (GSL_FN_EVAL(f, center + abscissa)
651 + GSL_FN_EVAL(f, center - abscissa));
652 }
653
654 /* test for convergence */
655
656 result_kronrod = res87 * half_length ;
657
658 err = rescale_error ((res87 - res43) * half_length, resabs, resasc);
659
660 if (err < epsabs || err < epsrel * fabs (result_kronrod))
661 {
662 * result = result_kronrod ;
663 * abserr = err ;
664 * neval = 87;
665 return GSL_SUCCESS;
666 }
667
668 /* failed to converge */
669
670 * result = result_kronrod ;
671 * abserr = err ;
672 * neval = 87;
673
674 // GSL_ERROR("failed to reach tolerance with highest-order rule", GSL_ETOL) ;
675 return GSL_ETOL ;
676}
double
Definition: Converters.cxx:939
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define GSL_ERROR(a, b)
static const double x2[5]
static double rescale_error(double err, const double result_abs, const double result_asc)
static const double w10[5]
#define GSL_DBL_EPSILON
static const double w43a[10]
static const double x4[22]
static const double w21b[6]
static const double x1[5]
double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
static const double w87a[21]
static const double w21a[5]
static const double w43b[12]
static const double x3[11]
#define GSL_SUCCESS
#define GSL_FN_EVAL(F, x)
static const double w87b[23]
int gsl_integration_qng(const gsl_function *f, double a, double b, double epsabs, double epsrel, double *result, double *abserr, size_t *neval)
#define GSL_DBL_MIN
#define oocoutE(o, a)
Definition: RooMsgService.h:48
#define oocoutI(o, a)
Definition: RooMsgService.h:45
unsigned int UInt_t
Definition: RtypesCore.h:46
const Bool_t kFALSE
Definition: RtypesCore.h:101
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:364
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
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
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:35
RooGaussKronrodIntegrator1D implements the Gauss-Kronrod integration algorithm.
friend double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
Bool_t initialize()
Perform one-time initialization of integrator.
RooGaussKronrodIntegrator1D()
coverity[UNINIT_CTOR] Default constructor
Double_t _xmax
Lower integration bound.
virtual RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const
Clone integrator with given function and configuration. Needed for RooNumIntFactory.
virtual Bool_t checkLimits() const
Check that our integration range is finite and otherwise return kFALSE.
virtual Double_t integral(const Double_t *yvec=0)
Calculate and return integral.
Bool_t setLimits(Double_t *xmin, Double_t *xmax)
Change our integration limits.
static void registerIntegrator(RooNumIntFactory &fact)
Register RooGaussKronrodIntegrator1D, its parameters and capabilities with RooNumIntConfig.
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
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 RooNumIntFactory & instance()
Static method returning reference to singleton instance of factory.
Mother of all ROOT objects.
Definition: TObject.h:37
Double_t x[n]
Definition: legend1.C:17
#define F(x, y, z)
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:152
static Roo_reg_GKInteg1D instance
static Roo_reg_AGKInteg1D instance
@ Integration
Definition: RooGlobalFunc.h:60
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:735
double(* function)(double x, void *params)
auto * a
Definition: textangle.C:12