diff --git a/petab/v2/petab1to2.py b/petab/v2/petab1to2.py index cd978e58..3201b94a 100644 --- a/petab/v2/petab1to2.py +++ b/petab/v2/petab1to2.py @@ -4,6 +4,7 @@ import shutil import warnings +from collections import defaultdict from contextlib import suppress from pathlib import Path from tempfile import TemporaryDirectory @@ -11,9 +12,12 @@ from uuid import uuid4 import pandas as pd +import sympy as sp from pandas.io.common import get_handle, is_url +from sbmlmath import TimeSymbol, sbml_math_to_sympy from .. import v1, v2 +from ..v1.math import petab_math_str, sympify_petab from ..v1.yaml import get_path_prefix, load_yaml from ..v1.yaml import validate as validate_yaml from ..versions import get_major_version @@ -25,6 +29,7 @@ def petab1to2( yaml_config: Path | str, output_dir: Path | str | None = None, + assignments_to_experiments: bool = False, *, validate: bool = True, ) -> v2.Problem | None: @@ -48,6 +53,10 @@ def petab1to2( :param output_dir: The output directory to save the converted PEtab problem, or ``None``, to return a :class:`petab.v2.Problem` instance. + :param assignments_to_experiments: + Whether to replace time-dependent piecewise assignments in the model, + such as ``k := dose * piecewise(0, time < t_dose, 1)``, by + experiment periods and conditions. :param validate: Whether to lint the input PEtab v1 problem, and the output PEtab v2 problem. @@ -60,16 +69,22 @@ def petab1to2( does not pass linting or if the generated files do not pass linting. """ if output_dir is not None: - return petab_files_1to2(yaml_config, output_dir, validate=validate) + return petab_files_1to2( + yaml_config, output_dir, assignments_to_experiments, + validate=validate + ) with TemporaryDirectory() as tmp_dir: - petab_files_1to2(yaml_config, tmp_dir, validate=validate) + petab_files_1to2( + yaml_config, tmp_dir, assignments_to_experiments, validate=validate + ) return v2.Problem.from_yaml(Path(tmp_dir, Path(yaml_config).name)) def petab_files_1to2( yaml_config: Path | str | dict, output_dir: Path | str, + assignments_to_experiments: bool = False, *, validate: bool = True, ): @@ -80,6 +95,8 @@ def petab_files_1to2( The PEtab problem as dictionary or YAML file name. :param output_dir: The output directory to save the converted PEtab problem. + :param assignments_to_experiments: + See :func:`petab1to2`. :param validate: Whether to lint the PEtab v1 problem before and the generated PEtab v2 problem after the conversion. See :func:`petab1to2`. @@ -112,6 +129,11 @@ def petab_files_1to2( ) if validate and v1.lint_problem(petab_problem): raise ValueError("Provided PEtab problem does not pass linting.") + new_periods = ( + _assignments_to_periods(petab_problem) + if assignments_to_experiments + else {} + ) output_dir = Path(output_dir) @@ -131,7 +153,10 @@ def petab_files_1to2( for file in ( model.location for model in new_yaml_config.model_files.values() ): - _copy_file(get_src_path(file), Path(get_dest_path(file))) + if new_periods: + petab_problem.model.to_file(get_dest_path(file)) + else: + _copy_file(get_src_path(file), Path(get_dest_path(file))) # Update observable table for observable_file in new_yaml_config.observable_files: @@ -144,6 +169,8 @@ def petab_files_1to2( # Update condition table for condition_file in new_yaml_config.condition_files: condition_df = v1.get_condition_df(get_src_path(condition_file)) + if new_periods: + condition_df = petab_problem.condition_df condition_df = v1v2_condition_df(condition_df, petab_problem.model) v2.write_condition_df(condition_df, get_dest_path(condition_file)) @@ -162,9 +189,11 @@ def v2_condition_id(cond_id: str) -> str: return "" def create_experiment_id(sim_cond_id: str, preeq_cond_id: str) -> str: - if not sim_cond_id and not preeq_cond_id: - return "" - if not preeq_cond_id and not v2_condition_id(sim_cond_id): + if ( + not preeq_cond_id + and not v2_condition_id(sim_cond_id) + and (sim_cond_id, preeq_cond_id) not in new_periods + ): # without pre-equilibration and without any changes, there is # nothing to be described by an experiment return "" @@ -208,6 +237,16 @@ def create_experiment_id(sim_cond_id: str, preeq_cond_id: str) -> str: v2.C.CONDITION_ID: v2_condition_id(sim_cond_id), } ) + experiments.extend( + {v2.C.EXPERIMENT_ID: exp_id, v2.C.TIME: t, v2.C.CONDITION_ID: c} + for t, c in new_periods.get((sim_cond_id, preeq_cond_id), []) + ) + if new_periods: + # drop periods of conditions whose changes were all replaced + has_changes = petab_problem.condition_df.notna().any(axis=1) + experiments = [ + e for e in experiments if has_changes[e[v2.C.CONDITION_ID]] + ] if experiments: exp_table_path = output_dir / "experiments.tsv" if exp_table_path.exists(): @@ -347,6 +386,101 @@ def _copy_file(src: Path | str, dest: Path): shutil.copy(str(src), str(dest)) +def _assignments_to_periods( + problem: v1.Problem, +) -> dict[tuple[str, str], list[tuple[float, str]]]: + """Replace time-dependent piecewise assignments by experiment periods. + + Assignment rules and initial assignments containing + ``piecewise(before, time < t_switch, after)`` are removed from the model, + which keeps the value before the switch. The condition table columns they + depend on are replaced by new conditions that set the assignment target + at the start of the respective periods. + + :returns: The new ``(time, condition ID)`` periods for each measured + ``(simulation condition ID, preequilibration condition ID)`` pair. + """ + sbml_model = problem.model.sbml_model + condition_df = problem.condition_df + overrides = { + cid: row.dropna().to_dict() for cid, row in condition_df.iterrows() + } + pairs = list( + problem.get_simulation_conditions_from_measurement_df() + .reindex( + columns=[ + v1.C.SIMULATION_CONDITION_ID, + v1.C.PREEQUILIBRATION_CONDITION_ID, + ], + fill_value="", + ) + .itertuples(index=False) + ) + parameter_ids = set(problem.parameter_df.index) + periods = defaultdict(list) + # (target ID, target value) -> new condition ID + new_conditions = {} + replaced_columns = 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 + (pw,) = pws + (time,) = pw.atoms(TimeSymbol) + (before, switch), _ = pw.args + t_switch = sp.solve(switch.lhs - switch.rhs, time)[0] + + target_id = assignment.getId() + target = sbml_model.getParameter(target_id) + target.setValue(float(expr.subs(pw, before))) + target.setConstant(False) + assignment.removeFromParentAndDelete() + replaced_columns |= {str(s) for s in expr.free_symbols} + + for sim_id, preeq_id in pairs: + t_sim = float(t_switch.subs(overrides[sim_id])) + values = { + t: expr.subs(overrides[sim_id]).subs(time, t) + for t in (0.0, t_sim) + } + if preeq_id: + values[v2.C.TIME_PREEQUILIBRATION] = expr.subs( + overrides[preeq_id] + ).subs(time, 0) + elif {str(s) for s in values[t_sim].free_symbols} <= parameter_ids: + # the model already has the pre-switch value, and the switch + # is valid as first period (only parameter table symbols) + del values[0.0] + for t, value in values.items(): + cond_id = new_conditions.setdefault( + (target_id, value), f"cond_{len(new_conditions)}" + ) + periods[sim_id, preeq_id].append((t, cond_id)) + + problem.condition_df = pd.concat( + [ + condition_df.drop(columns=replaced_columns, errors="ignore"), + pd.DataFrame.from_dict( + { + cid: { + target_id: float(value) + if value.is_number + else petab_math_str(value) + } + for (target_id, value), cid in new_conditions.items() + }, + orient="index", + ), + ] + ).rename_axis(v1.C.CONDITION_ID) + return {pair: sorted(p) for pair, p in periods.items()} + + def v1v2_condition_df( condition_df: pd.DataFrame, model: v1.Model ) -> pd.DataFrame: diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index d2d33bb1..5404b22e 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -5,8 +5,9 @@ import pytest from petab import v1, v2 -from petab.v1.models.sbml_model import SbmlModel +from petab.v1 import Problem as v1Problem from petab.v2 import Problem +from petab.v2.models.sbml_model import SbmlModel from petab.v2.petab1to2 import petab1to2, v1v2_observable_df @@ -383,3 +384,221 @@ def test_benchmark_collection(problem_id): pytest.skip(str(e)) assert isinstance(problem, Problem) assert len(problem.measurements) + + +def test_petab1to2_assignments_to_experiments(tmp_path): + problem = v1Problem() + 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_condition("c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", "c1", time=0.0, measurement=10.0) + problem.add_measurement("obs_a", "c1", time=20.0, measurement=5.0) + problem.add_parameter("k2", estimate=True, scale="lin", lb=1e-3, ub=1e3) + + yaml_path = problem.to_files_generic(prefix_path=tmp_path) + + converted = petab1to2(yaml_path, assignments_to_experiments=True) + + # 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: ["experiment____c1"], + 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() == [ + "experiment____c1", + "experiment____c1", + ] + + +def test_petab1to2_assignments_to_experiments_nested(tmp_path): + """Test conversion of piecewise assignments that are nested inside a + larger expression and switch at a time given by some expression.""" + problem = v1Problem() + 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_condition("c1") + problem.add_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", "c1", time=0.0, measurement=10.0) + problem.add_parameter("k2", estimate=True, scale="lin", lb=1e-3, ub=1e3) + + yaml_path = problem.to_files_generic(prefix_path=tmp_path) + + converted = petab1to2(yaml_path, assignments_to_experiments=True) + + 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: ["experiment____c1", "experiment____c1"], + 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_petab1to2_assignments_to_experiments_condition_specific(tmp_path): + """Test conversion of piecewise assignments whose switching times and + values are set by the condition table.""" + problem = v1Problem() + 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_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement("obs_a", "c1", time=30.0, measurement=10.0) + problem.add_measurement("obs_a", "c2", time=30.0, measurement=5.0) + problem.add_parameter("k2", estimate=True, scale="lin", lb=1e-3, ub=1e3) + + yaml_path = problem.to_files_generic(prefix_path=tmp_path) + + converted = petab1to2(yaml_path, assignments_to_experiments=True) + + # 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_petab1to2_assignments_to_experiments_preequilibration(tmp_path): + """Test that the converted assignments are applied during + preequilibration, using the preequilibration condition's values.""" + problem = v1Problem() + 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_observable("obs_a", formula="A", noise_formula=1.0) + problem.add_measurement( + "obs_a", "c1", time=30.0, measurement=10.0, preeq_cond_id="c2" + ) + problem.add_parameter("k2", estimate=True, scale="lin", lb=1e-3, ub=1e3) + + yaml_path = problem.to_files_generic(prefix_path=tmp_path) + + converted = petab1to2(yaml_path, assignments_to_experiments=True) + + # 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_2", "cond_0", "cond_1"], + } + ) + + 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: [0.0, 2.0, 3.0], + } + ) + + assert converted.experiment_df.equals(experiments_expected) + assert converted.condition_df.equals(conditions_expected)