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