add integrate_1d_gauss_kronrod and integrate_1d_double_exponential

stan/math/fwd/functor.hpp

@@ -6,6 +6,8 @@
 #include <stan/math/fwd/functor/finite_diff.hpp>
 #include <stan/math/fwd/functor/hessian.hpp>
 #include <stan/math/fwd/functor/integrate_1d.hpp>
+#include <stan/math/fwd/functor/integrate_1d_double_exponential.hpp>
+#include <stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp>
 #include <stan/math/fwd/functor/jacobian.hpp>
 #include <stan/math/fwd/functor/operands_and_partials.hpp>
 #include <stan/math/fwd/functor/partials_propagator.hpp>

stan/math/fwd/functor/integrate_1d_double_exponential.hpp

@@ -0,0 +1,97 @@
+#ifndef STAN_MATH_FWD_FUNCTOR_INTEGRATE_1D_DOUBLE_EXPONENTIAL_HPP
+#define STAN_MATH_FWD_FUNCTOR_INTEGRATE_1D_DOUBLE_EXPONENTIAL_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/functor/integrate_1d_double_exponential.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/meta/forward_as.hpp>
+#include <stan/math/prim/functor/apply.hpp>
+#include <stan/math/fwd/functor/finite_diff.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the integral of f from a to b using adaptive double-exponential
+ * quadrature, with tangents computed via finite differences over the
+ * integrand parameters.
+ *
+ * @tparam F Type of f
+ * @tparam T_a type of first limit
+ * @tparam T_b type of second limit
+ * @tparam Args types of parameter pack arguments
+ *
+ * @param f the functor to integrate
+ * @param a lower limit of integration
+ * @param b upper limit of integration
+ * @param relative_tolerance relative tolerance passed to Boost quadrature
+ * @param absolute_tolerance absolute-error floor on the convergence test
+ * @param max_refinements maximum refinement level passed to the Boost
+ *   quadrature class constructor
+ * @param[in, out] msgs the print stream for warning messages
+ * @param args additional arguments to pass to f
+ * @return numeric integral of function f
+ */
+template <typename F, typename T_a, typename T_b, typename... Args,
+          require_any_st_fvar<T_a, T_b, Args...> * = nullptr>
+inline return_type_t<T_a, T_b, Args...> integrate_1d_double_exponential_tol(
+    const F &f, const T_a &a, const T_b &b, double relative_tolerance,
+    double absolute_tolerance, int max_refinements, std::ostream *msgs,
+    const Args &... args) {
+  using FvarT = scalar_type_t<return_type_t<T_a, T_b, Args...>>;
+
+  auto a_val = value_of(a);
+  auto b_val = value_of(b);
+  auto func = [f, msgs, relative_tolerance, absolute_tolerance, max_refinements,
+               a_val, b_val](const auto &... args_var) {
+    return integrate_1d_double_exponential_tol(
+        f, a_val, b_val, relative_tolerance, absolute_tolerance,
+        max_refinements, msgs, args_var...);
+  };
+  FvarT ret = finite_diff(func, args...);
+  if constexpr (is_fvar<T_a>::value || is_fvar<T_b>::value) {
+    if constexpr (is_fvar<T_a>::value) {
+      ret.d_ += a.d_ * -f(a_val, 0.0, msgs, value_of(args)...);
+    }
+    if constexpr (is_fvar<T_b>::value) {
+      ret.d_ += b.d_ * f(b_val, 0.0, msgs, value_of(args)...);
+    }
+  }
+  return ret;
+}
+
+/**
+ * Compute the integral of the single variable function f from a to b using
+ * adaptive double-exponential quadrature. a and b can be finite or
+ * infinite.
+ *
+ * @tparam F Type of f
+ * @tparam T_a type of first limit
+ * @tparam T_b type of second limit
+ * @tparam Args types of parameter pack arguments
+ *
+ * @param f the functor to integrate
+ * @param a lower limit of integration
+ * @param b upper limit of integration
+ * @param relative_tolerance relative tolerance passed to Boost quadrature
+ * @param absolute_tolerance absolute-error floor on the convergence test
+ * @param max_refinements maximum refinement level passed to the Boost
+ *   quadrature class constructor
+ * @param[in, out] msgs the print stream for warning messages
+ * @param args additional arguments to pass to f
+ * @return numeric integral of function f
+ */
+template <typename F, typename T_a, typename T_b, typename... Args,
+          require_any_st_fvar<T_a, T_b, Args...> * = nullptr>
+inline return_type_t<T_a, T_b, Args...> integrate_1d_double_exponential_impl(
+    const F &f, const T_a &a, const T_b &b, double relative_tolerance,
+    double absolute_tolerance, int max_refinements, 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

stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp

@@ -0,0 +1,95 @@
+#ifndef STAN_MATH_FWD_FUNCTOR_INTEGRATE_1D_GAUSS_KRONROD_HPP
+#define STAN_MATH_FWD_FUNCTOR_INTEGRATE_1D_GAUSS_KRONROD_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/meta/forward_as.hpp>
+#include <stan/math/prim/functor/apply.hpp>
+#include <stan/math/fwd/functor/finite_diff.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the integral of f from a to b using adaptive Gauss-Kronrod (G21,K21)
+ * quadrature, with tangents computed via finite differences over the
+ * integrand parameters.
+ *
+ * @tparam F Type of f
+ * @tparam T_a type of first limit
+ * @tparam T_b type of second limit
+ * @tparam Args types of parameter pack arguments
+ *
+ * @param f the functor to integrate
+ * @param a lower limit of integration
+ * @param b upper limit of integration
+ * @param relative_tolerance relative tolerance passed to Boost quadrature
+ * @param absolute_tolerance absolute-error floor on the convergence test
+ * @param max_depth maximum recursive bisection depth passed to Boost
+ *   quadrature
+ * @param[in, out] msgs the print stream for warning messages
+ * @param args additional arguments to pass to f
+ * @return numeric integral of function f
+ */
+template <typename F, typename T_a, typename T_b, typename... Args,
+          require_any_st_fvar<T_a, T_b, Args...> * = nullptr>
+inline return_type_t<T_a, T_b, Args...> integrate_1d_gauss_kronrod_tol(
+    const F &f, const T_a &a, const T_b &b, double relative_tolerance,
+    double absolute_tolerance, int max_depth, std::ostream *msgs,
+    const Args &... args) {
+  using FvarT = scalar_type_t<return_type_t<T_a, T_b, Args...>>;
+
+  // Wrap integrate_1d_gauss_kronrod call in a functor where the input
+  // arguments are only those for which tangents are needed.
+  auto a_val = value_of(a);
+  auto b_val = value_of(b);
+  auto func = [f, msgs, relative_tolerance, absolute_tolerance, max_depth,
+               a_val, b_val](const auto &... args_var) {
+    return integrate_1d_gauss_kronrod_tol(f, a_val, b_val, relative_tolerance,
+                                          absolute_tolerance, max_depth, msgs,
+                                          args_var...);
+  };
+  FvarT ret = finite_diff(func, args...);
+  // Calculate tangents w.r.t. integration bounds if needed
+  if constexpr (is_fvar<T_a>::value || is_fvar<T_b>::value) {
+    if constexpr (is_fvar<T_a>::value) {
+      ret.d_ += a.d_ * -f(a_val, 0.0, msgs, value_of(args)...);
+    }
+    if constexpr (is_fvar<T_b>::value) {
+      ret.d_ += b.d_ * f(b_val, 0.0, msgs, value_of(args)...);
+    }
+  }
+  return ret;
+}
+
+/**
+ * Return the integral of f from a to b using adaptive Gauss-Kronrod (G21,K21)
+ * quadrature, with tangents computed via finite differences over the
+ * integrand parameters.
+ *
+ * @tparam F Type of f
+ * @tparam T_a type of first limit
+ * @tparam T_b type of second limit
+ * @tparam Args types of parameter pack arguments
+ *
+ * @param f the functor to integrate
+ * @param a lower limit of integration
+ * @param b upper limit of integration
+ * @param[in, out] msgs the print stream for warning messages
+ * @param args additional arguments to pass to f
+ * @return numeric integral of function f
+ */
+template <typename F, typename T_a, typename T_b, typename... Args,
+          require_any_st_fvar<T_a, T_b, Args...> * = nullptr>
+inline return_type_t<T_a, T_b, Args...> integrate_1d_gauss_kronrod(
+    const F &f, const T_a &a, const T_b &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

stan/math/mix/functor/laplace_likelihood.hpp

@@ -2,7 +2,7 @@
 #define STAN_MATH_MIX_FUNCTOR_LAPLACE_LIKELIHOOD_HPP
 
 #include <stan/math/mix/functor/hessian_block_diag.hpp>
-#include <stan/math/mix/functor/conditional_copy_and_promote.hpp>
+#include <stan/math/rev/functor/conditional_copy_and_promote.hpp>
 #include <stan/math/prim/functor.hpp>
 #include <stan/math/prim/fun.hpp>
 

stan/math/mix/functor/laplace_marginal_density.hpp

@@ -3,7 +3,7 @@
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/mix/functor/laplace_likelihood.hpp>
 #include <stan/math/mix/functor/laplace_marginal_density_estimator.hpp>
-#include <stan/math/mix/functor/conditional_copy_and_promote.hpp>
+#include <stan/math/rev/functor/conditional_copy_and_promote.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun.hpp>
@@ -135,40 +135,6 @@ inline constexpr void copy_compute_s2(Output&& output, Input&& input) {
       std::forward<Output>(output), std::forward<Input>(input));
 }
 
-/**
- * Collects adjoints from a tuple or std::vector of tuples
- * @tparam Output A tuple or std::vector of tuples where all scalar types are
- * `var` types
- * @tparam Input A tuple or std::vector of tuples where all scalar types are
- * `arithmetic` types
- * @param ret The vari object containing the adjoint to be added
- * @param output The output to which the adjoints will be added
- * @param input The input from which the adjoints will be collected
- */
-template <typename Output, typename Input>
-inline void reverse_pass_collect_adjoints(var ret, Output&& output,
-                                          Input&& input) {
-  if constexpr (is_tuple_v<Output>) {
-    stan::math::for_each(
-        [ret](auto&& inner_arg, auto&& inner_input) mutable {
-          reverse_pass_collect_adjoints(
-              ret, std::forward<decltype(inner_arg)>(inner_arg),
-              std::forward<decltype(inner_input)>(inner_input));
-        },
-        std::forward<Output>(output), std::forward<Input>(input));
-  } else if constexpr (is_std_vector_containing_tuple_v<Output>) {
-    for (std::size_t i = 0; i < output.size(); ++i) {
-      reverse_pass_collect_adjoints(ret, output[i], input[i]);
-    }
-  } else {
-    reverse_pass_callback(
-        [vi = ret.vi_, arg_arena = to_arena(std::forward<Output>(output)),
-         input_arena = to_arena(std::forward<Input>(input))]() mutable {
-          collect_adjoints(arg_arena, vi, input_arena);
-        });
-  }
-}
-
 }  // namespace internal
 
 /**

stan/math/prim/functor.hpp

@@ -14,6 +14,8 @@
 #include <stan/math/prim/functor/hcubature.hpp>
 #include <stan/math/prim/functor/integrate_1d.hpp>
 #include <stan/math/prim/functor/integrate_1d_adapter.hpp>
+#include <stan/math/prim/functor/integrate_1d_double_exponential.hpp>
+#include <stan/math/prim/functor/integrate_1d_gauss_kronrod.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/iter_tuple_nested.hpp>