From 8559b8e1f0f32099a6f02dc8c27f4ec6c2eee051 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 9 Oct 2024 16:06:08 +0200 Subject: [PATCH] simplify, stricter tolerances, dbg output --- python/sdist/amici/de_model.py | 41 +++++++++++++--------------------- python/tests/util.py | 37 +++++++++++++++++++++--------- 2 files changed, 42 insertions(+), 36 deletions(-) diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py index 8eecd0d2ba..2527b0916d 100644 --- a/python/sdist/amici/de_model.py +++ b/python/sdist/amici/de_model.py @@ -1719,13 +1719,13 @@ def _compute_equation(self, name: str) -> None: elif name == "deltaxB": event_eqs = [] for ie, event in enumerate(self._events): + # ==== 1st group of terms: Heaviside functions =========== + tmp_eq = smart_multiply( + self.sym("xdot") - self.sym("xdot_old"), + self.eq("dtaudx")[ie], + ) if event._state_update is not None: - # ==== 1st group of terms : Heaviside functions =========== - tmp_eq = smart_multiply( - self.sym("xdot") - self.sym("xdot_old"), - self.eq("dtaudx")[ie], - ) - # ==== 2nd group of terms : Derivatives of Dirac deltas === + # ==== 2nd group of terms: Derivatives of Dirac deltas === # Part 2a: explicit time dependence of bolus function tmp_eq -= smart_multiply( self.eq("ddeltaxdt")[ie], self.eq("dtaudx")[ie] @@ -1737,28 +1737,22 @@ def _compute_equation(self, name: str) -> None: ), self.eq("dtaudx")[ie], ) - # ==== 3rd group of terms : Dirac deltas ================== + # ==== 3rd group of terms: Dirac deltas ================== tmp_eq += self.eq("ddeltaxdx")[ie] - tmp_eq = smart_multiply(self.sym("xB").T, tmp_eq) - else: - tmp_eq = smart_multiply( - self.sym("xdot") - self.sym("xdot_old"), - self.eq("dtaudx")[ie], - ) - tmp_eq = smart_multiply(self.sym("xB").T, tmp_eq) + tmp_eq = smart_multiply(self.sym("xB").T, tmp_eq) event_eqs.append(tmp_eq) self._eqs[name] = event_eqs elif name == "deltaqB": event_eqs = [] for ie, event in enumerate(self._events): + # ==== 1st group of terms: Heaviside functions =========== + tmp_eq = smart_multiply( + self.sym("xdot") - self.sym("xdot_old"), + self.eq("dtaudp")[ie], + ) if event._state_update is not None: - # ==== 1st group of terms : Heaviside functions =========== - tmp_eq = smart_multiply( - self.sym("xdot") - self.sym("xdot_old"), - self.eq("dtaudp")[ie], - ) - # ==== 2nd group of terms : Derivatives of Dirac deltas === + # ==== 2nd group of terms: Derivatives of Dirac deltas === # Part 2a: explicit time dependence of bolus function tmp_eq -= smart_multiply( self.eq("ddeltaxdt")[ie], self.eq("dtaudp")[ie] @@ -1770,13 +1764,8 @@ def _compute_equation(self, name: str) -> None: ), self.eq("dtaudp")[ie], ) - # ==== 3rd group of terms : Dirac deltas ================== + # ==== 3rd group of terms: Dirac deltas ================== tmp_eq += self.eq("ddeltaxdp")[ie] - else: - tmp_eq = smart_multiply( - self.sym("xdot") - self.sym("xdot_old"), - self.eq("dtaudp")[ie], - ) event_eqs.append(smart_multiply(self.sym("xB").T, tmp_eq)) self._eqs[name] = event_eqs diff --git a/python/tests/util.py b/python/tests/util.py index 1a96fb375f..1cc542e1e1 100644 --- a/python/tests/util.py +++ b/python/tests/util.py @@ -6,6 +6,8 @@ import libsbml import numpy as np +import pandas as pd + from amici import ( AmiciModel, ExpData, @@ -185,30 +187,44 @@ def check_trajectories_with_adjoint_sensitivities(amici_model: AmiciModel): solver = amici_model.getSolver() solver.setSensitivityOrder(SensitivityOrder.first) solver.setSensitivityMethod(SensitivityMethod.forward) - solver.setAbsoluteTolerance(1e-15) - solver.setRelativeTolerance(1e-10) + solver.setAbsoluteTolerance(1e-16) + solver.setRelativeTolerance(1e-14) rdata_fsa = runAmiciSimulation(amici_model, solver=solver, edata=edata) # Show that we can do arbitrary precision here (test 8 digits) solver = amici_model.getSolver() solver.setSensitivityOrder(SensitivityOrder.first) solver.setSensitivityMethod(SensitivityMethod.adjoint) - solver.setAbsoluteTolerance(1e-15) - solver.setRelativeTolerance(1e-10) - solver.setAbsoluteToleranceB(1e-15) - solver.setRelativeToleranceB(1e-10) - solver.setAbsoluteToleranceQuadratures(1e-15) - solver.setRelativeToleranceQuadratures(1e-8) + solver.setAbsoluteTolerance(1e-16) + solver.setRelativeTolerance(1e-14) + solver.setAbsoluteToleranceB(1e-16) + solver.setRelativeToleranceB(1e-15) + solver.setAbsoluteToleranceQuadratures(1e-16) + solver.setRelativeToleranceQuadratures(1e-10) rdata_asa = runAmiciSimulation(amici_model, solver=solver, edata=edata) + assert_allclose(rdata_fsa.x, rdata_asa.x, atol=1e-14, rtol=1e-10) + assert_allclose(rdata_fsa.llh, rdata_asa.llh, atol=1e-14, rtol=1e-10) + df = pd.DataFrame( + { + "fsa": rdata_fsa["sllh"], + "asa": rdata_asa["sllh"], + "fd": np.nan, + }, + index=list(amici_model.getParameterIds()), + ) + df["abs_diff"] = df["fsa"] - df["asa"] + df["rel_diff"] = df["abs_diff"] / df["fsa"] + print(df) + # Also test against finite differences parameters = amici_model.getUnscaledParameters() sllh_fd = [] eps = 1e-5 for i_par, par in enumerate(parameters): solver = amici_model.getSolver() - solver.setSensitivityOrder(SensitivityOrder.none) - solver.setSensitivityMethod(SensitivityMethod.none) + solver.setSensitivityOrder(SensitivityOrder.first) + solver.setSensitivityMethod(SensitivityMethod.adjoint) solver.setAbsoluteTolerance(1e-15) solver.setRelativeTolerance(1e-13) tmp_par = np.array(parameters[:]) @@ -220,6 +236,7 @@ def check_trajectories_with_adjoint_sensitivities(amici_model: AmiciModel): amici_model.setParameters(tmp_par) rdata_m = runAmiciSimulation(amici_model, solver=solver, edata=edata) sllh_fd.append((rdata_p["llh"] - rdata_m["llh"]) / (2 * eps)) + df["fd"] = sllh_fd # test less strict in terms of absolute error, as the gradient are # typically in the order of 1e3