diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py index 0e48bf3af8..ba11687929 100644 --- a/python/sdist/amici/de_model.py +++ b/python/sdist/amici/de_model.py @@ -1617,15 +1617,12 @@ def _compute_equation(self, name: str) -> None: ] if name == "dzdx": for ie in range(self.num_events()): - dtaudx = ( - -self.eq("drootdx")[ie, :] - / self.eq("drootdt_total")[ie] - ) + dtaudx = self.eq("dtaudx") for iz in range(self.num_eventobs()): if ie != self._z2event[iz] - 1: continue dzdt = sp.diff(self.eq("z")[ie][iz], time_symbol) - self._eqs[name][ie][iz, :] += dzdt * dtaudx + self._eqs[name][ie][iz, :] += dzdt * -dtaudx[ie] elif name in ["rz", "drzdx", "drzdp"]: eq_events = [] @@ -1644,12 +1641,24 @@ def _compute_equation(self, name: str) -> None: elif name == "stau": self._eqs[name] = [ - -self.eq("sroot")[ie, :] / self.eq("drootdt_total")[ie] + self.eq("sroot")[ie, :] / self.eq("drootdt_total")[ie] if not self.eq("drootdt_total")[ie].is_zero else sp.zeros(*self.eq("sroot")[ie, :].shape) for ie in range(self.num_events()) ] + elif name == "dtaudx": + self._eqs[name] = [ + self.eq("drootdx")[ie, :] / self.eq("drootdt_total")[ie] + for ie in range(self.num_events()) + ] + + elif name == "dtaudp": + self._eqs[name] = [ + self.eq("drootdp")[ie, :] / self.eq("drootdt_total")[ie] + for ie in range(self.num_events()) + ] + elif name == "deltasx": if self.num_states_solver() * self.num_par() == 0: self._eqs[name] = [] @@ -1665,7 +1674,7 @@ def _compute_equation(self, name: str) -> None: self.eq("stau")[ie] ) and not smart_is_zero_matrix(self.eq("xdot")): tmp_eq += smart_multiply( - self.sym("xdot_old") - self.sym("xdot"), + self.sym("xdot") - self.sym("xdot_old"), self.sym("stau").T, ) @@ -1682,12 +1691,14 @@ def _compute_equation(self, name: str) -> None: if not smart_is_zero_matrix(self.eq("stau")[ie]): # chain rule for the time point tmp_eq += smart_multiply( - self.eq("ddeltaxdt")[ie], self.sym("stau").T + self.eq("ddeltaxdt")[ie], + -self.sym("stau").T, ) # additional part of chain rule state variables tmp_dxdp += smart_multiply( - self.sym("xdot_old"), self.sym("stau").T + self.sym("xdot_old"), + -self.sym("stau").T, ) # finish chain rule for the state variables @@ -1695,6 +1706,11 @@ def _compute_equation(self, name: str) -> None: self.eq("ddeltaxdx")[ie], tmp_dxdp ) + else: + tmp_eq = smart_multiply( + self.sym("xdot") - self.sym("xdot_old"), + self.eq("stau")[ie], + ) event_eqs.append(tmp_eq) self._eqs[name] = event_eqs diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index 1a0dc782db..793cade3e2 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -471,8 +471,8 @@ def _parse_piecewise_to_heaviside(args: Iterable[sp.Expr]) -> sp.Expr: symbolic expressions for arguments of the piecewise function """ # how many condition-expression pairs will we have? - formula = sp.Float(0.0) - not_condition = sp.Float(1.0) + formula = sp.Integer(0) + not_condition = sp.Integer(1) if all(isinstance(arg, ExprCondPair) for arg in args): # sympy piecewise @@ -483,7 +483,7 @@ def _parse_piecewise_to_heaviside(args: Iterable[sp.Expr]) -> sp.Expr: for coeff, trigger in grouped_args: if isinstance(coeff, BooleanAtom): - coeff = sp.Float(int(bool(coeff))) + coeff = sp.Integer(int(bool(coeff))) if trigger == sp.true: return formula + coeff * not_condition @@ -493,7 +493,7 @@ def _parse_piecewise_to_heaviside(args: Iterable[sp.Expr]) -> sp.Expr: tmp = _parse_heaviside_trigger(trigger) formula += coeff * sp.simplify(not_condition * tmp) - not_condition *= 1 - tmp + not_condition *= sp.Integer(1) - tmp return formula @@ -517,21 +517,24 @@ def _parse_heaviside_trigger(trigger: sp.Expr) -> sp.Expr: # step with H(0) = 1 if isinstance(trigger, sp.core.relational.StrictLessThan): # x < y => x - y < 0 => r < 0 - return 1 - sp.Heaviside(root) + return sp.Integer(1) - sp.Heaviside(root) if isinstance(trigger, sp.core.relational.LessThan): # x <= y => not(y < x) => not(y - x < 0) => not -r < 0 return sp.Heaviside(-root) if isinstance(trigger, sp.core.relational.StrictGreaterThan): # y > x => y - x < 0 => -r < 0 - return 1 - sp.Heaviside(-root) + return sp.Integer(1) - sp.Heaviside(-root) if isinstance(trigger, sp.core.relational.GreaterThan): # y >= x => not(x < y) => not(x - y < 0) => not r < 0 return sp.Heaviside(root) # or(x,y) = not(and(not(x),not(y)) if isinstance(trigger, sp.Or): - return 1 - sp.Mul( - *[1 - _parse_heaviside_trigger(arg) for arg in trigger.args] + return sp.Integer(1) - sp.Mul( + *[ + sp.Integer(1) - _parse_heaviside_trigger(arg) + for arg in trigger.args + ] ) if isinstance(trigger, sp.And): diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py index 2dc650d1b8..4892100877 100644 --- a/tests/benchmark-models/test_petab_benchmark.py +++ b/tests/benchmark-models/test_petab_benchmark.py @@ -19,6 +19,7 @@ from fiddy.extensions.amici import simulate_petab_to_cached_functions from fiddy.success import Consistency + repo_root = Path(__file__).parent.parent.parent # reuse compiled models from test_benchmark_collection.sh @@ -40,15 +41,9 @@ "Fujita_SciSignal2010", # FIXME: re-enable "Raia_CancerResearch2011", - "Zheng_PNAS2012", } models = list(sorted(models)) -debug = False -if debug: - debug_path = Path(__file__).parent / "debug" - debug_path.mkdir(exist_ok=True, parents=True) - @dataclass class GradientCheckSettings: @@ -58,6 +53,10 @@ class GradientCheckSettings: # 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 = [ 1e-1, @@ -68,50 +67,73 @@ class GradientCheckSettings: 1e-5, ] rng_seed: int = 0 + ss_sensitivity_mode: amici.SteadyStateSensitivityMode = ( + amici.SteadyStateSensitivityMode.integrateIfNewtonFails + ) + noise_level: float = 0.05 settings = defaultdict(GradientCheckSettings) -settings["Smith_BMCSystBiol2013"] = GradientCheckSettings( - atol_sim=1e-10, - rtol_sim=1e-10, +# 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=amici.SteadyStateSensitivityMode.integrationOnly, +) +settings["Borghans_BiophysChem1997"] = GradientCheckSettings( + rng_seed=2, + atol_check=1e-5, + rtol_check=1e-3, +) +settings["Brannmark_JBC2010"] = GradientCheckSettings( + ss_sensitivity_mode=amici.SteadyStateSensitivityMode.integrationOnly, ) +settings["Giordano_Nature2020"] = GradientCheckSettings( + atol_check=1e-6, rtol_check=1e-3, rng_seed=1 +) +settings["Okuonghae_ChaosSolitonsFractals2020"] = GradientCheckSettings( + atol_sim=1e-14, + rtol_sim=1e-14, + noise_level=0.01, + atol_consistency=1e-3, +) +settings["Okuonghae_ChaosSolitonsFractals2020"].step_sizes.extend([0.2, 0.005]) settings["Oliveira_NatCommun2021"] = GradientCheckSettings( # Avoid "root after reinitialization" atol_sim=1e-12, rtol_sim=1e-12, ) -settings["Okuonghae_ChaosSolitonsFractals2020"] = GradientCheckSettings( - atol_sim=1e-14, - rtol_sim=1e-14, +settings["Smith_BMCSystBiol2013"] = GradientCheckSettings( + atol_sim=1e-10, + rtol_sim=1e-10, ) -settings["Okuonghae_ChaosSolitonsFractals2020"].step_sizes.insert(0, 0.2) -settings["Giordano_Nature2020"] = GradientCheckSettings( - atol_check=1e-6, rtol_check=1e-3, rng_seed=1 +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=1, +) +settings["Zheng_PNAS2012"] = GradientCheckSettings( + rng_seed=1, + rtol_sim=1e-15, + atol_check=5e-4, + rtol_check=4e-3, + noise_level=0.01, ) -def assert_gradient_check_success( - derivative: fiddy.Derivative, - expected_derivative: np.array, - point: np.array, - atol: float, - rtol: float, -) -> 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: - return - - raise AssertionError(f"Gradient check failed:\n{check_result.df}") +debug = False +if debug: + debug_path = Path(__file__).parent / "debug" + debug_path.mkdir(exist_ok=True, parents=True) @pytest.mark.filterwarnings( @@ -135,7 +157,7 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): "Borghans_BiophysChem1997", "Sneyd_PNAS2002", "Bertozzi_PNAS2020", - "Okuonghae_ChaosSolitonsFractals2020", + "Zheng_PNAS2012", ): # not really worth the effort trying to fix these cases if they # only fail on linear scale @@ -144,7 +166,6 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): if ( model in ( - "Blasi_CellSystems2016", # events with parameter-dependent triggers # https://github.com/AMICI-dev/AMICI/issues/18 "Oliveira_NatCommun2021", @@ -154,25 +175,11 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): # FIXME pytest.skip() - if ( - model - in ( - "Weber_BMC2015", - "Sneyd_PNAS2002", - ) - and sensitivity_method == SensitivityMethod.forward - ): - # FIXME - pytest.skip() - petab_problem = benchmark_models_petab.get_problem(model) petab.flatten_timepoint_specific_output_overrides(petab_problem) # Only compute gradient for estimated parameters. - parameter_df_free = petab_problem.parameter_df.loc[ - petab_problem.x_free_ids - ] - parameter_ids = list(parameter_df_free.index) + parameter_ids = petab_problem.x_free_ids cur_settings = settings[model] # Setup AMICI objects. @@ -183,13 +190,10 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): amici_solver = amici_model.getSolver() amici_solver.setAbsoluteTolerance(cur_settings.atol_sim) amici_solver.setRelativeTolerance(cur_settings.rtol_sim) - amici_solver.setMaxSteps(int(1e5)) + amici_solver.setMaxSteps(2 * 10**5) amici_solver.setSensitivityMethod(sensitivity_method) - - if model in ("Brannmark_JBC2010",): - amici_model.setSteadyStateSensitivityMode( - amici.SteadyStateSensitivityMode.integrationOnly - ) + # TODO: we should probably test all sensitivity modes + amici_model.setSteadyStateSensitivityMode(cur_settings.ss_sensitivity_mode) amici_function, amici_derivative = simulate_petab_to_cached_functions( petab_problem=petab_problem, @@ -204,24 +208,20 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): # cache=not debug, cache=False, ) - - noise_level = 0.1 np.random.seed(cur_settings.rng_seed) # find a point where the derivative can be computed for _ in range(5): if scale: - point = np.asarray( - list( - petab_problem.scale_parameters( - dict(parameter_df_free.nominalValue) - ).values() - ) + point = petab_problem.x_nominal_free_scaled + point_noise = ( + np.random.randn(len(point)) * cur_settings.noise_level ) - point_noise = np.random.randn(len(point)) * noise_level else: - point = parameter_df_free.nominalValue.values - point_noise = np.random.randn(len(point)) * point * noise_level + point = petab_problem.x_nominal_free + point_noise = ( + np.random.randn(len(point)) * point * cur_settings.noise_level + ) point += point_noise # avoid small gradients at nominal value try: @@ -239,11 +239,19 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): sizes=cur_settings.step_sizes, direction_ids=parameter_ids, method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], - success_checker=Consistency(atol=1e-5, rtol=1e-1), + success_checker=Consistency( + 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:", expected_derivative) + print("Print actual derivative:", derivative.series.values) + if debug: write_debug_output( debug_path / f"{request.node.callspec.id}.tsv", @@ -258,9 +266,53 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): 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.array, + point: np.array, + 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() + with pd.option_context("display.max_columns", None, "display.width", None): + message = ( + f"Gradient check failed:\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 ): diff --git a/tests/conftest.py b/tests/conftest.py index 9b7dd7fb08..7642e73f1d 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -76,7 +76,7 @@ def pytest_generate_tests(metafunc): else: # Run all tests test_numbers = get_all_semantic_case_ids() - + test_numbers = map(format_test_id, test_numbers) metafunc.parametrize("test_number", test_numbers) @@ -89,8 +89,6 @@ def pytest_sessionfinish(session, exitstatus): terminalreporter.ensure_newline() # parse test names to get passed case IDs (don't know any better way to # access fixture values) - from testSBMLSuite import format_test_id - passed_ids = [format_test_id(_) for _ in passed_ids] if passed_ids: write_passed_tags(passed_ids, terminalreporter) @@ -157,3 +155,8 @@ def get_tags_for_test(test_id: str) -> tuple[set[str], set[str]]: return component_tags, test_tags print(f"No componentTags or testTags found for test case {test_id}.") return component_tags, test_tags + + +def format_test_id(test_id) -> str: + """Format numeric to 0-padded string""" + return f"{test_id:0>5}" diff --git a/tests/testSBMLSuite.py b/tests/testSBMLSuite.py index 2feb3ea7e8..d843efa5e8 100755 --- a/tests/testSBMLSuite.py +++ b/tests/testSBMLSuite.py @@ -24,6 +24,7 @@ from amici.constants import SymbolId from amici.gradient_check import check_derivatives from numpy.testing import assert_allclose +from conftest import format_test_id @pytest.fixture(scope="session") @@ -351,8 +352,3 @@ def read_settings_file(current_test_path: Path, test_id: str): (key, val) = line.split(":") settings[key] = val.strip() return settings - - -def format_test_id(test_id) -> str: - """Format numeric to 0-padded string""" - return f"{test_id:0>5}"