From 54ba6efedf22b40a00a325b12b0967d66ded6261 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Sun, 25 Oct 2026 12:26:49 +0000 Subject: [PATCH 1/8] assignments_to_experiments: fix lost switches, dropped changes and scope gaps - A switch exactly at the start of a later period was lost when that period did not change any symbol of the assignment, e.g. `piecewise(0, time < 10, 1)` with a period at t = 10, or `time <= 0` after preequilibration. A later period now also sets the value when it differs from the value at the end of the previous period. - The condition changes of every symbol of a converted assignment were dropped, even when the model still used the symbol elsewhere, e.g. a dose that also enters a reaction. They are now only dropped if nothing else (the model, an observable, a condition change) refers to the symbol. - Time-dependent initial assignments were converted into switches, but SBML evaluates them once, at the start of the simulation. They now stay in the model. - At preequilibration, which starts the experiment, a value that referred to model parameters outside the parameter table failed the lint. Those parameters now take their model values, as at the first period. - An assignment to a species or compartment crashed with an AttributeError; it now raises NotImplementedError. - A problem without time-dependent assignments was still changed: an experiment that only preequilibrates got an empty period at t = 0. `test_assignments_to_experiments` and `_unused_condition_id` used an initial assignment to mean a switching rule and now use an assignment rule. `_validate` relied on an estimated parameter with a time-dependent initial assignment; it now uses a ramp whose generated slope parameter collides with an observable parameter. The converted benchmark problems (Brannmark_JBC2010, Fujita_SciSignal2010, Isensee_JCB2018, Raimundez_PCB2020) are unchanged. Co-Authored-By: Claude Opus 5.5 --- petab/v2/converters.py | 115 ++++++++++++++++----- tests/v2/test_conversion.py | 195 ++++++++++++++++++++++++++++++++++-- 2 files changed, 278 insertions(+), 32 deletions(-) diff --git a/petab/v2/converters.py b/petab/v2/converters.py index 9eef1061..b0f5b0b5 100644 --- a/petab/v2/converters.py +++ b/petab/v2/converters.py @@ -757,9 +757,9 @@ def assignments_to_experiments(problem: Problem) -> Problem: if all of them were replaced, since a period's time overrides the model's initial time. - **Scope:** conversion is supported for SBML assignment rules that depend - on time through the conditions of piecewise functions, and at most - linearly through their values, where + **Scope:** conversion is supported for SBML assignment rules to + parameters that depend on time through the conditions of piecewise + functions, and at most linearly through their values, where * every such condition is an inequality that can be re-arranged explicitly into `t />= ...`, @@ -781,7 +781,11 @@ def assignments_to_experiments(problem: Problem) -> Problem: entity instead, with the rate rule `d target/dt = _slope`, and each period sets both the target and its slope; an assignment that only switches between constant values stays a constant parameter. Anything else - raises a :class:`NotImplementedError`. + raises a :class:`NotImplementedError`. Initial assignments stay in the + model, since they are evaluated once, at the start of the simulation. + + Condition changes of the replaced assignments' symbols are dropped, unless + the model, an observable or another change still refers to them. .. warning:: @@ -807,25 +811,29 @@ def assignments_to_experiments(problem: Problem) -> Problem: problem = deepcopy(problem) sbml_model = problem.model.sbml_model + # Initial assignments are evaluated once, at the start of the + # simulation, so they never switch and stay in the model. + time_dependent = [] + for rule in sbml_model.getListOfRules(): + if not rule.isAssignment(): + continue + expr = sbml_math_to_sympy(rule) + if any(pw.has(TimeSymbol) for pw in expr.atoms(sp.Piecewise)): + time_dependent.append((rule, expr)) + experiments = list(problem.experiments) if any(m.experiment_id is None for m in problem.measurements): default_experiment = Experiment(id=DEFAULT_EXPERIMENT_ID) experiments.append(default_experiment) else: default_experiment = None - droppable = _add_start_periods(experiments) + droppable = _add_start_periods(experiments) if time_dependent else set() # (experiment, period start time, target ID, target value) additions = [] replaced_symbols = set() converted_ids = [] - for rule in [ - *(r for r in sbml_model.getListOfRules() if r.isAssignment()), - *sbml_model.getListOfInitialAssignments(), - ]: - expr = sbml_math_to_sympy(rule) - if not any(pw.has(TimeSymbol) for pw in expr.atoms(sp.Piecewise)): - continue + for rule, expr in time_dependent: assignment = _PiecewiseAssignment(rule, expr, problem) assignment.replace_in_model() for experiment in experiments: @@ -846,7 +854,9 @@ def assignments_to_experiments(problem: Problem) -> Problem: stacklevel=2, ) - # replace the changes that are covered by the new conditions + # replace the changes that are covered by the new conditions, unless + # something else still refers to their targets + replaced_symbols -= _referenced_ids(problem, additions) _drop_changes(problem, experiments, replaced_symbols) _add_conditions(problem, additions) @@ -893,6 +903,11 @@ def __init__( self.problem = problem self.sbml_model = problem.model.sbml_model self.target_id = rule.getId() + if self.sbml_model.getParameter(self.target_id) is None: + raise NotImplementedError( + "Only assignments to parameters are supported, but " + f"`{self.target_id}` is not a parameter." + ) self.expr = expr (self.time,) = expr.atoms(TimeSymbol) self.symbol_ids = {str(s) for s in expr.free_symbols} @@ -977,6 +992,8 @@ def changes( changes = [] periods = experiment.sorted_periods overrides = {} + # the value at the end of the previous period + held = None for i_period, period in enumerate(periods): period_changes = self.problem.get_changes_for_period(period) overrides |= {c.target_id: c.target_value for c in period_changes} @@ -991,8 +1008,14 @@ def changes( # the value at time 0, which respects whether a condition # holds for a switch at time 0, held there: a ramp does not # run during preequilibration - value = self.expr.subs(subs).subs(self.time, 0) - changes.append((period.time, self.target_id, value)) + held = self.expr.subs(subs).subs(self.time, 0) + changes.append( + ( + period.time, + self.target_id, + self._with_model_values(held), + ) + ) if self.ramps: changes.append((period.time, self.slope_id, sp.S.Zero)) continue @@ -1024,31 +1047,38 @@ def changes( or {str(s) for s in self.default.free_symbols} & self.parameter_ids ): - # at the start of an experiment, a change may only refer - # to parameter table symbols; other symbols are - # constant model parameters, at their values - start_value = start_value.subs( - { - s: self._model_value(s) - for s in start_value.free_symbols - {self.time} - if str(s) not in self.parameter_ids - } + changes += self._start( + period.time, self._with_model_values(start_value) ) - changes += self._start(period.time, start_value) elif ( {c.target_id for c in period_changes} & self.symbol_ids # preequilibration held the slope at zero or (self.ramps and periods[i_period - 1].is_preequilibration) + # a switch at the start of the period + or not (start_value - held).is_zero ): changes += self._start(period.time, start_value) for t, value in switch_values: changes += self._start(t, value) + held = previous return changes def _model_value(self, symbol: sp.Symbol) -> float: """The value of the model parameter ``symbol``.""" return self.sbml_model.getParameter(str(symbol)).getValue() + def _with_model_values(self, value: sp.Expr) -> sp.Expr: + """``value`` for a change at the start of an experiment, which may + only refer to parameter table symbols: the other symbols are constant + model parameters, at their values.""" + return value.subs( + { + s: self._model_value(s) + for s in value.free_symbols - {self.time} + if str(s) not in self.parameter_ids + } + ) + def _switching_times(self, subs: dict, experiment_id: str) -> list[float]: """The switching times, for the symbol values ``subs``.""" times = [] @@ -1109,6 +1139,41 @@ def _add_start_periods(experiments: list[Experiment]) -> set[int]: return droppable +def _referenced_ids( + problem: Problem, + additions: list[tuple[Experiment, float, str, sp.Expr]], +) -> set[str]: + """The IDs that the math of the model, the observables, or the condition + changes, including those in ``additions``, refer to.""" + ids = set() + for element in problem.model.sbml_model.getListOfAllElements(): + if not (hasattr(element, "isSetMath") and element.isSetMath()): + continue + nodes = [element.getMath()] + while nodes: + node = nodes.pop() + if node.getType() == libsbml.AST_NAME: + ids.add(node.getName()) + nodes += [node.getChild(i) for i in range(node.getNumChildren())] + expressions = [ + *(o.formula for o in problem.observables), + *(o.noise_formula for o in problem.observables), + *(value for *_, value in additions), + ] + ids |= { + str(s) + for expression in expressions + if expression is not None + for s in expression.free_symbols + } + ids |= { + str(s) + for table in problem.condition_tables + for s in table.free_symbols + } + return ids + + def _drop_changes( problem: Problem, experiments: list[Experiment], target_ids: set[str] ) -> None: diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 201a6740..c558d9d7 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -670,7 +670,7 @@ def test_assignments_to_experiments(condition, before, after): species A, B; A = 10; B = 0; - k1 = piecewise(1, {condition}, 2); + k1 := piecewise(1, {condition}, 2); k2 = 0.5; R1: A -> B; k1 * A; R2: B -> A; k2 * B; @@ -719,7 +719,7 @@ def test_assignments_to_experiments_unused_condition_id(): problem = Problem() problem.model = SbmlModel.from_antimony(""" species A = 10; - k1 = piecewise(1, time < 10, 2); + k1 := piecewise(1, time < 10, 2); R1: A -> ; k1 * A; """) problem.add_observable("obs_a", formula="A", noise_formula="cond_0") @@ -1477,8 +1477,9 @@ def test_assignments_to_experiments_preequilibration_boundary( def test_assignments_to_experiments_validate(): """Test that the converted problem is linted. - Here, the input problem is valid, but the conversion makes the estimated - parameter `k1` a condition table target. + Here, the input problem is valid, but the conversion adds the model + parameter `k1_slope` for the ramp, which is also an estimated observable + parameter in the parameter table. """ problem = Problem() ant_model = """ @@ -1486,18 +1487,198 @@ def test_assignments_to_experiments_validate(): species A, B; A = 10; B = 0; - k1 = piecewise(0, time < 10, 1); + k1 := piecewise(time, time < 10, 1); k2 = 0.5; R1: A -> B; k1 * A; R2: B -> A; k2 * B; end """ problem.model = SbmlModel.from_antimony(ant_model) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_observable("obs_a", formula="k1_slope * A", noise_formula=1.0) problem.add_measurement("obs_a", time=20.0, measurement=5.0) - problem.add_parameter("k1", estimate=True, lb=1e-3, ub=1e3) + problem.add_parameter("k1_slope", estimate=True, lb=1e-3, ub=1e3) problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) assert not problem.validate().has_errors() with pytest.raises(ValueError, match="did not pass linting"): assignments_to_experiments(problem) + + +def _changes_by_time( + problem: Problem, experiment_id: str +) -> dict[float, dict[str, float]]: + """The changes of each period of an experiment, by the period's time.""" + return { + period.time: { + c.target_id: float(c.target_value) + for c in problem.get_changes_for_period(period) + } + for period in problem[experiment_id].sorted_periods + } + + +def _switch_model(assignment: str, extra: str = "") -> SbmlModel: + """A model with the assignment rule `k1 := `.""" + return SbmlModel.from_antimony(f""" + model conversion + species A; + A = 10; + other = 1; + {extra} + k1 := {assignment}; + R1: A -> ; k1 * A; + end + """) + + +def test_assignments_to_experiments_switch_at_period_start(): + """Test that a switch at the start of a later period is kept, although + that period does not change any symbol of the assignment.""" + problem = Problem() + problem.model = _switch_model("piecewise(0, time < 10, 1)") + problem.add_condition("c0", other=1.0) + problem.add_condition("c1", other=2.0) + problem.add_experiment("e1", 0, "c0", 10, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + assert _changes_by_time(converted, "e1") == { + 0.0: {"other": 1.0}, + 10.0: {"other": 2.0, "k1": 1.0}, + } + + +def test_assignments_to_experiments_switch_at_preequilibration_end(): + """Test that a switch at time 0 is kept after preequilibration, where the + condition still holds at time 0.""" + problem = Problem() + problem.model = _switch_model("piecewise(0, time <= 0, 1)") + problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, [], 0, []) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + with pytest.warns(UserWarning, match="during preequilibration"): + converted = assignments_to_experiments(problem) + + assert _changes_by_time(converted, "e1") == { + v2.C.TIME_PREEQUILIBRATION: {"k1": 0.0}, + 0.0: {"k1": 1.0}, + } + + +def test_assignments_to_experiments_shared_symbol(): + """Test that the changes of a symbol of a converted assignment stay if + the model still refers to that symbol elsewhere.""" + problem = Problem() + problem.model = _switch_model( + "dose * piecewise(0, time < 10, 1)", + extra="dose = 0; R2: -> A; dose;", + ) + problem.add_condition("c1", dose=2.0) + problem.add_experiment("e1", 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + assert _changes_by_time(converted, "e1") == { + 0.0: {"dose": 2.0}, + 10.0: {"k1": 2.0}, + } + + +def test_assignments_to_experiments_initial_assignment(): + """Test that a time-dependent initial assignment stays in the model, since + it is evaluated once, at the start of the simulation.""" + problem = Problem() + problem.model = SbmlModel.from_antimony(""" + model conversion + species A; + A = 10; + k1 = piecewise(0, time < 10, 1); + R1: A -> ; k1 * A; + end + """) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=1.0) + + converted = assignments_to_experiments(problem) + + assert converted.model.sbml_model.getInitialAssignment("k1") is not None + assert not converted.experiments + assert not converted.conditions + + +def test_assignments_to_experiments_preequilibration_model_value(): + """Test that a value at preequilibration, which starts the experiment, + only refers to parameter table symbols: other model parameters take their + values.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(level, time < 10, 0)", extra="level = 3;" + ) + problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, [], 0, []) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + with pytest.warns(UserWarning, match="during preequilibration"): + converted = assignments_to_experiments(problem) + + assert _changes_by_time(converted, "e1") == { + v2.C.TIME_PREEQUILIBRATION: {"k1": 3.0}, + 0.0: {}, + 10.0: {"k1": 0.0}, + } + + +def test_assignments_to_experiments_species_target(): + """Test that an assignment to a species is rejected.""" + problem = Problem() + problem.model = SbmlModel.from_antimony(""" + model conversion + species A, S; + A = 10; + S := piecewise(0, time < 10, 1); + R1: A -> ; S * A; + end + """) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=1.0) + + with pytest.raises(NotImplementedError, match="not a parameter"): + assignments_to_experiments(problem) + + +def test_assignments_to_experiments_nothing_to_convert(): + """Test that a problem without time-dependent assignments is unchanged, + including an experiment that only preequilibrates.""" + problem = Problem() + problem.model = SbmlModel.from_antimony(""" + model conversion + species A; + A = 10; + k2 = 0.5; + R1: A -> ; k2 * A; + end + """) + problem.add_condition("c1", k2=1.0) + problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + assert converted.experiment_df.equals(problem.experiment_df) + assert converted.condition_df.equals(problem.condition_df) From 14103608cdaf505a2784cc46468bdb7e6de5c12f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Sun, 25 Oct 2026 16:10:25 +0000 Subject: [PATCH 2/8] assignments_to_experiments: count conversion factors as uses, leave unconverted problems alone, pick an unused slope ID - A parameter used as a model or species conversion factor counted as unused, so a condition change of it was dropped when the parameter also appeared in a converted assignment. Conversion factors now count as references. - A problem without time-dependent assignments was still changed: its experiments' periods were re-sorted. It is now returned as it was, after the lint. - A ramp's slope parameter was only named so as not to clash with SBML IDs, and clashed with, e.g., an observable parameter `_slope`. It is now named with `Problem.get_unused_id`. `test_assignments_to_experiments_validate` relied on that slope clash to fail the lint; it now uses an existing experiment `default_experiment`, which the measurements without an experiment ID are assigned to as well. Co-Authored-By: Claude Opus 5.5 --- petab/v2/converters.py | 33 +++++++++---- tests/v2/test_conversion.py | 96 ++++++++++++++++++++++++++++--------- 2 files changed, 98 insertions(+), 31 deletions(-) diff --git a/petab/v2/converters.py b/petab/v2/converters.py index b0f5b0b5..dd14c366 100644 --- a/petab/v2/converters.py +++ b/petab/v2/converters.py @@ -820,6 +820,9 @@ def assignments_to_experiments(problem: Problem) -> Problem: expr = sbml_math_to_sympy(rule) if any(pw.has(TimeSymbol) for pw in expr.atoms(sp.Piecewise)): time_dependent.append((rule, expr)) + if not time_dependent: + _check_lint(problem) + return problem experiments = list(problem.experiments) if any(m.experiment_id is None for m in problem.measurements): @@ -827,7 +830,7 @@ def assignments_to_experiments(problem: Problem) -> Problem: experiments.append(default_experiment) else: default_experiment = None - droppable = _add_start_periods(experiments) if time_dependent else set() + droppable = _add_start_periods(experiments) # (experiment, period start time, target ID, target value) additions = [] @@ -876,6 +879,13 @@ def assignments_to_experiments(problem: Problem) -> Problem: if measurement.experiment_id is None: measurement.experiment_id = DEFAULT_EXPERIMENT_ID + _check_lint(problem) + return problem + + +def _check_lint(problem: Problem) -> None: + """Raise a :class:`ValueError` if the converted ``problem`` does not + pass linting.""" validation_issues = problem.validate() if validation_issues.has_errors(): messages = "\n".join(f"{issue.message}" for issue in validation_issues) @@ -883,8 +893,6 @@ def assignments_to_experiments(problem: Problem) -> Problem: f"The converted PEtab problem did not pass linting:\n{messages}" ) - return problem - class _PiecewiseAssignment: """A time-dependent piecewise assignment to a model parameter, and the @@ -964,9 +972,9 @@ def replace_in_model(self) -> None: target's, is set by every period that starts one. """ if self.ramps: - self.slope_id = f"{self.target_id}_slope" - while self.sbml_model.getElementBySId(self.slope_id) is not None: - self.slope_id += "_" + self.slope_id = self.problem.get_unused_id( + f"{self.target_id}_slope" + ) add_sbml_parameter( self.sbml_model, self.slope_id, value=0.0, constant=False ) @@ -1143,10 +1151,11 @@ def _referenced_ids( problem: Problem, additions: list[tuple[Experiment, float, str, sp.Expr]], ) -> set[str]: - """The IDs that the math of the model, the observables, or the condition - changes, including those in ``additions``, refer to.""" + """The IDs that the model, the observables, or the condition changes, + including those in ``additions``, refer to.""" + sbml_model = problem.model.sbml_model ids = set() - for element in problem.model.sbml_model.getListOfAllElements(): + for element in sbml_model.getListOfAllElements(): if not (hasattr(element, "isSetMath") and element.isSetMath()): continue nodes = [element.getMath()] @@ -1155,6 +1164,12 @@ def _referenced_ids( if node.getType() == libsbml.AST_NAME: ids.add(node.getName()) nodes += [node.getChild(i) for i in range(node.getNumChildren())] + # conversion factors refer to parameters by attribute, not in math + ids |= { + element.getConversionFactor() + for element in [sbml_model, *sbml_model.getListOfSpecies()] + if element.isSetConversionFactor() + } expressions = [ *(o.formula for o in problem.observables), *(o.noise_formula for o in problem.observables), diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index c558d9d7..beacc493 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -1477,27 +1477,22 @@ def test_assignments_to_experiments_preequilibration_boundary( def test_assignments_to_experiments_validate(): """Test that the converted problem is linted. - Here, the input problem is valid, but the conversion adds the model - parameter `k1_slope` for the ramp, which is also an estimated observable - parameter in the parameter table. + Here, the input problem is valid, but the measurement without an + experiment ID is assigned to a new experiment `default_experiment`, an ID + that an existing experiment already has. """ problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - k1 := piecewise(time, time < 10, 1); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - problem.add_observable("obs_a", formula="k1_slope * A", noise_formula=1.0) - problem.add_measurement("obs_a", time=20.0, measurement=5.0) - problem.add_parameter("k1_slope", estimate=True, lb=1e-3, ub=1e3) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + problem.model = _switch_model("piecewise(0, time < 10, 1)") + problem.add_condition("c1", other=2.0) + problem.add_experiment(DEFAULT_EXPERIMENT_ID, 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", + experiment_id=DEFAULT_EXPERIMENT_ID, + time=20.0, + measurement=5.0, + ) + problem.add_measurement("obs_a", time=30.0, measurement=5.0) assert not problem.validate().has_errors() with pytest.raises(ValueError, match="did not pass linting"): @@ -1661,7 +1656,8 @@ def test_assignments_to_experiments_species_target(): def test_assignments_to_experiments_nothing_to_convert(): """Test that a problem without time-dependent assignments is unchanged, - including an experiment that only preequilibrates.""" + including an experiment that only preequilibrates and periods that are + not in time order.""" problem = Problem() problem.model = SbmlModel.from_antimony(""" model conversion @@ -1672,7 +1668,43 @@ def test_assignments_to_experiments_nothing_to_convert(): end """) problem.add_condition("c1", k2=1.0) + problem.add_condition("c2", k2=2.0) problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, "c1") + # periods out of time order + problem.add_experiment("e2", 10, "c2", 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + for experiment_id in ("e1", "e2"): + problem.add_measurement( + "obs_a", experiment_id=experiment_id, time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + def periods(problem): + return [ + [(p.time, p.condition_ids) for p in e.periods] + for e in problem.experiments + ] + + assert periods(converted) == periods(problem) + assert converted.condition_df.equals(problem.condition_df) + + +@pytest.mark.parametrize("on", ["model", "species"]) +def test_assignments_to_experiments_conversion_factor(on): + """Test that the changes of a symbol of a converted assignment stay if + the model refers to that symbol as a conversion factor.""" + problem = Problem() + problem.model = _switch_model( + "cf * piecewise(0, time < 10, 1)", extra="cf = 1;" + ) + sbml_model = problem.model.sbml_model + element = sbml_model if on == "model" else sbml_model.getSpecies("A") + assert ( + element.setConversionFactor("cf") == libsbml.LIBSBML_OPERATION_SUCCESS + ) + problem.add_condition("c1", cf=2.0) + problem.add_experiment("e1", 0, "c1") problem.add_observable("obs_a", formula="A", noise_formula=1.0) problem.add_measurement( "obs_a", experiment_id="e1", time=20.0, measurement=1.0 @@ -1680,5 +1712,25 @@ def test_assignments_to_experiments_nothing_to_convert(): converted = assignments_to_experiments(problem) - assert converted.experiment_df.equals(problem.experiment_df) - assert converted.condition_df.equals(problem.condition_df) + assert _changes_by_time(converted, "e1") == { + 0.0: {"cf": 2.0}, + 10.0: {"k1": 2.0}, + } + + +def test_assignments_to_experiments_slope_id_taken(): + """Test that a ramp's slope parameter gets an ID that the problem does not + use yet, here because an observable parameter is called `k1_slope`.""" + problem = Problem() + problem.model = _switch_model("piecewise(time, time < 10, 1)") + problem.add_observable("obs_a", formula="k1_slope * A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=5.0) + problem.add_parameter("k1_slope", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + assert _rate_rule(converted, "k1") == "k1_slope_1" + assert _changes_by_time(converted, DEFAULT_EXPERIMENT_ID) == { + 0.0: {"k1": 0.0, "k1_slope_1": 1.0}, + 10.0: {"k1": 1.0, "k1_slope_1": 0.0}, + } From d3a085be28ef8b738caec305484839d820987f30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Sun, 25 Oct 2026 16:21:52 +0000 Subject: [PATCH 3/8] assignments_to_experiments: evaluate a condition change once, at the start of its period The converter kept a condition change's value as an expression and substituted it wherever a switch used it, so the value was re-evaluated at the switch with the values at that time. A condition change is evaluated once, at the start of its period, and all changes of a period simultaneously, with the values from before the period. E.g., with `d1 = b1 + 0.5` at t = 2 and `b1 = 2` at t = 3, a switch to `d1` at t = 5 used 2.5 instead of 1.7; likewise with `b1 = 2` in the same condition as `d1`. A change's value is now evaluated at the start of its period: symbols set by an earlier condition take that value, other constant model parameters their model values, and parameter table symbols stay. A change of a symbol of the assignment whose value still refers to anything else, e.g. a species, is rejected, since its value is not known there. Found by differential testing against simulations of the original and the converted problem. Co-Authored-By: Claude Opus 5.5 --- petab/v2/converters.py | 45 ++++++++++++++++++++++++++++- tests/v2/test_conversion.py | 57 +++++++++++++++++++++++++++++++++++++ 2 files changed, 101 insertions(+), 1 deletion(-) diff --git a/petab/v2/converters.py b/petab/v2/converters.py index dd14c366..b5666f04 100644 --- a/petab/v2/converters.py +++ b/petab/v2/converters.py @@ -1004,12 +1004,32 @@ def changes( held = None for i_period, period in enumerate(periods): period_changes = self.problem.get_changes_for_period(period) - overrides |= {c.target_id: c.target_value for c in period_changes} + # a condition's changes are evaluated once, at the start of their + # period, all with the values from before the period + overrides |= { + c.target_id: self._at_period_start(c.target_value, overrides) + for c in period_changes + } subs = self.fixed_values | { s: overrides[str(s)] for s in self.expr.free_symbols if str(s) in overrides } + if unknown := sorted( + { + str(x) + for s in self.expr.free_symbols + if str(s) in overrides + for x in overrides[str(s)].free_symbols + } + - self.parameter_ids + ): + raise NotImplementedError( + "Only condition changes whose values are known at the " + "start of their period are supported, but in experiment " + f"`{experiment.id}`, the changes of the symbols of the " + f"assignment to `{self.target_id}` refer to {unknown}." + ) times = self._switching_times(subs, experiment.id) if period.is_preequilibration: @@ -1075,6 +1095,29 @@ def _model_value(self, symbol: sp.Symbol) -> float: """The value of the model parameter ``symbol``.""" return self.sbml_model.getParameter(str(symbol)).getValue() + def _at_period_start(self, value: sp.Expr, overrides: dict) -> sp.Expr: + """The ``value`` of a change, evaluated at the start of its period. + + Symbols that a condition set before take that value (``overrides``), + and the other constant model parameters outside the parameter table + their model values. Parameter table symbols stay, since no condition + sets them; any other symbol stays as well, since its value is not + known here. + """ + if value is None: + return value + known = {} + for s in value.free_symbols: + if str(s) in overrides: + known[s] = overrides[str(s)] + elif ( + str(s) not in self.parameter_ids + and (p := self.sbml_model.getParameter(str(s))) is not None + and p.getConstant() + ): + known[s] = sp.Float(p.getValue()) + return value.xreplace(known) + def _with_model_values(self, value: sp.Expr) -> sp.Expr: """``value`` for a change at the start of an experiment, which may only refer to parameter table symbols: the other symbols are constant diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index beacc493..a3994cd4 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -1734,3 +1734,60 @@ def test_assignments_to_experiments_slope_id_taken(): 0.0: {"k1": 0.0, "k1_slope_1": 1.0}, 10.0: {"k1": 1.0, "k1_slope_1": 0.0}, } + + +@pytest.mark.parametrize( + "conditions", + [ + # the changes of a condition are simultaneous: `d1` takes the old `b1` + [(2.0, {"b1": 2.0, "d1": "b1 + 0.5"})], + # a change is evaluated once: `d1` keeps the `b1` of its period + [(2.0, {"d1": "b1 + 0.5"}), (3.0, {"b1": 2.0})], + ], +) +def test_assignments_to_experiments_change_evaluated_at_period_start( + conditions, +): + """Test that a condition change's value is evaluated at the start of its + period, with the values from before the period, also when a later switch + uses it.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(b0, time < 5, d1)", extra="b0 = 0.3; b1 = 1.2; d1 = 1;" + ) + args = [0, []] + for i, (t, changes) in enumerate(conditions): + problem.add_condition(f"c{i}", **changes) + args += [t, f"c{i}"] + problem.add_experiment("e1", *args) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + (switch,) = (p for p in converted["e1"].periods if p.time == 5.0) + (change,) = converted.get_changes_for_period(switch) + assert change.target_id == "k1" + # 1.2 + 0.5, not 2.0 + 0.5 + assert float(change.target_value) == pytest.approx(1.7) + + +def test_assignments_to_experiments_change_refers_to_species(): + """Test that a change of a symbol of the assignment whose value refers to + a species, which is not known at the start of the period, is + rejected.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(0, time < 5, d1)", extra="d1 = 1;" + ) + problem.add_condition("c1", d1="A") + problem.add_experiment("e1", 0, [], 2, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + with pytest.raises(NotImplementedError, match=r"refer to \['A'\]"): + assignments_to_experiments(problem) From 0bd974e7d6fc376ed63fbcb12a3c83f8e4512c8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Sun, 25 Oct 2026 16:49:35 +0000 Subject: [PATCH 4/8] Add a CHANGELOG entry for the assignments_to_experiments fixes Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0e818060..7a799186 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -59,6 +59,14 @@ 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) +* Several fixes to `assignments_to_experiments`: keep a switch at the start + of a later period, keep the condition changes of symbols that the model + still uses elsewhere, evaluate a condition change once at the start of its + period, leave time-dependent initial assignments in the model, reject + assignments to species or compartments, leave problems without + time-dependent assignments unchanged, and give a ramp's slope parameter an + unused ID + (by @FFroehlich in https://github.com/PEtab-dev/libpetab-python/pull/574) **Features** * `petab1to2` and `petab_files_1to2` can skip validation via `validate=False` From 02f33ee1036e7d2a4567a2062297f21f40d50b32 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Sun, 25 Oct 2026 17:31:23 +0000 Subject: [PATCH 5/8] assignments_to_experiments: use PEtab symbols in new conditions, and fixed parameter table values in condition changes - The values of the new condition changes used the symbols parsed from SBML, without the assumptions of the symbols PEtab creates (`petab_symbol`, real), so they compared unequal to the same expressions parsed by PEtab, e.g. after reading the condition table back. They now use PEtab's symbols. - A condition change evaluated at the start of its period kept parameter table symbols, also those of fixed parameters. A switching time set by a condition to, e.g., `p + 1` with `p` fixed by the parameter table was therefore rejected, although `t1 := p + 1` in the assignment itself is converted with `p`'s nominal value. Fixed parameters now take their nominal values there as well; estimated ones stay, and switching times that depend on them are still rejected. Both found by differential testing with a reference simulator. Co-Authored-By: Claude Opus 5.5 --- petab/v2/converters.py | 25 +++++++++++----- tests/v2/test_conversion.py | 57 +++++++++++++++++++++++++++++++++++++ 2 files changed, 75 insertions(+), 7 deletions(-) diff --git a/petab/v2/converters.py b/petab/v2/converters.py index b5666f04..2734151a 100644 --- a/petab/v2/converters.py +++ b/petab/v2/converters.py @@ -19,6 +19,7 @@ ExperimentPeriod, Problem, ) +from .math import petab_symbol from .models._sbml_utils import ( add_sbml_parameter, check, @@ -947,11 +948,11 @@ def __init__( ) # unless a condition sets them, the switching times' symbols keep # their value from the parameter table or else from the model - nominal_values = { + self.nominal_values = { p.id: p.nominal_value for p in problem.parameters if not p.estimate } self.fixed_values = { - s: nominal_values.get(str(s), self._model_value(s)) + s: self.nominal_values.get(str(s), self._model_value(s)) for t_switch, _ in self.switches.values() for s in t_switch.free_symbols } @@ -1099,10 +1100,10 @@ def _at_period_start(self, value: sp.Expr, overrides: dict) -> sp.Expr: """The ``value`` of a change, evaluated at the start of its period. Symbols that a condition set before take that value (``overrides``), - and the other constant model parameters outside the parameter table - their model values. Parameter table symbols stay, since no condition - sets them; any other symbol stays as well, since its value is not - known here. + parameters fixed by the parameter table their nominal values, and the + other constant model parameters outside the parameter table their + model values. Estimated parameters stay, since no condition sets them; + any other symbol stays as well, since its value is not known here. """ if value is None: return value @@ -1110,6 +1111,8 @@ def _at_period_start(self, value: sp.Expr, overrides: dict) -> sp.Expr: for s in value.free_symbols: if str(s) in overrides: known[s] = overrides[str(s)] + elif str(s) in self.nominal_values: + known[s] = sp.Float(self.nominal_values[str(s)]) elif ( str(s) not in self.parameter_ids and (p := self.sbml_model.getParameter(str(s))) is not None @@ -1260,10 +1263,18 @@ def _add_conditions( add a condition with that change to the period of ``experiment`` that starts at ``time``, which is added if needed. - Equal changes share a condition. + Equal changes share a condition. Their values use the symbols that PEtab + uses for the condition table, rather than those parsed from SBML. """ new_conditions = {} for experiment, t, target_id, value in additions: + value = value.xreplace( + { + s: petab_symbol(s.name) + for s in value.free_symbols + if not isinstance(s, TimeSymbol) + } + ) if (cond_id := new_conditions.get((target_id, value))) is None: cond_id = problem.get_unused_id(f"cond_{len(new_conditions)}") new_conditions[target_id, value] = cond_id diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index a3994cd4..72deade4 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -17,6 +17,7 @@ ExperimentsToSbmlConverter, assignments_to_experiments, ) +from petab.v2.math import sympify_petab from petab.v2.petab1to2 import petab1to2, v1v2_observable_df @@ -1791,3 +1792,59 @@ def test_assignments_to_experiments_change_refers_to_species(): with pytest.raises(NotImplementedError, match=r"refer to \['A'\]"): assignments_to_experiments(problem) + + +def test_assignments_to_experiments_petab_symbols(): + """Test that the values of new condition changes use PEtab's symbols, so + that they equal the same expressions parsed by PEtab.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(0, time < 5, d1)", extra="d1 = 1;" + ) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=1.0) + problem.add_parameter( + "d1", estimate=True, nominal_value=2.0, lb=0.1, ub=10 + ) + + converted = assignments_to_experiments(problem) + + assert [ + change.target_value + for condition in converted.conditions + for change in condition.changes + ] == [sympify_petab("d1")] + + +@pytest.mark.parametrize("estimate", [False, True]) +def test_assignments_to_experiments_switching_time_from_fixed_parameter( + estimate, +): + """Test that a condition may set a switching time to an expression in a + parameter that the parameter table fixes, which takes its nominal value, + but not in an estimated one.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(0, time < t1, 1)", extra="t1 = 0; p = 1;" + ) + problem.add_condition("c1", t1="p + 1") + problem.add_experiment("e1", 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + if estimate: + problem.add_parameter( + "p", estimate=True, nominal_value=4.0, lb=0.1, ub=10 + ) + with pytest.raises( + NotImplementedError, match="evaluate to a real number" + ): + assignments_to_experiments(problem) + return + problem.add_parameter("p", estimate=False, nominal_value=4.0) + + converted = assignments_to_experiments(problem) + + # the switch is at the nominal value 4 + 1, not the model value 1 + 1 + assert _changes_by_time(converted, "e1") == {0.0: {}, 5.0: {"k1": 1.0}} From ce6488cd91eadcdc8473d239123cc303d854b687 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Sun, 25 Oct 2026 17:31:50 +0000 Subject: [PATCH 6/8] CHANGELOG: mention switching times set from fixed parameter table values Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a799186..7a67b539 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -64,8 +64,9 @@ still uses elsewhere, evaluate a condition change once at the start of its period, leave time-dependent initial assignments in the model, reject assignments to species or compartments, leave problems without - time-dependent assignments unchanged, and give a ramp's slope parameter an - unused ID + time-dependent assignments unchanged, give a ramp's slope parameter an + unused ID, and accept switching times that conditions set from parameters + fixed by the parameter table (by @FFroehlich in https://github.com/PEtab-dev/libpetab-python/pull/574) **Features** From dbf03328cf055984551da824a11ecea3935d10c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Mon, 26 Oct 2026 07:59:17 +0000 Subject: [PATCH 7/8] Add SbmlModel.get_referenced_ids and ObservableTable.free_symbols `assignments_to_experiments` collected the IDs that the SBML model refers to, and the symbols in the observable and noise formulas, itself. These are now `SbmlModel.get_referenced_ids()` (IDs in any math of the model, and conversion factors) and `ObservableTable.free_symbols`, in line with `ConditionTable.free_symbols`. Co-Authored-By: Claude Opus 5.5 --- petab/v1/models/sbml_model.py | 31 +++++++++++++++++++++++++++++++ petab/v2/converters.py | 32 +++----------------------------- petab/v2/core.py | 16 ++++++++++++++++ tests/v1/test_sbml.py | 19 +++++++++++++++++++ tests/v2/test_core.py | 8 ++++++++ 5 files changed, 77 insertions(+), 29 deletions(-) diff --git a/petab/v1/models/sbml_model.py b/petab/v1/models/sbml_model.py index 030d0ab0..831a42e3 100644 --- a/petab/v1/models/sbml_model.py +++ b/petab/v1/models/sbml_model.py @@ -262,6 +262,37 @@ def get_valid_ids_for_condition_table(self) -> Iterable[str]: def symbol_allowed_in_observable_formula(self, id_: str) -> bool: return self.sbml_model.getElementBySId(id_) or id_ == "time" + def get_referenced_ids(self) -> set[str]: + """Get the IDs that the model refers to. + + These are the IDs in the math of any model element (e.g., kinetic + laws, rules, initial assignments, events), and those referred to as + conversion factors. IDs that are only defined, but not referred to, + are not included. + """ + ids = set() + for element in self.sbml_model.getListOfAllElements(): + if not (hasattr(element, "isSetMath") and element.isSetMath()): + continue + nodes = [element.getMath()] + while nodes: + node = nodes.pop() + if node.getType() == libsbml.AST_NAME: + ids.add(node.getName()) + nodes += [ + node.getChild(i) for i in range(node.getNumChildren()) + ] + # conversion factors refer to parameters by attribute, not in math + ids |= { + element.getConversionFactor() + for element in [ + self.sbml_model, + *self.sbml_model.getListOfSpecies(), + ] + if element.isSetConversionFactor() + } + return ids + def is_valid(self) -> bool: return is_sbml_consistent(self.sbml_model.getSBMLDocument()) diff --git a/petab/v2/converters.py b/petab/v2/converters.py index 2734151a..ea7fef41 100644 --- a/petab/v2/converters.py +++ b/petab/v2/converters.py @@ -1199,39 +1199,13 @@ def _referenced_ids( ) -> set[str]: """The IDs that the model, the observables, or the condition changes, including those in ``additions``, refer to.""" - sbml_model = problem.model.sbml_model - ids = set() - for element in sbml_model.getListOfAllElements(): - if not (hasattr(element, "isSetMath") and element.isSetMath()): - continue - nodes = [element.getMath()] - while nodes: - node = nodes.pop() - if node.getType() == libsbml.AST_NAME: - ids.add(node.getName()) - nodes += [node.getChild(i) for i in range(node.getNumChildren())] - # conversion factors refer to parameters by attribute, not in math - ids |= { - element.getConversionFactor() - for element in [sbml_model, *sbml_model.getListOfSpecies()] - if element.isSetConversionFactor() - } - expressions = [ - *(o.formula for o in problem.observables), - *(o.noise_formula for o in problem.observables), - *(value for *_, value in additions), - ] - ids |= { - str(s) - for expression in expressions - if expression is not None - for s in expression.free_symbols - } + ids = problem.model.get_referenced_ids() ids |= { str(s) - for table in problem.condition_tables + for table in [*problem.observable_tables, *problem.condition_tables] for s in table.free_symbols } + ids |= {str(s) for *_, value in additions for s in value.free_symbols} return ids diff --git a/petab/v2/core.py b/petab/v2/core.py index f69a16ed..b72aa187 100644 --- a/petab/v2/core.py +++ b/petab/v2/core.py @@ -484,6 +484,22 @@ def observables(self) -> list[Observable]: """List of observables.""" return self.elements + @property + def free_symbols(self) -> set[sp.Symbol]: + """Get all free symbols in the observable table. + + This includes all free symbols in the observable and noise formulas, + including the placeholders. + """ + return set( + chain.from_iterable( + formula.free_symbols + for observable in self.observables + for formula in (observable.formula, observable.noise_formula) + if formula is not None + ) + ) + @classmethod def from_df(cls, df: pd.DataFrame, **kwargs) -> ObservableTable: """Create an ObservableTable from a DataFrame.""" diff --git a/tests/v1/test_sbml.py b/tests/v1/test_sbml.py index c3e3e8de..fbc1a5b5 100644 --- a/tests/v1/test_sbml.py +++ b/tests/v1/test_sbml.py @@ -134,6 +134,25 @@ def test_get_condition_specific_models(): check_model(condition_model) +def test_sbml_model_get_referenced_ids(): + """Test that IDs referred to in any math or as a conversion factor are + found, and IDs that are only defined are not.""" + model = SbmlModel.from_antimony(""" + species S1, S2; + S1 = 1; S2 = 0; + k1 = 1; k2 = 2; k3 = 3; cf = 1; unused = 4; + R1: S1 -> S2; k1 * S1; + k4 := k2 * time; + E1: at time > k3: S1 = 0; + """) + model.sbml_model.getSpecies("S2").setConversionFactor("cf") + + ids = model.get_referenced_ids() + + assert {"S1", "k1", "k2", "k3", "cf"} <= ids + assert not {"unused", "k4", "S2"} & ids + + def test_sbml_model_repr(): sbml_document = libsbml.SBMLDocument() sbml_model = sbml_document.createModel() diff --git a/tests/v2/test_core.py b/tests/v2/test_core.py index 5f95ed15..67042bb1 100644 --- a/tests/v2/test_core.py +++ b/tests/v2/test_core.py @@ -482,6 +482,14 @@ def test_condition_table(): ).free_symbols == {x, y} +def test_observable_table_free_symbols(): + assert ObservableTable().free_symbols == set() + + assert ObservableTable( + [Observable(id="obs1", formula=x, noise_formula=y)] + ).free_symbols == {x, y} + + def test_to_df_math_expressions(): """Expressions are written in PEtab syntax and without precision loss. From 3e6327d60987d6720e467ca5778ae7fe97ed5f66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Mon, 26 Oct 2026 08:19:58 +0000 Subject: [PATCH 8/8] Move assignments_to_experiments to its own module `assignments_to_experiments` and its helpers move unchanged from `petab.v2.converters` to `petab.v2.assignments_to_experiments`, and its tests from `tests/v2/test_conversion.py` to `tests/v2/test_assignments_to_experiments.py`. `petab.v2.converters` still exports `assignments_to_experiments` and `DEFAULT_EXPERIMENT_ID`. Co-Authored-By: Claude Opus 5.5 --- doc/modules.rst | 1 + petab/v2/assignments_to_experiments.py | 618 ++++++++++ petab/v2/converters.py | 608 +--------- tests/v2/test_assignments_to_experiments.py | 1207 +++++++++++++++++++ tests/v2/test_conversion.py | 1198 ------------------ 5 files changed, 1836 insertions(+), 1796 deletions(-) create mode 100644 petab/v2/assignments_to_experiments.py create mode 100644 tests/v2/test_assignments_to_experiments.py diff --git a/doc/modules.rst b/doc/modules.rst index 1658aa1d..d461cebc 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.assignments_to_experiments petab.v2.converters petab.v2.core petab.v2.experiments diff --git a/petab/v2/assignments_to_experiments.py b/petab/v2/assignments_to_experiments.py new file mode 100644 index 00000000..8e83fe2a --- /dev/null +++ b/petab/v2/assignments_to_experiments.py @@ -0,0 +1,618 @@ +"""Conversion of time-dependent piecewise assignments to experiments.""" + +from __future__ import annotations + +import warnings +from copy import deepcopy + +import libsbml +import sympy as sp +from sbmlmath import TimeSymbol, sbml_math_to_sympy + +from .core import ( + Change, + Condition, + Experiment, + ExperimentPeriod, + Problem, +) +from .math import petab_symbol +from .models._sbml_utils import add_sbml_parameter +from .models.sbml_model import SbmlModel + +__all__ = ["DEFAULT_EXPERIMENT_ID", "assignments_to_experiments"] + + +DEFAULT_EXPERIMENT_ID = "default_experiment" + + +def assignments_to_experiments(problem: Problem) -> Problem: + """Replace time-dependent piecewise assignments by experiment periods. + + Fixed parameters may change over the course of a single simulation, e.g., + as described by a v2 experiment with multiple periods. As v1 did not have + experiments, then these changes were + sometimes directly encoded in an SBML model as assignment rules with + piecewise functions. This method takes a v2 problem + (possibly converted from a v1 problem) and replaces those SBML assignment + rules with v2 experiments. + + .. note:: + + Measurements without an experiment ID are assigned to a new experiment + :data:`DEFAULT_EXPERIMENT_ID`, if needed, which starts at time 0. + Every experiment keeps the period that starts it, without conditions + if all of them were replaced, since a period's time overrides the + model's initial time. + + **Scope:** conversion is supported for SBML assignment rules to + parameters that depend on time through the conditions of piecewise + functions, and at most linearly through their values, where + + * every such condition is an inequality that can be re-arranged + explicitly into `t />= ...`, + * every piecewise function that depends on time has an explicit + "otherwise" piece, + * every switching time evaluates to a real number in every experiment + (it may depend on parameters that are fixed or set by the condition + table, but not on estimated ones), + * and all other symbols are truly constant model parameters. + + Within that, the piecewise functions may be combined: a sum of doses + given at different times (`d1 * piecewise(0, t < t1, 1) + d2 * + piecewise(0, t < t2, 1)`), phases one after another as several pieces or + nested piecewise functions (`piecewise(b0, t < t1, piecewise(b1, t < t2, + b2))`). Each switch that changes the value becomes a period. A value may + also change linearly with time within a period (a ramp: `piecewise(r * t + + a, t <= t_end, 0)`). A condition change is applied once at the start of + its period, so the target of such an assignment becomes a differential + entity instead, with the rate rule `d target/dt = _slope`, and + each period sets both the target and its slope; an assignment that only + switches between constant values stays a constant parameter. Anything else + raises a :class:`NotImplementedError`. Initial assignments stay in the + model, since they are evaluated once, at the start of the simulation. + + Condition changes of the replaced assignments' symbols are dropped, unless + the model, an observable or another change still refers to them. + + .. warning:: + + As SBML assignment rules, the changes would occur during + pre-equilibration. As v2 experiments, the changes do not occur during + pre-equilibration. Hence, using this converter may affect the + correctness of your model's pre-equilibration. In this case, you could + e.g. keep the SBML assignment rules and/or add indicator variables to + trigger them correctly in the model for some/all pre-equilibration + conditions. + + + :param problem: The PEtab problem to convert. This will not be modified. + :returns: The converted PEtab problem. + :raises NotImplementedError: + If a time-dependent assignment is outside the scope described above. + :raises ValueError: + If the converted problem does not pass linting. + """ + if not isinstance(problem.model, SbmlModel): + raise TypeError("Only SBML models are supported.") + + problem = deepcopy(problem) + sbml_model = problem.model.sbml_model + + # Initial assignments are evaluated once, at the start of the + # simulation, so they never switch and stay in the model. + time_dependent = [] + for rule in sbml_model.getListOfRules(): + if not rule.isAssignment(): + continue + expr = sbml_math_to_sympy(rule) + if any(pw.has(TimeSymbol) for pw in expr.atoms(sp.Piecewise)): + time_dependent.append((rule, expr)) + if not time_dependent: + _check_lint(problem) + return problem + + experiments = list(problem.experiments) + if any(m.experiment_id is None for m in problem.measurements): + default_experiment = Experiment(id=DEFAULT_EXPERIMENT_ID) + experiments.append(default_experiment) + else: + default_experiment = None + droppable = _add_start_periods(experiments) + + # (experiment, period start time, target ID, target value) + additions = [] + replaced_symbols = set() + converted_ids = [] + for rule, expr in time_dependent: + assignment = _PiecewiseAssignment(rule, expr, problem) + assignment.replace_in_model() + for experiment in experiments: + additions += [ + (experiment, *change) + for change in assignment.changes(experiment) + ] + replaced_symbols |= assignment.symbol_ids + converted_ids.append(assignment.target_id) + + if converted_ids and any(e.has_preequilibration for e in experiments): + warnings.warn( + "Switches of time-dependent piecewise assignments during " + "preequilibration are not supported. The converted " + f"assignments to {converted_ids} keep their value at time 0 " + "throughout " + "preequilibration.", + stacklevel=2, + ) + + # replace the changes that are covered by the new conditions, unless + # something else still refers to their targets + replaced_symbols -= _referenced_ids(problem, additions) + _drop_changes(problem, experiments, replaced_symbols) + _add_conditions(problem, additions) + + for experiment in experiments: + experiment.sort_periods() + experiment.periods = [ + p + for p in experiment.periods + if p.condition_ids or id(p) not in droppable + ] + + if default_experiment and any( + p.condition_ids for p in default_experiment.periods + ): + problem += default_experiment + for measurement in problem.measurements: + if measurement.experiment_id is None: + measurement.experiment_id = DEFAULT_EXPERIMENT_ID + + _check_lint(problem) + return problem + + +def _check_lint(problem: Problem) -> None: + """Raise a :class:`ValueError` if the converted ``problem`` does not + pass linting.""" + validation_issues = problem.validate() + if validation_issues.has_errors(): + messages = "\n".join(f"{issue.message}" for issue in validation_issues) + raise ValueError( + f"The converted PEtab problem did not pass linting:\n{messages}" + ) + + +class _PiecewiseAssignment: + """A time-dependent piecewise assignment to a model parameter, and the + condition changes that replace it, see :func:`assignments_to_experiments`. + """ + + def __init__( + self, + rule: libsbml.Rule | libsbml.InitialAssignment, + expr: sp.Expr, + problem: Problem, + ): + """Check that the assignment ``rule``, with the value ``expr``, is + supported.""" + self.rule = rule + self.problem = problem + self.sbml_model = problem.model.sbml_model + self.target_id = rule.getId() + if self.sbml_model.getParameter(self.target_id) is None: + raise NotImplementedError( + "Only assignments to parameters are supported, but " + f"`{self.target_id}` is not a parameter." + ) + self.expr = expr + (self.time,) = expr.atoms(TimeSymbol) + self.symbol_ids = {str(s) for s in expr.free_symbols} + self.parameter_ids = {p.id for p in problem.parameters} + self.switches = _supported_time_switches( + expr, + self.time, + self.target_id, + {p.id for p in problem.parameters if p.estimate}, + ) + if non_constant := sorted( + str(s) + for s in expr.free_symbols - {self.time} + if not (p := self.sbml_model.getParameter(str(s))) + or not p.getConstant() + ): + raise NotImplementedError( + "Only constant model parameters are supported as symbols in " + f"the assignment to `{self.target_id}`, but it refers to " + f"{non_constant}." + ) + # the value before every switch: a condition whose two sides grow + # apart with time holds only after its switch + self.default = expr.xreplace( + { + c: sp.false if slope > 0 else sp.true + for c, (_, slope) in self.switches.items() + } + ) + # unless a condition sets them, the switching times' symbols keep + # their value from the parameter table or else from the model + self.nominal_values = { + p.id: p.nominal_value for p in problem.parameters if not p.estimate + } + self.fixed_values = { + s: self.nominal_values.get(str(s), self._model_value(s)) + for t_switch, _ in self.switches.values() + for s in t_switch.free_symbols + } + # whether the value changes with time within a period + self.ramps = expr.xreplace({c: sp.Dummy() for c in self.switches}).has( + self.time + ) + self.slope_id = None + + def replace_in_model(self) -> None: + """Replace the assignment in the model by a target that the + condition table sets. + + A value that changes with time within a period (a ramp) cannot be a + condition change on a constant target, which is applied once at the + start of its period. The model only makes the target a differential + entity, `d target/dt = _slope`; the slope's value, like the + target's, is set by every period that starts one. + """ + if self.ramps: + self.slope_id = self.problem.get_unused_id( + f"{self.target_id}_slope" + ) + add_sbml_parameter( + self.sbml_model, self.slope_id, value=0.0, constant=False + ) + rate_rule = self.sbml_model.createRateRule() + rate_rule.setVariable(self.target_id) + rate_rule.setMath(libsbml.parseL3Formula(self.slope_id)) + model_values = { + s: self._model_value(s) + for s in self.default.free_symbols - {self.time} + } + target = self.sbml_model.getParameter(self.target_id) + target.setValue( + float(self.default.subs(self.time, 0).subs(model_values)) + ) + target.setConstant(False) + self.rule.removeFromParentAndDelete() + + def changes( + self, experiment: Experiment + ) -> list[tuple[float, str, sp.Expr]]: + """The changes that replace the assignment in ``experiment``, each as + ``(period start time, target ID, target value)``.""" + changes = [] + periods = experiment.sorted_periods + overrides = {} + # the value at the end of the previous period + held = None + for i_period, period in enumerate(periods): + period_changes = self.problem.get_changes_for_period(period) + # a condition's changes are evaluated once, at the start of their + # period, all with the values from before the period + overrides |= { + c.target_id: self._at_period_start(c.target_value, overrides) + for c in period_changes + } + subs = self.fixed_values | { + s: overrides[str(s)] + for s in self.expr.free_symbols + if str(s) in overrides + } + if unknown := sorted( + { + str(x) + for s in self.expr.free_symbols + if str(s) in overrides + for x in overrides[str(s)].free_symbols + } + - self.parameter_ids + ): + raise NotImplementedError( + "Only condition changes whose values are known at the " + "start of their period are supported, but in experiment " + f"`{experiment.id}`, the changes of the symbols of the " + f"assignment to `{self.target_id}` refer to {unknown}." + ) + times = self._switching_times(subs, experiment.id) + + if period.is_preequilibration: + # the value at time 0, which respects whether a condition + # holds for a switch at time 0, held there: a ramp does not + # run during preequilibration + held = self.expr.subs(subs).subs(self.time, 0) + changes.append( + ( + period.time, + self.target_id, + self._with_model_values(held), + ) + ) + if self.ramps: + changes.append((period.time, self.slope_id, sp.S.Zero)) + continue + + start_value = self._value_from(subs, times, period.time) + # each switch within the period that changes the value + t_end = ( + periods[i_period + 1].time + if i_period + 1 < len(periods) + else sp.oo + ) + switch_values = [] + previous = start_value + for t in sorted({t for t in times if period.time < t < t_end}): + value = self._value_from(subs, times, t) + if not (value - previous).is_zero: + switch_values.append((t, value)) + previous = value + + if i_period == 0: + # The model already has the pre-switch value at the start of + # the first period, which stays even without conditions, + # unless that value refers to parameter table symbols, for + # which the model only has a number. A ramp's start and + # slope are always set by the tables. + if ( + self.ramps + or not (start_value - self.default).is_zero + or {str(s) for s in self.default.free_symbols} + & self.parameter_ids + ): + changes += self._start( + period.time, self._with_model_values(start_value) + ) + elif ( + {c.target_id for c in period_changes} & self.symbol_ids + # preequilibration held the slope at zero + or (self.ramps and periods[i_period - 1].is_preequilibration) + # a switch at the start of the period + or not (start_value - held).is_zero + ): + changes += self._start(period.time, start_value) + for t, value in switch_values: + changes += self._start(t, value) + held = previous + return changes + + def _model_value(self, symbol: sp.Symbol) -> float: + """The value of the model parameter ``symbol``.""" + return self.sbml_model.getParameter(str(symbol)).getValue() + + def _at_period_start(self, value: sp.Expr, overrides: dict) -> sp.Expr: + """The ``value`` of a change, evaluated at the start of its period. + + Symbols that a condition set before take that value (``overrides``), + parameters fixed by the parameter table their nominal values, and the + other constant model parameters outside the parameter table their + model values. Estimated parameters stay, since no condition sets them; + any other symbol stays as well, since its value is not known here. + """ + if value is None: + return value + known = {} + for s in value.free_symbols: + if str(s) in overrides: + known[s] = overrides[str(s)] + elif str(s) in self.nominal_values: + known[s] = sp.Float(self.nominal_values[str(s)]) + elif ( + str(s) not in self.parameter_ids + and (p := self.sbml_model.getParameter(str(s))) is not None + and p.getConstant() + ): + known[s] = sp.Float(p.getValue()) + return value.xreplace(known) + + def _with_model_values(self, value: sp.Expr) -> sp.Expr: + """``value`` for a change at the start of an experiment, which may + only refer to parameter table symbols: the other symbols are constant + model parameters, at their values.""" + return value.subs( + { + s: self._model_value(s) + for s in value.free_symbols - {self.time} + if str(s) not in self.parameter_ids + } + ) + + def _switching_times(self, subs: dict, experiment_id: str) -> list[float]: + """The switching times, for the symbol values ``subs``.""" + times = [] + for t_switch, _ in self.switches.values(): + t = t_switch.subs(subs) + if not (t.is_number and t.is_real and t.is_finite): + raise NotImplementedError( + "Only switching times that evaluate to a real number are " + "supported, but a switching time of the assignment to " + f"`{self.target_id}` in experiment `{experiment_id}` is " + f"`{t}`." + ) + times.append(float(t)) + return times + + def _value_from(self, subs: dict, times: list[float], t0: float): + """The value from ``t0`` up to its next switch, for the symbol values + ``subs``, every condition replaced by whether it holds in between.""" + later = [t for t in times if t > t0] + probe = (t0 + min(later)) / 2 if later else t0 + 1 + truth = {c: c.subs(subs).subs(self.time, probe) for c in self.switches} + return self.expr.xreplace(truth).subs(subs) + + def _start( + self, t: float, value: sp.Expr + ) -> list[tuple[float, str, sp.Expr]]: + """The changes that start ``value`` at time ``t``: the target set to + the value at ``t``, and, for a target that ramps, its slope set to + the value's rate of change.""" + changes = [(t, self.target_id, value.subs(self.time, t))] + if self.ramps: + changes.append((t, self.slope_id, sp.diff(value, self.time))) + return changes + + +def _add_start_periods(experiments: list[Experiment]) -> set[int]: + """Give every experiment a period that starts it, and return the + :func:`id` of the periods that are dropped if all their conditions are + replaced. + + The period that starts an experiment stays, even without conditions: its + time overrides the model's initial time, so the experiment would start at + its first switch instead. + """ + droppable = set() + for experiment in experiments: + if all(p.is_preequilibration for p in experiment.periods): + experiment += ExperimentPeriod(time=0) + start = min( + (p for p in experiment.periods if not p.is_preequilibration), + key=lambda p: p.time, + ) + droppable |= { + id(p) + for p in experiment.periods + if p.condition_ids and p is not start + } + return droppable + + +def _referenced_ids( + problem: Problem, + additions: list[tuple[Experiment, float, str, sp.Expr]], +) -> set[str]: + """The IDs that the model, the observables, or the condition changes, + including those in ``additions``, refer to.""" + ids = problem.model.get_referenced_ids() + ids |= { + str(s) + for table in [*problem.observable_tables, *problem.condition_tables] + for s in table.free_symbols + } + ids |= {str(s) for *_, value in additions for s in value.free_symbols} + return ids + + +def _drop_changes( + problem: Problem, experiments: list[Experiment], target_ids: set[str] +) -> None: + """Drop the changes of ``target_ids`` from the condition tables, and the + conditions left without changes from the periods of ``experiments``.""" + condition_ids = {c.id for c in problem.conditions} + for table in problem.condition_tables: + for condition in table.conditions: + condition.changes = [ + c for c in condition.changes if c.target_id not in target_ids + ] + table.elements = [c for c in table.conditions if c.changes] + removed_ids = condition_ids - {c.id for c in problem.conditions} + for experiment in experiments: + for period in experiment.periods: + period.condition_ids = [ + cid for cid in period.condition_ids if cid not in removed_ids + ] + + +def _add_conditions( + problem: Problem, + additions: list[tuple[Experiment, float, str, sp.Expr]], +) -> None: + """For each ``(experiment, time, target ID, value)`` in ``additions``, + add a condition with that change to the period of ``experiment`` that + starts at ``time``, which is added if needed. + + Equal changes share a condition. Their values use the symbols that PEtab + uses for the condition table, rather than those parsed from SBML. + """ + new_conditions = {} + for experiment, t, target_id, value in additions: + value = value.xreplace( + { + s: petab_symbol(s.name) + for s in value.free_symbols + if not isinstance(s, TimeSymbol) + } + ) + if (cond_id := new_conditions.get((target_id, value))) is None: + cond_id = problem.get_unused_id(f"cond_{len(new_conditions)}") + new_conditions[target_id, value] = cond_id + problem += Condition( + id=cond_id, + changes=[Change(target_id=target_id, target_value=value)], + ) + period = next((p for p in experiment.periods if p.time == t), None) + if period is None: + experiment += (period := ExperimentPeriod(time=t)) + period.condition_ids.append(cond_id) + + +def _supported_time_switches( + expr: sp.Expr, time: sp.Symbol, target_id: str, estimated_ids: set[str] +) -> dict[sp.Rel, tuple[sp.Expr, sp.Expr]]: + """Check that ``expr`` depends on ``time`` only in the supported way, and + return the conditions through which it does, each with the time at which + it switches and the slope of its two sides' difference. + + ``expr`` may combine several piecewise functions, each with several + pieces and nested in one another (a sum of doses, phases one after + another), as long as time enters only through their conditions, each a + single inequality linear in time, each piecewise function has an + explicit "otherwise" piece, and no switching time depends on the + parameters ``estimated_ids``. Anything else raises + :class:`NotImplementedError`. + """ + unsupported = NotImplementedError( + "Only time dependence through piecewise conditions, each a single " + "inequality linear in time, with an explicit otherwise piece, is " + f"supported, but the assignment to `{target_id}` is `{expr}`." + ) + switches = {} + for pw in expr.atoms(sp.Piecewise): + if not pw.has(time): + continue + # an explicit "otherwise" piece + if pw.args[-1].cond != sp.true: + raise unsupported + for piece in pw.args[:-1]: + condition = piece.cond + if not condition.has(time): + continue + if not ( + # a supported inequality + isinstance(condition, (sp.Gt, sp.Ge, sp.Lt, sp.Le)) + # that changes value at maximally a single time point + and ( + slope := sp.diff(condition.gts - condition.lts, time) + ).is_number + # and this time point is identifiable + and slope.is_nonzero + ): + raise unsupported + t_switch = sp.solve(condition.gts - condition.lts, time)[0] + switches[condition] = (t_switch, slope) + # outside the conditions, time enters at most linearly: the value may + # ramp, at a rate that does not itself change with time + rest = expr.xreplace({c: sp.Dummy() for c in switches}) + if rest.has(time) and sp.diff(rest, time).has(time): + raise NotImplementedError( + "Only values that change at most linearly with time are " + f"supported, but the assignment to `{target_id}` is `{expr}`." + ) + # no period can start at an estimated time + if estimated := sorted( + { + str(s) + for t_switch, _ in switches.values() + for s in t_switch.free_symbols + if str(s) in estimated_ids + } + ): + raise NotImplementedError( + "Only switching times that are not estimated are supported, but " + f"the switching times of the assignment to `{target_id}` depend " + f"on the estimated parameters {estimated}." + ) + return switches diff --git a/petab/v2/converters.py b/petab/v2/converters.py index ea7fef41..6d7d2d55 100644 --- a/petab/v2/converters.py +++ b/petab/v2/converters.py @@ -9,8 +9,12 @@ import libsbml import sympy as sp -from sbmlmath import TimeSymbol, sbml_math_to_sympy, set_math +from sbmlmath import sbml_math_to_sympy, set_math +from .assignments_to_experiments import ( + DEFAULT_EXPERIMENT_ID, + assignments_to_experiments, +) from .core import ( Change, Condition, @@ -19,7 +23,6 @@ ExperimentPeriod, Problem, ) -from .math import petab_symbol from .models._sbml_utils import ( add_sbml_parameter, check, @@ -27,7 +30,11 @@ ) from .models.sbml_model import SbmlModel -__all__ = ["ExperimentsToSbmlConverter", "assignments_to_experiments"] +__all__ = [ + "DEFAULT_EXPERIMENT_ID", + "ExperimentsToSbmlConverter", + "assignments_to_experiments", +] logger = logging.getLogger(__name__) @@ -734,598 +741,3 @@ def _add_indicators_to_conditions(self) -> None: ] ) ] - - -DEFAULT_EXPERIMENT_ID = "default_experiment" - - -def assignments_to_experiments(problem: Problem) -> Problem: - """Replace time-dependent piecewise assignments by experiment periods. - - Fixed parameters may change over the course of a single simulation, e.g., - as described by a v2 experiment with multiple periods. As v1 did not have - experiments, then these changes were - sometimes directly encoded in an SBML model as assignment rules with - piecewise functions. This method takes a v2 problem - (possibly converted from a v1 problem) and replaces those SBML assignment - rules with v2 experiments. - - .. note:: - - Measurements without an experiment ID are assigned to a new experiment - :data:`DEFAULT_EXPERIMENT_ID`, if needed, which starts at time 0. - Every experiment keeps the period that starts it, without conditions - if all of them were replaced, since a period's time overrides the - model's initial time. - - **Scope:** conversion is supported for SBML assignment rules to - parameters that depend on time through the conditions of piecewise - functions, and at most linearly through their values, where - - * every such condition is an inequality that can be re-arranged - explicitly into `t />= ...`, - * every piecewise function that depends on time has an explicit - "otherwise" piece, - * every switching time evaluates to a real number in every experiment - (it may depend on parameters that are fixed or set by the condition - table, but not on estimated ones), - * and all other symbols are truly constant model parameters. - - Within that, the piecewise functions may be combined: a sum of doses - given at different times (`d1 * piecewise(0, t < t1, 1) + d2 * - piecewise(0, t < t2, 1)`), phases one after another as several pieces or - nested piecewise functions (`piecewise(b0, t < t1, piecewise(b1, t < t2, - b2))`). Each switch that changes the value becomes a period. A value may - also change linearly with time within a period (a ramp: `piecewise(r * t - + a, t <= t_end, 0)`). A condition change is applied once at the start of - its period, so the target of such an assignment becomes a differential - entity instead, with the rate rule `d target/dt = _slope`, and - each period sets both the target and its slope; an assignment that only - switches between constant values stays a constant parameter. Anything else - raises a :class:`NotImplementedError`. Initial assignments stay in the - model, since they are evaluated once, at the start of the simulation. - - Condition changes of the replaced assignments' symbols are dropped, unless - the model, an observable or another change still refers to them. - - .. warning:: - - As SBML assignment rules, the changes would occur during - pre-equilibration. As v2 experiments, the changes do not occur during - pre-equilibration. Hence, using this converter may affect the - correctness of your model's pre-equilibration. In this case, you could - e.g. keep the SBML assignment rules and/or add indicator variables to - trigger them correctly in the model for some/all pre-equilibration - conditions. - - - :param problem: The PEtab problem to convert. This will not be modified. - :returns: The converted PEtab problem. - :raises NotImplementedError: - If a time-dependent assignment is outside the scope described above. - :raises ValueError: - If the converted problem does not pass linting. - """ - if not isinstance(problem.model, SbmlModel): - raise TypeError("Only SBML models are supported.") - - problem = deepcopy(problem) - sbml_model = problem.model.sbml_model - - # Initial assignments are evaluated once, at the start of the - # simulation, so they never switch and stay in the model. - time_dependent = [] - for rule in sbml_model.getListOfRules(): - if not rule.isAssignment(): - continue - expr = sbml_math_to_sympy(rule) - if any(pw.has(TimeSymbol) for pw in expr.atoms(sp.Piecewise)): - time_dependent.append((rule, expr)) - if not time_dependent: - _check_lint(problem) - return problem - - experiments = list(problem.experiments) - if any(m.experiment_id is None for m in problem.measurements): - default_experiment = Experiment(id=DEFAULT_EXPERIMENT_ID) - experiments.append(default_experiment) - else: - default_experiment = None - droppable = _add_start_periods(experiments) - - # (experiment, period start time, target ID, target value) - additions = [] - replaced_symbols = set() - converted_ids = [] - for rule, expr in time_dependent: - assignment = _PiecewiseAssignment(rule, expr, problem) - assignment.replace_in_model() - for experiment in experiments: - additions += [ - (experiment, *change) - for change in assignment.changes(experiment) - ] - replaced_symbols |= assignment.symbol_ids - converted_ids.append(assignment.target_id) - - if converted_ids and any(e.has_preequilibration for e in experiments): - warnings.warn( - "Switches of time-dependent piecewise assignments during " - "preequilibration are not supported. The converted " - f"assignments to {converted_ids} keep their value at time 0 " - "throughout " - "preequilibration.", - stacklevel=2, - ) - - # replace the changes that are covered by the new conditions, unless - # something else still refers to their targets - replaced_symbols -= _referenced_ids(problem, additions) - _drop_changes(problem, experiments, replaced_symbols) - _add_conditions(problem, additions) - - for experiment in experiments: - experiment.sort_periods() - experiment.periods = [ - p - for p in experiment.periods - if p.condition_ids or id(p) not in droppable - ] - - if default_experiment and any( - p.condition_ids for p in default_experiment.periods - ): - problem += default_experiment - for measurement in problem.measurements: - if measurement.experiment_id is None: - measurement.experiment_id = DEFAULT_EXPERIMENT_ID - - _check_lint(problem) - return problem - - -def _check_lint(problem: Problem) -> None: - """Raise a :class:`ValueError` if the converted ``problem`` does not - pass linting.""" - validation_issues = problem.validate() - if validation_issues.has_errors(): - messages = "\n".join(f"{issue.message}" for issue in validation_issues) - raise ValueError( - f"The converted PEtab problem did not pass linting:\n{messages}" - ) - - -class _PiecewiseAssignment: - """A time-dependent piecewise assignment to a model parameter, and the - condition changes that replace it, see :func:`assignments_to_experiments`. - """ - - def __init__( - self, - rule: libsbml.Rule | libsbml.InitialAssignment, - expr: sp.Expr, - problem: Problem, - ): - """Check that the assignment ``rule``, with the value ``expr``, is - supported.""" - self.rule = rule - self.problem = problem - self.sbml_model = problem.model.sbml_model - self.target_id = rule.getId() - if self.sbml_model.getParameter(self.target_id) is None: - raise NotImplementedError( - "Only assignments to parameters are supported, but " - f"`{self.target_id}` is not a parameter." - ) - self.expr = expr - (self.time,) = expr.atoms(TimeSymbol) - self.symbol_ids = {str(s) for s in expr.free_symbols} - self.parameter_ids = {p.id for p in problem.parameters} - self.switches = _supported_time_switches( - expr, - self.time, - self.target_id, - {p.id for p in problem.parameters if p.estimate}, - ) - if non_constant := sorted( - str(s) - for s in expr.free_symbols - {self.time} - if not (p := self.sbml_model.getParameter(str(s))) - or not p.getConstant() - ): - raise NotImplementedError( - "Only constant model parameters are supported as symbols in " - f"the assignment to `{self.target_id}`, but it refers to " - f"{non_constant}." - ) - # the value before every switch: a condition whose two sides grow - # apart with time holds only after its switch - self.default = expr.xreplace( - { - c: sp.false if slope > 0 else sp.true - for c, (_, slope) in self.switches.items() - } - ) - # unless a condition sets them, the switching times' symbols keep - # their value from the parameter table or else from the model - self.nominal_values = { - p.id: p.nominal_value for p in problem.parameters if not p.estimate - } - self.fixed_values = { - s: self.nominal_values.get(str(s), self._model_value(s)) - for t_switch, _ in self.switches.values() - for s in t_switch.free_symbols - } - # whether the value changes with time within a period - self.ramps = expr.xreplace({c: sp.Dummy() for c in self.switches}).has( - self.time - ) - self.slope_id = None - - def replace_in_model(self) -> None: - """Replace the assignment in the model by a target that the - condition table sets. - - A value that changes with time within a period (a ramp) cannot be a - condition change on a constant target, which is applied once at the - start of its period. The model only makes the target a differential - entity, `d target/dt = _slope`; the slope's value, like the - target's, is set by every period that starts one. - """ - if self.ramps: - self.slope_id = self.problem.get_unused_id( - f"{self.target_id}_slope" - ) - add_sbml_parameter( - self.sbml_model, self.slope_id, value=0.0, constant=False - ) - rate_rule = self.sbml_model.createRateRule() - rate_rule.setVariable(self.target_id) - rate_rule.setMath(libsbml.parseL3Formula(self.slope_id)) - model_values = { - s: self._model_value(s) - for s in self.default.free_symbols - {self.time} - } - target = self.sbml_model.getParameter(self.target_id) - target.setValue( - float(self.default.subs(self.time, 0).subs(model_values)) - ) - target.setConstant(False) - self.rule.removeFromParentAndDelete() - - def changes( - self, experiment: Experiment - ) -> list[tuple[float, str, sp.Expr]]: - """The changes that replace the assignment in ``experiment``, each as - ``(period start time, target ID, target value)``.""" - changes = [] - periods = experiment.sorted_periods - overrides = {} - # the value at the end of the previous period - held = None - for i_period, period in enumerate(periods): - period_changes = self.problem.get_changes_for_period(period) - # a condition's changes are evaluated once, at the start of their - # period, all with the values from before the period - overrides |= { - c.target_id: self._at_period_start(c.target_value, overrides) - for c in period_changes - } - subs = self.fixed_values | { - s: overrides[str(s)] - for s in self.expr.free_symbols - if str(s) in overrides - } - if unknown := sorted( - { - str(x) - for s in self.expr.free_symbols - if str(s) in overrides - for x in overrides[str(s)].free_symbols - } - - self.parameter_ids - ): - raise NotImplementedError( - "Only condition changes whose values are known at the " - "start of their period are supported, but in experiment " - f"`{experiment.id}`, the changes of the symbols of the " - f"assignment to `{self.target_id}` refer to {unknown}." - ) - times = self._switching_times(subs, experiment.id) - - if period.is_preequilibration: - # the value at time 0, which respects whether a condition - # holds for a switch at time 0, held there: a ramp does not - # run during preequilibration - held = self.expr.subs(subs).subs(self.time, 0) - changes.append( - ( - period.time, - self.target_id, - self._with_model_values(held), - ) - ) - if self.ramps: - changes.append((period.time, self.slope_id, sp.S.Zero)) - continue - - start_value = self._value_from(subs, times, period.time) - # each switch within the period that changes the value - t_end = ( - periods[i_period + 1].time - if i_period + 1 < len(periods) - else sp.oo - ) - switch_values = [] - previous = start_value - for t in sorted({t for t in times if period.time < t < t_end}): - value = self._value_from(subs, times, t) - if not (value - previous).is_zero: - switch_values.append((t, value)) - previous = value - - if i_period == 0: - # The model already has the pre-switch value at the start of - # the first period, which stays even without conditions, - # unless that value refers to parameter table symbols, for - # which the model only has a number. A ramp's start and - # slope are always set by the tables. - if ( - self.ramps - or not (start_value - self.default).is_zero - or {str(s) for s in self.default.free_symbols} - & self.parameter_ids - ): - changes += self._start( - period.time, self._with_model_values(start_value) - ) - elif ( - {c.target_id for c in period_changes} & self.symbol_ids - # preequilibration held the slope at zero - or (self.ramps and periods[i_period - 1].is_preequilibration) - # a switch at the start of the period - or not (start_value - held).is_zero - ): - changes += self._start(period.time, start_value) - for t, value in switch_values: - changes += self._start(t, value) - held = previous - return changes - - def _model_value(self, symbol: sp.Symbol) -> float: - """The value of the model parameter ``symbol``.""" - return self.sbml_model.getParameter(str(symbol)).getValue() - - def _at_period_start(self, value: sp.Expr, overrides: dict) -> sp.Expr: - """The ``value`` of a change, evaluated at the start of its period. - - Symbols that a condition set before take that value (``overrides``), - parameters fixed by the parameter table their nominal values, and the - other constant model parameters outside the parameter table their - model values. Estimated parameters stay, since no condition sets them; - any other symbol stays as well, since its value is not known here. - """ - if value is None: - return value - known = {} - for s in value.free_symbols: - if str(s) in overrides: - known[s] = overrides[str(s)] - elif str(s) in self.nominal_values: - known[s] = sp.Float(self.nominal_values[str(s)]) - elif ( - str(s) not in self.parameter_ids - and (p := self.sbml_model.getParameter(str(s))) is not None - and p.getConstant() - ): - known[s] = sp.Float(p.getValue()) - return value.xreplace(known) - - def _with_model_values(self, value: sp.Expr) -> sp.Expr: - """``value`` for a change at the start of an experiment, which may - only refer to parameter table symbols: the other symbols are constant - model parameters, at their values.""" - return value.subs( - { - s: self._model_value(s) - for s in value.free_symbols - {self.time} - if str(s) not in self.parameter_ids - } - ) - - def _switching_times(self, subs: dict, experiment_id: str) -> list[float]: - """The switching times, for the symbol values ``subs``.""" - times = [] - for t_switch, _ in self.switches.values(): - t = t_switch.subs(subs) - if not (t.is_number and t.is_real and t.is_finite): - raise NotImplementedError( - "Only switching times that evaluate to a real number are " - "supported, but a switching time of the assignment to " - f"`{self.target_id}` in experiment `{experiment_id}` is " - f"`{t}`." - ) - times.append(float(t)) - return times - - def _value_from(self, subs: dict, times: list[float], t0: float): - """The value from ``t0`` up to its next switch, for the symbol values - ``subs``, every condition replaced by whether it holds in between.""" - later = [t for t in times if t > t0] - probe = (t0 + min(later)) / 2 if later else t0 + 1 - truth = {c: c.subs(subs).subs(self.time, probe) for c in self.switches} - return self.expr.xreplace(truth).subs(subs) - - def _start( - self, t: float, value: sp.Expr - ) -> list[tuple[float, str, sp.Expr]]: - """The changes that start ``value`` at time ``t``: the target set to - the value at ``t``, and, for a target that ramps, its slope set to - the value's rate of change.""" - changes = [(t, self.target_id, value.subs(self.time, t))] - if self.ramps: - changes.append((t, self.slope_id, sp.diff(value, self.time))) - return changes - - -def _add_start_periods(experiments: list[Experiment]) -> set[int]: - """Give every experiment a period that starts it, and return the - :func:`id` of the periods that are dropped if all their conditions are - replaced. - - The period that starts an experiment stays, even without conditions: its - time overrides the model's initial time, so the experiment would start at - its first switch instead. - """ - droppable = set() - for experiment in experiments: - if all(p.is_preequilibration for p in experiment.periods): - experiment += ExperimentPeriod(time=0) - start = min( - (p for p in experiment.periods if not p.is_preequilibration), - key=lambda p: p.time, - ) - droppable |= { - id(p) - for p in experiment.periods - if p.condition_ids and p is not start - } - return droppable - - -def _referenced_ids( - problem: Problem, - additions: list[tuple[Experiment, float, str, sp.Expr]], -) -> set[str]: - """The IDs that the model, the observables, or the condition changes, - including those in ``additions``, refer to.""" - ids = problem.model.get_referenced_ids() - ids |= { - str(s) - for table in [*problem.observable_tables, *problem.condition_tables] - for s in table.free_symbols - } - ids |= {str(s) for *_, value in additions for s in value.free_symbols} - return ids - - -def _drop_changes( - problem: Problem, experiments: list[Experiment], target_ids: set[str] -) -> None: - """Drop the changes of ``target_ids`` from the condition tables, and the - conditions left without changes from the periods of ``experiments``.""" - condition_ids = {c.id for c in problem.conditions} - for table in problem.condition_tables: - for condition in table.conditions: - condition.changes = [ - c for c in condition.changes if c.target_id not in target_ids - ] - table.elements = [c for c in table.conditions if c.changes] - removed_ids = condition_ids - {c.id for c in problem.conditions} - for experiment in experiments: - for period in experiment.periods: - period.condition_ids = [ - cid for cid in period.condition_ids if cid not in removed_ids - ] - - -def _add_conditions( - problem: Problem, - additions: list[tuple[Experiment, float, str, sp.Expr]], -) -> None: - """For each ``(experiment, time, target ID, value)`` in ``additions``, - add a condition with that change to the period of ``experiment`` that - starts at ``time``, which is added if needed. - - Equal changes share a condition. Their values use the symbols that PEtab - uses for the condition table, rather than those parsed from SBML. - """ - new_conditions = {} - for experiment, t, target_id, value in additions: - value = value.xreplace( - { - s: petab_symbol(s.name) - for s in value.free_symbols - if not isinstance(s, TimeSymbol) - } - ) - if (cond_id := new_conditions.get((target_id, value))) is None: - cond_id = problem.get_unused_id(f"cond_{len(new_conditions)}") - new_conditions[target_id, value] = cond_id - problem += Condition( - id=cond_id, - changes=[Change(target_id=target_id, target_value=value)], - ) - period = next((p for p in experiment.periods if p.time == t), None) - if period is None: - experiment += (period := ExperimentPeriod(time=t)) - period.condition_ids.append(cond_id) - - -def _supported_time_switches( - expr: sp.Expr, time: sp.Symbol, target_id: str, estimated_ids: set[str] -) -> dict[sp.Rel, tuple[sp.Expr, sp.Expr]]: - """Check that ``expr`` depends on ``time`` only in the supported way, and - return the conditions through which it does, each with the time at which - it switches and the slope of its two sides' difference. - - ``expr`` may combine several piecewise functions, each with several - pieces and nested in one another (a sum of doses, phases one after - another), as long as time enters only through their conditions, each a - single inequality linear in time, each piecewise function has an - explicit "otherwise" piece, and no switching time depends on the - parameters ``estimated_ids``. Anything else raises - :class:`NotImplementedError`. - """ - unsupported = NotImplementedError( - "Only time dependence through piecewise conditions, each a single " - "inequality linear in time, with an explicit otherwise piece, is " - f"supported, but the assignment to `{target_id}` is `{expr}`." - ) - switches = {} - for pw in expr.atoms(sp.Piecewise): - if not pw.has(time): - continue - # an explicit "otherwise" piece - if pw.args[-1].cond != sp.true: - raise unsupported - for piece in pw.args[:-1]: - condition = piece.cond - if not condition.has(time): - continue - if not ( - # a supported inequality - isinstance(condition, (sp.Gt, sp.Ge, sp.Lt, sp.Le)) - # that changes value at maximally a single time point - and ( - slope := sp.diff(condition.gts - condition.lts, time) - ).is_number - # and this time point is identifiable - and slope.is_nonzero - ): - raise unsupported - t_switch = sp.solve(condition.gts - condition.lts, time)[0] - switches[condition] = (t_switch, slope) - # outside the conditions, time enters at most linearly: the value may - # ramp, at a rate that does not itself change with time - rest = expr.xreplace({c: sp.Dummy() for c in switches}) - if rest.has(time) and sp.diff(rest, time).has(time): - raise NotImplementedError( - "Only values that change at most linearly with time are " - f"supported, but the assignment to `{target_id}` is `{expr}`." - ) - # no period can start at an estimated time - if estimated := sorted( - { - str(s) - for t_switch, _ in switches.values() - for s in t_switch.free_symbols - if str(s) in estimated_ids - } - ): - raise NotImplementedError( - "Only switching times that are not estimated are supported, but " - f"the switching times of the assignment to `{target_id}` depend " - f"on the estimated parameters {estimated}." - ) - return switches diff --git a/tests/v2/test_assignments_to_experiments.py b/tests/v2/test_assignments_to_experiments.py new file mode 100644 index 00000000..54bc667b --- /dev/null +++ b/tests/v2/test_assignments_to_experiments.py @@ -0,0 +1,1207 @@ +"""Tests for `petab.v2.converters.assignments_to_experiments`.""" + +import libsbml +import numpy as np +import pandas as pd +import pytest + +from petab import v2 +from petab.v1.models.sbml_model import SbmlModel +from petab.v2 import Problem +from petab.v2.converters import ( + DEFAULT_EXPERIMENT_ID, + ExperimentsToSbmlConverter, + assignments_to_experiments, +) +from petab.v2.math import sympify_petab + + +@pytest.mark.parametrize( + "condition, before, after", + [("time < 10", 1.0, 2.0), ("10 < time", 2.0, 1.0)], +) +def test_assignments_to_experiments(condition, before, after): + """Test that switches in either direction are converted.""" + problem = Problem() + ant_model = f""" + model conversion + species A, B; + A = 10; + B = 0; + k1 := piecewise(1, {condition}, 2); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=0.0, measurement=10.0) + problem.add_measurement("obs_a", time=20.0, measurement=5.0) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + # the value before the switch is kept in the model, so that only the + # switch itself needs a condition; the experiment still starts at 0, + # where the first measurement is + assert converted.model.sbml_model.getParameter("k1").getValue() == before + + experiments_expected = pd.DataFrame( + data={ + v2.C.EXPERIMENT_ID: [DEFAULT_EXPERIMENT_ID] * 2, + v2.C.TIME: [0.0, 10.0], + v2.C.CONDITION_ID: ["", "cond_0"], + } + ) + + conditions_expected = pd.DataFrame( + data={ + v2.C.CONDITION_ID: ["cond_0"], + v2.C.TARGET_ID: ["k1"], + v2.C.TARGET_VALUE: [after], + } + ) + + assert converted.experiment_df.equals(experiments_expected) + assert converted.condition_df.equals(conditions_expected) + assert converted.measurement_df[v2.C.EXPERIMENT_ID].tolist() == [ + DEFAULT_EXPERIMENT_ID, + DEFAULT_EXPERIMENT_ID, + ] + + +def test_assignments_to_experiments_unused_condition_id(): + """Test that new condition IDs don't clash with other IDs.""" + problem = Problem() + problem.model = SbmlModel.from_antimony(""" + species A = 10; + k1 := piecewise(1, time < 10, 2); + R1: A -> ; k1 * A; + """) + problem.add_observable("obs_a", formula="A", noise_formula="cond_0") + problem.add_measurement("obs_a", time=20.0, measurement=5.0) + problem.add_parameter("cond_0", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + assert [c.id for c in converted.conditions] == ["cond_0_1"] + # the experiment starts at 0, where nothing changes + assert [ + (p.time, p.condition_ids) for p in converted.experiments[0].periods + ] == [(0.0, []), (10.0, ["cond_0_1"])] + + +def test_assignments_to_experiments_nested(): + """Test conversion of piecewise assignments that are nested inside a + larger expression and switch at a time given by some expression.""" + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + level = 3; + k1 := level * piecewise(0, time - 10 < 0, 1); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=0.0, measurement=10.0) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + assert converted.model.sbml_model.getParameter("k1").getValue() == 0.0 + + # `level` is not a parameter table entry, so it must not appear in the + # first period -- the switch may refer to it, since the experiment + # starts with an empty period at 0 + experiments_expected = pd.DataFrame( + data={ + v2.C.EXPERIMENT_ID: [DEFAULT_EXPERIMENT_ID] * 2, + v2.C.TIME: [0.0, 10.0], + v2.C.CONDITION_ID: ["", "cond_0"], + } + ) + + # the factor outside the piecewise expression is preserved + conditions_expected = pd.DataFrame( + data={ + v2.C.CONDITION_ID: ["cond_0"], + v2.C.TARGET_ID: ["k1"], + v2.C.TARGET_VALUE: ["level"], + } + ) + + assert converted.experiment_df.equals(experiments_expected) + assert converted.condition_df.equals(conditions_expected) + + +def test_assignments_to_experiments_parameter_table_default(): + """Test that a pre-switch value that refers to parameter table symbols is + set at the start, since the model only has a number for it.""" + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + d1 = 1; + d2 = 2; + k1 := piecewise(d1, time < 10, d2); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=5.0) + problem.add_parameter( + "d1", estimate=True, lb=1e-3, ub=1e3, nominal_value=3 + ) + problem.add_parameter( + "d2", estimate=True, lb=1e-3, ub=1e3, nominal_value=4 + ) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + experiments_expected = pd.DataFrame( + data={ + v2.C.EXPERIMENT_ID: [DEFAULT_EXPERIMENT_ID] * 2, + v2.C.TIME: [0.0, 10.0], + v2.C.CONDITION_ID: ["cond_0", "cond_1"], + } + ) + + conditions_expected = pd.DataFrame( + data={ + v2.C.CONDITION_ID: ["cond_0", "cond_1"], + v2.C.TARGET_ID: ["k1", "k1"], + v2.C.TARGET_VALUE: ["d1", "d2"], + } + ) + + assert converted.experiment_df.equals(experiments_expected) + assert converted.condition_df.equals(conditions_expected) + + +def test_assignments_to_experiments_condition_specific(): + """Test conversion of piecewise assignments whose switching times and + values are set by the condition table.""" + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + dose = 0; + dose_time = 0; + k1 := dose * piecewise(0, time - dose_time < 0, 1); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + problem.add_condition("c1", dose=2.0, dose_time=10.0) + problem.add_condition("c2", dose=3.0, dose_time=20.0) + problem.add_experiment("experiment____c1", 0, "c1") + problem.add_experiment("experiment____c2", 0, "c2") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="experiment____c1", time=30.0, measurement=10.0 + ) + problem.add_measurement( + "obs_a", experiment_id="experiment____c2", time=30.0, measurement=5.0 + ) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + # the pre-dose value is the same for both conditions and is kept in the + # model; each condition gets its own experiment and switching time, and + # each experiment still starts at 0, where its condition is now empty + assert converted.model.sbml_model.getParameter("k1").getValue() == 0.0 + + experiments_expected = pd.DataFrame( + data={ + v2.C.EXPERIMENT_ID: ["experiment____c1"] * 2 + + ["experiment____c2"] * 2, + v2.C.TIME: [0.0, 10.0, 0.0, 20.0], + v2.C.CONDITION_ID: ["", "cond_0", "", "cond_1"], + } + ) + + # `dose` and `dose_time` are covered by the new conditions + conditions_expected = pd.DataFrame( + data={ + v2.C.CONDITION_ID: ["cond_0", "cond_1"], + v2.C.TARGET_ID: ["k1", "k1"], + v2.C.TARGET_VALUE: [2.0, 3.0], + } + ) + + assert converted.experiment_df.equals(experiments_expected) + assert converted.condition_df.equals(conditions_expected) + assert converted.measurement_df[v2.C.EXPERIMENT_ID].tolist() == [ + "experiment____c1", + "experiment____c2", + ] + + +def _simulate(problem, experiment_id, timepoints): + """Simulate an experiment in BasiCO, after converting the experiments + to SBML events, and return the trajectory of `A`.""" + 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) + basico.load_model(converted.model.to_sbml_str()) + result = basico.run_time_course( + values=timepoints, r_tol=1e-10, a_tol=1e-12 + ) + return result["A"].to_numpy() + + +@pytest.mark.parametrize("by_condition", [False, True]) +def test_assignments_to_experiments_simulate_before_switch(by_condition): + """Test that a converted experiment starts at 0, not at its first switch, + if nothing changes at 0: without an experiment, or if the condition at + 0 only sets symbols of the converted assignment. Simulations before and + after the switch match the original problem's.""" + problem = Problem() + ant_model = f""" + model conversion + species A, B; + A = 10; + B = 0; + t_switch = 0; + k1 := piecewise(1, time < {"t_switch" if by_condition else 5}, 2); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + if by_condition: + experiment_id = "experiment____c1" + problem.add_condition("c1", t_switch=3.0) + problem.add_experiment(experiment_id, 0, "c1") + else: + experiment_id = DEFAULT_EXPERIMENT_ID + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + timepoints = [0.0, 1.0, 2.0, 4.0, 6.0, 10.0] + for t in timepoints: + problem.add_measurement( + "obs_a", + experiment_id=experiment_id if by_condition else None, + time=t, + measurement=5.0, + ) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + np.testing.assert_allclose( + _simulate(converted, experiment_id, timepoints), + _simulate(problem, experiment_id, timepoints), + rtol=1e-6, + ) + (experiment,) = converted.experiments + assert experiment.id == experiment_id + assert [(p.time, p.condition_ids) for p in experiment.periods] == [ + (0.0, []), + (3.0 if by_condition else 5.0, ["cond_0"]), + ] + + +def test_assignments_to_experiments_preequilibration(): + """Test that the converted assignments are applied during + preequilibration, using the preequilibration condition's values.""" + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + dose = 0; + dose_time = 0; + k1 := dose * piecewise(0, time - dose_time < 0, 1); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + problem.add_condition("c1", dose=2.0, dose_time=10.0) + # the dose is present from the start, i.e. also during preequilibration + problem.add_condition("c2", dose=3.0, dose_time=0.0) + problem.add_experiment( + "experiment__c2___c1", v2.C.TIME_PREEQUILIBRATION, "c2", 0, "c1" + ) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", + experiment_id="experiment__c2___c1", + time=30.0, + measurement=10.0, + ) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + with pytest.warns(UserWarning, match="during preequilibration"): + converted = assignments_to_experiments(problem) + + # the dose is applied during preequilibration, so the reset to the + # pre-dose value at t=0 is not redundant here and must be kept + experiments_expected = pd.DataFrame( + data={ + v2.C.EXPERIMENT_ID: ["experiment__c2___c1"] * 3, + v2.C.TIME: [v2.C.TIME_PREEQUILIBRATION, 0.0, 10.0], + v2.C.CONDITION_ID: ["cond_0", "cond_1", "cond_2"], + } + ) + + conditions_expected = pd.DataFrame( + data={ + v2.C.CONDITION_ID: ["cond_0", "cond_1", "cond_2"], + v2.C.TARGET_ID: ["k1", "k1", "k1"], + v2.C.TARGET_VALUE: [3.0, 0.0, 2.0], + } + ) + + assert converted.experiment_df.equals(experiments_expected) + assert converted.condition_df.equals(conditions_expected) + + +@pytest.mark.parametrize( + "assignment, match", + [ + ("piecewise(0, time == 10, 1)", "single"), + ("piecewise(0, time < 10 && time > 5, 1)", "single"), + ("piecewise(0, time^2 < 10, 1)", "single"), + ("piecewise(time^2, time < 10, 1)", "linearly"), + ("piecewise(exp(time), time < 10, 1)", "linearly"), + ( + "piecewise(0, time < t_dose, 1)", + r"estimated parameters \['t_dose'\]", + ), + ("piecewise(0, time < A, 1)", "constant model parameters"), + ("A * piecewise(0, time < 10, 1)", "constant model parameters"), + ], +) +def test_assignments_to_experiments_unsupported(assignment, match): + """Test that assignments outside the supported scope are rejected.""" + problem = Problem() + ant_model = f""" + model conversion + species A, B; + A = 10; + B = 0; + t_dose = 5; + k1 := {assignment}; + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=0.0, measurement=10.0) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + problem.add_parameter("t_dose", estimate=True, lb=0, ub=10) + + with pytest.raises(NotImplementedError, match=match): + assignments_to_experiments(problem) + + +@pytest.mark.parametrize( + "condition_value, parameter, t_switch", + [ + # a constant model parameter + (None, None, 5.0), + # a parameter fixed by the parameter table + (None, {"estimate": False, "nominal_value": 7.0}, 7.0), + # set by the condition table + (8.0, None, 8.0), + ], +) +def test_assignments_to_experiments_fixed_switching_time( + condition_value, parameter, t_switch +): + """Test that switching times may depend on parameters that are fixed or + set by the condition table.""" + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + t_dose = 5; + k1 := piecewise(0, time < t_dose, 1); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + experiment_id = None + if condition_value is not None: + experiment_id = "e1" + problem.add_condition("c1", t_dose=condition_value) + problem.add_experiment(experiment_id, 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id=experiment_id, time=20.0, measurement=5.0 + ) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + if parameter is not None: + problem.add_parameter("t_dose", **parameter) + + converted = assignments_to_experiments(problem) + + periods = [ + (row[v2.C.TIME], row[v2.C.CONDITION_ID]) + for _, row in converted.experiment_df.iterrows() + ] + # the experiment starts at 0, where nothing changes + assert periods[0] == (0.0, "") + conditions = converted.condition_df.set_index(v2.C.CONDITION_ID) + assert [ + (t, float(conditions.loc[condition_id, v2.C.TARGET_VALUE])) + for t, condition_id in periods[1:] + ] == [(t_switch, 1.0)] + + +def test_assignments_to_experiments_estimated_switching_time_condition(): + """Test that a switching time that the condition table sets to an + estimated parameter is rejected.""" + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + t_dose = 5; + t_est = 5; + k1 := piecewise(0, time < t_dose, 1); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + problem.add_condition("c1", t_dose="t_est") + problem.add_experiment("e1", 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=5.0 + ) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + problem.add_parameter("t_est", estimate=True, lb=0, ub=10) + + with pytest.raises(NotImplementedError, match="evaluate to a real number"): + assignments_to_experiments(problem) + + +@pytest.mark.parametrize( + "assignment, default, switches", + [ + # a sum of doses given at different times + ( + "piecewise(0, time < 10, 1) + 2 * piecewise(0, time < 20, 1)", + 0.0, + [(10.0, 1.0), (20.0, 3.0)], + ), + # phases one after another, as several pieces + ( + "piecewise(1, time < 10, 2, time < 20, 3)", + 1.0, + [(10.0, 2.0), (20.0, 3.0)], + ), + # and as nested piecewise functions + ( + "piecewise(1, time < 10, piecewise(2, time < 20, 3))", + 1.0, + [(10.0, 2.0), (20.0, 3.0)], + ), + # a switch that does not change the value is no period + ( + "piecewise(1, time < 10, piecewise(1, time < 20, 3))", + 1.0, + [(20.0, 3.0)], + ), + ], +) +def test_assignments_to_experiments_several_switches( + assignment, default, switches +): + """Test conversion of assignments that switch more than once.""" + problem = Problem() + ant_model = f""" + model conversion + species A, B; + A = 10; + B = 0; + k1 := {assignment}; + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=30.0, measurement=5.0) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + assert converted.model.sbml_model.getParameter("k1").getValue() == default + periods = [ + (row[v2.C.TIME], row[v2.C.CONDITION_ID]) + for _, row in converted.experiment_df.iterrows() + ] + # the experiment starts at 0, where nothing changes + assert periods[0] == (0.0, "") + conditions = converted.condition_df.set_index(v2.C.CONDITION_ID) + assert [ + (t, float(conditions.loc[condition_id, v2.C.TARGET_VALUE])) + for t, condition_id in periods[1:] + ] == switches + + +def test_assignments_to_experiments_condition_specific_doses(): + """Test conversion of two doses whose times and amounts are set by the + condition table, as in Brannmark_JBC2010.""" + problem = Problem() + doses = ( + "d1 * piecewise(0, time - t1 < 0, 1)" + " + d2 * piecewise(0, time - t2 < 0, 1)" + ) + ant_model = f""" + model conversion + species A, B; + A = 10; + B = 0; + d1 = 0; t1 = 0; d2 = 0; t2 = 0; + k1 := {doses}; + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + + problem.add_condition("c1", d1=2.0, t1=0.0, d2=0.0, t2=0.0) + problem.add_condition("c2", d1=1.0, t1=0.0, d2=3.0, t2=5.0) + problem.add_experiment("experiment____c1", 0, "c1") + problem.add_experiment("experiment____c2", 0, "c2") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + for experiment_id in ("experiment____c1", "experiment____c2"): + problem.add_measurement( + "obs_a", experiment_id=experiment_id, time=30.0, measurement=5.0 + ) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + conditions = converted.condition_df.set_index(v2.C.CONDITION_ID) + got = { + (row[v2.C.EXPERIMENT_ID], row[v2.C.TIME]): float( + conditions.loc[row[v2.C.CONDITION_ID], v2.C.TARGET_VALUE] + ) + for _, row in converted.experiment_df.iterrows() + } + # c1: the dose is there from the start; c2: the first dose from the + # start, the second added at t = 5 + assert got == { + ("experiment____c1", 0.0): 2.0, + ("experiment____c2", 0.0): 1.0, + ("experiment____c2", 5.0): 4.0, + } + + +def _rate_rule(problem, target_id): + rule = problem.model.sbml_model.getRateRuleByVariable(target_id) + return None if rule is None else libsbml.formulaToL3String(rule.getMath()) + + +def _changes_by_period(problem): + """{(experiment ID, time): {target ID: target value}}""" + conditions = problem.condition_df.set_index(v2.C.CONDITION_ID) + out = {} + for _, row in problem.experiment_df.iterrows(): + change = conditions.loc[[row[v2.C.CONDITION_ID]]] + out.setdefault((row[v2.C.EXPERIMENT_ID], row[v2.C.TIME]), {}).update( + { + r[v2.C.TARGET_ID]: float(r[v2.C.TARGET_VALUE]) + for _, r in change.iterrows() + } + ) + return out + + +def test_assignments_to_experiments_ramp(): + """Test that a value changing linearly with time becomes a rate rule + whose slope each period sets.""" + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + r = 0.5; a = 1; + k1 := piecewise(r * time + a, time <= 10, 0); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=5.0) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + sbml = converted.model.sbml_model + # the model only makes the target differential; the slope is a + # placeholder until the tables set it + assert _rate_rule(converted, "k1") == "k1_slope" + assert sbml.getParameter("k1_slope").getValue() == 0.0 + # the ramp starts at time 0, and at the switch, the target and its slope + # drop to zero + assert _changes_by_period(converted) == { + (DEFAULT_EXPERIMENT_ID, 0.0): {"k1": 1.0, "k1_slope": 0.5}, + (DEFAULT_EXPERIMENT_ID, 10.0): {"k1": 0.0, "k1_slope": 0.0}, + } + + +def test_assignments_to_experiments_condition_specific_ramps(): + """Test pulses and ramps whose start, rate and end are set by the + condition table, as in Fujita_SciSignal2010.""" + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + EGF_0 = 0; EGF_rate = 0; EGF_end = 0; + k1 := piecewise(EGF_rate * time + EGF_0, time <= EGF_end, 0); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + # a pulse: a constant dose withdrawn at t = 60 + problem.add_condition("pulse", EGF_0=1.0, EGF_rate=0.0, EGF_end=60.0) + # a ramp: growing from zero at 0.01 per unit time until t = 3600 + problem.add_condition("ramp", EGF_0=0.0, EGF_rate=0.01, EGF_end=3600.0) + for c in ("pulse", "ramp"): + problem.add_experiment(f"experiment____{c}", 0, c) + problem.add_measurement( + "obs_a", + experiment_id=f"experiment____{c}", + time=30.0, + measurement=5.0, + ) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + assert _rate_rule(converted, "k1") == "k1_slope" + assert _changes_by_period(converted) == { + ("experiment____pulse", 0.0): {"k1": 1.0, "k1_slope": 0.0}, + ("experiment____pulse", 60.0): {"k1": 0.0, "k1_slope": 0.0}, + ("experiment____ramp", 0.0): {"k1": 0.0, "k1_slope": 0.01}, + ("experiment____ramp", 3600.0): {"k1": 0.0, "k1_slope": 0.0}, + } + + +def test_assignments_to_experiments_ramp_preequilibration(): + """Test that a ramp is held at its value at time 0 during + preequilibration and starts at time 0.""" + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + r = 0; + k1 := piecewise(r * time + 1, time <= 10, 0); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + problem.add_condition("pre", r=0.0) + problem.add_condition("main", r=2.0) + problem.add_experiment( + "experiment__pre___main", v2.C.TIME_PREEQUILIBRATION, "pre", 0, "main" + ) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", + experiment_id="experiment__pre___main", + time=20.0, + measurement=5.0, + ) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + with pytest.warns(UserWarning, match="during preequilibration"): + converted = assignments_to_experiments(problem) + + assert _changes_by_period(converted) == { + ("experiment__pre___main", v2.C.TIME_PREEQUILIBRATION): { + "k1": 1.0, + "k1_slope": 0.0, + }, + ("experiment__pre___main", 0.0): {"k1": 1.0, "k1_slope": 2.0}, + ("experiment__pre___main", 10.0): {"k1": 0.0, "k1_slope": 0.0}, + } + + +@pytest.mark.parametrize("t_switch", [0, 10]) +@pytest.mark.parametrize( + "condition, holds_at_switch, holds_before", + [ + ("time < {t}", False, True), + ("time <= {t}", True, True), + ("time > {t}", False, False), + ("time >= {t}", True, False), + ("{t} >= time", True, True), + ], +) +def test_assignments_to_experiments_preequilibration_boundary( + t_switch, condition, holds_at_switch, holds_before +): + """Test that preequilibration uses the assignment's value at time 0, + including whether the condition holds for a switch at time 0.""" + condition = condition.format(t=t_switch) + problem = Problem() + ant_model = f""" + model conversion + species A, B; + A = 10; + B = 0; + k1 := piecewise(0, {condition}, 1); + k2 = 0.5; + R1: A -> B; k1 * A; + R2: B -> A; k2 * B; + end + """ + problem.model = SbmlModel.from_antimony(ant_model) + problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, [], 0, []) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=5.0 + ) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + with pytest.warns(UserWarning, match="during preequilibration"): + converted = assignments_to_experiments(problem) + + holds_at_0 = holds_at_switch if t_switch == 0 else holds_before + (preeq_period,) = ( + p for p in converted.experiments[0].periods if p.is_preequilibration + ) + assert [ + (c.target_id, float(c.target_value)) + for c in converted.get_changes_for_period(preeq_period) + ] == [("k1", 0.0 if holds_at_0 else 1.0)] + + +def test_assignments_to_experiments_validate(): + """Test that the converted problem is linted. + + Here, the input problem is valid, but the measurement without an + experiment ID is assigned to a new experiment `default_experiment`, an ID + that an existing experiment already has. + """ + problem = Problem() + problem.model = _switch_model("piecewise(0, time < 10, 1)") + problem.add_condition("c1", other=2.0) + problem.add_experiment(DEFAULT_EXPERIMENT_ID, 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", + experiment_id=DEFAULT_EXPERIMENT_ID, + time=20.0, + measurement=5.0, + ) + problem.add_measurement("obs_a", time=30.0, measurement=5.0) + assert not problem.validate().has_errors() + + with pytest.raises(ValueError, match="did not pass linting"): + assignments_to_experiments(problem) + + +def _changes_by_time( + problem: Problem, experiment_id: str +) -> dict[float, dict[str, float]]: + """The changes of each period of an experiment, by the period's time.""" + return { + period.time: { + c.target_id: float(c.target_value) + for c in problem.get_changes_for_period(period) + } + for period in problem[experiment_id].sorted_periods + } + + +def _switch_model(assignment: str, extra: str = "") -> SbmlModel: + """A model with the assignment rule `k1 := `.""" + return SbmlModel.from_antimony(f""" + model conversion + species A; + A = 10; + other = 1; + {extra} + k1 := {assignment}; + R1: A -> ; k1 * A; + end + """) + + +def test_assignments_to_experiments_switch_at_period_start(): + """Test that a switch at the start of a later period is kept, although + that period does not change any symbol of the assignment.""" + problem = Problem() + problem.model = _switch_model("piecewise(0, time < 10, 1)") + problem.add_condition("c0", other=1.0) + problem.add_condition("c1", other=2.0) + problem.add_experiment("e1", 0, "c0", 10, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + assert _changes_by_time(converted, "e1") == { + 0.0: {"other": 1.0}, + 10.0: {"other": 2.0, "k1": 1.0}, + } + + +def test_assignments_to_experiments_switch_at_preequilibration_end(): + """Test that a switch at time 0 is kept after preequilibration, where the + condition still holds at time 0.""" + problem = Problem() + problem.model = _switch_model("piecewise(0, time <= 0, 1)") + problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, [], 0, []) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + with pytest.warns(UserWarning, match="during preequilibration"): + converted = assignments_to_experiments(problem) + + assert _changes_by_time(converted, "e1") == { + v2.C.TIME_PREEQUILIBRATION: {"k1": 0.0}, + 0.0: {"k1": 1.0}, + } + + +def test_assignments_to_experiments_shared_symbol(): + """Test that the changes of a symbol of a converted assignment stay if + the model still refers to that symbol elsewhere.""" + problem = Problem() + problem.model = _switch_model( + "dose * piecewise(0, time < 10, 1)", + extra="dose = 0; R2: -> A; dose;", + ) + problem.add_condition("c1", dose=2.0) + problem.add_experiment("e1", 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + assert _changes_by_time(converted, "e1") == { + 0.0: {"dose": 2.0}, + 10.0: {"k1": 2.0}, + } + + +def test_assignments_to_experiments_initial_assignment(): + """Test that a time-dependent initial assignment stays in the model, since + it is evaluated once, at the start of the simulation.""" + problem = Problem() + problem.model = SbmlModel.from_antimony(""" + model conversion + species A; + A = 10; + k1 = piecewise(0, time < 10, 1); + R1: A -> ; k1 * A; + end + """) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=1.0) + + converted = assignments_to_experiments(problem) + + assert converted.model.sbml_model.getInitialAssignment("k1") is not None + assert not converted.experiments + assert not converted.conditions + + +def test_assignments_to_experiments_preequilibration_model_value(): + """Test that a value at preequilibration, which starts the experiment, + only refers to parameter table symbols: other model parameters take their + values.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(level, time < 10, 0)", extra="level = 3;" + ) + problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, [], 0, []) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + with pytest.warns(UserWarning, match="during preequilibration"): + converted = assignments_to_experiments(problem) + + assert _changes_by_time(converted, "e1") == { + v2.C.TIME_PREEQUILIBRATION: {"k1": 3.0}, + 0.0: {}, + 10.0: {"k1": 0.0}, + } + + +def test_assignments_to_experiments_species_target(): + """Test that an assignment to a species is rejected.""" + problem = Problem() + problem.model = SbmlModel.from_antimony(""" + model conversion + species A, S; + A = 10; + S := piecewise(0, time < 10, 1); + R1: A -> ; S * A; + end + """) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=1.0) + + with pytest.raises(NotImplementedError, match="not a parameter"): + assignments_to_experiments(problem) + + +def test_assignments_to_experiments_nothing_to_convert(): + """Test that a problem without time-dependent assignments is unchanged, + including an experiment that only preequilibrates and periods that are + not in time order.""" + problem = Problem() + problem.model = SbmlModel.from_antimony(""" + model conversion + species A; + A = 10; + k2 = 0.5; + R1: A -> ; k2 * A; + end + """) + problem.add_condition("c1", k2=1.0) + problem.add_condition("c2", k2=2.0) + problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, "c1") + # periods out of time order + problem.add_experiment("e2", 10, "c2", 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + for experiment_id in ("e1", "e2"): + problem.add_measurement( + "obs_a", experiment_id=experiment_id, time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + def periods(problem): + return [ + [(p.time, p.condition_ids) for p in e.periods] + for e in problem.experiments + ] + + assert periods(converted) == periods(problem) + assert converted.condition_df.equals(problem.condition_df) + + +@pytest.mark.parametrize("on", ["model", "species"]) +def test_assignments_to_experiments_conversion_factor(on): + """Test that the changes of a symbol of a converted assignment stay if + the model refers to that symbol as a conversion factor.""" + problem = Problem() + problem.model = _switch_model( + "cf * piecewise(0, time < 10, 1)", extra="cf = 1;" + ) + sbml_model = problem.model.sbml_model + element = sbml_model if on == "model" else sbml_model.getSpecies("A") + assert ( + element.setConversionFactor("cf") == libsbml.LIBSBML_OPERATION_SUCCESS + ) + problem.add_condition("c1", cf=2.0) + problem.add_experiment("e1", 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + assert _changes_by_time(converted, "e1") == { + 0.0: {"cf": 2.0}, + 10.0: {"k1": 2.0}, + } + + +def test_assignments_to_experiments_slope_id_taken(): + """Test that a ramp's slope parameter gets an ID that the problem does not + use yet, here because an observable parameter is called `k1_slope`.""" + problem = Problem() + problem.model = _switch_model("piecewise(time, time < 10, 1)") + problem.add_observable("obs_a", formula="k1_slope * A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=5.0) + problem.add_parameter("k1_slope", estimate=True, lb=1e-3, ub=1e3) + + converted = assignments_to_experiments(problem) + + assert _rate_rule(converted, "k1") == "k1_slope_1" + assert _changes_by_time(converted, DEFAULT_EXPERIMENT_ID) == { + 0.0: {"k1": 0.0, "k1_slope_1": 1.0}, + 10.0: {"k1": 1.0, "k1_slope_1": 0.0}, + } + + +@pytest.mark.parametrize( + "conditions", + [ + # the changes of a condition are simultaneous: `d1` takes the old `b1` + [(2.0, {"b1": 2.0, "d1": "b1 + 0.5"})], + # a change is evaluated once: `d1` keeps the `b1` of its period + [(2.0, {"d1": "b1 + 0.5"}), (3.0, {"b1": 2.0})], + ], +) +def test_assignments_to_experiments_change_evaluated_at_period_start( + conditions, +): + """Test that a condition change's value is evaluated at the start of its + period, with the values from before the period, also when a later switch + uses it.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(b0, time < 5, d1)", extra="b0 = 0.3; b1 = 1.2; d1 = 1;" + ) + args = [0, []] + for i, (t, changes) in enumerate(conditions): + problem.add_condition(f"c{i}", **changes) + args += [t, f"c{i}"] + problem.add_experiment("e1", *args) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + converted = assignments_to_experiments(problem) + + (switch,) = (p for p in converted["e1"].periods if p.time == 5.0) + (change,) = converted.get_changes_for_period(switch) + assert change.target_id == "k1" + # 1.2 + 0.5, not 2.0 + 0.5 + assert float(change.target_value) == pytest.approx(1.7) + + +def test_assignments_to_experiments_change_refers_to_species(): + """Test that a change of a symbol of the assignment whose value refers to + a species, which is not known at the start of the period, is + rejected.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(0, time < 5, d1)", extra="d1 = 1;" + ) + problem.add_condition("c1", d1="A") + problem.add_experiment("e1", 0, [], 2, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + + with pytest.raises(NotImplementedError, match=r"refer to \['A'\]"): + assignments_to_experiments(problem) + + +def test_assignments_to_experiments_petab_symbols(): + """Test that the values of new condition changes use PEtab's symbols, so + that they equal the same expressions parsed by PEtab.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(0, time < 5, d1)", extra="d1 = 1;" + ) + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", time=20.0, measurement=1.0) + problem.add_parameter( + "d1", estimate=True, nominal_value=2.0, lb=0.1, ub=10 + ) + + converted = assignments_to_experiments(problem) + + assert [ + change.target_value + for condition in converted.conditions + for change in condition.changes + ] == [sympify_petab("d1")] + + +@pytest.mark.parametrize("estimate", [False, True]) +def test_assignments_to_experiments_switching_time_from_fixed_parameter( + estimate, +): + """Test that a condition may set a switching time to an expression in a + parameter that the parameter table fixes, which takes its nominal value, + but not in an estimated one.""" + problem = Problem() + problem.model = _switch_model( + "piecewise(0, time < t1, 1)", extra="t1 = 0; p = 1;" + ) + problem.add_condition("c1", t1="p + 1") + problem.add_experiment("e1", 0, "c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", experiment_id="e1", time=20.0, measurement=1.0 + ) + if estimate: + problem.add_parameter( + "p", estimate=True, nominal_value=4.0, lb=0.1, ub=10 + ) + with pytest.raises( + NotImplementedError, match="evaluate to a real number" + ): + assignments_to_experiments(problem) + return + problem.add_parameter("p", estimate=False, nominal_value=4.0) + + converted = assignments_to_experiments(problem) + + # the switch is at the nominal value 4 + 1, not the model value 1 + 1 + assert _changes_by_time(converted, "e1") == {0.0: {}, 5.0: {"k1": 1.0}} diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 72deade4..f06ece64 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -1,7 +1,6 @@ import logging from pathlib import Path -import libsbml import numpy as np import pandas as pd import pytest @@ -12,12 +11,6 @@ from petab.v1.priors import Prior from petab.v2 import Problem from petab.v2.calculate import calculate_llh as calculate_llh_v2 -from petab.v2.converters import ( - DEFAULT_EXPERIMENT_ID, - ExperimentsToSbmlConverter, - assignments_to_experiments, -) -from petab.v2.math import sympify_petab from petab.v2.petab1to2 import petab1to2, v1v2_observable_df @@ -657,1194 +650,3 @@ def test_benchmark_collection(problem_id): pytest.skip(str(e)) assert isinstance(problem, Problem) assert len(problem.measurements) - - -@pytest.mark.parametrize( - "condition, before, after", - [("time < 10", 1.0, 2.0), ("10 < time", 2.0, 1.0)], -) -def test_assignments_to_experiments(condition, before, after): - """Test that switches in either direction are converted.""" - problem = Problem() - ant_model = f""" - model conversion - species A, B; - A = 10; - B = 0; - k1 := piecewise(1, {condition}, 2); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement("obs_a", time=0.0, measurement=10.0) - problem.add_measurement("obs_a", time=20.0, measurement=5.0) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - # the value before the switch is kept in the model, so that only the - # switch itself needs a condition; the experiment still starts at 0, - # where the first measurement is - assert converted.model.sbml_model.getParameter("k1").getValue() == before - - experiments_expected = pd.DataFrame( - data={ - v2.C.EXPERIMENT_ID: [DEFAULT_EXPERIMENT_ID] * 2, - v2.C.TIME: [0.0, 10.0], - v2.C.CONDITION_ID: ["", "cond_0"], - } - ) - - conditions_expected = pd.DataFrame( - data={ - v2.C.CONDITION_ID: ["cond_0"], - v2.C.TARGET_ID: ["k1"], - v2.C.TARGET_VALUE: [after], - } - ) - - assert converted.experiment_df.equals(experiments_expected) - assert converted.condition_df.equals(conditions_expected) - assert converted.measurement_df[v2.C.EXPERIMENT_ID].tolist() == [ - DEFAULT_EXPERIMENT_ID, - DEFAULT_EXPERIMENT_ID, - ] - - -def test_assignments_to_experiments_unused_condition_id(): - """Test that new condition IDs don't clash with other IDs.""" - problem = Problem() - problem.model = SbmlModel.from_antimony(""" - species A = 10; - k1 := piecewise(1, time < 10, 2); - R1: A -> ; k1 * A; - """) - problem.add_observable("obs_a", formula="A", noise_formula="cond_0") - problem.add_measurement("obs_a", time=20.0, measurement=5.0) - problem.add_parameter("cond_0", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - assert [c.id for c in converted.conditions] == ["cond_0_1"] - # the experiment starts at 0, where nothing changes - assert [ - (p.time, p.condition_ids) for p in converted.experiments[0].periods - ] == [(0.0, []), (10.0, ["cond_0_1"])] - - -def test_assignments_to_experiments_nested(): - """Test conversion of piecewise assignments that are nested inside a - larger expression and switch at a time given by some expression.""" - problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - level = 3; - k1 := level * piecewise(0, time - 10 < 0, 1); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement("obs_a", time=0.0, measurement=10.0) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - assert converted.model.sbml_model.getParameter("k1").getValue() == 0.0 - - # `level` is not a parameter table entry, so it must not appear in the - # first period -- the switch may refer to it, since the experiment - # starts with an empty period at 0 - experiments_expected = pd.DataFrame( - data={ - v2.C.EXPERIMENT_ID: [DEFAULT_EXPERIMENT_ID] * 2, - v2.C.TIME: [0.0, 10.0], - v2.C.CONDITION_ID: ["", "cond_0"], - } - ) - - # the factor outside the piecewise expression is preserved - conditions_expected = pd.DataFrame( - data={ - v2.C.CONDITION_ID: ["cond_0"], - v2.C.TARGET_ID: ["k1"], - v2.C.TARGET_VALUE: ["level"], - } - ) - - assert converted.experiment_df.equals(experiments_expected) - assert converted.condition_df.equals(conditions_expected) - - -def test_assignments_to_experiments_parameter_table_default(): - """Test that a pre-switch value that refers to parameter table symbols is - set at the start, since the model only has a number for it.""" - problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - d1 = 1; - d2 = 2; - k1 := piecewise(d1, time < 10, d2); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement("obs_a", time=20.0, measurement=5.0) - problem.add_parameter( - "d1", estimate=True, lb=1e-3, ub=1e3, nominal_value=3 - ) - problem.add_parameter( - "d2", estimate=True, lb=1e-3, ub=1e3, nominal_value=4 - ) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - experiments_expected = pd.DataFrame( - data={ - v2.C.EXPERIMENT_ID: [DEFAULT_EXPERIMENT_ID] * 2, - v2.C.TIME: [0.0, 10.0], - v2.C.CONDITION_ID: ["cond_0", "cond_1"], - } - ) - - conditions_expected = pd.DataFrame( - data={ - v2.C.CONDITION_ID: ["cond_0", "cond_1"], - v2.C.TARGET_ID: ["k1", "k1"], - v2.C.TARGET_VALUE: ["d1", "d2"], - } - ) - - assert converted.experiment_df.equals(experiments_expected) - assert converted.condition_df.equals(conditions_expected) - - -def test_assignments_to_experiments_condition_specific(): - """Test conversion of piecewise assignments whose switching times and - values are set by the condition table.""" - problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - dose = 0; - dose_time = 0; - k1 := dose * piecewise(0, time - dose_time < 0, 1); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - problem.add_condition("c1", dose=2.0, dose_time=10.0) - problem.add_condition("c2", dose=3.0, dose_time=20.0) - problem.add_experiment("experiment____c1", 0, "c1") - problem.add_experiment("experiment____c2", 0, "c2") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="experiment____c1", time=30.0, measurement=10.0 - ) - problem.add_measurement( - "obs_a", experiment_id="experiment____c2", time=30.0, measurement=5.0 - ) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - # the pre-dose value is the same for both conditions and is kept in the - # model; each condition gets its own experiment and switching time, and - # each experiment still starts at 0, where its condition is now empty - assert converted.model.sbml_model.getParameter("k1").getValue() == 0.0 - - experiments_expected = pd.DataFrame( - data={ - v2.C.EXPERIMENT_ID: ["experiment____c1"] * 2 - + ["experiment____c2"] * 2, - v2.C.TIME: [0.0, 10.0, 0.0, 20.0], - v2.C.CONDITION_ID: ["", "cond_0", "", "cond_1"], - } - ) - - # `dose` and `dose_time` are covered by the new conditions - conditions_expected = pd.DataFrame( - data={ - v2.C.CONDITION_ID: ["cond_0", "cond_1"], - v2.C.TARGET_ID: ["k1", "k1"], - v2.C.TARGET_VALUE: [2.0, 3.0], - } - ) - - assert converted.experiment_df.equals(experiments_expected) - assert converted.condition_df.equals(conditions_expected) - assert converted.measurement_df[v2.C.EXPERIMENT_ID].tolist() == [ - "experiment____c1", - "experiment____c2", - ] - - -def _simulate(problem, experiment_id, timepoints): - """Simulate an experiment in BasiCO, after converting the experiments - to SBML events, and return the trajectory of `A`.""" - 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) - basico.load_model(converted.model.to_sbml_str()) - result = basico.run_time_course( - values=timepoints, r_tol=1e-10, a_tol=1e-12 - ) - return result["A"].to_numpy() - - -@pytest.mark.parametrize("by_condition", [False, True]) -def test_assignments_to_experiments_simulate_before_switch(by_condition): - """Test that a converted experiment starts at 0, not at its first switch, - if nothing changes at 0: without an experiment, or if the condition at - 0 only sets symbols of the converted assignment. Simulations before and - after the switch match the original problem's.""" - problem = Problem() - ant_model = f""" - model conversion - species A, B; - A = 10; - B = 0; - t_switch = 0; - k1 := piecewise(1, time < {"t_switch" if by_condition else 5}, 2); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - if by_condition: - experiment_id = "experiment____c1" - problem.add_condition("c1", t_switch=3.0) - problem.add_experiment(experiment_id, 0, "c1") - else: - experiment_id = DEFAULT_EXPERIMENT_ID - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - timepoints = [0.0, 1.0, 2.0, 4.0, 6.0, 10.0] - for t in timepoints: - problem.add_measurement( - "obs_a", - experiment_id=experiment_id if by_condition else None, - time=t, - measurement=5.0, - ) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - np.testing.assert_allclose( - _simulate(converted, experiment_id, timepoints), - _simulate(problem, experiment_id, timepoints), - rtol=1e-6, - ) - (experiment,) = converted.experiments - assert experiment.id == experiment_id - assert [(p.time, p.condition_ids) for p in experiment.periods] == [ - (0.0, []), - (3.0 if by_condition else 5.0, ["cond_0"]), - ] - - -def test_assignments_to_experiments_preequilibration(): - """Test that the converted assignments are applied during - preequilibration, using the preequilibration condition's values.""" - problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - dose = 0; - dose_time = 0; - k1 := dose * piecewise(0, time - dose_time < 0, 1); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - problem.add_condition("c1", dose=2.0, dose_time=10.0) - # the dose is present from the start, i.e. also during preequilibration - problem.add_condition("c2", dose=3.0, dose_time=0.0) - problem.add_experiment( - "experiment__c2___c1", v2.C.TIME_PREEQUILIBRATION, "c2", 0, "c1" - ) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", - experiment_id="experiment__c2___c1", - time=30.0, - measurement=10.0, - ) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - with pytest.warns(UserWarning, match="during preequilibration"): - converted = assignments_to_experiments(problem) - - # the dose is applied during preequilibration, so the reset to the - # pre-dose value at t=0 is not redundant here and must be kept - experiments_expected = pd.DataFrame( - data={ - v2.C.EXPERIMENT_ID: ["experiment__c2___c1"] * 3, - v2.C.TIME: [v2.C.TIME_PREEQUILIBRATION, 0.0, 10.0], - v2.C.CONDITION_ID: ["cond_0", "cond_1", "cond_2"], - } - ) - - conditions_expected = pd.DataFrame( - data={ - v2.C.CONDITION_ID: ["cond_0", "cond_1", "cond_2"], - v2.C.TARGET_ID: ["k1", "k1", "k1"], - v2.C.TARGET_VALUE: [3.0, 0.0, 2.0], - } - ) - - assert converted.experiment_df.equals(experiments_expected) - assert converted.condition_df.equals(conditions_expected) - - -@pytest.mark.parametrize( - "assignment, match", - [ - ("piecewise(0, time == 10, 1)", "single"), - ("piecewise(0, time < 10 && time > 5, 1)", "single"), - ("piecewise(0, time^2 < 10, 1)", "single"), - ("piecewise(time^2, time < 10, 1)", "linearly"), - ("piecewise(exp(time), time < 10, 1)", "linearly"), - ( - "piecewise(0, time < t_dose, 1)", - r"estimated parameters \['t_dose'\]", - ), - ("piecewise(0, time < A, 1)", "constant model parameters"), - ("A * piecewise(0, time < 10, 1)", "constant model parameters"), - ], -) -def test_assignments_to_experiments_unsupported(assignment, match): - """Test that assignments outside the supported scope are rejected.""" - problem = Problem() - ant_model = f""" - model conversion - species A, B; - A = 10; - B = 0; - t_dose = 5; - k1 := {assignment}; - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement("obs_a", time=0.0, measurement=10.0) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - problem.add_parameter("t_dose", estimate=True, lb=0, ub=10) - - with pytest.raises(NotImplementedError, match=match): - assignments_to_experiments(problem) - - -@pytest.mark.parametrize( - "condition_value, parameter, t_switch", - [ - # a constant model parameter - (None, None, 5.0), - # a parameter fixed by the parameter table - (None, {"estimate": False, "nominal_value": 7.0}, 7.0), - # set by the condition table - (8.0, None, 8.0), - ], -) -def test_assignments_to_experiments_fixed_switching_time( - condition_value, parameter, t_switch -): - """Test that switching times may depend on parameters that are fixed or - set by the condition table.""" - problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - t_dose = 5; - k1 := piecewise(0, time < t_dose, 1); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - experiment_id = None - if condition_value is not None: - experiment_id = "e1" - problem.add_condition("c1", t_dose=condition_value) - problem.add_experiment(experiment_id, 0, "c1") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id=experiment_id, time=20.0, measurement=5.0 - ) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - if parameter is not None: - problem.add_parameter("t_dose", **parameter) - - converted = assignments_to_experiments(problem) - - periods = [ - (row[v2.C.TIME], row[v2.C.CONDITION_ID]) - for _, row in converted.experiment_df.iterrows() - ] - # the experiment starts at 0, where nothing changes - assert periods[0] == (0.0, "") - conditions = converted.condition_df.set_index(v2.C.CONDITION_ID) - assert [ - (t, float(conditions.loc[condition_id, v2.C.TARGET_VALUE])) - for t, condition_id in periods[1:] - ] == [(t_switch, 1.0)] - - -def test_assignments_to_experiments_estimated_switching_time_condition(): - """Test that a switching time that the condition table sets to an - estimated parameter is rejected.""" - problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - t_dose = 5; - t_est = 5; - k1 := piecewise(0, time < t_dose, 1); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - problem.add_condition("c1", t_dose="t_est") - problem.add_experiment("e1", 0, "c1") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=5.0 - ) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - problem.add_parameter("t_est", estimate=True, lb=0, ub=10) - - with pytest.raises(NotImplementedError, match="evaluate to a real number"): - assignments_to_experiments(problem) - - -@pytest.mark.parametrize( - "assignment, default, switches", - [ - # a sum of doses given at different times - ( - "piecewise(0, time < 10, 1) + 2 * piecewise(0, time < 20, 1)", - 0.0, - [(10.0, 1.0), (20.0, 3.0)], - ), - # phases one after another, as several pieces - ( - "piecewise(1, time < 10, 2, time < 20, 3)", - 1.0, - [(10.0, 2.0), (20.0, 3.0)], - ), - # and as nested piecewise functions - ( - "piecewise(1, time < 10, piecewise(2, time < 20, 3))", - 1.0, - [(10.0, 2.0), (20.0, 3.0)], - ), - # a switch that does not change the value is no period - ( - "piecewise(1, time < 10, piecewise(1, time < 20, 3))", - 1.0, - [(20.0, 3.0)], - ), - ], -) -def test_assignments_to_experiments_several_switches( - assignment, default, switches -): - """Test conversion of assignments that switch more than once.""" - problem = Problem() - ant_model = f""" - model conversion - species A, B; - A = 10; - B = 0; - k1 := {assignment}; - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement("obs_a", time=30.0, measurement=5.0) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - assert converted.model.sbml_model.getParameter("k1").getValue() == default - periods = [ - (row[v2.C.TIME], row[v2.C.CONDITION_ID]) - for _, row in converted.experiment_df.iterrows() - ] - # the experiment starts at 0, where nothing changes - assert periods[0] == (0.0, "") - conditions = converted.condition_df.set_index(v2.C.CONDITION_ID) - assert [ - (t, float(conditions.loc[condition_id, v2.C.TARGET_VALUE])) - for t, condition_id in periods[1:] - ] == switches - - -def test_assignments_to_experiments_condition_specific_doses(): - """Test conversion of two doses whose times and amounts are set by the - condition table, as in Brannmark_JBC2010.""" - problem = Problem() - doses = ( - "d1 * piecewise(0, time - t1 < 0, 1)" - " + d2 * piecewise(0, time - t2 < 0, 1)" - ) - ant_model = f""" - model conversion - species A, B; - A = 10; - B = 0; - d1 = 0; t1 = 0; d2 = 0; t2 = 0; - k1 := {doses}; - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - - problem.add_condition("c1", d1=2.0, t1=0.0, d2=0.0, t2=0.0) - problem.add_condition("c2", d1=1.0, t1=0.0, d2=3.0, t2=5.0) - problem.add_experiment("experiment____c1", 0, "c1") - problem.add_experiment("experiment____c2", 0, "c2") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - for experiment_id in ("experiment____c1", "experiment____c2"): - problem.add_measurement( - "obs_a", experiment_id=experiment_id, time=30.0, measurement=5.0 - ) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - conditions = converted.condition_df.set_index(v2.C.CONDITION_ID) - got = { - (row[v2.C.EXPERIMENT_ID], row[v2.C.TIME]): float( - conditions.loc[row[v2.C.CONDITION_ID], v2.C.TARGET_VALUE] - ) - for _, row in converted.experiment_df.iterrows() - } - # c1: the dose is there from the start; c2: the first dose from the - # start, the second added at t = 5 - assert got == { - ("experiment____c1", 0.0): 2.0, - ("experiment____c2", 0.0): 1.0, - ("experiment____c2", 5.0): 4.0, - } - - -def _rate_rule(problem, target_id): - rule = problem.model.sbml_model.getRateRuleByVariable(target_id) - return None if rule is None else libsbml.formulaToL3String(rule.getMath()) - - -def _changes_by_period(problem): - """{(experiment ID, time): {target ID: target value}}""" - conditions = problem.condition_df.set_index(v2.C.CONDITION_ID) - out = {} - for _, row in problem.experiment_df.iterrows(): - change = conditions.loc[[row[v2.C.CONDITION_ID]]] - out.setdefault((row[v2.C.EXPERIMENT_ID], row[v2.C.TIME]), {}).update( - { - r[v2.C.TARGET_ID]: float(r[v2.C.TARGET_VALUE]) - for _, r in change.iterrows() - } - ) - return out - - -def test_assignments_to_experiments_ramp(): - """Test that a value changing linearly with time becomes a rate rule - whose slope each period sets.""" - problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - r = 0.5; a = 1; - k1 := piecewise(r * time + a, time <= 10, 0); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement("obs_a", time=20.0, measurement=5.0) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - sbml = converted.model.sbml_model - # the model only makes the target differential; the slope is a - # placeholder until the tables set it - assert _rate_rule(converted, "k1") == "k1_slope" - assert sbml.getParameter("k1_slope").getValue() == 0.0 - # the ramp starts at time 0, and at the switch, the target and its slope - # drop to zero - assert _changes_by_period(converted) == { - (DEFAULT_EXPERIMENT_ID, 0.0): {"k1": 1.0, "k1_slope": 0.5}, - (DEFAULT_EXPERIMENT_ID, 10.0): {"k1": 0.0, "k1_slope": 0.0}, - } - - -def test_assignments_to_experiments_condition_specific_ramps(): - """Test pulses and ramps whose start, rate and end are set by the - condition table, as in Fujita_SciSignal2010.""" - problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - EGF_0 = 0; EGF_rate = 0; EGF_end = 0; - k1 := piecewise(EGF_rate * time + EGF_0, time <= EGF_end, 0); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - # a pulse: a constant dose withdrawn at t = 60 - problem.add_condition("pulse", EGF_0=1.0, EGF_rate=0.0, EGF_end=60.0) - # a ramp: growing from zero at 0.01 per unit time until t = 3600 - problem.add_condition("ramp", EGF_0=0.0, EGF_rate=0.01, EGF_end=3600.0) - for c in ("pulse", "ramp"): - problem.add_experiment(f"experiment____{c}", 0, c) - problem.add_measurement( - "obs_a", - experiment_id=f"experiment____{c}", - time=30.0, - measurement=5.0, - ) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - assert _rate_rule(converted, "k1") == "k1_slope" - assert _changes_by_period(converted) == { - ("experiment____pulse", 0.0): {"k1": 1.0, "k1_slope": 0.0}, - ("experiment____pulse", 60.0): {"k1": 0.0, "k1_slope": 0.0}, - ("experiment____ramp", 0.0): {"k1": 0.0, "k1_slope": 0.01}, - ("experiment____ramp", 3600.0): {"k1": 0.0, "k1_slope": 0.0}, - } - - -def test_assignments_to_experiments_ramp_preequilibration(): - """Test that a ramp is held at its value at time 0 during - preequilibration and starts at time 0.""" - problem = Problem() - ant_model = """ - model conversion - species A, B; - A = 10; - B = 0; - r = 0; - k1 := piecewise(r * time + 1, time <= 10, 0); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - problem.add_condition("pre", r=0.0) - problem.add_condition("main", r=2.0) - problem.add_experiment( - "experiment__pre___main", v2.C.TIME_PREEQUILIBRATION, "pre", 0, "main" - ) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", - experiment_id="experiment__pre___main", - time=20.0, - measurement=5.0, - ) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - with pytest.warns(UserWarning, match="during preequilibration"): - converted = assignments_to_experiments(problem) - - assert _changes_by_period(converted) == { - ("experiment__pre___main", v2.C.TIME_PREEQUILIBRATION): { - "k1": 1.0, - "k1_slope": 0.0, - }, - ("experiment__pre___main", 0.0): {"k1": 1.0, "k1_slope": 2.0}, - ("experiment__pre___main", 10.0): {"k1": 0.0, "k1_slope": 0.0}, - } - - -@pytest.mark.parametrize("t_switch", [0, 10]) -@pytest.mark.parametrize( - "condition, holds_at_switch, holds_before", - [ - ("time < {t}", False, True), - ("time <= {t}", True, True), - ("time > {t}", False, False), - ("time >= {t}", True, False), - ("{t} >= time", True, True), - ], -) -def test_assignments_to_experiments_preequilibration_boundary( - t_switch, condition, holds_at_switch, holds_before -): - """Test that preequilibration uses the assignment's value at time 0, - including whether the condition holds for a switch at time 0.""" - condition = condition.format(t=t_switch) - problem = Problem() - ant_model = f""" - model conversion - species A, B; - A = 10; - B = 0; - k1 := piecewise(0, {condition}, 1); - k2 = 0.5; - R1: A -> B; k1 * A; - R2: B -> A; k2 * B; - end - """ - problem.model = SbmlModel.from_antimony(ant_model) - problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, [], 0, []) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=5.0 - ) - problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) - - with pytest.warns(UserWarning, match="during preequilibration"): - converted = assignments_to_experiments(problem) - - holds_at_0 = holds_at_switch if t_switch == 0 else holds_before - (preeq_period,) = ( - p for p in converted.experiments[0].periods if p.is_preequilibration - ) - assert [ - (c.target_id, float(c.target_value)) - for c in converted.get_changes_for_period(preeq_period) - ] == [("k1", 0.0 if holds_at_0 else 1.0)] - - -def test_assignments_to_experiments_validate(): - """Test that the converted problem is linted. - - Here, the input problem is valid, but the measurement without an - experiment ID is assigned to a new experiment `default_experiment`, an ID - that an existing experiment already has. - """ - problem = Problem() - problem.model = _switch_model("piecewise(0, time < 10, 1)") - problem.add_condition("c1", other=2.0) - problem.add_experiment(DEFAULT_EXPERIMENT_ID, 0, "c1") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", - experiment_id=DEFAULT_EXPERIMENT_ID, - time=20.0, - measurement=5.0, - ) - problem.add_measurement("obs_a", time=30.0, measurement=5.0) - assert not problem.validate().has_errors() - - with pytest.raises(ValueError, match="did not pass linting"): - assignments_to_experiments(problem) - - -def _changes_by_time( - problem: Problem, experiment_id: str -) -> dict[float, dict[str, float]]: - """The changes of each period of an experiment, by the period's time.""" - return { - period.time: { - c.target_id: float(c.target_value) - for c in problem.get_changes_for_period(period) - } - for period in problem[experiment_id].sorted_periods - } - - -def _switch_model(assignment: str, extra: str = "") -> SbmlModel: - """A model with the assignment rule `k1 := `.""" - return SbmlModel.from_antimony(f""" - model conversion - species A; - A = 10; - other = 1; - {extra} - k1 := {assignment}; - R1: A -> ; k1 * A; - end - """) - - -def test_assignments_to_experiments_switch_at_period_start(): - """Test that a switch at the start of a later period is kept, although - that period does not change any symbol of the assignment.""" - problem = Problem() - problem.model = _switch_model("piecewise(0, time < 10, 1)") - problem.add_condition("c0", other=1.0) - problem.add_condition("c1", other=2.0) - problem.add_experiment("e1", 0, "c0", 10, "c1") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=1.0 - ) - - converted = assignments_to_experiments(problem) - - assert _changes_by_time(converted, "e1") == { - 0.0: {"other": 1.0}, - 10.0: {"other": 2.0, "k1": 1.0}, - } - - -def test_assignments_to_experiments_switch_at_preequilibration_end(): - """Test that a switch at time 0 is kept after preequilibration, where the - condition still holds at time 0.""" - problem = Problem() - problem.model = _switch_model("piecewise(0, time <= 0, 1)") - problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, [], 0, []) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=1.0 - ) - - with pytest.warns(UserWarning, match="during preequilibration"): - converted = assignments_to_experiments(problem) - - assert _changes_by_time(converted, "e1") == { - v2.C.TIME_PREEQUILIBRATION: {"k1": 0.0}, - 0.0: {"k1": 1.0}, - } - - -def test_assignments_to_experiments_shared_symbol(): - """Test that the changes of a symbol of a converted assignment stay if - the model still refers to that symbol elsewhere.""" - problem = Problem() - problem.model = _switch_model( - "dose * piecewise(0, time < 10, 1)", - extra="dose = 0; R2: -> A; dose;", - ) - problem.add_condition("c1", dose=2.0) - problem.add_experiment("e1", 0, "c1") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=1.0 - ) - - converted = assignments_to_experiments(problem) - - assert _changes_by_time(converted, "e1") == { - 0.0: {"dose": 2.0}, - 10.0: {"k1": 2.0}, - } - - -def test_assignments_to_experiments_initial_assignment(): - """Test that a time-dependent initial assignment stays in the model, since - it is evaluated once, at the start of the simulation.""" - problem = Problem() - problem.model = SbmlModel.from_antimony(""" - model conversion - species A; - A = 10; - k1 = piecewise(0, time < 10, 1); - R1: A -> ; k1 * A; - end - """) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement("obs_a", time=20.0, measurement=1.0) - - converted = assignments_to_experiments(problem) - - assert converted.model.sbml_model.getInitialAssignment("k1") is not None - assert not converted.experiments - assert not converted.conditions - - -def test_assignments_to_experiments_preequilibration_model_value(): - """Test that a value at preequilibration, which starts the experiment, - only refers to parameter table symbols: other model parameters take their - values.""" - problem = Problem() - problem.model = _switch_model( - "piecewise(level, time < 10, 0)", extra="level = 3;" - ) - problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, [], 0, []) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=1.0 - ) - - with pytest.warns(UserWarning, match="during preequilibration"): - converted = assignments_to_experiments(problem) - - assert _changes_by_time(converted, "e1") == { - v2.C.TIME_PREEQUILIBRATION: {"k1": 3.0}, - 0.0: {}, - 10.0: {"k1": 0.0}, - } - - -def test_assignments_to_experiments_species_target(): - """Test that an assignment to a species is rejected.""" - problem = Problem() - problem.model = SbmlModel.from_antimony(""" - model conversion - species A, S; - A = 10; - S := piecewise(0, time < 10, 1); - R1: A -> ; S * A; - end - """) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement("obs_a", time=20.0, measurement=1.0) - - with pytest.raises(NotImplementedError, match="not a parameter"): - assignments_to_experiments(problem) - - -def test_assignments_to_experiments_nothing_to_convert(): - """Test that a problem without time-dependent assignments is unchanged, - including an experiment that only preequilibrates and periods that are - not in time order.""" - problem = Problem() - problem.model = SbmlModel.from_antimony(""" - model conversion - species A; - A = 10; - k2 = 0.5; - R1: A -> ; k2 * A; - end - """) - problem.add_condition("c1", k2=1.0) - problem.add_condition("c2", k2=2.0) - problem.add_experiment("e1", v2.C.TIME_PREEQUILIBRATION, "c1") - # periods out of time order - problem.add_experiment("e2", 10, "c2", 0, "c1") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - for experiment_id in ("e1", "e2"): - problem.add_measurement( - "obs_a", experiment_id=experiment_id, time=20.0, measurement=1.0 - ) - - converted = assignments_to_experiments(problem) - - def periods(problem): - return [ - [(p.time, p.condition_ids) for p in e.periods] - for e in problem.experiments - ] - - assert periods(converted) == periods(problem) - assert converted.condition_df.equals(problem.condition_df) - - -@pytest.mark.parametrize("on", ["model", "species"]) -def test_assignments_to_experiments_conversion_factor(on): - """Test that the changes of a symbol of a converted assignment stay if - the model refers to that symbol as a conversion factor.""" - problem = Problem() - problem.model = _switch_model( - "cf * piecewise(0, time < 10, 1)", extra="cf = 1;" - ) - sbml_model = problem.model.sbml_model - element = sbml_model if on == "model" else sbml_model.getSpecies("A") - assert ( - element.setConversionFactor("cf") == libsbml.LIBSBML_OPERATION_SUCCESS - ) - problem.add_condition("c1", cf=2.0) - problem.add_experiment("e1", 0, "c1") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=1.0 - ) - - converted = assignments_to_experiments(problem) - - assert _changes_by_time(converted, "e1") == { - 0.0: {"cf": 2.0}, - 10.0: {"k1": 2.0}, - } - - -def test_assignments_to_experiments_slope_id_taken(): - """Test that a ramp's slope parameter gets an ID that the problem does not - use yet, here because an observable parameter is called `k1_slope`.""" - problem = Problem() - problem.model = _switch_model("piecewise(time, time < 10, 1)") - problem.add_observable("obs_a", formula="k1_slope * A", noise_formula=1.0) - problem.add_measurement("obs_a", time=20.0, measurement=5.0) - problem.add_parameter("k1_slope", estimate=True, lb=1e-3, ub=1e3) - - converted = assignments_to_experiments(problem) - - assert _rate_rule(converted, "k1") == "k1_slope_1" - assert _changes_by_time(converted, DEFAULT_EXPERIMENT_ID) == { - 0.0: {"k1": 0.0, "k1_slope_1": 1.0}, - 10.0: {"k1": 1.0, "k1_slope_1": 0.0}, - } - - -@pytest.mark.parametrize( - "conditions", - [ - # the changes of a condition are simultaneous: `d1` takes the old `b1` - [(2.0, {"b1": 2.0, "d1": "b1 + 0.5"})], - # a change is evaluated once: `d1` keeps the `b1` of its period - [(2.0, {"d1": "b1 + 0.5"}), (3.0, {"b1": 2.0})], - ], -) -def test_assignments_to_experiments_change_evaluated_at_period_start( - conditions, -): - """Test that a condition change's value is evaluated at the start of its - period, with the values from before the period, also when a later switch - uses it.""" - problem = Problem() - problem.model = _switch_model( - "piecewise(b0, time < 5, d1)", extra="b0 = 0.3; b1 = 1.2; d1 = 1;" - ) - args = [0, []] - for i, (t, changes) in enumerate(conditions): - problem.add_condition(f"c{i}", **changes) - args += [t, f"c{i}"] - problem.add_experiment("e1", *args) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=1.0 - ) - - converted = assignments_to_experiments(problem) - - (switch,) = (p for p in converted["e1"].periods if p.time == 5.0) - (change,) = converted.get_changes_for_period(switch) - assert change.target_id == "k1" - # 1.2 + 0.5, not 2.0 + 0.5 - assert float(change.target_value) == pytest.approx(1.7) - - -def test_assignments_to_experiments_change_refers_to_species(): - """Test that a change of a symbol of the assignment whose value refers to - a species, which is not known at the start of the period, is - rejected.""" - problem = Problem() - problem.model = _switch_model( - "piecewise(0, time < 5, d1)", extra="d1 = 1;" - ) - problem.add_condition("c1", d1="A") - problem.add_experiment("e1", 0, [], 2, "c1") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=1.0 - ) - - with pytest.raises(NotImplementedError, match=r"refer to \['A'\]"): - assignments_to_experiments(problem) - - -def test_assignments_to_experiments_petab_symbols(): - """Test that the values of new condition changes use PEtab's symbols, so - that they equal the same expressions parsed by PEtab.""" - problem = Problem() - problem.model = _switch_model( - "piecewise(0, time < 5, d1)", extra="d1 = 1;" - ) - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement("obs_a", time=20.0, measurement=1.0) - problem.add_parameter( - "d1", estimate=True, nominal_value=2.0, lb=0.1, ub=10 - ) - - converted = assignments_to_experiments(problem) - - assert [ - change.target_value - for condition in converted.conditions - for change in condition.changes - ] == [sympify_petab("d1")] - - -@pytest.mark.parametrize("estimate", [False, True]) -def test_assignments_to_experiments_switching_time_from_fixed_parameter( - estimate, -): - """Test that a condition may set a switching time to an expression in a - parameter that the parameter table fixes, which takes its nominal value, - but not in an estimated one.""" - problem = Problem() - problem.model = _switch_model( - "piecewise(0, time < t1, 1)", extra="t1 = 0; p = 1;" - ) - problem.add_condition("c1", t1="p + 1") - problem.add_experiment("e1", 0, "c1") - problem.add_observable("obs_a", formula="A", noise_formula=1.0) - problem.add_measurement( - "obs_a", experiment_id="e1", time=20.0, measurement=1.0 - ) - if estimate: - problem.add_parameter( - "p", estimate=True, nominal_value=4.0, lb=0.1, ub=10 - ) - with pytest.raises( - NotImplementedError, match="evaluate to a real number" - ): - assignments_to_experiments(problem) - return - problem.add_parameter("p", estimate=False, nominal_value=4.0) - - converted = assignments_to_experiments(problem) - - # the switch is at the nominal value 4 + 1, not the model value 1 + 1 - assert _changes_by_time(converted, "e1") == {0.0: {}, 5.0: {"k1": 1.0}}