diff --git a/.gitattributes b/.gitattributes index db32a13f211..4cc2852556f 100644 --- a/.gitattributes +++ b/.gitattributes @@ -5,3 +5,7 @@ *.tex text diff=tex lib/** binary lib/eigen_5.0.1/STAN_CHANGES.md -binary text diff=markdown +lib/boost_1.87.0/STAN_CHANGES -binary text +# Stan-patched vendored source: keep it diffable so the patch stays reviewable +# and so a future Boost upgrade that drops it is visible rather than silent. +lib/boost_1.87.0/boost/math/quadrature/gauss_kronrod.hpp -binary text diff=cpp diff --git a/lib/boost_1.87.0/STAN_CHANGES b/lib/boost_1.87.0/STAN_CHANGES new file mode 100644 index 00000000000..dc3f902b487 --- /dev/null +++ b/lib/boost_1.87.0/STAN_CHANGES @@ -0,0 +1,9 @@ +This file documents changes done for the stan-math project + +- Added an `abs_tol` parameter to + `boost/math/quadrature/gauss_kronrod.hpp`'s public + `gauss_kronrod::integrate()`, and forwarded it to the five + `recursive_adaptive_integrate()` entry points in place of the hard-coded + `Real(0)`. + + Used by `stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp`. diff --git a/lib/boost_1.87.0/boost/math/quadrature/gauss_kronrod.hpp b/lib/boost_1.87.0/boost/math/quadrature/gauss_kronrod.hpp index a81cb43c4f5..f9852119393 100644 --- a/lib/boost_1.87.0/boost/math/quadrature/gauss_kronrod.hpp +++ b/lib/boost_1.87.0/boost/math/quadrature/gauss_kronrod.hpp @@ -1227,7 +1227,7 @@ class gauss_kronrod : public detail::gauss_kronrod_detail - static auto integrate(F f, Real a, Real b, unsigned max_depth = 15, Real tol = tools::root_epsilon(), Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval()(std::declval())) + static auto integrate(F f, Real a, Real b, unsigned max_depth = 15, Real tol = tools::root_epsilon(), Real* error = nullptr, Real* pL1 = nullptr, Real abs_tol = Real(0))->decltype(std::declval()(std::declval())) { typedef decltype(f(a)) K; static_assert(!std::is_integral::value, @@ -1248,7 +1248,7 @@ class gauss_kronrod : public detail::gauss_kronrod_detail info = { u, tol }; - K res = recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1); + K res = recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, abs_tol, error, pL1); return res; } @@ -1263,7 +1263,7 @@ class gauss_kronrod : public detail::gauss_kronrod_detail info = { u, tol }; - K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1); + K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, abs_tol, error, pL1); if (pL1) { *pL1 *= 2; @@ -1280,7 +1280,7 @@ class gauss_kronrod : public detail::gauss_kronrod_detail info = { v, tol }; - K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1); + K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, abs_tol, error, pL1); if (pL1) { *pL1 *= 2; @@ -1297,9 +1297,9 @@ class gauss_kronrod : public detail::gauss_kronrod_detail info = { f, tol }; if (b < a) { - return -recursive_adaptive_integrate(&info, b, a, max_depth, Real(0), error, pL1); + return -recursive_adaptive_integrate(&info, b, a, max_depth, abs_tol, error, pL1); } - return recursive_adaptive_integrate(&info, a, b, max_depth, Real(0), error, pL1); + return recursive_adaptive_integrate(&info, a, b, max_depth, abs_tol, error, pL1); } } return static_cast(policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy())); diff --git a/stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp b/stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp index 1ddd29cc8e7..c8c6705868a 100644 --- a/stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp +++ b/stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp @@ -12,7 +12,7 @@ namespace stan { namespace math { /** - * Return the integral of f from a to b using adaptive Gauss-Kronrod (G21,K21) + * Return the integral of f from a to b using adaptive Gauss-Kronrod (G10,K21) * quadrature, with tangents computed via finite differences over the * integrand parameters. * @@ -64,7 +64,7 @@ inline return_type_t integrate_1d_gauss_kronrod_tol( } /** - * Return the integral of f from a to b using adaptive Gauss-Kronrod (G21,K21) + * Return the integral of f from a to b using adaptive Gauss-Kronrod (G10,K21) * quadrature, with tangents computed via finite differences over the * integrand parameters. * diff --git a/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp b/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp index 0f0013f5ee6..3a3ef788bbd 100644 --- a/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp +++ b/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp @@ -29,7 +29,7 @@ constexpr int INTEGRATE_1D_GAUSS_KRONROD_MAX_DEPTH = 15; /** * Integrate a single variable function f from a to b using Boost's adaptive - * Gauss-Kronrod (G21,K21) quadrature, with QUADPACK-style mixed convergence + * Gauss-Kronrod (G10,K21) quadrature, with QUADPACK-style mixed convergence * criterion. The integration succeeds (returns the Boost estimate Q) * whenever * error <= max(relative_tolerance * L1, absolute_tolerance) @@ -41,9 +41,18 @@ constexpr int INTEGRATE_1D_GAUSS_KRONROD_MAX_DEPTH = 15; * checking accumulated floating-point round-off against itself (this * happens routinely in nested integrate_1d_gauss_kronrod calls when the * outer integration probes the deep tail of the integrand and every - * inner evaluation sees an essentially-zero integrand). Setting it to - * zero (the default) reproduces the strict pure-relative-tolerance - * behaviour of integrate_1d. + * inner evaluation sees an essentially-zero integrand). + * + * absolute_tolerance is applied twice, in the same units: as a floor on + * refinement inside Boost's adaptive recursion (a panel whose error already + * sits below the floor is not bisected, which is what bounds the work in the + * round-off regime above) and as the floor on the convergence test below. + * + * Setting it to zero (the default) reproduces the strict + * pure-relative-tolerance behaviour of integrate_1d. Note that zero is + * Boost's sentinel for "derive the refinement budget from the root panel's + * own relative target", so a positive but negligible absolute_tolerance + * removes that derived budget and can refine slightly MORE than zero does. * * The signature for f should be: * double f(double x, double xc) @@ -86,7 +95,7 @@ inline double integrate_gk(const F& f, double a, double b, const unsigned int depth = max_depth < 0 ? 0u : static_cast(max_depth); double Q = gauss_kronrod::integrate( - f_wrap, a, b, depth, relative_tolerance, &error, &L1); + f_wrap, a, b, depth, relative_tolerance, &error, &L1, absolute_tolerance); // QUADPACK-style mixed convergence: throw only if the Boost error // exceeds both the relative-tolerance target (rel_tol * L1) and the @@ -106,7 +115,7 @@ inline double integrate_gk(const F& f, double a, double b, /** * Compute the integral of the single variable function f from a to b to within - * a specified relative tolerance using adaptive Gauss-Kronrod (G21,K21) + * a specified relative tolerance using adaptive Gauss-Kronrod (G10,K21) * quadrature. a and b can be finite or infinite. * * @tparam F type of function to integrate @@ -147,7 +156,7 @@ inline double integrate_1d_gauss_kronrod_tol(const F& f, double a, double b, /** * Compute the integral of the single variable function f from a to b using - * adaptive Gauss-Kronrod (G21,K21) quadrature. a and b can be finite or + * adaptive Gauss-Kronrod (G10,K21) quadrature. a and b can be finite or * infinite. * * The signature for f should be: diff --git a/stan/math/rev/functor/integrate_1d_gauss_kronrod.hpp b/stan/math/rev/functor/integrate_1d_gauss_kronrod.hpp index 135c4437d40..d5b75a323e6 100644 --- a/stan/math/rev/functor/integrate_1d_gauss_kronrod.hpp +++ b/stan/math/rev/functor/integrate_1d_gauss_kronrod.hpp @@ -16,7 +16,7 @@ namespace stan { namespace math { /** - * Return the integral of f from a to b using adaptive Gauss-Kronrod (G21,K21) + * Return the integral of f from a to b using adaptive Gauss-Kronrod (G10,K21) * quadrature. * * @tparam F Type of f @@ -57,7 +57,7 @@ inline return_type_t integrate_1d_gauss_kronrod_tol( /** * Compute the integral of the single variable function f from a to b using - * adaptive Gauss-Kronrod (G21,K21) quadrature. a and b can be finite or + * adaptive Gauss-Kronrod (G10,K21) quadrature. a and b can be finite or * infinite. * * f should be compatible with reverse mode autodiff and have the signature: diff --git a/test/unit/math/mix/functor/integrate_1d_gauss_kronrod_test.cpp b/test/unit/math/mix/functor/integrate_1d_gauss_kronrod_test.cpp index fa162ecd962..11b33ca7ce1 100644 --- a/test/unit/math/mix/functor/integrate_1d_gauss_kronrod_test.cpp +++ b/test/unit/math/mix/functor/integrate_1d_gauss_kronrod_test.cpp @@ -25,3 +25,20 @@ TEST(mixFunctor, integrate1DGaussKronrod) { // a correctness issue in the integrate_1d_gauss_kronrod wrapper, and so // the NaN-input case is intentionally omitted here. } + +// A positive absolute_tolerance must not disturb higher-order autodiff. The +// integrand is scaled to 1e-6 so the refinement floor is actually reached, +// which is the regime where abs_tol changes how many panels are evaluated. +TEST(mixFunctor, integrate1DGaussKronrodPositiveAbsoluteTolerance) { + auto f = [](const auto& theta, const auto& lb, const auto& ub) { + auto func = [](const auto& x, const auto& xc, std::ostream* msgs, + const auto& theta) { + return 1e-6 * theta * stan::math::sin(7.0 * x); + }; + std::ostringstream* msgs = nullptr; + return stan::math::integrate_1d_gauss_kronrod_tol(func, lb, ub, 1e-8, 1e-10, + 5, msgs, theta); + }; + + stan::test::expect_ad(f, 0.75, 0.0, 1.0); +} diff --git a/test/unit/math/prim/functor/integrate_1d_gauss_kronrod_test.cpp b/test/unit/math/prim/functor/integrate_1d_gauss_kronrod_test.cpp index e49ae98e4dd..f46f75e732d 100644 --- a/test/unit/math/prim/functor/integrate_1d_gauss_kronrod_test.cpp +++ b/test/unit/math/prim/functor/integrate_1d_gauss_kronrod_test.cpp @@ -444,3 +444,161 @@ TEST(StanMath_integrate_1d_gk_prim, abs_tol_argument_smoke) { EXPECT_NEAR(Q0, 1.0423499493102901, 1e-8); EXPECT_NEAR(Q1, Q0, 1e-12); } + +// --------------------------------------------------------------------------- +// absolute_tolerance during refinement +// +// absolute_tolerance is applied in two places, in the same units: as a floor +// on refinement inside Boost's adaptive recursion, and as the floor on the +// convergence test. The first of those depends on the abs_tol parameter added +// to Boost's gauss_kronrod::integrate by a Stan-local patch +// (lib/boost_1.87.0/STAN_CHANGES). +// +// --------------------------------------------------------------------------- + +// The motivating case: an integrand so small that the relative-tolerance test +// degenerates into comparing accumulated round-off against itself. A positive +// abs_tol stops the pointless refinement. +// +// ALSO the regression guard on the vendored-Boost patch -- see above. +TEST(StanMath_integrate_1d_gk_prim, + positive_abs_tol_reduces_work_on_negligible_integrand) { + constexpr double scale = 1e-12; + constexpr double frequency = 127.0; + constexpr double absolute_tolerance = 1e-14; + const double expected = scale * (1.0 - std::cos(frequency)) / frequency; + + auto run = [](double abs_tol, int &evaluations) { + auto integrand = [&evaluations](double x, double xc, std::ostream *msgs) { + ++evaluations; + return scale * std::sin(frequency * x); + }; + return stan::math::integrate_1d_gauss_kronrod_tol( + integrand, 0.0, 1.0, 1e-12, abs_tol, 5, integrate_1d_gk_test::msgs); + }; + + int relative_evaluations = 0, absolute_evaluations = 0; + const double relative_result = run(0.0, relative_evaluations); + const double absolute_result = run(absolute_tolerance, absolute_evaluations); + + // The contract: within the absolute tolerance requested. Asserting anything + // tighter would assert an accident of how much more accurate K21 happens to + // be than the caller asked for. + EXPECT_NEAR(absolute_result, expected, absolute_tolerance); + EXPECT_NEAR(relative_result, expected, absolute_tolerance); + EXPECT_NEAR(absolute_result, relative_result, absolute_tolerance); + EXPECT_LT(std::abs(absolute_result - expected) / std::abs(expected), 1e-10); + + // Materially less work. Stated as a ratio so it does not pin the exact + // recursion counts of the quadrature. + EXPECT_LT(absolute_evaluations, relative_evaluations); + EXPECT_LT(2 * absolute_evaluations, relative_evaluations); +} + +// Regression guard for unbounded refinement. x^{-0.9} has an endpoint +// singularity Gauss-Kronrod cannot resolve, so no attainable tolerance is met +// and every panel looks "not yet good enough". Refinement stays panel-local, +// so the cost stays at the level of the abs_tol == 0 call. +TEST(StanMath_integrate_1d_gk_prim, + positive_abs_tol_bounds_work_on_unresolvable_integrand) { + auto run = [](double absolute_tolerance, int &evaluations) { + auto integrand = [&evaluations](double x, double xc, std::ostream *msgs) { + ++evaluations; + return std::pow(x, -0.9); + }; + EXPECT_THROW(stan::math::integrate_1d_gauss_kronrod_tol( + integrand, 1e-300, 1.0, 1e-12, absolute_tolerance, 15, + integrate_1d_gk_test::msgs), + std::domain_error); + }; + + int relative_evaluations = 0, absolute_evaluations = 0; + run(0.0, relative_evaluations); + run(1e-14, absolute_evaluations); + + EXPECT_LE(absolute_evaluations, relative_evaluations); + EXPECT_LT(absolute_evaluations, 50000); +} + +// Among POSITIVE absolute tolerances, raising the floor never increases work. +// Zero is excluded because it is the derived-budget sentinel above. +TEST(StanMath_integrate_1d_gk_prim, work_is_monotone_in_positive_abs_tol) { + auto evaluations_for = [](double absolute_tolerance) { + int evaluations = 0; + auto integrand = [&evaluations](double x, double xc, std::ostream *msgs) { + ++evaluations; + return std::exp(-x * x) * std::cos(30 * x); + }; + try { + stan::math::integrate_1d_gauss_kronrod_tol(integrand, 0.0, 3.0, 1e-12, + absolute_tolerance, 12, + integrate_1d_gk_test::msgs); + } catch (const std::domain_error &) { + // Convergence failure is irrelevant here; only the cost is. + } + return evaluations; + }; + + int previous = evaluations_for(1e-300); + for (double absolute_tolerance : {1e-16, 1e-12, 1e-8, 1e-4, 1e-1}) { + const int current = evaluations_for(absolute_tolerance); + EXPECT_LE(current, previous) << "abs_tol = " << absolute_tolerance; + previous = current; + } +} + +// A positive abs_tol must not paper over a genuinely unconverged result. +TEST(StanMath_integrate_1d_gk_prim, positive_abs_tol_still_throws_when_needed) { + auto integrand = [](double x, double xc, std::ostream *msgs) { + return std::pow(x, -0.9); + }; + EXPECT_THROW( + stan::math::integrate_1d_gauss_kronrod_tol( + integrand, 1e-300, 1.0, 1e-12, 1e-14, 15, integrate_1d_gk_test::msgs), + std::domain_error); +} + +// max_depth = 0 means "one panel, no bisection"; a positive abs_tol must not +// disturb that (the refinement floor is never consulted). +TEST(StanMath_integrate_1d_gk_prim, positive_abs_tol_with_max_depth_zero) { + auto integrand + = [](double x, double xc, std::ostream *msgs) { return std::exp(x); }; + const double Q = stan::math::integrate_1d_gauss_kronrod_tol( + integrand, 0.0, 1.0, 1e-10, 1e-12, 0, integrate_1d_gk_test::msgs); + EXPECT_NEAR(Q, std::exp(1.0) - 1.0, 1e-12); +} + +// A positive abs_tol must agree with the abs_tol == 0 answer under every +// change of variable, not just the finite one. +TEST(StanMath_integrate_1d_gk_prim, positive_abs_tol_domain_transformations) { + constexpr double relative_tolerance = 1e-12; + constexpr double absolute_tolerance = 1e-10; + constexpr int max_depth = 8; + const double infinity = std::numeric_limits::infinity(); + + auto check_integral = [&](const auto &integrand, double lower, double upper, + double expected) { + const double legacy_result = stan::math::integrate_1d_gauss_kronrod_tol( + integrand, lower, upper, relative_tolerance, 0.0, max_depth, + integrate_1d_gk_test::msgs); + const double absolute_result = stan::math::integrate_1d_gauss_kronrod_tol( + integrand, lower, upper, relative_tolerance, absolute_tolerance, + max_depth, integrate_1d_gk_test::msgs); + EXPECT_NEAR(absolute_result, expected, absolute_tolerance); + EXPECT_NEAR(absolute_result, legacy_result, absolute_tolerance); + }; + + auto increasing_exponential + = [](double x, double xc, std::ostream *msgs) { return std::exp(x); }; + auto decreasing_exponential + = [](double x, double xc, std::ostream *msgs) { return std::exp(-x); }; + auto gaussian_kernel = [](double x, double xc, std::ostream *msgs) { + return std::exp(-x * x); + }; + + check_integral(increasing_exponential, 0.0, 1.0, std::exp(1.0) - 1.0); + check_integral(decreasing_exponential, 0.0, infinity, 1.0); + check_integral(increasing_exponential, -infinity, 0.0, 1.0); + check_integral(gaussian_kernel, -infinity, infinity, + std::sqrt(stan::math::pi())); +} diff --git a/test/unit/math/rev/functor/integrate_1d_gauss_kronrod_test.cpp b/test/unit/math/rev/functor/integrate_1d_gauss_kronrod_test.cpp index 7dfc53065d6..dfba6650b2b 100644 --- a/test/unit/math/rev/functor/integrate_1d_gauss_kronrod_test.cpp +++ b/test/unit/math/rev/functor/integrate_1d_gauss_kronrod_test.cpp @@ -481,4 +481,52 @@ TEST_F(AgradRev, StanMath_integrate_1d_gk_rev_TestUniform) { EXPECT_FLOAT_EQ(1, 1 + g[1]); } +// Reverse-mode counterpart of +// StanMath_integrate_1d_gk_prim.positive_abs_tol_reduces_work_on_negligible_ +// integrand. The gradient is computed by integrating the adjoint, so the +// refinement floor has to hold for that pass too -- and in practice that is +// where it matters most, because a parameter derivative is routinely orders +// of magnitude smaller than the value it differentiates. +TEST_F(AgradRev, StanMath_integrate_1d_gk_rev_positive_abs_tol_reduces_work) { + using stan::math::var; + + constexpr double scale = 1e-12; + constexpr double frequency = 127.0; + constexpr double absolute_tolerance = 1e-14; + const double expected_adjoint + = scale * (1.0 - std::cos(frequency)) / frequency; + const double theta_value = 1.7; + + auto run = [&](double abs_tol, int &evaluations, double &value, + double &adjoint) { + stan::math::nested_rev_autodiff nested; + auto integrand = [&evaluations](double x, double xc, std::ostream *msgs, + const auto &theta) { + ++evaluations; + return theta[0] * scale * stan::math::sin(frequency * x); + }; + std::vector theta{theta_value}; + var integral = stan::math::integrate_1d_gauss_kronrod_tol( + integrand, 0.0, 1.0, 1e-12, abs_tol, 5, msgs, theta); + integral.grad(); + value = integral.val(); + adjoint = theta[0].adj(); + }; + + int relative_evaluations = 0, absolute_evaluations = 0; + double relative_value = 0, relative_adjoint = 0; + double absolute_value = 0, absolute_adjoint = 0; + run(0.0, relative_evaluations, relative_value, relative_adjoint); + run(absolute_tolerance, absolute_evaluations, absolute_value, + absolute_adjoint); + + EXPECT_NEAR(absolute_value, theta_value * expected_adjoint, + absolute_tolerance); + EXPECT_NEAR(absolute_adjoint, expected_adjoint, absolute_tolerance); + EXPECT_NEAR(absolute_value, relative_value, absolute_tolerance); + EXPECT_NEAR(absolute_adjoint, relative_adjoint, absolute_tolerance); + + EXPECT_LT(absolute_evaluations, relative_evaluations); +} + } // namespace integrate_1d_gk_test