Auto-assign event priority to a single unprioritized event

CHANGELOG.md

@@ -54,6 +54,11 @@
 * `assignments_to_experiments` no longer creates condition IDs that clash
   with other IDs in the problem
   (by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/561)
+* `petab.v2.converters.ExperimentsToSbmlConverter` now assigns a priority to
+  a pre-existing event without priority if it is the only such event, so
+  that PEtab condition changes are applied before its assignments.
+  Previously, their order was undefined
+  (by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/562)
 
 **Features**
 * `petab1to2` and `petab_files_1to2` can skip validation via `validate=False`

petab/v2/converters.py

@@ -2,6 +2,8 @@
 
 from __future__ import annotations
 
+import logging
+import math
 import warnings
 from copy import deepcopy
 
@@ -22,18 +24,26 @@
 
 __all__ = ["ExperimentsToSbmlConverter", "assignments_to_experiments"]
 
+logger = logging.getLogger(__name__)
+
 
 class ExperimentsToSbmlConverter:
     """Convert PEtab experiments to SBML.
 
     For an SBML-model-based PEtab problem, this class converts the PEtab
     experiments to initial assignments and events as far as possible.
 
-    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.
+    If the model already contains events, the events implementing the PEtab
+    condition changes are added with a higher priority than the pre-existing
+    events, so that these changes are applied before any pre-existing
+    assignments.
     This requires that all event priorities in the original model are numeric
     constants.
+    If exactly one pre-existing event has no priority, it is assigned a
+    priority lower than all other priorities in the model.
+    If multiple events have no priority, their order relative to the
+    condition-change events is undefined, unless ``default_priority`` is
+    provided.
 
     The PEtab problem must not contain any identifiers starting with
     ``_petab``.
@@ -71,10 +81,15 @@ def __init__(
             In SBML, for event assignments that are to be applied at the same
             simulation time, the order of event execution is determined by the
             priority of the respective events.
-            If no priority is set, the order is undefined.
+            For events without priority, the order is undefined.
             See SBML specs for details.
-            To ensure that the PEtab condition-start-events are executed before
+            To ensure that the condition-change events are executed before
             any other events, all events should have a priority set.
+            Note that, according to the SBML specs, events with equal
+            priorities that are executed simultaneously must be executed in
+            random order.
+            If ``None``, a priority is only assigned automatically if exactly
+            one pre-existing event has no priority.
         """
         if len(problem.models) > 1:
             #  https://github.com/PEtab-dev/libpetab-python/issues/392
@@ -97,7 +112,7 @@ def __init__(
 
         # 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.
+        # The priority that will be used for the condition-change events.
         self._petab_event_priority = None
         self._default_priority = default_priority
         self._preprocess()
@@ -125,35 +140,30 @@ def _preprocess(self) -> None:
         # Apply default priority to all events that do not have a priority
         if self._default_priority is not None:
             for event in model.getListOfEvents():
-                if (
-                    not event.getPriority()
-                    or event.getPriority().getMath() is None
-                ):
+                if not self._has_priority(event):
                     priority = event.createPriority()
                     priority.setMath(
                         libsbml.parseL3Formula(str(self._default_priority))
                     )
+        else:
+            self._auto_assign_priority(self._get_event_priorities())
 
-        # 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."
+        # Track the maximum priority used so far
+        event_priorities = self._get_event_priorities()
+        for event_id, priority in event_priorities.items():
+            if math.isnan(priority) or priority == math.inf:
+                warnings.warn(
+                    f"Event `{event_id}` has priority {priority}. "
+                    "The condition-change events cannot be assigned a "
+                    "higher priority. Make sure that this event cannot "
+                    "trigger at the time of a PEtab condition change, "
+                    "otherwise simulation results may be incorrect.",
+                    stacklevel=1,
                 )
-            self._max_event_priority = max(
-                self._max_event_priority or 0, float(priority)
-            )
+        if non_nan := [
+            p for p in event_priorities.values() if not math.isnan(p)
+        ]:
+            self._max_event_priority = max(0, *non_nan)
 
         self._petab_event_priority = (
             self._max_event_priority + 1
@@ -190,6 +200,80 @@ def _preprocess(self) -> None:
                     stacklevel=1,
                 )
 
+    @staticmethod
+    def _has_priority(event: libsbml.Event) -> bool:
+        """Check whether the given event has a priority expression."""
+        return (
+            event.isSetPriority() and event.getPriority().getMath() is not None
+        )
+
+    def _get_event_priorities(self) -> dict[str, float]:
+        """Get the priorities of all events that have a priority.
+
+        :return: Mapping of event IDs to priorities.
+        :raises NotImplementedError: If any priority is not constant.
+        """
+        event_priorities = {}
+        for ev in self._model.getListOfEvents():
+            if not self._has_priority(ev):
+                continue
+            event_id = ev.getId() or str(ev)
+            priority = sbml_math_to_sympy(ev.getPriority())
+            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 `{event_id}` has a non-constant priority: "
+                    f"{priority}. This is currently not supported."
+                )
+            event_priorities[event_id] = float(priority)
+        return event_priorities
+
+    def _auto_assign_priority(
+        self, event_priorities: dict[str, float]
+    ) -> None:
+        """Assign a priority to a single unprioritized pre-existing event.
+
+        If exactly one pre-existing event has no priority, and no other
+        event has priority ``-inf`` or ``NaN``, the unprioritized event is
+        assigned a priority lower than any other finite priority in the
+        model.
+        Without priority, the order of this event with respect to all other
+        events is undefined, so assigning a priority only selects one of the
+        admissible orders. A distinct priority cannot create ties, which
+        would require random execution order (SBML L3V2 §4.12.3). This allows
+        executing the condition-change events before this event.
+
+        If more than one event has no priority, no priority is assigned,
+        since this would impose an arbitrary order on them if they are
+        executed simultaneously.
+
+        :param event_priorities: The priorities of all prioritized events.
+        """
+        unprioritized = [
+            ev
+            for ev in self._model.getListOfEvents()
+            if not self._has_priority(ev)
+        ]
+        if len(unprioritized) != 1:
+            return
+
+        priorities = event_priorities.values()
+        if any(math.isnan(p) or p == -math.inf for p in priorities):
+            # no distinct lower priority available
+            return
+
+        value = min((p for p in priorities if p != math.inf), default=1) - 1
+        event = unprioritized[0]
+        priority = event.createPriority()
+        set_math(priority, value)
+        logger.debug(
+            "Assigned priority %s to event `%s` to ensure that the "
+            "condition-change events are executed before this event.",
+            value,
+            event.getId(),
+        )
+
     def convert(self) -> Problem:
         """Convert the PEtab experiments to SBML events.
 

tests/v2/test_converters.py

@@ -1,6 +1,7 @@
 from math import inf
 
 import pandas as pd
+import pytest
 import sympy as sp
 from sbmlmath import sbml_math_to_sympy
 
@@ -286,3 +287,124 @@ def test_simulate_experiment_to_events():
 
     for col in df_expected.columns:
         assert (df_expected[col] == df_actual[col]).all()
+
+
+def _event_priorities(sbml_model) -> dict[str, float | None]:
+    """Get the numeric priorities of all events in the given model."""
+    return {
+        ev.getId(): float(sbml_math_to_sympy(ev.getPriority()))
+        if ev.isSetPriority() and ev.getPriority().getMath() is not None
+        else None
+        for ev in sbml_model.getListOfEvents()
+    }
+
+
+def _problem_0030(model_events: str) -> Problem:
+    """PEtab test suite case 0030 with the given model events."""
+    ant_model = f"""
+    compartment C = 4
+    species S in C = 3
+    p = 1
+    S' = p
+    {model_events}
+    """
+    problem = Problem()
+    problem.model = SbmlModel.from_antimony(ant_model)
+    problem.add_condition("condition1", S=2)
+    problem.add_condition("condition2", S="S + C", C=8)
+    problem.add_experiment("experiment1", 0, "condition1", 10, "condition2")
+    return problem
+
+
+E0 = "E0: at time >= 10, fromTrigger=false: C = C * 2"
+E0_PRIO = "E0: at time >= 10, priority=-2, fromTrigger=false: C = C * 2"
+E1 = "E1: at time >= 20, fromTrigger=false: C = C"
+E1_PRIO = "E1: at time >= 20, priority=5, fromTrigger=false: C = C"
+E1_PRIO_NEG_INF = "E1: at time >= 20, priority=-inf, fromTrigger=false: C = C"
+E1_PRIO_INF = "E1: at time >= 20, priority=inf, fromTrigger=false: C = C"
+E2_PRIO = "E2: at time >= 30, priority=5, fromTrigger=false: C = C"
+PETAB_EVENT = "_petab_event_experiment1_1"
+
+
+@pytest.mark.parametrize(
+    "model_events,expected",
+    [
+        # single unprioritized event
+        (E0, {"E0": 0, PETAB_EVENT: 1}),
+        # one unprioritized event, others prioritized
+        (f"{E0}\n{E1_PRIO}", {"E0": 4, "E1": 5, PETAB_EVENT: 6}),
+        # no distinct lower priority available: no automatic assignment for
+        #  the unprioritized event
+        (
+            f"{E0}\n{E1_PRIO_NEG_INF}",
+            {"E0": None, "E1": -inf, PETAB_EVENT: 1},
+        ),
+        # all events prioritized: unchanged
+        (f"{E0_PRIO}\n{E1_PRIO}", {"E0": -2, "E1": 5, PETAB_EVENT: 6}),
+        # multiple unprioritized events: no automatic assignment
+        (f"{E0}\n{E1}", {"E0": None, "E1": None, PETAB_EVENT: None}),
+        # no model events
+        ("", {PETAB_EVENT: None}),
+    ],
+)
+def test_auto_event_priority(model_events, expected):
+    """Test automatic priority assignment for pre-existing events."""
+    problem = _problem_0030(model_events)
+    orig_priorities = _event_priorities(problem.model.sbml_model)
+
+    converted = ExperimentsToSbmlConverter(problem).convert()
+
+    assert _event_priorities(converted.model.sbml_model) == expected
+    # the original problem is not modified
+    assert _event_priorities(problem.model.sbml_model) == orig_priorities
+
+
+def test_auto_event_priority_empty_priority():
+    """A priority without math is treated as unset."""
+    problem = _problem_0030(E0)
+    problem.model.sbml_model.getEvent("E0").createPriority()
+
+    converted = ExperimentsToSbmlConverter(problem).convert()
+
+    assert _event_priorities(converted.model.sbml_model) == {
+        "E0": 0,
+        PETAB_EVENT: 1,
+    }
+    # the original problem is not modified
+    assert _event_priorities(problem.model.sbml_model) == {"E0": None}
+
+
+def test_explicit_default_priority():
+    """An explicit default priority takes precedence.
+
+    It only applies to events without priority."""
+    problem = _problem_0030(f"{E0}\n{E1}\n{E2_PRIO}")
+
+    converted = ExperimentsToSbmlConverter(
+        problem, default_priority=3
+    ).convert()
+
+    assert _event_priorities(converted.model.sbml_model) == {
+        "E0": 3,
+        "E1": 3,
+        "E2": 5,
+        PETAB_EVENT: 6,
+    }
+
+
+def test_infinite_priority():
+    """Test handling of a pre-existing event with priority +inf.
+
+    The condition-change events cannot get a higher priority, but a lower
+    priority can still be assigned to an unprioritized event."""
+    problem = _problem_0030(f"{E0}\n{E1_PRIO_INF}\n{E2_PRIO}")
+
+    with pytest.warns(UserWarning, match="Event `E1` has priority inf"):
+        converted = ExperimentsToSbmlConverter(problem).convert()
+
+    assert _event_priorities(converted.model.sbml_model) == {
+        "E0": 4,
+        "E1": inf,
+        "E2": 5,
+        PETAB_EVENT: inf,
+    }