diff --git a/petab/v2/converters.py b/petab/v2/converters.py index 6d7d2d55..18a460d8 100644 --- a/petab/v2/converters.py +++ b/petab/v2/converters.py @@ -65,9 +65,34 @@ class ExperimentsToSbmlConverter: 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. + (if necessary), the experiment indicator parameter (if the experiment + has any changes encoded in the model), and the changes that set + *condition parameters*. + + A condition parameter is a constant SBML parameter that the model does + not assign (no rule, initial assignment or event assignment), that is + neither in the parameter table nor in the mapping table, nor an + observable or noise placeholder, and that is changed only in the first + period of an experiment, and only to numbers. + Its value is constant per experiment, like a column of a PEtab v1 + condition table, so its changes are kept as conditions of the first + period (and, for an experiment with pre-equilibration, of both kept + periods, since the value set at the start persists). + Encoding them in the model would add a branch per experiment to the + initial assignment of every such parameter. + + A species that meets the same criteria, and whose initial value in the + model is a number, is handled the same way through a new parameter + ``_petab_initial_value_``: it holds the species' initial + value, the species gets the initial assignment + `` = _petab_initial_value_``, and the changes of + the species are kept as changes of that parameter. """ + #: Prefix of the parameters that hold the initial value of a species + #: whose changes are kept in the condition table. + INITIAL_VALUE_PREFIX = "_petab_initial_value_" + #: ID of the parameter that indicates whether the model is in # the pre-equilibration phase (1) or not (0). PREEQ_INDICATOR = "_petab_preequilibration_indicator" @@ -122,6 +147,15 @@ def __init__( sciml._get_array_data_condition_ids() if sciml else set() ) + # The parameters and species whose changes stay in the condition + # table, see the class docstring. Set by `convert`. + self._condition_parameters: set[str] = set() + # species ID -> its initial value, for the species among them + self._initial_values: dict[str, float] = {} + # experiment ID -> the changes of its kept periods: the experiment + # indicator (if needed) and the changes of condition parameters + self._experiment_changes: dict[str, list[Change]] = {} + # The maximum event priority that was found in the unprocessed model. self._max_event_priority = None # The priority that will be used for the condition-change events. @@ -291,6 +325,8 @@ def convert(self) -> Problem: :return: The converted PEtab problem. """ + self._condition_parameters = self._get_condition_parameters() + self._add_initial_value_parameters() self._add_preequilibration_indicator() @@ -311,6 +347,117 @@ def convert(self) -> Problem: return self._new_problem + def _get_condition_parameters(self) -> set[str]: + """IDs of the parameters whose changes stay in the condition table. + + See the class docstring for the criteria. + """ + model = self._model + problem = self._new_problem + + assigned_in_model = { + rule.getVariable() for rule in model.getListOfRules() + } + assigned_in_model |= { + ia.getSymbol() for ia in model.getListOfInitialAssignments() + } + assigned_in_model |= { + ea.getVariable() + for event in model.getListOfEvents() + for ea in event.getListOfEventAssignments() + } + used_by_problem = {parameter.id for parameter in problem.parameters} + for mapping in problem.mappings: + used_by_problem.add(mapping.petab_id) + if mapping.model_id: + used_by_problem.add(mapping.model_id) + for observable in problem.observables: + used_by_problem |= { + str(placeholder) + for placeholder in ( + *observable.observable_placeholders, + *observable.noise_placeholders, + ) + } + + candidates = ( + ( + { + parameter.getId() + for parameter in model.getListOfParameters() + if parameter.getConstant() + } + | {species.getId() for species in model.getListOfSpecies()} + ) + - assigned_in_model + - used_by_problem + ) + + targets = set() + for experiment in problem.experiments: + for i_period, period in enumerate(experiment.sorted_periods): + if self._array_data_condition_ids and set( + period.condition_ids + ).issubset(self._array_data_condition_ids): + # set by array data, not converted + continue + for change in problem.get_changes_for_period(period): + targets.add(change.target_id) + if change.target_id in candidates and ( + i_period > 0 or not change.target_value.is_number + ): + candidates.discard(change.target_id) + + # only what conditions change + candidates &= targets + + # A species needs a numeric initial value for its parameter. The + # value of a concentration-based species given by an amount depends + # on its compartment, which must then have a fixed size. + fixed_sizes = { + sp.Symbol(compartment.getId()): sp.Float(compartment.getSize()) + for compartment in model.getListOfCompartments() + if compartment.isSetSize() + and compartment.getId() not in assigned_in_model + and compartment.getId() not in targets + } + self._initial_values = {} + for species_id in sorted(candidates): + if (species := model.getSpecies(species_id)) is None: + continue + value = self._initial_value_from_element(species).xreplace( + fixed_sizes + ) + if value.is_number: + self._initial_values[species_id] = float(value) + else: + candidates.discard(species_id) + return candidates + + def _add_initial_value_parameters(self) -> None: + """Add the parameters that hold the initial values of the species + whose changes stay in the condition table, and assign them to the + species.""" + model = self._model + for species_id, value in self._initial_values.items(): + par_id = self._condition_target(species_id) + if model.getElementBySId(par_id) is not None: + raise ValueError( + f"The model has entity with ID `{par_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_=par_id, constant=True, value=value) + set_initial_assignment(model, species_id, sp.Symbol(par_id)) + + def _condition_target(self, target_id: str) -> str: + """The ID that a kept change of `target_id` sets: the parameter + holding the initial value of a species, or the parameter itself.""" + if target_id in self._initial_values: + return f"{self.INITIAL_VALUE_PREFIX}{target_id}" + return target_id + def _convert_experiment( self, experiment: Experiment, @@ -330,7 +477,6 @@ def _convert_experiment( 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( @@ -339,7 +485,22 @@ def _convert_experiment( f"{self.__class__.__name__} and should not be used in the " "model." ) - add_sbml_parameter(model, id_=exp_ind_id, constant=False, value=0) + # the changes of the kept periods: the experiment indicator, added + # to the model only once a change is encoded there, and the changes + # of condition parameters + kept_changes: list[Change] = [] + indicator_change: list[Change] = [] + + def use_indicator() -> str: + if not indicator_change: + add_sbml_parameter( + model, id_=exp_ind_id, constant=False, value=0 + ) + indicator_change.append( + Change(target_id=exp_ind_id, target_value=1) + ) + return exp_ind_id + kept_periods: list[ExperimentPeriod] = [] for i_period, period in enumerate(experiment.sorted_periods): @@ -371,14 +532,25 @@ def _convert_experiment( # the changed values. # Additionally, tools that don't support events can still handle # single-period experiments. + # Changes of condition parameters are kept as conditions + # instead (see the class docstring). if i_period == 0: - exp_ind_id = self.get_experiment_indicator(experiment.id) - for change in self._new_problem.get_changes_for_period(period): + if change.target_id in self._condition_parameters: + kept_changes.append( + Change( + target_id=self._condition_target( + change.target_id + ), + target_value=change.target_value, + ) + ) + continue period0_assignments.setdefault( change.target_id, [] - ).append((exp_ind_id, change.target_value)) + ).append((use_indicator(), change.target_value)) else: + use_indicator() ev = self._create_period_start_event( experiment=experiment, i_period=i_period, @@ -392,10 +564,19 @@ def _convert_experiment( if len(kept_periods) > 2: raise AssertionError("Expected at most two periods to be kept.") - # add conditions that set the indicator parameters + # add conditions that set the indicator parameters and the + # condition parameters + changes = self._experiment_changes[experiment.id] = ( + indicator_change + kept_changes + ) + experiment_condition_ids = ( + [self._get_experiment_indicator_condition_id(experiment.id)] + if changes + else [] + ) for period in kept_periods: period.condition_ids = [ - self._get_experiment_indicator_condition_id(experiment.id), + *experiment_condition_ids, self.CONDITION_ID_PREEQ_ON if period.is_preequilibration else self.CONDITION_ID_PREEQ_OFF, @@ -714,24 +895,20 @@ def _add_indicators_to_conditions(self) -> None: changes=[Change(target_id=self._preeq_indicator, target_value=0)], ) - # add conditions for the experiment indicators + # add conditions for the experiment indicators and the condition + # parameters for experiment in problem.experiments: - cond_id = self._get_experiment_indicator_condition_id( - experiment.id - ) - changes = [ - Change( - target_id=self.get_experiment_indicator(experiment.id), - target_value=1, + if changes := self._experiment_changes.get(experiment.id): + problem += Condition( + id=self._get_experiment_indicator_condition_id( + experiment.id + ), + changes=changes, ) - ] - problem += 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. + # All other changes have been encoded in event assignments or + # initial assignments and can be removed. Only keep the conditions + # setting our indicators and the condition parameters. problem.condition_tables = [ ConditionTable( [ diff --git a/tests/v2/test_assignments_to_experiments.py b/tests/v2/test_assignments_to_experiments.py index 54bc667b..4b87a4c1 100644 --- a/tests/v2/test_assignments_to_experiments.py +++ b/tests/v2/test_assignments_to_experiments.py @@ -265,11 +265,19 @@ def _simulate(problem, experiment_id, timepoints): import basico converted = ExperimentsToSbmlConverter(problem).convert() - indicator = converted.model.sbml_model.getParameter( - ExperimentsToSbmlConverter.get_experiment_indicator(experiment_id) - ) - if indicator is not None: - indicator.setValue(1) + # apply what the converter leaves in the condition table: the experiment + # indicator and the condition parameters + sbml_model = converted.model.sbml_model + for experiment in converted.experiments: + if experiment.id != experiment_id: + continue + (period,) = ( + p for p in experiment.periods if not p.is_preequilibration + ) + for change in converted.get_changes_for_period(period): + sbml_model.getParameter(change.target_id).setValue( + float(change.target_value) + ) basico.load_model(converted.model.to_sbml_str()) result = basico.run_time_course( values=timepoints, r_tol=1e-10, a_tol=1e-12 diff --git a/tests/v2/test_converters.py b/tests/v2/test_converters.py index 31747237..22f840f2 100644 --- a/tests/v2/test_converters.py +++ b/tests/v2/test_converters.py @@ -101,9 +101,12 @@ def test_initial_assignments_many_experiments(): for i in range(n_experiments): problem.add_condition(f"c{i}", A=i + 2, B=i + 2) problem.add_experiment(f"e{i}", 0, f"c{i}") - # an experiment that does not change the initial values + # an experiment that does not change the initial values; it changes A + # later, so A's changes, like B's (which has an initial assignment), + # are encoded in the model rather than kept as conditions problem.add_condition("c_default", k=3) - problem.add_experiment("e_default", 0, "c_default") + problem.add_condition("c_late", A=7) + problem.add_experiment("e_default", 0, "c_default", 5, "c_late") converted = ExperimentsToSbmlConverter(problem).convert() sbml_model = converted.model.sbml_model @@ -134,6 +137,179 @@ def test_initial_assignments_many_experiments(): ) +def _kept_changes(problem: Problem) -> dict[str, list[dict[str, float]]]: + """Per experiment and kept period, the changes left in the condition + table, without the pre-equilibration indicator.""" + return { + experiment.id: [ + { + change.target_id: float(change.target_value) + for change in problem.get_changes_for_period(period) + if change.target_id + != ExperimentsToSbmlConverter.PREEQ_INDICATOR + } + for period in experiment.periods + ] + for experiment in problem.experiments + } + + +def test_condition_parameters_stay_in_conditions(): + """Numeric first-period changes of a constant parameter that nothing else + assigns stay in the condition table, instead of becoming a piecewise + initial assignment with a branch per experiment; those of a species stay + there as changes of a parameter holding its initial value. Only + experiments with changes encoded in the model get an indicator.""" + n_experiments = 300 + ant_model = """ + species A = 1 + species B = 0 + k = 3 + R1: A -> B; k * A + """ + problem = Problem() + problem.model = SbmlModel.from_antimony(ant_model) + for i in range(n_experiments): + problem.add_condition(f"c{i}", k=i + 1) + problem.add_experiment(f"e{i}", 0, f"c{i}") + problem.add_condition("c_species", k=7, A=5) + problem.add_experiment("e_species", 0, "c_species") + # a change after the start is an event, which needs the indicator + problem.add_condition("c_b", B=2) + problem.add_experiment("e_event", 0, "c0", 5, "c_b") + + converted = ExperimentsToSbmlConverter(problem).convert() + assert converted.validate().has_errors() is False + sbml_model = converted.model.sbml_model + + assert sbml_model.getInitialAssignmentBySymbol("k") is None + assert sbml_model.getParameter("k").getConstant() is True + initial_value_a = ExperimentsToSbmlConverter.INITIAL_VALUE_PREFIX + "A" + assert sbml_math_to_sympy( + sbml_model.getInitialAssignmentBySymbol("A").getMath() + ) == sp.Symbol(initial_value_a) + assert sbml_model.getParameter(initial_value_a).getValue() == 1 + assert sbml_model.getParameter(initial_value_a).getConstant() is True + assert [ + p.getId() + for p in sbml_model.getListOfParameters() + if p.getId().startswith("_petab_experiment_indicator_") + ] == ["_petab_experiment_indicator_e_event"] + assert sbml_model.getNumEvents() == 1 + + kept = _kept_changes(converted) + assert kept["e0"] == [{"k": 1.0}] + assert kept[f"e{n_experiments - 1}"] == [{"k": float(n_experiments)}] + assert kept["e_species"] == [{"k": 7.0, initial_value_a: 5.0}] + assert kept["e_event"] == [ + {"_petab_experiment_indicator_e_event": 1.0, "k": 1.0} + ] + + +@pytest.mark.parametrize( + "case", ["estimated", "later_period", "expression", "initial_assignment"] +) +def test_condition_parameters_encoded_otherwise(case): + """A parameter change that is not constant per experiment, or a parameter + that something else assigns, is still encoded in the model.""" + ant_model = """ + species A = 1 + k = 3 + p = 2 + R1: A -> ; k * A + """ + if case == "initial_assignment": + ant_model += "k = 2 * p\n" + problem = Problem() + problem.model = SbmlModel.from_antimony(ant_model) + problem.add_condition("c0", k="2 * p" if case == "expression" else 1) + if case == "later_period": + problem.add_condition("c1", k=2) + problem.add_experiment("e0", 0, "c0", 5, "c1") + else: + problem.add_experiment("e0", 0, "c0") + if case == "estimated": + problem.add_parameter("k", estimate=True, lb=0.1, ub=10) + + converted = ExperimentsToSbmlConverter(problem).convert() + sbml_model = converted.model.sbml_model + + assert sbml_model.getInitialAssignmentBySymbol("k") is not None + assert _kept_changes(converted)["e0"] == [ + {"_petab_experiment_indicator_e0": 1.0} + ] + assert sbml_model.getNumEvents() == int(case == "later_period") + # no condition sets A: it stays as it is + assert sbml_model.getInitialAssignmentBySymbol("A") is None + + +@pytest.mark.parametrize("case", ["initial_assignment", "compartment_target"]) +def test_condition_parameters_species_encoded_otherwise(case): + """A species with an initial assignment, or whose initial value depends + on a compartment that conditions resize, keeps its piecewise initial + assignment.""" + ant_model = """ + compartment C = 2 + species A in C = 1 + p = 3 + """ + if case == "initial_assignment": + ant_model += "A = 2 * p\n" + problem = Problem() + problem.model = SbmlModel.from_antimony(ant_model) + if case == "compartment_target": + # concentration-based, given as an amount: the initial + # concentration depends on the size of C + species = problem.model.sbml_model.getSpecies("A") + species.unsetInitialConcentration() + species.setInitialAmount(4) + problem.add_condition("c1", C=3) + problem.add_experiment("e1", 0, "c1") + problem.add_condition("c0", A=5) + problem.add_experiment("e0", 0, "c0") + + converted = ExperimentsToSbmlConverter(problem).convert() + sbml_model = converted.model.sbml_model + + expr = sbml_math_to_sympy( + sbml_model.getInitialAssignmentBySymbol("A").getMath() + ) + assert isinstance(expr, sp.Piecewise) + assert ( + sbml_model.getParameter( + ExperimentsToSbmlConverter.INITIAL_VALUE_PREFIX + "A" + ) + is None + ) + + +def test_condition_parameters_preequilibration(): + """A condition parameter set in the pre-equilibration period keeps its + value in the main period.""" + ant_model = """ + species A = 1 + k = 3 + R1: A -> ; k * A + """ + problem = Problem() + problem.model = SbmlModel.from_antimony(ant_model) + problem.add_condition("c_preeq", k=2) + problem.add_condition("c_main", A=5) + problem.add_experiment("e1", -inf, "c_preeq", 0, "c_main") + + converted = ExperimentsToSbmlConverter(problem).convert() + assert converted.validate().has_errors() is False + sbml_model = converted.model.sbml_model + + assert sbml_model.getInitialAssignmentBySymbol("k") is None + # the event that sets A needs e1's indicator; k is set in both periods + assert sbml_model.getNumEvents() == 1 + assert _kept_changes(converted)["e1"] == [ + {"_petab_experiment_indicator_e1": 1.0, "k": 2.0}, + {"_petab_experiment_indicator_e1": 1.0, "k": 2.0}, + ] + + def test_simulate_experiment_to_events(): """ Convert PEtab experiment to SBML events and compare BasiCO simulation