NumericalIntegration.h

Go to the documentation of this file.

#ifndef  AMPGEN_NUMERICALINTEGRATION_H 
#define  AMPGEN_NUMERICALINTEGRATION_H 1 
#include <iostream>
#include <queue>
 
#include <gsl/gsl_integration.h>
 
#include "AmpGen/MetaUtils.h"
#include "AmpGen/simd/utils.h"
#include "AmpGen/simd/integrate_fp.h"
 
namespace AmpGen {
  template <typename function_type> double integrate_1d( const function_type& fcn, const double& min, const double& max, gsl_integration_workspace* ws = nullptr)
  {
    gsl_integration_workspace* w = (ws == nullptr) ? gsl_integration_workspace_alloc (1000) : ws;
    double result = 0;
    double error  = 0;
    gsl_function F;
    //  static_assert( is_functor<function_type, double( const double& )>::value == std::true_type , "function is of wrong type");
    if constexpr( is_functor<function_type, double( const double& )>::value )
    { 
      F.function = [](double x, void* p) -> double { return (*static_cast<function_type*>(p))(x); };
    }
    else if constexpr( is_functor<function_type, double( const std::array<double, 1>& )>::value )
    {
      F.function = [](double x, void* p) -> double { return (*static_cast<function_type*>(p))( std::array<double,1>{x} ); };
    }
    else static_assert( true, "function matches no signature!");
 
    F.params   = const_cast<function_type*>(&fcn);
    gsl_integration_qags (&F, min, max, 0, 1e-5, 1000, w, &result, &error);
    if( ws == nullptr ) gsl_integration_workspace_free (w);
    return result;
  }
  template <typename function_type> double integrate_1d_inf( const function_type& fcn, gsl_integration_workspace* ws = nullptr)
  {
    gsl_integration_workspace* w = (ws == nullptr) ? gsl_integration_workspace_alloc (1000) : ws;
    double result = 0;
    double error  = 0;
    gsl_function F;
    F.function = [](double x, void* p) -> double { return (*static_cast<function_type*>(p))(x); };
    F.params   = const_cast<function_type*>(&fcn);
    gsl_integration_qagi(&F, 0, 1e-5, 1000, w, &result, &error);
    if( ws == nullptr ) gsl_integration_workspace_free (w);
    //std::cout << result << " +/- " << error << std::endl; 
    return result;
  }
 
  template <typename function_type> double integrate_1d_cauchy( const function_type& fcn, const double& x0, const double& min, const double& max, gsl_integration_workspace* ws = nullptr)
  {
    gsl_integration_workspace* w = (ws == nullptr) ? gsl_integration_workspace_alloc (1000) : ws;
    double result = 0;
    double error  = 0;
    gsl_function F;
    F.function = [](double x, void* p) -> double { return (*static_cast<function_type*>(p))(x); };
    F.params   = const_cast<function_type*>(&fcn);
    gsl_integration_qawc(&F, min, max, x0, 0, 1e-8, 1000, w, &result, &error);
    if( ws == nullptr ) gsl_integration_workspace_free (w);
    return result;
  }
 
  template <unsigned dim> struct integral {
    double value = {0};
    double var   = {0};
    int index = {0};
    std::array<double, dim> a     = {0};
    std::array<double, dim> b     = {0};
    bool operator<( const integral<dim>& other) const { return var < other.var; }
  };  
 
  template<unsigned dim, typename fcn> std::tuple<double, double, unsigned> integrate_fp( const fcn& F, const std::array<double,dim>& ctr, const std::array<double, dim>& wth )
  {
    #if INSTRUCTION_SET == INSTRUCTION_SET_AVX2d
      return integrate_fp<dim, 4>(F, ctr, wth);
    #endif
    if constexpr( is_functor<fcn, double(const std::array<double, dim>&)>::value )
    {
      return integrate_fp_scalar<dim>(F, ctr, wth);
    }
  }
 
  template <unsigned dim, typename fcn> double integrate(const fcn& F, const std::array<double, dim> xmin, const std::array<double, dim>&  xmax)
  {
    // Based off of ROOT Math/IntegratorMultiDim
    // With improved speed and safety:
    // - No dynamic memory allocation
    // - Static dimension of integrals 
    // - Supports vectorised integrands  
    // References to actual method [again, from ROOT Math/IntegratorMultiDim]
    //   1.A.C. Genz and A.A. Malik, Remarks on algorithm 006:
    //     An adaptive algorithm for numerical integration over
    //     an N-dimensional rectangular region, J. Comput. Appl. Math. 6 (1980) 295-302.
    //   2.A. van Doren and L. de Ridder, An adaptive algorithm for numerical
    //     integration over an n-dimensional cube, J.Comput. Appl. Math. 2 (1976) 207-217.
    if constexpr( dim == 1 ){
      if constexpr( is_functor<fcn, double(const double&)>::value ) { integrate_1d(F, xmin[0], xmax[0] ); }
      else if constexpr( is_functor<fcn, double(const std::array<double, 1>& ) >::value ){
        return integrate_1d( [&F](const double& x){ return F( std::array<double,1>{x} ); }, xmin[0], xmax[0] );
      }
      else static_assert(true, "1D function doesn't have recognised signature");
      return 0;
    }
    else {
 
      double epsrel = 1e-10;  //specified relative accuracy
      double epsabs = 0.; //specified relative accuracy
      //output parameters
      double relerr = 0 ; //an estimation of the relative accuracy of the result
 
      double result = 0;
      double abserr = 0;
      auto status  = 3;
 
      unsigned int ifncls = 0; 
      bool  ldv   = false;
      unsigned isbrgn = 1;
      unsigned isbrgs = 1;
 
      constexpr unsigned irlcls = get_power<2,dim>::value +2*dim*(dim+1)+1; // number of function evaluations per iteration
      constexpr unsigned minpts = get_power<2,dim>::value +2*dim*(dim+1)+1; // minimum number of function evaluations
      constexpr unsigned maxpts = 1000000;                                   // maximum number of function evaluations
 
      std::array<integral<dim>, ( 1 + maxpts/irlcls)/2 > partial_integrals; 
 
      integral<dim>* current_integral = &(partial_integrals[0]); 
 
      for (unsigned j=0; j<dim; j++) {
        current_integral->a[j] = (xmax[j] + xmin[j])*0.5;
        current_integral->b[j] = (xmax[j] - xmin[j])*0.5;
      }
 
      unsigned int idvax0=0; 
      integral<dim> tmp; 
 
      do {
        auto [rgnval, rgnerr, idvaxn] = integrate_fp<dim>( F, current_integral->a, current_integral->b );
        result += rgnval;
        abserr += rgnerr;
        ifncls += irlcls;
 
        double aresult = std::abs(result);
        tmp.var   = rgnerr;
        tmp.value = rgnval;
        tmp.index = idvaxn;
        tmp.a     = current_integral->a;
        tmp.b     = current_integral->b;
 
        if (ldv) {
          unsigned isbtmp =0;  
          while( true )
          {
            isbtmp = 2*isbrgn; 
            if (isbtmp > isbrgs) break;
            if (isbtmp < isbrgs && partial_integrals[isbtmp].var < partial_integrals[isbtmp+1].var ) isbtmp++;
            if (rgnerr >= partial_integrals[isbtmp].var ) break;
            partial_integrals[isbrgn] = partial_integrals[isbtmp] ;
            isbrgn = isbtmp;
          }
        }
        else {
          unsigned isbtmp =0;  
          do {
            isbtmp = isbrgn/2;
            if ( isbtmp >= 1 && rgnerr > partial_integrals[isbtmp].var ) {
              partial_integrals[isbrgn] = partial_integrals[isbtmp]; 
              isbrgn = isbtmp;
            }
          } while ( isbtmp >= 1 && rgnerr > partial_integrals[isbtmp].var );
        }
 
        partial_integrals[isbrgn] = tmp; 
        if (ldv) {//divison along chosen coordinate
          ldv = false; 
          current_integral = &(partial_integrals[isbrgn]);
          isbrgs += 1; 
          isbrgn  = isbrgs;
          partial_integrals[isbrgn].a     = current_integral->a;
          partial_integrals[isbrgn].b     = current_integral->b;
          partial_integrals[isbrgn].a[idvax0] += 2*current_integral->b[idvax0];
          current_integral = &(partial_integrals[isbrgn]);
          continue; 
        }
        //if no divisions to be made..
        relerr = std::abs(result) == 0 ? abserr : abserr/std::abs(result);
        if ((relerr < epsrel && aresult < epsabs) or 
            ( ( relerr < epsrel || abserr < epsabs ) && ifncls > minpts) ){ status = 0; break ; }
        if (( isbrgs >= partial_integrals.size()-1) or ( ifncls+2*irlcls > maxpts ) )  { status = 2; break; }
       // std::cout << "#calls: " << ifncls << ", #cycles: " << isbrgs << " / " << ( 1 + maxpts/irlcls)/2 << " " << abserr << std::endl; 
        ldv = true;
        isbrgn  = 1;
        current_integral = &( partial_integrals[isbrgn] );
 
        abserr -= current_integral->var;
        result -= current_integral->value;
        idvax0  = current_integral->index;
 
        current_integral -> b[idvax0] *= 0.5;
        current_integral -> a[idvax0] -= current_integral->b[idvax0];
 
      } while( status == 3 ); 
//      std::cout << "IFNCLS: " << ifncls << " ISBRGN: " << isbrgn  << " ISBRGS: " << isbrgs << " " << abserr << " " << relerr << std::endl; 
      return result;         //an approximate value of the integral
    }
  }
}
 
 
#endif