GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GIntegral.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GIntegral.cpp - Integration class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2015 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GIntegral.cpp
23  * @brief Integration class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #include <cmath> // std::abs()
29 #include <vector>
30 #include <algorithm> // std::sort
31 #include "GIntegral.hpp"
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 
35 /* __ Method name definitions ____________________________________________ */
36 #define G_ROMBERG "GIntegral::romberg(double&, double&, int&)"
37 #define G_TRAPZD "GIntegral::trapzd(double&, double&, int&, double)"
38 #define G_POLINT "GIntegral::polint(double*, double*, int, double, double*)"
39 
40 /* __ Macros _____________________________________________________________ */
41 
42 /* __ Coding definitions _________________________________________________ */
43 
44 /* __ Debug definitions __________________________________________________ */
45 
46 /* __ Constants __________________________________________________________ */
47 namespace gammalib {
48 
49  // Gauss-Kronrod abscissae, common to the 10-, 21-, 43- and 87-point rule
50  const double gkx1[5] = {
51  0.973906528517171720077964012084452,
52  0.865063366688984510732096688423493,
53  0.679409568299024406234327365114874,
54  0.433395394129247190799265943165784,
55  0.148874338981631210884826001129720
56  };
57 
58  // Gauss-Kronrod weights of the 10-point rule
59  const double gkw10[5] = {
60  0.066671344308688137593568809893332,
61  0.149451349150580593145776339657697,
62  0.219086362515982043995534934228163,
63  0.269266719309996355091226921569469,
64  0.295524224714752870173892994651338
65  };
66 
67  // Gauss-Kronrod abscissae, common to the 21-, 43- and 87-point rule
68  const double gkx2[5] = {
69  0.995657163025808080735527280689003,
70  0.930157491355708226001207180059508,
71  0.780817726586416897063717578345042,
72  0.562757134668604683339000099272694,
73  0.294392862701460198131126603103866
74  };
75 
76  // Gauss-Kronrod weights of the 21-point rule for abscissae gkx1
77  const double gkw21a[5] = {
78  0.032558162307964727478818972459390,
79  0.075039674810919952767043140916190,
80  0.109387158802297641899210590325805,
81  0.134709217311473325928054001771707,
82  0.147739104901338491374841515972068
83  };
84 
85  // Gauss-Kronrod weights of the 21-point rule for abscissae gkx2
86  const double gkw21b[6] = {
87  0.011694638867371874278064396062192,
88  0.054755896574351996031381300244580,
89  0.093125454583697605535065465083366,
90  0.123491976262065851077958109831074,
91  0.142775938577060080797094273138717,
92  0.149445554002916905664936468389821
93  };
94 
95  // Gauss-Kronrod abscissae, common to the 43- and 87-point rule
96  const double gkx3[11] = {
97  0.999333360901932081394099323919911,
98  0.987433402908088869795961478381209,
99  0.954807934814266299257919200290473,
100  0.900148695748328293625099494069092,
101  0.825198314983114150847066732588520,
102  0.732148388989304982612354848755461,
103  0.622847970537725238641159120344323,
104  0.499479574071056499952214885499755,
105  0.364901661346580768043989548502644,
106  0.222254919776601296498260928066212,
107  0.074650617461383322043914435796506
108  };
109 
110  // Gauss-Kronrod weights of the 43-point rule for abscissae gkx1, gkx3
111  const double gkw43a[10] = {
112  0.016296734289666564924281974617663,
113  0.037522876120869501461613795898115,
114  0.054694902058255442147212685465005,
115  0.067355414609478086075553166302174,
116  0.073870199632393953432140695251367,
117  0.005768556059769796184184327908655,
118  0.027371890593248842081276069289151,
119  0.046560826910428830743339154433824,
120  0.061744995201442564496240336030883,
121  0.071387267268693397768559114425516
122  };
123 
124  // Gauss-Kronrod weights of the 43-point formula for abscissae gkx3
125  const double gkw43b[12] = {
126  0.001844477640212414100389106552965,
127  0.010798689585891651740465406741293,
128  0.021895363867795428102523123075149,
129  0.032597463975345689443882222526137,
130  0.042163137935191811847627924327955,
131  0.050741939600184577780189020092084,
132  0.058379395542619248375475369330206,
133  0.064746404951445885544689259517511,
134  0.069566197912356484528633315038405,
135  0.072824441471833208150939535192842,
136  0.074507751014175118273571813842889,
137  0.074722147517403005594425168280423
138  };
139 
140  // Gauss-Kronrod abscissae, of the 87-point rule
141  const double gkx4[22] = {
142  0.999902977262729234490529830591582,
143  0.997989895986678745427496322365960,
144  0.992175497860687222808523352251425,
145  0.981358163572712773571916941623894,
146  0.965057623858384619128284110607926,
147  0.943167613133670596816416634507426,
148  0.915806414685507209591826430720050,
149  0.883221657771316501372117548744163,
150  0.845710748462415666605902011504855,
151  0.803557658035230982788739474980964,
152  0.757005730685495558328942793432020,
153  0.706273209787321819824094274740840,
154  0.651589466501177922534422205016736,
155  0.593223374057961088875273770349144,
156  0.531493605970831932285268948562671,
157  0.466763623042022844871966781659270,
158  0.399424847859218804732101665817923,
159  0.329874877106188288265053371824597,
160  0.258503559202161551802280975429025,
161  0.185695396568346652015917141167606,
162  0.111842213179907468172398359241362,
163  0.037352123394619870814998165437704
164  };
165 
166  // Gauss-Kronrod weights of the 87-point rule for abscissae gkx1, gkx2, gkx3
167  const double gkw87a[21] = {
168  0.008148377384149172900002878448190,
169  0.018761438201562822243935059003794,
170  0.027347451050052286161582829741283,
171  0.033677707311637930046581056957588,
172  0.036935099820427907614589586742499,
173  0.002884872430211530501334156248695,
174  0.013685946022712701888950035273128,
175  0.023280413502888311123409291030404,
176  0.030872497611713358675466394126442,
177  0.035693633639418770719351355457044,
178  0.000915283345202241360843392549948,
179  0.005399280219300471367738743391053,
180  0.010947679601118931134327826856808,
181  0.016298731696787335262665703223280,
182  0.021081568889203835112433060188190,
183  0.025370969769253827243467999831710,
184  0.029189697756475752501446154084920,
185  0.032373202467202789685788194889595,
186  0.034783098950365142750781997949596,
187  0.036412220731351787562801163687577,
188  0.037253875503047708539592001191226
189  };
190 
191  // Gauss-Kronrod weights of the 87-point formula for abscissae gkx4
192  const double gkw87b[23] = {
193  0.000274145563762072350016527092881,
194  0.001807124155057942948341311753254,
195  0.004096869282759164864458070683480,
196  0.006758290051847378699816577897424,
197  0.009549957672201646536053581325377,
198  0.012329447652244853694626639963780,
199  0.015010447346388952376697286041943,
200  0.017548967986243191099665352925900,
201  0.019938037786440888202278192730714,
202  0.022194935961012286796332102959499,
203  0.024339147126000805470360647041454,
204  0.026374505414839207241503786552615,
205  0.028286910788771200659968002987960,
206  0.030052581128092695322521110347341,
207  0.031646751371439929404586051078883,
208  0.033050413419978503290785944862689,
209  0.034255099704226061787082821046821,
210  0.035262412660156681033782717998428,
211  0.036076989622888701185500318003895,
212  0.036698604498456094498018047441094,
213  0.037120549269832576114119958413599,
214  0.037334228751935040321235449094698,
215  0.037361073762679023410321241766599
216  };
217 
218 } // end gammalib namespace
219 
220 
221 /*==========================================================================
222  = =
223  = Constructors/destructors =
224  = =
225  ==========================================================================*/
226 
227 /***********************************************************************//**
228  * @brief Void constructor
229  ***************************************************************************/
231 {
232  // Initialise members
233  init_members();
234 
235  // Return
236  return;
237 }
238 
239 
240 /***********************************************************************//**
241  * @brief Function kernel constructor
242  *
243  * @param[in] kernel Pointer to function kernel.
244  *
245  * The function kernel constructor assigns the function kernel pointer in
246  * constructing the object.
247  ***************************************************************************/
249 {
250  // Initialise members
251  init_members();
252 
253  // Set function kernel
254  m_kernel = kernel;
255 
256  // Return
257  return;
258 }
259 
260 
261 /***********************************************************************//**
262  * @brief Copy constructor
263  *
264  * @param[in] integral Integral.
265  ***************************************************************************/
267 {
268  // Initialise members
269  init_members();
270 
271  // Copy members
272  copy_members(integral);
273 
274  // Return
275  return;
276 }
277 
278 
279 /***********************************************************************//**
280  * @brief Destructor
281  ***************************************************************************/
283 {
284  // Free members
285  free_members();
286 
287  // Return
288  return;
289 }
290 
291 
292 /*==========================================================================
293  = =
294  = Operators =
295  = =
296  ==========================================================================*/
297 
298 /***********************************************************************//**
299  * @brief Assignment operator
300  *
301  * @param[in] integral Integral.
302  * @return Integral.
303  ***************************************************************************/
305 {
306  // Execute only if object is not identical
307  if (this != &integral) {
308 
309  // Free members
310  free_members();
311 
312  // Initialise integral
313  init_members();
314 
315  // Copy members
316  copy_members(integral);
317 
318  } // endif: object was not identical
319 
320  // Return
321  return *this;
322 }
323 
324 
325 /*==========================================================================
326  = =
327  = Public methods =
328  = =
329  ==========================================================================*/
330 
331 /***********************************************************************//**
332  * @brief Clear integral
333  ***************************************************************************/
335 {
336  // Free members
337  free_members();
338 
339  // Initialise private members
340  init_members();
341 
342  // Return
343  return;
344 }
345 
346 
347 /***********************************************************************//**
348  * @brief Clone integral
349  *
350  * @return Pointer to deep copy of integral.
351  ***************************************************************************/
353 {
354  return new GIntegral(*this);
355 }
356 
357 
358 /***********************************************************************//**
359  * @brief Perform Romberg integration
360  *
361  * @param[in] bounds Integration boundaries.
362  * @param[in] order Integration order (default: 5)
363  *
364  * @exception GException::invalid_argument
365  * Integration order incompatible with number of iterations.
366  *
367  * Returns the integral of the integrand, computed over a number of
368  * intervals [a0,a1], [a1,a2], ... that are given as an unordered vector
369  * by the @p bounds argument.
370  *
371  * Integration is performed by Romberg's method of order 2*order, where
372  *
373  * order=1 is equivalent to the trapezoidal rule,
374  * order=2 is equivalent to Simpson's rule, and
375  * order=3 is equivalent to Boole's rule.
376  *
377  * The number of iterations is limited by m_max_iter. m_eps specifies the
378  * requested fractional accuracy. By default it is set to 1e-6.
379  ***************************************************************************/
380 double GIntegral::romberg(std::vector<double> bounds, const int& order)
381 {
382  // Sort integration boundaries in ascending order
383  std::sort(bounds.begin(), bounds.end());
384 
385  // Initialise integral
386  double value = 0.0;
387 
388  // Initialise integration status information
389  int calls = 0;
390 
391  // Add integral of all intervals
392  for (int i = 0; i < bounds.size()-1; ++i) {
393  value += romberg(bounds[i], bounds[i+1], order);
394  calls += m_calls;
395  }
396 
397  // Set integration status information
398  m_calls = calls;
399 
400  // Return value
401  return value;
402 }
403 
404 
405 /***********************************************************************//**
406  * @brief Perform Romberg integration
407  *
408  * @param[in] a Left integration boundary.
409  * @param[in] b Right integration boundary.
410  * @param[in] order Integration order (default: 5)
411  *
412  * @exception GException::invalid_argument
413  * Integration order incompatible with number of iterations.
414  *
415  * Returns the integral of the integrand from a to b. Integration is
416  * performed by Romberg's method of order 2*order, where
417  *
418  * order=1 is equivalent to the trapezoidal rule,
419  * order=2 is equivalent to Simpson's rule, and
420  * order=3 is equivalent to Boole's rule.
421  *
422  * The number of iterations is limited by m_max_iter. m_eps specifies the
423  * requested fractional accuracy. By default it is set to 1e-6.
424  ***************************************************************************/
425 double GIntegral::romberg(const double& a, const double& b, const int& order)
426 {
427  // Initialise result and status
428  double result = 0.0;
429 
430  // Initialise integration status information
431  m_isvalid = true;
432  m_calls = 0;
433  m_has_abserr = false;
434  m_relerr = false;
435 
436  // Continue only if integration range is valid
437  if (b > a) {
438 
439  // Initialise variables
440  bool converged = false;
441  double dss = 0.0;
442 
443  // Determine (maximum) number of iterations
444  int max_iter = (m_fix_iter > 0) ? m_fix_iter : m_max_iter;
445 
446  // Check whether maximum number of iterations is compliant with
447  // order
448  if (order > max_iter) {
449  std::string msg = "Requested integration order "+
450  gammalib::str(order)+" is larger than the "
451  "maximum number of iterations "+
452  gammalib::str(max_iter)+". Either reduce the "
453  "integration order or increase the (maximum) "
454  "number of iterations.";
456  }
457 
458  // Allocate temporal storage
459  double* s = new double[max_iter+2];
460  double* h = new double[max_iter+2];
461 
462  // Initialise step size
463  h[1] = 1.0;
464  s[0] = 0.0;
465 
466  // Iterative loop
467  for (m_iter = 1; m_iter <= max_iter; ++m_iter) {
468 
469  // Integration using Trapezoid rule
470  s[m_iter] = trapzd(a, b, m_iter, s[m_iter-1]);
471 
472  // Compile option: Check for NaN/Inf
473  #if defined(G_NAN_CHECK)
474  if (is_notanumber(s[m_iter]) || is_infinite(s[m_iter])) {
475  m_message = "*** ERROR: GIntegral::romberg"
476  "(a="+gammalib::str(a)+", b="+gammalib::str(b)+""
477  ", k="+gammalib::str(k)+"): NaN/Inf encountered"
478  " (s["+gammalib::str(m_iter)+"]="
479  ""+gammalib::str(s[m_iter])+")";
480  std::cout << m_message << std::endl;
481  m_isvalid = false;
482  }
483  #endif
484 
485  // Starting from iteration order on, use polynomial interpolation
486  if (m_iter >= order) {
487 
488  // Compute result using polynom interpolation
489  result = polint(&h[m_iter-order], &s[m_iter-order],
490  order, 0.0, &dss);
491 
492  // If a fixed number of iterations has been requested then
493  // check whether we reached the final one; otherwise check
494  // whether we reached the requested precision.
495  if (m_fix_iter > 0) {
496  if (m_iter == max_iter) {
497  converged = true;
498  }
499  }
500  else if (std::abs(dss) <= m_eps * std::abs(result)) {
501  converged = true;
502  }
503  if (converged) {
504  m_has_abserr = true;
505  m_abserr = std::abs(dss);
506  if (std::abs(result) > 0) {
507  m_has_relerr = true;
508  m_relerr = m_abserr / std::abs(result);
509  }
510  break;
511  }
512 
513  } // endif: polynomial interpolation performed
514 
515  // Reduce step size
516  h[m_iter+1]= 0.25 * h[m_iter];
517 
518  } // endfor: iterative loop
519 
520  // Free temporal storage
521  delete [] s;
522  delete [] h;
523 
524  // Set status and optionally dump warning
525  if (!converged) {
526  m_isvalid = false;
527  m_message = "Integration uncertainty "+
528  gammalib::str(std::abs(dss))+
529  " exceeds absolute tolerance of "+
530  gammalib::str(m_eps * std::abs(result))+
531  " after "+gammalib::str(m_iter)+
532  " iterations. Result "+
533  gammalib::str(result)+
534  " is inaccurate.";
535  if (!m_silent) {
536  std::string origin = "GIntegral::romberg("+
537  gammalib::str(a)+", "+
538  gammalib::str(b)+", "+
539  gammalib::str(order)+")";
540  gammalib::warning(origin, m_message);
541  }
542  }
543 
544  } // endif: integration range was valid
545 
546  // Return result
547  return result;
548 }
549 
550 
551 /***********************************************************************//**
552  * @brief Perform Trapezoidal integration
553  *
554  * @param[in] a Left integration boundary.
555  * @param[in] b Right integration boundary.
556  * @param[in] n Number of steps.
557  * @param[in] result Result from a previous trapezoidal integration step.
558  *
559  * The original Numerical Recipes function had result declared as a static
560  * variable, yet this led to some untrackable integration problems. For this
561  * reason, previous results are now passed using an argument.
562  * Result initialisation is done if n=1.
563  ***************************************************************************/
564 double GIntegral::trapzd(const double& a, const double& b, const int& n,
565  double result)
566 {
567  // Handle case of identical boundaries
568  if (a == b) {
569  result = 0.0;
570  }
571 
572  // ... otherwise use trapeziodal rule
573  else {
574 
575  // Case A: Only a single step is requested
576  if (n == 1) {
577 
578  // Evaluate integrand at boundaries
579  double y_a = m_kernel->eval(a);
580  double y_b = m_kernel->eval(b);
581  m_calls += 2;
582 
583  // Compute result
584  result = 0.5*(b-a)*(y_a + y_b);
585 
586  } // endif: only a single step was requested
587 
588  // Case B: More than a single step is requested
589  else {
590 
591  // Compute step level 2^(n-1)
592  int it = 1;
593  for (int j = 1; j < n-1; ++j) {
594  it <<= 1;
595  }
596 
597  // Verify that step level is valid
598  if (it == 0) {
599  m_isvalid = false;
600  m_message = "Invalid step level "+gammalib::str(it)+
601  " encountered for"
602  " a="+gammalib::str(a)+
603  ", b="+gammalib::str(b)+
604  ", n="+gammalib::str(n)+
605  ", result="+gammalib::str(result)+
606  ". Looks like n is too large.";
608  }
609 
610  // Set step size
611  double tnm = double(it);
612  double del = (b-a)/tnm;
613 
614  // Verify that step is >0
615  if (del == 0) {
616  m_isvalid = false;
617  m_message = "Invalid step size "+gammalib::str(del)+
618  " encountered for"
619  " a="+gammalib::str(a)+
620  ", b="+gammalib::str(b)+
621  ", n="+gammalib::str(n)+
622  ", result="+gammalib::str(result)+
623  ". Step is too small to make sense.";
625  }
626 
627  // Sum up values
628  double x = a + 0.5*del;
629  double sum = 0.0;
630  for (int j = 0; j < it; ++j, x+=del) {
631 
632  // Evaluate integrand
633  double y = m_kernel->eval(x);
634  m_calls++;
635 
636  // Add integrand
637  sum += y;
638 
639  } // endfor: looped over steps
640 
641  // Set result
642  result = 0.5*(result + (b-a)*sum/tnm);
643  }
644 
645  } // endelse: trapeziodal rule was applied
646 
647  // Return result
648  return result;
649 }
650 
651 
652 /***********************************************************************//**
653  * @brief Adaptive Simpson's integration
654  *
655  * @param[in] a Left integration boundary.
656  * @param[in] b Right integration boundary.
657  *
658  * Integrates the function using an adaptive Simpson's rule. The initial
659  * interval [a,b] is split into two sub-intervals [a,c] and [c,b] for which
660  * the integral is computed using
661  *
662  * \f[
663  * \frac{b-a}{6} f(a) + 4f(c) + f(b)
664  * \f]
665  *
666  * where \f$c=(a+b)/2\f$ is the mid-point of interval [a,b]. Each
667  * sub-interval is then recursively divided into sub-interval and the process
668  * is repeated. Dividing of sub-intervals is stopped when the difference
669  * between subsequent intervals falls below the relative tolerance specified
670  * by eps(). The maximum recursion depth is set by the max_iter() method.
671  *
672  * I almost do not dare to confess, but the code has been taken from
673  * http://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method
674  * It's really pretty simple ...
675  ***************************************************************************/
676 double GIntegral::adaptive_simpson(const double& a, const double& b) const
677 {
678  // Initialise integration status information
679  m_isvalid = true;
680  m_calls = 0;
681  m_iter = m_max_iter;
682  m_has_abserr = false;
683  m_relerr = false;
684 
685  // Compute mid-point c
686  double c = 0.5*(a + b); //!< Mid-point of interval [a,b]
687  double h = b - a; //!< Length of interval [a,b]
688 
689  // Evaluate function at boundaries and mid-point c
690  double fa = m_kernel->eval(a);
691  double fb = m_kernel->eval(b);
692  double fc = m_kernel->eval(c);
693  m_calls += 3;
694 
695  // Compute integral using Simpson's rule
696  double S = (h/6.0) * (fa + 4.0*fc + fb);
697 
698  // Initialise absolute precision
699  double epsilon = (std::abs(S) > 0) ? m_eps * S : m_eps;
700 
701  // Call recursive auxiliary function
702  double value = adaptive_simpson_aux(a, b, epsilon, S, fa, fb, fc, m_max_iter);
703 
704  // Deduce the number of iterations from the iteration counter
706 
707  // If result is not valid, set and output status message
708  if (!m_isvalid) {
709  m_message = "Integration uncertainty exceeds relative tolerance "
710  "of "+gammalib::str(m_eps)+" and absolute tolerance of "
711  ""+gammalib::str(epsilon)+" after "+gammalib::str(m_iter)+
712  " iterations and "+gammalib::str(m_calls)+" function "
713  "calls. Result "+gammalib::str(value)+" inaccurate.";
714  if (!m_silent) {
715  std::string origin = "GIntegral::adaptive_simpson("+
716  gammalib::str(a)+", "+
717  gammalib::str(b)+")";
718  gammalib::warning(origin, m_message);
719  }
720  }
721 
722  // Return result
723  return value;
724 }
725 
726 
727 /***********************************************************************//**
728  * @brief Gauss-Kronrod integration
729  *
730  * @param[in] a Left integration boundary.
731  * @param[in] b Right integration boundary.
732  *
733  ***************************************************************************/
734 double GIntegral::gauss_kronrod(const double& a, const double& b) const
735 {
736  // Initialise integration status information
737  m_isvalid = true;
738  m_iter = 0;
739  m_calls = 0;
740  m_has_abserr = false;
741  m_relerr = false;
742 
743  // Initialise integration result
744  double result = 0.0;
745  double error = 0.0;
746 
747  // Allocate some arrays
748  double fv1[5];
749  double fv2[5];
750  double fv3[5];
751  double fv4[5];
752  double savfun[21];
753 
754  // Main code loop (so that we can exit using break)
755  do {
756 
757  // Tolerance check
758  if (m_eps < 1.12e-14) {
759  m_isvalid = false;
760  m_message = "Requested relative tolerance of "+gammalib::str(m_eps)+
761  " cannot be acheived. Please relax the integration "
762  "precision.";
763  if (!m_silent) {
764  std::string origin = "GIntegral::gauss_kronrod("+
765  gammalib::str(a)+", "+
766  gammalib::str(b)+")";
767  gammalib::warning(origin, m_message);
768  }
769  }
770 
771  // Compute function at mid-point
772  double h = 0.5 * (b - a);
773  double abs_h = std::abs(h);
774  double c = 0.5 * (b + a);
775  double f_c = m_kernel->eval(c);
776  m_calls++;
777 
778  // Compute the integral using the 10- and 21-point formulae
779  m_iter++;
780  double res10 = 0;
781  double res21 = gammalib::gkw21b[5] * f_c;
782  double resabs = gammalib::gkw21b[5] * std::abs(f_c);
783  for (int k = 0; k < 5; ++k) {
784  double x = h * gammalib::gkx1[k];
785  double fval1 = m_kernel->eval(c+x);
786  double fval2 = m_kernel->eval(c-x);
787  double fval = fval1 + fval2;
788  m_calls += 2;
789  res10 += gammalib::gkw10[k] * fval;
790  res21 += gammalib::gkw21a[k] * fval;
791  resabs += gammalib::gkw21a[k] * (std::abs(fval1) + std::abs(fval2));
792  savfun[k] = fval;
793  fv1[k] = fval1;
794  fv2[k] = fval2;
795  }
796  for (int k = 0; k < 5; ++k) {
797  double x = h * gammalib::gkx2[k];
798  double fval1 = m_kernel->eval(c+x);
799  double fval2 = m_kernel->eval(c-x);
800  double fval = fval1 + fval2;
801  m_calls += 2;
802  res21 += gammalib::gkw21b[k] * fval;
803  resabs += gammalib::gkw21b[k] * (std::abs(fval1) + std::abs(fval2));
804  savfun[k+5] = fval;
805  fv3[k] = fval1;
806  fv4[k] = fval2;
807  }
808  resabs *= abs_h;
809  double mean = 0.5 * res21;
810  double resasc = gammalib::gkw21b[5] * std::abs(f_c - mean);
811  for (int k = 0; k < 5; ++k) {
812  resasc += (gammalib::gkw21a[k] *
813  (std::abs(fv1[k] - mean) + std::abs(fv2[k] - mean)) +
814  gammalib::gkw21b[k] *
815  (std::abs(fv3[k] - mean) + std::abs(fv4[k] - mean)));
816  }
817  resasc *= abs_h ;
818  result = res21 * h;
819  error = rescale_error((res21 - res10) * h, resabs, resasc);
820 
821  // Test for convergence */
822  //if (err < epsabs || err < epsrel * fabs (result))
823  if (error < m_eps * std::abs(result)) {
824  m_has_abserr = true;
825  m_abserr = error;
826  if (std::abs(result) > 0) {
827  m_has_relerr = true;
828  m_relerr = error / std::abs(result);
829  }
830  break;
831  }
832 
833  // Compute the integral using the 43-point formula
834  m_iter++;
835  double res43 = gammalib::gkw43b[11] * f_c;
836  for (int k = 0; k < 10; ++k) {
837  res43 += savfun[k] * gammalib::gkw43a[k];
838  }
839  for (int k = 0; k < 11; ++k) {
840  double x = h * gammalib::gkx3[k];
841  double fval = (m_kernel->eval(c+x) +
842  m_kernel->eval(c-x));
843  m_calls += 2;
844  res43 += fval * gammalib::gkw43b[k];
845  savfun[k+10] = fval;
846  }
847  result = res43 * h;
848 
849  // Test for convergence */
850  error = rescale_error((res43 - res21) * h, resabs, resasc);
851  //if (err < epsabs || err < epsrel * fabs (result))
852  if (error < m_eps * std::abs(result)) {
853  m_has_abserr = true;
854  m_abserr = error;
855  if (std::abs(result) > 0) {
856  m_has_relerr = true;
857  m_relerr = error / std::abs(result);
858  }
859  break;
860  }
861 
862  // Compute the integral using the 87-point formula
863  m_iter++;
864  double res87 = gammalib::gkw87b[22] * f_c;
865  for (int k = 0; k < 21; ++k) {
866  res87 += savfun[k] * gammalib::gkw87a[k];
867  }
868  for (int k = 0; k < 22; ++k) {
869  double x = h * gammalib::gkx4[k];
870  res87 += gammalib::gkw87b[k] *
871  (m_kernel->eval(c+x) +
872  m_kernel->eval(c-x));
873  m_calls += 2;
874  }
875  result = res87 * h ;
876 
877  // Test for convergence */
878  error = rescale_error ((res87 - res43) * h, resabs, resasc);
879  //if (err < epsabs || err < epsrel * fabs (result))
880  if (error < m_eps * std::abs(result)) {
881  m_has_abserr = true;
882  m_abserr = error;
883  if (std::abs(result) > 0) {
884  m_has_relerr = true;
885  m_relerr = error / std::abs(result);
886  }
887  break;
888  }
889 
890  // Failed to converge
891  m_isvalid = false;
892  m_message = "Integration uncertainty "+gammalib::str(error)+" exceeds "
893  "absolute tolerance of "+
894  gammalib::str(m_eps * std::abs(result))+" after "+
895  gammalib::str(m_iter)+" iterations and "+
896  gammalib::str(m_calls)+" function calls. Result "+
897  gammalib::str(result)+" inaccurate.";
898  if (!m_silent) {
899  std::string origin = "GIntegral::gauss_kronrod("+
900  gammalib::str(a)+", "+
901  gammalib::str(b)+")";
902  gammalib::warning(origin, m_message);
903  }
904 
905  } while (false); // end of main loop
906 
907  // Return result
908  return result;
909 }
910 
911 
912 /***********************************************************************//**
913  * @brief Print integral information
914  *
915  * @param[in] chatter Chattiness (defaults to NORMAL).
916  * @return String containing integral information.
917  ***************************************************************************/
918 std::string GIntegral::print(const GChatter& chatter) const
919 {
920  // Initialise result string
921  std::string result;
922 
923  // Continue only if chatter is not silent
924  if (chatter != SILENT) {
925 
926  // Append header
927  result.append("=== GIntegral ===");
928 
929  // Append information
930  result.append("\n"+gammalib::parformat("Relative precision"));
931  result.append(gammalib::str(eps()));
932  if (m_has_abserr) {
933  result.append("\n"+gammalib::parformat("Absolute error"));
934  result.append(gammalib::str(m_abserr));
935  }
936  if (m_has_relerr) {
937  result.append("\n"+gammalib::parformat("Relative error"));
938  result.append(gammalib::str(m_relerr));
939  }
940  result.append("\n"+gammalib::parformat("Function calls"));
941  result.append(gammalib::str(calls()));
942  result.append("\n"+gammalib::parformat("Iterations"));
943  result.append(gammalib::str(iter()));
944  if (m_fix_iter > 0) {
945  result.append(" (fixed: ");
946  result.append(gammalib::str(fixed_iter()));
947  result.append(")");
948  }
949  else {
950  result.append(" (maximum: ");
951  result.append(gammalib::str(max_iter()));
952  result.append(")");
953  }
954 
955  // Append status information
956  result.append("\n"+gammalib::parformat("Status"));
957  if (is_valid()) {
958  result.append("Result accurate.");
959  }
960  else {
961  result.append(message());
962  }
963  if (silent()) {
964  result.append("\n"+gammalib::parformat("Warnings")+"suppressed");
965  }
966  else {
967  result.append("\n"+gammalib::parformat("Warnings"));
968  result.append("in standard output");
969  }
970 
971  } // endif: chatter was not silent
972 
973  // Return result
974  return result;
975 }
976 
977 
978 /*==========================================================================
979  = =
980  = Protected methods =
981  = =
982  ==========================================================================*/
983 
984 /***********************************************************************//**
985  * @brief Initialise class members
986  ***************************************************************************/
988 {
989  // Initialise members
990  m_kernel = NULL;
991  m_eps = 1.0e-6;
992  m_max_iter = 20;
993  m_fix_iter = 0;
994  m_message.clear();
995  m_silent = false;
996 
997  // Initialise results
998  m_iter = 0;
999  m_calls = 0;
1000  m_isvalid = true;
1001  m_has_abserr = false;
1002  m_has_relerr = false;
1003  m_abserr = 0.0;
1004  m_relerr = 0.0;
1005 
1006  // Return
1007  return;
1008 }
1009 
1010 
1011 /***********************************************************************//**
1012  * @brief Copy class members
1013  *
1014  * @param[in] integral Integral.
1015  ***************************************************************************/
1016 void GIntegral::copy_members(const GIntegral& integral)
1017 {
1018  // Copy attributes
1019  m_kernel = integral.m_kernel;
1020  m_eps = integral.m_eps;
1021  m_max_iter = integral.m_max_iter;
1022  m_fix_iter = integral.m_fix_iter;
1023  m_message = integral.m_message;
1024  m_silent = integral.m_silent;
1025 
1026  // Copy results
1027  m_iter = integral.m_iter;
1028  m_calls = integral.m_calls;
1029  m_isvalid = integral.m_isvalid;
1030  m_has_abserr = integral.m_has_abserr;
1031  m_has_relerr = integral.m_has_relerr;
1032  m_abserr = integral.m_abserr;
1033  m_relerr = integral.m_relerr;
1034 
1035  // Return
1036  return;
1037 }
1038 
1039 
1040 /***********************************************************************//**
1041  * @brief Delete class members
1042  ***************************************************************************/
1044 {
1045  // Return
1046  return;
1047 }
1048 
1049 
1050 /***********************************************************************//**
1051  * @brief Perform Polynomial interpolation
1052  *
1053  * @param[in] xa Pointer to array of X values.
1054  * @param[in] ya Pointer to array of Y values.
1055  * @param[in] n Number of elements in arrays.
1056  * @param[in] x X value at which interpolations should be performed.
1057  * @param[out] dy Error estimate for interpolated values.
1058  *
1059  * Given arrays xa[1,..,n] and ya[1,..,n], and given a value x, this
1060  * method returns a value y, and an error estimate dy. If P(x) is the
1061  * polynomial of degree n-1, then the returned value y=P(x).
1062  *
1063  * @todo Implement exceptions instead of screen dump.
1064  * @todo Use std::vector for xa and ya and start at 0
1065  ***************************************************************************/
1066 double GIntegral::polint(double* xa, double* ya, int n, double x, double* dy)
1067 {
1068  // Initialise result
1069  double y = 0.0;
1070 
1071  // Allocate temporary memory
1072  std::vector<double> c(n, 0.0);
1073  std::vector<double> d(n, 0.0);
1074 
1075  // Compute initial distance to first node
1076  double dif = std::abs(x-xa[1]);
1077 
1078  // Find index ns of the closest table entry
1079  int ns = 0;
1080  for (int i = 0; i < n; ++i) {
1081  double dift = std::abs(x-xa[i+1]);
1082  if (dift < dif) {
1083  ns = i;
1084  dif = dift;
1085  }
1086  c[i] = ya[i+1];
1087  d[i] = ya[i+1];
1088  }
1089 
1090  // Get initial approximation to y
1091  y = ya[ns+1];
1092  ns--;
1093 
1094  // Loop over each column of the tableau
1095  for (int m = 1; m < n; ++m) {
1096 
1097  // Update current c's and d's
1098  for (int i = 0; i < n-m; ++i) {
1099  double ho = xa[i+1] - x;
1100  double hp = xa[i+m+1] - x;
1101  double w = c[i+1] - d[i];
1102  double den = ho - hp;
1103  if (den == 0.0) {
1104  m_isvalid = false;
1105  m_message = "Invalid step size "+gammalib::str(den)+
1106  " encountered. Two values in xa array are"
1107  " identical.";
1109  }
1110  den = w/den;
1111  d[i] = hp*den;
1112  c[i] = ho*den;
1113  }
1114 
1115  // Compute y correction
1116  *dy = (2*(ns+1) < (n-m)) ? c[ns+1] : d[ns--];
1117 
1118  // Update y
1119  y += *dy;
1120 
1121  } // endfor: looped over columns of tableau
1122 
1123  // Return
1124  return y;
1125 }
1126 
1127 
1128 /***********************************************************************//**
1129  * @brief Auxiliary function for adaptive Simpson's method
1130  *
1131  * @param[in] a Left integration boundary.
1132  * @param[in] b Right integration boundary.
1133  * @param[in] eps Precision.
1134  * @param[in] S Integral of last computation.
1135  * @param[in] fa Function value at left integration boundary.
1136  * @param[in] fb Function value at right integration boundary.
1137  * @param[in] fc Function value at mid-point of interval [a,b]
1138  * @param[in] bottom Iteration counter (stop when 0)
1139  *
1140  * Implements a recursive auxiliary method for the adative_simpson()
1141  * integrator.
1142  ***************************************************************************/
1143 double GIntegral::adaptive_simpson_aux(const double& a, const double& b,
1144  const double& eps, const double& S,
1145  const double& fa, const double& fb,
1146  const double& fc,
1147  const int& bottom) const
1148 {
1149  // Store the iteration counter
1150  if (bottom < m_iter) {
1151  m_iter = bottom;
1152  }
1153 
1154  // Compute mid-point c bet
1155  double c = 0.5*(a + b); //!< Mid-point of interval [a,b]
1156  double h = b - a; //!< Length of interval [a,b]
1157  double d = 0.5*(a + c); //!< Mid-point of interval [a,c]
1158  double e = 0.5*(c + b); //!< Mid-point of interval [c,b]
1159 
1160  // Evaluate function at mid-points d and e
1161  double fd = m_kernel->eval(d);
1162  double fe = m_kernel->eval(e);
1163  m_calls += 2;
1164 
1165  // Compute integral using Simpson's rule for the left and right interval
1166  double h12 = h / 12.0;
1167  double Sleft = h12 * (fa + 4.0*fd + fc);
1168  double Sright = h12 * (fc + 4.0*fe + fb);
1169  double S2 = Sleft + Sright;
1170 
1171  // Allocate result
1172  double value;
1173 
1174  // If converged then compute the result ...
1175  if (std::abs(S2 - S) <= 15.0 * eps) {
1176 // if (std::abs(S2 - S) <= 15.0 * m_eps * std::abs(S2)) {
1177  value = S2 + (S2 - S)/15.0;
1178  }
1179 
1180  // ... else if the maximum recursion depth was reached then compute the
1181  // result and signal result invalidity
1182  else if (bottom <= 0) {
1183  value = S2 + (S2 - S)/15.0;
1184  m_isvalid = false;
1185  }
1186 
1187  // ... otherwise call this method recursively
1188  else {
1189  value = adaptive_simpson_aux(a, c, 0.5*eps, Sleft, fa, fc, fd, bottom-1) +
1190  adaptive_simpson_aux(c, b, 0.5*eps, Sright, fc, fb, fe, bottom-1);
1191  }
1192 
1193  // Return result
1194  return value;
1195 }
1196 
1197 
1198 //#define GSL_DBL_EPSILON 2.2204460492503131e-16
1199 //#define GSL_DBL_MIN 2.2250738585072014e-308
1200 /***********************************************************************//**
1201  * @brief Rescale errors for Gauss-Kronrod integration
1202  *
1203  * @param[in] err Error estimate.
1204  * @param[in] result_abs ???.
1205  * @param[in] result_asc ???.
1206  * @return Rescaled error estimate.
1207  ***************************************************************************/
1208 double GIntegral::rescale_error(double err, const double& result_abs, const double& result_asc) const
1209 {
1210  // Take absolute value of error
1211  err = std::abs(err);
1212 
1213  // ...
1214  if (result_asc != 0.0 && err != 0.0) {
1215  double scale = std::pow((200.0 * err / result_asc), 1.5);
1216  if (scale < 1.0) {
1217  err = result_asc * scale ;
1218  }
1219  else {
1220  err = result_asc ;
1221  }
1222  }
1223  if (result_abs > 2.2250738585072014e-308 / (50.0 * 2.2204460492503131e-16)) {
1224  double min_err = 50.0 * 2.2204460492503131e-16 * result_abs;
1225  if (min_err > err) {
1226  err = min_err ;
1227  }
1228  }
1229 
1230  // Return error
1231  return err;
1232 }
const double gkx2[5]
Definition: GIntegral.cpp:68
const double & eps(void) const
Get relative precision.
Definition: GIntegral.hpp:219
const double gkx1[5]
Definition: GIntegral.cpp:50
const int & fixed_iter(void) const
Return fixed number of iterations.
Definition: GIntegral.hpp:194
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
int m_iter
Number of iterations used.
Definition: GIntegral.hpp:109
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1263
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
const bool & silent(void) const
Get silence flag.
Definition: GIntegral.hpp:256
void init_members(void)
Initialise class members.
Definition: GIntegral.cpp:987
const double gkw43a[10]
Definition: GIntegral.cpp:111
#define G_ROMBERG
Definition: GIntegral.cpp:36
double polint(double *xa, double *ya, int n, double x, double *dy)
Perform Polynomial interpolation.
Definition: GIntegral.cpp:1066
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:901
double rescale_error(double err, const double &result_abs, const double &result_asc) const
Rescale errors for Gauss-Kronrod integration.
Definition: GIntegral.cpp:1208
int m_max_iter
Maximum number of iterations.
Definition: GIntegral.hpp:104
const std::string & message(void) const
Return integration status message.
Definition: GIntegral.hpp:307
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
double gauss_kronrod(const double &a, const double &b) const
Gauss-Kronrod integration.
Definition: GIntegral.cpp:734
int m_fix_iter
Fixed number of iterations.
Definition: GIntegral.hpp:105
GIntegral & operator=(const GIntegral &integral)
Assignment operator.
Definition: GIntegral.cpp:304
void clear(void)
Clear integral.
Definition: GIntegral.cpp:334
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:178
#define G_TRAPZD
Definition: GIntegral.cpp:37
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:161
const int & iter(void) const
Return number of iterations.
Definition: GIntegral.hpp:139
bool m_isvalid
Integration result valid (true=yes)
Definition: GIntegral.hpp:111
GIntegral(void)
Void constructor.
Definition: GIntegral.cpp:230
const double gkw87a[21]
Definition: GIntegral.cpp:167
const double gkw43b[12]
Definition: GIntegral.cpp:125
double trapzd(const double &a, const double &b, const int &n=1, double result=0.0)
Perform Trapezoidal integration.
Definition: GIntegral.cpp:564
bool m_has_abserr
Has absolute integration error.
Definition: GIntegral.hpp:112
std::string m_message
Status message (if result is invalid)
Definition: GIntegral.hpp:116
const double gkw87b[23]
Definition: GIntegral.cpp:192
GChatter
Definition: GTypemaps.hpp:33
bool m_has_relerr
Has relative integration error.
Definition: GIntegral.hpp:113
const int & calls(void) const
Get number of function calls.
Definition: GIntegral.hpp:231
virtual double eval(const double &x)=0
double adaptive_simpson_aux(const double &a, const double &b, const double &eps, const double &S, const double &fa, const double &fb, const double &fc, const int &bottom) const
Auxiliary function for adaptive Simpson&#39;s method.
Definition: GIntegral.cpp:1143
const GFunction * kernel(void) const
Get kernel.
Definition: GIntegral.hpp:283
double m_relerr
Absolute integration error.
Definition: GIntegral.hpp:115
const double gkw21a[5]
Definition: GIntegral.cpp:77
void copy_members(const GIntegral &integral)
Copy class members.
Definition: GIntegral.cpp:1016
double m_eps
Requested relative integration precision.
Definition: GIntegral.hpp:103
bool m_silent
Suppress integration warnings in console.
Definition: GIntegral.hpp:106
std::string print(const GChatter &chatter=NORMAL) const
Print integral information.
Definition: GIntegral.cpp:918
virtual ~GIntegral(void)
Destructor.
Definition: GIntegral.cpp:282
const double gkw10[5]
Definition: GIntegral.cpp:59
Single parameter function abstract base class.
Definition: GFunction.hpp:44
const int & max_iter(void) const
Return maximum number of iterations.
Definition: GIntegral.hpp:164
#define G_POLINT
Definition: GIntegral.cpp:38
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
Exception handler interface definition.
GFunction * m_kernel
Pointer to function kernel.
Definition: GIntegral.hpp:102
int m_calls
Number of function calls used.
Definition: GIntegral.hpp:110
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1020
const double gkw21b[6]
Definition: GIntegral.cpp:86
Integration class interface definition.
const double gkx4[22]
Definition: GIntegral.cpp:141
const double gkx3[11]
Definition: GIntegral.cpp:96
const bool & is_valid(void) const
Signal if integration result is valid.
Definition: GIntegral.hpp:295
double m_abserr
Absolute integration error.
Definition: GIntegral.hpp:114
void free_members(void)
Delete class members.
Definition: GIntegral.cpp:1043
GIntegral * clone(void) const
Clone integral.
Definition: GIntegral.cpp:352
double adaptive_simpson(const double &a, const double &b) const
Adaptive Simpson&#39;s integration.
Definition: GIntegral.cpp:676
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:411