diff --git a/doc/modules.rst b/doc/modules.rst index 627ba9d8..6dacba5a 100644 --- a/doc/modules.rst +++ b/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 diff --git a/petab/v2/__init__.py b/petab/v2/__init__.py index 4f8d28ea..68069010 100644 --- a/petab/v2/__init__.py +++ b/petab/v2/__init__.py @@ -37,6 +37,7 @@ models, # noqa: F401, E402 ) from .conditions import * # noqa: F403, F401, E402 +from .core import * # noqa: F401, E402 from .experiments import ( # noqa: F401, E402 get_experiment_df, write_experiment_df, diff --git a/petab/v2/converters.py b/petab/v2/converters.py new file mode 100644 index 00000000..f0736087 --- /dev/null +++ b/petab/v2/converters.py @@ -0,0 +1,408 @@ +"""Conversion of PEtab problems.""" + +from __future__ import annotations + +import warnings +from copy import deepcopy + +import libsbml +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. + + 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`` + and one period with a finite start time. The associated changes with + these periods are only the pre-equilibration indicator + (if necessary), and the experiment indicator parameter. + """ + + #: ID of the parameter that indicates whether the model is in + # the pre-equilibration phase (1) or not (0). + PREEQ_INDICATOR = "_petab_preequilibration_indicator" + + #: The condition ID of the condition that sets the + #: pre-equilibration indicator to 1. + CONDITION_ID_PREEQ_ON = "_petab_preequilibration_on" + + #: The condition ID of the condition that sets the + #: pre-equilibration indicator to 0. + CONDITION_ID_PREEQ_OFF = "_petab_preequilibration_off" + + def __init__(self, problem: Problem, default_priority: float = None): + """Initialize the converter. + + :param problem: The PEtab problem to convert. + This will not be modified. + :param default_priority: The priority value to apply to any events that + preexist in the model and do not have a priority set. + + 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. + See SBML specs for details. + To ensure that the PEtab condition-start-events are executed before + any other events, all events should have a priority set. + """ + if not isinstance(problem.model, SbmlModel): + raise ValueError("Only SBML models are supported.") + + self._original_problem = problem + self._new_problem = deepcopy(self._original_problem) + + self._model: libsbml.Model = self._new_problem.model.sbml_model + self._preeq_indicator = self.PREEQ_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._default_priority = default_priority + self._preprocess() + + def _get_experiment_indicator_condition_id( + self, experiment_id: str + ) -> str: + """Get the condition ID for the experiment indicator parameter.""" + return f"_petab_experiment_condition_{experiment_id}" + + def _preprocess(self) -> None: + """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 of the model failed." + ) + + # 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 + ): + priority = event.createPriority() + priority.setMath( + libsbml.parseL3Formula(str(self._default_priority)) + ) + + # 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 + ) + + for event in model.getListOfEvents(): + # Check for undefined event priorities and warn + 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 " + "a PEtab condition change, otherwise the behavior is " + "undefined. To avoid this warning, see the " + "`default_priority` parameter of " + f"{self.__class__.__name__}.", + stacklevel=1, + ) + + # Check for useValuesFromTrigger time + if event.getUseValuesFromTriggerTime(): + # Non-PEtab-condition-change events must be executed *after* + # PEtab condition changes have been applied, based on the + # updated model state. This would be violated by + # useValuesFromTriggerTime=true. + warnings.warn( + f"Event `{event.getId()}` has " + "`useValuesFromTriggerTime=true'. " + "Make sure that this event cannot trigger at the time of " + "a PEtab condition change, or consider changing " + "`useValuesFromTriggerTime' to `false'. Otherwise " + "simulation results may be incorrect.", + stacklevel=1, + ) + + def convert(self) -> Problem: + """Convert the PEtab experiments to SBML events. + + :return: The converted PEtab problem. + """ + + self._add_preequilibration_indicator() + + for experiment in self._new_problem.experiment_table.experiments: + self._convert_experiment(experiment) + + self._add_indicators_to_conditions() + + validation_results = self._new_problem.validate() + validation_results.log() + + return self._new_problem + + def _convert_experiment(self, experiment: Experiment) -> None: + """Convert a single experiment to SBML events.""" + model = self._model + experiment.sort_periods() + has_preequilibration = experiment.has_preequilibration + + # add experiment indicator + exp_ind_id = self.get_experiment_indicator(experiment.id) + if model.getElementBySId(exp_ind_id) is not None: + raise ValueError( + f"The model has entity with ID `{exp_ind_id}`. " + "IDs starting with `petab_` are reserved for " + f"{self.__class__.__name__} and should not be used in the " + "model." + ) + add_sbml_parameter(model, id_=exp_ind_id, constant=False, value=0) + kept_periods = [] + for i_period, period in enumerate(experiment.periods): + if period.is_preequilibration: + # pre-equilibration cannot be represented in SBML, + # so we need to keep this period in the Problem. + kept_periods.append(period) + elif i_period == int(has_preequilibration): + # we always keep the first non-pre-equilibration period + # to set the indicator parameters + kept_periods.append(period) + elif not period.condition_ids: + # no condition, no changes, no need for an event, + # no need to keep the period unless it's the pre-equilibration + # or the only non-equilibration period (handled above) + continue + + ev = self._create_period_start_event( + experiment=experiment, + i_period=i_period, + period=period, + ) + self._create_event_assignments_for_period( + ev, + [ + self._new_problem.condition_table[condition_id] + for condition_id in period.condition_ids + ], + ) + + if len(kept_periods) > 2: + raise AssertionError("Expected at most two periods to be kept.") + + # add conditions that set the indicator parameters + for period in kept_periods: + period.condition_ids = [ + self._get_experiment_indicator_condition_id(experiment.id), + self.CONDITION_ID_PREEQ_ON + if period.is_preequilibration + else self.CONDITION_ID_PREEQ_OFF, + ] + + experiment.periods = kept_periods + + def _create_period_start_event( + self, experiment: Experiment, i_period: int, period: ExperimentPeriod + ) -> libsbml.Event: + """Create an event that triggers at the start of a period.""" + + # TODO: for now, add separate events for each experiment x period, + # this could be optimized to reuse events + + 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) + + # Create trigger expressions + # Since handling of == and !=, and distinguishing < and <= + # (and > and >=), is a bit tricky in terms of root-finding, + # we use these slightly more convoluted expressions. + # (assuming that the indicator parameters are {0, 1}) + if period.is_preequilibration: + trig_math = libsbml.parseL3Formula( + f"({exp_ind_id} > 0.5) && ({self._preeq_indicator} > 0.5)" + ) + else: + trig_math = libsbml.parseL3Formula( + f"({exp_ind_id} > 0.5) " + f"&& ({self._preeq_indicator} < 0.5) " + f"&& (time >= {period.time})" + ) + check(trigger.setMath(trig_math)) + + return ev + + def _add_preequilibration_indicator( + self, + ) -> None: + """Add an indicator parameter for the pre-equilibration to the SBML + model.""" + par_id = self._preeq_indicator + if self._model.getElementBySId(par_id) is not None: + raise ValueError( + f"Entity with ID {par_id} already exists in the SBML model." + ) + + # add the pre-steady-state indicator parameter + add_sbml_parameter(self._model, id_=par_id, value=0, constant=False) + + @staticmethod + def get_experiment_indicator(experiment_id: str) -> str: + """The ID of the experiment indicator parameter. + + The experiment indicator parameter is used to identify the + experiment in the SBML model. It is a parameter that is set + to 1 for the current experiment and 0 for all other + experiments. The parameter is used in the event trigger + to determine whether the event should be triggered. + + :param experiment_id: The ID of the experiment for which to create + the experiment indicator parameter ID. + """ + return f"_petab_experiment_indicator_{experiment_id}" + + @staticmethod + def _create_event_assignments_for_period( + event: libsbml.Event, conditions: list[Condition] + ) -> None: + """Create an event assignments for a given period.""" + for condition in conditions: + for change in condition.changes: + ExperimentsToEventsConverter._change_to_event_assignment( + change, event + ) + + @staticmethod + def _change_to_event_assignment( + change: Change, event: libsbml.Event + ) -> None: + """Convert a PEtab ``Change`` to an SBML event assignment.""" + sbml_model = event.getModel() + + ea = event.createEventAssignment() + ea.setVariable(change.target_id) + set_math(ea, change.target_value) + + # target needs const=False, and target may not exist yet + # (e.g., in case of output parameters added in the observable + # table) + target = sbml_model.getElementBySId(change.target_id) + if target is None: + add_sbml_parameter( + sbml_model, id_=change.target_id, constant=False, value=0 + ) + else: + # We can safely change the `constant` attribute of the target. + # "Constant" does not imply "boundary condition" in SBML. + target.setConstant(False) + + # the target value may depend on parameters that are only + # introduced in the PEtab parameter table - those need + # to be added to the model + for sym in change.target_value.free_symbols: + if sbml_model.getElementBySId(sym.name) is None: + add_sbml_parameter( + sbml_model, id_=sym.name, constant=True, value=0 + ) + + def _add_indicators_to_conditions(self) -> None: + """After converting the experiments to events, add the indicator + parameters for the pre-equilibration period and for the different + experiments to the remaining conditions. + Then remove all other conditions.""" + problem = self._new_problem + + # create conditions for indicator parameters + problem.condition_table.conditions.append( + Condition( + id=self.CONDITION_ID_PREEQ_ON, + changes=[ + Change(target_id=self._preeq_indicator, target_value=1) + ], + ) + ) + problem.condition_table.conditions.append( + Condition( + id=self.CONDITION_ID_PREEQ_OFF, + changes=[ + Change(target_id=self._preeq_indicator, target_value=0) + ], + ) + ) + # add conditions for the experiment indicators + for experiment in problem.experiment_table.experiments: + cond_id = self._get_experiment_indicator_condition_id( + experiment.id + ) + changes = [ + Change( + target_id=self.get_experiment_indicator(experiment.id), + target_value=1, + ) + ] + problem.condition_table.conditions.append( + Condition( + id=cond_id, + changes=changes, + ) + ) + + # All changes have been encoded in event assignments and can be + # removed. Only keep the conditions setting our indicators. + problem.condition_table.conditions = [ + condition + for condition in problem.condition_table.conditions + if condition.id.startswith("_petab") + ] diff --git a/petab/v2/core.py b/petab/v2/core.py index a847b196..193be335 100644 --- a/petab/v2/core.py +++ b/petab/v2/core.py @@ -114,6 +114,8 @@ class NoiseDistribution(str, Enum): LOG_NORMAL = C.LOG_NORMAL #: Log-Laplace distribution LOG_LAPLACE = C.LOG_LAPLACE + #: Log10-Normal + LOG10_NORMAL = C.LOG10_NORMAL class PriorDistribution(str, Enum): @@ -519,6 +521,11 @@ def _validate_ids(cls, condition_ids): raise ValueError(f"Invalid {C.CONDITION_ID}: `{condition_id}'") return condition_ids + @property + def is_preequilibration(self) -> bool: + """Check if this period is a preequilibration period.""" + return self.time == C.TIME_PREEQUILIBRATION + class Experiment(BaseModel): """An experiment or a timecourse defined by an ID and a set of different @@ -553,6 +560,15 @@ def __iadd__(self, other: ExperimentPeriod) -> Experiment: self.periods.append(other) return self + @property + def has_preequilibration(self) -> bool: + """Check if the experiment has preequilibration enabled.""" + return any(period.is_preequilibration for period in self.periods) + + def sort_periods(self) -> None: + """Sort the periods of the experiment by time.""" + self.periods.sort(key=lambda period: period.time) + class ExperimentTable(BaseModel): """PEtab experiments table.""" diff --git a/petab/v2/models/_sbml_utils.py b/petab/v2/models/_sbml_utils.py new file mode 100644 index 00000000..cbccde2b --- /dev/null +++ b/petab/v2/models/_sbml_utils.py @@ -0,0 +1,51 @@ +"""Private utility functions for SBML handling.""" + +import libsbml + +retval_to_str = { + getattr(libsbml, attr): attr + for attr in ( + "LIBSBML_DUPLICATE_OBJECT_ID", + "LIBSBML_INDEX_EXCEEDS_SIZE", + "LIBSBML_INVALID_ATTRIBUTE_VALUE", + "LIBSBML_INVALID_OBJECT", + "LIBSBML_INVALID_XML_OPERATION", + "LIBSBML_LEVEL_MISMATCH", + "LIBSBML_NAMESPACES_MISMATCH", + "LIBSBML_OPERATION_FAILED", + "LIBSBML_UNEXPECTED_ATTRIBUTE", + "LIBSBML_PKG_UNKNOWN", + "LIBSBML_PKG_VERSION_MISMATCH", + "LIBSBML_PKG_CONFLICTED_VERSION", + ) +} + + +def check(res: int): + """Check the return value of a libsbml function that returns a status code. + + :param res: The return value to check. + :raises RuntimeError: If the return value indicates an error. + """ + if res != libsbml.LIBSBML_OPERATION_SUCCESS: + raise RuntimeError(f"libsbml error: {retval_to_str.get(res, res)}") + + +def add_sbml_parameter( + model: libsbml.Model, + id_: str, + value: float = None, + constant: bool = None, +) -> libsbml.Parameter: + """Add a parameter to the SBML model.""" + param = model.createParameter() + + check(param.setId(id_)) + + if value is not None: + check(param.setValue(value)) + + if constant is not None: + check(param.setConstant(constant)) + + return param diff --git a/petab/v2/problem.py b/petab/v2/problem.py index a191942f..97684241 100644 --- a/petab/v2/problem.py +++ b/petab/v2/problem.py @@ -838,7 +838,7 @@ def validate( ) validation_results = ValidationResultList() - if self.config.extensions: + if self.config and self.config.extensions: extensions = ",".join(self.config.extensions.keys()) validation_results.append( ValidationIssue( diff --git a/tests/v2/test_converters.py b/tests/v2/test_converters.py new file mode 100644 index 00000000..76ba6a86 --- /dev/null +++ b/tests/v2/test_converters.py @@ -0,0 +1,76 @@ +from math import inf + +from petab.v2 import Change, Condition, Experiment, ExperimentPeriod, Problem +from petab.v2.converters import ExperimentsToEventsConverter +from petab.v2.models.sbml_model import SbmlModel + + +def test_experiments_to_events_converter(): + """Test the ExperimentsToEventsConverter.""" + ant_model = """ + species X = 0 + X' = 1 + """ + problem = Problem() + problem.model = SbmlModel.from_antimony(ant_model) + problem.add_condition("c1", X=1) + problem.add_condition("c2", X=2) + problem.add_experiment("e1", -inf, "c1", 10, "c2") + + converter = ExperimentsToEventsConverter(problem) + converted = converter.convert() + assert converted.validate().has_errors() is False + + assert isinstance(converted.model, SbmlModel) + sbml_model = converted.model.sbml_model + + assert sbml_model.getNumEvents() == 2 + assert converted.condition_table.conditions == [ + Condition( + id="_petab_preequilibration_on", + changes=[ + Change( + target_id="_petab_preequilibration_indicator", + target_value=1, + ) + ], + ), + Condition( + id="_petab_preequilibration_off", + changes=[ + Change( + target_id="_petab_preequilibration_indicator", + target_value=0, + ) + ], + ), + Condition( + id="_petab_experiment_condition_e1", + changes=[ + Change( + target_id="_petab_experiment_indicator_e1", target_value=1 + ) + ], + ), + ] + assert converted.experiment_table.experiments == [ + Experiment( + id="e1", + periods=[ + ExperimentPeriod( + time=-inf, + condition_ids=[ + "_petab_experiment_condition_e1", + "_petab_preequilibration_on", + ], + ), + ExperimentPeriod( + time=10.0, + condition_ids=[ + "_petab_experiment_condition_e1", + "_petab_preequilibration_off", + ], + ), + ], + ), + ]