Skip to content

Commit a514bfd

Browse files
dweindldweindl
committed
v2: Add experiments->events converter
Add functionality to convert PEtab v2 experiments/conditions to SBML events. This should make it easier to implement v2 support in other tools.
1 parent ef79523 commit a514bfd

File tree

3 files changed

+310
-0
lines changed

3 files changed

+310
-0
lines changed

‎petab/v1/calculate.py‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -408,6 +408,17 @@ def calculate_single_llh(
408408
Returns:
409409
The computed likelihood for the given values.
410410
"""
411+
412+
# PEtab v2:
413+
from petab.v2.C import LOG10_NORMAL, LOG_NORMAL
414+
415+
if noise_distribution == LOG10_NORMAL and scale == LIN:
416+
noise_distribution = NORMAL
417+
scale = LOG10
418+
elif noise_distribution == LOG_NORMAL and scale == LIN:
419+
noise_distribution = NORMAL
420+
scale = LOG
421+
411422
# short-hand
412423
m, s, sigma = measurement, simulation, noise_value
413424
pi, log, log10 = np.pi, np.log, np.log10

‎petab/v2/converters.py‎

Lines changed: 289 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,289 @@
1+
"""Conversion of PEtab problems."""
2+
3+
from __future__ import annotations
4+
5+
from copy import deepcopy
6+
from math import inf
7+
from typing import TYPE_CHECKING
8+
9+
from .core import Change, Condition
10+
from .problem import Problem
11+
12+
if TYPE_CHECKING:
13+
import libsbml
14+
15+
16+
PRE_STEADY_STATE_INDICATOR = "_petab_pre_steady_state_indicator"
17+
18+
# TODO to class?
19+
20+
21+
def get_experiment_indicator(experiment_id: str) -> str:
22+
"""The ID of the experiment indicator parameter.
23+
24+
The experiment indicator parameter is used to identify the
25+
experiment in the SBML model. It is a parameter that is set
26+
to 1 for the current experiment and 0 for all other
27+
experiments. The parameter is used in the event trigger
28+
to determine when the event should be triggered.
29+
30+
:param experiment_id: The ID of the experiment for which to create
31+
the experiment indicator parameter ID.
32+
"""
33+
return f"_petab_experiment_indicator_{experiment_id}"
34+
35+
36+
def experiments_to_events(problem: Problem) -> Problem:
37+
"""
38+
Convert experiments to events.
39+
40+
For an SBML-model-based PEtab problem, replace the PEtab experiments
41+
by events as far as possible.
42+
43+
Currently, this assumes that there is no other event in the model
44+
that could trigger at the same time as the events created here.
45+
I.e., the events responsible for applying PEtab condition changes
46+
don't have a priority assigned that would guarantee that they are executed
47+
before any pre-existing events.
48+
49+
The PEtab problem must not contain any identifiers starting with
50+
``_petab``.
51+
52+
:param problem: The PEtab problem to convert. This will not be modified.
53+
:return: The converted PEtab problem.
54+
All periods and condition changes that are represented by events
55+
are removed from the condition table.
56+
Each experiment will have at most one period with a start time of -inf
57+
and one period with a finite start time. The associated changes with
58+
these periods are only the steady-state pre-simulation indicator
59+
(if necessary), and the experiment indicator parameter.
60+
"""
61+
import libsbml
62+
63+
from .models.sbml_model import SbmlModel
64+
65+
if not isinstance(problem.model, SbmlModel):
66+
raise ValueError("Only SBML models are supported.")
67+
68+
# TODO move elsewhere
69+
retval_to_str = {
70+
getattr(libsbml, attr): attr
71+
for attr in (
72+
"LIBSBML_DUPLICATE_OBJECT_ID",
73+
"LIBSBML_INDEX_EXCEEDS_SIZE",
74+
"LIBSBML_INVALID_ATTRIBUTE_VALUE",
75+
"LIBSBML_INVALID_OBJECT",
76+
"LIBSBML_INVALID_XML_OPERATION",
77+
"LIBSBML_LEVEL_MISMATCH",
78+
"LIBSBML_NAMESPACES_MISMATCH",
79+
"LIBSBML_OPERATION_FAILED",
80+
"LIBSBML_UNEXPECTED_ATTRIBUTE",
81+
"LIBSBML_PKG_UNKNOWN",
82+
"LIBSBML_PKG_VERSION_MISMATCH",
83+
"LIBSBML_PKG_CONFLICTED_VERSION",
84+
)
85+
}
86+
87+
def _check(res):
88+
if res != libsbml.LIBSBML_OPERATION_SUCCESS:
89+
raise RuntimeError(f"libsbml error: {retval_to_str.get(res, res)}")
90+
91+
problem = deepcopy(problem)
92+
sbml_mod = problem.model.sbml_model
93+
94+
# add pre-steady-state indicator parameter
95+
pre_steady_state_ind_id = PRE_STEADY_STATE_INDICATOR
96+
if sbml_mod.getElementBySId(pre_steady_state_ind_id) is not None:
97+
raise AssertionError(
98+
f"Entity with ID {pre_steady_state_ind_id} exists already."
99+
)
100+
pre_steady_state_ind = sbml_mod.createParameter()
101+
pre_steady_state_ind.setId(pre_steady_state_ind_id)
102+
pre_steady_state_ind.setValue(0)
103+
pre_steady_state_ind.setConstant(False)
104+
105+
for experiment in problem.experiment_table.experiments:
106+
experiment.sort_periods()
107+
108+
# add experiment indicator
109+
exp_ind_id = get_experiment_indicator(experiment.id)
110+
if sbml_mod.getElementBySId(exp_ind_id) is not None:
111+
raise AssertionError(
112+
f"Entity with ID {exp_ind_id} exists already."
113+
)
114+
_add_sbml_parameter(sbml_mod, id_=exp_ind_id, constant=False, value=0)
115+
116+
kept_periods = []
117+
for i_period, period in enumerate(experiment.periods):
118+
if (
119+
i_period == 0
120+
or (i_period == 1 and experiment.periods[0].time == -inf)
121+
) and period.time != 0:
122+
raise NotImplementedError(
123+
"Cannot represent non-zero initial time in SBML."
124+
)
125+
126+
if is_preeq_period := (period.time == -inf):
127+
# steady-state pre-simulation cannot be represented in SBML,
128+
# so we need to keep this period in the Problem.
129+
kept_periods.append(period)
130+
elif not period.condition_ids:
131+
# no condition, no changes, no need for an event,
132+
# no need to keep the period unless it's the initial
133+
# steady-state simulation
134+
continue
135+
elif period.time == 0:
136+
kept_periods.append(period)
137+
138+
if sbml_mod.getLevel() < 3:
139+
if not sbml_mod.getSBMLDocument().setLevelAndVersion(3, 2):
140+
raise ValueError(
141+
"Cannot handle SBML models with SBML level < 3, "
142+
"because they do not support initial values for event "
143+
"triggers and automatic upconversion failed."
144+
)
145+
146+
# TODO: for now, add separate events for each experiment x period,
147+
# this could be optimized to reuse events
148+
149+
# TODO if there is already some event that could trigger
150+
# at this time, we need event priorities. This is difficult to
151+
# implement, though, since in general, we can't know the maximum
152+
# priority of the other events, unless they are static.
153+
ev = sbml_mod.createEvent()
154+
_check(ev.setId(f"_petab_event_{experiment.id}_{i_period}"))
155+
_check(ev.setUseValuesFromTriggerTime(True))
156+
trigger = ev.createTrigger()
157+
_check(trigger.setInitialValue(False)) # may trigger at t=0
158+
_check(trigger.setPersistent(True))
159+
160+
if is_preeq_period:
161+
trig_math = libsbml.parseL3Formula(
162+
f"({exp_ind_id} == 1) && ({pre_steady_state_ind_id} == 1)"
163+
)
164+
else:
165+
trig_math = libsbml.parseL3Formula(
166+
f"({exp_ind_id} == 1) && ({pre_steady_state_ind_id} != 1) "
167+
f"&& (time >= {period.time})"
168+
)
169+
_check(trigger.setMath(trig_math))
170+
171+
_create_event_assignments_for_period(
172+
ev,
173+
[
174+
problem.condition_table[condition_id]
175+
for condition_id in period.condition_ids
176+
],
177+
)
178+
179+
if len(kept_periods) > 2:
180+
raise AssertionError("Expected at most two periods to be kept.")
181+
182+
# add conditions that set the indicator parameters
183+
for period in kept_periods:
184+
period.condition_ids = [
185+
f"_petab_experiment_condition_{experiment.id}",
186+
"_petab_steady_state_pre_simulation"
187+
if period.time == -inf
188+
else "_petab_no_steady_state_pre_simulation",
189+
]
190+
experiment.periods = kept_periods
191+
192+
# create conditions for indicator parameters
193+
problem.condition_table.conditions.append(
194+
Condition(
195+
id="_petab_steady_state_pre_simulation",
196+
changes=[
197+
Change(target_id=pre_steady_state_ind_id, target_value=1)
198+
],
199+
)
200+
)
201+
problem.condition_table.conditions.append(
202+
Condition(
203+
id="_petab_no_steady_state_pre_simulation",
204+
changes=[
205+
Change(target_id=pre_steady_state_ind_id, target_value=0)
206+
],
207+
)
208+
)
209+
for experiment in problem.experiment_table.experiments:
210+
problem.condition_table.conditions.append(
211+
Condition(
212+
id=f"_petab_experiment_condition_{experiment.id}",
213+
changes=[
214+
Change(
215+
target_id=get_experiment_indicator(experiment.id),
216+
target_value=1,
217+
)
218+
],
219+
)
220+
)
221+
222+
# all changes have been encoded in event assignments and can be
223+
# removed
224+
problem.condition_table.conditions = [
225+
condition
226+
for condition in problem.condition_table.conditions
227+
if condition.id.startswith("_petab")
228+
]
229+
230+
validation_results = problem.validate()
231+
validation_results.log()
232+
233+
return problem
234+
235+
236+
def _create_event_assignments_for_period(
237+
event: libsbml.Event, conditions: list[Condition]
238+
) -> None:
239+
"""Create an event assignment for a given period."""
240+
from sbmlmath import set_math
241+
242+
sbml_model = event.getModel()
243+
for condition in conditions:
244+
for change in condition.changes:
245+
ea = event.createEventAssignment()
246+
ea.setVariable(change.target_id)
247+
set_math(ea, change.target_value)
248+
249+
# target needs const=False, and target may not exist yet
250+
# (e.g., in case of output parameters added in the observable
251+
# table)
252+
target = sbml_model.getElementBySId(change.target_id)
253+
if target is None:
254+
_add_sbml_parameter(
255+
sbml_model, id_=change.target_id, constant=False, value=0
256+
)
257+
else:
258+
# TODO: can that break models??
259+
target.setConstant(False)
260+
261+
# the target value may depend on parameters that are only
262+
# introduced in the PEtab parameter table - those need
263+
# to be added to the model
264+
for sym in change.target_value.free_symbols:
265+
if sbml_model.getElementBySId(sym.name) is None:
266+
_add_sbml_parameter(
267+
sbml_model, id_=sym.name, constant=True, value=0
268+
)
269+
270+
271+
def _add_sbml_parameter(
272+
model: libsbml.Model,
273+
id_: str = None,
274+
value: float = None,
275+
constant: bool = None,
276+
) -> libsbml.Parameter:
277+
"""Add a parameter to the SBML model."""
278+
param = model.createParameter()
279+
280+
if id_ is not None:
281+
param.setId(id_)
282+
283+
if value is not None:
284+
param.setValue(value)
285+
286+
if constant is not None:
287+
param.setConstant(constant)
288+
289+
return param

‎petab/v2/core.py‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -114,6 +114,8 @@ class NoiseDistribution(str, Enum):
114114
LOG_NORMAL = C.LOG_NORMAL
115115
#: Log-Laplace distribution
116116
LOG_LAPLACE = C.LOG_LAPLACE
117+
#: Log10-Normal
118+
LOG10_NORMAL = C.LOG10_NORMAL
117119

118120

119121
class PriorDistribution(str, Enum):
@@ -553,6 +555,14 @@ def __iadd__(self, other: ExperimentPeriod) -> Experiment:
553555
self.periods.append(other)
554556
return self
555557

558+
def has_steady_state_presimulation(self) -> bool:
559+
"""Check if the experiment has a steady-state pre-simulation."""
560+
return any(period.time == -np.inf for period in self.periods)
561+
562+
def sort_periods(self) -> None:
563+
"""Sort the periods of the experiment by time."""
564+
self.periods.sort(key=lambda period: period.time)
565+
556566

557567
class ExperimentTable(BaseModel):
558568
"""PEtab experiments table."""

0 commit comments

Comments
 (0)