diff --git a/stan/math/prim/functor.hpp b/stan/math/prim/functor.hpp index 0591d4b65dc..0ec5c343ff7 100644 --- a/stan/math/prim/functor.hpp +++ b/stan/math/prim/functor.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index 8e76694dd4b..e4494c81ed8 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/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 #include #include +#include #include #include #include @@ -50,17 +51,53 @@ namespace math { template 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 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 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 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 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 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 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 * = 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 -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& theta, const std::vector& x_r, const std::vector& 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(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), a, b, relative_tolerance, + msgs, theta, x_r, x_i); } } // namespace math diff --git a/stan/math/prim/functor/integrate_1d_adapter.hpp b/stan/math/prim/functor/integrate_1d_adapter.hpp new file mode 100644 index 00000000000..ecfcaaaa9a6 --- /dev/null +++ b/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 +#include + +/** + * Adapt the non-variadic integrate_1d arguments to the variadic + * integrate_1d_impl interface + * + * @tparam F type of function to adapt + */ +template +struct integrate_1d_adapter { + const F& f_; + + explicit integrate_1d_adapter(const F& f) : f_(f) {} + + template + auto operator()(const T_a& x, const T_b& xc, std::ostream* msgs, + const std::vector& theta, + const std::vector& x_r, + const std::vector& x_i) const { + return f_(x, xc, theta, x_r, x_i, msgs); + } +}; + +#endif diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index a02e36b57ad..be8ff106393 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -20,43 +20,139 @@ namespace stan { namespace math { /** - * Calculate first derivative of f(x, param, std::ostream&) - * with respect to the nth parameter. Uses nested reverse mode autodiff + * Return the integral of f from a to b to the given relative tolerance * - * Gradients that evaluate to NaN are set to zero if the function itself - * evaluates to zero. If the function is not zero and the gradient evaluates to - * NaN, a std::domain_error is thrown + * @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 * - * @tparam F type of f + * @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[in, out] msgs the print stream for warning messages + * @param args additional arguments to pass to f + * @return numeric integral of function f */ -template -inline double gradient_of_f(const F &f, const double &x, const double &xc, - const std::vector &theta_vals, - const std::vector &x_r, - const std::vector &x_i, size_t n, - std::ostream *msgs) { - double gradient = 0.0; - - // Run nested autodiff in this scope - nested_rev_autodiff nested; - - std::vector theta_var(theta_vals.size()); - for (size_t i = 0; i < theta_vals.size(); i++) { - theta_var[i] = theta_vals[i]; - } - var fx = f(x, xc, theta_var, x_r, x_i, msgs); - fx.grad(); - gradient = theta_var[n].adj(); - if (is_nan(gradient)) { - if (fx.val() == 0) { - gradient = 0; - } else { - throw_domain_error("gradient_of_f", "The gradient of f", n, - "is nan for parameter ", ""); +template * = nullptr> +inline return_type_t integrate_1d_impl( + const F &f, const T_a &a, const T_b &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); + + double a_val = value_of(a); + double b_val = value_of(b); + + if (unlikely(a_val == b_val)) { + if (is_inf(a_val)) { + throw_domain_error(function, "Integration endpoints are both", a_val, "", + ""); + } + return var(0.0); + } else { + auto args_val_tuple = std::make_tuple(value_of(args)...); + + double integral = integrate( + [&](const auto &x, const auto &xc) { + return apply( + [&](auto &&... val_args) { return f(x, xc, msgs, val_args...); }, + args_val_tuple); + }, + a_val, b_val, relative_tolerance); + + constexpr size_t num_vars_ab = is_var::value + is_var::value; + size_t num_vars_args = count_vars(args...); + vari **varis = ChainableStack::instance_->memalloc_.alloc_array( + num_vars_ab + num_vars_args); + double *partials = ChainableStack::instance_->memalloc_.alloc_array( + num_vars_ab + num_vars_args); + // We move this pointer up based on whether we a or b is a var type. + double *partials_ptr = partials; + + save_varis(varis, a, b, args...); + + for (size_t i = 0; i < num_vars_ab + num_vars_args; ++i) { + partials[i] = 0.0; + } + + if (is_var::value && !is_inf(a)) { + *partials_ptr = apply( + [&f, a_val, msgs](auto &&... val_args) { + return -f(a_val, 0.0, msgs, val_args...); + }, + args_val_tuple); + partials_ptr++; + } + + if (!is_inf(b) && is_var::value) { + *partials_ptr = apply( + [&f, b_val, msgs](auto &&... val_args) { + return f(b_val, 0.0, msgs, val_args...); + }, + args_val_tuple); + partials_ptr++; + } + + { + nested_rev_autodiff argument_nest; + // The arguments copy is used multiple times in the following nests, so + // do it once in a separate nest for efficiency + auto args_tuple_local_copy = std::make_tuple(deep_copy_vars(args)...); + + // Save the varis so it's easy to efficiently access the nth adjoint + std::vector local_varis(num_vars_args); + apply( + [&](const auto &... args) { + save_varis(local_varis.data(), args...); + }, + args_tuple_local_copy); + + for (size_t n = 0; n < num_vars_args; ++n) { + // This computes the integral of the gradient of f with respect to the + // nth parameter in args using a nested nested reverse mode autodiff + *partials_ptr = integrate( + [&](const auto &x, const auto &xc) { + argument_nest.set_zero_all_adjoints(); + + nested_rev_autodiff gradient_nest; + var fx = apply( + [&f, &x, &xc, msgs](auto &&... local_args) { + return f(x, xc, msgs, local_args...); + }, + args_tuple_local_copy); + fx.grad(); + + double gradient = local_varis[n]->adj(); + + // Gradients that evaluate to NaN are set to zero if the function + // itself evaluates to zero. If the function is not zero and the + // gradient evaluates to NaN, a std::domain_error is thrown + if (is_nan(gradient)) { + if (fx.val() == 0) { + gradient = 0; + } else { + throw_domain_error("gradient_of_f", "The gradient of f", n, + "is nan for parameter ", ""); + } + } + return gradient; + }, + a_val, b_val, relative_tolerance); + partials_ptr++; + } } - } - return gradient; + return make_callback_var( + integral, + [total_vars = num_vars_ab + num_vars_args, varis, partials](auto &vi) { + for (size_t i = 0; i < total_vars; ++i) { + varis[i]->adj_ += partials[i] * vi.adj(); + } + }); + } } /** @@ -120,52 +216,8 @@ inline return_type_t integrate_1d( const F &f, const T_a &a, const T_b &b, const std::vector &theta, const std::vector &x_r, const std::vector &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 (value_of(a) == value_of(b)) { - if (is_inf(a)) { - throw_domain_error(function, "Integration endpoints are both", - value_of(a), "", ""); - } - return var(0.0); - } else { - double integral = integrate( - std::bind(f, std::placeholders::_1, std::placeholders::_2, - value_of(theta), x_r, x_i, msgs), - value_of(a), value_of(b), relative_tolerance); - - size_t N_theta_vars = is_var::value ? theta.size() : 0; - std::vector dintegral_dtheta(N_theta_vars); - std::vector theta_concat(N_theta_vars); - - if (N_theta_vars > 0) { - std::vector theta_vals = value_of(theta); - - for (size_t n = 0; n < N_theta_vars; ++n) { - dintegral_dtheta[n] = integrate( - std::bind(gradient_of_f, f, std::placeholders::_1, - std::placeholders::_2, theta_vals, x_r, x_i, n, - msgs), - value_of(a), value_of(b), relative_tolerance); - theta_concat[n] = theta[n]; - } - } - - if (!is_inf(a) && is_var::value) { - theta_concat.push_back(a); - dintegral_dtheta.push_back( - -value_of(f(value_of(a), 0.0, theta, x_r, x_i, msgs))); - } - - if (!is_inf(b) && is_var::value) { - theta_concat.push_back(b); - dintegral_dtheta.push_back( - value_of(f(value_of(b), 0.0, theta, x_r, x_i, msgs))); - } - - return precomputed_gradients(integral, theta_concat, dintegral_dtheta); - } + return integrate_1d_impl(integrate_1d_adapter(f), a, b, relative_tolerance, + msgs, theta, x_r, x_i); } } // namespace math diff --git a/test/unit/math/prim/functor/integrate_1d_impl_test.cpp b/test/unit/math/prim/functor/integrate_1d_impl_test.cpp new file mode 100644 index 00000000000..985a4cd6392 --- /dev/null +++ b/test/unit/math/prim/functor/integrate_1d_impl_test.cpp @@ -0,0 +1,431 @@ +#include +#include +#include +#include +#include +#include +#include + +namespace integrate_1d_impl_test { + +std::ostringstream *msgs = nullptr; + +struct f1 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { + return exp(-x) / sqrt(x); + } +}; + +struct f2 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { + if (x <= 0.5) { + return sqrt(x) / sqrt(1 - x * x); + } else { + return sqrt(x / ((x + 1) * (xc))); + } + } +}; + +struct f3 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { + return exp(-x); + } +}; + +struct f4 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs, + const std::vector &theta) const { + return exp(x) + theta[0]; + } +}; + +struct f5 { + template + inline stan::return_type_t operator()( + const T1 &x, const T1 &xc, std::ostream *msgs, + const std::vector &theta) const { + return exp(x) + pow(theta[0], 2) + pow(theta[1], 3); + } +}; + +struct f6 { + template + inline stan::return_type_t operator()( + const T1 &x, const T1 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_i) const { + return exp(x) + pow(x_i[0], 2) + pow(theta[0], 4) + 3 * theta[1]; + } +}; + +struct f7 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs, + const std::vector &x_r) const { + return exp(x) + pow(x_r[0], 2) + pow(x_r[1], 5) + 3 * x_r[2]; + } +}; + +struct f8 { + template + inline stan::return_type_t operator()( + const T1 &x, const T1 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return exp(-pow(x - theta[0], x_i[0]) / pow(x_r[0], x_i[0])); + } +}; + +struct f9 { + template + inline stan::return_type_t operator()( + const T1 &x, const T1 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_i) const { + return 1.0 / (1.0 + pow(x, x_i[0]) / theta[0]); + } +}; + +struct f10 { + template + inline stan::return_type_t operator()( + const T1 &x, const T1 &xc, std::ostream *msgs, + const std::vector &theta) const { + return pow(x, theta[0] - 1.0) + * pow((x > 0.5) ? xc : (1 - x), theta[1] - 1.0); + } +}; + +struct f11 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { + return (std::isnan(xc)) ? xc : 0.0; + } +}; + +struct f12 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { + T1 out = stan::math::modified_bessel_second_kind(0, x); + if (out > 0) + return 2 * x * out; + return out; + } +}; + +struct f13 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { + T1 out = stan::math::modified_bessel_second_kind(0, x); + if (out > 0) + return 2 * x * stan::math::square(out); + return out; + } +}; + +struct f14 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { + return exp(x) * stan::math::inv_sqrt(x > 0.5 ? xc : 1 - x); + } +}; + +struct f15 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { + T1 x2 = x * x; + T1 numer = x2 * log(x); + T1 denom = x < 0.5 ? (x + 1) * (x - 1) : (x + 1) * -xc; + denom *= x2 * x2 + 1; + return numer / denom; + } +}; + +struct f16 { + template + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { + return x * sin(x) / (1 + stan::math::square(cos(x))); + } +}; + +double lbaX_pdf(double X, double t, double A, double v, double s, + std::ostream *pstream__) { + double b_A_tv_ts; + double b_tv_ts; + double term_1; + double term_2; + double pdf; + + b_A_tv_ts = (((X - A) - (t * v)) / (t * s)); + b_tv_ts = ((X - (t * v)) / (t * s)); + term_1 = stan::math::Phi(b_A_tv_ts); + term_2 = stan::math::Phi(b_tv_ts); + pdf = ((1 / A) * (-term_1 + term_2)); + return pdf; +} + +double lbaX_cdf(double X, double t, double A, double v, double s, + std::ostream *pstream__) { + double b_A_tv; + double b_tv; + double ts; + double term_1; + double term_2; + double term_3; + double term_4; + double cdf; + + b_A_tv = ((X - A) - (t * v)); + b_tv = (X - (t * v)); + ts = (t * s); + term_1 = (b_A_tv * stan::math::Phi((b_A_tv / ts))); + term_2 = (b_tv * stan::math::Phi((b_tv / ts))); + term_3 = (ts + * stan::math::exp( + stan::math::normal_lpdf((b_A_tv / ts), 0, 1))); + term_4 + = (ts + * stan::math::exp(stan::math::normal_lpdf((b_tv / ts), 0, 1))); + cdf = ((1 / A) * (((-term_1 + term_2) - term_3) + term_4)); + return cdf; +} + +double rank_density(double x, double xc, const std::vector &theta, + const std::vector &x_r, const std::vector &x_i, + std::ostream *pstream__) { + double t = theta[0]; + double A = theta[1]; + double v1 = theta[2]; + double v2 = theta[3]; + double s = theta[4]; + double v = (lbaX_pdf(x, t, A, v1, s, pstream__) + * lbaX_cdf(x, t, A, v2, s, pstream__)); + return v; +} + +struct rank_density_functor__ { + double operator()(double x, double xc, std::ostream *pstream__, + const std::vector &theta, + const std::vector &x_r, + const std::vector &x_i) const { + return rank_density(x, xc, theta, x_r, x_i, pstream__); + } +}; + +double order(double down, double up, const std::vector &theta, + const std::vector &x_r, std::ostream *pstream__) { + std::vector x_i; + + double v; + + v = stan::math::integrate_1d_impl(rank_density_functor__(), down, up, 1e-8, + pstream__, theta, x_r, x_i); + return v; +} +} // namespace integrate_1d_impl_test +/* + * test_integration is a helper function to make it easy to test the + * integrate_1d_impl function. + * + * It takes in a callable function object, integration limits, parameters, real + * and integer data. It integrates the provided function and compares the + * computed integral against the provided integral (val). + * + * The prototype for f is: + * struct f10 { + * inline double operator()(const double& x, const double& xc, const + * std::vector& theta, const std::vector& x_r, const + * std::vector& x_i, std::ostream& msgs) const { + * } + * }; + * + * @tparam F Type of f + * @param f a functor with signature given above + * @param a lower limit of integration + * @param b upper limit of integration + * @param param parameters to be passed to f (should be + * std::vector) + * @param x_r real data to be passed to f (should be std::vector) + * @param x_i integer data to be passed to f (should be std::vector) + * @param val correct value of integral + */ +template +void test_integration(const F &f, double a, double b, double val, + const Args &... args) { + using stan::math::integrate_1d_impl; + + std::vector tolerances = {1e-4, 1e-6, 1e-8}; + + for (auto tolerance : tolerances) { + EXPECT_LE(std::abs(integrate_1d_impl(f, a, b, tolerance, + integrate_1d_impl_test::msgs, args...) + - val), + tolerance); + // Flip the domain of integration and check that the integral is working + auto flipped + = [&](double x, double xc, std::ostream *msgs, const auto &... args) { + return f(-x, -xc, msgs, args...); + }; + + EXPECT_LE(std::abs(integrate_1d_impl(flipped, -b, -a, tolerance, + integrate_1d_impl_test::msgs, args...) + - val), + tolerance); + } +} + +TEST(StanMath_integrate_1d_impl_prim, TestThrows) { + // Left limit of integration must be less than or equal to right limit + EXPECT_THROW( + stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, 1.0, 0.0, 0.0, + integrate_1d_impl_test::msgs), + std::domain_error); + // NaN limits not okay + EXPECT_THROW( + stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, 0.0, + std::numeric_limits::quiet_NaN(), + 0.0, integrate_1d_impl_test::msgs), + std::domain_error); + EXPECT_THROW( + stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, + std::numeric_limits::quiet_NaN(), + 0.0, 0.0, integrate_1d_impl_test::msgs), + std::domain_error); + EXPECT_THROW( + stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, + std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN(), + 0.0, integrate_1d_impl_test::msgs), + std::domain_error); + // Two of the same inf limits not okay + EXPECT_THROW( + stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, + -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + 0.0, integrate_1d_impl_test::msgs), + std::domain_error); + + EXPECT_THROW( + stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + 0.0, integrate_1d_impl_test::msgs), + std::domain_error); + // xc should be nan if there are infinite limits + EXPECT_THROW( + stan::math::integrate_1d_impl(integrate_1d_impl_test::f11{}, 0.0, + std::numeric_limits::infinity(), + 0.0, integrate_1d_impl_test::msgs), + std::runtime_error); + EXPECT_THROW( + stan::math::integrate_1d_impl(integrate_1d_impl_test::f11{}, + std::numeric_limits::infinity(), + 0.0, 0.0, integrate_1d_impl_test::msgs), + std::domain_error); + EXPECT_THROW( + stan::math::integrate_1d_impl(integrate_1d_impl_test::f11{}, + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + 0.0, integrate_1d_impl_test::msgs), + std::domain_error); + // But not otherwise + EXPECT_NO_THROW(stan::math::integrate_1d_impl(integrate_1d_impl_test::f11{}, + 0.0, 1.0, 0.0, + integrate_1d_impl_test::msgs)); +} + +TEST(StanMath_integrate_1d_impl_prim, test_integer_arguments) { + double v; + EXPECT_NO_THROW( + v = stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, 0, 1, 0.0, + integrate_1d_impl_test::msgs)); + EXPECT_NO_THROW( + v = stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, 0.0, 1, + 0.0, integrate_1d_impl_test::msgs)); + EXPECT_NO_THROW( + v = stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, 0, 1.0, + 0.0, integrate_1d_impl_test::msgs)); +} + +TEST(StanMath_integrate_1d_impl_prim, test1) { + // Tricky integral from Boost docs + limit at infinity + test_integration(integrate_1d_impl_test::f1{}, 0.0, + std::numeric_limits::infinity(), 1.772453850905516); + // Tricky integral from Boost 1d integration docs + test_integration(integrate_1d_impl_test::f2{}, 0.0, 1.0, 1.198140234735592); + // Tricky integral from Boost 1d integration docs + test_integration(integrate_1d_impl_test::f2{}, 0, 1, 1.198140234735592); + // Zero crossing integral + limit at infinity + test_integration(integrate_1d_impl_test::f3{}, -2.0, + std::numeric_limits::infinity(), 7.38905609893065); + // Easy integrals + test_integration(integrate_1d_impl_test::f4{}, 0.2, 0.7, 1.0423499493102901, + std::vector({0.5})); + test_integration(integrate_1d_impl_test::f5{}, -0.2, 0.7, 1.396621954392482, + std::vector({0.4, 0.4})); + test_integration(integrate_1d_impl_test::f4{}, 0.0, 0.0, 0.0, + std::vector({0.5})); + test_integration(integrate_1d_impl_test::f5{}, 1.0, 1.0, 0.0, + std::vector({0.4, 0.4})); + // Test x_i + test_integration(integrate_1d_impl_test::f6{}, -0.2, 2.9, 4131.985414616364, + std::vector({6.0, 5.1}), std::vector({4})); + // Test x_r + test_integration(integrate_1d_impl_test::f7{}, -0.2, 2.9, 24219.985414616367, + std::vector({4.0, 6.0, 5.1})); + // Both limits at infinity + test x_r/x_i + test_integration(integrate_1d_impl_test::f8{}, + -std::numeric_limits::infinity(), + std::numeric_limits::infinity(), 3.013171546539377, + std::vector({5.0}), std::vector({1.7}), + std::vector({2})); + // Both limits at infinity + test x_i + test_integration(integrate_1d_impl_test::f9{}, + -std::numeric_limits::infinity(), + std::numeric_limits::infinity(), 2.372032924895055, + std::vector({1.3}), std::vector({4})); + // Various integrals of beta function + test_integration(integrate_1d_impl_test::f10{}, 0.0, 1.0, 19.71463948905016, + std::vector({0.1, 0.1})); + test_integration(integrate_1d_impl_test::f10{}, 0.0, 1.0, 11.32308697521577, + std::vector({0.1, 0.5})); + test_integration(integrate_1d_impl_test::f10{}, 0.0, 1.0, 11.32308697521577, + std::vector({0.5, 0.1})); + test_integration(integrate_1d_impl_test::f10{}, 0.0, 1.0, 0.00952380952380952, + std::vector({5.0, 3.0})); + + // Integrals from + // http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/dhb-tanh-sinh.pdf + test_integration(integrate_1d_impl_test::f12{}, 0.0, + std::numeric_limits::infinity(), 2.0); + test_integration(integrate_1d_impl_test::f13{}, 0.0, + std::numeric_limits::infinity(), 1.0); + test_integration(integrate_1d_impl_test::f14{}, 0.0, 1.0, + exp(1) * sqrt(stan::math::pi()) * stan::math::erf(1.0)); + + // Integrals from http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf + // works normally but not to tolerance when limits of integration are flipped + // test_integration(f15{}, 0.0, 1.0, {}, {}, {}, + // stan::math::square(stan::math::pi()) * (2 - sqrt(2.0)) / + // 32); + test_integration(integrate_1d_impl_test::f16{}, 0.0, stan::math::pi(), + stan::math::square(stan::math::pi()) / 4); +} + +TEST(StanMath_integrate_1d_impl_prim, TestTolerance) { + std::ostringstream *msgs = nullptr; + + double t = 0.5; + double b = 1.0; + double A = 0.5; + double v1 = 1.0; + double v2 = 1.0; + double s = 1.0; + + std::vector theta = {t, A, v1, v2, s}; + std::vector x_r; + + EXPECT_NO_THROW(integrate_1d_impl_test::order(-10, 0.67, theta, x_r, msgs)); +} diff --git a/test/unit/math/rev/functor/integrate_1d_impl_test.cpp b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp new file mode 100644 index 00000000000..8244a2638ad --- /dev/null +++ b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp @@ -0,0 +1,939 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +std::ostringstream *msgs = nullptr; + +struct f1 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return exp(x) + theta[0]; + } +}; + +struct f2 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return exp(theta[0] * cos(2 * 3.141593 * x)) + theta[0]; + } +}; + +struct f3 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return exp(x) + pow(theta[0], x_r[0]) + 2 * pow(theta[1], x_r[1]) + + 2 * theta[2]; + } +}; + +struct f4 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return exp(-x) / sqrt(x); + } +}; + +struct f5 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return exp(-theta[0] * x) / sqrt(theta[1] * x); + } +}; + +struct f6 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return sqrt(x / (1 - theta[0] * x * x)); + } +}; + +struct f7 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return exp(-theta[0] * x); + } +}; + +struct f8 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return exp(theta[0] * x); + } +}; + +struct f10 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return 1 / (1 + pow(x, x_i[0]) / x_r[0]); + } +}; + +struct f11 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return pow(x, theta[0] - 1.0) + * pow((x > 0.5) ? xc : (1 - x), theta[1] - 1.0); + } +}; + +struct f12 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + T3 mu = theta[0]; + T3 sigma = theta[1]; + return exp(-0.5 * stan::math::square((x - mu) / sigma)) + / (sigma * sqrt(2.0 * stan::math::pi())); + } +}; + +struct f13 { + template + inline stan::return_type_t operator()( + const T1 &x, const T2 &xc, std::ostream *msgs, + const std::vector &theta, const std::vector &x_r, + const std::vector &x_i) const { + return x + theta[0] + theta[1]; + } +}; + +/* + * Return the adjoint if the argument is a var + * + * @param v variable + * @return adjoint of var + */ +double get_adjoint_if_var(stan::math::var v) { return v.adj(); } + +/* + * If the argument is not a var, return a NaN + * + * @param v variable + * @return NaN + */ +double get_adjoint_if_var(double v) { + return std::numeric_limits::quiet_NaN(); +} + +/* + * test_derivatives is a helper function to make it easy to test the + * integrate_1d function. + * + * It takes in a callable function object, integration limits, parameters, real + * and integer data. It integrates the provided function and compares the + * computed integral and gradients against the provided integral value (val) and + * gradients (grad, d_a, d_b). + * + * The first three template arguments specify how the left limit, right limit, + * and parameters are tested. If the template arguments are vars, then the + * gradients are checked against the provided gradients, otherwise the provided + * gradients aren't used. + * + * The prototype for f is: + * struct f10 { + * template + * inline stan::return_type_t operator()( + * const T1& x, const T2& xc, const std::vector& theta, const + * std::vector& x_r, const std::vector& x_i, std::ostream& msgs) + * const { + * } + * }; + * + * @tparam T_a Type of integration left limit + * @tparam T_b Type of integration right limit + * @tparam T_theta Type of parameters + * @tparam F Type of f + * @param f a functor with signature given above + * @param a lower limit of integration + * @param b upper limit of integration + * @param param parameters to be passed to f (should be + * std::vector) + * @param x_r real data to be passed to f (should be std::vector) + * @param x_i integer data to be passed to f (should be std::vector) + * @param val correct value of integral + * @param grad correct value of gradients (not used if T_theta is not var) + * @param d_a correct value of gradient of integral with respect to left hand + * limit (not used if T_a is not var) + * @param d_b correct value of gradient of integral with respect to right hand + * limit (not used if T_b is not var) + */ +template +void test_derivatives(const F &f, double a, double b, + std::vector thetas, + const std::vector &x_r, + const std::vector &x_i, double val, + std::vector grad, double d_a = 0.0, + double d_b = 0.0) { + using stan::math::value_of; + using stan::math::var; + + std::vector tolerances = {1e-4, 1e-6, 1e-8}; + + for (auto tolerance : tolerances) { + // Reset autodiff stack + stan::math::recover_memory(); + // Convert endpoints into given template types + T_a a_(a); + T_b b_(b); + std::vector thetas_(thetas.size()); + for (size_t i = 0; i < thetas.size(); ++i) + thetas_[i] = thetas[i]; + + var integral = stan::math::integrate_1d_impl(f, a_, b_, tolerance, msgs, + thetas_, x_r, x_i); + integral.grad(); + EXPECT_LE(std::abs(val - integral.val()), tolerance); + if (stan::is_var::value) { + for (size_t i = 0; i < grad.size(); ++i) + EXPECT_LE(std::abs(grad[i] - get_adjoint_if_var(thetas_[i])), + tolerance); + } + if (stan::is_var::value) { + EXPECT_LE(std::abs(d_a - get_adjoint_if_var(a_)), tolerance); + } + if (stan::is_var::value) { + EXPECT_LE(std::abs(d_b - get_adjoint_if_var(b_)), tolerance); + } + } +} + +TEST(StanMath_integrate_1d_impl_rev, test_integer_arguments) { + stan::math::var v; + std::vector theta = {0.5}; + std::vector x_r; + std::vector x_i; + EXPECT_NO_THROW(v = stan::math::integrate_1d_impl(f2{}, 0, 1, 1e-6, msgs, + theta, x_r, x_i)); + EXPECT_NO_THROW(v = stan::math::integrate_1d_impl(f2{}, 0.0, 1, 1e-6, msgs, + theta, x_r, x_i)); + EXPECT_NO_THROW(v = stan::math::integrate_1d_impl(f2{}, 0, 1.0, 1e-6, msgs, + theta, x_r, x_i)); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_easy) { + // Easy integrals + using stan::math::var; + test_derivatives(f1{}, 0.2, 0.7, {0.75}, {}, {}, + 0.7923499493102901 + 0.5 * 0.75, {0.5}); + test_derivatives(f2{}, 0.0, 1.0, {0.5}, {}, {}, + 1.56348343527304, {1.25789445875152}, + -2.148721270700128, 2.14872127069993); + test_derivatives(f2{}, 0.0, 1.0, {0.5}, {}, {}, + 1.56348343527304, {}, -2.148721270700128, + 2.14872127069993); + test_derivatives(f1{}, 0.0, 0.0, {0.75}, {}, {}, 0.0, + {0.0}); + test_derivatives(f2{}, 1.0, 1.0, {0.5}, {}, {}, 0.0, + {0.0}); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_zero_crossing) { + // Zero crossing integral + test x_r + vars at endpoints + using stan::math::var; + test_derivatives(f3{}, -1.0, 1.0, {0.5, 1.75, 3.9}, {2.5, 3.0}, + {}, + 2.350402387287579 + 2.0 * pow(0.5, 2.5) + + 4.0 * pow(1.75, 3.0) + 4.0 * 3.9, + {5 * pow(0.5, 1.5), 12 * 1.75 * 1.75, 4.0}, + -19.06340613646808, 21.41380852375568); +} + +TEST(StanMath_integrate_1d_impl_rev, + TestDerivatives_var_right_endpoint_var_params) { + // Zero crossing integral + test x_r + vars at right endpoint + using stan::math::var; + test_derivatives( + f3{}, -1.0, 1.0, {0.5, 1.75, 3.9}, {2.5, 3.0}, {}, + 2.350402387287579 + 2.0 * pow(0.5, 2.5) + 4.0 * pow(1.75, 3.0) + + 4.0 * 3.9, + {5 * pow(0.5, 1.5), 12 * 1.75 * 1.75, 4.0}, 0.0, 21.41380852375568); +} + +TEST(StanMath_integrate_1d_impl_rev, + TestDerivatives_var_left_endpoint_var_params) { + // Zero crossing integral + test x_r + var at left endpoint + using stan::math::var; + test_derivatives( + f3{}, -1.0, 1.0, {0.5, 1.75, 3.9}, {2.5, 3.0}, {}, + 2.350402387287579 + 2.0 * pow(0.5, 2.5) + 4.0 * pow(1.75, 3.0) + + 4.0 * 3.9, + {5 * pow(0.5, 1.5), 12 * 1.75 * 1.75, 4.0}, -19.06340613646808, 0.0); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_no_param_vars) { + // No param vars + using stan::math::var; + test_derivatives(f3{}, -1.0, 1.0, {0.5, 1.75, 3.9}, + {2.5, 3.0}, {}, + 2.350402387287579 + 2.0 * pow(0.5, 2.5) + + 4.0 * pow(1.75, 3.0) + 4.0 * 3.9, + {}, -19.06340613646808, 21.41380852375568); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_left_limit_var) { + // No param vars, only left limit var + using stan::math::var; + test_derivatives(f3{}, -1.0, 1.0, {0.5, 1.75, 3.9}, + {2.5, 3.0}, {}, + 2.350402387287579 + 2.0 * pow(0.5, 2.5) + + 4.0 * pow(1.75, 3.0) + 4.0 * 3.9, + {}, -19.06340613646808, 0.0); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_right_limit_var) { + // No param vars, only right limit var + using stan::math::var; + test_derivatives(f3{}, -1.0, 1.0, {0.5, 1.75, 3.9}, + {2.5, 3.0}, {}, + 2.350402387287579 + 2.0 * pow(0.5, 2.5) + + 4.0 * pow(1.75, 3.0) + 4.0 * 3.9, + {}, 0.0, 21.41380852375568); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_tricky1) { + // Tricky integral from Boost docs + limit at infinity + no gradients + using stan::math::var; + test_derivatives(f4{}, 0.0, + std::numeric_limits::infinity(), + {}, {}, {}, 1.772453850905516, {}); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_tricky2) { + // Tricky integral from Boost docs + limit at infinity with gradients + using stan::math::var; + test_derivatives( + f5{}, 0.0, std::numeric_limits::infinity(), {0.5, 3.0}, {}, {}, + 1.772453850905516 / sqrt(0.5 * 3.0), + {-1.772453850905516 * 3.0 / (2 * pow(0.5 * 3.0, 1.5)), + -1.772453850905516 * 0.5 / (2 * pow(0.5 * 3.0, 1.5))}); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_tricky3) { + // Tricky integral from Boost docs + using stan::math::var; + test_derivatives( + f6{}, 0.0, 1.0, {0.75}, {}, {}, 0.851926727945904, {0.4814066053874294}); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_zero_crossing2) { + // Zero crossing integral + limit at infinity + var at left limit + using stan::math::var; + test_derivatives( + f7{}, -5.0, std::numeric_limits::infinity(), {1.5}, {}, {}, + 1205.361609637375, {5223.23364176196}, -1808.042414456063, + std::numeric_limits::quiet_NaN()); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_zero_crossing3) { + // Zero crossing integral + limit at negative infinity + var at right limit + using stan::math::var; + test_derivatives( + f8{}, -std::numeric_limits::infinity(), 5.0, {1.5}, {}, {}, + 1205.361609637375, {5223.23364176196}, + std::numeric_limits::quiet_NaN(), 1808.042414456063); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_indefinite) { + // Both limits at infinity + test x_r/x_i + no gradients + using stan::math::var; + test_derivatives( + f10{}, -std::numeric_limits::infinity(), + std::numeric_limits::infinity(), {}, {1.7}, {4}, + 2.536571480364399, {}); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_endpoint_precision) { + // Various integrals of beta function + using stan::math::var; + test_derivatives(f11{}, 0.0, 1.0, {0.1, 0.1}, {}, {}, + 19.71463948905016, + {-101.229055967892, -101.229055967892}); + test_derivatives( + f11{}, 0.0, 1.0, {0.5, 0.51}, {}, {}, 3.098843783331868, + {-4.346514423368675, -4.196150770134913}); + test_derivatives( + f11{}, 0.0, 1.0, {0.51, 0.5}, {}, {}, 3.098843783331868, + {-4.196150770134913, -4.346514423368675}); + test_derivatives( + f11{}, 0.0, 1.0, {5.0, 3.0}, {}, {}, 0.00952380952380952, + {-0.004852607709750566, -0.01040816326530613}); + test_derivatives( + f11{}, 0.0, 1.0, {3.0, 5.0}, {}, {}, 0.00952380952380952, + {-0.01040816326530613, -0.004852607709750566}); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_gaussian) { + // Check Gaussian integrates to 1.0 always + using stan::math::var; + test_derivatives( + f12{}, -std::numeric_limits::infinity(), + std::numeric_limits::infinity(), {5.7, 1}, {}, {}, 1.0, + {0.0, 0.0}); +} + +TEST(StanMath_integrate_1d_impl_rev, + TestDerivativesSameVarAtEndpointAndInParams) { + using stan::math::var; + + var a = 2.0; + var b = 4.0; + std::vector thetas = {a, b}; + std::vector x_r = {}; + std::vector x_i = {}; + + var integral = stan::math::integrate_1d_impl(f13{}, a, b, 1e-8, msgs, thetas, + x_r, x_i); + integral.grad(); + + EXPECT_LT(std::abs(18.0 - integral.val()), 1e-8); + EXPECT_LT(std::abs(-6.0 - a.adj()), 1e-8); + EXPECT_LT(std::abs(12.0 - b.adj()), 1e-8); +} + +TEST(StanMath_integrate_1d_impl_rev, TestBeta) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var alpha = 9.0 / 5; + var beta = 13.0 / 7; + AVEC theta = {alpha, beta}; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::beta_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, 0.0, 1.0, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(alpha, beta); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestCauchy) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var mu = 9.0 / 5; + var sigma = 13.0 / 7; + AVEC theta = {mu, sigma}; + double b = std::numeric_limits::infinity(); + double a = -b; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::cauchy_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(mu, sigma); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestChiSquare) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var nu = 9.0 / 5; + AVEC theta = {nu}; + double b = std::numeric_limits::infinity(); + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf + = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { return exp(stan::math::chi_square_lpdf(x, theta[0])); }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(nu); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestDoubleExponential) { + using stan::math::exp; + using stan::math::integrate_1d; + using stan::math::var; + + var mu = 9.0 / 5; + var sigma = 13.0 / 7; + AVEC theta = {mu, sigma}; + double a = -std::numeric_limits::infinity(); + double b = mu.val(); + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::double_exponential_lpdf(x, theta[0], theta[1])); + }; + // requires two subintervals to achieve numerical accuracy + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i) + + integrate_1d_impl(pdf, b, -a, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(mu, sigma); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestExponential) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var beta = 9.0 / 5; + AVEC theta = {beta}; + double b = std::numeric_limits::infinity(); + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf + = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { return exp(stan::math::exponential_lpdf(x, theta[0])); }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(beta); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestFrechet) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var alpha = 9.0 / 5; + var sigma = 13.0 / 7; + AVEC theta = {alpha, sigma}; + double b = std::numeric_limits::infinity(); + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::frechet_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(alpha, sigma); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestGamma) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var alpha = 9.0 / 5; + var beta = 13.0 / 7; + AVEC theta = {alpha, beta}; + double b = std::numeric_limits::infinity(); + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::gamma_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(alpha, beta); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestGumbel) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var mu = 9.0 / 5; + var beta = 13.0 / 7; + AVEC theta = {mu, beta}; + double b = std::numeric_limits::infinity(); + double a = -b; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::gumbel_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(mu, beta); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestInvChiSquared) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var nu = 9.0 / 5; + AVEC theta = {nu}; + double b = std::numeric_limits::infinity(); + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::inv_chi_square_lpdf(x, theta[0])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(nu); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestLogistic) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var mu = 9.0 / 5; + var sigma = 13.0 / 7; + AVEC theta = {mu, sigma}; + double b = std::numeric_limits::infinity(); + double a = -b; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::logistic_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(mu, sigma); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestLogNormal) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var mu = 9.0 / 5; + var sigma = 13.0 / 7; + AVEC theta = {mu, sigma}; + double b = std::numeric_limits::infinity(); + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::lognormal_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(mu, sigma); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestNormal) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var mu = 9.0 / 5; + var sigma = 13.0 / 7; + AVEC theta = {mu, sigma}; + double b = std::numeric_limits::infinity(); + double a = -b; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::normal_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(mu, sigma); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestPareto) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var m = 9.0 / 5; + var alpha = 13.0 / 7; + AVEC theta = {m, alpha}; + double b = std::numeric_limits::infinity(); + var a = m; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::pareto_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(m, alpha); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestPareto2) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var mu = 9.0 / 5; + var lambda = 13.0 / 7; + var alpha = 11.0 / 3; + AVEC theta = {mu, lambda, alpha}; + double b = std::numeric_limits::infinity(); + var a = mu; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::pareto_type_2_lpdf(x, theta[0], theta[1], theta[2])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(mu, lambda, alpha); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); + EXPECT_FLOAT_EQ(1, 1 + g[2]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestRayleigh) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var sigma = 13.0 / 7; + AVEC theta = {sigma}; + double b = std::numeric_limits::infinity(); + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf + = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { return exp(stan::math::rayleigh_lpdf(x, theta[0])); }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(sigma); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestScaledInvChiSquare) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var nu = 9.0 / 5; + var s = 13.0 / 7; + AVEC theta = {nu, s}; + double b = std::numeric_limits::infinity(); + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::scaled_inv_chi_square_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(nu, s); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestStudentT) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var nu = 11.0 / 3; + var mu = 9.0 / 5; + var sigma = 13.0 / 7; + AVEC theta = {nu, mu, sigma}; + double b = std::numeric_limits::infinity(); + double a = -b; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::student_t_lpdf(x, theta[0], theta[1], theta[2])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(nu, mu, sigma); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); + EXPECT_FLOAT_EQ(1, 1 + g[2]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestUniform) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var a = 9.0 / 5; + var b = 13.0 / 7; + std::vector x_r = {}; + std::vector x_i = {}; + AVEC theta = {a, b}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::uniform_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(a, b); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestVonMises) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var mu = 9.0 / 5; + var kappa = 13.0 / 7; + AVEC theta = {mu, kappa}; + double b = stan::math::pi() * 2; + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::von_mises_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(mu, kappa); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +} + +TEST(StanMath_integrate_1d_impl_rev, TestWeibull) { + using stan::math::exp; + using stan::math::integrate_1d_impl; + using stan::math::var; + + var alpha = 9.0 / 5; + var sigma = 13.0 / 7; + AVEC theta = {alpha, sigma}; + double b = std::numeric_limits::infinity(); + double a = 0; + std::vector x_r = {}; + std::vector x_i = {}; + auto pdf = [](auto x, auto xc, std::ostream *msgs, auto theta, auto x_r, + auto x_i) { + return exp(stan::math::weibull_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d_impl(pdf, a, b, 1e-8, msgs, theta, x_r, x_i); + EXPECT_FLOAT_EQ(1, I.val()); + + AVEC x = createAVEC(alpha, sigma); + VEC g; + I.grad(x, g); + EXPECT_FLOAT_EQ(1, 1 + g[0]); + EXPECT_FLOAT_EQ(1, 1 + g[1]); +}