    /*****************************************************************************
     * Project: RooFit                                                           *
     * Package: RooFitCore                                                       *
     *    File: $Id: RooGaussKronrodIntegrator1D.cc,v 1.6 2005/06/20 15:44:52 wverkerke Exp $
     * Authors:                                                                  *
     *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
     *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
     *                                                                           *
     * Copyright (c) 2000-2005, Regents of the University of California          *
     *                          and Stanford University. All rights reserved.    *
     *                                                                           *
     * Redistribution and use in source and binary forms,                        *
     * with or without modification, are permitted according to the terms        *
     * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
     *****************************************************************************/

    // -- CLASS DESCRIPTION [AUX] --
    // RooGaussKronrodIntegrator1D implements the Gauss-Kronrod integration algorithm.
    //
    // An Gaussian quadrature method for numerical integration in which
    // error is estimation based on evaluation at special points known as
    // "Kronrod points."  By suitably picking these points, abscissas from
    // previous iterations can be reused as part of the new set of points,
    // whereas usual Gaussian quadrature would require recomputation of
    // all abscissas at each iteration.
    //
    // This class automatically handles (-inf,+inf) integrals by dividing
    // the integration in three regions (-inf,-1), (-1,1), (1,inf) and
    // calculating the 1st and 3rd term using a x -> 1/x coordinate
    // transformation
    //
    // This class embeds the Gauss-Kronrod integrator from the GNU
    // Scientific Library version 1.5 and applies the 10-, 21-, 43- and
    // 87-point rule in succession until the required target precision is
    // reached

    #include "RooFit.h"

    #include <assert.h>
    #include <assert.h>
    #include "Riostream.h"
    #include "TMath.h"
    #include "RooGaussKronrodIntegrator1D.h"
    #include "RooArgSet.h"
    #include "RooRealVar.h"
    #include "RooNumber.h"
    #include "RooNumIntFactory.h"
    #include "RooIntegratorBinding.h"

    ClassImp(RooGaussKronrodIntegrator1D)
    ;

    // --- From GSL_MATH.h -------------------------------------------
    struct gsl_function_struct
    {
      double (* function) (double x, void * params);
      void * params;
    };
    typedef struct gsl_function_struct gsl_function ;
    #define GSL_FN_EVAL(F,x) (*((F)->function))(x,(F)->params)
    //----From GSL_INTEGRATION.h ---------------------------------------
    int gsl_integration_qng (const gsl_function * f,
                             double a, double b,
                             double epsabs, double epsrel,
                             double *result, double *abserr,
                             size_t * neval);
    //-------------------------------------------------------------------

    // Register this class with RooNumIntConfig
    static void registerGaussKronrodIntegrator1D(RooNumIntFactory& fact)
    {
      fact.storeProtoIntegrator(new RooGaussKronrodIntegrator1D(),RooArgSet()) ;
    }
    static Bool_t dummy = RooNumIntFactory::instance().registerInitializer(&registerGaussKronrodIntegrator1D) ;

     RooGaussKronrodIntegrator1D::RooGaussKronrodIntegrator1D()
    {
    }

     RooGaussKronrodIntegrator1D::RooGaussKronrodIntegrator1D(const RooAbsFunc& function, const RooNumIntConfig& config) :
      RooAbsIntegrator(function),
      _epsAbs(config.epsRel()),
      _epsRel(config.epsAbs())
    {
      // Use this form of the constructor to integrate over the function's default range.
      _useIntegrandLimits= kTRUE;
      _valid= initialize();
    }

     RooGaussKronrodIntegrator1D::RooGaussKronrodIntegrator1D(const RooAbsFunc& function,
    							 Double_t xmin, Double_t xmax, const RooNumIntConfig& config) :
      RooAbsIntegrator(function),
      _epsAbs(config.epsRel()),
      _epsRel(config.epsAbs()),
      _xmin(xmin),
      _xmax(xmax)
    {
      // Use this form of the constructor to integrate over the function's default range.
      _useIntegrandLimits= kFALSE;
      _valid= initialize();
    }

     RooAbsIntegrator* RooGaussKronrodIntegrator1D::clone(const RooAbsFunc& function, const RooNumIntConfig& config) const
    {
      return new RooGaussKronrodIntegrator1D(function,config) ;
    }

     Bool_t RooGaussKronrodIntegrator1D::initialize()
    {
      // Allocate coordinate buffer size after number of function dimensions
      _x = new Double_t[_function->getDimension()] ;

      return checkLimits();
    }

     RooGaussKronrodIntegrator1D::~RooGaussKronrodIntegrator1D()
    {
      delete[] _x ;
    }

     Bool_t RooGaussKronrodIntegrator1D::setLimits(Double_t xmin, Double_t xmax) {
      // Change our integration limits. Return kTRUE if the new limits are
      // ok, or otherwise kFALSE. Always returns kFALSE and does nothing
      // if this object was constructed to always use our integrand's limits.

      if(_useIntegrandLimits) {
        cout << "RooGaussKronrodIntegrator1D::setLimits: cannot override integrand's limits" << endl;
        return kFALSE;
      }
      _xmin= xmin;
      _xmax= xmax;
      return checkLimits();
    }

     Bool_t RooGaussKronrodIntegrator1D::checkLimits() const {
      // Check that our integration range is finite and otherwise return kFALSE.
      // Update the limits from the integrand if requested.

      if(_useIntegrandLimits) {
        assert(0 != integrand() && integrand()->isValid());
        _xmin= integrand()->getMinLimit(0);
        _xmax= integrand()->getMaxLimit(0);
      }
      return kTRUE ;
    }

    double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
    {
      RooGaussKronrodIntegrator1D* instance = (RooGaussKronrodIntegrator1D*) data ;
      return instance->integrand(instance->xvec(x)) ;
    }

     Double_t RooGaussKronrodIntegrator1D::integral(const Double_t *yvec)
    {
      assert(isValid());

      // Copy yvec to xvec if provided
      if (yvec) {
        UInt_t i ; for (i=0 ; i<_function->getDimension()-1 ; i++) {
          _x[i+1] = yvec[i] ;
        }
      }

      // Setup glue function
      gsl_function F;
      F.function = &RooGaussKronrodIntegrator1D_GSL_GlueFunction ;
      F.params = this ;

      // Return values
      double result, error;
      size_t neval = 0 ;

      // Call GSL implementation of integeator
      gsl_integration_qng (&F, _xmin, _xmax, _epsAbs, _epsRel, &result, &error, &neval);

      return result;
    }

    // ----------------------------------------------------------------------------
    // ---------- Code below imported from GSL version 1.5 ------------------------
    // ----------------------------------------------------------------------------

    /*
     *
     * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
     *
     * This program is free software; you can redistribute it and/or modify
     * it under the terms of the GNU General Public License as published by
     * the Free Software Foundation; either version 2 of the License, or (at
     * your option) any later version.
     *
     * This program is distributed in the hope that it will be useful, but
     * WITHOUT ANY WARRANTY; without even the implied warranty of
     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     * General Public License for more details.
     *
     * You should have received a copy of the GNU General Public License
     * along with this program; if not, write to the Free Software
     * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
     */

    #define GSL_SUCCESS 0
    #define GSL_EBADTOL 13  /* user specified an invalid tolerance */
    #define GSL_ETOL    14  /* failed to reach the specified tolerance */
    #define GSL_ERROR(a,b) cout << "RooGaussKronrodIntegrator1D::integral() ERROR: " << a << endl ; return b ;
    #define GSL_DBL_MIN        2.2250738585072014e-308
    #define GSL_DBL_EPSILON    2.2204460492503131e-16

    // INCLUDED BELOW #include "qng.c"
    #include <config.h>
    #include <math.h>
    #include <float.h>
    #include <stdlib.h>

    int gsl_integration_qng (const gsl_function * f,
                             double a, double b,
                             double epsabs, double epsrel,
                             double *result, double *abserr,
                             size_t * neval);

    // INCLUDED BELOW #include "err.c"
    static double rescale_error (double err, const double result_abs, const double result_asc) ;

    static double
    rescale_error (double err, const double result_abs, const double result_asc)
    {
      err = fabs(err) ;

      if (result_asc != 0 && err != 0)
          {
            double scale = TMath::Power((200 * err / result_asc), 1.5) ;

            if (scale < 1)
              {
                err = result_asc * scale ;
              }
            else
              {
                err = result_asc ;
              }
          }
      if (result_abs > GSL_DBL_MIN / (50 * GSL_DBL_EPSILON))
        {
          double min_err = 50 * GSL_DBL_EPSILON * result_abs ;

          if (min_err > err)
            {
              err = min_err ;
            }
        }

      return err ;
    }

    // INCLUDED BELOW #include "qng.h"
    /* Gauss-Kronrod-Patterson quadrature coefficients for use in
       quadpack routine qng. These coefficients were calculated with
       101 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov
       1981. */

    /* x1, abscissae common to the 10-, 21-, 43- and 87-point rule */
    static const double x1[5] = {
      0.973906528517171720077964012084452,
      0.865063366688984510732096688423493,
      0.679409568299024406234327365114874,
      0.433395394129247190799265943165784,
      0.148874338981631210884826001129720
    } ;

    /* w10, weights of the 10-point formula */
    static const double w10[5] = {
      0.066671344308688137593568809893332,
      0.149451349150580593145776339657697,
      0.219086362515982043995534934228163,
      0.269266719309996355091226921569469,
      0.295524224714752870173892994651338
    } ;

    /* x2, abscissae common to the 21-, 43- and 87-point rule */
    static const double x2[5] = {
      0.995657163025808080735527280689003,
      0.930157491355708226001207180059508,
      0.780817726586416897063717578345042,
      0.562757134668604683339000099272694,
      0.294392862701460198131126603103866
    } ;

    /* w21a, weights of the 21-point formula for abscissae x1 */
    static const double w21a[5] = {
      0.032558162307964727478818972459390,
      0.075039674810919952767043140916190,
      0.109387158802297641899210590325805,
      0.134709217311473325928054001771707,
      0.147739104901338491374841515972068
    } ;

    /* w21b, weights of the 21-point formula for abscissae x2 */
    static const double w21b[6] = {
      0.011694638867371874278064396062192,
      0.054755896574351996031381300244580,
      0.093125454583697605535065465083366,
      0.123491976262065851077958109831074,
      0.142775938577060080797094273138717,
      0.149445554002916905664936468389821
    } ;

    /* x3, abscissae common to the 43- and 87-point rule */
    static const double x3[11] = {
      0.999333360901932081394099323919911,
      0.987433402908088869795961478381209,
      0.954807934814266299257919200290473,
      0.900148695748328293625099494069092,
      0.825198314983114150847066732588520,
      0.732148388989304982612354848755461,
      0.622847970537725238641159120344323,
      0.499479574071056499952214885499755,
      0.364901661346580768043989548502644,
      0.222254919776601296498260928066212,
      0.074650617461383322043914435796506
    } ;

    /* w43a, weights of the 43-point formula for abscissae x1, x3 */
    static const double w43a[10] = {
      0.016296734289666564924281974617663,
      0.037522876120869501461613795898115,
      0.054694902058255442147212685465005,
      0.067355414609478086075553166302174,
      0.073870199632393953432140695251367,
      0.005768556059769796184184327908655,
      0.027371890593248842081276069289151,
      0.046560826910428830743339154433824,
      0.061744995201442564496240336030883,
      0.071387267268693397768559114425516
    } ;

    /* w43b, weights of the 43-point formula for abscissae x3 */
    static const double w43b[12] = {
      0.001844477640212414100389106552965,
      0.010798689585891651740465406741293,
      0.021895363867795428102523123075149,
      0.032597463975345689443882222526137,
      0.042163137935191811847627924327955,
      0.050741939600184577780189020092084,
      0.058379395542619248375475369330206,
      0.064746404951445885544689259517511,
      0.069566197912356484528633315038405,
      0.072824441471833208150939535192842,
      0.074507751014175118273571813842889,
      0.074722147517403005594425168280423
    } ;

    /* x4, abscissae of the 87-point rule */
    static const double x4[22] = {
      0.999902977262729234490529830591582,
      0.997989895986678745427496322365960,
      0.992175497860687222808523352251425,
      0.981358163572712773571916941623894,
      0.965057623858384619128284110607926,
      0.943167613133670596816416634507426,
      0.915806414685507209591826430720050,
      0.883221657771316501372117548744163,
      0.845710748462415666605902011504855,
      0.803557658035230982788739474980964,
      0.757005730685495558328942793432020,
      0.706273209787321819824094274740840,
      0.651589466501177922534422205016736,
      0.593223374057961088875273770349144,
      0.531493605970831932285268948562671,
      0.466763623042022844871966781659270,
      0.399424847859218804732101665817923,
      0.329874877106188288265053371824597,
      0.258503559202161551802280975429025,
      0.185695396568346652015917141167606,
      0.111842213179907468172398359241362,
      0.037352123394619870814998165437704
    } ;

    /* w87a, weights of the 87-point formula for abscissae x1, x2, x3 */
    static const double w87a[21] = {
      0.008148377384149172900002878448190,
      0.018761438201562822243935059003794,
      0.027347451050052286161582829741283,
      0.033677707311637930046581056957588,
      0.036935099820427907614589586742499,
      0.002884872430211530501334156248695,
      0.013685946022712701888950035273128,
      0.023280413502888311123409291030404,
      0.030872497611713358675466394126442,
      0.035693633639418770719351355457044,
      0.000915283345202241360843392549948,
      0.005399280219300471367738743391053,
      0.010947679601118931134327826856808,
      0.016298731696787335262665703223280,
      0.021081568889203835112433060188190,
      0.025370969769253827243467999831710,
      0.029189697756475752501446154084920,
      0.032373202467202789685788194889595,
      0.034783098950365142750781997949596,
      0.036412220731351787562801163687577,
      0.037253875503047708539592001191226
    } ;

    /* w87b, weights of the 87-point formula for abscissae x4    */
    static const double w87b[23] = {
      0.000274145563762072350016527092881,
      0.001807124155057942948341311753254,
      0.004096869282759164864458070683480,
      0.006758290051847378699816577897424,
      0.009549957672201646536053581325377,
      0.012329447652244853694626639963780,
      0.015010447346388952376697286041943,
      0.017548967986243191099665352925900,
      0.019938037786440888202278192730714,
      0.022194935961012286796332102959499,
      0.024339147126000805470360647041454,
      0.026374505414839207241503786552615,
      0.028286910788771200659968002987960,
      0.030052581128092695322521110347341,
      0.031646751371439929404586051078883,
      0.033050413419978503290785944862689,
      0.034255099704226061787082821046821,
      0.035262412660156681033782717998428,
      0.036076989622888701185500318003895,
      0.036698604498456094498018047441094,
      0.037120549269832576114119958413599,
      0.037334228751935040321235449094698,
      0.037361073762679023410321241766599
    } ;

    int
    gsl_integration_qng (const gsl_function *f,
                         double a, double b,
                         double epsabs, double epsrel,
                         double * result, double * abserr, size_t * neval)
    {
      double fv1[5], fv2[5], fv3[5], fv4[5];
      double savfun[21];  /* array of function values which have been computed */
      double res10, res21, res43, res87;    /* 10, 21, 43 and 87 point results */
      double result_kronrod, err ;
      double resabs; /* approximation to the integral of abs(f) */
      double resasc; /* approximation to the integral of abs(f-i/(b-a)) */

      const double half_length =  0.5 * (b - a);
      const double abs_half_length = fabs (half_length);
      const double center = 0.5 * (b + a);
      const double f_center = GSL_FN_EVAL(f, center);

      int k ;

      if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
        {
          * result = 0;
          * abserr = 0;
          * neval = 0;
          GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
    		 GSL_EBADTOL);
        };

      /* Compute the integral using the 10- and 21-point formula. */

      res10 = 0;
      res21 = w21b[5] * f_center;
      resabs = w21b[5] * fabs (f_center);

      for (k = 0; k < 5; k++)
        {
          const double abscissa = half_length * x1[k];
          const double fval1 = GSL_FN_EVAL(f, center + abscissa);
          const double fval2 = GSL_FN_EVAL(f, center - abscissa);
          const double fval = fval1 + fval2;
          res10 += w10[k] * fval;
          res21 += w21a[k] * fval;
          resabs += w21a[k] * (fabs (fval1) + fabs (fval2));
          savfun[k] = fval;
          fv1[k] = fval1;
          fv2[k] = fval2;
        }

      for (k = 0; k < 5; k++)
        {
          const double abscissa = half_length * x2[k];
          const double fval1 = GSL_FN_EVAL(f, center + abscissa);
          const double fval2 = GSL_FN_EVAL(f, center - abscissa);
          const double fval = fval1 + fval2;
          res21 += w21b[k] * fval;
          resabs += w21b[k] * (fabs (fval1) + fabs (fval2));
          savfun[k + 5] = fval;
          fv3[k] = fval1;
          fv4[k] = fval2;
        }

      resabs *= abs_half_length ;

      {
        const double mean = 0.5 * res21;

        resasc = w21b[5] * fabs (f_center - mean);

        for (k = 0; k < 5; k++)
          {
            resasc +=
              (w21a[k] * (fabs (fv1[k] - mean) + fabs (fv2[k] - mean))
              + w21b[k] * (fabs (fv3[k] - mean) + fabs (fv4[k] - mean)));
          }
        resasc *= abs_half_length ;
      }

      result_kronrod = res21 * half_length;

      err = rescale_error ((res21 - res10) * half_length, resabs, resasc) ;

      /*   test for convergence. */

      if (err < epsabs || err < epsrel * fabs (result_kronrod))
        {
          * result = result_kronrod ;
          * abserr = err ;
          * neval = 21;
          return GSL_SUCCESS;
        }

      /* compute the integral using the 43-point formula. */

      res43 = w43b[11] * f_center;

      for (k = 0; k < 10; k++)
        {
          res43 += savfun[k] * w43a[k];
        }

      for (k = 0; k < 11; k++)
        {
          const double abscissa = half_length * x3[k];
          const double fval = (GSL_FN_EVAL(f, center + abscissa)
                               + GSL_FN_EVAL(f, center - abscissa));
          res43 += fval * w43b[k];
          savfun[k + 10] = fval;
        }

      /*  test for convergence */

      result_kronrod = res43 * half_length;
      err = rescale_error ((res43 - res21) * half_length, resabs, resasc);

      if (err < epsabs || err < epsrel * fabs (result_kronrod))
        {
          * result = result_kronrod ;
          * abserr = err ;
          * neval = 43;
          return GSL_SUCCESS;
        }

      /* compute the integral using the 87-point formula. */

      res87 = w87b[22] * f_center;

      for (k = 0; k < 21; k++)
        {
          res87 += savfun[k] * w87a[k];
        }

      for (k = 0; k < 22; k++)
        {
          const double abscissa = half_length * x4[k];
          res87 += w87b[k] * (GSL_FN_EVAL(f, center + abscissa)
                              + GSL_FN_EVAL(f, center - abscissa));
        }

      /*  test for convergence */

      result_kronrod = res87 * half_length ;

      err = rescale_error ((res87 - res43) * half_length, resabs, resasc);

      if (err < epsabs || err < epsrel * fabs (result_kronrod))
        {
          * result = result_kronrod ;
          * abserr = err ;
          * neval = 87;
          return GSL_SUCCESS;
        }

      /* failed to converge */

      * result = result_kronrod ;
      * abserr = err ;
      * neval = 87;

      // GSL_ERROR("failed to reach tolerance with highest-order rule", GSL_ETOL) ;
      return GSL_ETOL ;
    }