Samuel Flis / gammalib

Commit 01e2aa98f04047722eda22196e732b94499fe4a6

Authored by Jürgen Knödlseder
1 parent 07a039a9

Implement Gauss-Kronrod quadrature in GIntegral class.

Inline Side-by-side
Showing 5 changed files with 503 additions and 20 deletions
include/GIntegral.hpp
  View file @01e2aa9
... ... @@ -39,7 +39,7 @@
39 39 * @brief GIntegral class interface definition.
40 40 *
41 41 * This class allows to perform integration using various methods. The
42   - * integrand is implemented by a derived class of GIntegrand.
  42 + * integrand is implemented by a derived class of GFunction.
43 43 ***************************************************************************/
44 44 class GIntegral : public GBase {
45 45  
... ... @@ -73,6 +73,7 @@ public:
73 73 double trapzd(const double& a, const double& b, const int& n = 1,
74 74 double result = 0.0);
75 75 double adaptive_simpson(const double& a, const double& b) const;
  76 + double gauss_kronrod(const double& a, const double& b) const;
76 77 std::string print(const GChatter& chatter = NORMAL) const;
77 78  
78 79 protected:
... ... @@ -86,6 +87,9 @@ protected:
86 87 const double& fa, const double& fb,
87 88 const double& fc,
88 89 const int& bottom) const;
  90 + double rescale_error(double err,
  91 + const double& result_abs,
  92 + const double& result_asc) const;
89 93  
90 94 // Protected data area
91 95 GFunction* m_kernel; //!< Pointer to function kernel
... ... @@ -94,10 +98,14 @@ protected:
94 98 bool m_silent; //!< Suppress integration warnings in console
95 99  
96 100 // Integrator results
97   - mutable int m_iter; //!< Number of iterations used
98   - mutable int m_calls; //!< Number of function calls used
99   - mutable bool m_isvalid; //!< Integration result valid (true=yes)
100   - mutable std::string m_message; //!< Status message (if result is invalid)
  101 + mutable int m_iter; //!< Number of iterations used
  102 + mutable int m_calls; //!< Number of function calls used
  103 + mutable bool m_isvalid; //!< Integration result valid (true=yes)
  104 + mutable bool m_has_abserr; //!< Has absolute integration error
  105 + mutable bool m_has_relerr; //!< Has relative integration error
  106 + mutable double m_abserr; //!< Absolute integration error
  107 + mutable double m_relerr; //!< Absolute integration error
  108 + mutable std::string m_message; //!< Status message (if result is invalid)
101 109  
102 110 };
103 111  
... ...
pyext/GIntegral.i
  View file @01e2aa9
... ... @@ -36,7 +36,7 @@
36 36 * @brief Integration class Python interface definition.
37 37 *
38 38 * This class allows to perform integration using various methods. The
39   - * integrand is implemented by a derived class of GIntegrand.
  39 + * integrand is implemented by a derived class of GFunction.
40 40 ***************************************************************************/
41 41 class GIntegral : public GBase {
42 42 public:
... ... @@ -54,16 +54,19 @@ public:
54 54 void eps(const double& eps);
55 55 void silent(const bool& silent);
56 56 const int& iter(void) const;
  57 + const int& calls(void) const;
57 58 const int& max_iter(void) const;
58 59 const double& eps(void) const;
59 60 const bool& silent(void) const;
60   - const bool& isvalid(void) const;
  61 + const bool& is_valid(void) const;
61 62 const std::string& message(void) const;
62 63 void kernel(GFunction* kernel);
63 64 const GFunction* kernel(void) const;
64 65 double romb(const double& a, const double& b, const int& k = 5);
65 66 double trapzd(const double& a, const double& b, const int& n = 1,
66 67 double result = 0.0);
  68 + double adaptive_simpson(const double& a, const double& b) const;
  69 + double gauss_kronrod(const double& a, const double& b) const;
67 70 };
68 71  
69 72  
... ...
src/numerics/GIntegral.cpp
  View file @01e2aa9
... ... @@ -41,6 +41,180 @@
41 41  
42 42 /* __ Debug definitions __________________________________________________ */
43 43  
  44 +/* __ Constants __________________________________________________________ */
  45 +namespace gammalib {
  46 +
  47 + // Gauss-Kronrod abscissae, common to the 10-, 21-, 43- and 87-point rule
  48 + const double gkx1[5] = {
  49 + 0.973906528517171720077964012084452,
  50 + 0.865063366688984510732096688423493,
  51 + 0.679409568299024406234327365114874,
  52 + 0.433395394129247190799265943165784,
  53 + 0.148874338981631210884826001129720
  54 + };
  55 +
  56 + // Gauss-Kronrod weights of the 10-point rule
  57 + const double gkw10[5] = {
  58 + 0.066671344308688137593568809893332,
  59 + 0.149451349150580593145776339657697,
  60 + 0.219086362515982043995534934228163,
  61 + 0.269266719309996355091226921569469,
  62 + 0.295524224714752870173892994651338
  63 + };
  64 +
  65 + // Gauss-Kronrod abscissae, common to the 21-, 43- and 87-point rule
  66 + const double gkx2[5] = {
  67 + 0.995657163025808080735527280689003,
  68 + 0.930157491355708226001207180059508,
  69 + 0.780817726586416897063717578345042,
  70 + 0.562757134668604683339000099272694,
  71 + 0.294392862701460198131126603103866
  72 + };
  73 +
  74 + // Gauss-Kronrod weights of the 21-point rule for abscissae gkx1
  75 + const double gkw21a[5] = {
  76 + 0.032558162307964727478818972459390,
  77 + 0.075039674810919952767043140916190,
  78 + 0.109387158802297641899210590325805,
  79 + 0.134709217311473325928054001771707,
  80 + 0.147739104901338491374841515972068
  81 + };
  82 +
  83 + // Gauss-Kronrod weights of the 21-point rule for abscissae gkx2
  84 + const double gkw21b[6] = {
  85 + 0.011694638867371874278064396062192,
  86 + 0.054755896574351996031381300244580,
  87 + 0.093125454583697605535065465083366,
  88 + 0.123491976262065851077958109831074,
  89 + 0.142775938577060080797094273138717,
  90 + 0.149445554002916905664936468389821
  91 + };
  92 +
  93 + // Gauss-Kronrod abscissae, common to the 43- and 87-point rule
  94 + const double gkx3[11] = {
  95 + 0.999333360901932081394099323919911,
  96 + 0.987433402908088869795961478381209,
  97 + 0.954807934814266299257919200290473,
  98 + 0.900148695748328293625099494069092,
  99 + 0.825198314983114150847066732588520,
  100 + 0.732148388989304982612354848755461,
  101 + 0.622847970537725238641159120344323,
  102 + 0.499479574071056499952214885499755,
  103 + 0.364901661346580768043989548502644,
  104 + 0.222254919776601296498260928066212,
  105 + 0.074650617461383322043914435796506
  106 + };
  107 +
  108 + // Gauss-Kronrod weights of the 43-point rule for abscissae gkx1, gkx3
  109 + const double gkw43a[10] = {
  110 + 0.016296734289666564924281974617663,
  111 + 0.037522876120869501461613795898115,
  112 + 0.054694902058255442147212685465005,
  113 + 0.067355414609478086075553166302174,
  114 + 0.073870199632393953432140695251367,
  115 + 0.005768556059769796184184327908655,
  116 + 0.027371890593248842081276069289151,
  117 + 0.046560826910428830743339154433824,
  118 + 0.061744995201442564496240336030883,
  119 + 0.071387267268693397768559114425516
  120 + };
  121 +
  122 + // Gauss-Kronrod weights of the 43-point formula for abscissae gkx3
  123 + const double gkw43b[12] = {
  124 + 0.001844477640212414100389106552965,
  125 + 0.010798689585891651740465406741293,
  126 + 0.021895363867795428102523123075149,
  127 + 0.032597463975345689443882222526137,
  128 + 0.042163137935191811847627924327955,
  129 + 0.050741939600184577780189020092084,
  130 + 0.058379395542619248375475369330206,
  131 + 0.064746404951445885544689259517511,
  132 + 0.069566197912356484528633315038405,
  133 + 0.072824441471833208150939535192842,
  134 + 0.074507751014175118273571813842889,
  135 + 0.074722147517403005594425168280423
  136 + };
  137 +
  138 + // Gauss-Kronrod abscissae, of the 87-point rule
  139 + const double gkx4[22] = {
  140 + 0.999902977262729234490529830591582,
  141 + 0.997989895986678745427496322365960,
  142 + 0.992175497860687222808523352251425,
  143 + 0.981358163572712773571916941623894,
  144 + 0.965057623858384619128284110607926,
  145 + 0.943167613133670596816416634507426,
  146 + 0.915806414685507209591826430720050,
  147 + 0.883221657771316501372117548744163,
  148 + 0.845710748462415666605902011504855,
  149 + 0.803557658035230982788739474980964,
  150 + 0.757005730685495558328942793432020,
  151 + 0.706273209787321819824094274740840,
  152 + 0.651589466501177922534422205016736,
  153 + 0.593223374057961088875273770349144,
  154 + 0.531493605970831932285268948562671,
  155 + 0.466763623042022844871966781659270,
  156 + 0.399424847859218804732101665817923,
  157 + 0.329874877106188288265053371824597,
  158 + 0.258503559202161551802280975429025,
  159 + 0.185695396568346652015917141167606,
  160 + 0.111842213179907468172398359241362,
  161 + 0.037352123394619870814998165437704
  162 + };
  163 +
  164 + // Gauss-Kronrod weights of the 87-point rule for abscissae gkx1, gkx2, gkx3
  165 + const double gkw87a[21] = {
  166 + 0.008148377384149172900002878448190,
  167 + 0.018761438201562822243935059003794,
  168 + 0.027347451050052286161582829741283,
  169 + 0.033677707311637930046581056957588,
  170 + 0.036935099820427907614589586742499,
  171 + 0.002884872430211530501334156248695,
  172 + 0.013685946022712701888950035273128,
  173 + 0.023280413502888311123409291030404,
  174 + 0.030872497611713358675466394126442,
  175 + 0.035693633639418770719351355457044,
  176 + 0.000915283345202241360843392549948,
  177 + 0.005399280219300471367738743391053,
  178 + 0.010947679601118931134327826856808,
  179 + 0.016298731696787335262665703223280,
  180 + 0.021081568889203835112433060188190,
  181 + 0.025370969769253827243467999831710,
  182 + 0.029189697756475752501446154084920,
  183 + 0.032373202467202789685788194889595,
  184 + 0.034783098950365142750781997949596,
  185 + 0.036412220731351787562801163687577,
  186 + 0.037253875503047708539592001191226
  187 + };
  188 +
  189 + // Gauss-Kronrod weights of the 87-point formula for abscissae gkx4
  190 + const double gkw87b[23] = {
  191 + 0.000274145563762072350016527092881,
  192 + 0.001807124155057942948341311753254,
  193 + 0.004096869282759164864458070683480,
  194 + 0.006758290051847378699816577897424,
  195 + 0.009549957672201646536053581325377,
  196 + 0.012329447652244853694626639963780,
  197 + 0.015010447346388952376697286041943,
  198 + 0.017548967986243191099665352925900,
  199 + 0.019938037786440888202278192730714,
  200 + 0.022194935961012286796332102959499,
  201 + 0.024339147126000805470360647041454,
  202 + 0.026374505414839207241503786552615,
  203 + 0.028286910788771200659968002987960,
  204 + 0.030052581128092695322521110347341,
  205 + 0.031646751371439929404586051078883,
  206 + 0.033050413419978503290785944862689,
  207 + 0.034255099704226061787082821046821,
  208 + 0.035262412660156681033782717998428,
  209 + 0.036076989622888701185500318003895,
  210 + 0.036698604498456094498018047441094,
  211 + 0.037120549269832576114119958413599,
  212 + 0.037334228751935040321235449094698,
  213 + 0.037361073762679023410321241766599
  214 + };
  215 +
  216 +} // end gammalib namespace
  217 +
44 218  
45 219 /*==========================================================================
46 220 = =
... ... @@ -431,8 +605,11 @@ double GIntegral::adaptive_simpson(const double&amp; a, const double&amp; b) const
431 605 // Compute integral using Simpson's rule
432 606 double S = (h/6.0) * (fa + 4.0*fc + fb);
433 607  
  608 + // Initialise absolute precision
  609 + double epsilon = (std::abs(S) > 0) ? m_eps * S : m_eps;
  610 +
434 611 // Call recursive auxiliary function
435   - double value = adaptive_simpson_aux(a, b, m_eps, S, fa, fb, fc, m_max_iter);
  612 + double value = adaptive_simpson_aux(a, b, epsilon, S, fa, fb, fc, m_max_iter);
436 613  
437 614 // Deduce the number of iterations from the iteration counter
438 615 m_iter = m_max_iter - m_iter;
... ... @@ -440,8 +617,10 @@ double GIntegral::adaptive_simpson(const double&amp; a, const double&amp; b) const
440 617 // If result is not valid, set and output status message
441 618 if (!m_isvalid) {
442 619 m_message = "Integration uncertainty exceeds relative tolerance "
443   - "of "+gammalib::str(m_eps)+" after "+gammalib::str(m_iter)+
444   - " iterations. Result "+gammalib::str(value)+" inaccurate.";
  620 + "of "+gammalib::str(m_eps)+" and absolute tolerance of "
  621 + ""+gammalib::str(epsilon)+" after "+gammalib::str(m_iter)+
  622 + " iterations and "+gammalib::str(m_calls)+" function "
  623 + "calls. Result "+gammalib::str(value)+" inaccurate.";
445 624 if (!m_silent) {
446 625 std::string origin = "GIntegral::adaptive_simpson("+
447 626 gammalib::str(a)+", "+
... ... @@ -456,6 +635,179 @@ double GIntegral::adaptive_simpson(const double&amp; a, const double&amp; b) const
456 635  
457 636  
458 637 /***********************************************************************//**
  638 + * @brief Gauss-Kronrod integration
  639 + *
  640 + * @param[in] a Left integration boundary.
  641 + * @param[in] b Right integration boundary.
  642 + *
  643 + ***************************************************************************/
  644 +double GIntegral::gauss_kronrod(const double& a, const double& b) const
  645 +{
  646 + // Initialise integration status information
  647 + m_isvalid = true;
  648 + m_iter = 0;
  649 + m_calls = 0;
  650 +
  651 + // Initialise integration result
  652 + double result = 0.0;
  653 + double error = 0.0;
  654 + //double resabs = 0.0; // Approximation to the integral of abs(f)
  655 + //double resasc = 0.0; // Approximation to the integral of abs(f-i/(b-a))
  656 +
  657 + // Allocate some arrays
  658 + double fv1[5];
  659 + double fv2[5];
  660 + double fv3[5];
  661 + double fv4[5];
  662 + double savfun[21];
  663 +
  664 + // Main code loop (so that we can exit using break)
  665 + do {
  666 +
  667 + // Tolerance check
  668 + if (m_eps < 1.12e-14) {
  669 + m_isvalid = false;
  670 + m_message = "Requested relative tolerance of "+gammalib::str(m_eps)+
  671 + " cannot be acheived. Please relax the integration "
  672 + "precision.";
  673 + if (!m_silent) {
  674 + std::string origin = "GIntegral::gauss_kronrod("+
  675 + gammalib::str(a)+", "+
  676 + gammalib::str(b)+")";
  677 + gammalib::warning(origin, m_message);
  678 + }
  679 + }
  680 +
  681 + // Compute function at mid-point
  682 + double h = 0.5 * (b - a);
  683 + double abs_h = std::abs(h);
  684 + double c = 0.5 * (b + a);
  685 + double f_c = m_kernel->eval(c);
  686 + m_calls++;
  687 +
  688 + // Compute the integral using the 10- and 21-point formulae
  689 + m_iter++;
  690 + double res10 = 0;
  691 + double res21 = gammalib::gkw21b[5] * f_c;
  692 + double resabs = gammalib::gkw21b[5] * std::abs(f_c);
  693 + for (int k = 0; k < 5; ++k) {
  694 + double x = h * gammalib::gkx1[k];
  695 + double fval1 = m_kernel->eval(c+x);
  696 + double fval2 = m_kernel->eval(c-x);
  697 + double fval = fval1 + fval2;
  698 + m_calls += 2;
  699 + res10 += gammalib::gkw10[k] * fval;
  700 + res21 += gammalib::gkw21a[k] * fval;
  701 + resabs += gammalib::gkw21a[k] * (std::abs(fval1) + std::abs(fval2));
  702 + savfun[k] = fval;
  703 + fv1[k] = fval1;
  704 + fv2[k] = fval2;
  705 + }
  706 + for (int k = 0; k < 5; ++k) {
  707 + double x = h * gammalib::gkx2[k];
  708 + double fval1 = m_kernel->eval(c+x);
  709 + double fval2 = m_kernel->eval(c-x);
  710 + double fval = fval1 + fval2;
  711 + m_calls += 2;
  712 + res21 += gammalib::gkw21b[k] * fval;
  713 + resabs += gammalib::gkw21b[k] * (std::abs(fval1) + std::abs(fval2));
  714 + savfun[k+5] = fval;
  715 + fv3[k] = fval1;
  716 + fv4[k] = fval2;
  717 + }
  718 + resabs *= abs_h;
  719 + double mean = 0.5 * res21;
  720 + double resasc = gammalib::gkw21b[5] * std::abs(f_c - mean);
  721 + for (int k = 0; k < 5; ++k) {
  722 + resasc += (gammalib::gkw21a[k] *
  723 + (std::abs(fv1[k] - mean) + std::abs(fv2[k] - mean)) +
  724 + gammalib::gkw21b[k] *
  725 + (std::abs(fv3[k] - mean) + std::abs(fv4[k] - mean)));
  726 + }
  727 + resasc *= abs_h ;
  728 + result = res21 * h;
  729 + error = rescale_error((res21 - res10) * h, resabs, resasc);
  730 +
  731 + // Test for convergence */
  732 + //if (err < epsabs || err < epsrel * fabs (result))
  733 + if (error < m_eps * std::abs(result)) {
  734 + m_has_abserr = true;
  735 + m_abserr = error;
  736 + break;
  737 + }
  738 +
  739 + // Compute the integral using the 43-point formula
  740 + m_iter++;
  741 + double res43 = gammalib::gkw43b[11] * f_c;
  742 + for (int k = 0; k < 10; ++k) {
  743 + res43 += savfun[k] * gammalib::gkw43a[k];
  744 + }
  745 + for (int k = 0; k < 11; ++k) {
  746 + double x = h * gammalib::gkx3[k];
  747 + double fval = (m_kernel->eval(c+x) +
  748 + m_kernel->eval(c-x));
  749 + m_calls += 2;
  750 + res43 += fval * gammalib::gkw43b[k];
  751 + savfun[k+10] = fval;
  752 + }
  753 + result = res43 * h;
  754 +
  755 + // Test for convergence */
  756 + error = rescale_error((res43 - res21) * h, resabs, resasc);
  757 + //if (err < epsabs || err < epsrel * fabs (result))
  758 + if (error < m_eps * std::abs(result)) {
  759 + m_has_abserr = true;
  760 + m_abserr = error;
  761 + break;
  762 + }
  763 +
  764 + // Compute the integral using the 87-point formula
  765 + m_iter++;
  766 + double res87 = gammalib::gkw87b[22] * f_c;
  767 + for (int k = 0; k < 21; ++k) {
  768 + res87 += savfun[k] * gammalib::gkw87a[k];
  769 + }
  770 + for (int k = 0; k < 22; ++k) {
  771 + double x = h * gammalib::gkx4[k];
  772 + res87 += gammalib::gkw87b[k] *
  773 + (m_kernel->eval(c+x) +
  774 + m_kernel->eval(c-x));
  775 + m_calls += 2;
  776 + }
  777 + result = res87 * h ;
  778 +
  779 + // Test for convergence */
  780 + error = rescale_error ((res87 - res43) * h, resabs, resasc);
  781 + //if (err < epsabs || err < epsrel * fabs (result))
  782 + if (error < m_eps * std::abs(result)) {
  783 + m_has_abserr = true;
  784 + m_abserr = error;
  785 + break;
  786 + }
  787 +
  788 + // Failed to converge
  789 + m_isvalid = false;
  790 + m_message = "Integration uncertainty "+gammalib::str(error)+" exceeds "
  791 + "absolute tolerance of "+
  792 + gammalib::str(m_eps * std::abs(result))+" after "+
  793 + gammalib::str(m_iter)+" iterations and "+
  794 + gammalib::str(m_calls)+" function calls. Result "+
  795 + gammalib::str(result)+" inaccurate.";
  796 + if (!m_silent) {
  797 + std::string origin = "GIntegral::gauss_kronrod("+
  798 + gammalib::str(a)+", "+
  799 + gammalib::str(b)+")";
  800 + gammalib::warning(origin, m_message);
  801 + }
  802 +
  803 + } while (false); // end of main loop
  804 +
  805 + // Return result
  806 + return result;
  807 +}
  808 +
  809 +
  810 +/***********************************************************************//**
459 811 * @brief Print integral information
460 812 *
461 813 * @param[in] chatter Chattiness (defaults to NORMAL).
... ... @@ -475,6 +827,14 @@ std::string GIntegral::print(const GChatter&amp; chatter) const
475 827 // Append information
476 828 result.append("\n"+gammalib::parformat("Relative precision"));
477 829 result.append(gammalib::str(eps()));
  830 + if (m_has_abserr) {
  831 + result.append("\n"+gammalib::parformat("Absolute error"));
  832 + result.append(gammalib::str(m_abserr));
  833 + }
  834 + if (m_has_relerr) {
  835 + result.append("\n"+gammalib::parformat("Relative error"));
  836 + result.append(gammalib::str(m_relerr));
  837 + }
478 838 result.append("\n"+gammalib::parformat("Function calls"));
479 839 result.append(gammalib::str(calls()));
480 840 result.append("\n"+gammalib::parformat("Iterations"));
... ... @@ -521,11 +881,17 @@ void GIntegral::init_members(void)
521 881 m_kernel = NULL;
522 882 m_eps = 1.0e-6;
523 883 m_max_iter = 20;
524   - m_iter = 0;
525   - m_isvalid = true;
526 884 m_message.clear();
527 885 m_silent = false;
528   - m_calls = 0;
  886 +
  887 + // Initialise results
  888 + m_iter = 0;
  889 + m_calls = 0;
  890 + m_isvalid = true;
  891 + m_has_abserr = false;
  892 + m_has_relerr = false;
  893 + m_abserr = 0.0;
  894 + m_relerr = 0.0;
529 895  
530 896 // Return
531 897 return;
... ... @@ -543,11 +909,17 @@ void GIntegral::copy_members(const GIntegral&amp; integral)
543 909 m_kernel = integral.m_kernel;
544 910 m_eps = integral.m_eps;
545 911 m_max_iter = integral.m_max_iter;
546   - m_iter = integral.m_iter;
547   - m_isvalid = integral.m_isvalid;
548 912 m_message = integral.m_message;
549 913 m_silent = integral.m_silent;
550   - m_calls = integral.m_calls;
  914 +
  915 + // Copy results
  916 + m_iter = integral.m_iter;
  917 + m_calls = integral.m_calls;
  918 + m_isvalid = integral.m_isvalid;
  919 + m_has_abserr = integral.m_has_abserr;
  920 + m_has_relerr = integral.m_has_relerr;
  921 + m_abserr = integral.m_abserr;
  922 + m_relerr = integral.m_relerr;
551 923  
552 924 // Return
553 925 return;
... ... @@ -643,7 +1015,7 @@ double GIntegral::polint(double* xa, double* ya, int n, double x, double* dy)
643 1015  
644 1016  
645 1017 /***********************************************************************//**
646   - * @brief Auxiliary function for adaptive Simpson's method.
  1018 + * @brief Auxiliary function for adaptive Simpson's method
647 1019 *
648 1020 * @param[in] a Left integration boundary.
649 1021 * @param[in] b Right integration boundary.
... ... @@ -663,7 +1035,9 @@ double GIntegral::adaptive_simpson_aux(const double&amp; a, const double&amp; b,
663 1035 const int& bottom) const
664 1036 {
665 1037 // Store the iteration counter
666   - m_iter = bottom;
  1038 + if (bottom < m_iter) {
  1039 + m_iter = bottom;
  1040 + }
667 1041  
668 1042 // Compute mid-point c bet
669 1043 double c = 0.5*(a + b); //!< Mid-point of interval [a,b]
... ... @@ -686,8 +1060,8 @@ double GIntegral::adaptive_simpson_aux(const double&amp; a, const double&amp; b,
686 1060 double value;
687 1061  
688 1062 // If converged then compute the result ...
689   -// if (std::abs(S2 - S) <= 15.0 * eps * std::abs(S2)) {
690   - if (std::abs(S2 - S) <= 15.0 * m_eps * std::abs(S2)) {
  1063 + if (std::abs(S2 - S) <= 15.0 * eps) {
  1064 +// if (std::abs(S2 - S) <= 15.0 * m_eps * std::abs(S2)) {
691 1065 value = S2 + (S2 - S)/15.0;
692 1066 }
693 1067  
... ... @@ -707,3 +1081,40 @@ double GIntegral::adaptive_simpson_aux(const double&amp; a, const double&amp; b,
707 1081 // Return result
708 1082 return value;
709 1083 }
  1084 +
  1085 +
  1086 +//#define GSL_DBL_EPSILON 2.2204460492503131e-16
  1087 +//#define GSL_DBL_MIN 2.2250738585072014e-308
  1088 +/***********************************************************************//**
  1089 + * @brief Rescale errors for Gauss-Kronrod integration
  1090 + *
  1091 + * @param[in] err Error estimate.
  1092 + * @param[in] result_abs ???.
  1093 + * @param[in] result_asc ???.
  1094 + * @return Rescaled error estimate.
  1095 + ***************************************************************************/
  1096 +double GIntegral::rescale_error(double err, const double& result_abs, const double& result_asc) const
  1097 +{
  1098 + // Take absolute value of error
  1099 + err = std::abs(err);
  1100 +
  1101 + // ...
  1102 + if (result_asc != 0.0 && err != 0.0) {
  1103 + double scale = std::pow((200.0 * err / result_asc), 1.5);
  1104 + if (scale < 1.0) {
  1105 + err = result_asc * scale ;
  1106 + }
  1107 + else {
  1108 + err = result_asc ;
  1109 + }
  1110 + }
  1111 + if (result_abs > 2.2250738585072014e-308 / (50.0 * 2.2204460492503131e-16)) {
  1112 + double min_err = 50.0 * 2.2204460492503131e-16 * result_abs;
  1113 + if (min_err > err) {
  1114 + err = min_err ;
  1115 + }
  1116 + }
  1117 +
  1118 + // Return error
  1119 + return err;
  1120 +}
... ...
test/test_GNumerics.cpp
  View file @01e2aa9
... ... @@ -53,6 +53,7 @@ void TestGNumerics::set(void)
53 53 append(static_cast<pfunction>(&TestGNumerics::test_integral),"Test GIntegral");
54 54 append(static_cast<pfunction>(&TestGNumerics::test_romberg_integration),"Test Romberg integration");
55 55 append(static_cast<pfunction>(&TestGNumerics::test_adaptive_simpson_integration),"Test adaptive Simpson integration");
  56 + append(static_cast<pfunction>(&TestGNumerics::test_gauss_kronrod_integration),"Test Gauss-Kronrod integration");
56 57  
57 58 // Return
58 59 return;
... ... @@ -95,9 +96,11 @@ void TestGNumerics::test_romberg_integration(void)
95 96  
96 97 // Integrate over the entire Gaussian
97 98 double result = integral.romb(-10.0*m_sigma, 10.0*m_sigma);
  99 +/*
98 100 std::cout << result << std::endl;
99 101 std::cout << result-1.0 << std::endl;
100 102 std::cout << integral << std::endl;
  103 +*/
101 104 test_value(result,1.0,1.0e-6,"","Gaussian integral is not 1.0 (integral="+gammalib::str(result)+")");
102 105  
103 106 // Test [-1sigma, 1sigma]
... ... @@ -124,20 +127,77 @@ void TestGNumerics::test_adaptive_simpson_integration(void)
124 127  
125 128 // Integrate over the entire Gaussian
126 129 double result = integral.adaptive_simpson(-10.0*m_sigma, 10.0*m_sigma);
  130 +/*
  131 +std::cout << result << std::endl;
  132 +std::cout << result-1.0 << std::endl;
  133 +std::cout << integral << std::endl;
  134 +*/
127 135 test_value(result,1.0,1.0e-6,"","Gaussian integral is not 1.0 (integral="+gammalib::str(result)+")");
128 136  
129 137 // Test [-1sigma, 1sigma]
130 138 result = integral.adaptive_simpson(-m_sigma, m_sigma);
  139 +/*
  140 +std::cout << result << std::endl;
  141 +std::cout << result-0.68268948130801355 << std::endl;
  142 +std::cout << integral << std::endl;
  143 +*/
131 144 test_value(result,0.68268948130801355,1.0e-6,"","Gaussian integral is not 0.682689 (difference="+gammalib::str((result-0.68268948130801355))+")");
132 145  
133 146 // Test [0.0, 1sigma]
134 147 result = integral.adaptive_simpson(0.0, m_sigma);
  148 +/*
  149 +std::cout << result << std::endl;
  150 +std::cout << result-0.3413447460687748 << std::endl;
  151 +std::cout << integral << std::endl;
  152 +*/
  153 + test_value(result,0.3413447460687748,1.0e-6,"","Gaussian integral is not 0.341345 (difference="+gammalib::str((result-0.3413447460687748))+")");
  154 +
  155 + // Return
  156 + return;
  157 +}
  158 +
  159 +
  160 +/***********************************************************************//**
  161 + * @brief Test Gauss-Kronrod integration
  162 + ***************************************************************************/
  163 +void TestGNumerics::test_gauss_kronrod_integration(void)
  164 +{
  165 + // Set-up integral
  166 + Gauss integrand(m_sigma);
  167 + GIntegral integral(&integrand);
  168 +
  169 + // Integrate over the entire Gaussian
  170 + double result = integral.gauss_kronrod(-10.0*m_sigma, 10.0*m_sigma);
  171 +/*
  172 +std::cout << result << std::endl;
  173 +std::cout << result-1.0 << std::endl;
  174 +std::cout << integral << std::endl;
  175 +*/
  176 + test_value(result,1.0,1.0e-6,"","Gaussian integral is not 1.0 (integral="+gammalib::str(result)+")");
  177 +
  178 + // Test [-1sigma, 1sigma]
  179 + result = integral.gauss_kronrod(-m_sigma, m_sigma);
  180 +/*
  181 +std::cout << result << std::endl;
  182 +std::cout << result-0.68268948130801355 << std::endl;
  183 +std::cout << integral << std::endl;
  184 +*/
  185 + test_value(result,0.68268948130801355,1.0e-6,"","Gaussian integral is not 0.682689 (difference="+gammalib::str((result-0.68268948130801355))+")");
  186 +
  187 + // Test [0.0, 1sigma]
  188 + result = integral.gauss_kronrod(0.0, m_sigma);
  189 +/*
  190 +std::cout << result << std::endl;
  191 +std::cout << result-0.3413447460687748 << std::endl;
  192 +std::cout << integral << std::endl;
  193 +*/
135 194 test_value(result,0.3413447460687748,1.0e-6,"","Gaussian integral is not 0.341345 (difference="+gammalib::str((result-0.3413447460687748))+")");
136 195  
137 196 // Return
138 197 return;
139 198 }
140 199  
  200 +
141 201 /***********************************************************************//**
142 202 * @brief Main test function
143 203 ***************************************************************************/
... ...
test/test_GNumerics.hpp
  View file @01e2aa9
... ... @@ -70,6 +70,7 @@ public:
70 70 void test_integral(void);
71 71 void test_romberg_integration(void);
72 72 void test_adaptive_simpson_integration(void);
  73 + void test_gauss_kronrod_integration(void);
73 74  
74 75 private:
75 76 // Private members
... ...