diff --git a/petab/v2/converters.py b/petab/v2/converters.py index 33f6f550..91e3a811 100644 --- a/petab/v2/converters.py +++ b/petab/v2/converters.py @@ -4,10 +4,11 @@ import warnings from copy import deepcopy +from itertools import count import libsbml import sympy as sp -from sbmlmath import sbml_math_to_sympy, set_math +from sbmlmath import TimeSymbol, sbml_math_to_sympy, set_math from .core import ( Change, @@ -20,7 +21,7 @@ from .models._sbml_utils import add_sbml_parameter, check from .models.sbml_model import SbmlModel -__all__ = ["ExperimentsToSbmlConverter"] +__all__ = ["ExperimentsToSbmlConverter", "assignments_to_experiments"] class ExperimentsToSbmlConverter: @@ -632,3 +633,263 @@ 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. + + **Scope:** conversion is supported for SBML assignment rules that contain + a single piecewise function with exactly 3 arguments, where the 2nd + argument is an inequality that can be + re-arranged explicitly into `t />= ...`, and whose other symbols + are all truly constant model parameters. Anything else raises a + :class:`NotImplementedError`. + + .. 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 + parameter_ids = {p.id for p in problem.parameters} + + 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 + # periods that are dropped if all their conditions were replaced + droppable = { + id(p) for e in experiments for p in e.periods if p.condition_ids + } + for experiment in experiments: + if all(p.is_preequilibration for p in experiment.periods): + experiment += (period := ExperimentPeriod(time=0)) + droppable.add(id(period)) + + # (experiment, period start time, target ID, target value) + additions = [] + replaced_symbols = set() + converted_ids = [] + + for assignment in [ + *(r for r in sbml_model.getListOfRules() if r.isAssignment()), + *sbml_model.getListOfInitialAssignments(), + ]: + expr = sbml_math_to_sympy(assignment) + pws = [pw for pw in expr.atoms(sp.Piecewise) if pw.has(TimeSymbol)] + if not pws: + continue + target_id = assignment.getId() + pw = pws[0] + (time,) = pw.atoms(TimeSymbol) + switch = pw.args[0].cond + if not ( + # A single piecewise + len(pws) == 1 + # with "conditional" and "otherwise" pieces + and len(pw.args) == 2 + # and explicitly an "otherwise" piece, + and pw.args[1].cond == sp.true + # with a supported inequality condition + and isinstance(switch, (sp.Gt, sp.Ge, sp.Lt, sp.Le)) + # that involves the time variable, + and not expr.subs(switch, sp.Dummy()).has(time) + # and that changes value at maximally a single time point + and (slope := sp.diff(switch.gts - switch.lts, time)).is_number + # and this time point is identifiable. + and slope.is_nonzero + ): + raise NotImplementedError( + "Only time dependence via a single 3-argument piecewise " + "function with a condition linear in time is supported, but " + f"the assignment to `{target_id}` is `{expr}`." + ) + if non_constant := sorted( + str(s) + for s in expr.free_symbols - {time} + if not (p := sbml_model.getParameter(str(s))) + or not p.getConstant() + ): + raise NotImplementedError( + "Only constant model parameters are supported as symbols in " + f"the assignment to `{target_id}`, but it refers to " + f"{non_constant}." + ) + # `switch` holds after the switch iff `gts - lts` increases with time + value, otherwise = (piece.expr for piece in pw.args) + before, after = (otherwise, value) if slope > 0 else (value, otherwise) + t_switch = sp.solve(switch.gts - switch.lts, time)[0] + + default = expr.subs(pw, before) + after_expr = expr.subs(pw, after) + target = sbml_model.getParameter(target_id) + target.setValue(float(default)) + target.setConstant(False) + assignment.removeFromParentAndDelete() + symbol_ids = {str(s) for s in expr.free_symbols} + replaced_symbols |= symbol_ids + converted_ids.append(target_id) + + for experiment in experiments: + periods = experiment.sorted_periods + overrides = {} + for i, period in enumerate(periods): + changes = problem.get_changes_for_period(period) + overrides |= {c.target_id: c.target_value for c in changes} + subs = { + s: overrides[str(s)] + for s in expr.free_symbols + if str(s) in overrides + } + 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 the switching time of the " + f"assignment to `{target_id}` in experiment " + f"`{experiment.id}` is `{t}`." + ) + t = float(t) + + if period.is_preequilibration: + # the value at time 0, which respects whether the + # condition holds at a switch at time 0 + value = expr.subs(subs).subs(time, 0) + additions.append( + (experiment, period.time, target_id, value) + ) + continue + + t_end = periods[i + 1].time if i + 1 < len(periods) else sp.oo + switch_value = ( + after_expr.subs(subs) if period.time < t < t_end else None + ) + start_value = ( + default if period.time < t else after_expr + ).subs(subs) + # The model already has the pre-switch value at the start of + # the first period, and the switch can become the first + # period if it only refers to parameter table symbols. + redundant = ( + i == 0 + and (start_value - default).is_zero + and ( + switch_value is None + or {str(s) for s in switch_value.free_symbols} + <= parameter_ids + ) + ) + if not redundant and ( + i == 0 or {c.target_id for c in changes} & symbol_ids + ): + additions.append( + (experiment, period.time, target_id, start_value) + ) + if switch_value is not None: + additions.append((experiment, t, target_id, switch_value)) + + 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, + ) + + # drop changes that are covered by the new conditions + 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 replaced_symbols + ] + table.elements = [c for c in table.conditions if c.changes] + removed_condition_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_condition_ids + ] + + # add the new conditions to new or existing periods + new_conditions = {} + cond_counter = count() + for experiment, t, target_id, value in additions: + if (cond_id := new_conditions.get((target_id, value))) is None: + while (cond_id := f"cond_{next(cond_counter)}") in condition_ids: + pass + 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) + + 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 + ] or experiment.periods[:1] + + 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 + + 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}" + ) + + return problem diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 68685c37..86ee3f69 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -11,6 +11,10 @@ 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, + assignments_to_experiments, +) from petab.v2.petab1to2 import petab1to2, v1v2_observable_df @@ -611,3 +615,347 @@ 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 + assert converted.model.sbml_model.getParameter("k1").getValue() == before + + experiments_expected = pd.DataFrame( + data={ + v2.C.EXPERIMENT_ID: [DEFAULT_EXPERIMENT_ID], + v2.C.TIME: [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_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 pre-dose value is kept as a condition here + 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"], + } + ) + + # the factor outside the piecewise expression is preserved + conditions_expected = pd.DataFrame( + data={ + v2.C.CONDITION_ID: ["cond_0", "cond_1"], + v2.C.TARGET_ID: ["k1", "k1"], + v2.C.TARGET_VALUE: [0.0, "level"], + } + ) + + 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 + assert converted.model.sbml_model.getParameter("k1").getValue() == 0.0 + + experiments_expected = pd.DataFrame( + data={ + v2.C.EXPERIMENT_ID: ["experiment____c1", "experiment____c2"], + v2.C.TIME: [10.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 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) + piecewise(0, time < 20, 1)", "single"), + ("piecewise(0, time < 10, 1, time < 20, 2)", "single"), + ("piecewise(0, time == 10, 1)", "single"), + ("piecewise(0, time^2 < 10, 1)", "single"), + ("piecewise(time, time < 10, 1)", "single"), + ("piecewise(0, time < t_dose, 1)", "evaluate to a real number"), + ("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("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 conversion makes the estimated + parameter `k1` a condition table target. + """ + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + k1 = piecewise(0, 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_measurement("obs_a", time=20.0, measurement=5.0) + problem.add_parameter("k1", 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)