From 6e0f680319644346cec05e9b61e899c5b3233668 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 22 Nov 2020 14:24:37 -0500 Subject: [PATCH 01/20] Saving work --- stan/math/prim/functor/integrate_1d.hpp | 21 + .../prim/functor/integrate_1d_new_test.cpp | 450 ++++++++++++++++++ 2 files changed, 471 insertions(+) create mode 100644 test/unit/math/prim/functor/integrate_1d_new_test.cpp diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index cc564fb6dd6..537a92e31ce 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -205,6 +205,27 @@ inline double integrate_1d(const F& f, const double a, const double b, } } +template +inline double integrate_1d_new(const F& f, double a, double b, + double relative_tolerance, + std::ostream* msgs, const Args&... args) { + //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, + msgs, args...), + a, b, relative_tolerance); + } +} + } // namespace math } // namespace stan diff --git a/test/unit/math/prim/functor/integrate_1d_new_test.cpp b/test/unit/math/prim/functor/integrate_1d_new_test.cpp new file mode 100644 index 00000000000..bc84975673f --- /dev/null +++ b/test/unit/math/prim/functor/integrate_1d_new_test.cpp @@ -0,0 +1,450 @@ +#include +#include +#include +#include +#include +#include +#include + +namespace integrate_1d_new_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 stan::return_type_t operator()(const T1 &x, const T1 &xc, + std::ostream *msgs, + const std::vector &theta, + 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_new(rank_density_functor__(), + down, up, + 1e-8, pstream__, + theta, x_r, x_i); + return v; +} +} // namespace integrate_1d_new_test +/* + * test_integration is a helper function to make it easy to test the + * integrate_1d_new 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_new; + + std::vector tolerances = {1e-4, 1e-6, 1e-8}; + + for (auto tolerance : tolerances) { + EXPECT_LE(std::abs(integrate_1d_new(f, a, b, tolerance, + integrate_1d_new_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_new(flipped, -b, -a, tolerance, + integrate_1d_new_test::msgs, args...) + - val), + tolerance); + } +} + +/*TEST(StanMath_integrate_1d_new_prim, TestThrows) { + // Left limit of integration must be less than or equal to right limit + EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 1.0, 0.0, + std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6), + std::domain_error); + // NaN limits not okay + EXPECT_THROW( + stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0.0, + std::numeric_limits::quiet_NaN(), + std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6), + std::domain_error); + EXPECT_THROW( + stan::math::integrate_1d_new( + integrate_1d_new_test::f2{}, std::numeric_limits::quiet_NaN(), + 0.0, std::vector(), {}, {}, integrate_1d_new_test::msgs, 1e-6), + std::domain_error); + EXPECT_THROW( + stan::math::integrate_1d_new( + integrate_1d_new_test::f2{}, std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN(), std::vector(), {}, + {}, integrate_1d_new_test::msgs, 1e-6), + std::domain_error); + // Two of the same inf limits not okay + EXPECT_THROW( + stan::math::integrate_1d_new( + integrate_1d_new_test::f2{}, -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), std::vector(), {}, + {}, integrate_1d_new_test::msgs, 1e-6), + std::domain_error); + + EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6), + std::domain_error); + // xc should be nan if there are infinite limits + EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, 0.0, + std::numeric_limits::infinity(), + std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6), + std::runtime_error); + EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, + std::numeric_limits::infinity(), + 0.0, std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6), + std::domain_error); + EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6), + std::domain_error); + // But not otherwise + EXPECT_NO_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, 0.0, 1.0, + std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6)); +} + +TEST(StanMath_integrate_1d_new_prim, test_integer_arguments) { + double v; + EXPECT_NO_THROW(v = stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0, 1, + std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6)); + EXPECT_NO_THROW(v = stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0.0, 1, + std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6)); + EXPECT_NO_THROW(v = stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0, 1.0, + std::vector(), {}, {}, + integrate_1d_new_test::msgs, 1e-6)); + }*/ + +TEST(StanMath_integrate_1d_new_prim, test1) { + // Tricky integral from Boost docs + limit at infinity + test_integration(integrate_1d_new_test::f1{}, 0.0, + std::numeric_limits::infinity(), + 1.772453850905516); + // Tricky integral from Boost 1d integration docs + test_integration(integrate_1d_new_test::f2{}, 0.0, 1.0, 1.198140234735592); + // Tricky integral from Boost 1d integration docs + test_integration(integrate_1d_new_test::f2{}, 0, 1, 1.198140234735592); + // Zero crossing integral + limit at infinity + test_integration(integrate_1d_new_test::f3{}, -2.0, + std::numeric_limits::infinity(), 7.38905609893065); + // Easy integrals + /*test_integration(integrate_1d_new_test::f4{}, 0.2, 0.7, {0.5}, {}, {}, + 1.0423499493102901); + test_integration(integrate_1d_new_test::f5{}, -0.2, 0.7, {0.4, 0.4}, {}, {}, + 1.396621954392482); + test_integration(integrate_1d_new_test::f4{}, 0.0, 0.0, {0.5}, {}, {}, 0.0); + test_integration(integrate_1d_new_test::f5{}, 1.0, 1.0, {0.4, 0.4}, {}, {}, 0.0); + // Test x_i + test_integration(integrate_1d_new_test::f6{}, -0.2, 2.9, {6.0, 5.1}, {}, {4}, + 4131.985414616364); + // Test x_r + test_integration(integrate_1d_new_test::f7{}, -0.2, 2.9, {}, {4.0, 6.0, 5.1}, {}, + 24219.985414616367); + // Both limits at infinity + test x_r/x_i + test_integration(integrate_1d_new_test::f8{}, + -std::numeric_limits::infinity(), + std::numeric_limits::infinity(), {5.0}, {1.7}, {2}, + 3.013171546539377); + // Both limits at infinity + test x_i + test_integration(integrate_1d_new_test::f9{}, + -std::numeric_limits::infinity(), + std::numeric_limits::infinity(), {1.3}, {}, {4}, + 2.372032924895055); + // Various integrals of beta function + test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, {0.1, 0.1}, {}, {}, + 19.71463948905016); + test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, {0.1, 0.5}, {}, {}, + 11.32308697521577); + test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, {0.5, 0.1}, {}, {}, + 11.32308697521577); + test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, {5.0, 3.0}, {}, {}, + 0.00952380952380952); + + // Integrals from + // http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/dhb-tanh-sinh.pdf + test_integration(integrate_1d_new_test::f12{}, 0.0, + std::numeric_limits::infinity(), {}, {}, {}, 2.0); + test_integration(integrate_1d_new_test::f13{}, 0.0, + std::numeric_limits::infinity(), {}, {}, {}, 1.0); + test_integration(integrate_1d_new_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_new_test::f16{}, 0.0, stan::math::pi(), {}, {}, {}, + stan::math::square(stan::math::pi()) / 4);*/ +} + +/*TEST(StanMath_integrate_1d_new_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_new_test::order(-10, 0.67, theta, x_r, msgs)); + }*/ + From 340e62ca812039401efeaa17366fc9326de80f08 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 23 Nov 2020 16:46:24 -0500 Subject: [PATCH 02/20] integrate_1d_new working (variadic integrate_1d) (Issue #2110) --- stan/math/prim/functor.hpp | 1 + stan/math/prim/functor/integrate_1d.hpp | 67 +- stan/math/rev/functor/integrate_1d.hpp | 163 +++- .../prim/functor/integrate_1d_new_test.cpp | 140 ++- .../rev/functor/integrate_1d_new_test.cpp | 897 ++++++++++++++++++ 5 files changed, 1113 insertions(+), 155 deletions(-) create mode 100644 test/unit/math/rev/functor/integrate_1d_new_test.cpp diff --git a/stan/math/prim/functor.hpp b/stan/math/prim/functor.hpp index aab4147aa96..d64b2666617 100644 --- a/stan/math/prim/functor.hpp +++ b/stan/math/prim/functor.hpp @@ -12,6 +12,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 537a92e31ce..94115e7d17b 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 @@ -136,6 +137,28 @@ inline double integrate(const F& f, double a, double b, return Q; } +template * = nullptr> +inline double integrate_1d_new(const F& f, double a, double b, + double relative_tolerance, + std::ostream* msgs, const Args&... args) { + //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, + msgs, args...), + a, b, relative_tolerance); + } +} + /** * 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. @@ -183,47 +206,15 @@ 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); - } -} - -template -inline double integrate_1d_new(const F& f, double a, double b, - double relative_tolerance, - std::ostream* msgs, const Args&... args) { - //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, - msgs, args...), - a, b, relative_tolerance); - } + return integrate_1d_new(integrate_1d_adapter(f), a, b, + relative_tolerance, msgs, + theta, x_r, x_i); } } // namespace math diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index a02e36b57ad..cd3beaf236c 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -59,6 +59,120 @@ inline double gradient_of_f(const F &f, const double &x, const double &xc, return gradient; } +template +inline double gradient_of_f_new(const F &f, const double &x, const double &xc, + size_t n, std::ostream *msgs, + const Args&... args) { + double gradient = 0.0; + + // Run nested autodiff in this scope + nested_rev_autodiff nested; + + std::tuple + args_tuple_local_copy(deep_copy_vars(args)...); + + Eigen::VectorXd adjoints = Eigen::VectorXd::Zero(count_vars(args...)); + + var fx = apply([&f, &x, &xc, msgs](auto&&... args) { + return f(x, xc, msgs, args...); + }, args_tuple_local_copy); + + fx.grad(); + + apply([&](auto&&... args) { + accumulate_adjoints(adjoints.data(), + std::forward(args)...); + }, + std::move(args_tuple_local_copy)); + + gradient = adjoints.coeff(n); + 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; +} + + +template * = nullptr> +inline return_type_t integrate_1d_new( + const F &f, const T_a &a, const T_b &b, + double relative_tolerance, + std::ostream *msgs, const Args&... args) { + static 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 (a_val == b_val) { + if (is_inf(a_val)) { + throw_domain_error(function, "Integration endpoints are both", + a_val, "", ""); + } + return var(0.0); + } else { + std::tuple + args_val_tuple(value_of(args)...); + + double integral = integrate(apply([&f, msgs](auto&&... args) { + return std::bind(f, std::placeholders::_1, std::placeholders::_2, + msgs, args...); + }, args_val_tuple), + a_val, b_val, relative_tolerance); + + size_t num_vars_ab = count_vars(a, b); + 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); + 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_inf(a) && is_var::value) { + *partials_ptr = apply([&f, a_val, msgs](auto&&... args) { + return -f(a_val, 0.0, msgs, args...); + }, args_val_tuple); + partials_ptr++; + } + + if (!is_inf(b) && is_var::value) { + *partials_ptr = apply([&f, b_val, msgs](auto&&... args) { + return f(b_val, 0.0, msgs, args...); + }, args_val_tuple); + partials_ptr++; + } + + for (size_t n = 0; n < num_vars_args; ++n) { + *partials_ptr = integrate( + std::bind(gradient_of_f_new, f, + std::placeholders::_1, std::placeholders::_2, + n, msgs, args...), + a_val, b_val, relative_tolerance); + partials_ptr++; + } + + return var(new precomputed_gradients_vari(integral, + num_vars_ab + num_vars_args, + varis, partials)); + } +} + /** * 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. @@ -120,52 +234,9 @@ 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_new(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_new_test.cpp b/test/unit/math/prim/functor/integrate_1d_new_test.cpp index bc84975673f..7daadabbb68 100644 --- a/test/unit/math/prim/functor/integrate_1d_new_test.cpp +++ b/test/unit/math/prim/functor/integrate_1d_new_test.cpp @@ -67,11 +67,10 @@ struct f6 { }; struct f7 { - 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 { + 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]; } }; @@ -140,7 +139,7 @@ struct f13 { struct f14 { template inline T1 operator()(const T1 &x, const T1 &xc, - std::ostream *msgs) const { + std::ostream *msgs) const { return exp(x) * stan::math::inv_sqrt(x > 0.5 ? xc : 1 - x); } }; @@ -158,7 +157,7 @@ struct f15 { }; struct f16 { - template + template inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { return x * sin(x) / (1 + stan::math::square(cos(x))); @@ -295,79 +294,70 @@ void test_integration(const F &f, double a, double b, } } -/*TEST(StanMath_integrate_1d_new_prim, TestThrows) { +TEST(StanMath_integrate_1d_new_prim, TestThrows) { // Left limit of integration must be less than or equal to right limit EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 1.0, 0.0, - std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6), + 0.0, integrate_1d_new_test::msgs), std::domain_error); // NaN limits not okay EXPECT_THROW( stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0.0, - std::numeric_limits::quiet_NaN(), - std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6), + std::numeric_limits::quiet_NaN(), + 0.0, integrate_1d_new_test::msgs), std::domain_error); EXPECT_THROW( stan::math::integrate_1d_new( - integrate_1d_new_test::f2{}, std::numeric_limits::quiet_NaN(), - 0.0, std::vector(), {}, {}, integrate_1d_new_test::msgs, 1e-6), + integrate_1d_new_test::f2{}, + std::numeric_limits::quiet_NaN(), + 0.0, 0.0, integrate_1d_new_test::msgs), std::domain_error); EXPECT_THROW( stan::math::integrate_1d_new( integrate_1d_new_test::f2{}, std::numeric_limits::quiet_NaN(), - std::numeric_limits::quiet_NaN(), std::vector(), {}, - {}, integrate_1d_new_test::msgs, 1e-6), + std::numeric_limits::quiet_NaN(), 0.0, + integrate_1d_new_test::msgs), std::domain_error); // Two of the same inf limits not okay EXPECT_THROW( stan::math::integrate_1d_new( integrate_1d_new_test::f2{}, -std::numeric_limits::infinity(), - -std::numeric_limits::infinity(), std::vector(), {}, - {}, integrate_1d_new_test::msgs, 1e-6), + -std::numeric_limits::infinity(), 0.0, + integrate_1d_new_test::msgs), std::domain_error); EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, - std::numeric_limits::infinity(), - std::numeric_limits::infinity(), - std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6), + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + 0.0, integrate_1d_new_test::msgs), std::domain_error); // xc should be nan if there are infinite limits EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, 0.0, - std::numeric_limits::infinity(), - std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6), + std::numeric_limits::infinity(), + 0.0, integrate_1d_new_test::msgs), std::runtime_error); EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, std::numeric_limits::infinity(), - 0.0, std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6), + 0.0, 0.0, integrate_1d_new_test::msgs), std::domain_error); EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, - std::numeric_limits::infinity(), - std::numeric_limits::infinity(), - std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6), + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + 0.0, integrate_1d_new_test::msgs), std::domain_error); // But not otherwise EXPECT_NO_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, 0.0, 1.0, - std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6)); + 0.0, integrate_1d_new_test::msgs)); } TEST(StanMath_integrate_1d_new_prim, test_integer_arguments) { double v; EXPECT_NO_THROW(v = stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0, 1, - std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6)); + 0.0, integrate_1d_new_test::msgs)); EXPECT_NO_THROW(v = stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0.0, 1, - std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6)); + 0.0, integrate_1d_new_test::msgs)); EXPECT_NO_THROW(v = stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0, 1.0, - std::vector(), {}, {}, - integrate_1d_new_test::msgs, 1e-6)); - }*/ + 0.0, integrate_1d_new_test::msgs)); +} TEST(StanMath_integrate_1d_new_prim, test1) { // Tricky integral from Boost docs + limit at infinity @@ -382,45 +372,54 @@ TEST(StanMath_integrate_1d_new_prim, test1) { test_integration(integrate_1d_new_test::f3{}, -2.0, std::numeric_limits::infinity(), 7.38905609893065); // Easy integrals - /*test_integration(integrate_1d_new_test::f4{}, 0.2, 0.7, {0.5}, {}, {}, - 1.0423499493102901); - test_integration(integrate_1d_new_test::f5{}, -0.2, 0.7, {0.4, 0.4}, {}, {}, - 1.396621954392482); - test_integration(integrate_1d_new_test::f4{}, 0.0, 0.0, {0.5}, {}, {}, 0.0); - test_integration(integrate_1d_new_test::f5{}, 1.0, 1.0, {0.4, 0.4}, {}, {}, 0.0); + test_integration(integrate_1d_new_test::f4{}, 0.2, 0.7, + 1.0423499493102901, std::vector({0.5})); + test_integration(integrate_1d_new_test::f5{}, -0.2, 0.7, + 1.396621954392482, std::vector({0.4, 0.4})); + test_integration(integrate_1d_new_test::f4{}, 0.0, 0.0, 0.0, + std::vector({0.5})); + test_integration(integrate_1d_new_test::f5{}, 1.0, 1.0, 0.0, + std::vector({0.4, 0.4})); // Test x_i - test_integration(integrate_1d_new_test::f6{}, -0.2, 2.9, {6.0, 5.1}, {}, {4}, - 4131.985414616364); + test_integration(integrate_1d_new_test::f6{}, -0.2, 2.9, + 4131.985414616364, std::vector({6.0, 5.1}), + std::vector({4})); // Test x_r - test_integration(integrate_1d_new_test::f7{}, -0.2, 2.9, {}, {4.0, 6.0, 5.1}, {}, - 24219.985414616367); + test_integration(integrate_1d_new_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_new_test::f8{}, -std::numeric_limits::infinity(), - std::numeric_limits::infinity(), {5.0}, {1.7}, {2}, - 3.013171546539377); + 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_new_test::f9{}, -std::numeric_limits::infinity(), - std::numeric_limits::infinity(), {1.3}, {}, {4}, - 2.372032924895055); + std::numeric_limits::infinity(), + 2.372032924895055, + std::vector({1.3}), + std::vector({4})); // Various integrals of beta function - test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, {0.1, 0.1}, {}, {}, - 19.71463948905016); - test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, {0.1, 0.5}, {}, {}, - 11.32308697521577); - test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, {0.5, 0.1}, {}, {}, - 11.32308697521577); - test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, {5.0, 3.0}, {}, {}, - 0.00952380952380952); + test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, + 19.71463948905016, std::vector({0.1, 0.1})); + test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, + 11.32308697521577, std::vector({0.1, 0.5})); + test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, + 11.32308697521577, std::vector({0.5, 0.1})); + test_integration(integrate_1d_new_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_new_test::f12{}, 0.0, - std::numeric_limits::infinity(), {}, {}, {}, 2.0); + std::numeric_limits::infinity(), 2.0); test_integration(integrate_1d_new_test::f13{}, 0.0, - std::numeric_limits::infinity(), {}, {}, {}, 1.0); - test_integration(integrate_1d_new_test::f14{}, 0.0, 1.0, {}, {}, {}, + std::numeric_limits::infinity(), 1.0); + test_integration(integrate_1d_new_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 @@ -428,11 +427,11 @@ TEST(StanMath_integrate_1d_new_prim, test1) { // test_integration(f15{}, 0.0, 1.0, {}, {}, {}, // stan::math::square(stan::math::pi()) * (2 - sqrt(2.0)) / // 32); - test_integration(integrate_1d_new_test::f16{}, 0.0, stan::math::pi(), {}, {}, {}, - stan::math::square(stan::math::pi()) / 4);*/ + test_integration(integrate_1d_new_test::f16{}, 0.0, stan::math::pi(), + stan::math::square(stan::math::pi()) / 4); } -/*TEST(StanMath_integrate_1d_new_prim, TestTolerance) { +TEST(StanMath_integrate_1d_new_prim, TestTolerance) { std::ostringstream *msgs = nullptr; double t = 0.5; @@ -446,5 +445,4 @@ TEST(StanMath_integrate_1d_new_prim, test1) { std::vector x_r; EXPECT_NO_THROW(integrate_1d_new_test::order(-10, 0.67, theta, x_r, msgs)); - }*/ - +} diff --git a/test/unit/math/rev/functor/integrate_1d_new_test.cpp b/test/unit/math/rev/functor/integrate_1d_new_test.cpp new file mode 100644 index 00000000000..4f417fee777 --- /dev/null +++ b/test/unit/math/rev/functor/integrate_1d_new_test.cpp @@ -0,0 +1,897 @@ +#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_new(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_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_new(f2{}, 0, 1, 1e-6, msgs, theta, x_r, x_i)); + EXPECT_NO_THROW( + v = stan::math::integrate_1d_new(f2{}, 0.0, 1, 1e-6, msgs, theta, x_r, x_i)); + EXPECT_NO_THROW( + v = stan::math::integrate_1d_new(f2{}, 0, 1.0, 1e-6, msgs, theta, x_r, x_i)); +} + +TEST(StanMath_integrate_1d_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_rev, TestDerivativesSameVarAtEndpointAndInParams) { + using stan::math::var; + + var a = 2.0; + var b = 4.0; + std::vector thetas = {a, b}; + + var integral + = stan::math::integrate_1d(f13{}, a, b, thetas, {}, {}, msgs, 1e-8); + 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_rev, TestBeta) { + using stan::math::exp; + using stan::math::integrate_1d; + using stan::math::var; + + var alpha = 9.0 / 5; + var beta = 13.0 / 7; + AVEC theta = {alpha, beta}; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::beta_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, 0.0, 1.0, theta, {}, {}, msgs, 1e-8); + 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_rev, TestCauchy) { + 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 b = std::numeric_limits::infinity(); + double a = -b; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::cauchy_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestChiSquare) { + using stan::math::exp; + using stan::math::integrate_1d; + using stan::math::var; + + var nu = 9.0 / 5; + AVEC theta = {nu}; + double b = std::numeric_limits::infinity(); + double a = 0; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::chi_square_lpdf(x, theta[0])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_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(); + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::double_exponential_lpdf(x, theta[0], theta[1])); + }; + // requires two subintervals to achieve numerical accuracy + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8) + + integrate_1d(pdf, b, -a, theta, {}, {}, msgs, 1e-8); + 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_rev, TestExponential) { + using stan::math::exp; + using stan::math::integrate_1d; + using stan::math::var; + + var beta = 9.0 / 5; + AVEC theta = {beta}; + double b = std::numeric_limits::infinity(); + double a = 0; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::exponential_lpdf(x, theta[0])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestFrechet) { + using stan::math::exp; + using stan::math::integrate_1d; + 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; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::frechet_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestGamma) { + using stan::math::exp; + using stan::math::integrate_1d; + 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; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::gamma_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestGumbel) { + using stan::math::exp; + using stan::math::integrate_1d; + 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; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::gumbel_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestInvChiSquared) { + using stan::math::exp; + using stan::math::integrate_1d; + using stan::math::var; + + var nu = 9.0 / 5; + AVEC theta = {nu}; + double b = std::numeric_limits::infinity(); + double a = 0; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::inv_chi_square_lpdf(x, theta[0])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestLogistic) { + 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 b = std::numeric_limits::infinity(); + double a = -b; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::logistic_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestLogNormal) { + 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 b = std::numeric_limits::infinity(); + double a = 0; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::lognormal_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestNormal) { + 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 b = std::numeric_limits::infinity(); + double a = -b; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::normal_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestPareto) { + using stan::math::exp; + using stan::math::integrate_1d; + 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; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::pareto_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestPareto2) { + using stan::math::exp; + using stan::math::integrate_1d; + 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; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::pareto_type_2_lpdf(x, theta[0], theta[1], theta[2])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestRayleigh) { + using stan::math::exp; + using stan::math::integrate_1d; + using stan::math::var; + + var sigma = 13.0 / 7; + AVEC theta = {sigma}; + double b = std::numeric_limits::infinity(); + double a = 0; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::rayleigh_lpdf(x, theta[0])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestScaledInvChiSquare) { + using stan::math::exp; + using stan::math::integrate_1d; + 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; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::scaled_inv_chi_square_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestStudentT) { + using stan::math::exp; + using stan::math::integrate_1d; + 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; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::student_t_lpdf(x, theta[0], theta[1], theta[2])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestUniform) { + using stan::math::exp; + using stan::math::integrate_1d; + using stan::math::var; + + var a = 9.0 / 5; + var b = 13.0 / 7; + AVEC theta = {a, b}; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::uniform_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestVonMises) { + using stan::math::exp; + using stan::math::integrate_1d; + 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; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::von_mises_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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_rev, TestWeibull) { + using stan::math::exp; + using stan::math::integrate_1d; + 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; + auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, + std::ostream *msgs) { + return exp(stan::math::weibull_lpdf(x, theta[0], theta[1])); + }; + var I = integrate_1d(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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]); +} +*/ From f84f7eb2c37f1adef86947c72796cb2e0b3ba7bb Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 24 Feb 2021 15:29:03 -0500 Subject: [PATCH 03/20] Bit of cleanup and added adapter file (Issue #2197) --- stan/math/prim/functor/integrate_1d.hpp | 26 +++++-- .../prim/functor/integrate_1d_adapter.hpp | 29 ++++++++ stan/math/rev/functor/integrate_1d.hpp | 69 ++++++++----------- 3 files changed, 78 insertions(+), 46 deletions(-) create mode 100644 stan/math/prim/functor/integrate_1d_adapter.hpp diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index 803f98c99a8..31c428c54c8 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -132,12 +132,24 @@ inline double integrate(const F& f, double a, double b, 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 * = nullptr> -inline double integrate_1d_new(const F& f, double a, double b, - double relative_tolerance, - std::ostream* msgs, const Args&... args) { - //const double relative_tolerance = std::sqrt(EPSILON); +inline double integrate_1d_impl(const F& f, double a, double b, + double relative_tolerance, + std::ostream* msgs, const Args&... args) { static const char* function = "integrate_1d"; check_less_or_equal(function, "lower limit", a, b); @@ -207,9 +219,9 @@ inline double integrate_1d(const F& f, double a, double b, const std::vector& x_i, std::ostream* msgs, const double relative_tolerance = std::sqrt(EPSILON)) { - return integrate_1d_new(integrate_1d_adapter(f), a, b, - relative_tolerance, msgs, - theta, x_r, x_i); + 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..82033a70823 --- /dev/null +++ b/stan/math/prim/functor/integrate_1d_adapter.hpp @@ -0,0 +1,29 @@ +#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 cd3beaf236c..60cc3eb4b9f 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -28,41 +28,17 @@ namespace math { * NaN, a std::domain_error is thrown * * @tparam F type of f + * @tparam Args types of arguments to f + * @param f function to compute gradients of + * @param x location at which to evaluate gradients + * @param xc complement of location (if bounded domain of integration) + * @param n compute gradient with respect to nth parameter + * @param msgs stream for messages + * @param args other arguments to pass to 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 ", ""); - } - } - - return gradient; -} - template -inline double gradient_of_f_new(const F &f, const double &x, const double &xc, - size_t n, std::ostream *msgs, - const Args&... args) { +inline double gradient_of_f(const F &f, const double &x, const double &xc, + size_t n, std::ostream *msgs, const Args&... args) { double gradient = 0.0; // Run nested autodiff in this scope @@ -98,11 +74,26 @@ inline double gradient_of_f_new(const F &f, const double &x, const double &xc, return gradient; } - +/** + * Return the integral of f from a to b to the given relative tolerance + * + * @tparam T_a type of first limit + * @tparam T_b type of second limit + * @tparam T_theta type of parameters + * @tparam T 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 * = nullptr> -inline return_type_t integrate_1d_new( +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) { @@ -160,7 +151,7 @@ inline return_type_t integrate_1d_new( for (size_t n = 0; n < num_vars_args; ++n) { *partials_ptr = integrate( - std::bind(gradient_of_f_new, f, + std::bind(gradient_of_f, f, std::placeholders::_1, std::placeholders::_2, n, msgs, args...), a_val, b_val, relative_tolerance); @@ -234,9 +225,9 @@ 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)) { - return integrate_1d_new(integrate_1d_adapter(f), a, b, - relative_tolerance, msgs, - theta, x_r, x_i); + return integrate_1d_impl(integrate_1d_adapter(f), a, b, + relative_tolerance, msgs, + theta, x_r, x_i); } } // namespace math From 5019637932eeebde5f92754f2d11d2c5c90d559e Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 24 Feb 2021 15:36:10 -0500 Subject: [PATCH 04/20] Renamed tests (Issue #2197) --- ...ew_test.cpp => integrate_1d_impl_test.cpp} | 126 +++++++++--------- ...ew_test.cpp => integrate_1d_impl_test.cpp} | 8 +- 2 files changed, 67 insertions(+), 67 deletions(-) rename test/unit/math/prim/functor/{integrate_1d_new_test.cpp => integrate_1d_impl_test.cpp} (76%) rename test/unit/math/rev/functor/{integrate_1d_new_test.cpp => integrate_1d_impl_test.cpp} (98%) diff --git a/test/unit/math/prim/functor/integrate_1d_new_test.cpp b/test/unit/math/prim/functor/integrate_1d_impl_test.cpp similarity index 76% rename from test/unit/math/prim/functor/integrate_1d_new_test.cpp rename to test/unit/math/prim/functor/integrate_1d_impl_test.cpp index 7daadabbb68..9f20b385744 100644 --- a/test/unit/math/prim/functor/integrate_1d_new_test.cpp +++ b/test/unit/math/prim/functor/integrate_1d_impl_test.cpp @@ -6,7 +6,7 @@ #include #include -namespace integrate_1d_new_test { +namespace integrate_1d_impl_test { std::ostringstream *msgs = nullptr; @@ -235,16 +235,16 @@ double order(double down, double up, const std::vector &theta, double v; - v = stan::math::integrate_1d_new(rank_density_functor__(), + v = stan::math::integrate_1d_impl(rank_density_functor__(), down, up, 1e-8, pstream__, theta, x_r, x_i); return v; } -} // namespace integrate_1d_new_test +} // namespace integrate_1d_impl_test /* * test_integration is a helper function to make it easy to test the - * integrate_1d_new function. + * 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 @@ -271,13 +271,13 @@ double order(double down, double up, const std::vector &theta, template void test_integration(const F &f, double a, double b, double val, const Args&... args) { - using stan::math::integrate_1d_new; + 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_new(f, a, b, tolerance, - integrate_1d_new_test::msgs, args...) + 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 @@ -287,109 +287,109 @@ void test_integration(const F &f, double a, double b, return f(-x, -xc, msgs, args...); }; - EXPECT_LE(std::abs(integrate_1d_new(flipped, -b, -a, tolerance, - integrate_1d_new_test::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_new_prim, TestThrows) { +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_new(integrate_1d_new_test::f2{}, 1.0, 0.0, - 0.0, integrate_1d_new_test::msgs), + 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_new(integrate_1d_new_test::f2{}, 0.0, + stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, 0.0, std::numeric_limits::quiet_NaN(), - 0.0, integrate_1d_new_test::msgs), + 0.0, integrate_1d_impl_test::msgs), std::domain_error); EXPECT_THROW( - stan::math::integrate_1d_new( - integrate_1d_new_test::f2{}, + stan::math::integrate_1d_impl( + integrate_1d_impl_test::f2{}, std::numeric_limits::quiet_NaN(), - 0.0, 0.0, integrate_1d_new_test::msgs), + 0.0, 0.0, integrate_1d_impl_test::msgs), std::domain_error); EXPECT_THROW( - stan::math::integrate_1d_new( - integrate_1d_new_test::f2{}, std::numeric_limits::quiet_NaN(), + stan::math::integrate_1d_impl( + integrate_1d_impl_test::f2{}, std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN(), 0.0, - integrate_1d_new_test::msgs), + integrate_1d_impl_test::msgs), std::domain_error); // Two of the same inf limits not okay EXPECT_THROW( - stan::math::integrate_1d_new( - integrate_1d_new_test::f2{}, -std::numeric_limits::infinity(), + stan::math::integrate_1d_impl( + integrate_1d_impl_test::f2{}, -std::numeric_limits::infinity(), -std::numeric_limits::infinity(), 0.0, - integrate_1d_new_test::msgs), + integrate_1d_impl_test::msgs), std::domain_error); - EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, + EXPECT_THROW(stan::math::integrate_1d_impl(integrate_1d_impl_test::f2{}, std::numeric_limits::infinity(), std::numeric_limits::infinity(), - 0.0, integrate_1d_new_test::msgs), + 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_new(integrate_1d_new_test::f11{}, 0.0, + EXPECT_THROW(stan::math::integrate_1d_impl(integrate_1d_impl_test::f11{}, 0.0, std::numeric_limits::infinity(), - 0.0, integrate_1d_new_test::msgs), + 0.0, integrate_1d_impl_test::msgs), std::runtime_error); - EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, + EXPECT_THROW(stan::math::integrate_1d_impl(integrate_1d_impl_test::f11{}, std::numeric_limits::infinity(), - 0.0, 0.0, integrate_1d_new_test::msgs), + 0.0, 0.0, integrate_1d_impl_test::msgs), std::domain_error); - EXPECT_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, + EXPECT_THROW(stan::math::integrate_1d_impl(integrate_1d_impl_test::f11{}, std::numeric_limits::infinity(), std::numeric_limits::infinity(), - 0.0, integrate_1d_new_test::msgs), + 0.0, integrate_1d_impl_test::msgs), std::domain_error); // But not otherwise - EXPECT_NO_THROW(stan::math::integrate_1d_new(integrate_1d_new_test::f11{}, 0.0, 1.0, - 0.0, integrate_1d_new_test::msgs)); + 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_new_prim, test_integer_arguments) { +TEST(StanMath_integrate_1d_impl_prim, test_integer_arguments) { double v; - EXPECT_NO_THROW(v = stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0, 1, - 0.0, integrate_1d_new_test::msgs)); - EXPECT_NO_THROW(v = stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0.0, 1, - 0.0, integrate_1d_new_test::msgs)); - EXPECT_NO_THROW(v = stan::math::integrate_1d_new(integrate_1d_new_test::f2{}, 0, 1.0, - 0.0, integrate_1d_new_test::msgs)); + 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_new_prim, test1) { +TEST(StanMath_integrate_1d_impl_prim, test1) { // Tricky integral from Boost docs + limit at infinity - test_integration(integrate_1d_new_test::f1{}, 0.0, + 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_new_test::f2{}, 0.0, 1.0, 1.198140234735592); + test_integration(integrate_1d_impl_test::f2{}, 0.0, 1.0, 1.198140234735592); // Tricky integral from Boost 1d integration docs - test_integration(integrate_1d_new_test::f2{}, 0, 1, 1.198140234735592); + test_integration(integrate_1d_impl_test::f2{}, 0, 1, 1.198140234735592); // Zero crossing integral + limit at infinity - test_integration(integrate_1d_new_test::f3{}, -2.0, + test_integration(integrate_1d_impl_test::f3{}, -2.0, std::numeric_limits::infinity(), 7.38905609893065); // Easy integrals - test_integration(integrate_1d_new_test::f4{}, 0.2, 0.7, + test_integration(integrate_1d_impl_test::f4{}, 0.2, 0.7, 1.0423499493102901, std::vector({0.5})); - test_integration(integrate_1d_new_test::f5{}, -0.2, 0.7, + test_integration(integrate_1d_impl_test::f5{}, -0.2, 0.7, 1.396621954392482, std::vector({0.4, 0.4})); - test_integration(integrate_1d_new_test::f4{}, 0.0, 0.0, 0.0, + test_integration(integrate_1d_impl_test::f4{}, 0.0, 0.0, 0.0, std::vector({0.5})); - test_integration(integrate_1d_new_test::f5{}, 1.0, 1.0, 0.0, + 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_new_test::f6{}, -0.2, 2.9, + 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_new_test::f7{}, -0.2, 2.9, + 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_new_test::f8{}, + test_integration(integrate_1d_impl_test::f8{}, -std::numeric_limits::infinity(), std::numeric_limits::infinity(), 3.013171546539377, @@ -397,29 +397,29 @@ TEST(StanMath_integrate_1d_new_prim, test1) { std::vector({1.7}), std::vector({2})); // Both limits at infinity + test x_i - test_integration(integrate_1d_new_test::f9{}, + 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_new_test::f10{}, 0.0, 1.0, + test_integration(integrate_1d_impl_test::f10{}, 0.0, 1.0, 19.71463948905016, std::vector({0.1, 0.1})); - test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, + test_integration(integrate_1d_impl_test::f10{}, 0.0, 1.0, 11.32308697521577, std::vector({0.1, 0.5})); - test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, + test_integration(integrate_1d_impl_test::f10{}, 0.0, 1.0, 11.32308697521577, std::vector({0.5, 0.1})); - test_integration(integrate_1d_new_test::f10{}, 0.0, 1.0, + 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_new_test::f12{}, 0.0, + test_integration(integrate_1d_impl_test::f12{}, 0.0, std::numeric_limits::infinity(), 2.0); - test_integration(integrate_1d_new_test::f13{}, 0.0, + test_integration(integrate_1d_impl_test::f13{}, 0.0, std::numeric_limits::infinity(), 1.0); - test_integration(integrate_1d_new_test::f14{}, 0.0, 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 @@ -427,11 +427,11 @@ TEST(StanMath_integrate_1d_new_prim, test1) { // test_integration(f15{}, 0.0, 1.0, {}, {}, {}, // stan::math::square(stan::math::pi()) * (2 - sqrt(2.0)) / // 32); - test_integration(integrate_1d_new_test::f16{}, 0.0, stan::math::pi(), + test_integration(integrate_1d_impl_test::f16{}, 0.0, stan::math::pi(), stan::math::square(stan::math::pi()) / 4); } -TEST(StanMath_integrate_1d_new_prim, TestTolerance) { +TEST(StanMath_integrate_1d_impl_prim, TestTolerance) { std::ostringstream *msgs = nullptr; double t = 0.5; @@ -444,5 +444,5 @@ TEST(StanMath_integrate_1d_new_prim, TestTolerance) { std::vector theta = {t, A, v1, v2, s}; std::vector x_r; - EXPECT_NO_THROW(integrate_1d_new_test::order(-10, 0.67, theta, x_r, msgs)); + 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_new_test.cpp b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp similarity index 98% rename from test/unit/math/rev/functor/integrate_1d_new_test.cpp rename to test/unit/math/rev/functor/integrate_1d_impl_test.cpp index 4f417fee777..2288b4baa01 100644 --- a/test/unit/math/rev/functor/integrate_1d_new_test.cpp +++ b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp @@ -216,7 +216,7 @@ void test_derivatives(const F &f, double a, double b, for (size_t i = 0; i < thetas.size(); ++i) thetas_[i] = thetas[i]; - var integral = stan::math::integrate_1d_new(f, a_, b_, tolerance, msgs, thetas_, x_r, x_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) { @@ -239,11 +239,11 @@ TEST(StanMath_integrate_1d_rev, test_integer_arguments) { std::vector x_r; std::vector x_i; EXPECT_NO_THROW( - v = stan::math::integrate_1d_new(f2{}, 0, 1, 1e-6, msgs, theta, x_r, x_i)); + 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_new(f2{}, 0.0, 1, 1e-6, msgs, theta, x_r, x_i)); + 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_new(f2{}, 0, 1.0, 1e-6, msgs, theta, x_r, x_i)); + v = stan::math::integrate_1d_impl(f2{}, 0, 1.0, 1e-6, msgs, theta, x_r, x_i)); } TEST(StanMath_integrate_1d_rev, TestDerivatives_easy) { From 403987fdd44df09bd41b973f059cbf1e98975496 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 24 Feb 2021 16:23:09 -0500 Subject: [PATCH 05/20] Turned on reverse mode tests (Issue #2197) --- .../rev/functor/integrate_1d_impl_test.cpp | 279 ++++++++++-------- 1 file changed, 160 insertions(+), 119 deletions(-) diff --git a/test/unit/math/rev/functor/integrate_1d_impl_test.cpp b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp index 2288b4baa01..6b46a7d428a 100644 --- a/test/unit/math/rev/functor/integrate_1d_impl_test.cpp +++ b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp @@ -233,7 +233,7 @@ void test_derivatives(const F &f, double a, double b, } } -TEST(StanMath_integrate_1d_rev, test_integer_arguments) { +TEST(StanMath_integrate_1d_impl_rev, test_integer_arguments) { stan::math::var v; std::vector theta = {0.5}; std::vector x_r; @@ -246,7 +246,7 @@ TEST(StanMath_integrate_1d_rev, test_integer_arguments) { v = stan::math::integrate_1d_impl(f2{}, 0, 1.0, 1e-6, msgs, theta, x_r, x_i)); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_easy) { +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_easy) { // Easy integrals using stan::math::var; test_derivatives(f1{}, 0.2, 0.7, {0.75}, {}, {}, @@ -263,7 +263,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_easy) { {0.0}); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_zero_crossing) { +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}, @@ -274,7 +274,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_zero_crossing) { -19.06340613646808, 21.41380852375568); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_var_right_endpoint_var_params) { +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( @@ -284,7 +284,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_var_right_endpoint_var_params) { {5 * pow(0.5, 1.5), 12 * 1.75 * 1.75, 4.0}, 0.0, 21.41380852375568); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_var_left_endpoint_var_params) { +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( @@ -294,7 +294,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_var_left_endpoint_var_params) { {5 * pow(0.5, 1.5), 12 * 1.75 * 1.75, 4.0}, -19.06340613646808, 0.0); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_no_param_vars) { +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}, @@ -304,7 +304,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_no_param_vars) { {}, -19.06340613646808, 21.41380852375568); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_left_limit_var) { +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}, @@ -314,7 +314,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_left_limit_var) { {}, -19.06340613646808, 0.0); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_right_limit_var) { +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}, @@ -324,7 +324,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_right_limit_var) { {}, 0.0, 21.41380852375568); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_tricky1) { +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, @@ -332,7 +332,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_tricky1) { {}, {}, {}, 1.772453850905516, {}); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_tricky2) { +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_tricky2) { // Tricky integral from Boost docs + limit at infinity with gradients using stan::math::var; test_derivatives( @@ -342,14 +342,14 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_tricky2) { -1.772453850905516 * 0.5 / (2 * pow(0.5 * 3.0, 1.5))}); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_tricky3) { +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_rev, TestDerivatives_zero_crossing2) { +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( @@ -358,7 +358,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_zero_crossing2) { std::numeric_limits::quiet_NaN()); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_zero_crossing3) { +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( @@ -367,7 +367,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_zero_crossing3) { std::numeric_limits::quiet_NaN(), 1808.042414456063); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_indefinite) { +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( @@ -376,7 +376,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_indefinite) { 2.536571480364399, {}); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_endpoint_precision) { +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}, {}, {}, @@ -396,7 +396,7 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_endpoint_precision) { {-0.01040816326530613, -0.004852607709750566}); } -TEST(StanMath_integrate_1d_rev, TestDerivatives_gaussian) { +TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_gaussian) { // Check Gaussian integrates to 1.0 always using stan::math::var; test_derivatives( @@ -405,15 +405,17 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_gaussian) { {0.0, 0.0}); } -/*TEST(StanMath_integrate_1d_rev, TestDerivativesSameVarAtEndpointAndInParams) { +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(f13{}, a, b, thetas, {}, {}, msgs, 1e-8); + = 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); @@ -421,19 +423,21 @@ TEST(StanMath_integrate_1d_rev, TestDerivatives_gaussian) { EXPECT_LT(std::abs(12.0 - b.adj()), 1e-8); } -TEST(StanMath_integrate_1d_rev, TestBeta) { +TEST(StanMath_integrate_1d_impl_rev, TestBeta) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var alpha = 9.0 / 5; var beta = 13.0 / 7; AVEC theta = {alpha, beta}; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, 0.0, 1.0, theta, {}, {}, msgs, 1e-8); + 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); @@ -443,9 +447,9 @@ TEST(StanMath_integrate_1d_rev, TestBeta) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestCauchy) { +TEST(StanMath_integrate_1d_impl_rev, TestCauchy) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var mu = 9.0 / 5; @@ -453,11 +457,13 @@ TEST(StanMath_integrate_1d_rev, TestCauchy) { AVEC theta = {mu, sigma}; double b = std::numeric_limits::infinity(); double a = -b; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -467,20 +473,22 @@ TEST(StanMath_integrate_1d_rev, TestCauchy) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestChiSquare) { +TEST(StanMath_integrate_1d_impl_rev, TestChiSquare) { using stan::math::exp; - using stan::math::integrate_1d; + 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; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -489,7 +497,7 @@ TEST(StanMath_integrate_1d_rev, TestChiSquare) { EXPECT_FLOAT_EQ(1, 1 + g[0]); } -TEST(StanMath_integrate_1d_rev, TestDoubleExponential) { +TEST(StanMath_integrate_1d_impl_rev, TestDoubleExponential) { using stan::math::exp; using stan::math::integrate_1d; using stan::math::var; @@ -499,13 +507,15 @@ TEST(StanMath_integrate_1d_rev, TestDoubleExponential) { AVEC theta = {mu, sigma}; double a = -std::numeric_limits::infinity(); double b = mu.val(); - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8) - + integrate_1d(pdf, b, -a, theta, {}, {}, msgs, 1e-8); + 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); @@ -515,20 +525,22 @@ TEST(StanMath_integrate_1d_rev, TestDoubleExponential) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestExponential) { +TEST(StanMath_integrate_1d_impl_rev, TestExponential) { using stan::math::exp; - using stan::math::integrate_1d; + 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; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -537,9 +549,9 @@ TEST(StanMath_integrate_1d_rev, TestExponential) { EXPECT_FLOAT_EQ(1, 1 + g[0]); } -TEST(StanMath_integrate_1d_rev, TestFrechet) { +TEST(StanMath_integrate_1d_impl_rev, TestFrechet) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var alpha = 9.0 / 5; @@ -547,11 +559,13 @@ TEST(StanMath_integrate_1d_rev, TestFrechet) { AVEC theta = {alpha, sigma}; double b = std::numeric_limits::infinity(); double a = 0; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -561,9 +575,9 @@ TEST(StanMath_integrate_1d_rev, TestFrechet) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestGamma) { +TEST(StanMath_integrate_1d_impl_rev, TestGamma) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var alpha = 9.0 / 5; @@ -571,11 +585,13 @@ TEST(StanMath_integrate_1d_rev, TestGamma) { AVEC theta = {alpha, beta}; double b = std::numeric_limits::infinity(); double a = 0; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -585,9 +601,9 @@ TEST(StanMath_integrate_1d_rev, TestGamma) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestGumbel) { +TEST(StanMath_integrate_1d_impl_rev, TestGumbel) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var mu = 9.0 / 5; @@ -595,11 +611,13 @@ TEST(StanMath_integrate_1d_rev, TestGumbel) { AVEC theta = {mu, beta}; double b = std::numeric_limits::infinity(); double a = -b; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -609,20 +627,22 @@ TEST(StanMath_integrate_1d_rev, TestGumbel) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestInvChiSquared) { +TEST(StanMath_integrate_1d_impl_rev, TestInvChiSquared) { using stan::math::exp; - using stan::math::integrate_1d; + 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; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -631,9 +651,9 @@ TEST(StanMath_integrate_1d_rev, TestInvChiSquared) { EXPECT_FLOAT_EQ(1, 1 + g[0]); } -TEST(StanMath_integrate_1d_rev, TestLogistic) { +TEST(StanMath_integrate_1d_impl_rev, TestLogistic) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var mu = 9.0 / 5; @@ -641,11 +661,13 @@ TEST(StanMath_integrate_1d_rev, TestLogistic) { AVEC theta = {mu, sigma}; double b = std::numeric_limits::infinity(); double a = -b; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -655,9 +677,9 @@ TEST(StanMath_integrate_1d_rev, TestLogistic) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestLogNormal) { +TEST(StanMath_integrate_1d_impl_rev, TestLogNormal) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var mu = 9.0 / 5; @@ -665,11 +687,13 @@ TEST(StanMath_integrate_1d_rev, TestLogNormal) { AVEC theta = {mu, sigma}; double b = std::numeric_limits::infinity(); double a = 0; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -679,9 +703,9 @@ TEST(StanMath_integrate_1d_rev, TestLogNormal) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestNormal) { +TEST(StanMath_integrate_1d_impl_rev, TestNormal) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var mu = 9.0 / 5; @@ -689,11 +713,13 @@ TEST(StanMath_integrate_1d_rev, TestNormal) { AVEC theta = {mu, sigma}; double b = std::numeric_limits::infinity(); double a = -b; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -703,9 +729,9 @@ TEST(StanMath_integrate_1d_rev, TestNormal) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestPareto) { +TEST(StanMath_integrate_1d_impl_rev, TestPareto) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var m = 9.0 / 5; @@ -713,11 +739,13 @@ TEST(StanMath_integrate_1d_rev, TestPareto) { AVEC theta = {m, alpha}; double b = std::numeric_limits::infinity(); var a = m; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -727,9 +755,9 @@ TEST(StanMath_integrate_1d_rev, TestPareto) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestPareto2) { +TEST(StanMath_integrate_1d_impl_rev, TestPareto2) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var mu = 9.0 / 5; @@ -738,11 +766,13 @@ TEST(StanMath_integrate_1d_rev, TestPareto2) { AVEC theta = {mu, lambda, alpha}; double b = std::numeric_limits::infinity(); var a = mu; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -753,20 +783,22 @@ TEST(StanMath_integrate_1d_rev, TestPareto2) { EXPECT_FLOAT_EQ(1, 1 + g[2]); } -TEST(StanMath_integrate_1d_rev, TestRayleigh) { +TEST(StanMath_integrate_1d_impl_rev, TestRayleigh) { using stan::math::exp; - using stan::math::integrate_1d; + 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; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -775,9 +807,9 @@ TEST(StanMath_integrate_1d_rev, TestRayleigh) { EXPECT_FLOAT_EQ(1, 1 + g[0]); } -TEST(StanMath_integrate_1d_rev, TestScaledInvChiSquare) { +TEST(StanMath_integrate_1d_impl_rev, TestScaledInvChiSquare) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var nu = 9.0 / 5; @@ -785,11 +817,13 @@ TEST(StanMath_integrate_1d_rev, TestScaledInvChiSquare) { AVEC theta = {nu, s}; double b = std::numeric_limits::infinity(); double a = 0; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -799,9 +833,9 @@ TEST(StanMath_integrate_1d_rev, TestScaledInvChiSquare) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestStudentT) { +TEST(StanMath_integrate_1d_impl_rev, TestStudentT) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var nu = 11.0 / 3; @@ -810,11 +844,13 @@ TEST(StanMath_integrate_1d_rev, TestStudentT) { AVEC theta = {nu, mu, sigma}; double b = std::numeric_limits::infinity(); double a = -b; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -825,19 +861,21 @@ TEST(StanMath_integrate_1d_rev, TestStudentT) { EXPECT_FLOAT_EQ(1, 1 + g[2]); } -TEST(StanMath_integrate_1d_rev, TestUniform) { +TEST(StanMath_integrate_1d_impl_rev, TestUniform) { using stan::math::exp; - using stan::math::integrate_1d; + 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, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -847,9 +885,9 @@ TEST(StanMath_integrate_1d_rev, TestUniform) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestVonMises) { +TEST(StanMath_integrate_1d_impl_rev, TestVonMises) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var mu = 9.0 / 5; @@ -857,11 +895,13 @@ TEST(StanMath_integrate_1d_rev, TestVonMises) { AVEC theta = {mu, kappa}; double b = stan::math::pi() * 2; double a = 0; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -871,9 +911,9 @@ TEST(StanMath_integrate_1d_rev, TestVonMises) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } -TEST(StanMath_integrate_1d_rev, TestWeibull) { +TEST(StanMath_integrate_1d_impl_rev, TestWeibull) { using stan::math::exp; - using stan::math::integrate_1d; + using stan::math::integrate_1d_impl; using stan::math::var; var alpha = 9.0 / 5; @@ -881,11 +921,13 @@ TEST(StanMath_integrate_1d_rev, TestWeibull) { AVEC theta = {alpha, sigma}; double b = std::numeric_limits::infinity(); double a = 0; - auto pdf = [](auto x, auto xc, auto theta, auto x_r, auto x_i, - std::ostream *msgs) { + 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(pdf, a, b, theta, {}, {}, msgs, 1e-8); + 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); @@ -894,4 +936,3 @@ TEST(StanMath_integrate_1d_rev, TestWeibull) { EXPECT_FLOAT_EQ(1, 1 + g[0]); EXPECT_FLOAT_EQ(1, 1 + g[1]); } -*/ From 132cd2c2f0e7a44e3873a8257ec1f9854ce9b4f6 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 27 Feb 2021 13:58:49 +0000 Subject: [PATCH 06/20] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/functor/integrate_1d.hpp | 18 +- .../prim/functor/integrate_1d_adapter.hpp | 9 +- stan/math/rev/functor/integrate_1d.hpp | 102 +++---- .../prim/functor/integrate_1d_impl_test.cpp | 257 ++++++++---------- .../rev/functor/integrate_1d_impl_test.cpp | 191 ++++++------- 5 files changed, 279 insertions(+), 298 deletions(-) diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index 31c428c54c8..14cbf015be8 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -146,10 +146,10 @@ inline double integrate(const F& f, double a, double b, * @return numeric integral of function f */ template * = nullptr> + require_all_not_st_var* = nullptr> inline double integrate_1d_impl(const F& f, double a, double b, - double relative_tolerance, - std::ostream* msgs, const Args&... args) { + double relative_tolerance, std::ostream* msgs, + const Args&... args) { static const char* function = "integrate_1d"; check_less_or_equal(function, "lower limit", a, b); @@ -159,10 +159,9 @@ inline double integrate_1d_impl(const F& f, double a, double b, } return 0.0; } else { - return integrate( - std::bind(f, std::placeholders::_1, std::placeholders::_2, - msgs, args...), - a, b, relative_tolerance); + return integrate(std::bind(f, std::placeholders::_1, + std::placeholders::_2, msgs, args...), + a, b, relative_tolerance); } } @@ -219,9 +218,8 @@ inline double integrate_1d(const F& f, double a, double b, const std::vector& x_i, std::ostream* msgs, const double relative_tolerance = std::sqrt(EPSILON)) { - return integrate_1d_impl(integrate_1d_adapter(f), a, b, - relative_tolerance, msgs, - theta, x_r, x_i); + 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 index 82033a70823..ecfcaaaa9a6 100644 --- a/stan/math/prim/functor/integrate_1d_adapter.hpp +++ b/stan/math/prim/functor/integrate_1d_adapter.hpp @@ -17,11 +17,10 @@ struct integrate_1d_adapter { 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 { + 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); } }; diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index 60cc3eb4b9f..e0568ea8ed3 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -38,28 +38,30 @@ namespace math { */ template inline double gradient_of_f(const F &f, const double &x, const double &xc, - size_t n, std::ostream *msgs, const Args&... args) { + size_t n, std::ostream *msgs, + const Args &... args) { double gradient = 0.0; // Run nested autodiff in this scope nested_rev_autodiff nested; - std::tuple - args_tuple_local_copy(deep_copy_vars(args)...); + std::tuple args_tuple_local_copy( + deep_copy_vars(args)...); Eigen::VectorXd adjoints = Eigen::VectorXd::Zero(count_vars(args...)); - var fx = apply([&f, &x, &xc, msgs](auto&&... args) { - return f(x, xc, msgs, args...); - }, args_tuple_local_copy); + var fx = apply( + [&f, &x, &xc, msgs](auto &&... args) { return f(x, xc, msgs, args...); }, + args_tuple_local_copy); fx.grad(); - apply([&](auto&&... args) { - accumulate_adjoints(adjoints.data(), - std::forward(args)...); - }, - std::move(args_tuple_local_copy)); + apply( + [&](auto &&... args) { + accumulate_adjoints(adjoints.data(), + std::forward(args)...); + }, + std::move(args_tuple_local_copy)); gradient = adjoints.coeff(n); if (is_nan(gradient)) { @@ -90,13 +92,11 @@ inline double gradient_of_f(const F &f, const double &x, const double &xc, * @param args additional arguments to pass to f * @return numeric integral of function f */ -template * = nullptr> +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) { + const F &f, const T_a &a, const T_b &b, double relative_tolerance, + std::ostream *msgs, const Args &... args) { static const char *function = "integrate_1d"; check_less_or_equal(function, "lower limit", a, b); @@ -105,65 +105,66 @@ inline return_type_t integrate_1d_impl( if (a_val == b_val) { if (is_inf(a_val)) { - throw_domain_error(function, "Integration endpoints are both", - a_val, "", ""); + throw_domain_error(function, "Integration endpoints are both", a_val, "", + ""); } return var(0.0); } else { - std::tuple - args_val_tuple(value_of(args)...); + std::tuple args_val_tuple(value_of(args)...); - double integral = integrate(apply([&f, msgs](auto&&... args) { - return std::bind(f, std::placeholders::_1, std::placeholders::_2, - msgs, args...); - }, args_val_tuple), - a_val, b_val, relative_tolerance); + double integral = integrate(apply( + [&f, msgs](auto &&... args) { + return std::bind( + f, std::placeholders::_1, + std::placeholders::_2, msgs, args...); + }, + args_val_tuple), + a_val, b_val, relative_tolerance); size_t num_vars_ab = count_vars(a, b); 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); - double* partials_ptr = partials; + 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); + double *partials_ptr = partials; save_varis(varis, a, b, args...); - for(size_t i = 0; i < num_vars_ab + num_vars_args; ++i) { + for (size_t i = 0; i < num_vars_ab + num_vars_args; ++i) { partials[i] = 0.0; } if (!is_inf(a) && is_var::value) { - *partials_ptr = apply([&f, a_val, msgs](auto&&... args) { - return -f(a_val, 0.0, msgs, args...); - }, args_val_tuple); + *partials_ptr = apply( + [&f, a_val, msgs](auto &&... args) { + return -f(a_val, 0.0, msgs, args...); + }, + args_val_tuple); partials_ptr++; } if (!is_inf(b) && is_var::value) { - *partials_ptr = apply([&f, b_val, msgs](auto&&... args) { - return f(b_val, 0.0, msgs, args...); - }, args_val_tuple); + *partials_ptr + = apply([&f, b_val, msgs]( + auto &&... args) { return f(b_val, 0.0, msgs, args...); }, + args_val_tuple); partials_ptr++; } for (size_t n = 0; n < num_vars_args; ++n) { *partials_ptr = integrate( - std::bind(gradient_of_f, f, - std::placeholders::_1, std::placeholders::_2, - n, msgs, args...), - a_val, b_val, relative_tolerance); + std::bind(gradient_of_f, f, std::placeholders::_1, + std::placeholders::_2, n, msgs, args...), + a_val, b_val, relative_tolerance); partials_ptr++; } - return var(new precomputed_gradients_vari(integral, - num_vars_ab + num_vars_args, - varis, partials)); + return var(new precomputed_gradients_vari( + integral, num_vars_ab + num_vars_args, varis, partials)); } } - + /** * 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. @@ -225,9 +226,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)) { - return integrate_1d_impl(integrate_1d_adapter(f), a, b, - relative_tolerance, msgs, - theta, x_r, x_i); + 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 index 9f20b385744..985a4cd6392 100644 --- a/test/unit/math/prim/functor/integrate_1d_impl_test.cpp +++ b/test/unit/math/prim/functor/integrate_1d_impl_test.cpp @@ -12,16 +12,14 @@ std::ostringstream *msgs = nullptr; struct f1 { template - inline T1 operator()(const T1 &x, const T1 &xc, - std::ostream *msgs) const { + 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 { + 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 { @@ -32,92 +30,84 @@ struct f2 { struct f3 { template - inline T1 operator()(const T1 &x, const T1 &xc, - std::ostream *msgs) const { + 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 { + 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 { + 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 { + 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 { + 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 { + 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 { + 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 { + 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); + * 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 { + 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 { + 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; @@ -127,8 +117,7 @@ struct f12 { struct f13 { template - inline T1 operator()(const T1 &x, const T1 &xc, - std::ostream *msgs) const { + 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); @@ -138,16 +127,14 @@ struct f13 { struct f14 { template - inline T1 operator()(const T1 &x, const T1 &xc, - std::ostream *msgs) const { + 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 { + 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; @@ -158,8 +145,7 @@ struct f15 { struct f16 { template - inline T1 operator()(const T1 &x, const T1 &xc, - std::ostream *msgs) const { + inline T1 operator()(const T1 &x, const T1 &xc, std::ostream *msgs) const { return x * sin(x) / (1 + stan::math::square(cos(x))); } }; @@ -220,11 +206,10 @@ double rank_density(double x, double xc, const std::vector &theta, } struct rank_density_functor__ { - double operator()(double x, double xc, - std::ostream *pstream__, - const std::vector &theta, + double operator()(double x, double xc, std::ostream *pstream__, + const std::vector &theta, const std::vector &x_r, - const std::vector &x_i) const { + const std::vector &x_i) const { return rank_density(x, xc, theta, x_r, x_i, pstream__); } }; @@ -235,10 +220,8 @@ double order(double down, double up, const std::vector &theta, double v; - v = stan::math::integrate_1d_impl(rank_density_functor__(), - down, up, - 1e-8, pstream__, - theta, x_r, x_i); + 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 @@ -269,26 +252,25 @@ double order(double down, double up, const std::vector &theta, * @param val correct value of integral */ template -void test_integration(const F &f, double a, double b, - double val, const Args&... args) { +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...) + 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...); - }; + 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...) + integrate_1d_impl_test::msgs, args...) - val), tolerance); } @@ -296,74 +278,81 @@ void test_integration(const F &f, double a, double b, 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); + 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::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), + 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), + 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), + 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); + 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); + 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)); + 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)); + 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); + 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 @@ -372,53 +361,47 @@ TEST(StanMath_integrate_1d_impl_prim, test1) { 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.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})); + std::vector({0.5})); test_integration(integrate_1d_impl_test::f5{}, 1.0, 1.0, 0.0, - std::vector({0.4, 0.4})); + 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_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})); + 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})); + 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})); + 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})); + 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); + 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)); @@ -428,7 +411,7 @@ TEST(StanMath_integrate_1d_impl_prim, test1) { // 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); + stan::math::square(stan::math::pi()) / 4); } TEST(StanMath_integrate_1d_impl_prim, TestTolerance) { diff --git a/test/unit/math/rev/functor/integrate_1d_impl_test.cpp b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp index 6b46a7d428a..8244a2638ad 100644 --- a/test/unit/math/rev/functor/integrate_1d_impl_test.cpp +++ b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp @@ -12,9 +12,9 @@ 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 { + 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]; } }; @@ -22,9 +22,9 @@ struct f1 { 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 { + 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]; } }; @@ -32,9 +32,9 @@ struct f2 { 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 { + 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]; } @@ -43,9 +43,9 @@ struct f3 { 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 { + 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); } }; @@ -53,9 +53,9 @@ struct f4 { 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 { + 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); } }; @@ -63,9 +63,9 @@ struct f5 { 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 { + 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)); } }; @@ -73,9 +73,9 @@ struct f6 { 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 { + 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); } }; @@ -83,9 +83,9 @@ struct f7 { 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 { + 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); } }; @@ -93,9 +93,9 @@ struct f8 { 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 { + 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]); } }; @@ -103,9 +103,9 @@ struct f10 { 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 { + 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); } @@ -114,9 +114,9 @@ struct f11 { 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 { + 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)) @@ -127,9 +127,9 @@ struct f12 { 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 { + 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]; } }; @@ -216,7 +216,8 @@ void test_derivatives(const F &f, double a, double b, 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); + 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) { @@ -238,12 +239,12 @@ TEST(StanMath_integrate_1d_impl_rev, test_integer_arguments) { 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)); + 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) { @@ -260,7 +261,7 @@ TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_easy) { 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}); + {0.0}); } TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_zero_crossing) { @@ -274,7 +275,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_zero_crossing) { -19.06340613646808, 21.41380852375568); } -TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_var_right_endpoint_var_params) { +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( @@ -284,7 +286,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_var_right_endpoint_var_para {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) { +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( @@ -405,7 +408,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestDerivatives_gaussian) { {0.0, 0.0}); } -TEST(StanMath_integrate_1d_impl_rev, TestDerivativesSameVarAtEndpointAndInParams) { +TEST(StanMath_integrate_1d_impl_rev, + TestDerivativesSameVarAtEndpointAndInParams) { using stan::math::var; var a = 2.0; @@ -414,8 +418,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestDerivativesSameVarAtEndpointAndInParams 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); + 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); @@ -433,8 +437,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestBeta) { 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) { + 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); @@ -459,8 +463,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestCauchy) { 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) { + 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); @@ -484,10 +488,9 @@ TEST(StanMath_integrate_1d_impl_rev, TestChiSquare) { 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])); - }; + 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()); @@ -509,8 +512,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestDoubleExponential) { 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) { + 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 @@ -536,10 +539,9 @@ TEST(StanMath_integrate_1d_impl_rev, TestExponential) { 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])); - }; + 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()); @@ -561,8 +563,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestFrechet) { 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) { + 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); @@ -587,8 +589,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestGamma) { 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) { + 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); @@ -613,8 +615,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestGumbel) { 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) { + 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); @@ -638,8 +640,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestInvChiSquared) { 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) { + 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); @@ -663,8 +665,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestLogistic) { 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) { + 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); @@ -689,8 +691,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestLogNormal) { 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) { + 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); @@ -715,8 +717,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestNormal) { 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) { + 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); @@ -741,8 +743,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestPareto) { 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) { + 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); @@ -768,8 +770,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestPareto2) { 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) { + 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); @@ -794,10 +796,9 @@ TEST(StanMath_integrate_1d_impl_rev, TestRayleigh) { 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])); - }; + 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()); @@ -819,8 +820,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestScaledInvChiSquare) { 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) { + 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); @@ -846,8 +847,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestStudentT) { 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) { + 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); @@ -871,8 +872,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestUniform) { 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) { + 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); @@ -897,8 +898,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestVonMises) { 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) { + 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); @@ -923,8 +924,8 @@ TEST(StanMath_integrate_1d_impl_rev, TestWeibull) { 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) { + 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); From 59d1099b05f65c1de73334808344b4200cd9b4f2 Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 25 Mar 2021 15:36:52 -0400 Subject: [PATCH 07/20] Switch binds to lambdas --- stan/math/prim/functor/integrate_1d.hpp | 7 ++++--- stan/math/rev/functor/integrate_1d.hpp | 20 ++++++++------------ 2 files changed, 12 insertions(+), 15 deletions(-) diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index 14cbf015be8..e520e798bb1 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -58,6 +58,7 @@ inline double integrate(const F& f, double a, double b, bool used_two_integrals = false; size_t levels; double Q = 0.0; + // 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 integrator; @@ -159,9 +160,9 @@ inline double integrate_1d_impl(const F& f, double a, double b, } return 0.0; } else { - return integrate(std::bind(f, std::placeholders::_1, - std::placeholders::_2, msgs, args...), - a, b, relative_tolerance); + return integrate([&](const auto& x, const auto& xc) { + return f(x, xc, msgs, args...); + }, a, b, relative_tolerance); } } diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index e0568ea8ed3..663e487c578 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -112,14 +112,11 @@ inline return_type_t integrate_1d_impl( } else { std::tuple args_val_tuple(value_of(args)...); - double integral = integrate(apply( - [&f, msgs](auto &&... args) { - return std::bind( - f, std::placeholders::_1, - std::placeholders::_2, msgs, args...); - }, - args_val_tuple), - a_val, b_val, relative_tolerance); + double integral = integrate([&](const auto& x, const auto& xc) { + return apply([&](auto &&... args) { + return f(x, xc, msgs, args...); + }, args_val_tuple); + }, a_val, b_val, relative_tolerance); size_t num_vars_ab = count_vars(a, b); size_t num_vars_args = count_vars(args...); @@ -153,10 +150,9 @@ inline return_type_t integrate_1d_impl( } for (size_t n = 0; n < num_vars_args; ++n) { - *partials_ptr = integrate( - std::bind(gradient_of_f, f, std::placeholders::_1, - std::placeholders::_2, n, msgs, args...), - a_val, b_val, relative_tolerance); + *partials_ptr = integrate([&](const auto& x, const auto& xc) { + return gradient_of_f(f, x, xc, n, msgs, args...); + }, a_val, b_val, relative_tolerance); partials_ptr++; } From a6cebfc94e8b17cc38fe8c166753ca755f65c305 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 25 Mar 2021 19:41:25 +0000 Subject: [PATCH 08/20] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/functor/integrate_1d.hpp | 9 +++++---- stan/math/rev/functor/integrate_1d.hpp | 19 +++++++++++-------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index e520e798bb1..e0fb50c1687 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -58,7 +58,8 @@ inline double integrate(const F& f, double a, double b, bool used_two_integrals = false; size_t levels; double Q = 0.0; - // if a or b is infinite, set xc argument to NaN (see docs above for user function for xc info) + // 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 integrator; @@ -160,9 +161,9 @@ inline double integrate_1d_impl(const F& f, double a, double b, } return 0.0; } else { - return integrate([&](const auto& x, const auto& xc) { - return f(x, xc, msgs, args...); - }, a, b, relative_tolerance); + return integrate( + [&](const auto& x, const auto& xc) { return f(x, xc, msgs, args...); }, + a, b, relative_tolerance); } } diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index 663e487c578..be355d4001d 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -112,11 +112,12 @@ inline return_type_t integrate_1d_impl( } else { std::tuple args_val_tuple(value_of(args)...); - double integral = integrate([&](const auto& x, const auto& xc) { - return apply([&](auto &&... args) { - return f(x, xc, msgs, args...); - }, args_val_tuple); - }, a_val, b_val, relative_tolerance); + double integral = integrate( + [&](const auto &x, const auto &xc) { + return apply([&](auto &&... args) { return f(x, xc, msgs, args...); }, + args_val_tuple); + }, + a_val, b_val, relative_tolerance); size_t num_vars_ab = count_vars(a, b); size_t num_vars_args = count_vars(args...); @@ -150,9 +151,11 @@ inline return_type_t integrate_1d_impl( } for (size_t n = 0; n < num_vars_args; ++n) { - *partials_ptr = integrate([&](const auto& x, const auto& xc) { - return gradient_of_f(f, x, xc, n, msgs, args...); - }, a_val, b_val, relative_tolerance); + *partials_ptr = integrate( + [&](const auto &x, const auto &xc) { + return gradient_of_f(f, x, xc, n, msgs, args...); + }, + a_val, b_val, relative_tolerance); partials_ptr++; } From abaa60298ab000e0ddf1ebdc2b7c5a19cf5f85f1 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Thu, 25 Mar 2021 19:52:37 -0400 Subject: [PATCH 09/20] use double nested reverse pass to save making N tuple copies --- stan/math/rev/functor/integrate_1d.hpp | 58 +++++++++++++------------- 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index be355d4001d..f0bf3a3e4c7 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -39,30 +39,21 @@ namespace math { template inline double gradient_of_f(const F &f, const double &x, const double &xc, size_t n, std::ostream *msgs, - const Args &... args) { + const std::tuple &args_tuple, + size_t num_var_args) { double gradient = 0.0; // Run nested autodiff in this scope + Eigen::VectorXd adjoints = Eigen::VectorXd::Zero(num_var_args); nested_rev_autodiff nested; - - std::tuple args_tuple_local_copy( - deep_copy_vars(args)...); - - Eigen::VectorXd adjoints = Eigen::VectorXd::Zero(count_vars(args...)); - var fx = apply( [&f, &x, &xc, msgs](auto &&... args) { return f(x, xc, msgs, args...); }, - args_tuple_local_copy); + args_tuple); fx.grad(); - apply( - [&](auto &&... args) { - accumulate_adjoints(adjoints.data(), - std::forward(args)...); - }, - std::move(args_tuple_local_copy)); - + apply([&](auto &&... args) { accumulate_adjoints(adjoints.data(), args...); }, + args_tuple); gradient = adjoints.coeff(n); if (is_nan(gradient)) { if (fx.val() == 0) { @@ -97,7 +88,7 @@ template integrate_1d_impl( const F &f, const T_a &a, const T_b &b, double relative_tolerance, std::ostream *msgs, const Args &... args) { - static const char *function = "integrate_1d"; + static constexpr const char* function = "integrate_1d"; check_less_or_equal(function, "lower limit", a, b); double a_val = value_of(a); @@ -110,7 +101,7 @@ inline return_type_t integrate_1d_impl( } return var(0.0); } else { - std::tuple args_val_tuple(value_of(args)...); + auto args_val_tuple = std::make_tuple(value_of(args)...); double integral = integrate( [&](const auto &x, const auto &xc) { @@ -119,7 +110,7 @@ inline return_type_t integrate_1d_impl( }, a_val, b_val, relative_tolerance); - size_t num_vars_ab = count_vars(a, b); + 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); @@ -149,18 +140,29 @@ inline return_type_t integrate_1d_impl( args_val_tuple); partials_ptr++; } - - for (size_t n = 0; n < num_vars_args; ++n) { - *partials_ptr = integrate( - [&](const auto &x, const auto &xc) { - return gradient_of_f(f, x, xc, n, msgs, args...); - }, - a_val, b_val, relative_tolerance); - partials_ptr++; + { + nested_rev_autodiff nested; + auto args_tuple_local_copy = std::make_tuple(deep_copy_vars(args)...); + for (size_t n = 0; n < num_vars_args; ++n) { + *partials_ptr = integrate( + [&](const auto &x, const auto &xc) { + auto blah = gradient_of_f(f, x, xc, n, msgs, + args_tuple_local_copy, num_vars_args); + nested.set_zero_all_adjoints(); + return blah; + }, + a_val, b_val, relative_tolerance); + partials_ptr++; + } } - return var(new precomputed_gradients_vari( - integral, num_vars_ab + num_vars_args, varis, partials)); + 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(); + } + }); } } From 82624e8dff82a390acb28a1b399b1e5feef306ee Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Mon, 29 Mar 2021 10:30:53 -0400 Subject: [PATCH 10/20] Update stan/math/rev/functor/integrate_1d.hpp Co-authored-by: Steve Bronder --- stan/math/rev/functor/integrate_1d.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index be355d4001d..588cc4bbc25 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -45,8 +45,7 @@ inline double gradient_of_f(const F &f, const double &x, const double &xc, // Run nested autodiff in this scope nested_rev_autodiff nested; - std::tuple args_tuple_local_copy( - deep_copy_vars(args)...); + auto args_tuple_local_copy = std::make_tuple(deep_copy_vars(args)...); Eigen::VectorXd adjoints = Eigen::VectorXd::Zero(count_vars(args...)); From 3469fa1e96ffb90076c6250d1ebd6a6c1263203e Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Mon, 29 Mar 2021 10:31:08 -0400 Subject: [PATCH 11/20] Update stan/math/prim/functor/integrate_1d.hpp Co-authored-by: Steve Bronder --- stan/math/prim/functor/integrate_1d.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index e0fb50c1687..3f63c8edca0 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -152,7 +152,7 @@ template Date: Mon, 29 Mar 2021 10:42:15 -0400 Subject: [PATCH 12/20] Updated docs --- stan/math/rev/functor/integrate_1d.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index 588cc4bbc25..3ab5df24d86 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -20,8 +20,8 @@ namespace stan { namespace math { /** - * Calculate first derivative of f(x, param, std::ostream&) - * with respect to the nth parameter. Uses nested reverse mode autodiff + * Calculate first derivative of f(x, xc, msgs, args...) + * with respect to the nth parameter in args. Uses nested reverse mode autodiff * * 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 @@ -78,10 +78,10 @@ inline double gradient_of_f(const F &f, const double &x, const double &xc, /** * Return the integral of f from a to b to the given relative tolerance * + * @tparam F Type of f * @tparam T_a type of first limit * @tparam T_b type of second limit - * @tparam T_theta type of parameters - * @tparam T Type of f + * @tparam Args types of parameter pack arguments * * @param f the functor to integrate * @param a lower limit of integration From 95b4032071c8080850c86089d1a0b9ad0647d4a8 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Mon, 29 Mar 2021 10:42:52 -0400 Subject: [PATCH 13/20] Update stan/math/rev/functor/integrate_1d.hpp Co-authored-by: Steve Bronder --- stan/math/rev/functor/integrate_1d.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index 588cc4bbc25..f6cf0993d06 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -96,7 +96,7 @@ template integrate_1d_impl( const F &f, const T_a &a, const T_b &b, double relative_tolerance, std::ostream *msgs, const Args &... args) { - static const char *function = "integrate_1d"; + static constexpr const char *function = "integrate_1d"; check_less_or_equal(function, "lower limit", a, b); double a_val = value_of(a); From c4540c572e7678a3485960456d67961c68f25513 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 29 Mar 2021 10:46:08 -0400 Subject: [PATCH 14/20] Reordered if --- stan/math/rev/functor/integrate_1d.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index 71a3967aece..f9bad1d8a15 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -132,7 +132,7 @@ inline return_type_t integrate_1d_impl( partials[i] = 0.0; } - if (!is_inf(a) && is_var::value) { + if (is_var::value && !is_inf(a)) { *partials_ptr = apply( [&f, a_val, msgs](auto &&... args) { return -f(a_val, 0.0, msgs, args...); From aba98b607dccfba47636e677295c991549ff2e54 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 29 Mar 2021 11:34:46 -0400 Subject: [PATCH 15/20] Put error checks into functions --- stan/math/prim/functor/integrate_1d.hpp | 57 ++++++++++++++----------- 1 file changed, 32 insertions(+), 25 deletions(-) diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index 3f63c8edca0..85491da47b2 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -51,19 +51,44 @@ 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 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 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 integrator; /** @@ -74,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 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 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 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); }; @@ -103,34 +130,14 @@ 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"); - } - } 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 Q; } From 77164ccf4dd275f96a9d58d47bb594be231eb1f3 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 29 Mar 2021 16:57:50 +0000 Subject: [PATCH 16/20] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/functor/integrate_1d.hpp | 14 +++++++------- stan/math/rev/functor/integrate_1d.hpp | 23 ++++++++++++++--------- 2 files changed, 21 insertions(+), 16 deletions(-) diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index 85491da47b2..09c58753e59 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -70,18 +70,18 @@ inline double integrate(const F& f, double a, double b, 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"); + 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"); + 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); }; diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index 96e3a06154d..0701a161ab1 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -99,8 +99,8 @@ inline return_type_t integrate_1d_impl( // do it once in a separate nest for efficiency auto args_tuple_local_copy = std::make_tuple(deep_copy_vars(args)...); for (size_t n = 0; n < num_vars_args; ++n) { - - // This computes the integral of the gradient of f with respect to the nth + // This computes the integral of the gradient of f with respect to the + // nth // parameter in args. Uses nested reverse mode autodiff *partials_ptr = integrate( [&](const auto &x, const auto &xc) { @@ -111,25 +111,30 @@ inline return_type_t integrate_1d_impl( nested_rev_autodiff gradient_nest; var fx = apply( - [&f, &x, &xc, msgs](auto &&... args) { return f(x, xc, msgs, args...); }, + [&f, &x, &xc, msgs](auto &&... args) { + return f(x, xc, msgs, args...); + }, args_tuple_local_copy); fx.grad(); - apply([&](auto &&... args) { accumulate_adjoints(adjoints.data(), args...); }, - args_tuple_local_copy); + apply( + [&](auto &&... args) { + accumulate_adjoints(adjoints.data(), args...); + }, + args_tuple_local_copy); double gradient = adjoints.coeff(n); - // 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 + // 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 ", ""); + "is nan for parameter ", ""); } } From 4e882f21ff895c9c3df98434942ae3172212bc20 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Mon, 29 Mar 2021 19:27:40 -0400 Subject: [PATCH 17/20] small cleanups and use get the adjoint in each integration by a lookup instead of copying all the adjoints --- stan/math/prim/fun/get.hpp | 24 ++++---- stan/math/prim/functor/integrate_1d.hpp | 77 ++++++++++++++----------- stan/math/rev/functor/integrate_1d.hpp | 64 +++++++++++--------- 3 files changed, 93 insertions(+), 72 deletions(-) diff --git a/stan/math/prim/fun/get.hpp b/stan/math/prim/fun/get.hpp index 8a466561278..0552df50c64 100644 --- a/stan/math/prim/fun/get.hpp +++ b/stan/math/prim/fun/get.hpp @@ -25,27 +25,27 @@ inline T get(const T& x, size_t n) { } /** \ingroup type_trait - * Returns the n-th element of the provided std::vector. + * Returns the n-th element of the provided Eigen matrix. * - * @param x input vector + * @param m input \c Eigen \c Matrix or expression * @param n index of the element to return - * @return n-th element of the input vector + * @return n-th element of the \c Eigen \c Matrix or expression */ -template -inline T get(const std::vector& x, size_t n) { - return x[n]; +template > +inline scalar_type_t get(const T& m, size_t n) { + return m(static_cast(n)); } /** \ingroup type_trait - * Returns the n-th element of the provided Eigen matrix. + * Returns the n-th element of the provided std::vector. * - * @param m input \c Eigen \c Matrix or expression + * @param x input vector * @param n index of the element to return - * @return n-th element of the \c Eigen \c Matrix or expression + * @return n-th element of the input vector */ -template > -inline scalar_type_t get(const T& m, size_t n) { - return m(static_cast(n)); +template +inline auto get(const std::vector& x, size_t n) { + return get(x[n], n - x.size()); } } // namespace stan diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index 09c58753e59..f3a00a95503 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -57,38 +57,44 @@ inline double integrate(const F& f, double a, double b, double L1 = 0.0; double L2 = 0.0; size_t levels; - double Q = 0.0; - auto one_integral_convergence_check = [&]() { - if (error1 > relative_tolerance * L1) { + 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 = [&]() { - 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"); + 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 > 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 (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 = [&](double x) { return f(x, NOT_A_NUMBER); }; + 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); - one_integral_convergence_check(); + 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; /** @@ -97,48 +103,52 @@ 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, + double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1, &levels); - one_integral_convergence_check(); + 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, + 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(); + 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, + double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1, &levels); - one_integral_convergence_check(); + 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, + 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(); + 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, + 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(); + two_integral_convergence_check(error1, error2, relative_tolerance, L1, L2); + return Q; } else { - Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1, + double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1, &levels); - one_integral_convergence_check(); + one_integral_convergence_check(error1, relative_tolerance, L1); + return Q; } } - - return Q; } /** @@ -161,8 +171,7 @@ inline double integrate_1d_impl(const F& f, double a, double b, const Args&... args) { static constexpr const char* function = "integrate_1d"; check_less_or_equal(function, "lower limit", a, b); - - if (a == b) { + if (unlikely(a == b)) { if (std::isinf(a)) { throw_domain_error(function, "Integration endpoints are both", a, "", ""); } diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index 0701a161ab1..5d71b5c2250 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -46,7 +47,7 @@ inline return_type_t integrate_1d_impl( double a_val = value_of(a); double b_val = value_of(b); - if (a_val == b_val) { + if (unlikely(a_val == b_val)) { if (is_inf(a_val)) { throw_domain_error(function, "Integration endpoints are both", a_val, "", ""); @@ -57,8 +58,9 @@ inline return_type_t integrate_1d_impl( double integral = integrate( [&](const auto &x, const auto &xc) { - return apply([&](auto &&... args) { return f(x, xc, msgs, args...); }, - args_val_tuple); + return apply( + [&](auto &&... val_args) { return f(x, xc, msgs, val_args...); }, + args_val_tuple); }, a_val, b_val, relative_tolerance); @@ -68,6 +70,7 @@ inline return_type_t integrate_1d_impl( 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...); @@ -78,18 +81,19 @@ inline return_type_t integrate_1d_impl( if (is_var::value && !is_inf(a)) { *partials_ptr = apply( - [&f, a_val, msgs](auto &&... args) { - return -f(a_val, 0.0, msgs, args...); + [&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 &&... args) { return f(b_val, 0.0, msgs, args...); }, - args_val_tuple); + *partials_ptr = apply( + [&f, b_val, msgs](auto &&... val_args) { + return f(b_val, 0.0, msgs, val_args...); + }, + args_val_tuple); partials_ptr++; } @@ -100,32 +104,41 @@ inline return_type_t integrate_1d_impl( auto args_tuple_local_copy = std::make_tuple(deep_copy_vars(args)...); 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. Uses nested reverse mode autodiff + // 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(); - - Eigen::VectorXd adjoints = Eigen::VectorXd::Zero(num_vars_args); - nested_rev_autodiff gradient_nest; - var fx = apply( - [&f, &x, &xc, msgs](auto &&... args) { - return f(x, xc, msgs, args...); + [&f, &x, &xc, msgs](auto &&... local_args) { + return f(x, xc, msgs, local_args...); }, args_tuple_local_copy); - fx.grad(); - - apply( - [&](auto &&... args) { - accumulate_adjoints(adjoints.data(), args...); + size_t adjoint_count = 0; + double gradient = 0; + bool not_found = true; + // for_each is guaranteed to go off from left to right. + // So for var argument we count the number of previous vars + // until we go past n, then index into that argument to get + // the correct adjoint. + stan::math::for_each( + [&](auto &arg) { + using arg_t = decltype(arg); + using scalar_arg_t = scalar_type_t; + if (is_var::value) { + size_t var_count = count_vars(arg); + if (((adjoint_count + var_count) < n) && not_found) { + adjoint_count += var_count; + } else if (not_found) { + not_found = false; + gradient + = forward_as(stan::get(arg, n - adjoint_count)) + .adj(); + } + } }, args_tuple_local_copy); - - double gradient = adjoints.coeff(n); - // 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 @@ -137,7 +150,6 @@ inline return_type_t integrate_1d_impl( "is nan for parameter ", ""); } } - return gradient; }, a_val, b_val, relative_tolerance); From cd1b0cb29a9e27067dd7f09e4be116ed0fff234a Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 30 Mar 2021 13:57:07 -0700 Subject: [PATCH 18/20] Use varis to get nth gradient --- stan/math/rev/functor/integrate_1d.hpp | 35 ++++++++------------------ 1 file changed, 11 insertions(+), 24 deletions(-) diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index 5d71b5c2250..1197a247e21 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -102,12 +102,20 @@ inline return_type_t integrate_1d_impl( // 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) { @@ -115,30 +123,9 @@ inline return_type_t integrate_1d_impl( }, args_tuple_local_copy); fx.grad(); - size_t adjoint_count = 0; - double gradient = 0; - bool not_found = true; - // for_each is guaranteed to go off from left to right. - // So for var argument we count the number of previous vars - // until we go past n, then index into that argument to get - // the correct adjoint. - stan::math::for_each( - [&](auto &arg) { - using arg_t = decltype(arg); - using scalar_arg_t = scalar_type_t; - if (is_var::value) { - size_t var_count = count_vars(arg); - if (((adjoint_count + var_count) < n) && not_found) { - adjoint_count += var_count; - } else if (not_found) { - not_found = false; - gradient - = forward_as(stan::get(arg, n - adjoint_count)) - .adj(); - } - } - }, - args_tuple_local_copy); + + 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 From 93600842ce0b92457cf6b37b9b6daa3e49cd85aa Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 30 Mar 2021 21:54:20 +0000 Subject: [PATCH 19/20] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/functor/integrate_1d.hpp | 95 +++++++++++++------------ stan/math/rev/functor/integrate_1d.hpp | 10 +-- 2 files changed, 57 insertions(+), 48 deletions(-) diff --git a/stan/math/prim/functor/integrate_1d.hpp b/stan/math/prim/functor/integrate_1d.hpp index f3a00a95503..e4494c81ed8 100644 --- a/stan/math/prim/functor/integrate_1d.hpp +++ b/stan/math/prim/functor/integrate_1d.hpp @@ -58,41 +58,44 @@ inline double integrate(const F& f, double a, double b, double L2 = 0.0; size_t levels; - auto one_integral_convergence_check = [](auto& error1, auto& rel_tol, auto& L1) { + 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"); + 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"); - }(); - } - }; + 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; - double 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)) { @@ -103,48 +106,52 @@ 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) { - double 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; - 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); + 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) { - double 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; - 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); + 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 = [&f](double x, double xc) { return f(x, xc); }; boost::math::quadrature::tanh_sinh integrator; 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, L1, L2); + 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 { - double 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; } diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index 1197a247e21..ecd81887a42 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -104,10 +104,12 @@ inline return_type_t integrate_1d_impl( 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); + 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 From 49ba9553caaaf2aaaabba28d8bae10a48f64f1c3 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 30 Mar 2021 19:07:13 -0400 Subject: [PATCH 20/20] remove changes to math::get() now that we don't need it in integrate1d --- stan/math/prim/fun/get.hpp | 24 ++++++++++++------------ stan/math/rev/functor/integrate_1d.hpp | 1 - 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/stan/math/prim/fun/get.hpp b/stan/math/prim/fun/get.hpp index 0552df50c64..8a466561278 100644 --- a/stan/math/prim/fun/get.hpp +++ b/stan/math/prim/fun/get.hpp @@ -25,27 +25,27 @@ inline T get(const T& x, size_t n) { } /** \ingroup type_trait - * Returns the n-th element of the provided Eigen matrix. + * Returns the n-th element of the provided std::vector. * - * @param m input \c Eigen \c Matrix or expression + * @param x input vector * @param n index of the element to return - * @return n-th element of the \c Eigen \c Matrix or expression + * @return n-th element of the input vector */ -template > -inline scalar_type_t get(const T& m, size_t n) { - return m(static_cast(n)); +template +inline T get(const std::vector& x, size_t n) { + return x[n]; } /** \ingroup type_trait - * Returns the n-th element of the provided std::vector. + * Returns the n-th element of the provided Eigen matrix. * - * @param x input vector + * @param m input \c Eigen \c Matrix or expression * @param n index of the element to return - * @return n-th element of the input vector + * @return n-th element of the \c Eigen \c Matrix or expression */ -template -inline auto get(const std::vector& x, size_t n) { - return get(x[n], n - x.size()); +template > +inline scalar_type_t get(const T& m, size_t n) { + return m(static_cast(n)); } } // namespace stan diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp index ecd81887a42..be8ff106393 100644 --- a/stan/math/rev/functor/integrate_1d.hpp +++ b/stan/math/rev/functor/integrate_1d.hpp @@ -8,7 +8,6 @@ #include #include #include -#include #include #include #include