diff --git a/CHANGELOG.md b/CHANGELOG.md index 230e675a..46ad2449 100644 --- a/CHANGELOG.md +++ b/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` diff --git a/petab/v2/converters.py b/petab/v2/converters.py index 37c5d001..5317d270 100644 --- a/petab/v2/converters.py +++ b/petab/v2/converters.py @@ -2,6 +2,8 @@ from __future__ import annotations +import logging +import math import warnings from copy import deepcopy @@ -22,6 +24,8 @@ __all__ = ["ExperimentsToSbmlConverter", "assignments_to_experiments"] +logger = logging.getLogger(__name__) + class ExperimentsToSbmlConverter: """Convert PEtab experiments to SBML. @@ -29,11 +33,17 @@ class ExperimentsToSbmlConverter: 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. diff --git a/tests/v2/test_converters.py b/tests/v2/test_converters.py index c38238f0..31747237 100644 --- a/tests/v2/test_converters.py +++ b/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, + }