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