Add sllh computation back to petab_objective.simulate_petab (#1548)

* add sllh to `simulate_petab`

* use ExpData for plist if available

* use plist if set; breakpoint at error; preliminary test

* add sensi test to preexisting benchmark model test

* remove outdated test and test cases

* add FD check for petab problems; add petab problem test case for SLLH; clean up

* typo, RTOL

* fix paths

* docs

* Apply suggestions from code review

Co-authored-by: Yannik Schälte <31767307+yannikschaelte@users.noreply.github.com>
Co-authored-by: Daniel Weindl <dweindl@users.noreply.github.com>

* reviews: docstring;error handling;optional sllh aggregation

* Apply suggestions from code review

Co-authored-by: Fabian Fröhlich <fabian_froehlich@hms.harvard.edu>

* review: clean up

* ignore parameters missing due to estimate=0

* add test for petab benchmark models

* use fiddy

* develop merge

* fix after merge

* merge conflicts

* restore changes

* make scaled gradients optional, default to linear

* test scaled and unscaled

* add fiddy to installation

* fix deps

* update for fiddy redesign

* fix fiddy version

* update benchmark settings

* move tests

* add fiddy

* Update test_benchmark_collection_models.yml

* fix paths

* Update test_petab_benchmark.py

* restrict test set

* Update test_petab_benchmark.py

* Update test_petab_benchmark.py

* Update test_petab_benchmark.py

* Update test_petab_benchmark.py

* use fiddy main

* Update test_petab_benchmark.py

* oopsie

* address code review

* remove custom grad check methods; adjust test

* Update python/sdist/amici/petab_objective.py

Co-authored-by: Dilan Pathirana <59329744+dilpath@users.noreply.github.com>

---------

Co-authored-by: Daniel Weindl <dweindl@users.noreply.github.com>
Co-authored-by: Yannik Schälte <31767307+yannikschaelte@users.noreply.github.com>
Co-authored-by: Fabian Fröhlich <fabian_froehlich@hms.harvard.edu>
Co-authored-by: Daniel Weindl <daniel.weindl@helmholtz-muenchen.de>
Co-authored-by: Fabian Fröhlich <fabian@schaluck.com>
Co-authored-by: Fabian Fröhlich <fabian.frohlich@crick.ac.uk>

.github/workflows/test_benchmark_collection_models.yml

@@ -61,6 +61,12 @@ jobs:
         git clone --depth 1 https://github.com/benchmarking-initiative/Benchmark-Models-PEtab.git \
           && export BENCHMARK_COLLECTION="$(pwd)/Benchmark-Models-PEtab/Benchmark-Models/" \
           && AMICI_PARALLEL_COMPILE=2 tests/benchmark-models/test_benchmark_collection.sh
+          
+    # run gradient checks
+    - name: Run Gradient Checks
+      run: |
+        pip install git+https://github.com/ICB-DCM/fiddy.git \
+          && cd tests/benchmark-models && pytest ./test_petab_benchmark.py
 
     # upload results
     - uses: actions/upload-artifact@v3

.gitignore

@@ -198,3 +198,6 @@ Benchmark-Models-PEtab/
 /test_bmc/
 
 CS_Signalling_ERBB_RAS_AKT/
+cache_fiddy/*
+debug/*
+tests/benchmark-models/cache_fiddy/*

python/sdist/amici/petab_objective.py

@@ -24,7 +24,10 @@
 from .logging import get_logger, log_execution_time
 from .petab_import import PREEQ_INDICATOR_ID, element_is_state
 from .parameter_mapping import (
-    fill_in_parameters, ParameterMappingForCondition, ParameterMapping)
+    fill_in_parameters,
+    ParameterMappingForCondition,
+    ParameterMapping,
+)
 
 logger = get_logger(__name__)
 
@@ -53,9 +56,14 @@ def simulate_petab(
         log_level: int = logging.WARNING,
         num_threads: int = 1,
         failfast: bool = True,
+        scaled_gradients: bool = False,
 ) -> Dict[str, Any]:
     """Simulate PEtab model.
 
+    .. note::
+        Regardless of `scaled_parameters`, unscaled sensitivities are returned,
+        unless `scaled_gradients=True`.
+
     :param petab_problem:
         PEtab problem to work on.
     :param amici_model:
@@ -87,6 +95,9 @@ def simulate_petab(
     :param failfast:
         Returns as soon as an integration failure is encountered, skipping
         any remaining simulations.
+    :param scaled_gradients:
+        Whether to compute gradients on parameter scale (``True``) or not
+        (``False``).
 
     :return:
         Dictionary of
@@ -104,15 +115,13 @@ def simulate_petab(
     if solver is None:
         solver = amici_model.getSolver()
 
-    # Get parameters
-    if problem_parameters is None:
-        # Use PEtab nominal values as default
-        problem_parameters = {t.Index: getattr(t, NOMINAL_VALUE) for t in
-                              petab_problem.parameter_df.itertuples()}
-        if scaled_parameters:
-            raise NotImplementedError(
-                "scaled_parameters=True in combination with "
-                "problem_parameters=None is currently not supported.")
+    # Switch to scaled parameters.
+    problem_parameters = _default_scaled_parameters(
+        petab_problem=petab_problem,
+        problem_parameters=problem_parameters,
+        scaled_parameters=scaled_parameters,
+    )
+    scaled_parameters = True
 
     # number of amici simulations will be number of unique
     # (preequilibrationConditionId, simulationConditionId) pairs.
@@ -153,6 +162,29 @@ def simulate_petab(
 
     # Compute total llh
     llh = sum(rdata['llh'] for rdata in rdatas)
+    # Compute total sllh
+    sllh = None
+    if solver.getSensitivityOrder() != amici.SensitivityOrder.none:
+        sllh = aggregate_sllh(
+            amici_model=amici_model,
+            rdatas=rdatas,
+            parameter_mapping=parameter_mapping,
+            petab_scale=scaled_parameters,
+            petab_problem=petab_problem,
+            edatas=edatas,
+        )
+        if not scaled_gradients and sllh is not None:
+            sllh = {
+                parameter_id: rescale_sensitivity(
+                    sensitivity=sensitivity,
+                    parameter_value=problem_parameters[parameter_id],
+                    old_scale=petab_problem.parameter_df.loc[
+                        parameter_id, PARAMETER_SCALE
+                    ],
+                    new_scale=LIN,
+                )
+                for parameter_id, sensitivity in sllh.items()
+            }
 
     # Log results
     sim_cond = petab_problem.get_simulation_conditions_from_measurement_df()
@@ -165,11 +197,158 @@ def simulate_petab(
 
     return {
         LLH: llh,
+        SLLH: sllh,
         RDATAS: rdatas,
         EDATAS: edatas,
     }
 
 
+def aggregate_sllh(
+        amici_model: AmiciModel,
+        rdatas: Sequence[amici.ReturnDataView],
+        parameter_mapping: Optional[ParameterMapping],
+        edatas: List[AmiciExpData],
+        petab_scale: bool = True,
+        petab_problem: petab.Problem = None,
+) -> Union[None, Dict[str, float]]:
+    """
+    Aggregate likelihood gradient for all conditions, according to PEtab
+    parameter mapping.
+
+    :param amici_model:
+        AMICI model from which ``rdatas`` were obtained.
+    :param rdatas:
+        Simulation results.
+    :param parameter_mapping:
+        PEtab parameter mapping to condition-specific simulation parameters.
+    :param edatas:
+        Experimental data used for simulation.
+    :param petab_scale:
+        Whether to check that sensitivities were computed with parameters on
+        the scales provided in the PEtab parameters table.
+    :param petab_problem:
+        The PEtab problem that defines the parameter scales.
+
+    :return:
+        Aggregated likelihood sensitivities.
+    """
+    accumulated_sllh = {}
+    model_parameter_ids = amici_model.getParameterIds()
+
+    if petab_scale and petab_problem is None:
+        raise ValueError(
+            'Please provide the PEtab problem, when using '
+            '`petab_scale=True`.'
+        )
+
+    # Check for issues in all condition simulation results.
+    for rdata in rdatas:
+        # Condition failed during simulation.
+        if rdata.status != amici.AMICI_SUCCESS:
+            return None
+        # Condition simulation result does not provide SLLH.
+        if rdata.sllh is None:
+            raise ValueError(
+                'The sensitivities of the likelihood for a condition were '
+                'not computed.'
+            )
+
+    for condition_parameter_mapping, edata, rdata in \
+            zip(parameter_mapping, edatas, rdatas):
+        for sllh_parameter_index, condition_parameter_sllh in \
+                enumerate(rdata.sllh):
+            # Get PEtab parameter ID
+            # Use ExpData if it provides a parameter list, else default to
+            # Model.
+            if edata.plist:
+                model_parameter_index = edata.plist[sllh_parameter_index]
+            else:
+                model_parameter_index = amici_model.plist(sllh_parameter_index)
+            model_parameter_id = model_parameter_ids[model_parameter_index]
+            petab_parameter_id = (
+                condition_parameter_mapping
+                .map_sim_var[model_parameter_id]
+            )
+
+            # Initialize
+            if petab_parameter_id not in accumulated_sllh:
+                accumulated_sllh[petab_parameter_id] = 0
+
+            # Check that the scale is consistent
+            if petab_scale:
+                # `ParameterMappingForCondition` objects provide the scale in
+                # terms of `petab.C` constants already, not AMICI equivalents.
+                model_parameter_scale = (
+                    condition_parameter_mapping
+                    .scale_map_sim_var[model_parameter_id]
+                )
+                petab_parameter_scale = (
+                    petab_problem
+                    .parameter_df
+                    .loc[petab_parameter_id, PARAMETER_SCALE]
+                )
+                if model_parameter_scale != petab_parameter_scale:
+                    raise ValueError(
+                        f'The scale of the parameter `{petab_parameter_id}` '
+                        'differs between the AMICI model '
+                        f'({model_parameter_scale}) and the PEtab problem '
+                        f'({petab_parameter_scale}).'
+                    )
+
+            # Accumulate
+            accumulated_sllh[petab_parameter_id] += condition_parameter_sllh
+
+    return accumulated_sllh
+
+
+def rescale_sensitivity(
+    sensitivity: float,
+    parameter_value: float,
+    old_scale: str,
+    new_scale: str,
+) -> float:
+    """Rescale a sensitivity between parameter scales.
+
+    :param sensitivity:
+        The sensitivity corresponding to the parameter value.
+    :param parameter_value:
+        The parameter vector element, on ``old_scale``.
+    :param old_scale:
+        The scale of the parameter value.
+    :param new_scale:
+        The parameter scale on which to rescale the sensitivity.
+
+    :return:
+        The rescaled sensitivity.
+    """
+    LOG_E_10 = np.log(10)
+
+    if old_scale == new_scale:
+        return sensitivity
+
+    unscaled_parameter_value = petab.parameters.unscale(
+        parameter=parameter_value,
+        scale_str=old_scale,
+    )
+
+    scale = {
+        (LIN, LOG): lambda s: s * unscaled_parameter_value,
+        (LOG, LIN): lambda s: s / unscaled_parameter_value,
+        (LIN, LOG10): lambda s: s * (unscaled_parameter_value * LOG_E_10),
+        (LOG10, LIN): lambda s: s / (unscaled_parameter_value * LOG_E_10),
+    }
+
+    scale[(LOG, LOG10)] = lambda s: scale[(LIN, LOG10)](scale[(LOG, LIN)](s))
+    scale[(LOG10, LOG)] = lambda s: scale[(LIN, LOG)](scale[(LOG10, LIN)](s))
+
+    if (old_scale, new_scale) not in scale:
+        raise NotImplementedError(
+            f"Old scale: {old_scale}. New scale: {new_scale}."
+        )
+
+    return scale[(old_scale, new_scale)](sensitivity)
+
+
 def create_parameterized_edatas(
         amici_model: AmiciModel,
         petab_problem: petab.Problem,
@@ -811,3 +990,39 @@ def rdatas_to_simulation_df(
                                   measurement_df=measurement_df)
 
     return df.rename(columns={MEASUREMENT: SIMULATION})
+
+
+def _default_scaled_parameters(
+    petab_problem: petab.Problem,
+    problem_parameters: Optional[Dict[str, float]] = None,
+    scaled_parameters: bool = False,
+) -> Optional[Dict[str, float]]:
+    """
+    Helper method to handle an unscaled or unspecified parameter vector.
+
+    The parameter vector defaults to the nominal values in the PEtab
+    parameter table.
+
+    Unscaled parameter values are scaled.
+
+    :param petab_problem:
+        The PEtab problem.
+    :param problem_parameters:
+        Keys are PEtab parameter IDs, values are parameter values on the scale
+        defined in the PEtab parameter table. Defaults to the nominal values in
+        the PEtab parameter table.
+    :param scaled_parameters:
+        Whether `problem_parameters` are on the scale defined in the PEtab
+        parameter table.
+
+    :return:
+        The scaled parameter vector.
+    """
+    if problem_parameters is None:
+        problem_parameters = dict(zip(
+            petab_problem.x_ids,
+            petab_problem.x_nominal_scaled,
+        ))
+    elif not scaled_parameters:
+        problem_parameters = petab_problem.scale_parameters(problem_parameters)
+    return problem_parameters

python/tests/petab_test_problems/lotka_volterra/README.md

@@ -0,0 +1,5 @@
+A Lotka-Volterra model, based on the model provided as an example in the `yaml2sbml` package: https://github.com/yaml2sbml-dev/yaml2sbml
+
+The PEtab problem can be created by running `bash model/convert_to_petab.sh` (test on Ubuntu 20.04) inside the `model` directory, in a Python 3 environment with `yaml2sbml`: https://pypi.org/project/yaml2sbml/
+
+The model is augmented with new parameters `departure_prey` and `arrival_predator`, to allow for a steady-state under certain conditions. This results in a model that can pre-equilibrate, then switch to the expected oscillations of a Lotka-Volterra model.

python/tests/petab_test_problems/lotka_volterra/model/convert_to_petab.sh

@@ -0,0 +1,24 @@
+#!/usr/bin/env bash
+set -eou pipefail
+script_dir=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )
+petab_path=$script_dir/../petab
+
+mkdir -p "${petab_path}"
+
+# Validate input yaml2sbml model
+yaml2sbml_validate lotka_volterra.yaml
+
+# Copy measurements to PEtab directory
+cp "${script_dir}/measurements.tsv" "$petab_path/measurements.tsv"
+
+# Write the PEtab problem
+python "${script_dir}/writer.py"
+
+# Remove condition parameters from PEtab parameters table
+for condition_parameter in beta delta departure_prey arrival_predator
+do
+	sed -i "/^${condition_parameter}/d" "${petab_path}/parameter*"
+done
+
+# Validate the PEtab problem
+petablint -vy "${petab_path}/problem.yaml"