Fix zero sensitivity for conservation laws with parameter-dependent initial values

.github/workflows/test_sbml_semantic_test_suite.yml

@@ -46,6 +46,10 @@ jobs:
       uses: ./.github/actions/install-apt-dependencies
 
     - run: AMICI_PARALLEL_COMPILE="" ./scripts/installAmiciSource.sh
+
+    - name: Install test dependencies
+      run: python3 -m pip install --upgrade git+https://github.com/ICB-DCM/fiddy.git@main
+
     - run: AMICI_PARALLEL_COMPILE="" ./scripts/run-SBMLTestsuite.sh ${{ matrix.cases }}
     - run: mv tests/sbml/amici-semantic-results/results.json tests/sbml/amici-semantic-results/results_${{ matrix.cases }}.json
 

python/sdist/amici/_symbolic/de_model.py

@@ -1166,22 +1166,7 @@ def _generate_symbol(self, name: str) -> None:
             )
             return
         elif name == "dtcldp":
-            # check, whether the CL consists of only one state. Then,
-            # sensitivities drop out, otherwise generate symbols
-            self._syms[name] = sp.Matrix(
-                [
-                    [
-                        sp.Symbol(
-                            f"s{tcl.get_id()}__{par.get_id()}",
-                            real=True,
-                        )
-                        for par in self._free_parameters
-                    ]
-                    if self.conservation_law_has_multispecies(tcl)
-                    else [0] * self.num_par()
-                    for tcl in self._conservation_laws
-                ]
-            )
+            self._syms[name] = self._dtcldp_symbols()
             return
         elif name == "x_old":
             length = len(self.eq("xdot"))
@@ -2440,20 +2425,34 @@ def state_is_constant(self, ix: int) -> bool:
 
         return state.get_dt().is_zero
 
-    def conservation_law_has_multispecies(self, tcl: ConservationLaw) -> bool:
+    def _dtcldp_symbols(self) -> sp.Matrix:
         """
-        Checks whether a conservation law has multiple species or it just
-        defines one constant species
+        Builds the symbol matrix for ``dtcldp``, the sensitivity of each
+        conservation law's total abundance w.r.t. the free parameters.
 
-        :param tcl:
-            conservation law
+        This is always a symbol, never a literal zero: whether a
+        conservation law's total abundance actually depends on a free
+        parameter can't be decided from the model's own symbolic initial
+        conditions alone. Preequilibration (steady-state Newton solve) and
+        state reinitialization compute/override initial states and
+        sensitivities at runtime, entirely bypassing the model's compiled
+        ``fx0``/``fsx0``, so a species' initial value can depend on a
+        parameter even where the static formula shows no such dependence.
+        The correct numeric value is always computed at runtime via
+        ``fdtotal_cldp``/``fdtotal_cldx_rdata``.
 
         :return:
-            boolean indicating if conservation_law is not None
+            symbol matrix, one row per conservation law
         """
-        state_set = set(self.sym("x_rdata"))
-        n_species = len(state_set.intersection(tcl.get_val().free_symbols))
-        return n_species > 1
+        return sp.Matrix(
+            [
+                [
+                    symbol_with_assumptions(f"s{tcl.get_id()}__{par.get_id()}")
+                    for par in self._free_parameters
+                ]
+                for tcl in self._conservation_laws
+            ]
+        )
 
     def _expr_is_time_dependent(self, expr: sp.Expr) -> bool:
         """Determine whether an expression is time-dependent.

python/sdist/amici/adapters/fiddy.py

@@ -179,7 +179,7 @@ def output_labels_for_derivatives(
 def run_simulation_to_function_and_derivative(
     amici_model: AmiciModel,
     *,
-    cache: bool = True,
+    cache: bool = False,
     free_parameter_ids: list[str] = None,
     amici_solver: AmiciSolver = None,
     amici_edata: AmiciExpData = None,
@@ -293,7 +293,7 @@ def simulate_petab_to_function_and_derivative(
     *,
     amici_model: Model,
     free_parameter_ids: list[str] = None,
-    cache: bool = True,
+    cache: bool = False,
     precreate_edatas: bool = True,
     precreate_parameter_mapping: bool = True,
     simulate_petab: Callable[[Any], str] = None,
@@ -416,7 +416,7 @@ def simulate_petab_v2_to_function_and_derivative(
     petab_simulator: PetabSimulator,
     *,
     free_parameter_ids: list[str] = None,
-    cache: bool = True,
+    cache: bool = False,
 ) -> tuple[Type.FUNCTION, Type.FUNCTION]:
     r"""Create a fiddy-checkable ``(function, derivative)`` pair for a
     `PetabSimulator`, e.g. for :func:`fiddy.check_gradient`.

tests/sbml/conftest.py

@@ -16,22 +16,35 @@
 if str(script_dir) not in sys.path:
     sys.path.insert(0, str(script_dir))
 
-# the two independently-reported checks per SBML semantic test suite case
+# the independently-reported checks per SBML semantic test suite case
 SIMULATION_CHECK = "test_sbml_testsuite_case"
-SENSITIVITY_CHECK = "test_sbml_testsuite_case_sensitivity"
+SENSITIVITY_FORWARD_CHECK = "test_sbml_testsuite_case_sensitivity_forward"
+SENSITIVITY_ADJOINT_CHECK = "test_sbml_testsuite_case_sensitivity_adjoint"
+SENSITIVITY_CONSISTENCY_CHECK = (
+    "test_sbml_testsuite_case_sensitivity_consistency"
+)
+CHECKS = (
+    SIMULATION_CHECK,
+    SENSITIVITY_FORWARD_CHECK,
+    SENSITIVITY_ADJOINT_CHECK,
+    SENSITIVITY_CONSISTENCY_CHECK,
+)
+# short suffix per check, used for the `results.json` field names
+_CHECK_SUFFIXES = {
+    SIMULATION_CHECK: "simulation",
+    SENSITIVITY_FORWARD_CHECK: "sensitivity_forward",
+    SENSITIVITY_ADJOINT_CHECK: "sensitivity_adjoint",
+    SENSITIVITY_CONSISTENCY_CHECK: "sensitivity_consistency",
+}
 
 # stores passed SBML semantic test suite IDs, by check
-passed_ids: dict[str, list[str]] = {
-    SIMULATION_CHECK: [],
-    SENSITIVITY_CHECK: [],
-}
+passed_ids: dict[str, list[str]] = {check: [] for check in CHECKS}
 # test tags we encountered (from the simulation check only -- that's what
 # the SBML test suite's own tag-support semantics are about)
 encountered_tags: set[str] = set()
 # failed/skipped tests with error message, by check
 failed_or_skipped_ids: dict[str, dict[str, str]] = {
-    SIMULATION_CHECK: {},
-    SENSITIVITY_CHECK: {},
+    check: {} for check in CHECKS
 }
 
 SBML_SEMANTIC_CASES_DIR = (
@@ -171,37 +184,34 @@ def write_passed_tags(passed_simulation_ids, out=sys.stdout):
         )
         out.write("  " + "\n  ".join(sorted(passed_test_tags)))
 
-    with open(RESULT_PATH / "results.json", "w") as f:
-        json.dump(
-            {
-                "supported_tags": sorted(
-                    passed_test_tags | passed_component_tags
-                ),
-                "encountered_tags": sorted(encountered_tags),
-                "passed_tests_simulation": sorted(passed_simulation_ids),
-                "passed_tests_sensitivity": sorted(
-                    passed_ids[SENSITIVITY_CHECK]
-                ),
-                "failed_or_skipped_simulation": {
-                    k: failed_or_skipped_ids[SIMULATION_CHECK][k]
-                    for k in sorted(failed_or_skipped_ids[SIMULATION_CHECK])
-                },
-                "failed_or_skipped_sensitivity": {
-                    k: failed_or_skipped_ids[SENSITIVITY_CHECK][k]
-                    for k in sorted(failed_or_skipped_ids[SENSITIVITY_CHECK])
-                },
-            },
-            f,
-            indent=2,
+    result = {
+        "supported_tags": sorted(passed_test_tags | passed_component_tags),
+        "encountered_tags": sorted(encountered_tags),
+    }
+    for check in CHECKS:
+        suffix = _CHECK_SUFFIXES[check]
+        ids = (
+            passed_simulation_ids
+            if check == SIMULATION_CHECK
+            else passed_ids[check]
         )
+        result[f"passed_tests_{suffix}"] = sorted(ids)
+        result[f"failed_or_skipped_{suffix}"] = {
+            k: failed_or_skipped_ids[check][k]
+            for k in sorted(failed_or_skipped_ids[check])
+        }
+
+    with open(RESULT_PATH / "results.json", "w") as f:
+        json.dump(result, f, indent=2)
 
 
 def pytest_runtest_logreport(report: "TestReport") -> None:
     """Collect test case IDs of passed SBML semantic test suite cases"""
     if report.when != "call":
         return
     match = re.search(
-        r"::(test_sbml_testsuite_case(?:_sensitivity)?)\[(\d+)\]",
+        r"::(test_sbml_testsuite_case"
+        r"(?:_sensitivity_(?:forward|adjoint|consistency))?)\[(\d+)\]",
         report.nodeid,
     )
     if not match:

tests/sbml/consolidate_results.py

@@ -13,29 +13,33 @@
 # where all result artifacts have been unpacked to
 result_dir = Path("combined")
 
+CHECKS = (
+    "simulation",
+    "sensitivity_forward",
+    "sensitivity_adjoint",
+    "sensitivity_consistency",
+)
+
 # tags encountered across all tests (from the simulation check only)
 encountered_tags: set[str] = set()
 # tags for which at least one test passed
 supported_tags: set[str] = set()
 
 # test IDs of passed tests, by check
-passed_ids: dict[str, set[str]] = {"simulation": set(), "sensitivity": set()}
+passed_ids: dict[str, set[str]] = {check: set() for check in CHECKS}
 # failed or skipped tests with error message, by check
-failed_or_skipped: dict[str, dict[str, str]] = {
-    "simulation": dict(),
-    "sensitivity": dict(),
-}
+failed_or_skipped: dict[str, dict[str, str]] = {check: {} for check in CHECKS}
 
 for tag_file in result_dir.glob("results_*.json"):
     with open(tag_file) as f:
         cur_tags = json.load(f)
     encountered_tags |= set(cur_tags["encountered_tags"])
     supported_tags |= set(cur_tags["supported_tags"])
-    for check in ("simulation", "sensitivity"):
+    for check in CHECKS:
         passed_ids[check] |= set(cur_tags[f"passed_tests_{check}"])
         failed_or_skipped[check] |= cur_tags[f"failed_or_skipped_{check}"]
 
-for check in ("simulation", "sensitivity"):
+for check in CHECKS:
     num_tests_success = len(passed_ids[check])
     num_tests_total = num_tests_success + len(failed_or_skipped[check])
     frac_tests_passed = num_tests_success / num_tests_total
@@ -57,7 +61,7 @@
 print()
 print(",".join(sorted(list(unsupported_tags))))
 print()
-for check in ("simulation", "sensitivity"):
+for check in CHECKS:
     print(f"Failed or-skipped tests [{check}]")
     print("-----------------------" + "-" * len(check))
     print()