Add gradient checks for SBML test suite

.github/workflows/test_sbml_semantic_test_suite.yml

@@ -2,6 +2,7 @@ name: SBML
 on:
   push:
     branches:
+      - check_grad_sbml_test_suite
       - main
   pull_request:
     paths:
@@ -43,6 +44,9 @@ jobs:
     - name: Install apt dependencies
       uses: ./.github/actions/install-apt-dependencies
 
+    - run: |
+        source build/venv/bin/activate && pip install git+https://github.com/ICB-DCM/fiddy.git
+
     - run: AMICI_PARALLEL_COMPILE="" ./scripts/installAmiciSource.sh
     - run: AMICI_PARALLEL_COMPILE="" ./scripts/run-SBMLTestsuite.sh ${{ matrix.cases }}
 

python/sdist/pyproject.toml

@@ -65,6 +65,7 @@ test = [
     "pytest",
     "pytest-cov",
     "pytest-rerunfailures",
+    "pytest-xdist",
     "coverage",
     "shyaml",
     "antimony>=2.13",

scripts/installAmiciSource.sh

@@ -37,6 +37,7 @@ python -m pip install --upgrade pip wheel
 #  --no-build-isolation below.
 #  The latter is necessary for code coverage to work.
 python -m pip install --upgrade pip setuptools cmake_build_extension==0.6.0 numpy petab swig
+pip install --upgrade git+https://github.com/ICB-DCM/fiddy.git
 python -m pip install git+https://github.com/pysb/pysb@master # for SPM with compartments
 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

tests/sbml/conftest.py

@@ -76,6 +76,13 @@ def pytest_generate_tests(metafunc):
         # Get CLI option
         cases = metafunc.config.getoption("cases")
         if cases:
+            # iff specific case IDs are given and the SBML semantic test suite is not there, we should fail.
+            if not SBML_SEMANTIC_CASES_DIR.exists():
+                raise ValueError(
+                    "The SBML semantic cases are missing. You can install them with "
+                    "'AMICI/scripts/run-SBMLTestsuite.sh'."
+                )
+
             # Run selected tests
             last_id = int(list(get_all_semantic_case_ids())[-1])
             test_numbers = sorted(set(parse_selection(cases, last_id)))

tests/sbml/testSBMLSuite.py

@@ -16,15 +16,21 @@
 import amici
 import pandas as pd
 import pytest
-from amici.gradient_check import check_derivatives
-
+from fiddy import MethodId, get_derivative
+from fiddy.extensions.amici import (
+    reshape,
+    run_amici_simulation_to_cached_functions,
+)
+from fiddy.success import Consistency
 from utils import (
     verify_results,
     write_result_file,
     find_model_file,
     read_settings_file,
     apply_settings,
 )
+import libsbml
+import numpy as np
 
 
 @pytest.fixture(scope="session")
@@ -37,12 +43,12 @@ def test_sbml_testsuite_case(test_id, result_path, sbml_semantic_cases_dir):
 
     # test cases for which sensitivities are to be checked
     #  key: case ID; value: epsilon for finite differences
-    sensitivity_check_cases = {
-        # parameter-dependent conservation laws
-        "00783": 1.5e-2,
-        # initial events
-        "00995": 1e-3,
-    }
+    # sensitivity_check_cases = {
+    #     # parameter-dependent conservation laws
+    #     "00783": 1.5e-2,
+    #     # initial events
+    #     "00995": 1e-3,
+    # }
 
     try:
         current_test_path = sbml_semantic_cases_dir / test_id
@@ -55,13 +61,29 @@ def test_sbml_testsuite_case(test_id, result_path, sbml_semantic_cases_dir):
             inplace=True,
         )
 
+        # TODO remove after https://github.com/AMICI-dev/AMICI/pull/2101
+        #   and https://github.com/AMICI-dev/AMICI/issues/2106
+        # Don't attempt to generate sensitivity code for models with events+algebraic rules, which will fail
+        sbml_file = find_model_file(current_test_path, test_id)
+        sbml_document = libsbml.SBMLReader().readSBMLFromFile(str(sbml_file))
+        sbml_model = sbml_document.getModel()
+        has_events = sbml_model.getNumEvents() > 0
+        has_algebraic_rules = any(
+            rule.getTypeCode() == libsbml.SBML_ALGEBRAIC_RULE
+            for rule in sbml_model.getListOfRules()
+        )
+        generate_sensitivity_code = not (has_events and has_algebraic_rules)
+        # TODO https://github.com/AMICI-dev/AMICI/issues/2109
+        generate_sensitivity_code &= test_id not in {"01240"}
+        # ^^^^^^^^
+
         # setup model
         model_dir = Path(__file__).parent / "SBMLTestModels" / test_id
         model, solver, wrapper = compile_model(
             current_test_path,
             test_id,
             model_dir,
-            generate_sensitivity_code=test_id in sensitivity_check_cases,
+            generate_sensitivity_code=generate_sensitivity_code,
         )
         settings = read_settings_file(current_test_path, test_id)
 
@@ -75,19 +97,78 @@ def test_sbml_testsuite_case(test_id, result_path, sbml_semantic_cases_dir):
             else:
                 raise RuntimeError("Simulation failed unexpectedly")
 
-        # verify
+        # verify simulation results
         simulated = verify_results(
             settings, rdata, results, wrapper, model, atol, rtol
         )
 
         # record results
         write_result_file(simulated, test_id, result_path)
 
-        # check sensitivities for selected models
-        if epsilon := sensitivity_check_cases.get(test_id):
-            solver.setSensitivityOrder(amici.SensitivityOrder.first)
-            solver.setSensitivityMethod(amici.SensitivityMethod.forward)
-            check_derivatives(model, solver, epsilon=epsilon)
+        # test sensitivities
+        if not model.getParameters():
+            pytest.skip("No parameters -> no sensitivities to check")
+
+        # TODO see https://github.com/AMICI-dev/AMICI/pull/2101
+        if not generate_sensitivity_code:
+            pytest.skip("Sensitivity analysis is known to fail.")
+        if any(id_ == 0 for id_ in model.idlist):
+            pytest.skip("Sensitivity analysis for DAE is known to fail.")
+
+        solver.setSensitivityOrder(amici.SensitivityOrder.first)
+        solver.setSensitivityMethod(amici.SensitivityMethod.forward)
+        # currently only checking "x"/"sx" for FSA
+        (
+            amici_function_f,
+            amici_derivative_f,
+            structures_f,
+        ) = run_amici_simulation_to_cached_functions(
+            amici_model=model,
+            amici_solver=solver,
+            derivative_variables=["x"],
+            cache=False,
+        )
+        rdata_f = amici.runAmiciSimulation(model, solver)
+
+        # solver.setSensitivityMethod(amici.SensitivityMethod.adjoint)
+        # (
+        #    amici_function_a,
+        #    amici_derivative_a,
+        #    structures_a,
+        # ) = run_amici_simulation_to_cached_functions(
+        #    amici_model=model,
+        #    amici_solver=solver,
+        #    derivative_variables=["x"],
+        #    cache=False,
+        # )
+        # rdata_a = amici.runAmiciSimulation(model, solver)
+
+        point = np.asarray(model.getParameters())
+
+        derivative = get_derivative(
+            # can use `_f` or `_a` here, should be no difference
+            function=amici_function_f,
+            point=point,
+            sizes=[1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1],
+            direction_ids=model.getParameterIds(),
+            method_ids=[MethodId.FORWARD, MethodId.BACKWARD, MethodId.CENTRAL],
+            relative_sizes=True,
+            success_checker=Consistency(rtol=1e-2, atol=1e-4),
+        )
+
+        derivative_fd = reshape(
+            derivative.value.flat,
+            structures_f["derivative"],
+            sensitivities=True,
+        )["x"]
+        derivative_fsa = rdata_f.sx
+        # derivative_asa = rdata_a.sllh  # currently None, define some objective?
+
+        # could alternatively use a `fiddy.DerivativeCheck` class
+        if not np.isclose(
+            derivative_fd, derivative_fsa, rtol=5e-2, atol=5e-2
+        ).all():
+            raise ValueError("Gradients were not validated.")
 
     except amici.sbml_import.SBMLException as err:
         pytest.skip(str(err))

tests/sbml/utils.py

@@ -4,6 +4,7 @@
 import numpy as np
 import pandas as pd
 from amici.constants import SymbolId
+
 from numpy.testing import assert_allclose