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

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,44 @@ 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 one_integral_convergence_check = [&]() {
+    if (error1 > relative_tolerance * L1) {
+      throw_domain_error(
+          function, "error estimate of integral", error1, "",
+          " exceeds the given relative tolerance times norm of integral");
+    }
+  };
+
+  auto two_integral_convergence_check = [&]() {
+    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");
+    }
+  };
+
+  // if a or b is infinite, set xc argument to NaN (see docs above for user
+  // function for xc info)
   auto f_wrap = [&](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);
+    one_integral_convergence_check();
   } else if (std::isinf(a)) {
     boost::math::quadrature::exp_sinh<double> integrator;
     /**
@@ -71,26 +99,28 @@ inline double integrate(const F& f, double a, double b,
     if (b <= 0.0) {
       Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1,
                                &levels);
+      one_integral_convergence_check();
     } 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;
+      two_integral_convergence_check();
     }
   } 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);
+      one_integral_convergence_check();
     } 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;
+      two_integral_convergence_check();
     }
   } else {
     auto f_wrap = [&](double x, double xc) { return f(x, xc); };
@@ -100,35 +130,48 @@ inline double integrate(const F& f, double a, double b,
                                &levels)
           + integrator.integrate(f_wrap, 0.0, b, relative_tolerance, &error2,
                                  &L2, &levels);
-      used_two_integrals = true;
+      two_integral_convergence_check();
     } else {
       Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1,
                                &levels);
+      one_integral_convergence_check();
     }
   }
 
-  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");
+  return Q;
+}
+
+/**
+ * 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 (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 +221,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