From a20447a0982f014f171985dbd3b0c58e7632e026 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Wed, 23 Oct 2024 19:08:20 +0200
Subject: [PATCH 01/16] Add missing simulation status codes (#2560)

Closes #2559.
---
 include/amici/defines.h |  6 +++
 src/amici.cpp           |  6 +++
 src/solver_cvodes.cpp   | 83 ++++++++++++-----------------------------
 3 files changed, 36 insertions(+), 59 deletions(-)

diff --git a/include/amici/defines.h b/include/amici/defines.h
index 1c6741f8d2..b49ebc6a83 100644
--- a/include/amici/defines.h
+++ b/include/amici/defines.h
@@ -74,7 +74,13 @@ constexpr int AMICI_CONSTR_FAIL=             -15;
 constexpr int AMICI_CVODES_CONSTR_FAIL=      -15;
 constexpr int AMICI_IDAS_CONSTR_FAIL=        -11;
 constexpr int AMICI_ILL_INPUT=               -22;
+constexpr int AMICI_BAD_T=                   -25;
+constexpr int AMICI_BAD_DKY=                 -26;
 constexpr int AMICI_FIRST_QRHSFUNC_ERR=      -32;
+constexpr int AMICI_SRHSFUNC_FAIL=           -41;
+constexpr int AMICI_FIRST_SRHSFUNC_ERR=      -42;
+constexpr int AMICI_REPTD_SRHSFUNC_ERR=      -43;
+constexpr int AMICI_UNREC_SRHSFUNC_ERR=      -44;
 constexpr int AMICI_ERROR=                   -99;
 constexpr int AMICI_NO_STEADY_STATE=         -81;
 constexpr int AMICI_DAMPING_FACTOR_ERROR=    -86;
diff --git a/src/amici.cpp b/src/amici.cpp
index bba344358f..8312200e0b 100644
--- a/src/amici.cpp
+++ b/src/amici.cpp
@@ -51,6 +51,12 @@ std::map<int, std::string> simulation_status_to_str_map = {
     {AMICI_LSETUP_FAIL, "AMICI_LSETUP_FAIL"},
     {AMICI_FIRST_QRHSFUNC_ERR, "AMICI_FIRST_QRHSFUNC_ERR"},
     {AMICI_WARNING, "AMICI_WARNING"},
+    {AMICI_BAD_T, "AMICI_BAD_T"},
+    {AMICI_BAD_DKY, "AMICI_BAD_DKY"},
+    {AMICI_FIRST_SRHSFUNC_ERR, "AMICI_FIRST_SRHSFUNC_ERR"},
+    {AMICI_SRHSFUNC_FAIL, "AMICI_SRHSFUNC_FAIL"},
+    {AMICI_REPTD_SRHSFUNC_ERR, "AMICI_REPTD_SRHSFUNC_ERR"},
+    {AMICI_UNREC_SRHSFUNC_ERR, "AMICI_UNREC_SRHSFUNC_ERR"},
 };
 
 std::unique_ptr<ReturnData> runAmiciSimulation(
diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp
index f777011329..dfccbb2389 100644
--- a/src/solver_cvodes.cpp
+++ b/src/solver_cvodes.cpp
@@ -34,65 +34,30 @@ static_assert((int)InterpolationType::polynomial == CV_POLYNOMIAL, "");
 static_assert((int)LinearMultistepMethod::adams == CV_ADAMS, "");
 static_assert((int)LinearMultistepMethod::BDF == CV_BDF, "");
 
-static_assert(AMICI_ROOT_RETURN == CV_ROOT_RETURN, "");
-
-static_assert(
-    amici::AMICI_SUCCESS == CV_SUCCESS, "AMICI_SUCCESS != CV_SUCCESS"
-);
-static_assert(
-    amici::AMICI_DATA_RETURN == CV_TSTOP_RETURN,
-    "AMICI_DATA_RETURN != CV_TSTOP_RETURN"
-);
-static_assert(
-    amici::AMICI_ROOT_RETURN == CV_ROOT_RETURN,
-    "AMICI_ROOT_RETURN != CV_ROOT_RETURN"
-);
-static_assert(
-    amici::AMICI_ILL_INPUT == CV_ILL_INPUT, "AMICI_ILL_INPUT != CV_ILL_INPUT"
-);
-static_assert(amici::AMICI_NORMAL == CV_NORMAL, "AMICI_NORMAL != CV_NORMAL");
-static_assert(
-    amici::AMICI_ONE_STEP == CV_ONE_STEP, "AMICI_ONE_STEP != CV_ONE_STEP"
-);
-static_assert(
-    amici::AMICI_TOO_MUCH_ACC == CV_TOO_MUCH_ACC,
-    "AMICI_TOO_MUCH_ACC != CV_TOO_MUCH_ACC"
-);
-static_assert(
-    amici::AMICI_TOO_MUCH_WORK == CV_TOO_MUCH_WORK,
-    "AMICI_TOO_MUCH_WORK != CV_TOO_MUCH_WORK"
-);
-static_assert(
-    amici::AMICI_ERR_FAILURE == CV_ERR_FAILURE,
-    "AMICI_ERR_FAILURE != CV_ERR_FAILURE"
-);
-static_assert(
-    amici::AMICI_CONV_FAILURE == CV_CONV_FAILURE,
-    "AMICI_CONV_FAILURE != CV_CONV_FAILURE"
-);
-static_assert(
-    amici::AMICI_LSETUP_FAIL == CV_LSETUP_FAIL,
-    "AMICI_LSETUP_FAIL != CV_LSETUP_FAIL"
-);
-static_assert(
-    amici::AMICI_CONSTR_FAIL == CV_CONSTR_FAIL,
-    "AMICI_CONSTR_FAIL != CV_CONSTR_FAIL"
-);
-static_assert(
-    amici::AMICI_RHSFUNC_FAIL == CV_RHSFUNC_FAIL,
-    "AMICI_RHSFUNC_FAIL != CV_RHSFUNC_FAIL"
-);
-static_assert(
-    amici::AMICI_FIRST_RHSFUNC_ERR == CV_FIRST_RHSFUNC_ERR,
-    "AMICI_FIRST_RHSFUNC_ERR != CV_FIRST_RHSFUNC_ERR"
-);
-static_assert(
-    amici::AMICI_FIRST_QRHSFUNC_ERR == CV_FIRST_QRHSFUNC_ERR,
-    "AMICI_FIRST_QRHSFUNC_ERR != CV_FIRST_QRHSFUNC_ERR"
-);
-static_assert(
-    amici::AMICI_WARNING == CV_WARNING, "AMICI_WARNING != CV_WARNING"
-);
+#define STATIC_ASSERT_EQUAL(amici_constant, cv_constant) \
+static_assert(amici_constant == cv_constant, #amici_constant " != " #cv_constant)
+
+STATIC_ASSERT_EQUAL(amici::AMICI_SUCCESS, CV_SUCCESS);
+STATIC_ASSERT_EQUAL(amici::AMICI_ROOT_RETURN, CV_ROOT_RETURN);
+STATIC_ASSERT_EQUAL(amici::AMICI_DATA_RETURN, CV_TSTOP_RETURN);
+STATIC_ASSERT_EQUAL(amici::AMICI_ILL_INPUT, CV_ILL_INPUT);
+STATIC_ASSERT_EQUAL(amici::AMICI_NORMAL, CV_NORMAL);
+STATIC_ASSERT_EQUAL(amici::AMICI_ONE_STEP, CV_ONE_STEP);
+STATIC_ASSERT_EQUAL(amici::AMICI_TOO_MUCH_ACC, CV_TOO_MUCH_ACC);
+STATIC_ASSERT_EQUAL(amici::AMICI_TOO_MUCH_WORK, CV_TOO_MUCH_WORK);
+STATIC_ASSERT_EQUAL(amici::AMICI_ERR_FAILURE, CV_ERR_FAILURE);
+STATIC_ASSERT_EQUAL(amici::AMICI_CONV_FAILURE, CV_CONV_FAILURE);
+STATIC_ASSERT_EQUAL(amici::AMICI_LSETUP_FAIL, CV_LSETUP_FAIL);
+STATIC_ASSERT_EQUAL(amici::AMICI_CONSTR_FAIL, CV_CONSTR_FAIL);
+STATIC_ASSERT_EQUAL(amici::AMICI_RHSFUNC_FAIL, CV_RHSFUNC_FAIL);
+STATIC_ASSERT_EQUAL(amici::AMICI_FIRST_RHSFUNC_ERR, CV_FIRST_RHSFUNC_ERR);
+STATIC_ASSERT_EQUAL(amici::AMICI_FIRST_QRHSFUNC_ERR, CV_FIRST_QRHSFUNC_ERR);
+STATIC_ASSERT_EQUAL(amici::AMICI_BAD_T, CV_BAD_T);
+STATIC_ASSERT_EQUAL(amici::AMICI_BAD_DKY, CV_BAD_DKY);
+STATIC_ASSERT_EQUAL(amici::AMICI_SRHSFUNC_FAIL, CV_SRHSFUNC_FAIL);
+STATIC_ASSERT_EQUAL(amici::AMICI_FIRST_SRHSFUNC_ERR, CV_FIRST_SRHSFUNC_ERR);
+STATIC_ASSERT_EQUAL(amici::AMICI_REPTD_SRHSFUNC_ERR, CV_REPTD_SRHSFUNC_ERR);
+STATIC_ASSERT_EQUAL(amici::AMICI_UNREC_SRHSFUNC_ERR, CV_UNREC_SRHSFUNC_ERR);
 
 /*
  * The following static members are callback function to CVODES.

From 9c603e1015ff0f233078f4e741b564bb9e9dc6c8 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Thu, 24 Oct 2024 10:53:08 +0200
Subject: [PATCH 02/16] PEtab import: Fix potentially incorrect sensitivities
 with state-dependent sigmas (#2562)

Due to PEtab not allowing and observables in noiseFormula and
AMICI not allowing state variables in sigma expressions, we are trying
to replace any state variables by matching observables during PEtab import.

However, if there were multiple observables with the same observableFormula,
or one observableFormula was contained in another one, the substition
of observableFormula for observableId could - depending on the ordering -
result in noiseFormulas depending on other observables, which could lead to
incorrect sensitivities due to https://github.com/AMICI-dev/AMICI/issues/2561.

This change ensures that we are only using the observableId from the same
row of the observable table as the noiseFormula.

---------

Co-authored-by: Dilan Pathirana <59329744+dilpath@users.noreply.github.com>
---
 python/sdist/amici/petab/import_helpers.py | 26 ++++++++++------------
 1 file changed, 12 insertions(+), 14 deletions(-)

diff --git a/python/sdist/amici/petab/import_helpers.py b/python/sdist/amici/petab/import_helpers.py
index b52d74004d..70af87c3b3 100644
--- a/python/sdist/amici/petab/import_helpers.py
+++ b/python/sdist/amici/petab/import_helpers.py
@@ -47,8 +47,8 @@ def get_observation_model(
 
     observables = {}
     sigmas = {}
-
     nan_pat = r"^[nN]a[nN]$"
+
     for _, observable in observable_df.iterrows():
         oid = str(observable.name)
         # need to sanitize due to https://github.com/PEtab-dev/PEtab/issues/447
@@ -56,20 +56,18 @@ def get_observation_model(
         formula_obs = re.sub(nan_pat, "", str(observable[OBSERVABLE_FORMULA]))
         formula_noise = re.sub(nan_pat, "", str(observable[NOISE_FORMULA]))
         observables[oid] = {"name": name, "formula": formula_obs}
-        sigmas[oid] = formula_noise
-
-    # PEtab does currently not allow observables in noiseFormula and AMICI
-    #  cannot handle states in sigma expressions. Therefore, where possible,
-    #  replace species occurring in error model definition by observableIds.
-    replacements = {
-        sp.sympify(observable["formula"], locals=_clash): sp.Symbol(
-            observable_id
+        formula_noise_sym = sp.sympify(formula_noise, locals=_clash)
+        formula_obs_sym = sp.sympify(formula_obs, locals=_clash)
+        # PEtab does currently not allow observables in noiseFormula, and
+        #  AMICI cannot handle state variables in sigma expressions.
+        #  Therefore, where possible, replace state variables occurring in
+        #  the error model definition by equivalent observableIds.
+        #  Do not use observableIds from other rows
+        #  (https://github.com/AMICI-dev/AMICI/issues/2561).
+        formula_noise_sym = formula_noise_sym.subs(
+            formula_obs_sym, sp.Symbol(oid)
         )
-        for observable_id, observable in observables.items()
-    }
-    for observable_id, formula in sigmas.items():
-        repl = sp.sympify(formula, locals=_clash).subs(replacements)
-        sigmas[observable_id] = str(repl)
+        sigmas[oid] = str(formula_noise_sym)
 
     noise_distrs = petab_noise_distributions_to_amici(observable_df)
 

From 56e69568ad4c5de3f074b0d9143eeadf5e82367b Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Thu, 24 Oct 2024 15:05:55 +0200
Subject: [PATCH 03/16] Fix benchmark collection gradient check (#2564)

* Fix instance vs class attribute for step sizes
* Only flatten problems where necessary
---
 .../test_benchmark_collection.sh              | 15 +++++++++-
 .../benchmark-models/test_petab_benchmark.py  | 29 ++++++++++++-------
 tests/benchmark-models/test_petab_model.py    |  7 ++++-
 3 files changed, 38 insertions(+), 13 deletions(-)

diff --git a/tests/benchmark-models/test_benchmark_collection.sh b/tests/benchmark-models/test_benchmark_collection.sh
index 581b8db028..d3d1cad712 100755
--- a/tests/benchmark-models/test_benchmark_collection.sh
+++ b/tests/benchmark-models/test_benchmark_collection.sh
@@ -93,9 +93,22 @@ for model in $models; do
     yaml="${model_dir}"/"${model}"/problem.yaml
   fi
 
+  # problems we need to flatten
+  to_flatten=(
+    "Bruno_JExpBot2016" "Chen_MSB2009" "Crauste_CellSystems2017"
+    "Fiedler_BMCSystBiol2016" "Fujita_SciSignal2010" "SalazarCavazos_MBoC2020"
+  )
+  flatten=""
+  for item in "${to_flatten[@]}"; do
+    if [[ "$item" == "$model" ]]; then
+      flatten="--flatten"
+      break
+    fi
+  done
+
   amici_model_dir=test_bmc/"${model}"
   mkdir -p "$amici_model_dir"
-  cmd_import="amici_import_petab ${yaml} -o ${amici_model_dir} -n ${model} --flatten"
+  cmd_import="amici_import_petab ${yaml} -o ${amici_model_dir} -n ${model} ${flatten}"
   cmd_run="$script_path/test_petab_model.py -y ${yaml} -d ${amici_model_dir} -m ${model} -c"
 
   printf '=%.0s' {1..40}
diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py
index 4892100877..69df16f181 100644
--- a/tests/benchmark-models/test_petab_benchmark.py
+++ b/tests/benchmark-models/test_petab_benchmark.py
@@ -12,8 +12,9 @@
 from amici.petab.petab_import import import_petab_problem
 import benchmark_models_petab
 from collections import defaultdict
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 from amici import SensitivityMethod
+from petab.v1.lint import measurement_table_has_timepoint_specific_mappings
 from fiddy import MethodId, get_derivative
 from fiddy.derivative_check import NumpyIsCloseDerivativeCheck
 from fiddy.extensions.amici import simulate_petab_to_cached_functions
@@ -58,14 +59,18 @@ class GradientCheckSettings:
     atol_consistency: float = 1e-5
     rtol_consistency: float = 1e-1
     # Step sizes for finite difference gradient checks.
-    step_sizes = [
-        1e-1,
-        5e-2,
-        1e-2,
-        1e-3,
-        1e-4,
-        1e-5,
-    ]
+    step_sizes: list[float] = field(
+        default_factory=lambda: [
+            2e-1,
+            1e-1,
+            5e-2,
+            1e-2,
+            5e-1,
+            1e-3,
+            1e-4,
+            1e-5,
+        ]
+    )
     rng_seed: int = 0
     ss_sensitivity_mode: amici.SteadyStateSensitivityMode = (
         amici.SteadyStateSensitivityMode.integrateIfNewtonFails
@@ -97,7 +102,6 @@ class GradientCheckSettings:
     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,
@@ -176,7 +180,10 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request):
         pytest.skip()
 
     petab_problem = benchmark_models_petab.get_problem(model)
-    petab.flatten_timepoint_specific_output_overrides(petab_problem)
+    if measurement_table_has_timepoint_specific_mappings(
+        petab_problem.measurement_df,
+    ):
+        petab.flatten_timepoint_specific_output_overrides(petab_problem)
 
     # Only compute gradient for estimated parameters.
     parameter_ids = petab_problem.x_free_ids
diff --git a/tests/benchmark-models/test_petab_model.py b/tests/benchmark-models/test_petab_model.py
index c4ec2f5dd2..125a046a5e 100755
--- a/tests/benchmark-models/test_petab_model.py
+++ b/tests/benchmark-models/test_petab_model.py
@@ -25,6 +25,7 @@
     simulate_petab,
 )
 from petab.v1.visualize import plot_problem
+from petab.v1.lint import measurement_table_has_timepoint_specific_mappings
 
 logger = get_logger(f"amici.{__name__}", logging.WARNING)
 
@@ -115,7 +116,11 @@ def main():
 
     # load PEtab files
     problem = petab.Problem.from_yaml(args.yaml_file_name)
-    petab.flatten_timepoint_specific_output_overrides(problem)
+
+    if measurement_table_has_timepoint_specific_mappings(
+        problem.measurement_df
+    ):
+        petab.flatten_timepoint_specific_output_overrides(problem)
 
     # load model
     if args.model_directory:

From 3b7c6d1de820fa1b7d536d3a45e2014333541e25 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Thu, 24 Oct 2024 18:55:40 +0200
Subject: [PATCH 04/16] GHA: PYTHONFAULTHANDLER=1 (#2566)

Set PYTHONFAULTHANDLER=1 for more output on crashes (#2565). So far, this was only set for Ubuntu tests.
---
 .github/workflows/test_python_cplusplus.yml | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml
index 23337986db..284a170031 100644
--- a/.github/workflows/test_python_cplusplus.yml
+++ b/.github/workflows/test_python_cplusplus.yml
@@ -247,6 +247,8 @@ jobs:
   macos_cpp_py:
     name: Tests MacOS C++/Python
     runs-on: macos-latest
+    env:
+      PYTHONFAULTHANDLER: "1"
 
     steps:
     - name: Set up Python
@@ -300,6 +302,8 @@ jobs:
   macos_python:
     name: Tests MacOS Python
     runs-on: macos-latest
+    env:
+      PYTHONFAULTHANDLER: "1"
 
     steps:
     - name: Cache

From f5fa6cdfd5813bdaf2e4799df10eee603bec9cd3 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Thu, 24 Oct 2024 20:28:51 +0200
Subject: [PATCH 05/16] Refactor benchmark collection tests (#2569)

* Refactor tests/benchmark-models/test_benchmark_collection.sh into pytest cases
* Unify with gradient checks in `tests/benchmark-models/test_petab_benchmark.py`
* Test some additional benchmark problems

For now, just consolidating the old tests. To be cleaned up further later on.

Closes #2510.
---
 .../test_benchmark_collection_models.yml      |  27 +-
 .../test_benchmark_collection.sh              | 149 ---------
 .../benchmark-models/test_petab_benchmark.py  | 304 ++++++++++++++++--
 tests/benchmark-models/test_petab_model.py    | 268 ---------------
 4 files changed, 299 insertions(+), 449 deletions(-)
 delete mode 100755 tests/benchmark-models/test_benchmark_collection.sh
 delete mode 100755 tests/benchmark-models/test_petab_model.py

diff --git a/.github/workflows/test_benchmark_collection_models.yml b/.github/workflows/test_benchmark_collection_models.yml
index 0b14d73716..3e46af8992 100644
--- a/.github/workflows/test_benchmark_collection_models.yml
+++ b/.github/workflows/test_benchmark_collection_models.yml
@@ -52,25 +52,28 @@ jobs:
         AMICI_PARALLEL_COMPILE="" pip3 install -v --user \
             $(ls -t python/sdist/dist/amici-*.tar.gz | head -1)[petab,test,vis]
 
-    - run: |
+    - name: Install test dependencies
+      run: |
         python3 -m pip uninstall -y petab && python3 -m pip install git+https://github.com/petab-dev/libpetab-python.git@develop \
-        &&  python3 -m pip install -U sympy
+        &&  python3 -m pip install -U sympy \
+        &&  python3 -m pip install git+https://github.com/ICB-DCM/fiddy.git
 
-    # retrieve test models
-    - name: Download and test benchmark collection
+    - name: Download benchmark collection
       run: |
         git clone --depth 1 https://github.com/benchmarking-initiative/Benchmark-Models-PEtab.git \
-          && export BENCHMARK_COLLECTION="$(pwd)/Benchmark-Models-PEtab/Benchmark-Models/" \
-          && pip3 install -e $BENCHMARK_COLLECTION/../src/python \
-          && AMICI_PARALLEL_COMPILE="" tests/benchmark-models/test_benchmark_collection.sh
+          && python3 -m pip install -e Benchmark-Models-PEtab/src/python
+
+    - name: Run tests
+      env:
+        AMICI_PARALLEL_COMPILE: ""
+      run: |
+          cd tests/benchmark-models && pytest --durations=10
 
-    # run gradient checks
-    - name: Run Gradient Checks
+    # collect & upload results
+    - name: Aggregate results
       run: |
-        pip install git+https://github.com/ICB-DCM/fiddy.git \
-          && cd tests/benchmark-models && pytest --durations=10 ./test_petab_benchmark.py
+          cd tests/benchmark-models && python3 evaluate_benchmark.py
 
-    # upload results
     - uses: actions/upload-artifact@v4
       with:
         name: computation-times-${{ matrix.python-version }}-${{ matrix.extract_subexpressions }}
diff --git a/tests/benchmark-models/test_benchmark_collection.sh b/tests/benchmark-models/test_benchmark_collection.sh
deleted file mode 100755
index d3d1cad712..0000000000
--- a/tests/benchmark-models/test_benchmark_collection.sh
+++ /dev/null
@@ -1,149 +0,0 @@
-#!/bin/bash
-# Import and run selected benchmark models with nominal parameters and check
-# agreement with reference values
-#
-# Expects environment variable BENCHMARK_COLLECTION to provide path to
-# benchmark collection model directory
-
-# Confirmed to be working
-models="
-Bachmann_MSB2011
-Beer_MolBioSystems2014
-Boehm_JProteomeRes2014
-Borghans_BiophysChem1997
-Brannmark_JBC2010
-Bruno_JExpBot2016
-Crauste_CellSystems2017
-Elowitz_Nature2000
-Fiedler_BMCSystBiol2016
-Fujita_SciSignal2010
-Isensee_JCB2018
-Lucarelli_CellSystems2018
-Schwen_PONE2014
-Smith_BMCSystBiol2013
-Sneyd_PNAS2002
-Weber_BMC2015
-Zheng_PNAS2012"
-
-#
-# not merged:
-# Becker_Science2010 (multiple models)             https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Becker_Science2010
-# Hass_PONE2017 (???)                              https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Hass_PONE2017
-# Korkut_eLIFE2015 (???)                           https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Korkut_eLIFE2015
-# Casaletto_PNAS2019 (yaml missing)                https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Casaletto_PNAS2019
-# Merkle_PCB2016 (model missing)                   https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Merkle_PCB2016
-# Parmar_PCB2019 (SBML extensions)                 https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Parmar_PCB2019
-# Swameye_PNAS2003 (splines)                       https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Swameye_PNAS2003
-# Sobotta_Frontiers2017 (???)                      https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Sobotta_Frontiers2017
-# Raia_CancerResearch2011 (state dependent sigmas) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Raia_CancerResearch2011
-#
-# no reference value:
-# Alkan_SciSignal2018 (d2d: Alkan_DDRP_SciSignal2018)
-# Bertozzi_PNAS2020 (gh: vanako, doi: https://doi.org/10.1073/pnas.2006520117, code/data: https://github.com/gomohler/pnas2020 (R))
-# Blasi_CellSystems2016 (gh: Leonard Schmiester, doi: https://doi.org/10.1016/j.cels.2016.01.002, code/data: not available)
-# Giordano_Nature2020 (gh: Paul Jonas Jost, doi: https://doi.org/10.1038/s41591-020-0883-7, code/data: http://users.dimi.uniud.it/~giulia.giordano/docs/papers/SIDARTHEcode.zip (MATLAB))
-# Laske_PLOSComputBiol2019 (gh: Clemens Peiter, doi: https://doi.org/10.1128/JVI.00080-12 (?), code/data: ???)
-# Okuonghae_ChaosSolitonsFractals2020 (gh: Paul Jonas Jost, doi: https://doi.org/10.1016/j.chaos.2020.110032, code/data: ???)
-# Oliveira_NatCommun2021 (gh: lorenapohl, doi: https://doi.org/10.1038/s41467-020-19798-3, code: https://github.com/cidacslab/Mathematical-and-Statistical-Modeling-of-COVID19-in-Brazil (python) data: https://infovis.sei.ba.gov.br/covid19/ )
-# Perelson_Science1996 (gh: Philipp Staedter, doi: https://doi.org/10.1126/science.271.5255.1582, code/data: ???)
-# Rahman_MBS2016 (gh: Yannik Schaelte, doi: https://doi.org/10.1016/j.mbs.2016.07.009, code: not available, data: table in paper ...)
-# Raimundez_PCB2020 (gh: Elba Raimundez, doi: https://doi.org/10.1371/journal.pcbi.1007147, code/data: https://zenodo.org/record/2908234#.Y5hUUS8w3yw (d2d))
-# SalazarCavazos_MBoC2020 (gh: Dilan Pathirana, doi: https://doi.org/10.1091/mbc.E19-09-0548, code/data: supplement (BNGL))
-# Zhao_QuantBiol2020 (gh: Iva Ewert, doi: https://doi.org/10.1007/s40484-020-0199-0, code: not available, data: table in supp)
-#
-# covered by performance test:
-# Froehlich_CellSystems2018
-#
-# Unknown reasons:
-# Chen_MSB2009
-#
-
-set -e
-
-[[ -n "${BENCHMARK_COLLECTION}" ]] && model_dir="${BENCHMARK_COLLECTION}"
-
-function show_help() {
-  echo "-h: this help; -n: dry run, print commands; -b path_to_models_dir"
-}
-
-OPTIND=1
-while getopts "h?nb:" opt; do
-  case "$opt" in
-  h | \?)
-    show_help
-    exit 0
-    ;;
-  n)
-    dry_run=1
-    ;;
-  b)
-    model_dir=$OPTARG
-    ;;
-  esac
-done
-
-script_path=$(dirname "$BASH_SOURCE")
-script_path=$(cd "$script_path" && pwd)
-
-for model in $models; do
-  yaml="${model_dir}"/"${model}"/"${model}".yaml
-
-  # different naming scheme
-  if [[ "$model" == "Bertozzi_PNAS2020" ]]; then
-    yaml="${model_dir}"/"${model}"/problem.yaml
-  fi
-
-  # problems we need to flatten
-  to_flatten=(
-    "Bruno_JExpBot2016" "Chen_MSB2009" "Crauste_CellSystems2017"
-    "Fiedler_BMCSystBiol2016" "Fujita_SciSignal2010" "SalazarCavazos_MBoC2020"
-  )
-  flatten=""
-  for item in "${to_flatten[@]}"; do
-    if [[ "$item" == "$model" ]]; then
-      flatten="--flatten"
-      break
-    fi
-  done
-
-  amici_model_dir=test_bmc/"${model}"
-  mkdir -p "$amici_model_dir"
-  cmd_import="amici_import_petab ${yaml} -o ${amici_model_dir} -n ${model} ${flatten}"
-  cmd_run="$script_path/test_petab_model.py -y ${yaml} -d ${amici_model_dir} -m ${model} -c"
-
-  printf '=%.0s' {1..40}
-  printf "   %s   " "${model}"
-  printf '=%.0s' {1..40}
-  echo
-
-  if [[ -z "$dry_run" ]]; then
-    $cmd_import
-    $cmd_run
-  else
-    echo "$cmd_import"
-    echo "$cmd_run"
-  fi
-
-  printf '=%.0s' {1..100}
-  echo
-  echo
-done
-
-cd "$script_path" && python evaluate_benchmark.py
-
-# Test deprecated import from individual PEtab files
-model="Zheng_PNAS2012"
-problem_dir="${model_dir}/${model}"
-amici_model_dir=test_bmc/"${model}-deprecated"
-cmd_import="amici_import_petab -s "${problem_dir}/model_${model}.xml" \
-  -m "${problem_dir}/measurementData_${model}.tsv" \
-  -c "${problem_dir}/experimentalCondition_${model}.tsv" \
-  -p "${problem_dir}/parameters_${model}.tsv" \
-  -b "${problem_dir}/observables_${model}.tsv" \
-  -o ${amici_model_dir} -n ${model} --no-compile"
-
-if [[ -z "$dry_run" ]]; then
-  $cmd_import
-else
-  echo "$cmd_import"
-fi
diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py
index 69df16f181..b4e1f50e68 100644
--- a/tests/benchmark-models/test_petab_benchmark.py
+++ b/tests/benchmark-models/test_petab_benchmark.py
@@ -1,9 +1,12 @@
-"""Tests for simulate_petab on PEtab benchmark problems."""
+"""Tests based on the PEtab benchmark problems.
 
-from pathlib import Path
+Tests simulate_petab, correctness of the log-likelihood computation at nominal
+parameters, correctness of the gradient computation, and simulation times
+for a subset of the benchmark problems.
+"""
 
+from pathlib import Path
 import fiddy
-
 import amici
 import numpy as np
 import pandas as pd
@@ -19,13 +22,33 @@
 from fiddy.derivative_check import NumpyIsCloseDerivativeCheck
 from fiddy.extensions.amici import simulate_petab_to_cached_functions
 from fiddy.success import Consistency
+import contextlib
+import logging
+import yaml
+from amici.logging import get_logger
+from amici.petab.simulations import (
+    LLH,
+    RDATAS,
+    rdatas_to_measurement_df,
+    simulate_petab,
+)
+from petab.v1.visualize import plot_problem
 
 
-repo_root = Path(__file__).parent.parent.parent
+logger = get_logger(f"amici.{__name__}", logging.WARNING)
 
-# reuse compiled models from test_benchmark_collection.sh
+script_dir = Path(__file__).parent.absolute()
+repo_root = script_dir.parent.parent
 benchmark_outdir = repo_root / "test_bmc"
-models = set(benchmark_models_petab.MODELS) - {
+
+# reference values for simulation times and log-likelihoods
+references_yaml = script_dir / "benchmark_models.yaml"
+with open(references_yaml) as f:
+    reference_values = yaml.full_load(f)
+
+# problem IDs for which to check the gradient
+# TODO: extend
+problems_for_gradient_check = set(benchmark_models_petab.MODELS) - {
     # excluded due to excessive runtime
     "Bachmann_MSB2011",
     "Chen_MSB2009",
@@ -39,15 +62,75 @@
     "Smith_BMCSystBiol2013",
     # excluded due to excessive numerical failures
     "Crauste_CellSystems2017",
-    "Fujita_SciSignal2010",
-    # FIXME: re-enable
-    "Raia_CancerResearch2011",
 }
-models = list(sorted(models))
+problems_for_gradient_check = list(sorted(problems_for_gradient_check))
+
+
+# Problems for checking the log-likelihood computation at nominal parameters
+#
+# not merged:
+# Becker_Science2010 (multiple models)             https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Becker_Science2010
+# Hass_PONE2017 (???)                              https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Hass_PONE2017
+# Korkut_eLIFE2015 (???)                           https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Korkut_eLIFE2015
+# Casaletto_PNAS2019 (yaml missing)                https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Casaletto_PNAS2019
+# Merkle_PCB2016 (model missing)                   https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Merkle_PCB2016
+# Parmar_PCB2019 (SBML extensions)                 https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Parmar_PCB2019
+# Swameye_PNAS2003 (splines)                       https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Swameye_PNAS2003
+# Sobotta_Frontiers2017 (???)                      https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Sobotta_Frontiers2017
+# Raia_CancerResearch2011 (state dependent sigmas) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Raia_CancerResearch2011
+#
+# no reference value:
+# Alkan_SciSignal2018 (d2d: Alkan_DDRP_SciSignal2018)
+# Bertozzi_PNAS2020 (gh: vanako, doi: https://doi.org/10.1073/pnas.2006520117, code/data: https://github.com/gomohler/pnas2020 (R))
+# Blasi_CellSystems2016 (gh: Leonard Schmiester, doi: https://doi.org/10.1016/j.cels.2016.01.002, code/data: not available)
+# Giordano_Nature2020 (gh: Paul Jonas Jost, doi: https://doi.org/10.1038/s41591-020-0883-7, code/data: http://users.dimi.uniud.it/~giulia.giordano/docs/papers/SIDARTHEcode.zip (MATLAB))
+# Laske_PLOSComputBiol2019 (gh: Clemens Peiter, doi: https://doi.org/10.1128/JVI.00080-12 (?), code/data: ???)
+# Okuonghae_ChaosSolitonsFractals2020 (gh: Paul Jonas Jost, doi: https://doi.org/10.1016/j.chaos.2020.110032, code/data: ???)
+# Oliveira_NatCommun2021 (gh: lorenapohl, doi: https://doi.org/10.1038/s41467-020-19798-3, code: https://github.com/cidacslab/Mathematical-and-Statistical-Modeling-of-COVID19-in-Brazil (python) data: https://infovis.sei.ba.gov.br/covid19/ )
+# Perelson_Science1996 (gh: Philipp Staedter, doi: https://doi.org/10.1126/science.271.5255.1582, code/data: ???)
+# Rahman_MBS2016 (gh: Yannik Schaelte, doi: https://doi.org/10.1016/j.mbs.2016.07.009, code: not available, data: table in paper ...)
+# Raimundez_PCB2020 (gh: Elba Raimundez, doi: https://doi.org/10.1371/journal.pcbi.1007147, code/data: https://zenodo.org/record/2908234#.Y5hUUS8w3yw (d2d))
+# SalazarCavazos_MBoC2020 (gh: Dilan Pathirana, doi: https://doi.org/10.1091/mbc.E19-09-0548, code/data: supplement (BNGL))
+# Zhao_QuantBiol2020 (gh: Iva Ewert, doi: https://doi.org/10.1007/s40484-020-0199-0, code: not available, data: table in supp)
+#
+# covered by performance test:
+# Froehlich_CellSystems2018
+#
+# Unknown reasons:
+# Chen_MSB2009
+#
+# Confirmed to be working:
+# TODO: extend
+problems_for_llh_check = [
+    "Bachmann_MSB2011",
+    "Beer_MolBioSystems2014",
+    "Boehm_JProteomeRes2014",
+    "Borghans_BiophysChem1997",
+    "Brannmark_JBC2010",
+    "Bruno_JExpBot2016",
+    "Crauste_CellSystems2017",
+    "Elowitz_Nature2000",
+    "Fiedler_BMCSystBiol2016",
+    "Fujita_SciSignal2010",
+    "Isensee_JCB2018",
+    "Lucarelli_CellSystems2018",
+    "Schwen_PONE2014",
+    "Smith_BMCSystBiol2013",
+    "Sneyd_PNAS2002",
+    "Weber_BMC2015",
+    "Zheng_PNAS2012",
+]
+
+# all PEtab problems we need to import
+problems = list(
+    sorted(set(problems_for_gradient_check + problems_for_llh_check))
+)
 
 
 @dataclass
 class GradientCheckSettings:
+    """Problem-specific settings for gradient checks."""
+
     # Absolute and relative tolerances for simulation
     atol_sim: float = 1e-16
     rtol_sim: float = 1e-12
@@ -93,6 +176,10 @@ class GradientCheckSettings:
 settings["Brannmark_JBC2010"] = GradientCheckSettings(
     ss_sensitivity_mode=amici.SteadyStateSensitivityMode.integrationOnly,
 )
+settings["Fujita_SciSignal2010"] = GradientCheckSettings(
+    atol_check=1e-7,
+    rtol_check=5e-4,
+)
 settings["Giordano_Nature2020"] = GradientCheckSettings(
     atol_check=1e-6, rtol_check=1e-3, rng_seed=1
 )
@@ -107,6 +194,10 @@ class GradientCheckSettings:
     atol_sim=1e-12,
     rtol_sim=1e-12,
 )
+settings["Raia_CancerResearch2011"] = GradientCheckSettings(
+    atol_check=1e-10,
+    rtol_check=1e-3,
+)
 settings["Smith_BMCSystBiol2013"] = GradientCheckSettings(
     atol_sim=1e-10,
     rtol_sim=1e-10,
@@ -140,6 +231,175 @@ class GradientCheckSettings:
     debug_path.mkdir(exist_ok=True, parents=True)
 
 
+@pytest.fixture(scope="session", params=problems, ids=problems)
+def benchmark_problem(request):
+    """Fixture providing model and PEtab problem for a problem from
+    the benchmark problem collection."""
+    problem_id = request.param
+    petab_problem = benchmark_models_petab.get_problem(problem_id)
+    if measurement_table_has_timepoint_specific_mappings(
+        petab_problem.measurement_df,
+    ):
+        petab.flatten_timepoint_specific_output_overrides(petab_problem)
+
+    # Setup AMICI objects.
+    amici_model = import_petab_problem(
+        petab_problem,
+        model_output_dir=benchmark_outdir / problem_id,
+    )
+    return problem_id, petab_problem, amici_model
+
+
+@pytest.mark.filterwarnings(
+    "ignore:divide by zero encountered in log",
+    # https://github.com/AMICI-dev/AMICI/issues/18
+    "ignore:Adjoint sensitivity analysis for models with discontinuous "
+    "right hand sides .*:UserWarning",
+)
+def test_nominal_parameters_llh(benchmark_problem):
+    """Test the log-likelihood computation at nominal parameters.
+
+    Also check that the simulation time is within the reference range.
+    """
+    problem_id, petab_problem, amici_model = benchmark_problem
+    if problem_id not in problems_for_llh_check:
+        pytest.skip("Excluded from log-likelihood check.")
+
+    amici_solver = amici_model.getSolver()
+    amici_solver.setAbsoluteTolerance(1e-8)
+    amici_solver.setRelativeTolerance(1e-8)
+    amici_solver.setMaxSteps(10_000)
+    if problem_id in ("Brannmark_JBC2010", "Isensee_JCB2018"):
+        amici_model.setSteadyStateSensitivityMode(
+            amici.SteadyStateSensitivityMode.integrationOnly
+        )
+
+    times = dict()
+
+    for label, sensi_mode in {
+        "t_sim": amici.SensitivityMethod.none,
+        "t_fwd": amici.SensitivityMethod.forward,
+        "t_adj": amici.SensitivityMethod.adjoint,
+    }.items():
+        amici_solver.setSensitivityMethod(sensi_mode)
+        if sensi_mode == amici.SensitivityMethod.none:
+            amici_solver.setSensitivityOrder(amici.SensitivityOrder.none)
+        else:
+            amici_solver.setSensitivityOrder(amici.SensitivityOrder.first)
+
+        res_repeats = [
+            simulate_petab(
+                petab_problem=petab_problem,
+                amici_model=amici_model,
+                solver=amici_solver,
+                log_level=logging.DEBUG,
+            )
+            for _ in range(3)  # repeat to get more stable timings
+        ]
+        res = res_repeats[0]
+
+        times[label] = np.min(
+            [
+                sum(r.cpu_time + r.cpu_timeB for r in res[RDATAS]) / 1000
+                # only forwards/backwards simulation
+                for res in res_repeats
+            ]
+        )
+
+        if sensi_mode == amici.SensitivityMethod.none:
+            rdatas = res[RDATAS]
+            llh = res[LLH]
+    # TODO: check that all llh match, check that all sllh match
+    times["np"] = sum(petab_problem.parameter_df[petab.ESTIMATE])
+
+    pd.Series(times).to_csv(script_dir / f"{problem_id}_benchmark.csv")
+    for rdata in rdatas:
+        assert (
+            rdata.status == amici.AMICI_SUCCESS
+        ), f"Simulation failed for {rdata.id}"
+
+    if debug:
+        # create simulation PEtab table
+        sim_df = rdatas_to_measurement_df(
+            rdatas=rdatas,
+            model=amici_model,
+            measurement_df=petab_problem.measurement_df,
+        )
+        sim_df.rename(
+            columns={petab.MEASUREMENT: petab.SIMULATION}, inplace=True
+        )
+        sim_df.to_csv(
+            benchmark_outdir / problem_id / f"{problem_id}_simulation.tsv",
+            index=False,
+            sep="\t",
+        )
+
+        # visualize fit, save to file
+        with contextlib.suppress(NotImplementedError):
+            axs = plot_problem(
+                petab_problem=petab_problem, simulations_df=sim_df
+            )
+            for plot_id, ax in axs.items():
+                fig_path = (
+                    benchmark_outdir
+                    / problem_id
+                    / f"{problem_id}_{plot_id}_vis.png"
+                )
+                logger.info(f"Saving figure to {fig_path}")
+                ax.get_figure().savefig(fig_path, dpi=150)
+
+    try:
+        ref_llh = reference_values[problem_id]["llh"]
+        rdiff = np.abs((llh - ref_llh) / ref_llh)
+        rtol = 1e-3
+        adiff = np.abs(llh - ref_llh)
+        atol = 1e-3
+        tolstr = (
+            f" Absolute difference is {adiff:.2e} "
+            f"(tol {atol:.2e}) and relative difference is "
+            f"{rdiff:.2e} (tol {rtol:.2e})."
+        )
+
+        if np.isclose(llh, ref_llh, rtol=rtol, atol=atol):
+            logger.info(
+                f"Computed llh {llh:.4e} matches reference {ref_llh:.4e}."
+                + tolstr
+            )
+        else:
+            pytest.fail(
+                f"Computed llh {llh:.4e} does not match reference "
+                f"{ref_llh:.4e}." + tolstr
+            )
+    except KeyError:
+        logger.error(
+            "No reference likelihood found for "
+            f"{problem_id} in {references_yaml}"
+        )
+
+    for label, key in {
+        "simulation": "t_sim",
+        "adjoint sensitivity": "t_adj",
+        "forward sensitivity": "t_fwd",
+    }.items():
+        try:
+            ref = reference_values[problem_id][key]
+            if times[key] > ref:
+                pytest.fail(
+                    f"Computation time for {label} ({times[key]:.2e}) "
+                    f"exceeds reference ({ref:.2e})."
+                )
+            else:
+                logger.info(
+                    f"Computation time for {label} ({times[key]:.2e}) "
+                    f"within reference ({ref:.2e})."
+                )
+        except KeyError:
+            logger.error(
+                f"No reference time for {label} found for "
+                f"{problem_id} in {references_yaml}"
+            )
+
+
 @pytest.mark.filterwarnings(
     "ignore:divide by zero encountered in log",
     # https://github.com/AMICI-dev/AMICI/issues/18
@@ -152,9 +412,14 @@ class GradientCheckSettings:
     (SensitivityMethod.forward, SensitivityMethod.adjoint),
     ids=["forward", "adjoint"],
 )
-@pytest.mark.parametrize("model", models)
-def test_benchmark_gradient(model, scale, sensitivity_method, request):
-    if not scale and model in (
+def test_benchmark_gradient(
+    benchmark_problem, scale, sensitivity_method, request
+):
+    problem_id, petab_problem, amici_model = benchmark_problem
+    if problem_id not in problems_for_gradient_check:
+        pytest.skip("Excluded from gradient check.")
+
+    if not scale and problem_id in (
         "Smith_BMCSystBiol2013",
         "Brannmark_JBC2010",
         "Elowitz_Nature2000",
@@ -165,10 +430,10 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request):
     ):
         # not really worth the effort trying to fix these cases if they
         # only fail on linear scale
-        pytest.skip()
+        pytest.skip("scale=False disabled for this problem")
 
     if (
-        model
+        problem_id
         in (
             # events with parameter-dependent triggers
             #  https://github.com/AMICI-dev/AMICI/issues/18
@@ -176,10 +441,9 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request):
         )
         and sensitivity_method == SensitivityMethod.adjoint
     ):
-        # FIXME
-        pytest.skip()
+        pytest.skip("Unsupported ASA+events")
 
-    petab_problem = benchmark_models_petab.get_problem(model)
+    petab_problem = benchmark_models_petab.get_problem(problem_id)
     if measurement_table_has_timepoint_specific_mappings(
         petab_problem.measurement_df,
     ):
@@ -187,12 +451,12 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request):
 
     # Only compute gradient for estimated parameters.
     parameter_ids = petab_problem.x_free_ids
-    cur_settings = settings[model]
+    cur_settings = settings[problem_id]
 
     # Setup AMICI objects.
     amici_model = import_petab_problem(
         petab_problem,
-        model_output_dir=benchmark_outdir / model,
+        model_output_dir=benchmark_outdir / problem_id,
     )
     amici_solver = amici_model.getSolver()
     amici_solver.setAbsoluteTolerance(cur_settings.atol_sim)
diff --git a/tests/benchmark-models/test_petab_model.py b/tests/benchmark-models/test_petab_model.py
deleted file mode 100755
index 125a046a5e..0000000000
--- a/tests/benchmark-models/test_petab_model.py
+++ /dev/null
@@ -1,268 +0,0 @@
-#!/usr/bin/env python3
-
-"""
-Simulate a PEtab problem and compare results to reference values
-"""
-
-import argparse
-import contextlib
-import importlib
-import logging
-import os
-import sys
-from pathlib import Path
-
-import amici
-import numpy as np
-import pandas as pd
-import petab.v1 as petab
-import yaml
-from amici.logging import get_logger
-from amici.petab.simulations import (
-    LLH,
-    RDATAS,
-    rdatas_to_measurement_df,
-    simulate_petab,
-)
-from petab.v1.visualize import plot_problem
-from petab.v1.lint import measurement_table_has_timepoint_specific_mappings
-
-logger = get_logger(f"amici.{__name__}", logging.WARNING)
-
-
-def parse_cli_args():
-    """Parse command line arguments
-
-    Returns:
-        Parsed CLI arguments from ``argparse``.
-    """
-
-    parser = argparse.ArgumentParser(
-        description="Simulate PEtab-format model using AMICI."
-    )
-
-    # General options:
-    parser.add_argument(
-        "-v",
-        "--verbose",
-        dest="verbose",
-        action="store_true",
-        help="More verbose output",
-    )
-    parser.add_argument(
-        "-c",
-        "--check",
-        dest="check",
-        action="store_true",
-        help="Compare to reference value",
-    )
-    parser.add_argument(
-        "-p",
-        "--plot",
-        dest="plot",
-        action="store_true",
-        help="Plot measurement and simulation results",
-    )
-
-    # PEtab problem
-    parser.add_argument(
-        "-y",
-        "--yaml",
-        dest="yaml_file_name",
-        required=True,
-        help="PEtab YAML problem filename",
-    )
-
-    # Corresponding AMICI model
-    parser.add_argument(
-        "-m",
-        "--model-name",
-        dest="model_name",
-        help="Name of the AMICI module of the model to " "simulate.",
-        required=True,
-    )
-    parser.add_argument(
-        "-d",
-        "--model-dir",
-        dest="model_directory",
-        help="Directory containing the AMICI module of the "
-        "model to simulate. Required if model is not "
-        "in python path.",
-    )
-
-    parser.add_argument(
-        "-o",
-        "--simulation-file",
-        dest="simulation_file",
-        help="File to write simulation result to, in PEtab"
-        "measurement table format.",
-    )
-
-    return parser.parse_args()
-
-
-def main():
-    """Simulate the model specified on the command line"""
-    script_dir = Path(__file__).parent.absolute()
-    args = parse_cli_args()
-    loglevel = logging.DEBUG if args.verbose else logging.INFO
-    logger.setLevel(loglevel)
-
-    logger.info(
-        f"Simulating '{args.model_name}' "
-        f"({args.model_directory}) using PEtab data from "
-        f"{args.yaml_file_name}"
-    )
-
-    # load PEtab files
-    problem = petab.Problem.from_yaml(args.yaml_file_name)
-
-    if measurement_table_has_timepoint_specific_mappings(
-        problem.measurement_df
-    ):
-        petab.flatten_timepoint_specific_output_overrides(problem)
-
-    # load model
-    if args.model_directory:
-        sys.path.insert(0, args.model_directory)
-    model_module = importlib.import_module(args.model_name)
-    amici_model = model_module.getModel()
-    amici_solver = amici_model.getSolver()
-
-    amici_solver.setAbsoluteTolerance(1e-8)
-    amici_solver.setRelativeTolerance(1e-8)
-    amici_solver.setMaxSteps(int(1e4))
-    if args.model_name in ("Brannmark_JBC2010", "Isensee_JCB2018"):
-        amici_model.setSteadyStateSensitivityMode(
-            amici.SteadyStateSensitivityMode.integrationOnly
-        )
-
-    times = dict()
-
-    for label, sensi_mode in {
-        "t_sim": amici.SensitivityMethod.none,
-        "t_fwd": amici.SensitivityMethod.forward,
-        "t_adj": amici.SensitivityMethod.adjoint,
-    }.items():
-        amici_solver.setSensitivityMethod(sensi_mode)
-        if sensi_mode == amici.SensitivityMethod.none:
-            amici_solver.setSensitivityOrder(amici.SensitivityOrder.none)
-        else:
-            amici_solver.setSensitivityOrder(amici.SensitivityOrder.first)
-
-        res_repeats = [
-            simulate_petab(
-                petab_problem=problem,
-                amici_model=amici_model,
-                solver=amici_solver,
-                log_level=loglevel,
-            )
-            for _ in range(3)  # repeat to get more stable timings
-        ]
-        res = res_repeats[0]
-
-        times[label] = np.min(
-            [
-                sum(r.cpu_time + r.cpu_timeB for r in res[RDATAS]) / 1000
-                # only forwards/backwards simulation
-                for res in res_repeats
-            ]
-        )
-
-        if sensi_mode == amici.SensitivityMethod.none:
-            rdatas = res[RDATAS]
-            llh = res[LLH]
-
-    times["np"] = sum(problem.parameter_df[petab.ESTIMATE])
-
-    pd.Series(times).to_csv(script_dir / f"{args.model_name}_benchmark.csv")
-    for rdata in rdatas:
-        assert (
-            rdata.status == amici.AMICI_SUCCESS
-        ), f"Simulation failed for {rdata.id}"
-
-    # create simulation PEtab table
-    sim_df = rdatas_to_measurement_df(
-        rdatas=rdatas, model=amici_model, measurement_df=problem.measurement_df
-    )
-    sim_df.rename(columns={petab.MEASUREMENT: petab.SIMULATION}, inplace=True)
-
-    if args.simulation_file:
-        sim_df.to_csv(args.simulation_file, index=False, sep="\t")
-
-    if args.plot:
-        with contextlib.suppress(NotImplementedError):
-            # visualize fit
-            axs = plot_problem(petab_problem=problem, simulations_df=sim_df)
-
-            # save figure
-            for plot_id, ax in axs.items():
-                fig_path = os.path.join(
-                    args.model_directory,
-                    f"{args.model_name}_{plot_id}_vis.png",
-                )
-                logger.info(f"Saving figure to {fig_path}")
-                ax.get_figure().savefig(fig_path, dpi=150)
-
-    if args.check:
-        references_yaml = script_dir / "benchmark_models.yaml"
-        with open(references_yaml) as f:
-            refs = yaml.full_load(f)
-
-        try:
-            ref_llh = refs[args.model_name]["llh"]
-
-            rdiff = np.abs((llh - ref_llh) / ref_llh)
-            rtol = 1e-3
-            adiff = np.abs(llh - ref_llh)
-            atol = 1e-3
-            tolstr = (
-                f" Absolute difference is {adiff:.2e} "
-                f"(tol {atol:.2e}) and relative difference is "
-                f"{rdiff:.2e} (tol {rtol:.2e})."
-            )
-
-            if np.isclose(llh, ref_llh, rtol=rtol, atol=atol):
-                logger.info(
-                    f"Computed llh {llh:.4e} matches reference {ref_llh:.4e}."
-                    + tolstr
-                )
-            else:
-                logger.error(
-                    f"Computed llh {llh:.4e} does not match reference "
-                    f"{ref_llh:.4e}." + tolstr
-                )
-                sys.exit(1)
-        except KeyError:
-            logger.error(
-                "No reference likelihood found for "
-                f"{args.model_name} in {references_yaml}"
-            )
-
-        for label, key in {
-            "simulation": "t_sim",
-            "adjoint sensitivity": "t_adj",
-            "forward sensitivity": "t_fwd",
-        }.items():
-            try:
-                ref = refs[args.model_name][key]
-                if times[key] > ref:
-                    logger.error(
-                        f"Computation time for {label} ({times[key]:.2e}) "
-                        f"exceeds reference ({ref:.2e})."
-                    )
-                    sys.exit(1)
-                else:
-                    logger.info(
-                        f"Computation time for {label} ({times[key]:.2e}) "
-                        f"within reference ({ref:.2e})."
-                    )
-            except KeyError:
-                logger.error(
-                    f"No reference time for {label} found for "
-                    f"{args.model_name} in {references_yaml}"
-                )
-
-
-if __name__ == "__main__":
-    main()

From ff6b3f1f37fcb36bcd2ff3e065225a4bd891a14c Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Fri, 25 Oct 2024 09:37:55 +0200
Subject: [PATCH 06/16] Fix problem ID and output in benchmark collection tests
 (#2570)

* Reference value wasn't found for Fiedler_BMCSystBiol2016, because it was recently renamed
* Print all rows from gradient check results
---
 tests/benchmark-models/benchmark_models.yaml  |  2 +-
 .../benchmark-models/test_petab_benchmark.py  | 29 +++++++++++++------
 2 files changed, 21 insertions(+), 10 deletions(-)

diff --git a/tests/benchmark-models/benchmark_models.yaml b/tests/benchmark-models/benchmark_models.yaml
index 8550c9ca38..2fed74cf3f 100644
--- a/tests/benchmark-models/benchmark_models.yaml
+++ b/tests/benchmark-models/benchmark_models.yaml
@@ -57,7 +57,7 @@ Elowitz_Nature2000:
   t_adj: 0.11
   note: benchmark collection reference value matches up to sign when applying log10-correction +sum(log(meas*log(10))) / 2
 
-Fiedler_BMC2016:
+Fiedler_BMCSystBiol2016:
   llh: 58.58390161681
   t_sim: 0.005
   t_fwd: 0.05
diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py
index b4e1f50e68..b055f8961f 100644
--- a/tests/benchmark-models/test_petab_benchmark.py
+++ b/tests/benchmark-models/test_petab_benchmark.py
@@ -35,11 +35,19 @@
 from petab.v1.visualize import plot_problem
 
 
-logger = get_logger(f"amici.{__name__}", logging.WARNING)
+# Enable various debug output
+debug = False
+
+logger = get_logger(
+    f"amici.{__name__}", logging.DEBUG if debug else logging.INFO
+)
 
 script_dir = Path(__file__).parent.absolute()
 repo_root = script_dir.parent.parent
 benchmark_outdir = repo_root / "test_bmc"
+debug_path = script_dir / "debug"
+if debug:
+    debug_path.mkdir(exist_ok=True, parents=True)
 
 # reference values for simulation times and log-likelihoods
 references_yaml = script_dir / "benchmark_models.yaml"
@@ -225,12 +233,6 @@ class GradientCheckSettings:
 )
 
 
-debug = False
-if debug:
-    debug_path = Path(__file__).parent / "debug"
-    debug_path.mkdir(exist_ok=True, parents=True)
-
-
 @pytest.fixture(scope="session", params=problems, ids=problems)
 def benchmark_problem(request):
     """Fixture providing model and PEtab problem for a problem from
@@ -570,9 +572,18 @@ def assert_gradient_check_success(
     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):
+
+    success_fail = "succeeded" if check_result.success else "failed"
+    with pd.option_context(
+        "display.max_columns",
+        None,
+        "display.width",
+        None,
+        "display.max_rows",
+        None,
+    ):
         message = (
-            f"Gradient check failed:\n{df}\n\n"
+            f"Gradient check {success_fail}:\n{df}\n\n"
             f"Maximum absolute difference: {max_adiff} (tolerance: {atol})\n"
             f"Maximum relative difference: {max_rdiff} (tolerance: {rtol})"
         )

From 2b05a19952a81f409336a9cff9364535bc6e62ff Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Wed, 6 Nov 2024 15:20:05 +0100
Subject: [PATCH 07/16] GHA: Workaround for macos segfaults (#2572)

Ignore warnings during macos pytest runs.
---
 .github/workflows/test_python_cplusplus.yml | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml
index 284a170031..fcb1067b85 100644
--- a/.github/workflows/test_python_cplusplus.yml
+++ b/.github/workflows/test_python_cplusplus.yml
@@ -293,7 +293,8 @@ jobs:
 
     - name: Python tests
       run: |
-        scripts/run-python-tests.sh \
+        # ignore warnings until https://github.com/swig/swig/issues/3061 is resolved
+        scripts/run-python-tests.sh -W ignore:: \
           test_pregenerated_models.py \
           test_splines_short.py \
           test_misc.py
@@ -355,7 +356,8 @@ jobs:
 
     - name: Python tests
       run: |
-        scripts/run-python-tests.sh \
+        # ignore warnings until https://github.com/swig/swig/issues/3061 is resolved
+        scripts/run-python-tests.sh -W ignore:: \
           --ignore=test_pregenerated_models.py \
           --ignore=test_splines_short.py \
           --ignore=test_misc.py

From 0009c4ad8f8c6e330df165e8a3f0e4e128e640a4 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Wed, 6 Nov 2024 18:04:57 +0100
Subject: [PATCH 08/16] Check for unsupported observable IDs in sigma
 expressions (#2563)

Closes #2561.
---
 python/sdist/amici/sbml_import.py | 20 +++++++++++++++++---
 1 file changed, 17 insertions(+), 3 deletions(-)

diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py
index 61ce9a0ee1..fcaa1ed752 100644
--- a/python/sdist/amici/sbml_import.py
+++ b/python/sdist/amici/sbml_import.py
@@ -2052,12 +2052,26 @@ def _process_log_likelihood(
                 obs["reg_symbol"] = generate_regularization_symbol(obs_id)
 
         if not event_reg:
+            sigmas = {
+                obs_id: self._sympy_from_sbml_math(sigma)
+                for obs_id, sigma in sigmas.items()
+            }
+            obs_syms = set(self.symbols[obs_symbol].keys())
+            for obs_id in self.symbols[obs_symbol]:
+                if (sigma := sigmas.get(str(obs_id))) and sigma.has(
+                    *(obs_syms - {obs_id})
+                ):
+                    raise ValueError(
+                        f"Sigma expression for {obs_id} ({sigma}) must not "
+                        f"contain any observable symbols other than {obs_id}."
+                    )
+            # check that only the corresponding observable ID is used in the
+            #  sigma formula, but not any other observable ID
+            #  https://github.com/AMICI-dev/AMICI/issues/2561
             self.symbols[sigma_symbol] = {
                 symbol_with_assumptions(f"sigma_{obs_id}"): {
                     "name": f'sigma_{obs["name"]}',
-                    "value": self._sympy_from_sbml_math(
-                        sigmas.get(str(obs_id), "1.0")
-                    ),
+                    "value": sigmas.get(str(obs_id), sp.Float(1.0)),
                 }
                 for obs_id, obs in self.symbols[obs_symbol].items()
             }

From 94b616fd44214ac1f86d8d0eca0f4a948246442b Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Sun, 10 Nov 2024 09:54:10 +0100
Subject: [PATCH 09/16] GHA: require Python3.13 tests to pass (#2577)

It looks like all our dependencies are now Python3.13-ready
(https://github.com/AMICI-dev/AMICI/actions/runs/11753436377/job/32746292107),
so we should require Python3.13 tests to pass.
---
 .github/workflows/test_python_ver_matrix.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test_python_ver_matrix.yml b/.github/workflows/test_python_ver_matrix.yml
index 8ccf2c1558..a7d2d0e514 100644
--- a/.github/workflows/test_python_ver_matrix.yml
+++ b/.github/workflows/test_python_ver_matrix.yml
@@ -33,7 +33,7 @@ jobs:
           - python-version: '3.12'
             experimental: false
           - python-version: '3.13'
-            experimental: true
+            experimental: false
 
     steps:
     - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV

From 9e7f06ba02117bb1dc36a9717bbaf07fccb64809 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Sun, 10 Nov 2024 09:54:23 +0100
Subject: [PATCH 10/16] Fix potentially incorrect disabling of Newton's method
 by `SteadyStateSensitivityMode::integrationOnly` (#2576)

Fixes a bug that would result in disabling Newton's method for
steady-state computation when steady-state sensitivity mode is set
to `integrationOnly`, even if no senitivities are computed.

This was because only the sensitivity method was checked, but not the sensitivity order.

Fixes #2575.
---
 src/steadystateproblem.cpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp
index d2923d5bcd..ce8096dee8 100644
--- a/src/steadystateproblem.cpp
+++ b/src/steadystateproblem.cpp
@@ -133,7 +133,8 @@ void SteadystateProblem::findSteadyState(
         = model.getSteadyStateComputationMode()
               == SteadyStateComputationMode::integrationOnly
           || solver.getNewtonMaxSteps() == 0
-          || (model.getSteadyStateSensitivityMode()
+          || (solver.getSensitivityOrder() >= SensitivityOrder::first
+              && model.getSteadyStateSensitivityMode()
                   == SteadyStateSensitivityMode::integrationOnly
               && ((it == -1
                    && solver.getSensitivityMethodPreequilibration()

From a1cdcb8a4ca32702b01b958329fe62c021585f4b Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Sun, 10 Nov 2024 10:38:19 +0100
Subject: [PATCH 11/16] Change default steady-state method to `integrationOnly`
 (#2574)

**Breaking change**

Change the default mode for computing steady states and sensitivities
at steady state to `integrationOnly` (from previously `integrateIfNewtonFails`).

This is done for a more robust default behaviour. For example, the
evaluation in https://doi.org/10.1371/journal.pone.0312148 shows
that - at least for some models - Newton's method may easily lead to
physically impossible solutions.

To keep the previous behaviour, use:
```python
amici_model.setSteadyStateComputationMode(amici.SteadyStateComputationMode.integrateIfNewtonFails)
amici_model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.integrateIfNewtonFails)
```

Closes #2571.
---
 include/amici/model.h                               | 4 ++--
 python/tests/test_compare_conservation_laws_sbml.py | 2 ++
 python/tests/test_preequilibration.py               | 7 ++++++-
 python/tests/test_pregenerated_models.py            | 7 +++++++
 python/tests/test_swig_interface.py                 | 8 ++++----
 5 files changed, 21 insertions(+), 7 deletions(-)

diff --git a/include/amici/model.h b/include/amici/model.h
index 533f0139dc..57591dd2f7 100644
--- a/include/amici/model.h
+++ b/include/amici/model.h
@@ -2063,12 +2063,12 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /** method for steady-state computation */
     SteadyStateComputationMode steadystate_computation_mode_{
-        SteadyStateComputationMode::integrateIfNewtonFails
+        SteadyStateComputationMode::integrationOnly
     };
 
     /** method for steadystate sensitivities computation */
     SteadyStateSensitivityMode steadystate_sensitivity_mode_{
-        SteadyStateSensitivityMode::integrateIfNewtonFails
+        SteadyStateSensitivityMode::integrationOnly
     };
 
     /**
diff --git a/python/tests/test_compare_conservation_laws_sbml.py b/python/tests/test_compare_conservation_laws_sbml.py
index 640d2dd988..a19a7565d2 100644
--- a/python/tests/test_compare_conservation_laws_sbml.py
+++ b/python/tests/test_compare_conservation_laws_sbml.py
@@ -102,6 +102,7 @@ def get_results(
     sensi_order=0,
     sensi_meth=amici.SensitivityMethod.forward,
     sensi_meth_preeq=amici.SensitivityMethod.forward,
+    stst_mode=amici.SteadyStateComputationMode.integrateIfNewtonFails,
     stst_sensi_mode=amici.SteadyStateSensitivityMode.newtonOnly,
     reinitialize_states=False,
 ):
@@ -115,6 +116,7 @@ def get_results(
     solver.setSensitivityMethodPreequilibration(sensi_meth_preeq)
     solver.setSensitivityMethod(sensi_meth)
     model.setSteadyStateSensitivityMode(stst_sensi_mode)
+    model.setSteadyStateComputationMode(stst_mode)
     if edata is None:
         model.setTimepoints(np.linspace(0, 5, 101))
     else:
diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py
index bf81e0256b..c7e254114e 100644
--- a/python/tests/test_preequilibration.py
+++ b/python/tests/test_preequilibration.py
@@ -18,7 +18,12 @@
 def preeq_fixture(pysb_example_presimulation_module):
     model = pysb_example_presimulation_module.getModel()
     model.setReinitializeFixedParameterInitialStates(True)
-
+    model.setSteadyStateComputationMode(
+        amici.SteadyStateComputationMode.integrateIfNewtonFails
+    )
+    model.setSteadyStateSensitivityMode(
+        amici.SteadyStateSensitivityMode.integrateIfNewtonFails
+    )
     solver = model.getSolver()
     solver.setSensitivityOrder(amici.SensitivityOrder.first)
     solver.setSensitivityMethod(amici.SensitivityMethod.forward)
diff --git a/python/tests/test_pregenerated_models.py b/python/tests/test_pregenerated_models.py
index 5a110cdfc2..e9daf13ba8 100755
--- a/python/tests/test_pregenerated_models.py
+++ b/python/tests/test_pregenerated_models.py
@@ -63,6 +63,13 @@ def test_pregenerated_model(sub_test, case):
     amici.readModelDataFromHDF5(
         options_file, model.get(), f"/{sub_test}/{case}/options"
     )
+    if model_name == "model_steadystate":
+        model.setSteadyStateComputationMode(
+            amici.SteadyStateComputationMode.integrateIfNewtonFails
+        )
+        model.setSteadyStateSensitivityMode(
+            amici.SteadyStateSensitivityMode.integrateIfNewtonFails
+        )
     amici.readSolverSettingsFromHDF5(
         options_file, solver.get(), f"/{sub_test}/{case}/options"
     )
diff --git a/python/tests/test_swig_interface.py b/python/tests/test_swig_interface.py
index f61b28ea4d..9686a25d94 100644
--- a/python/tests/test_swig_interface.py
+++ b/python/tests/test_swig_interface.py
@@ -105,12 +105,12 @@ def test_copy_constructors(pysb_example_presimulation_module):
     # `pysb_example_presimulation_module.getModel()`.
     "StateIsNonNegative": None,
     "SteadyStateComputationMode": [
-        2,
-        1,
+        amici.SteadyStateComputationMode.integrationOnly,
+        amici.SteadyStateComputationMode.integrateIfNewtonFails,
     ],
     "SteadyStateSensitivityMode": [
-        2,
-        1,
+        amici.SteadyStateSensitivityMode.integrationOnly,
+        amici.SteadyStateSensitivityMode.integrateIfNewtonFails,
     ],
     ("t0", "setT0"): [
         0.0,

From 2c362b432b6da97c3f612be6ca4a0f5d692186b5 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@uni-bonn.de>
Date: Sun, 10 Nov 2024 14:52:40 +0100
Subject: [PATCH 12/16] clang-format

---
 src/solver_cvodes.cpp      | 6 ++++--
 src/steadystateproblem.cpp | 2 +-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp
index dfccbb2389..0b60304c55 100644
--- a/src/solver_cvodes.cpp
+++ b/src/solver_cvodes.cpp
@@ -34,8 +34,10 @@ static_assert((int)InterpolationType::polynomial == CV_POLYNOMIAL, "");
 static_assert((int)LinearMultistepMethod::adams == CV_ADAMS, "");
 static_assert((int)LinearMultistepMethod::BDF == CV_BDF, "");
 
-#define STATIC_ASSERT_EQUAL(amici_constant, cv_constant) \
-static_assert(amici_constant == cv_constant, #amici_constant " != " #cv_constant)
+#define STATIC_ASSERT_EQUAL(amici_constant, cv_constant)                       \
+    static_assert(                                                             \
+        amici_constant == cv_constant, #amici_constant " != " #cv_constant     \
+    )
 
 STATIC_ASSERT_EQUAL(amici::AMICI_SUCCESS, CV_SUCCESS);
 STATIC_ASSERT_EQUAL(amici::AMICI_ROOT_RETURN, CV_ROOT_RETURN);
diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp
index ce8096dee8..7c3ac9fab3 100644
--- a/src/steadystateproblem.cpp
+++ b/src/steadystateproblem.cpp
@@ -135,7 +135,7 @@ void SteadystateProblem::findSteadyState(
           || solver.getNewtonMaxSteps() == 0
           || (solver.getSensitivityOrder() >= SensitivityOrder::first
               && model.getSteadyStateSensitivityMode()
-                  == SteadyStateSensitivityMode::integrationOnly
+                     == SteadyStateSensitivityMode::integrationOnly
               && ((it == -1
                    && solver.getSensitivityMethodPreequilibration()
                           == SensitivityMethod::forward)

From 40834389e826e8b38f253820438f375bb1971a20 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Sun, 10 Nov 2024 17:53:26 +0100
Subject: [PATCH 13/16] Optional warning in `fill_in_parameters` (#2578)

`amici.petab.conditions.fill_in_parameters` emits a warnings if parameters are supplied that don't occur in the parameter mapping.
This is to point out potential issues with the parameter mapping.
However, sometimes it's more convenient to silently ignore those extra parameters. Therefore, make this warning optional.
---
 python/sdist/amici/petab/conditions.py | 22 +++++++++++++++++++---
 1 file changed, 19 insertions(+), 3 deletions(-)

diff --git a/python/sdist/amici/petab/conditions.py b/python/sdist/amici/petab/conditions.py
index 08c2f90302..ab06e8850d 100644
--- a/python/sdist/amici/petab/conditions.py
+++ b/python/sdist/amici/petab/conditions.py
@@ -41,6 +41,7 @@ def fill_in_parameters(
     scaled_parameters: bool,
     parameter_mapping: ParameterMapping,
     amici_model: AmiciModel,
+    warn_unused: bool = True,
 ) -> None:
     """Fill fixed and dynamic parameters into the edatas (in-place).
 
@@ -59,9 +60,15 @@ def fill_in_parameters(
         Parameter mapping for all conditions.
     :param amici_model:
         AMICI model.
+    :param warn_unused:
+        Whether a warning should be emitted if not all problem parameters
+        were used. I.e., if there are parameters in `problem_parameters`
+        that are not in `parameter_mapping`.
     """
-    if unused_parameters := (
-        set(problem_parameters.keys()) - parameter_mapping.free_symbols
+    if warn_unused and (
+        unused_parameters := (
+            set(problem_parameters.keys()) - parameter_mapping.free_symbols
+        )
     ):
         warnings.warn(
             "The following problem parameters were not used: "
@@ -223,6 +230,7 @@ def create_parameterized_edatas(
     scaled_parameters: bool = False,
     parameter_mapping: ParameterMapping = None,
     simulation_conditions: pd.DataFrame | dict = None,
+    warn_unused: bool = True,
 ) -> list[amici.ExpData]:
     """Create list of :class:amici.ExpData objects with parameters filled in.
 
@@ -244,6 +252,11 @@ def create_parameterized_edatas(
     :param simulation_conditions:
         Result of :func:`petab.get_simulation_conditions`. Can be provided to
         save time if this has been obtained before.
+    :param warn_unused:
+        Whether a warning should be emitted if not all problem parameters
+        were used. I.e., if there are parameters in `problem_parameters`
+        that are not in `parameter_mapping` or in the generated parameter
+        mapping.
 
     :return:
         List with one :class:`amici.amici.ExpData` per simulation condition,
@@ -282,6 +295,7 @@ def create_parameterized_edatas(
         scaled_parameters=scaled_parameters,
         parameter_mapping=parameter_mapping,
         amici_model=amici_model,
+        warn_unused=warn_unused,
     )
 
     return edatas
@@ -387,7 +401,9 @@ def create_edatas(
 
     :return:
         List with one :class:`amici.amici.ExpData` per simulation condition,
-        with filled in timepoints and data.
+        with filled in timepoints and data, but without parameter values
+        (see :func:`create_parameterized_edatas` or
+        :func:`fill_in_parameters` for that).
     """
     if simulation_conditions is None:
         simulation_conditions = (

From e84a04284bd066f72370b88c735ad508d63c5207 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Sun, 10 Nov 2024 17:55:38 +0100
Subject: [PATCH 14/16] Update codecov config (#2580)

Update `fixes`, ignore test models.
---
 codecov.yml | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/codecov.yml b/codecov.yml
index 15beea9dfc..2779b486d1 100644
--- a/codecov.yml
+++ b/codecov.yml
@@ -1,8 +1,10 @@
-fixes:
-  - "build/venv/lib/python3.9/site-packages/::python/"
-  - "build/venv/lib/python3.10/site-packages/::python/"
-  - "build/venv/lib/python3.11/site-packages/::python/"
+# see https://docs.codecov.com/docs/codecovyml-reference
 
+fixes:
+  # https://docs.codecov.com/docs/fixing-paths
+  - "build/venv/lib/python[0-9.]+/site-packages/::python/"
+  - "python/sdist/build/temp_amici/CMakeFiles/amici.dir/src/::src/"
+  - "build/CMakeFiles/amici.dir/src/::src/"
 codecov:
   require_ci_to_pass: yes
 
@@ -27,6 +29,8 @@ comment:
 ignore:
   - "tests/*"
   - "tests/**/*"
+  - "build/tests/**"
+  - "amici_models/**"
 
 flags:
   python:

From 5c7079489d97f220bd7d0f85cca4b8865b55a125 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Sun, 10 Nov 2024 17:56:38 +0100
Subject: [PATCH 15/16] Fix SUNContext in ModelStateDerived copy ctor (#2583)

When copying SUNDIALS objects, we need to make sure they have a valid `SUNContext` set.
When copying `ModelStateDerived`, the `SUNContext` needs to be updated. Otherwise,
dangling SUNContext pointers likely result in segfaults later on (#2579).

This was done only for a subset of objects. In particular the ones added only during
https://github.com/AMICI-dev/AMICI/pull/2505 were missing, because of other changes
happening in parallel.

To make this less messy: Add `set_ctx(.)` functions all our sundials vector/matrix
wrappers.

Closes #2579.
---
 include/amici/model_state.h             | 90 ++++++++-----------------
 include/amici/sundials_matrix_wrapper.h | 13 ++++
 include/amici/vector.h                  | 14 ++++
 3 files changed, 55 insertions(+), 62 deletions(-)

diff --git a/include/amici/model_state.h b/include/amici/model_state.h
index 149c5757b8..defb12d4c0 100644
--- a/include/amici/model_state.h
+++ b/include/amici/model_state.h
@@ -146,71 +146,37 @@ struct ModelStateDerived {
         , dwdw_(other.dwdw_)
         , dwdx_hierarchical_(other.dwdx_hierarchical_)
         , dJydy_dense_(other.dJydy_dense_) {
-        // Update the SUNContext of all matrices
-        if (J_.data()) {
-            J_.get()->sunctx = sunctx_;
+        // Update the SUNContext of all SUNDIALS objects
+        J_.set_ctx(sunctx_);
+        JB_.set_ctx(sunctx_);
+        dxdotdw_.set_ctx(sunctx_);
+        dwdx_.set_ctx(sunctx_);
+        dwdp_.set_ctx(sunctx_);
+        M_.set_ctx(sunctx_);
+        MSparse_.set_ctx(sunctx_);
+        dfdx_.set_ctx(sunctx_);
+        dxdotdp_full.set_ctx(sunctx_);
+        dxdotdp_explicit.set_ctx(sunctx_);
+        dxdotdp_implicit.set_ctx(sunctx_);
+        dxdotdx_explicit.set_ctx(sunctx_);
+        dxdotdx_implicit.set_ctx(sunctx_);
+        dx_rdatadx_solver.set_ctx(sunctx_);
+        dx_rdatadtcl.set_ctx(sunctx_);
+        dtotal_cldx_rdata.set_ctx(sunctx_);
+        dxdotdp.set_ctx(sunctx_);
+
+        for (auto& dJydy : dJydy_) {
+            dJydy.set_ctx(sunctx_);
         }
-        if (JB_.data()) {
-            JB_.get()->sunctx = sunctx_;
+        for (auto& dwdp : dwdp_hierarchical_) {
+            dwdp.set_ctx(sunctx_);
         }
-        if (dxdotdw_.data()) {
-            dxdotdw_.get()->sunctx = sunctx_;
-        }
-        if (dwdx_.data()) {
-            dwdx_.get()->sunctx = sunctx_;
-        }
-        if (dwdp_.data()) {
-            dwdp_.get()->sunctx = sunctx_;
-        }
-        if (M_.data()) {
-            M_.get()->sunctx = sunctx_;
-        }
-        if (MSparse_.data()) {
-            MSparse_.get()->sunctx = sunctx_;
-        }
-        if (dfdx_.data()) {
-            dfdx_.get()->sunctx = sunctx_;
-        }
-        if (dxdotdp_full.data()) {
-            dxdotdp_full.get()->sunctx = sunctx_;
-        }
-        if (dxdotdp_explicit.data()) {
-            dxdotdp_explicit.get()->sunctx = sunctx_;
-        }
-        if (dxdotdp_implicit.data()) {
-            dxdotdp_implicit.get()->sunctx = sunctx_;
-        }
-        if (dxdotdx_explicit.data()) {
-            dxdotdx_explicit.get()->sunctx = sunctx_;
-        }
-        if (dxdotdx_implicit.data()) {
-            dxdotdx_implicit.get()->sunctx = sunctx_;
-        }
-        if (dx_rdatadx_solver.data()) {
-            dx_rdatadx_solver.get()->sunctx = sunctx_;
-        }
-        if (dx_rdatadtcl.data()) {
-            dx_rdatadtcl.get()->sunctx = sunctx_;
-        }
-        if (dtotal_cldx_rdata.data()) {
-            dtotal_cldx_rdata.get()->sunctx = sunctx_;
-        }
-        for (auto const& dwdp : dwdp_hierarchical_) {
-            if (dwdp.data()) {
-                dwdp.get()->sunctx = sunctx_;
-            }
-        }
-        for (auto const& dwdx : dwdx_hierarchical_) {
-            if (dwdx.data()) {
-                dwdx.get()->sunctx = sunctx_;
-            }
-        }
-        if (dwdw_.data()) {
-            dwdw_.get()->sunctx = sunctx_;
-        }
-        if (dJydy_dense_.data()) {
-            dJydy_dense_.get()->sunctx = sunctx_;
+        for (auto& dwdx : dwdx_hierarchical_) {
+            dwdx.set_ctx(sunctx_);
         }
+        sspl_.set_ctx(sunctx_);
+        dwdw_.set_ctx(sunctx_);
+        dJydy_dense_.set_ctx(sunctx_);
     }
 
     /**
diff --git a/include/amici/sundials_matrix_wrapper.h b/include/amici/sundials_matrix_wrapper.h
index 942069a803..a7ad361acb 100644
--- a/include/amici/sundials_matrix_wrapper.h
+++ b/include/amici/sundials_matrix_wrapper.h
@@ -506,6 +506,19 @@ class SUNMatrixWrapper {
      */
     SUNContext get_ctx() const;
 
+    /**
+     * @brief Set SUNContext
+     *
+     * Update the SUNContext of the wrapped SUNMatrix.
+     *
+     * @param ctx SUNDIALS context to set
+     */
+    void set_ctx(SUNContext ctx) {
+        if (matrix_) {
+            matrix_->sunctx = ctx;
+        }
+    }
+
   private:
     /**
      * @brief SUNMatrix to which all methods are applied
diff --git a/include/amici/vector.h b/include/amici/vector.h
index 32c436fbda..0a7648e460 100644
--- a/include/amici/vector.h
+++ b/include/amici/vector.h
@@ -413,6 +413,20 @@ class AmiVectorArray {
      */
     void copy(AmiVectorArray const& other);
 
+    /**
+     * @brief Set SUNContext
+     *
+     * If any AmiVector is non-empty, this changes the current SUNContext of the
+     * associated N_Vector. If empty, do nothing.
+     *
+     * @param ctx SUNDIALS context to set
+     */
+    void set_ctx(SUNContext ctx) {
+        for (auto& vec : vec_array_) {
+            vec.set_ctx(ctx);
+        }
+    }
+
   private:
     /** main data storage */
     std::vector<AmiVector> vec_array_;

From 307fb2380608b7a315e4102a5327512d47efa90a Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@uni-bonn.de>
Date: Sun, 10 Nov 2024 11:35:34 +0100
Subject: [PATCH 16/16] Update changelog, bump version

---
 CHANGELOG.md | 57 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 version.txt  |  2 +-
 2 files changed, 58 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8eb0abd8ae..d4c667b7c7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,63 @@ See also our [versioning policy](https://amici.readthedocs.io/en/latest/versioni
 
 ## v0.X Series
 
+### v0.28.0 (2024-11-11)
+
+**Breaking changes**
+
+* Changed the default steady-state method to `integrationOnly`
+  (by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2574)
+
+  The default mode for computing steady states and sensitivities at steady
+  state was changed to `integrationOnly`
+  (from previously `integrateIfNewtonFails`).
+
+  This was done for a more robust default behavior. For example, the evaluation
+  in https://doi.org/10.1371/journal.pone.0312148 shows that - at least for
+  some models - Newton's method may easily lead to physically impossible
+  solutions.
+
+  To keep the previous behavior, use:
+  ```python
+  amici_model.setSteadyStateComputationMode(amici.SteadyStateComputationMode.integrateIfNewtonFails)
+  amici_model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.integrateIfNewtonFails)
+  ```
+
+**Fixes**
+
+* PEtab import: **Fixed potentially incorrect sensitivities** with
+  observable/state-dependent sigmas.
+  This was fixed for all cases amici can handle, others cases will now result
+  in `ValueError`s (https://github.com/AMICI-dev/AMICI/pull/2563).
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2562
+
+* Fixed potentially incorrect disabling of Newton's method
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2576
+
+* Fixed `ModelStateDerived` copy ctor, where previously dangling pointers
+  could lead to crashes in some situations
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2583
+
+* Added missing simulation status codes
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2560
+
+* Check for unsupported observable IDs in sigma expressions
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2563
+
+
+**Features**
+
+* Optional warning in `fill_in_parameters`
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2578
+
+**Full Changelog**: https://github.com/AMICI-dev/AMICI/compare/v0.27.0...v0.28.0
+
 ### v0.27.0 (2024-10-21)
 
 This release comes with an **updated version of the SUNDIALS package (7.1.1)** (https://github.com/AMICI-dev/AMICI/pull/2513).
diff --git a/version.txt b/version.txt
index 1b58cc1018..697f087f37 100644
--- a/version.txt
+++ b/version.txt
@@ -1 +1 @@
-0.27.0
+0.28.0