simplify, stricter tolerances, dbg output

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
 

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