diff --git a/python/sdist/amici/adapters/fiddy.py b/python/sdist/amici/adapters/fiddy.py index d91841ff8d..e07ac49106 100644 --- a/python/sdist/amici/adapters/fiddy.py +++ b/python/sdist/amici/adapters/fiddy.py @@ -10,6 +10,7 @@ from __future__ import annotations +import warnings from collections.abc import Callable from functools import partial from inspect import signature @@ -18,6 +19,8 @@ import numpy as np import petab.v1 as petab from fiddy import CachedFunction, Type, fiddy_array +from fiddy.directional_derivative import DirectionalDerivative +from fiddy.success import Consistency from petab.v1.C import LIN, LOG, LOG10 from amici.sim.sundials import ( @@ -38,6 +41,7 @@ from amici.sim.sundials.petab import PetabSimulationResult, PetabSimulator __all__ = [ + "RobustConsistency", "run_simulation_to_cached_functions", "simulate_petab_to_cached_functions", "simulate_petab_v2_to_cached_functions", @@ -46,6 +50,219 @@ LOG_E_10 = np.log(10) +class RobustConsistency(Consistency): + """`Consistency`, plus rejection of step sizes that are self-consistent + but inconsistent with the majority of other step sizes. + + `Consistency` checks whether the requested methods (e.g. + forward/backward/central) agree with each other at each step size + ("self-consistent"), then blends every self-consistent size's mean into + the final value. Self-consistency alone is not a strong guarantee on its + own: a step size can be small enough that all methods sample points + within the target function's floating-point noise floor and become + correlated (affected by the same rounding/cancellation error) -- + self-consistent, yet biased away from the truth. Symmetrically, a step + size can also be large enough that all methods are biased the same way + by higher-order/truncation effects. + + To guard against this, self-consistent step sizes are additionally + required to agree with the majority of other self-consistent step sizes, + via iterative outlier rejection (order-independent; step size magnitude + is not used as a proxy for trustworthiness): repeatedly compute the + median and a robust (MAD-based) spread of the current candidates, and + drop the single worst-deviating one if it exceeds ``trend_n_sigma`` + scaled MADs from the median, until nothing looks anomalous. This only + activates once there are at least ``min_trend_samples`` self-consistent + step sizes; below that, there isn't enough data to estimate a spread, and + all self-consistent step sizes are used, as in `Consistency`. A + `UserWarning` is emitted whenever one or more step sizes are rejected + this way. + + This addresses a long-standing intermittent CI failure in AMICI's PEtab + benchmark gradient test + (``test_benchmark_gradient[Weber_BMC2015-*-unscaled]``, see + https://github.com/AMICI-dev/AMICI/issues/3078): that test uses + `Consistency` to finite-difference-check an analytically computed + gradient for a model parameter (``a32``) several orders of magnitude + smaller than the model's other free parameters, and a small step size + could become spuriously self-consistent while biased away from the true + derivative. + + Note that this is a majority-vote style method: like any check based + purely on the agreement of the values themselves (no independent ground + truth), it has a breakdown point of roughly 50% (a property of the + underlying median/MAD statistics) -- if close to half (or more) of the + self-consistent step sizes are corrupted, this check cannot reliably + tell which subset is trustworthy. This is a fundamental limitation of + any purely data-driven consistency check, not something this + implementation can detect or work around; sufficient step sizes with a + real chance of being individually trustworthy should be provided. + + This was originally proposed upstream, in fiddy, as + https://github.com/ICB-DCM/fiddy/pull/77, but was not merged; it lives + here instead. + """ + + id = "robust_consistency" + + def __init__( + self, + *args, + trend_n_sigma: float = 5.0, + min_trend_samples: int = 3, + **kwargs, + ): + """Construct. + + :param trend_n_sigma: + The number of scaled median-absolute-deviations a + self-consistent step size's estimate may deviate from the + median of the other trusted step sizes' estimates, before it + is rejected as an outlier. + :param min_trend_samples: + The minimum number of self-consistent step sizes required + before the cross-step-size outlier rejection is attempted. + Below this, all self-consistent step sizes are trusted, same + as in `Consistency`. + :param args: + Positional arguments passed to `Consistency.__init__`. + :param kwargs: + Keyword arguments passed to `Consistency.__init__` + (e.g. ``rtol``, ``atol``, ``equal_nan``). + """ + super().__init__(*args, **kwargs) + self.trend_n_sigma = trend_n_sigma + self.min_trend_samples = min_trend_samples + + def _self_consistent_means( + self, directional_derivative: DirectionalDerivative + ) -> list[Type.DIRECTIONAL_DERIVATIVE]: + """Group results by step size, and return the per-size mean for + every step size whose requested methods agree with each other + ("self-consistent") within ``rtol/2``, ``atol/2``.""" + computer_results = directional_derivative.get_computer_results() + analysis_results = directional_derivative.get_analysis_results() + results_by_size = {} + for result in [*computer_results, *analysis_results]: + size = result.metadata.get("size_absolute", None) + if size is None: + continue + if size not in results_by_size: + results_by_size[size] = {} + if result.method_id in results_by_size[size]: + raise ValueError( + f"Duplicate, and possibly conflicting, results for method " + f'"{result.method_id}" and size "{size}".', + ) + results_by_size[size][result.method_id] = result.value + + self_consistent_means = [] + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", "Mean of empty slice", RuntimeWarning + ) + for results in results_by_size.values(): + values = list(results.values()) + mean = np.nanmean(values, axis=0) + is_self_consistent = np.isclose( + values, + mean, + rtol=self.rtol / 2, + atol=self.atol / 2, + equal_nan=self.equal_nan, + ).all() + if is_self_consistent: + self_consistent_means.append(mean) + return self_consistent_means + + def method( + self, directional_derivative: DirectionalDerivative + ) -> tuple[bool, float]: + self_consistent_means = self._self_consistent_means( + directional_derivative + ) + + if not self_consistent_means: + return False, np.nan + + trusted_means = self._reject_outliers(self_consistent_means) + + if not trusted_means: + return False, np.nan + + n_rejected = len(self_consistent_means) - len(trusted_means) + if n_rejected: + warnings.warn( + f"{n_rejected} step size(s) were self-consistent (the " + "requested methods agreed with each other) but were " + "rejected as inconsistent with the majority of other step " + "sizes; see `RobustConsistency`'s docstring.", + stacklevel=2, + ) + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", "Mean of empty slice", RuntimeWarning + ) + value = np.nanmean(trusted_means, axis=0) + + success = ( + np.isclose( + trusted_means, + value, + rtol=self.rtol, + atol=self.atol, + equal_nan=self.equal_nan, + ).all() + and not np.isnan(trusted_means).all() + ) + return success, value + + def _reject_outliers( + self, means: list[Type.DIRECTIONAL_DERIVATIVE] + ) -> list[Type.DIRECTIONAL_DERIVATIVE]: + """Iteratively reject step sizes whose estimate is an outlier. + + See the class docstring for the rationale. Order-independent: does + not assume larger (or smaller) step sizes are inherently more + trustworthy. + + :param means: + The per-step-size mean estimates that passed the + within-step-size self-consistency check. + :return: + The subset of `means` that are also mutually consistent with + each other. + """ + trusted = list(means) + if len(trusted) < self.min_trend_samples: + return trusted + + floor = max(self.atol / 2, np.finfo(float).tiny) + while len(trusted) >= self.min_trend_samples: + stacked = np.asarray(trusted, dtype=float) + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", "All-NaN", RuntimeWarning) + center = np.nanmedian(stacked, axis=0) + mad = np.nanmedian(np.abs(stacked - center), axis=0) + scale = np.maximum(mad * 1.4826, floor) + # One badness score per candidate, reduced across all + # output dimensions (a candidate is an outlier if it + # deviates too much in *any* output element). + badness = np.nanmax( + (np.abs(stacked - center) / scale).reshape( + len(trusted), -1 + ), + axis=1, + ) + worst = int(np.nanargmax(badness)) + if badness[worst] > self.trend_n_sigma: + trusted.pop(worst) + else: + break + return trusted + + def _transform_gradient_lin_to_lin(gradient_value, _): return gradient_value diff --git a/python/tests/adapters/test_fiddy.py b/python/tests/adapters/test_fiddy.py index fa5750d4c0..b60c042799 100644 --- a/python/tests/adapters/test_fiddy.py +++ b/python/tests/adapters/test_fiddy.py @@ -5,6 +5,7 @@ import numpy as np import pytest from amici.adapters.fiddy import ( + RobustConsistency, run_simulation_to_cached_functions, simulate_petab_to_cached_functions, ) @@ -12,9 +13,10 @@ from amici.sim.sundials import SensitivityOrder, SteadyStateSensitivityMode from fiddy import MethodId, Type, get_derivative from fiddy.derivative_check import NumpyIsCloseDerivativeCheck -from fiddy.success import Consistency +from fiddy.directional_derivative import ComputerResult from numpy.testing import assert_allclose from petab import v1 +from scipy.optimize import rosen # Absolute and relative tolerances for finite difference gradient checks. ATOL: float = 1e-3 @@ -77,7 +79,7 @@ def test_run_amici_simulation_to_functions(problem_generator): # analysis_classes=[ # lambda: TransformByDirectionScale(scales=parameter_scales), # ], - success_checker=Consistency(atol=1e-2), + success_checker=RobustConsistency(atol=1e-2), ) test_derivative = derivative.value @@ -157,7 +159,7 @@ def test_simulate_petab_to_functions(problem_generator, scaled_parameters): sizes=[1e-10, 1e-5, 1e-3, 1e-1], direction_ids=free_parameter_ids, method_ids=[MethodId.FORWARD, MethodId.BACKWARD, MethodId.CENTRAL], - success_checker=Consistency(), + success_checker=RobustConsistency(), ) check = NumpyIsCloseDerivativeCheck( @@ -167,3 +169,246 @@ def test_simulate_petab_to_functions(problem_generator, scaled_parameters): ) result = check(rtol=1e-2) assert result.success + + +class FakeDirectionalDerivative: + """Minimal stand-in exposing only what `RobustConsistency.method` calls.""" + + def __init__(self, computer_results, analysis_results=None): + self._computer_results = computer_results + self._analysis_results = analysis_results or [] + + def get_computer_results(self): + return self._computer_results + + def get_analysis_results(self): + return self._analysis_results + + +def test_robust_consistency_rejects_rounding_noise_dominated_step_sizes(): + """Regression test for the mechanism behind the flaky + `test_benchmark_gradient[Weber_BMC2015-*-unscaled]` failures + (AMICI-dev/AMICI#3078). + + A step size can become small enough that forward/backward/central all + sample points within the target function's floating-point noise floor. + They then become correlated (affected by the same rounding/cancellation + error), and can spuriously agree with each other ("self-consistent") + while being biased away from the true derivative. `RobustConsistency` + must not blend such a step size into the final value while reporting + `success=True`. + """ + true_slope = 872.68 + noise_floor = 2e-7 + + def f(point): + x0 = point[0] + value = -1023.447 + true_slope * (x0 - 1e-4) + value += noise_floor * np.sin(1e8 * x0) + return np.array(value) + + point = np.array([9.579126317171899e-05]) + step_sizes = [5e-1, 2e-1, 1e-1, 5e-2, 1e-2, 1e-3, 1e-4, 1e-5] + + with pytest.warns(UserWarning, match="rejected as inconsistent"): + derivative = get_derivative( + function=f, + point=point, + sizes=step_sizes, + direction_ids=["x0"], + method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], + success_checker=RobustConsistency(rtol=0.1, atol=1e-5), + relative_sizes=True, + ) + + success = bool(derivative.df["success"].values[0]) + value = float(np.squeeze(derivative.series.values[0])) + + # Reporting `success=True` is only acceptable if the value is actually + # accurate; silently returning a significantly biased value (as + # `Consistency` does: ~868.5, a ~0.5% error) is the bug being fixed. + if success: + assert np.isclose(value, true_slope, rtol=1e-2) + + +def test_robust_consistency_averages_all_trustworthy_step_sizes(): + """A wide, but genuinely well-behaved, range of step sizes should not + trigger spurious outlier rejection (and thus no rejection warning).""" + + def f(point): + return np.array([rosen(point)]) + + point = np.array([1.3, 0.7]) + # Chosen to have comparable precision across the whole range (see + # `test_robust_consistency_narrows_to_the_most_precise_step_sizes` below + # for what happens once the range gets wide enough that the smallest + # steps are far more precise than the largest). + step_sizes = [1e-2, 1e-3, 1e-4] + + derivative = get_derivative( + function=f, + point=point, + sizes=step_sizes, + direction_ids=["x0"], + directions=[np.array([1.0, 0.0])], + method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], + success_checker=RobustConsistency(rtol=1e-2, atol=1e-8), + ) + + assert bool(derivative.df["success"].values[0]) + value = float(np.squeeze(derivative.series.values[0])) + h = 1e-6 + expected = ( + rosen(point + np.array([h, 0.0])) - rosen(point - np.array([h, 0.0])) + ) / (2 * h) + assert np.isclose(value, expected, rtol=1e-3) + + +def test_robust_consistency_narrows_to_the_most_precise_step_sizes(): + """Rejection isn't only about *biased* step sizes (the motivating bug): + a genuinely wide, noise-free step-size range can legitimately narrow + down to just the handful of smallest, most precise steps, even though + the larger, excluded ones weren't wrong -- just comparatively imprecise + (ordinary, shrinking-with-h truncation error) next to a cluster that + happens to already be near machine precision. The blended value must + stay accurate either way. + """ + + def f(point): + return np.array([rosen(point)]) + + point = np.array([1.3, 0.7]) + step_sizes = [1e-2, 1e-3, 1e-4, 1e-5, 1e-6] + + with pytest.warns(UserWarning, match="rejected as inconsistent"): + derivative = get_derivative( + function=f, + point=point, + sizes=step_sizes, + direction_ids=["x0"], + directions=[np.array([1.0, 0.0])], + method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], + success_checker=RobustConsistency(rtol=1e-2, atol=1e-8), + ) + + assert bool(derivative.df["success"].values[0]) + value = float(np.squeeze(derivative.series.values[0])) + h = 1e-6 + expected = ( + rosen(point + np.array([h, 0.0])) - rosen(point - np.array([h, 0.0])) + ) / (2 * h) + assert np.isclose(value, expected, rtol=1e-6) + + +def test_robust_consistency_warns_when_rejecting_step_sizes(): + """`RobustConsistency` should tell the user when it rejects a step size + that looked self-consistent on its own -- otherwise a legitimate-seeming + result could silently vanish from the blend without a trace.""" + results = [ + ComputerResult( + method_id="central", value=10.0, metadata={"size_absolute": 1.0} + ), + ComputerResult( + method_id="central", value=10.01, metadata={"size_absolute": 0.5} + ), + ComputerResult( + method_id="central", value=9.99, metadata={"size_absolute": 0.1} + ), + ComputerResult( + method_id="central", + value=500.0, + metadata={"size_absolute": 0.01}, + ), + ] + + checker = RobustConsistency() + with pytest.warns(UserWarning, match="1 step size"): + success, value = checker.method(FakeDirectionalDerivative(results)) + + assert success + assert np.isclose(value, np.mean([10.0, 10.01, 9.99])) + + +class TestRejectOutliers: + """Unit tests for `RobustConsistency._reject_outliers`, the + order-independent iterative outlier-rejection pass over step sizes' + per-size means.""" + + def test_below_min_trend_samples_keeps_everything(self): + # An outlier (500.0) is present, but there are fewer candidates than + # `min_trend_samples` -- too little data to estimate a spread, so no + # trimming is attempted at all. + checker = RobustConsistency(min_trend_samples=5) + means = [10.0, 10.01, 500.0] + assert checker._reject_outliers(means) == means + + def test_no_outliers_keeps_all(self): + checker = RobustConsistency() + means = [10.0, 10.01, 9.99, 10.02] + assert checker._reject_outliers(means) == means + + def test_removes_single_outlier(self): + checker = RobustConsistency() + means = [10.0, 10.01, 9.99, 10.02, 500.0] + trusted = checker._reject_outliers(means) + assert trusted == [10.0, 10.01, 9.99, 10.02] + + def test_removes_multiple_outliers_iteratively(self): + # Two outliers on opposite sides of the trustworthy cluster; both + # must be dropped, one per iteration, worst-first. + checker = RobustConsistency() + means = [10.0, 10.01, 9.99, 10.02, 500.0, -500.0] + trusted = checker._reject_outliers(means) + assert trusted == [10.0, 10.01, 9.99, 10.02] + + def test_order_independent(self): + # Dropping is based on value, not position: shuffling the input + # must not change which candidates survive. + checker = RobustConsistency() + means = [500.0, 10.0, 10.01, 9.99, 10.02] + trusted = checker._reject_outliers(means) + assert sorted(trusted) == [9.99, 10.0, 10.01, 10.02] + + def test_respects_trend_n_sigma(self): + means = [10.0, 10.01, 9.99, 10.02, 500.0] + lenient_checker = RobustConsistency(trend_n_sigma=1e6) + assert lenient_checker._reject_outliers(means) == means + + strict_checker = RobustConsistency(trend_n_sigma=5.0) + assert strict_checker._reject_outliers(means) == [ + 10.0, + 10.01, + 9.99, + 10.02, + ] + + def test_vector_valued_drops_whole_candidate_on_any_element_outlier(self): + # A candidate that's fine in one output element but a severe + # outlier in another must still be dropped entirely (not just + # masked in the bad element) -- "badness" is reduced across all + # output dimensions before picking the worst candidate. + checker = RobustConsistency() + means = [ + np.array([10.0, 5.0]), + np.array([10.01, 5.01]), + np.array([9.99, 500.0]), # fine in element 0, an outlier in 1 + ] + trusted = checker._reject_outliers(means) + assert len(trusted) == 2 + assert all( + np.array_equal(t, means[i]) + for t, i in zip(trusted, [0, 1], strict=True) + ) + + def test_nan_candidate_is_never_flagged_as_worst(self): + # Known, documented limitation: `nanargmax` ignores NaNs, so a + # candidate whose mean is entirely NaN can never be selected as + # "the worst" and is left in the trusted set untouched (harmless in + # practice: it doesn't shift `np.nanmean` of the final value, and + # `RobustConsistency.method`'s final blanket `isclose` check against + # a non-NaN blended value still reports `success=False` overall). + checker = RobustConsistency() + means = [10.0, 10.01, np.nan] + trusted = checker._reject_outliers(means) + assert len(trusted) == 3 + assert np.isnan(trusted[-1]) diff --git a/tests/benchmark_models/test_petab_benchmark.py b/tests/benchmark_models/test_petab_benchmark.py index fb58aa8a4b..9d3a4246ef 100644 --- a/tests/benchmark_models/test_petab_benchmark.py +++ b/tests/benchmark_models/test_petab_benchmark.py @@ -21,6 +21,7 @@ import yaml from amici import get_model_root_dir from amici.adapters.fiddy import ( + RobustConsistency, simulate_petab_to_cached_functions, simulate_petab_v2_to_cached_functions, ) @@ -43,7 +44,6 @@ ) from fiddy import MethodId, get_derivative from fiddy.derivative_check import NumpyIsCloseDerivativeCheck -from fiddy.success import Consistency from petab.v1.lint import measurement_table_has_timepoint_specific_mappings from petab.v1.visualize import plot_problem @@ -388,6 +388,13 @@ def test_nominal_parameters_llh(benchmark_problem): # https://github.com/AMICI-dev/AMICI/issues/18 "ignore:Adjoint sensitivity analysis for models with discontinuous " "right hand sides .*:UserWarning", + # RobustConsistency deliberately warns when it rejects a step size that + # was self-consistent on its own but inconsistent with the majority of + # other step sizes -- this is the intended corrective behavior, not a + # test failure (see https://github.com/ICB-DCM/fiddy/pull/77, fixes + # AMICI-dev/AMICI#3078). + "ignore:.*were rejected as inconsistent with the majority of other " + "step sizes.*:UserWarning", ) @pytest.mark.parametrize("scale", (True, False), ids=["scaled", "unscaled"]) @pytest.mark.parametrize( @@ -405,11 +412,6 @@ def test_benchmark_gradient( if not scale and problem_id in ( "Smith_BMCSystBiol2013", "Brannmark_JBC2010", - "Elowitz_Nature2000", - "Borghans_BiophysChem1997", - "Sneyd_PNAS2002", - "Bertozzi_PNAS2020", - "Zheng_PNAS2012", ): # not really worth the effort trying to fix these cases if they # only fail on linear scale @@ -485,7 +487,7 @@ def test_benchmark_gradient( sizes=cur_settings.step_sizes, direction_ids=parameter_ids, method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], - success_checker=Consistency( + success_checker=RobustConsistency( rtol=cur_settings.rtol_consistency, atol=cur_settings.atol_consistency, ), @@ -601,6 +603,13 @@ def write_debug_output( "right hand sides .*:UserWarning", "ignore:.*has `useValuesFromTriggerTime=true'.*:UserWarning", "ignore:.*Using `log-normal` instead.*:UserWarning", + # RobustConsistency deliberately warns when it rejects a step size that + # was self-consistent on its own but inconsistent with the majority of + # other step sizes -- this is the intended corrective behavior, not a + # test failure (see https://github.com/ICB-DCM/fiddy/pull/77, fixes + # AMICI-dev/AMICI#3078). + "ignore:.*were rejected as inconsistent with the majority of other " + "step sizes.*:UserWarning", ) @pytest.mark.parametrize("problem_id", problems_for_llh_check) def test_nominal_parameters_llh_v2(problem_id): @@ -811,7 +820,7 @@ def test_nominal_parameters_llh_v2(problem_id): sizes=cur_settings.step_sizes, direction_ids=parameter_ids, method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], - success_checker=Consistency( + success_checker=RobustConsistency( rtol=cur_settings.rtol_consistency, atol=cur_settings.atol_consistency, ),