diff --git a/.github/workflows/test_benchmark_collection_models.yml b/.github/workflows/test_benchmark_collection_models.yml index 0f6054b1f0..ac4926cf18 100644 --- a/.github/workflows/test_benchmark_collection_models.yml +++ b/.github/workflows/test_benchmark_collection_models.yml @@ -163,7 +163,7 @@ jobs: run: | python3 -m pip uninstall -y petab && python3 -m pip install git+https://github.com/petab-dev/libpetab-python.git@1b8599dd1eb9bda74853255b4cc4baf75b7e4d63 \ && python3 -m pip install -U sympy \ - && python3 -m pip install git+https://github.com/ICB-DCM/fiddy.git@amici100 + && python3 -m pip install git+https://github.com/ICB-DCM/fiddy.git@main - run: pip uninstall -y diffrax && pip install git+https://github.com/patrick-kidger/diffrax@main # TODO FIXME https://github.com/patrick-kidger/diffrax/issues/654 + event dirs diff --git a/.github/workflows/test_windows.yml b/.github/workflows/test_windows.yml index 9b11c2e98a..117139e5d9 100644 --- a/.github/workflows/test_windows.yml +++ b/.github/workflows/test_windows.yml @@ -58,6 +58,11 @@ jobs: shell: bash run: pip install -v $(ls -t dist/amici-*.tar.gz | head -1)[petab,test,jax] + # TODO: switch back to PyPI once the fiddy redesign is released + - name: Install fiddy from GitHub + shell: bash + run: pip install git+https://github.com/ICB-DCM/fiddy.git@main + - run: python -m amici - name: Get Pooch Cache Directory diff --git a/CHANGELOG.md b/CHANGELOG.md index 829f1e1303..c9c357d4e3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,12 @@ See also our [versioning policy](https://amici.readthedocs.io/en/latest/versioni ### v1.1.1 (unreleased) +**Features** + +* The `fiddy` adapter (`amici.adapters.fiddy`) and the `fiddy` package itself + were redesigned, making finite-difference gradient checks much more + robust and requiring few, if any, hyperparameters. + **Fixes** * There are no more reserved names: previously, model import or diff --git a/python/sdist/amici/adapters/fiddy.py b/python/sdist/amici/adapters/fiddy.py index e07ac49106..f2a1453787 100644 --- a/python/sdist/amici/adapters/fiddy.py +++ b/python/sdist/amici/adapters/fiddy.py @@ -10,7 +10,6 @@ from __future__ import annotations -import warnings from collections.abc import Callable from functools import partial from inspect import signature @@ -18,9 +17,7 @@ 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 fiddy import CachedFunction, Type from petab.v1.C import LIN, LOG, LOG10 from amici.sim.sundials import ( @@ -41,228 +38,15 @@ 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", + "run_simulation_to_function_and_derivative", + "simulate_petab_to_function_and_derivative", + "simulate_petab_v2_to_function_and_derivative", + "output_labels_for_derivatives", ] 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 @@ -326,8 +110,73 @@ def _rdata_array_transpose(array: np.ndarray, variable: str) -> tuple[int]: if v not in ["sz", "srz", "ssigmaz", "s2llh"] } +# Entities to id type mapping +_entity_ids_by_variable = { + "x": "state", + "x0": "state", + "x_ss": "state", + "y": "observable", + "sigmay": "observable", + "res": "observable", +} +# Entities that have a time +_has_timepoint_axis = {"x", "y", "sigmay", "res"} -def run_simulation_to_cached_functions( + +def output_labels_for_derivatives( + amici_model: AmiciModel, + derivative_variables: list[str] = None, + timepoints: list[float] = None, +) -> list[str]: + """Per-flat-row labels for fiddy's `function`/`derivative`'s bundled output. + + :param amici_model: The AMICI model (for state/observable IDs). + :param derivative_variables: Same meaning/default as + :func:`run_simulation_to_function_and_derivative`. + :param timepoints: Output timepoints, for variables with a timepoint + axis. Defaults to `amici_model.get_timepoints()`. + :return: One label per flat output row, in bundling order. + :raises NotImplementedError: For a variable with no label source + (``z``, ``rz``, ``sigmaz``, or second-order ``sllh``). + """ + variables = list( + default_derivatives + if derivative_variables is None + else derivative_variables + ) + unsupported = [v for v in variables if v not in default_derivatives] + if unsupported: + raise NotImplementedError( + f"No output labels available for {unsupported} -- only " + f"{list(default_derivatives)} are supported." + ) + if timepoints is None: + timepoints = list(amici_model.get_timepoints()) + ids_by_kind = { + "state": list(amici_model.get_state_ids()), + "observable": list(amici_model.get_observable_ids()), + } + + labels = [] + for variable in variables: + if variable == "llh": + labels.append("llh") + continue + entity_ids = ids_by_kind[_entity_ids_by_variable[variable]] + if variable in _has_timepoint_axis: + labels.extend( + f"{variable}[t={t:g}, id={entity_id}]" + for t in timepoints + for entity_id in entity_ids + ) + else: + labels.extend( + f"{variable}[id={entity_id}]" for entity_id in entity_ids + ) + return labels + + +def run_simulation_to_function_and_derivative( amici_model: AmiciModel, *, cache: bool = True, @@ -336,7 +185,15 @@ def run_simulation_to_cached_functions( amici_edata: AmiciExpData = None, derivative_variables: list[str] = None, ): - """Convert `run_simulation` to fiddy functions. + """Convert `run_simulation` to a fiddy-checkable ``(function, + derivative)`` pair, e.g. for :func:`fiddy.check_jacobian`. + + Both `function` and `derivative` return a dict keyed by + `derivative_variables` (or `default_derivatives`' keys, if not given) + -- one simulation output per key for `function` (`x`, `y`, `llh`, ...), + its forward-sensitivity counterpart for `derivative` (`sx`, `sy`, + `sllh`, ..., with the parameter axis moved last, and sliced/reordered to + `free_parameter_ids` from each simulation's `rdata.plist`. :param amici_model: The AMICI model to simulate. @@ -349,17 +206,18 @@ def run_simulation_to_cached_functions( The variables that derivatives will be computed or approximated for. See the keys of `all_rdata_derivatives` for options. :param free_parameter_ids: - The IDs that correspond to the values in the free parameter vector that is - simulated. + IDs for the values in the simulated free parameter vector. Each + must be in the resolved `plist` (`amici_model` or `amici_edata`), + or `derivative` raises `ValueError`. :param cache: Whether to cache the function calls. - :returns: function, derivatives and structure + :returns: A tuple of `(function, derivative)`. """ if amici_solver is None: amici_solver = amici_model.create_solver() if free_parameter_ids is None: free_parameter_ids = amici_model.get_free_parameter_ids() - if amici_edata is not None and amici_edata.free_parameters is not None: + if amici_edata is not None and amici_edata.free_parameters: raise NotImplementedError( "Customization of parameter values inside AMICI ExpData." ) @@ -368,6 +226,7 @@ def run_simulation_to_cached_functions( chosen_derivatives = { k: all_rdata_derivatives[k] for k in derivative_variables } + amici_free_parameter_ids = amici_model.get_free_parameter_ids() def run_amici_simulation( point: Type.POINT, order: SensitivityOrder @@ -380,84 +239,56 @@ def run_amici_simulation( ) return rdata - def function(point: Type.POINT): + def function(point: Type.POINT) -> dict[str, np.ndarray]: rdata = run_amici_simulation(point=point, order=SensitivityOrder.none) - outputs = { - variable: fiddy_array(getattr(rdata, variable)) - for variable in chosen_derivatives - } - rdata_flat = np.concatenate( - [output.flat for output in outputs.values()] - ) - return rdata_flat - - def derivative(point: Type.POINT, return_dict: bool = False): + outputs = {} + for variable in chosen_derivatives: + value = getattr(rdata, variable) + # AMICI represents a structurally empty field (e.g. `x` for a + # model with zero states) as `None`, not an empty array -- + # `np.asarray(None, dtype=float)` would silently produce a 0-d + # NaN scalar instead, which is both the wrong shape and would + # spuriously fail fiddy's non-finite-value check. + if value is not None: + outputs[variable] = np.asarray(value, dtype=float) + return outputs + + def derivative(point: Type.POINT) -> dict[str, np.ndarray]: rdata = run_amici_simulation(point=point, order=SensitivityOrder.first) - outputs = { - variable: _rdata_array_transpose( - array=fiddy_array(getattr(rdata, derivative_variable)), - variable=derivative_variable, - ) - for variable, derivative_variable in chosen_derivatives.items() - } - rdata_flat = np.concatenate( - [ - output_array.reshape(-1, output_array.shape[-1]) - for output_array in outputs.values() - ], - axis=0, - ) - if return_dict: - return outputs - return rdata_flat + rdata_free_parameter_ids = [ + amici_free_parameter_ids[i] for i in rdata.plist + ] + try: + parameter_indices = [ + rdata_free_parameter_ids.index(parameter_id) + for parameter_id in free_parameter_ids + ] + except ValueError as error: + raise ValueError( + f"{error}. `free_parameter_ids` requested a parameter " + "whose sensitivity was not computed by this simulation " + "-- check `amici_model.get_parameter_list()` and " + "`amici_edata.plist` (if `amici_edata` is given, its own " + "`plist` takes priority over the model's whenever it is " + "non-empty)." + ) from error + outputs = {} + for variable, derivative_variable in chosen_derivatives.items(): + value = getattr(rdata, derivative_variable) + if value is not None: # see `function`'s comment above + outputs[variable] = _rdata_array_transpose( + array=np.asarray(value, dtype=float), + variable=derivative_variable, + )[..., parameter_indices] + return outputs if cache: function = CachedFunction(function) - derivative = CachedFunction(derivative) - # Get structure - dummy_point = fiddy_array( - [ - amici_model.get_free_parameter_by_id(par_id) - for par_id in free_parameter_ids - ] - ) - dummy_rdata = run_amici_simulation( - point=dummy_point, order=SensitivityOrder.first - ) - - structures = { - "function": {variable: None for variable in chosen_derivatives}, - "derivative": {variable: None for variable in chosen_derivatives}, - } - function_position = 0 - derivative_position = 0 - for variable, derivative_variable in chosen_derivatives.items(): - function_array = fiddy_array(getattr(dummy_rdata, variable)) - derivative_array = fiddy_array( - getattr(dummy_rdata, derivative_variable) - ) - structures["function"][variable] = ( - function_position, - function_position + function_array.size, - function_array.shape, - ) - structures["derivative"][variable] = ( - derivative_position, - derivative_position + derivative_array.size, - derivative_array.shape, - ) - function_position += function_array.size - derivative_position += derivative_array.size - - return function, derivative, structures - - -# (start, stop, shape) -TYPE_STRUCTURE = tuple[int, int, tuple[int, ...]] + return function, derivative -def simulate_petab_to_cached_functions( +def simulate_petab_to_function_and_derivative( petab_problem: petab.Problem, *, amici_model: Model, @@ -470,7 +301,8 @@ def simulate_petab_to_cached_functions( ) -> tuple[Type.FUNCTION, Type.FUNCTION]: """ Convert :func:`amici.sim.sundials.petab.v1.simulate_petab` - (PEtab v1 simulations) to fiddy functions. + (PEtab v1 simulations) to a fiddy-checkable ``(function, derivative)`` + pair, e.g. for :func:`fiddy.check_gradient`. Note that all gradients are provided on linear scale. The correction from `'log10'` scale is automatically done. @@ -576,18 +408,18 @@ def derivative(point: Type.POINT) -> Type.POINT: if cache: function = CachedFunction(function) - derivative = CachedFunction(derivative) return function, derivative -def simulate_petab_v2_to_cached_functions( +def simulate_petab_v2_to_function_and_derivative( petab_simulator: PetabSimulator, *, free_parameter_ids: list[str] = None, cache: bool = True, ) -> tuple[Type.FUNCTION, Type.FUNCTION]: - r"""Create fiddy functions for PetabSimulator. + r"""Create a fiddy-checkable ``(function, derivative)`` pair for a + `PetabSimulator`, e.g. for :func:`fiddy.check_gradient`. :param petab_simulator: The PEtab simulator to use. @@ -634,6 +466,5 @@ def derivative(point: Type.POINT) -> Type.POINT: if cache: function = CachedFunction(function) - derivative = CachedFunction(derivative) return function, derivative diff --git a/python/tests/adapters/test_fiddy.py b/python/tests/adapters/test_fiddy.py index b60c042799..0febfb11d8 100644 --- a/python/tests/adapters/test_fiddy.py +++ b/python/tests/adapters/test_fiddy.py @@ -1,26 +1,27 @@ """Tests for `amici.adapters.fiddy`.""" +import sys from pathlib import Path +import amici import numpy as np import pytest from amici.adapters.fiddy import ( - RobustConsistency, - run_simulation_to_cached_functions, - simulate_petab_to_cached_functions, + output_labels_for_derivatives, + run_simulation_to_function_and_derivative, + simulate_petab_to_function_and_derivative, ) from amici.importers.petab.v1 import import_petab_problem from amici.sim.sundials import SensitivityOrder, SteadyStateSensitivityMode -from fiddy import MethodId, Type, get_derivative -from fiddy.derivative_check import NumpyIsCloseDerivativeCheck -from fiddy.directional_derivative import ComputerResult -from numpy.testing import assert_allclose +from fiddy import ( + JoblibExecutor, + SequentialExecutor, + Type, + check_gradient, + check_jacobian, + estimate_gradient, +) from petab import v1 -from scipy.optimize import rosen - -# Absolute and relative tolerances for finite difference gradient checks. -ATOL: float = 1e-3 -RTOL: float = 1e-3 def lotka_volterra() -> tuple[v1.Problem, np.ndarray]: @@ -37,11 +38,23 @@ def lotka_volterra() -> tuple[v1.Problem, np.ndarray]: return petab_problem, point -@pytest.mark.parametrize("problem_generator", [lotka_volterra]) -def test_run_amici_simulation_to_functions(problem_generator): - petab_problem, point = problem_generator() +@pytest.fixture(scope="session") +def lotka_volterra_model_module(): + """Imports `lotka_volterra` model module.""" + petab_problem, _ = lotka_volterra() + import_petab_problem(petab_problem) + model_name = petab_problem.model.model_id + return amici.import_model_module( + model_name, amici.get_model_dir(model_name) + ) + + +def test_run_amici_simulation_to_function_and_derivative( + lotka_volterra_model_module, +): + petab_problem, point = lotka_volterra() timepoints = sorted(set(petab_problem.measurement_df.time)) - amici_model = import_petab_problem(petab_problem) + amici_model = lotka_volterra_model_module.get_model() amici_model.set_timepoints(timepoints) amici_solver = amici_model.create_solver() @@ -52,61 +65,129 @@ def test_run_amici_simulation_to_functions(problem_generator): petab_problem.parameter_df.estimate == 1 ].index ) - parameter_indices = [ - amici_model.get_free_parameter_ids().index(parameter_id) - for parameter_id in parameter_ids - ] - - ( - amici_function, - amici_derivative, - structures, - ) = run_simulation_to_cached_functions( + + derivative_variables = ["x", "x0", "y", "sigmay"] + function, derivative = run_simulation_to_function_and_derivative( free_parameter_ids=parameter_ids, amici_model=amici_model, amici_solver=amici_solver, + derivative_variables=derivative_variables, + ) + + expected = derivative(point) + output_labels = output_labels_for_derivatives( + amici_model, + derivative_variables=derivative_variables, + timepoints=timepoints, + ) + result = check_jacobian( + function, + point, + expected, + direction_labels=parameter_ids, + output_labels=output_labels, + ) + assert len(output_labels) == len(result.output_results) + result.assert_success(always_print=True) + + +def test_run_simulation_respects_a_customized_parameter_list( + lotka_volterra_model_module, +): + """Regression test for proper `plist` handling.""" + petab_problem, _ = lotka_volterra() + timepoints = sorted(set(petab_problem.measurement_df.time)) + amici_model = lotka_volterra_model_module.get_model() + amici_model.set_timepoints(timepoints) + amici_solver = amici_model.create_solver() + amici_solver.set_sensitivity_order(SensitivityOrder.first) + + free_parameter_ids = list(amici_model.get_free_parameter_ids()) + alpha_id = "alpha" + sigma_id = "noiseParameter1_observable_prey" + assert set(free_parameter_ids) == {alpha_id, "gamma", sigma_id} + alpha_index = free_parameter_ids.index(alpha_id) + sigma_index = free_parameter_ids.index(sigma_id) + # Confirm the natural order actually puts alpha before sigma, so the + # `[sigma_index, alpha_index]` plist below is really a reversal. + assert alpha_index < sigma_index + + values = {alpha_id: 2.0, "gamma": 3.0, sigma_id: 1.0} + full_point = np.array([values[pid] for pid in free_parameter_ids]) + + # Baseline: default (identity) plist + _, baseline_derivative = run_simulation_to_function_and_derivative( + free_parameter_ids=free_parameter_ids, + amici_model=amici_model, + amici_solver=amici_solver, + derivative_variables=["llh"], + ) + baseline = baseline_derivative(full_point)["llh"] + + # A subset, reversed relative to the model's own natural order. + amici_model.set_parameter_list([sigma_index, alpha_index]) + subset_ids = [alpha_id, sigma_id] + subset_point = np.array([full_point[alpha_index], full_point[sigma_index]]) + _, subset_derivative = run_simulation_to_function_and_derivative( + free_parameter_ids=subset_ids, + amici_model=amici_model, + amici_solver=amici_solver, + derivative_variables=["llh"], ) + restricted = subset_derivative(subset_point)["llh"] - expected_derivative = amici_derivative(point)[..., parameter_indices] - - derivative = get_derivative( - function=amici_function, - point=point, - sizes=[1e-10, 1e-5], - direction_ids=parameter_ids, - method_ids=[MethodId.FORWARD, MethodId.BACKWARD, MethodId.CENTRAL], - # analysis_classes=[], - # analysis_classes=[ - # lambda: TransformByDirectionScale(scales=parameter_scales), - # ], - success_checker=RobustConsistency(atol=1e-2), + np.testing.assert_allclose(restricted[0], baseline[alpha_index]) + np.testing.assert_allclose(restricted[1], baseline[sigma_index]) + + +@pytest.mark.skipif( + sys.platform == "win32", + reason="Parallelization/pickling requires HDF5 support -- unavailable on Windows builds.", +) +def test_joblib_executor_agrees_with_sequential_executor( + lotka_volterra_model_module, +): + """Results from `SequentialExecutor` and `JoblibExecutor` + must agree exactly. + """ + petab_problem, point = lotka_volterra() + timepoints = sorted(set(petab_problem.measurement_df.time)) + amici_model = lotka_volterra_model_module.get_model() + amici_model.set_timepoints(timepoints) + amici_solver = amici_model.create_solver() + amici_solver.set_sensitivity_order(SensitivityOrder.first) + + parameter_ids = list( + petab_problem.parameter_df[ + petab_problem.parameter_df.estimate == 1 + ].index ) - test_derivative = derivative.value - - # The test derivative is close to the expected derivative. - assert_allclose( - test_derivative, - expected_derivative, - rtol=1e-1, - atol=1e-1, - equal_nan=True, + + function, _ = run_simulation_to_function_and_derivative( + free_parameter_ids=parameter_ids, + amici_model=amici_model, + amici_solver=amici_solver, + derivative_variables=["x", "x0", "y", "sigmay"], ) - # Same as above assert. - check = NumpyIsCloseDerivativeCheck( - derivative=derivative, - expectation=expected_derivative, - point=point, + sequential = estimate_gradient( + function, point, executor=SequentialExecutor() ) - result = check(rtol=1e-1, atol=1e-1, equal_nan=True) - assert result.success + parallel = estimate_gradient( + function, point, executor=JoblibExecutor(n_jobs=4) + ) + + for s, p in zip(sequential, parallel, strict=True): + assert s.value == p.value + assert s.status == p.status -@pytest.mark.parametrize("problem_generator", [lotka_volterra]) @pytest.mark.parametrize("scaled_parameters", (False, True)) -def test_simulate_petab_to_functions(problem_generator, scaled_parameters): - petab_problem, point = problem_generator() - amici_model = import_petab_problem(petab_problem) +def test_simulate_petab_to_function_and_derivative( + scaled_parameters, lotka_volterra_model_module +): + petab_problem, point = lotka_volterra() + amici_model = lotka_volterra_model_module.get_model() amici_solver = amici_model.create_solver() if amici_model.get_name() == "simple": @@ -131,7 +212,7 @@ def test_simulate_petab_to_functions(problem_generator, scaled_parameters): ) ) - amici_function, amici_derivative = simulate_petab_to_cached_functions( + function, derivative = simulate_petab_to_function_and_derivative( free_parameter_ids=petab_problem.parameter_df.index, petab_problem=petab_problem, amici_model=amici_model, @@ -140,275 +221,6 @@ def test_simulate_petab_to_functions(problem_generator, scaled_parameters): scaled_parameters=scaled_parameters, ) - expected_derivative = amici_derivative(point) - - free_parameter_ids = list( - petab_problem.parameter_df[ - petab_problem.parameter_df.estimate == 1 - ].index - ) - # parameter_scales = dict( - # petab_problem.parameter_df[ - # petab_problem.parameter_df.estimate == 1 - # ].parameterScale - # ) - - derivative = get_derivative( - function=amici_function, - point=point, - sizes=[1e-10, 1e-5, 1e-3, 1e-1], - direction_ids=free_parameter_ids, - method_ids=[MethodId.FORWARD, MethodId.BACKWARD, MethodId.CENTRAL], - success_checker=RobustConsistency(), - ) - - check = NumpyIsCloseDerivativeCheck( - derivative=derivative, - expectation=expected_derivative, - point=point, - ) - 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]) + expected = derivative(point) + result = check_gradient(function, point, expected) + result.assert_success(always_print=True) diff --git a/scripts/installAmiciSource.sh b/scripts/installAmiciSource.sh index 7881eecc2a..8b99f79bef 100755 --- a/scripts/installAmiciSource.sh +++ b/scripts/installAmiciSource.sh @@ -45,4 +45,6 @@ python -m pip uninstall petab -y python -m pip install git+https://github.com/petab-dev/libpetab-python.git@main AMICI_BUILD_TEMP="${AMICI_PATH}/python/sdist/build/temp" \ python -m pip install --verbose -e "${AMICI_PATH}/python/sdist[petab,test,vis,jax]" --no-build-isolation +# TODO: switch back to PyPI once the fiddy redesign is released +python -m pip install git+https://github.com/ICB-DCM/fiddy.git@main deactivate diff --git a/swig/solver.i b/swig/solver.i index 4c04cefdb1..dfbc9b8805 100644 --- a/swig/solver.i +++ b/swig/solver.i @@ -109,6 +109,13 @@ def _solver_reduce(self: "Solver"): reboots and will not work in distributed (MPI) settings. This requires that amici was compiled with HDF5 support. """ + from amici.sim.sundials import hdf5_enabled + if not hdf5_enabled: + raise RuntimeError( + "Cannot pickle this Solver: pickling stores solver settings " + "via HDF5, but this AMICI installation was built without " + "HDF5 support." + ) from amici.sim.sundials._swig_wrappers import restore_solver, write_solver_settings_to_hdf5 from tempfile import NamedTemporaryFile import os diff --git a/tests/benchmark_models/test_petab_benchmark.py b/tests/benchmark_models/test_petab_benchmark.py index 9d3a4246ef..9759af307a 100644 --- a/tests/benchmark_models/test_petab_benchmark.py +++ b/tests/benchmark_models/test_petab_benchmark.py @@ -9,11 +9,10 @@ import logging import os from collections import defaultdict -from dataclasses import dataclass, field +from dataclasses import dataclass from pathlib import Path import benchmark_models_petab -import fiddy import numpy as np import pandas as pd import petab.v1 as petab @@ -21,9 +20,8 @@ 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, + simulate_petab_to_function_and_derivative, + simulate_petab_v2_to_function_and_derivative, ) from amici.importers.petab.v1 import ( import_petab_problem, @@ -42,8 +40,7 @@ rdatas_to_measurement_df, simulate_petab, ) -from fiddy import MethodId, get_derivative -from fiddy.derivative_check import NumpyIsCloseDerivativeCheck +from fiddy import check_gradient from petab.v1.lint import measurement_table_has_timepoint_specific_mappings from petab.v1.visualize import plot_problem @@ -153,26 +150,6 @@ class GradientCheckSettings: # Absolute and relative tolerances for simulation atol_sim: float = 1e-16 rtol_sim: float = 1e-12 - # Absolute and relative tolerances for finite difference gradient checks. - atol_check: float = 1e-3 - rtol_check: float = 1e-2 - # Absolute and relative tolerances for fiddy consistency check between - # forward/backward/central differences. - atol_consistency: float = 1e-5 - rtol_consistency: float = 1e-1 - # Step sizes for finite difference gradient checks. - step_sizes: list[float] = field( - default_factory=lambda: [ - 2e-1, - 1e-1, - 5e-2, - 1e-2, - 5e-1, - 1e-3, - 1e-4, - 1e-5, - ] - ) rng_seed: int = 0 ss_computation_mode: SteadyStateComputationMode = ( SteadyStateComputationMode.integrationOnly @@ -186,39 +163,26 @@ class GradientCheckSettings: settings = defaultdict(GradientCheckSettings) # NOTE: Newton method fails badly with ASA for Blasi_CellSystems2016 settings["Blasi_CellSystems2016"] = GradientCheckSettings( - atol_check=1e-12, - rtol_check=1e-4, ss_sensitivity_mode=SteadyStateSensitivityMode.integrationOnly, ) -settings["Borghans_BiophysChem1997"] = GradientCheckSettings( - rng_seed=2, - atol_check=1e-5, - rtol_check=1e-3, -) settings["Brannmark_JBC2010"] = GradientCheckSettings( + rtol_sim=1e-14, ss_sensitivity_mode=SteadyStateSensitivityMode.integrationOnly, ) -settings["Fujita_SciSignal2010"] = GradientCheckSettings( - atol_check=1e-7, - rtol_check=5e-4, -) -settings["Giordano_Nature2020"] = GradientCheckSettings( - atol_check=1e-6, rtol_check=1e-3, rng_seed=1 -) +settings["Giordano_Nature2020"] = GradientCheckSettings(rng_seed=1) settings["Okuonghae_ChaosSolitonsFractals2020"] = GradientCheckSettings( atol_sim=1e-14, rtol_sim=1e-14, noise_level=0.01, - atol_consistency=1e-3, ) settings["Oliveira_NatCommun2021"] = GradientCheckSettings( # Avoid "root after reinitialization" atol_sim=1e-12, rtol_sim=1e-12, ) -settings["Raia_CancerResearch2011"] = GradientCheckSettings( - atol_check=1e-10, - rtol_check=1e-3, +settings["SalazarCavazos_MBoC2020"] = GradientCheckSettings( + atol_sim=1e-12, + rtol_sim=1e-12, ) settings["Smith_BMCSystBiol2013"] = GradientCheckSettings( atol_sim=1e-10, @@ -227,38 +191,34 @@ class GradientCheckSettings: settings["Sneyd_PNAS2002"] = GradientCheckSettings( atol_sim=1e-15, rtol_sim=1e-12, - atol_check=1e-5, - rtol_check=1e-4, rng_seed=7, ) settings["Weber_BMC2015"] = GradientCheckSettings( - atol_sim=1e-12, - rtol_sim=1e-12, - atol_check=1e-6, - rtol_check=1e-2, - rng_seed=2, + atol_sim=1e-13, + rtol_sim=1e-13, + rng_seed=1, ) settings["Zheng_PNAS2012"] = GradientCheckSettings( - rng_seed=1, + rng_seed=2, rtol_sim=1e-15, - atol_check=5e-4, - rtol_check=4e-3, noise_level=0.01, ss_sensitivity_mode=SteadyStateSensitivityMode.integrationOnly, - step_sizes=[ - 3e-1, - 2e-1, - 1e-1, - 5e-2, - 1e-2, - 5e-1, - 1e-3, - 1e-4, - 1e-5, - ], ) +def assert_gradient_check_confirms_something(result) -> None: + """`check_gradient`'s `success` is `True` as long as no direction is + confidently *wrong* -- a check where every direction came back + "inconclusive" (noise-dominated/discontinuity-suspected) would still + report success, having actually confirmed nothing. Require at least + one direction to have been confirmed converged. + """ + assert any(r.outcome == "passed" for r in result.direction_results), ( + "check_gradient reported success, but every direction was " + "inconclusive -- nothing was actually confirmed correct." + ) + + @pytest.mark.filterwarnings( "ignore:divide by zero encountered in log", # https://github.com/AMICI-dev/AMICI/issues/18 @@ -388,13 +348,6 @@ 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( @@ -402,21 +355,11 @@ def test_nominal_parameters_llh(benchmark_problem): (SensitivityMethod.forward, SensitivityMethod.adjoint), ids=["forward", "adjoint"], ) -def test_benchmark_gradient( - benchmark_problem, scale, sensitivity_method, request -): +def test_benchmark_gradient(benchmark_problem, scale, sensitivity_method): problem_id, petab_problem, _, amici_model = benchmark_problem if problem_id not in problems_for_gradient_check: pytest.skip("Excluded from gradient check.") - if not scale and problem_id in ( - "Smith_BMCSystBiol2013", - "Brannmark_JBC2010", - ): - # not really worth the effort trying to fix these cases if they - # only fail on linear scale - pytest.skip("scale=False disabled for this problem") - petab_problem = benchmark_models_petab.get_problem(problem_id) if measurement_table_has_timepoint_specific_mappings( petab_problem.measurement_df, @@ -442,22 +385,24 @@ def test_benchmark_gradient( cur_settings.ss_sensitivity_mode ) - amici_function, amici_derivative = simulate_petab_to_cached_functions( - petab_problem=petab_problem, - free_parameter_ids=parameter_ids, - amici_model=amici_model, - solver=amici_solver, - scaled_parameters=scale, - scaled_gradients=scale, - # FIXME: there is some issue with caching in fiddy - # e.g. Elowitz_Nature2000-True fails with cache=True, - # but not with cache=False - # cache=not debug, - cache=False, - num_threads=os.cpu_count(), + amici_function, amici_derivative = ( + simulate_petab_to_function_and_derivative( + petab_problem=petab_problem, + free_parameter_ids=parameter_ids, + amici_model=amici_model, + solver=amici_solver, + scaled_parameters=scale, + scaled_gradients=scale, + num_threads=os.cpu_count(), + ) ) np.random.seed(cur_settings.rng_seed) + bounds = ( + np.array(petab_problem.get_lb(free=True, fixed=False, scaled=scale)), + np.array(petab_problem.get_ub(free=True, fixed=False, scaled=scale)), + ) + # find a point where the derivative can be computed for _ in range(5): if scale: @@ -471,6 +416,9 @@ def test_benchmark_gradient( np.random.randn(len(point)) * point * cur_settings.noise_level ) point += point_noise # avoid small gradients at nominal value + # Jittering can push a point outside its own bounds; clip before + # passing both to fiddy. + point = np.clip(point, bounds[0], bounds[1]) try: expected_derivative = amici_derivative(point) @@ -481,119 +429,15 @@ def test_benchmark_gradient( else: raise RuntimeError("Could not compute expected derivative.") - derivative = get_derivative( - function=amici_function, - point=point, - sizes=cur_settings.step_sizes, - direction_ids=parameter_ids, - method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], - success_checker=RobustConsistency( - rtol=cur_settings.rtol_consistency, - atol=cur_settings.atol_consistency, - ), - expected_result=expected_derivative, - relative_sizes=not scale, - ) - print() print("Testing at:", point) print("Expected derivative (amici):", expected_derivative) - print("Print actual derivative (fiddy):", derivative.series.values) - - if debug: - write_debug_output( - debug_path / f"{request.node.callspec.id}.tsv", - derivative, - expected_derivative, - parameter_ids, - ) - - assert_gradient_check_success( - derivative, - expected_derivative, - point, - rtol=cur_settings.rtol_check, - atol=cur_settings.atol_check, - always_print=True, - ) - - -def assert_gradient_check_success( - derivative: fiddy.Derivative, - expected_derivative: np.ndarray, - point: np.ndarray, - atol: float, - rtol: float, - always_print: bool = False, -) -> None: - if not derivative.df.success.all(): - raise AssertionError( - f"Failed to compute finite differences:\n{derivative.df}" - ) - check = NumpyIsCloseDerivativeCheck( - derivative=derivative, - expectation=expected_derivative, - point=point, - ) - check_result = check(rtol=rtol, atol=atol) - - if check_result.success is True and not always_print: - return - df = check_result.df - df["abs_diff"] = np.abs(df["expectation"] - df["test"]) - df["rel_diff"] = df["abs_diff"] / np.abs(df["expectation"]) - df["atol_success"] = df["abs_diff"] <= atol - df["rtol_success"] = df["rel_diff"] <= rtol - max_adiff = df["abs_diff"].max() - max_rdiff = df["rel_diff"].max() - - success_fail = "succeeded" if check_result.success else "failed" - with pd.option_context( - "display.max_columns", - None, - "display.width", - None, - "display.max_rows", - None, - ): - message = ( - f"Gradient check {success_fail}:\n{df}\n\n" - f"Maximum absolute difference: {max_adiff} (tolerance: {atol})\n" - f"Maximum relative difference: {max_rdiff} (tolerance: {rtol})" - ) - - if check_result.success is False: - raise AssertionError(message) - - if always_print: - print(message) - - -def write_debug_output( - file_name, derivative, expected_derivative, parameter_ids -): - df = pd.DataFrame( - [ - { - ( - "fd", - r.metadata["size_absolute"], - str(r.method_id), - ): r.value - for c in d.computers - for r in c.results - } - for d in derivative.directional_derivatives - ], - index=parameter_ids, + result = check_gradient( + amici_function, point, expected_derivative, bounds=bounds ) - df[("fd", "full", "")] = derivative.series.values - df[("amici", "", "")] = expected_derivative - df["abs_diff"] = np.abs(df[("fd", "full", "")] - df[("amici", "", "")]) - df["rel_diff"] = df["abs_diff"] / np.abs(df[("amici", "", "")]) - - df.to_csv(file_name, sep="\t") + result.assert_success(always_print=True) + assert_gradient_check_confirms_something(result) @pytest.mark.filterwarnings( @@ -603,13 +447,6 @@ 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): @@ -753,20 +590,6 @@ def test_nominal_parameters_llh_v2(problem_id): # TODO scale = False - # also excluded from v1 test - 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 - pytest.skip("scale=False disabled for this problem") - cur_settings = settings[problem_id] ps.solver.set_absolute_tolerance(cur_settings.atol_sim) ps.solver.set_relative_tolerance(cur_settings.rtol_sim) @@ -784,13 +607,19 @@ def test_nominal_parameters_llh_v2(problem_id): ) parameter_ids = ps._petab_problem.x_free_ids - amici_function, amici_derivative = simulate_petab_v2_to_cached_functions( - ps, - free_parameter_ids=parameter_ids, - cache=False, + amici_function, amici_derivative = ( + simulate_petab_v2_to_function_and_derivative( + ps, + free_parameter_ids=parameter_ids, + ) ) np.random.seed(cur_settings.rng_seed) + bounds = ( + np.array(ps._petab_problem.get_lb(free=True, fixed=False)), + np.array(ps._petab_problem.get_ub(free=True, fixed=False)), + ) + # find a point where the derivative can be computed for _ in range(5): if scale: @@ -804,6 +633,9 @@ def test_nominal_parameters_llh_v2(problem_id): np.random.randn(len(point)) * point * cur_settings.noise_level ) point += point_noise # avoid small gradients at nominal value + # Jittering can push a point outside its own bounds; clip before + # passing both to fiddy. + point = np.clip(point, bounds[0], bounds[1]) try: expected_derivative = amici_derivative(point) @@ -814,41 +646,15 @@ def test_nominal_parameters_llh_v2(problem_id): else: raise RuntimeError("Could not compute expected derivative.") - derivative = get_derivative( - function=amici_function, - point=point, - sizes=cur_settings.step_sizes, - direction_ids=parameter_ids, - method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], - success_checker=RobustConsistency( - rtol=cur_settings.rtol_consistency, - atol=cur_settings.atol_consistency, - ), - expected_result=expected_derivative, - relative_sizes=not scale, - ) - print() print("Testing at:", point) print("Expected derivative (amici):", expected_derivative) - print("Print actual derivative (fiddy):", derivative.series.values) - - # if debug: - # write_debug_output( - # debug_path / f"{request.node.callspec.id}.tsv", - # derivative, - # expected_derivative, - # parameter_ids, - # ) - - assert_gradient_check_success( - derivative, - expected_derivative, - point, - rtol=cur_settings.rtol_check, - atol=cur_settings.atol_check, - always_print=True, + + result = check_gradient( + amici_function, point, expected_derivative, bounds=bounds ) + result.assert_success(always_print=True) + assert_gradient_check_confirms_something(result) def compare_to_reference(problem_id: str, llh: float):