..

Author: dweindl

doc/modules.rst

@@ -32,6 +32,7 @@ API Reference
    petab.v1.yaml
    petab.v2
    petab.v2.C
+   petab.v2.converters
    petab.v2.core
    petab.v2.experiments
    petab.v2.lint

petab/v2/converters.py

@@ -2,44 +2,45 @@
 
 from __future__ import annotations
 
+import warnings
 from copy import deepcopy
 from math import inf
 
 import libsbml
-from sbmlmath import set_math
+from sbmlmath import sbml_math_to_sympy, set_math
 
 from .core import Change, Condition, Experiment, ExperimentPeriod
 from .models._sbml_utils import add_sbml_parameter, check
 from .models.sbml_model import SbmlModel
 from .problem import Problem
 
+__all__ = ["ExperimentsToEventsConverter"]
+
 
 class ExperimentsToEventsConverter:
     """Convert PEtab experiments to SBML events.
 
     For an SBML-model-based PEtab problem, this class converts the PEtab
     experiments to events as far as possible.
 
-    Currently, this assumes that there is no other event in the model
-    that could trigger at the same time as the events created here.
-    I.e., the events responsible for applying PEtab condition changes
-    don't have a priority assigned that would guarantee that they are executed
-    before any pre-existing events.
+    If the model already contains events, PEtab events are added with a higher
+    priority than the existing events to guarantee that PEtab condition changes
+    are applied before any pre-existing assignments.
 
     The PEtab problem must not contain any identifiers starting with
     ``_petab``.
 
     All periods and condition changes that are represented by events
     will be removed from the condition table.
-    Each experiment will have at most one period with a start time of -inf
+    Each experiment will have at most one period with a start time of ``-inf``
     and one period with a finite start time. The associated changes with
     these periods are only the steady-state pre-simulation indicator
     (if necessary), and the experiment indicator parameter.
     """
 
     #: ID of the parameter that indicates whether the model is in
-    #   the steady-state pre-simulation phase (1) or not (0).
-    PRE_STEADY_STATE_INDICATOR = "_petab_pre_steady_state_indicator"
+    #  the steady-state pre-simulation phase (1) or not (0).
+    PRE_SIM_INDICATOR = "_petab_pre_simulation_indicator"
 
     def __init__(self, problem: Problem):
         """Initialize the converter.
@@ -50,30 +51,81 @@ def __init__(self, problem: Problem):
         if not isinstance(problem.model, SbmlModel):
             raise ValueError("Only SBML models are supported.")
 
-        self._problem = problem
-        self._model = problem.model.sbml_model
-        self._presim_indicator = self.PRE_STEADY_STATE_INDICATOR
+        self._original_problem = problem
+        self._new_problem = deepcopy(self._original_problem)
+
+        self._model = self._new_problem.model.sbml_model
+        self._presim_indicator = self.PRE_SIM_INDICATOR
+
+        # The maximum event priority that was found in the unprocessed model.
+        self._max_event_priority = None
+        # The priority that will be used for the PEtab events.
+        self._petab_event_priority = None
+
+        self._preprocess()
 
+    def _preprocess(self):
+        """Check whether we can handle the given problem and store some model
+        information."""
         model = self._model
         if model.getLevel() < 3:
+            # try to upgrade the SBML model
             if not model.getSBMLDocument().setLevelAndVersion(3, 2):
                 raise ValueError(
                     "Cannot handle SBML models with SBML level < 3, "
                     "because they do not support initial values for event "
                     "triggers and automatic upconversion failed."
                 )
 
+        # Collect event priorities
+        event_priorities = {
+            ev.getId() or str(ev): sbml_math_to_sympy(ev.getPriority())
+            for ev in model.getListOfEvents()
+            if ev.getPriority() and ev.getPriority().getMath() is not None
+        }
+
+        # Check for non-constant event priorities and track the maximum
+        #  priority used so far.
+        for e, priority in event_priorities.items():
+            if priority.free_symbols:
+                # We'd need to find the maximum priority of all events,
+                #  which is challenging/impossible to do in general.
+                raise NotImplementedError(
+                    f"Event `{e}` has a non-constant priority: {priority}. "
+                    "This is currently not supported."
+                )
+            self._max_event_priority = max(
+                self._max_event_priority or 0, float(priority)
+            )
+
+        self._petab_event_priority = (
+            self._max_event_priority + 1
+            if self._max_event_priority is not None
+            else None
+        )
+        # Check for undefined event priorities and warn
+        for event in model.getListOfEvents():
+            if (prio := event.getPriority()) and prio.getMath() is None:
+                warnings.warn(
+                    f"Event `{event.getId()}` has no priority set. "
+                    "Make sure that this event cannot trigger at the time of "
+                    "PEtab condition change, otherwise the behavior is "
+                    "undefined.",
+                    stacklevel=1,
+                )
+
     def convert(self) -> Problem:
         """Convert the PEtab experiments to SBML events.
 
         :return: The converted PEtab problem.
         """
-        problem = deepcopy(self._problem)
 
         self._add_presimulation_indicator()
 
+        problem = self._new_problem
         for experiment in problem.experiment_table.experiments:
             self._convert_experiment(problem, experiment)
+
         self._add_indicators_to_conditions(problem)
 
         validation_results = problem.validate()
@@ -99,9 +151,7 @@ def _convert_experiment(self, problem: Problem, experiment: Experiment):
         kept_periods = []
         for i_period, period in enumerate(experiment.periods):
             # check for non-zero initial times of the first period
-            if (
-                i_period == 0 or (i_period == 1 and has_presimulation)
-            ) and period.time != 0:
+            if (i_period == int(has_presimulation)) and period.time != 0:
                 # TODO: we could address that by offsetting all occurrences of
                 #  the SBML time in the model (except for the newly added
                 #  events triggers). Or we better just leave it to the
@@ -161,17 +211,15 @@ def _create_period_begin_event(
         # TODO: for now, add separate events for each experiment x period,
         #  this could be optimized to reuse events
 
-        # TODO if there is already some event that could trigger
-        #  at this time, we need event priorities. This is difficult to
-        #  implement, though, since in general, we can't know the maximum
-        #  priority of the other events, unless they are static.
-
         ev = self._model.createEvent()
         check(ev.setId(f"_petab_event_{experiment.id}_{i_period}"))
         check(ev.setUseValuesFromTriggerTime(True))
         trigger = ev.createTrigger()
         check(trigger.setInitialValue(False))  # may trigger at t=0
         check(trigger.setPersistent(True))
+        if self._petab_event_priority is not None:
+            priority = ev.createPriority()
+            set_math(priority, self._petab_event_priority)
 
         exp_ind_id = self.get_experiment_indicator(experiment.id)