ExperimentsToSbmlConverter: keep per-experiment constants in the condition table

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_<species ID>``: it holds the species' initial
+    value, the species gets the initial assignment
+    ``<species ID> = _petab_initial_value_<species ID>``, 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(
                 [

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