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