Merge pull request #2397 from stan-dev/feature/variadic-integrate-1d

Update integrate_1d to use variadic autodiff stuff internally in preparation for closures

Author: SteveBronder

stan/math/prim/functor.hpp

@@ -10,6 +10,7 @@
 #include <stan/math/prim/functor/finite_diff_gradient_auto.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
 #include <stan/math/prim/functor/integrate_1d.hpp>
+#include <stan/math/prim/functor/integrate_1d_adapter.hpp>
 #include <stan/math/prim/functor/integrate_ode_rk45.hpp>
 #include <stan/math/prim/functor/integrate_ode_std_vector_interface_adapter.hpp>
 #include <stan/math/prim/functor/ode_ckrk.hpp>

stan/math/prim/functor/integrate_1d.hpp

@@ -1,9 +1,10 @@
-#ifndef STAN_MATH_PRIM_FUNCTOR_integrate_1d_HPP
-#define STAN_MATH_PRIM_FUNCTOR_integrate_1d_HPP
+#ifndef STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_HPP
+#define STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_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>
@@ -50,17 +51,53 @@ namespace math {
 template <typename F>
 inline double integrate(const F& f, double a, double b,
                         double relative_tolerance) {
+  static constexpr const char* function = "integrate";
   double error1 = 0.0;
   double error2 = 0.0;
   double L1 = 0.0;
   double L2 = 0.0;
-  bool used_two_integrals = false;
   size_t levels;
-  double Q = 0.0;
-  auto f_wrap = [&](double x) { return f(x, NOT_A_NUMBER); };
+
+  auto one_integral_convergence_check = [](auto& error1, auto& rel_tol,
+                                           auto& L1) {
+    if (error1 > rel_tol * L1) {
+      [error1]() STAN_COLD_PATH {
+        throw_domain_error(
+            function, "error estimate of integral", error1, "",
+            " exceeds the given relative tolerance times norm of integral");
+      }();
+    }
+  };
+
+  auto two_integral_convergence_check
+      = [](auto& error1, auto& error2, auto& rel_tol, auto& L1, auto& L2) {
+          if (error1 > rel_tol * L1) {
+            [error1]() STAN_COLD_PATH {
+              throw_domain_error(
+                  function, "error estimate of integral below zero", error1, "",
+                  " exceeds the given relative tolerance times norm of "
+                  "integral below zero");
+            }();
+          }
+          if (error2 > rel_tol * L2) {
+            [error2]() STAN_COLD_PATH {
+              throw_domain_error(
+                  function, "error estimate of integral above zero", error2, "",
+                  " exceeds the given relative tolerance times norm of "
+                  "integral above zero");
+            }();
+          }
+        };
+
+  // if a or b is infinite, set xc argument to NaN (see docs above for user
+  // function for xc info)
+  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;
-    Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1, &levels);
+    double Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
+                                    &levels);
+    one_integral_convergence_check(error1, relative_tolerance, L1);
+    return Q;
   } else if (std::isinf(a)) {
     boost::math::quadrature::exp_sinh<double> integrator;
     /**
@@ -69,66 +106,88 @@ inline double integrate(const F& f, double a, double b,
      * https://www.boost.org/doc/libs/1_66_0/libs/math/doc/html/math_toolkit/double_exponential/de_caveats.html)
      */
     if (b <= 0.0) {
-      Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1,
-                               &levels);
+      double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1,
+                                      &L1, &levels);
+      one_integral_convergence_check(error1, relative_tolerance, L1);
+      return Q;
     } else {
       boost::math::quadrature::tanh_sinh<double> integrator_right;
-      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);
-      used_two_integrals = true;
+      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, L1,
+                                     L2);
+      return Q;
     }
   } else if (std::isinf(b)) {
     boost::math::quadrature::exp_sinh<double> integrator;
     if (a >= 0.0) {
-      Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1,
-                               &levels);
+      double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1,
+                                      &L1, &levels);
+      one_integral_convergence_check(error1, relative_tolerance, L1);
+      return Q;
     } else {
       boost::math::quadrature::tanh_sinh<double> integrator_left;
-      Q = integrator_left.integrate(f_wrap, a, 0, relative_tolerance, &error1,
-                                    &L1, &levels)
-          + integrator.integrate(f_wrap, relative_tolerance, &error2, &L2,
-                                 &levels);
-      used_two_integrals = true;
+      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, L1,
+                                     L2);
+      return Q;
     }
   } else {
-    auto f_wrap = [&](double x, double xc) { return f(x, xc); };
+    auto f_wrap = [&f](double x, double xc) { return f(x, xc); };
     boost::math::quadrature::tanh_sinh<double> integrator;
     if (a < 0.0 && b > 0.0) {
-      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);
-      used_two_integrals = true;
+      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, L1,
+                                     L2);
+      return Q;
     } else {
-      Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1,
-                               &levels);
+      double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1,
+                                      &L1, &levels);
+      one_integral_convergence_check(error1, relative_tolerance, L1);
+      return Q;
     }
   }
+}
 
-  static const char* function = "integrate";
-  if (used_two_integrals) {
-    if (error1 > relative_tolerance * L1) {
-      throw_domain_error(function, "error estimate of integral below zero",
-                         error1, "",
-                         " exceeds the given relative tolerance times norm of "
-                         "integral below zero");
-    }
-    if (error2 > relative_tolerance * L2) {
-      throw_domain_error(function, "error estimate of integral above zero",
-                         error2, "",
-                         " exceeds the given relative tolerance times norm of "
-                         "integral above zero");
+/**
+ * Compute the integral of the single variable function f from a to b to within
+ * a specified relative tolerance. a and b can be finite or infinite.
+ *
+ * @tparam T Type of f
+ * @param f the function to be integrated
+ * @param a lower limit of integration
+ * @param b upper limit of integration
+ * @param relative_tolerance tolerance passed to Boost quadrature
+ * @param[in, out] msgs the print stream for warning messages
+ * @param args additional arguments passed to f
+ * @return numeric integral of function f
+ */
+template <typename F, typename... Args,
+          require_all_not_st_var<Args...>* = nullptr>
+inline double integrate_1d_impl(const F& f, double a, double b,
+                                double relative_tolerance, std::ostream* msgs,
+                                const Args&... args) {
+  static constexpr const char* function = "integrate_1d";
+  check_less_or_equal(function, "lower limit", a, b);
+  if (unlikely(a == b)) {
+    if (std::isinf(a)) {
+      throw_domain_error(function, "Integration endpoints are both", a, "", "");
     }
+    return 0.0;
   } else {
-    if (error1 > relative_tolerance * L1) {
-      throw_domain_error(
-          function, "error estimate of integral", error1, "",
-          " exceeds the given relative tolerance times norm of integral");
-    }
+    return integrate(
+        [&](const auto& x, const auto& xc) { return f(x, xc, msgs, args...); },
+        a, b, relative_tolerance);
   }
-  return Q;
 }
 
 /**
@@ -178,26 +237,14 @@ inline double integrate(const F& f, double a, double b,
  * @return numeric integral of function f
  */
 template <typename F>
-inline double integrate_1d(const F& f, const double a, const double b,
+inline double integrate_1d(const F& f, double a, double b,
                            const std::vector<double>& theta,
                            const std::vector<double>& x_r,
                            const std::vector<int>& x_i, std::ostream* msgs,
                            const double relative_tolerance
                            = std::sqrt(EPSILON)) {
-  static const char* function = "integrate_1d";
-  check_less_or_equal(function, "lower limit", a, b);
-
-  if (a == b) {
-    if (std::isinf(a)) {
-      throw_domain_error(function, "Integration endpoints are both", a, "", "");
-    }
-    return 0.0;
-  } else {
-    return integrate(
-        std::bind<double>(f, std::placeholders::_1, std::placeholders::_2,
-                          theta, x_r, x_i, msgs),
-        a, b, relative_tolerance);
-  }
+  return integrate_1d_impl(integrate_1d_adapter<F>(f), a, b, relative_tolerance,
+                           msgs, theta, x_r, x_i);
 }
 
 }  // namespace math

stan/math/prim/functor/integrate_1d_adapter.hpp

@@ -0,0 +1,28 @@
+#ifndef STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_ADAPTER_HPP
+#define STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_ADAPTER_HPP
+
+#include <ostream>
+#include <vector>
+
+/**
+ * Adapt the non-variadic integrate_1d arguments to the variadic
+ * integrate_1d_impl interface
+ *
+ * @tparam F type of function to adapt
+ */
+template <typename F>
+struct integrate_1d_adapter {
+  const F& f_;
+
+  explicit integrate_1d_adapter(const F& f) : f_(f) {}
+
+  template <typename T_a, typename T_b, typename T_theta>
+  auto operator()(const T_a& x, const T_b& xc, std::ostream* msgs,
+                  const std::vector<T_theta>& theta,
+                  const std::vector<double>& x_r,
+                  const std::vector<int>& x_i) const {
+    return f_(x, xc, theta, x_r, x_i, msgs);
+  }
+};
+
+#endif