integrate_1d_double_exponential.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_DOUBLE_EXPONENTIAL_HPP
#define STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_DOUBLE_EXPONENTIAL_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/exp_sinh.hpp>
#include <boost/math/quadrature/sinh_sinh.hpp>
#include <boost/math/quadrature/tanh_sinh.hpp>
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <ostream>
 
namespace stan {
namespace math {
 
constexpr int INTEGRATE_1D_DOUBLE_EXPONENTIAL_MAX_REFINEMENTS = 15;
 
template <typename F>
inline double integrate_de(const F& f, double a, double b,
                           double relative_tolerance, double absolute_tolerance,
                           int max_refinements) {
  static constexpr const char* function = "integrate_1d_double_exponential";
  double error1 = 0.0;
  double error2 = 0.0;
  double L1 = 0.0;
  double L2 = 0.0;
  size_t levels = 0;
 
  const size_t mr
      = max_refinements < 0 ? 0u : static_cast<size_t>(max_refinements);
 
  auto one_integral_convergence_check
      = [&](double e, double rel_tol, double abs_tol, double L) {
          const double threshold = std::max(rel_tol * L, abs_tol);
          if (e > threshold) {
            [e]() STAN_COLD_PATH {
              throw_domain_error(
                  function, "error estimate of integral", e, "",
                  " exceeds max(relative_tolerance * L1, absolute_tolerance)");
            }();
          }
        };
  auto two_integral_convergence_check = [&](double e1, double e2,
                                            double rel_tol, double abs_tol,
                                            double La, double Lb) {
    const double threshold_a = std::max(rel_tol * La, abs_tol);
    const double threshold_b = std::max(rel_tol * Lb, abs_tol);
    if (e1 > threshold_a) {
      [e1]() STAN_COLD_PATH {
        throw_domain_error(
            function, "error estimate of integral below zero", e1, "",
            " exceeds max(relative_tolerance * L1, absolute_tolerance) for "
            "the lower piece of the split");
      }();
    }
    if (e2 > threshold_b) {
      [e2]() STAN_COLD_PATH {
        throw_domain_error(
            function, "error estimate of integral above zero", e2, "",
            " exceeds max(relative_tolerance * L1, absolute_tolerance) for "
            "the upper piece of the split");
      }();
    }
  };
 
  // If a or b is infinite, set xc to NaN (see docs above)
  auto f_wrap = [&f](double x) { return f(x, NOT_A_NUMBER); };
  if (std::isinf(a) && std::isinf(b)) {
    boost::math::quadrature::sinh_sinh<double> integrator(mr);
    double Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
                                    &levels);
    one_integral_convergence_check(error1, relative_tolerance,
                                   absolute_tolerance, L1);
    return Q;
  } else if (std::isinf(a)) {
    boost::math::quadrature::exp_sinh<double> integrator(mr);
    // If the integral crosses zero, break it into two (advice from the
    // Boost implementation).
    if (b <= 0.0) {
      double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1,
                                      &L1, &levels);
      one_integral_convergence_check(error1, relative_tolerance,
                                     absolute_tolerance, L1);
      return Q;
    } else {
      boost::math::quadrature::tanh_sinh<double> integrator_right(mr);
      double Q = integrator.integrate(f_wrap, a, 0.0, relative_tolerance,
                                      &error1, &L1, &levels)
                 + integrator_right.integrate(
                     f_wrap, 0.0, b, relative_tolerance, &error2, &L2, &levels);
      two_integral_convergence_check(error1, error2, relative_tolerance,
                                     absolute_tolerance, L1, L2);
      return Q;
    }
  } else if (std::isinf(b)) {
    boost::math::quadrature::exp_sinh<double> integrator(mr);
    if (a >= 0.0) {
      double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1,
                                      &L1, &levels);
      one_integral_convergence_check(error1, relative_tolerance,
                                     absolute_tolerance, L1);
      return Q;
    } else {
      boost::math::quadrature::tanh_sinh<double> integrator_left(mr);
      double Q = integrator_left.integrate(f_wrap, a, 0, relative_tolerance,
                                           &error1, &L1, &levels)
                 + integrator.integrate(f_wrap, relative_tolerance, &error2,
                                        &L2, &levels);
      two_integral_convergence_check(error1, error2, relative_tolerance,
                                     absolute_tolerance, L1, L2);
      return Q;
    }
  } else {
    auto f_wrap = [&f](double x, double xc) { return f(x, xc); };
    boost::math::quadrature::tanh_sinh<double> integrator(mr);
    if (a < 0.0 && b > 0.0) {
      double Q = integrator.integrate(f_wrap, a, 0.0, relative_tolerance,
                                      &error1, &L1, &levels)
                 + integrator.integrate(f_wrap, 0.0, b, relative_tolerance,
                                        &error2, &L2, &levels);
      two_integral_convergence_check(error1, error2, relative_tolerance,
                                     absolute_tolerance, L1, L2);
      return Q;
    } else {
      double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1,
                                      &L1, &levels);
      one_integral_convergence_check(error1, relative_tolerance,
                                     absolute_tolerance, L1);
      return Q;
    }
  }
}
 
template <typename F, typename... Args,
          require_all_st_arithmetic<Args...>* = nullptr>
inline double integrate_1d_double_exponential_tol(
    const F& f, double a, double b, double relative_tolerance,
    double absolute_tolerance, int max_refinements, std::ostream* msgs,
    const Args&... args) {
  static constexpr const char* function = "integrate_1d_double_exponential";
  check_less_or_equal(function, "lower limit", a, b);
  check_nonnegative(function, "max_refinements", max_refinements);
  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_de(
        [&](auto&& x, auto&& xc) { return f(x, xc, msgs, args...); }, a, b,
        relative_tolerance, absolute_tolerance, max_refinements);
  }
}
 
template <typename F, typename... Args,
          require_all_st_arithmetic<Args...>* = nullptr>
inline double integrate_1d_double_exponential(const F& f, double a, double b,
                                              std::ostream* msgs,
                                              const Args&... args) {
  return integrate_1d_double_exponential_tol(
      f, a, b, std::sqrt(EPSILON), 0.0,
      INTEGRATE_1D_DOUBLE_EXPONENTIAL_MAX_REFINEMENTS, msgs, args...);
}
 
}  // namespace math
}  // namespace stan
 
#endif