diff --git a/petab/v2/converters.py b/petab/v2/converters.py index 33f6f550..99b2ed7a 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,182 @@ 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. + + Measurements without an experiment ID are assigned to a new experiment + :data:`DEFAULT_EXPERIMENT_ID`, if needed. + + :param problem: The PEtab problem to convert. This will not be modified. + :returns: The converted PEtab problem. + """ + 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() + + 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 + if len(pws) > 1 or len(pws[0].args) != 2: + raise NotImplementedError( + "Only assignments with a single time-dependent piecewise " + "function with a single condition are supported, but " + f"the assignment to `{assignment.getId()}` is `{expr}`." + ) + (pw,) = pws + (time,) = pw.atoms(TimeSymbol) + (before, switch), _ = pw.args + t_switch = sp.solve(switch.lhs - switch.rhs, time)[0] + + target_id = assignment.getId() + default = expr.subs(pw, before) + 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 + + 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 + } + cur_expr = expr.subs(subs) + if period.is_preequilibration: + additions.append( + ( + experiment, + period.time, + target_id, + cur_expr.subs(time, 0), + ) + ) + continue + + t_end = periods[i + 1].time if i + 1 < len(periods) else sp.oo + t = t_switch.subs(subs) + if not (t.is_number and t.is_real): + 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) + switch_value = ( + cur_expr.subs(time, t) if period.time < t < t_end else None + ) + start_value = cur_expr.subs(time, period.time) + # 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)) + + # 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 + + return problem diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 68685c37..762c6667 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -11,6 +11,7 @@ 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 assignments_to_experiments from petab.v2.petab1to2 import petab1to2, v1v2_observable_df @@ -611,3 +612,343 @@ def test_benchmark_collection(problem_id): pytest.skip(str(e)) assert isinstance(problem, Problem) assert len(problem.measurements) + + +def test_assignments_to_experiments(): + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + k1 = piecewise(1, time < 10, 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() == 1.0 + + experiments_expected = pd.DataFrame( + data={ + v2.C.EXPERIMENT_ID: ["default_experiment"], + 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: [2.0], + } + ) + + 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", + "default_experiment", + ] + + +def test_assignments_to_experiments_reversed(): + problem = Problem() + ant_model = """ + model conversion + species A, B; + A = 10; + B = 0; + k1 = piecewise(1, 10 > time, 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() == 1.0 + + experiments_expected = pd.DataFrame( + data={ + v2.C.EXPERIMENT_ID: ["default_experiment"], + 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: [2.0], + } + ) + + 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", + "default_experiment", + ] + + +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"] * 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) + + 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", + [ + "piecewise(0, time < 10, 1) + piecewise(0, time < 20, 1)", + "piecewise(0, time < 10, 1, time < 20, 2)", + ], +) +def test_assignments_to_experiments_multiple_switches(assignment): + """Test that assignments with more than one switch are rejected.""" + 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=0.0, measurement=10.0) + problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3) + + with pytest.raises(NotImplementedError, match="single time-dependent"): + assignments_to_experiments(problem) + + +@pytest.mark.parametrize( + "assignment, estimated", + [ + # depends on a species + ("piecewise(0, time < A, 1)", None), + ("piecewise(0, A * time < 5, 1)", None), + # depends on an estimated parameter + ("piecewise(0, time < t_dose, 1)", "t_dose"), + # no real solution + ("piecewise(0, time^2 + 1 < 0, 1)", None), + ], +) +def test_assignments_to_experiments_non_numeric_switch(assignment, estimated): + """Test that switching times that are not real numbers 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) + if estimated: + problem.add_parameter(estimated, estimate=True, lb=0, ub=10) + + with pytest.raises(NotImplementedError, match="evaluate to a real number"): + assignments_to_experiments(problem)