Skip to content

Commit d2fb35e

Browse files
dilpathdweindlyannikschaelteFabian FröhlichFFroehlich
authored
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 <[email protected]> Co-authored-by: Daniel Weindl <[email protected]> * reviews: docstring;error handling;optional sllh aggregation * Apply suggestions from code review Co-authored-by: Fabian Fröhlich <[email protected]> * 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 <[email protected]> --------- Co-authored-by: Daniel Weindl <[email protected]> Co-authored-by: Yannik Schälte <[email protected]> Co-authored-by: Fabian Fröhlich <[email protected]> Co-authored-by: Daniel Weindl <[email protected]> Co-authored-by: Fabian Fröhlich <[email protected]> Co-authored-by: Fabian Fröhlich <[email protected]>
1 parent 255be19 commit d2fb35e

17 files changed

+793
-11
lines changed

‎.github/workflows/test_benchmark_collection_models.yml‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,12 @@ jobs:
6161
git clone --depth 1 https://github.com/benchmarking-initiative/Benchmark-Models-PEtab.git \
6262
&& export BENCHMARK_COLLECTION="$(pwd)/Benchmark-Models-PEtab/Benchmark-Models/" \
6363
&& AMICI_PARALLEL_COMPILE=2 tests/benchmark-models/test_benchmark_collection.sh
64+
65+
# run gradient checks
66+
- name: Run Gradient Checks
67+
run: |
68+
pip install git+https://github.com/ICB-DCM/fiddy.git \
69+
&& cd tests/benchmark-models && pytest ./test_petab_benchmark.py
6470
6571
# upload results
6672
- uses: actions/upload-artifact@v3

‎.gitignore‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -198,3 +198,6 @@ Benchmark-Models-PEtab/
198198
/test_bmc/
199199

200200
CS_Signalling_ERBB_RAS_AKT/
201+
cache_fiddy/*
202+
debug/*
203+
tests/benchmark-models/cache_fiddy/*

‎python/sdist/amici/petab_objective.py‎

Lines changed: 225 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,10 @@
2424
from .logging import get_logger, log_execution_time
2525
from .petab_import import PREEQ_INDICATOR_ID, element_is_state
2626
from .parameter_mapping import (
27-
fill_in_parameters, ParameterMappingForCondition, ParameterMapping)
27+
fill_in_parameters,
28+
ParameterMappingForCondition,
29+
ParameterMapping,
30+
)
2831

2932
logger = get_logger(__name__)
3033

@@ -53,9 +56,14 @@ def simulate_petab(
5356
log_level: int = logging.WARNING,
5457
num_threads: int = 1,
5558
failfast: bool = True,
59+
scaled_gradients: bool = False,
5660
) -> Dict[str, Any]:
5761
"""Simulate PEtab model.
5862
63+
.. note::
64+
Regardless of `scaled_parameters`, unscaled sensitivities are returned,
65+
unless `scaled_gradients=True`.
66+
5967
:param petab_problem:
6068
PEtab problem to work on.
6169
:param amici_model:
@@ -87,6 +95,9 @@ def simulate_petab(
8795
:param failfast:
8896
Returns as soon as an integration failure is encountered, skipping
8997
any remaining simulations.
98+
:param scaled_gradients:
99+
Whether to compute gradients on parameter scale (``True``) or not
100+
(``False``).
90101
91102
:return:
92103
Dictionary of
@@ -104,15 +115,13 @@ def simulate_petab(
104115
if solver is None:
105116
solver = amici_model.getSolver()
106117

107-
# Get parameters
108-
if problem_parameters is None:
109-
# Use PEtab nominal values as default
110-
problem_parameters = {t.Index: getattr(t, NOMINAL_VALUE) for t in
111-
petab_problem.parameter_df.itertuples()}
112-
if scaled_parameters:
113-
raise NotImplementedError(
114-
"scaled_parameters=True in combination with "
115-
"problem_parameters=None is currently not supported.")
118+
# Switch to scaled parameters.
119+
problem_parameters = _default_scaled_parameters(
120+
petab_problem=petab_problem,
121+
problem_parameters=problem_parameters,
122+
scaled_parameters=scaled_parameters,
123+
)
124+
scaled_parameters = True
116125

117126
# number of amici simulations will be number of unique
118127
# (preequilibrationConditionId, simulationConditionId) pairs.
@@ -153,6 +162,29 @@ def simulate_petab(
153162

154163
# Compute total llh
155164
llh = sum(rdata['llh'] for rdata in rdatas)
165+
# Compute total sllh
166+
sllh = None
167+
if solver.getSensitivityOrder() != amici.SensitivityOrder.none:
168+
sllh = aggregate_sllh(
169+
amici_model=amici_model,
170+
rdatas=rdatas,
171+
parameter_mapping=parameter_mapping,
172+
petab_scale=scaled_parameters,
173+
petab_problem=petab_problem,
174+
edatas=edatas,
175+
)
176+
if not scaled_gradients and sllh is not None:
177+
sllh = {
178+
parameter_id: rescale_sensitivity(
179+
sensitivity=sensitivity,
180+
parameter_value=problem_parameters[parameter_id],
181+
old_scale=petab_problem.parameter_df.loc[
182+
parameter_id, PARAMETER_SCALE
183+
],
184+
new_scale=LIN,
185+
)
186+
for parameter_id, sensitivity in sllh.items()
187+
}
156188

157189
# Log results
158190
sim_cond = petab_problem.get_simulation_conditions_from_measurement_df()
@@ -165,11 +197,158 @@ def simulate_petab(
165197

166198
return {
167199
LLH: llh,
200+
SLLH: sllh,
168201
RDATAS: rdatas,
169202
EDATAS: edatas,
170203
}
171204

172205

206+
def aggregate_sllh(
207+
amici_model: AmiciModel,
208+
rdatas: Sequence[amici.ReturnDataView],
209+
parameter_mapping: Optional[ParameterMapping],
210+
edatas: List[AmiciExpData],
211+
petab_scale: bool = True,
212+
petab_problem: petab.Problem = None,
213+
) -> Union[None, Dict[str, float]]:
214+
"""
215+
Aggregate likelihood gradient for all conditions, according to PEtab
216+
parameter mapping.
217+
218+
:param amici_model:
219+
AMICI model from which ``rdatas`` were obtained.
220+
:param rdatas:
221+
Simulation results.
222+
:param parameter_mapping:
223+
PEtab parameter mapping to condition-specific simulation parameters.
224+
:param edatas:
225+
Experimental data used for simulation.
226+
:param petab_scale:
227+
Whether to check that sensitivities were computed with parameters on
228+
the scales provided in the PEtab parameters table.
229+
:param petab_problem:
230+
The PEtab problem that defines the parameter scales.
231+
232+
:return:
233+
Aggregated likelihood sensitivities.
234+
"""
235+
accumulated_sllh = {}
236+
model_parameter_ids = amici_model.getParameterIds()
237+
238+
if petab_scale and petab_problem is None:
239+
raise ValueError(
240+
'Please provide the PEtab problem, when using '
241+
'`petab_scale=True`.'
242+
)
243+
244+
# Check for issues in all condition simulation results.
245+
for rdata in rdatas:
246+
# Condition failed during simulation.
247+
if rdata.status != amici.AMICI_SUCCESS:
248+
return None
249+
# Condition simulation result does not provide SLLH.
250+
if rdata.sllh is None:
251+
raise ValueError(
252+
'The sensitivities of the likelihood for a condition were '
253+
'not computed.'
254+
)
255+
256+
for condition_parameter_mapping, edata, rdata in \
257+
zip(parameter_mapping, edatas, rdatas):
258+
for sllh_parameter_index, condition_parameter_sllh in \
259+
enumerate(rdata.sllh):
260+
# Get PEtab parameter ID
261+
# Use ExpData if it provides a parameter list, else default to
262+
# Model.
263+
if edata.plist:
264+
model_parameter_index = edata.plist[sllh_parameter_index]
265+
else:
266+
model_parameter_index = amici_model.plist(sllh_parameter_index)
267+
model_parameter_id = model_parameter_ids[model_parameter_index]
268+
petab_parameter_id = (
269+
condition_parameter_mapping
270+
.map_sim_var[model_parameter_id]
271+
)
272+
273+
# Initialize
274+
if petab_parameter_id not in accumulated_sllh:
275+
accumulated_sllh[petab_parameter_id] = 0
276+
277+
# Check that the scale is consistent
278+
if petab_scale:
279+
# `ParameterMappingForCondition` objects provide the scale in
280+
# terms of `petab.C` constants already, not AMICI equivalents.
281+
model_parameter_scale = (
282+
condition_parameter_mapping
283+
.scale_map_sim_var[model_parameter_id]
284+
)
285+
petab_parameter_scale = (
286+
petab_problem
287+
.parameter_df
288+
.loc[petab_parameter_id, PARAMETER_SCALE]
289+
)
290+
if model_parameter_scale != petab_parameter_scale:
291+
raise ValueError(
292+
f'The scale of the parameter `{petab_parameter_id}` '
293+
'differs between the AMICI model '
294+
f'({model_parameter_scale}) and the PEtab problem '
295+
f'({petab_parameter_scale}).'
296+
)
297+
298+
# Accumulate
299+
accumulated_sllh[petab_parameter_id] += condition_parameter_sllh
300+
301+
return accumulated_sllh
302+
303+
304+
def rescale_sensitivity(
305+
sensitivity: float,
306+
parameter_value: float,
307+
old_scale: str,
308+
new_scale: str,
309+
) -> float:
310+
"""Rescale a sensitivity between parameter scales.
311+
312+
:param sensitivity:
313+
The sensitivity corresponding to the parameter value.
314+
:param parameter_value:
315+
The parameter vector element, on ``old_scale``.
316+
:param old_scale:
317+
The scale of the parameter value.
318+
:param new_scale:
319+
The parameter scale on which to rescale the sensitivity.
320+
321+
:return:
322+
The rescaled sensitivity.
323+
"""
324+
LOG_E_10 = np.log(10)
325+
326+
if old_scale == new_scale:
327+
return sensitivity
328+
329+
unscaled_parameter_value = petab.parameters.unscale(
330+
parameter=parameter_value,
331+
scale_str=old_scale,
332+
)
333+
334+
scale = {
335+
(LIN, LOG): lambda s: s * unscaled_parameter_value,
336+
(LOG, LIN): lambda s: s / unscaled_parameter_value,
337+
(LIN, LOG10): lambda s: s * (unscaled_parameter_value * LOG_E_10),
338+
(LOG10, LIN): lambda s: s / (unscaled_parameter_value * LOG_E_10),
339+
}
340+
341+
scale[(LOG, LOG10)] = lambda s: scale[(LIN, LOG10)](scale[(LOG, LIN)](s))
342+
scale[(LOG10, LOG)] = lambda s: scale[(LIN, LOG)](scale[(LOG10, LIN)](s))
343+
344+
if (old_scale, new_scale) not in scale:
345+
raise NotImplementedError(
346+
f"Old scale: {old_scale}. New scale: {new_scale}."
347+
)
348+
349+
return scale[(old_scale, new_scale)](sensitivity)
350+
351+
173352
def create_parameterized_edatas(
174353
amici_model: AmiciModel,
175354
petab_problem: petab.Problem,
@@ -811,3 +990,39 @@ def rdatas_to_simulation_df(
811990
measurement_df=measurement_df)
812991

813992
return df.rename(columns={MEASUREMENT: SIMULATION})
993+
994+
995+
def _default_scaled_parameters(
996+
petab_problem: petab.Problem,
997+
problem_parameters: Optional[Dict[str, float]] = None,
998+
scaled_parameters: bool = False,
999+
) -> Optional[Dict[str, float]]:
1000+
"""
1001+
Helper method to handle an unscaled or unspecified parameter vector.
1002+
1003+
The parameter vector defaults to the nominal values in the PEtab
1004+
parameter table.
1005+
1006+
Unscaled parameter values are scaled.
1007+
1008+
:param petab_problem:
1009+
The PEtab problem.
1010+
:param problem_parameters:
1011+
Keys are PEtab parameter IDs, values are parameter values on the scale
1012+
defined in the PEtab parameter table. Defaults to the nominal values in
1013+
the PEtab parameter table.
1014+
:param scaled_parameters:
1015+
Whether `problem_parameters` are on the scale defined in the PEtab
1016+
parameter table.
1017+
1018+
:return:
1019+
The scaled parameter vector.
1020+
"""
1021+
if problem_parameters is None:
1022+
problem_parameters = dict(zip(
1023+
petab_problem.x_ids,
1024+
petab_problem.x_nominal_scaled,
1025+
))
1026+
elif not scaled_parameters:
1027+
problem_parameters = petab_problem.scale_parameters(problem_parameters)
1028+
return problem_parameters

‎python/tests/petab_test_problems/lotka_volterra/README.md‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
A Lotka-Volterra model, based on the model provided as an example in the `yaml2sbml` package: https://github.com/yaml2sbml-dev/yaml2sbml
2+
3+
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/
4+
5+
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‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
#!/usr/bin/env bash
2+
set -eou pipefail
3+
script_dir=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )
4+
petab_path=$script_dir/../petab
5+
6+
mkdir -p "${petab_path}"
7+
8+
# Validate input yaml2sbml model
9+
yaml2sbml_validate lotka_volterra.yaml
10+
11+
# Copy measurements to PEtab directory
12+
cp "${script_dir}/measurements.tsv" "$petab_path/measurements.tsv"
13+
14+
# Write the PEtab problem
15+
python "${script_dir}/writer.py"
16+
17+
# Remove condition parameters from PEtab parameters table
18+
for condition_parameter in beta delta departure_prey arrival_predator
19+
do
20+
sed -i "/^${condition_parameter}/d" "${petab_path}/parameter*"
21+
done
22+
23+
# Validate the PEtab problem
24+
petablint -vy "${petab_path}/problem.yaml"

0 commit comments

Comments
 (0)