From 3bef8848b0b1560b646f0a4e33c023a76339a489 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Tue, 15 Sep 2026 21:35:35 +0300 Subject: [PATCH 1/5] Add abs_tol parameter to vendored Boost gauss_kronrod::integrate --- .gitattributes | 4 ++++ lib/boost_1.87.0/STAN_CHANGES | 21 +++++++++++++++++++ .../boost/math/quadrature/gauss_kronrod.hpp | 12 +++++------ 3 files changed, 31 insertions(+), 6 deletions(-) create mode 100644 lib/boost_1.87.0/STAN_CHANGES 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..603297decb4 --- /dev/null +++ b/lib/boost_1.87.0/STAN_CHANGES @@ -0,0 +1,21 @@ +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)`. + + The recursion already accepts an absolute-error budget and honours it in its + leaf criterion; only the public entry point was missing a way to set it. + Without this, a panel whose integrand is orders of magnitude below the value + being integrated -- e.g. the derivative of a marginal likelihood with respect + to a concentration parameter in a saturated tail -- can never satisfy the + leaf-relative test and bisects to max_depth, costing ~2^15 * 21 evaluations + per panel for mass that is numerically negligible. + + The parameter is trailing and defaults to `Real(0)`, which is the value + previously hard-coded, so every existing caller is unaffected and the + quadrature is bit-for-bit unchanged for them. + + 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())); From a9c5f69959b6d7f9baacac26bb8b057c021530ad Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Tue, 15 Sep 2026 21:35:47 +0300 Subject: [PATCH 2/5] Fix Gauss-Kronrod order in docs: the pair is (G10,K21), not (G21,K21) --- stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp | 4 ++-- stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp | 6 +++--- stan/math/rev/functor/integrate_1d_gauss_kronrod.hpp | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) 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..462a472bd2e 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) @@ -106,7 +106,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 +147,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: From 27fd485c4477c34c8d72212afdf6daac93aa66db Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Tue, 15 Sep 2026 21:36:28 +0300 Subject: [PATCH 3/5] Apply absolute tolerance during Gauss-Kronrod refinement --- .../functor/integrate_1d_gauss_kronrod.hpp | 20 +- .../integrate_1d_gauss_kronrod_test.cpp | 17 ++ .../integrate_1d_gauss_kronrod_test.cpp | 188 ++++++++++++++++++ .../integrate_1d_gauss_kronrod_test.cpp | 48 +++++ 4 files changed, 269 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp b/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp index 462a472bd2e..7d6c4e1ffcc 100644 --- a/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp +++ b/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp @@ -41,9 +41,21 @@ 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. + * Reaching the first of those needs the abs_tol parameter added to Boost's + * gauss_kronrod::integrate by a Stan-local patch; see + * lib/boost_1.87.0/STAN_CHANGES. + * + * 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 +98,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 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..866191c3db1 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,191 @@ 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). +// +// positive_abs_tol_reduces_work_on_negligible_integrand below is the guard on +// that patch: if a Boost upgrade drops it, absolute_tolerance stops affecting +// refinement, the two evaluation counts become equal, and the test fails. +// --------------------------------------------------------------------------- + +// 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); +} + +// absolute_tolerance == 0 is Boost's sentinel for "derive the refinement +// budget from the root panel's own relative target", NOT "no absolute floor". +// A positive but negligible value therefore removes that derived budget and +// refines at least as much as zero does. This is a wart of Boost's interface, +// not of the patch; it is pinned here so the behaviour is documented rather +// than discovered. +TEST(StanMath_integrate_1d_gk_prim, abs_tol_zero_is_a_derived_budget_sentinel) { + auto run = [](double absolute_tolerance, int *evaluations) { + auto integrand = [evaluations](double x, double xc, std::ostream *msgs) { + ++*evaluations; + return std::exp(-x * x) * std::cos(30 * x); + }; + return stan::math::integrate_1d_gauss_kronrod_tol( + integrand, 0.0, 3.0, 1e-10, absolute_tolerance, 15, + integrate_1d_gk_test::msgs); + }; + + int zero_evaluations = 0, tiny_evaluations = 0; + const double zero_result = run(0.0, &zero_evaluations); + const double tiny_result = run(1e-300, &tiny_evaluations); + + // Both are valid answers to the requested relative tolerance ... + EXPECT_NEAR(zero_result, tiny_result, 1e-10 * std::abs(zero_result)); + // ... but the negligible floor costs at least as much work as no floor. + EXPECT_GE(tiny_evaluations, zero_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..9aafe8a931b 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 From 58c35e4785f16cca614bde9e2132ff7bccc766e4 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Tue, 15 Sep 2026 23:22:01 +0300 Subject: [PATCH 4/5] Trim abs_tol documentation and drop the sentinel test --- lib/boost_1.87.0/STAN_CHANGES | 12 -------- .../functor/integrate_1d_gauss_kronrod.hpp | 3 -- .../integrate_1d_gauss_kronrod_test.cpp | 30 ------------------- 3 files changed, 45 deletions(-) diff --git a/lib/boost_1.87.0/STAN_CHANGES b/lib/boost_1.87.0/STAN_CHANGES index 603297decb4..dc3f902b487 100644 --- a/lib/boost_1.87.0/STAN_CHANGES +++ b/lib/boost_1.87.0/STAN_CHANGES @@ -6,16 +6,4 @@ This file documents changes done for the stan-math project `recursive_adaptive_integrate()` entry points in place of the hard-coded `Real(0)`. - The recursion already accepts an absolute-error budget and honours it in its - leaf criterion; only the public entry point was missing a way to set it. - Without this, a panel whose integrand is orders of magnitude below the value - being integrated -- e.g. the derivative of a marginal likelihood with respect - to a concentration parameter in a saturated tail -- can never satisfy the - leaf-relative test and bisects to max_depth, costing ~2^15 * 21 evaluations - per panel for mass that is numerically negligible. - - The parameter is trailing and defaults to `Real(0)`, which is the value - previously hard-coded, so every existing caller is unaffected and the - quadrature is bit-for-bit unchanged for them. - Used by `stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp`. diff --git a/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp b/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp index 7d6c4e1ffcc..3a3ef788bbd 100644 --- a/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp +++ b/stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp @@ -47,9 +47,6 @@ constexpr int INTEGRATE_1D_GAUSS_KRONROD_MAX_DEPTH = 15; * 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. - * Reaching the first of those needs the abs_tol parameter added to Boost's - * gauss_kronrod::integrate by a Stan-local patch; see - * lib/boost_1.87.0/STAN_CHANGES. * * Setting it to zero (the default) reproduces the strict * pure-relative-tolerance behaviour of integrate_1d. Note that zero is 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 866191c3db1..124fc3e2851 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 @@ -454,9 +454,6 @@ TEST(StanMath_integrate_1d_gk_prim, abs_tol_argument_smoke) { // to Boost's gauss_kronrod::integrate by a Stan-local patch // (lib/boost_1.87.0/STAN_CHANGES). // -// positive_abs_tol_reduces_work_on_negligible_integrand below is the guard on -// that patch: if a Boost upgrade drops it, absolute_tolerance stops affecting -// refinement, the two evaluation counts become equal, and the test fails. // --------------------------------------------------------------------------- // The motivating case: an integrand so small that the relative-tolerance test @@ -498,33 +495,6 @@ TEST(StanMath_integrate_1d_gk_prim, EXPECT_LT(2 * absolute_evaluations, relative_evaluations); } -// absolute_tolerance == 0 is Boost's sentinel for "derive the refinement -// budget from the root panel's own relative target", NOT "no absolute floor". -// A positive but negligible value therefore removes that derived budget and -// refines at least as much as zero does. This is a wart of Boost's interface, -// not of the patch; it is pinned here so the behaviour is documented rather -// than discovered. -TEST(StanMath_integrate_1d_gk_prim, abs_tol_zero_is_a_derived_budget_sentinel) { - auto run = [](double absolute_tolerance, int *evaluations) { - auto integrand = [evaluations](double x, double xc, std::ostream *msgs) { - ++*evaluations; - return std::exp(-x * x) * std::cos(30 * x); - }; - return stan::math::integrate_1d_gauss_kronrod_tol( - integrand, 0.0, 3.0, 1e-10, absolute_tolerance, 15, - integrate_1d_gk_test::msgs); - }; - - int zero_evaluations = 0, tiny_evaluations = 0; - const double zero_result = run(0.0, &zero_evaluations); - const double tiny_result = run(1e-300, &tiny_evaluations); - - // Both are valid answers to the requested relative tolerance ... - EXPECT_NEAR(zero_result, tiny_result, 1e-10 * std::abs(zero_result)); - // ... but the negligible floor costs at least as much work as no floor. - EXPECT_GE(tiny_evaluations, zero_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, From 36fd07225933c7088b41565b26ae0e4fd4cd2117 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Tue, 15 Sep 2026 23:22:07 +0300 Subject: [PATCH 5/5] Pass evaluation counters by reference in Gauss-Kronrod tests --- .../integrate_1d_gauss_kronrod_test.cpp | 20 +++++++++---------- .../integrate_1d_gauss_kronrod_test.cpp | 20 +++++++++---------- 2 files changed, 20 insertions(+), 20 deletions(-) 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 124fc3e2851..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 @@ -468,9 +468,9 @@ TEST(StanMath_integrate_1d_gk_prim, 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; + 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( @@ -478,8 +478,8 @@ TEST(StanMath_integrate_1d_gk_prim, }; 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); + 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 @@ -501,9 +501,9 @@ TEST(StanMath_integrate_1d_gk_prim, // 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; + 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( @@ -513,8 +513,8 @@ TEST(StanMath_integrate_1d_gk_prim, }; int relative_evaluations = 0, absolute_evaluations = 0; - run(0.0, &relative_evaluations); - run(1e-14, &absolute_evaluations); + run(0.0, relative_evaluations); + run(1e-14, absolute_evaluations); EXPECT_LE(absolute_evaluations, relative_evaluations); EXPECT_LT(absolute_evaluations, 50000); 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 9aafe8a931b..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 @@ -497,28 +497,28 @@ TEST_F(AgradRev, StanMath_integrate_1d_gk_rev_positive_abs_tol_reduces_work) { = scale * (1.0 - std::cos(frequency)) / frequency; const double theta_value = 1.7; - auto run = [&](double abs_tol, int *evaluations, double *value, - double *adjoint) { + 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; + 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(); + 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); + 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);