From 7fb286f11511bdee8fda496af326b96a626196b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Fri, 7 Aug 2026 13:48:13 +0100 Subject: [PATCH] Fix flaky Weber_BMC2015 gradient check `test_benchmark_gradient[Weber_BMC2015-*-unscaled]` has been flip-flopping between passing and failing on `Benchmark Collection CPP (3.14, false)` -- two runs of the same tree (differing only in CHANGELOG.md, with identical dependency versions) produced one pass and one failure. The forward and adjoint AMICI gradients agree to nine significant digits, so the gradients are fine; the finite differences are not. The failing direction is `a32`, where fiddy reports 860.83 against AMICI's 872.68 (1.36%, just over `rtol_check=1e-2`). Root cause: fiddy's `Consistency` checker averages the estimates from every step size whose forward/backward/central differences agree with each other, and Weber has no step size window that suits all directions at once. The noise parameters (e.g. `std_yPKDpN25`) have such a large second derivative that forward and backward differences only come within 5% of the central one for relative steps <= 5e-5. At those steps the differences for parameters with a small nominal value (`a32` ~ 1e-4) are integrator-noise dominated -- the central difference is off by 0.8% at 1e-4 and 4.2% at 1e-5, while being off by only 2e-6 at 2e-1. Those noisy estimates are internally consistent (the bias comes from a shifted solver trajectory, not from per-evaluation jitter), so `Consistency` accepts and averages them in. Whether a given step size clears the internal agreement check depends on last-ulp details of the runner, which is what makes the outcome vary between runs. Use central differences only for this problem, over the step size window in which every direction is converged. Measured worst-case relative deviation from the AMICI gradient over {forward, adjoint} x {scaled, unscaled} x {AMICI_EXTRACT_CSE=true, false} drops from 5.5e-3 to 1.9e-3, against an unchanged `rtol_check` of 1e-2. `method_ids` becomes a `GradientCheckSettings` field, defaulting to the previous forward/backward/central triple, so all other problems are unaffected. Co-Authored-By: Claude Opus 5 --- .../benchmark_models/test_petab_benchmark.py | 25 +++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/tests/benchmark_models/test_petab_benchmark.py b/tests/benchmark_models/test_petab_benchmark.py index fb58aa8a4b..52b044baef 100644 --- a/tests/benchmark_models/test_petab_benchmark.py +++ b/tests/benchmark_models/test_petab_benchmark.py @@ -173,6 +173,15 @@ class GradientCheckSettings: 1e-5, ] ) + # Finite difference schemes to use. All of them have to agree within + # (atol|rtol)_consistency for a step size to be considered. + method_ids: list[MethodId] = field( + default_factory=lambda: [ + MethodId.CENTRAL, + MethodId.FORWARD, + MethodId.BACKWARD, + ] + ) rng_seed: int = 0 ss_computation_mode: SteadyStateComputationMode = ( SteadyStateComputationMode.integrationOnly @@ -237,6 +246,18 @@ class GradientCheckSettings: atol_check=1e-6, rtol_check=1e-2, rng_seed=2, + # No single step size works for all directions with the default + # forward/backward/central agreement requirement: the noise parameters + # (e.g. std_yPKDpN25) are so strongly curved that forward and backward + # differences only agree with the central one for steps <= 5e-5, while at + # those steps the finite differences of parameters with a small nominal + # value (e.g. a32 ~ 1e-4) are dominated by integrator noise. Since + # `Consistency` averages over all step sizes it accepts, the noisy small + # steps then bias the result by ~1%, right at `rtol_check` -- which is what + # made this test flip-flop between runs. Use central differences only, over + # the step size window in which every direction is converged. + method_ids=[MethodId.CENTRAL], + step_sizes=[5e-3, 2e-3, 1e-3, 5e-4], ) settings["Zheng_PNAS2012"] = GradientCheckSettings( rng_seed=1, @@ -484,7 +505,7 @@ def test_benchmark_gradient( point=point, sizes=cur_settings.step_sizes, direction_ids=parameter_ids, - method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], + method_ids=cur_settings.method_ids, success_checker=Consistency( rtol=cur_settings.rtol_consistency, atol=cur_settings.atol_consistency, @@ -810,7 +831,7 @@ def test_nominal_parameters_llh_v2(problem_id): point=point, sizes=cur_settings.step_sizes, direction_ids=parameter_ids, - method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], + method_ids=cur_settings.method_ids, success_checker=Consistency( rtol=cur_settings.rtol_consistency, atol=cur_settings.atol_consistency,