integrate_1d_gauss_kronrod.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_GAUSS_KRONROD_HPP
#define STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_GAUSS_KRONROD_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/functor/integrate_1d_adapter.hpp>
#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <algorithm>
#include <cmath>
#include <ostream>
 
namespace stan {
namespace math {
 
constexpr unsigned int INTEGRATE_1D_GAUSS_KRONROD_ORDER = 21;
 
constexpr int INTEGRATE_1D_GAUSS_KRONROD_MAX_DEPTH = 15;
 
template <typename F>
inline double integrate_gk(const F& f, double a, double b,
                           double relative_tolerance, double absolute_tolerance,
                           int max_depth) {
  static constexpr const char* function = "integrate_1d_gauss_kronrod";
  double error = 0.0;
  double L1 = 0.0;
 
  // Gauss-Kronrod does not pass a distance-to-boundary; the user functor
  // still takes (x, xc) for signature compatibility with integrate_1d, but
  // xc is unused here.
  auto f_wrap = [&f](double x) { return f(x, NOT_A_NUMBER); };
 
  using boost::math::quadrature::gauss_kronrod;
  const unsigned int depth
      = max_depth < 0 ? 0u : static_cast<unsigned int>(max_depth);
  double Q = gauss_kronrod<double, INTEGRATE_1D_GAUSS_KRONROD_ORDER>::integrate(
      f_wrap, a, b, depth, relative_tolerance, &error, &L1, absolute_tolerance);
 
  // QUADPACK-style mixed convergence: throw only if the Boost error
  // exceeds both the relative-tolerance target (rel_tol * L1) and the
  // user-supplied absolute floor. With absolute_tolerance = 0 (default)
  // this reduces to the strict pure-relative test.
  const double convergence_threshold
      = std::max(relative_tolerance * L1, absolute_tolerance);
  if (error > convergence_threshold) {
    [error]() STAN_COLD_PATH {
      throw_domain_error(
          function, "error estimate of integral", error, "",
          " exceeds max(relative_tolerance * L1, absolute_tolerance)");
    }();
  }
  return Q;
}
 
template <typename F, typename... Args,
          require_all_st_arithmetic<Args...>* = nullptr>
inline double integrate_1d_gauss_kronrod_tol(const F& f, double a, double b,
                                             double relative_tolerance,
                                             double absolute_tolerance,
                                             int max_depth, std::ostream* msgs,
                                             const Args&... args) {
  static constexpr const char* function = "integrate_1d_gauss_kronrod";
  check_less_or_equal(function, "lower limit", a, b);
  check_nonnegative(function, "max_depth", max_depth);
  check_nonnegative(function, "absolute_tolerance", absolute_tolerance);
  if (unlikely(a == b)) {
    if (std::isinf(a)) {
      throw_domain_error(function, "Integration endpoints are both", a, "", "");
    }
    return 0.0;
  } else {
    return integrate_gk(
        [&](auto&& x, auto&& xc) { return f(x, xc, msgs, args...); }, a, b,
        relative_tolerance, absolute_tolerance, max_depth);
  }
}
 
template <typename F, typename... Args,
          require_all_st_arithmetic<Args...>* = nullptr>
inline double integrate_1d_gauss_kronrod(const F& f, double a, double b,
                                         std::ostream* msgs,
                                         const Args&... args) {
  return integrate_1d_gauss_kronrod_tol(f, a, b, std::sqrt(EPSILON), 0.0,
                                        INTEGRATE_1D_GAUSS_KRONROD_MAX_DEPTH,
                                        msgs, args...);
}
 
}  // namespace math
}  // namespace stan
 
#endif