diff --git a/CHANGELOG.md b/CHANGELOG.md index 79db7e5b..4bc4f5e5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -33,6 +33,9 @@ https://github.com/PEtab-dev/libpetab-python/pull/533, https://github.com/PEtab-dev/libpetab-python/pull/535, and https://github.com/PEtab-dev/libpetab-python/pull/536) +* `petab1to2` no longer adds default `uniform` priors to non-estimated + parameters, or if no estimated parameter has an objective prior + (by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/539) * `petab_math_str` now produces valid PEtab math expressions for boolean operators, (in)equality, `min`/`max`, powers with negative bases, and `pi`/`e`, and prints floats without loss of precision. PEtab v2 condition @@ -98,6 +101,10 @@ * `petab.v2.reparametrize_log` estimates parameters on log scale by replacing them with new log-transformed parameters (by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/560) +* `petab1to2` and `petab_files_1to2` can preserve the objective function of + log-scaled parameters with `parameterScale*` priors via + `reparametrize="priors"` or `reparametrize="all"` + (by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/539) ## 0.9 series diff --git a/petab/v2/petab1to2.py b/petab/v2/petab1to2.py index b03cf714..cea4a73d 100644 --- a/petab/v2/petab1to2.py +++ b/petab/v2/petab1to2.py @@ -8,9 +8,11 @@ from contextlib import suppress from pathlib import Path from tempfile import TemporaryDirectory +from typing import Literal, get_args from urllib.parse import urlparse import pandas as pd +import sympy as sp from pandas.io.common import get_handle, is_url from .. import v1, v2 @@ -19,15 +21,25 @@ from ..v1.yaml import validate as validate_yaml from ..versions import get_major_version from .models import MODEL_TYPE_SBML +from .reparametrize import ( + _get_override_slots, + _log, + _log10, + reparametrize_log, +) __all__ = ["petab1to2"] +#: Options for reparametrizing log-scaled parameters +Reparametrize = Literal["none", "priors", "all"] + def petab1to2( yaml_config: Path | str, output_dir: Path | str | None = None, *, validate: bool = True, + reparametrize: Reparametrize = "none", ) -> v2.Problem | None: """Convert from PEtab 1.0 to PEtab 2.0 format. @@ -55,16 +67,51 @@ def petab1to2( If ``True`` (default), both problems are validated. If ``False``, neither problem is validated. The ``yaml_config`` content is validated regardless. + :param reparametrize: + Which ``log``/``log10``-scaled parameters to estimate as + ``log10_``/``log_`` (see :func:`petab.v2.reparametrize_log`). + PEtab v2 has no ``parameterScale``, so without reparametrization, + ``parameterScale*`` priors are converted to densities over the + unscaled parameter, which changes the maximum a posteriori estimate. + + * ``"none"`` (default): don't reparametrize. + * ``"priors"``: estimated parameters with ``parameterScale*`` priors. + * ``"all"``: all estimated log-scaled parameters, except those with + priors over the unscaled parameter (e.g., ``logNormal``). + This also preserves the v1 default prior (uniform on parameter + scale). Parameters without ``parameterScale*`` priors are not + reparametrized (with a warning) if they override a placeholder + that cannot be log-transformed. + + Parameters with a lower bound of ``0`` get a lower bound of ``-inf``. + They are not reparametrized if their prior on log scale would be + an improper uniform prior. + Further log-scaled parameters are reparametrized if they share a + placeholder with reparametrized ones in the measurement table. :raises ValueError: If the input is invalid or, if ``validate`` is ``True``, if the input does not pass linting or if the generated files do not pass linting. + :raises NotImplementedError: + If a parameter with a ``parameterScale*`` prior cannot be + reparametrized, e.g., because it overrides a placeholder that + cannot be log-transformed. """ if output_dir is not None: - return petab_files_1to2(yaml_config, output_dir, validate=validate) + return petab_files_1to2( + yaml_config, + output_dir, + validate=validate, + reparametrize=reparametrize, + ) with TemporaryDirectory() as tmp_dir: - petab_files_1to2(yaml_config, tmp_dir, validate=validate) + petab_files_1to2( + yaml_config, + tmp_dir, + validate=validate, + reparametrize=reparametrize, + ) return v2.Problem.from_yaml(Path(tmp_dir, Path(yaml_config).name)) @@ -73,6 +120,7 @@ def petab_files_1to2( output_dir: Path | str, *, validate: bool = True, + reparametrize: Reparametrize = "none", ): """Convert PEtab files from PEtab 1.0 to PEtab 2.0. @@ -84,6 +132,9 @@ def petab_files_1to2( :param validate: Whether to lint the PEtab v1 problem before and the generated PEtab v2 problem after the conversion. See :func:`petab1to2`. + :param reparametrize: + Which log-scaled parameters to reparametrize. + See :func:`petab1to2`. :raises ValueError: If the input is invalid or, if ``validate`` is ``True``, if the input @@ -101,6 +152,12 @@ def petab_files_1to2( get_dest_path = lambda filename: f"{output_dir}/{filename}" # noqa: E731 + if reparametrize not in get_args(Reparametrize): + raise ValueError( + f"Invalid value for `reparametrize': {reparametrize}. " + f"Must be one of {get_args(Reparametrize)}." + ) + # Validate the original PEtab problem validate_yaml(yaml_config, path_prefix=path_prefix) if get_major_version(yaml_config) != 1: @@ -123,7 +180,9 @@ def petab_files_1to2( # Update tables # parameter table - parameter_df = v1v2_parameter_df(petab_problem.parameter_df.copy()) + parameter_df = v1v2_parameter_df( + petab_problem.parameter_df.copy(), reparametrize=reparametrize + ) v2.write_parameter_df( parameter_df, get_dest_path(new_yaml_config.parameter_files[0]) ) @@ -289,6 +348,15 @@ def create_experiment_id(sim_cond_id: str, preeq_cond_id: str) -> str: new_yaml_file = output_dir / Path(yaml_file).name new_yaml_config.to_yaml(new_yaml_file) + if reparametrize != "none" and ( + problem := _reparametrize_log_scaled( + v2.Problem.from_yaml(new_yaml_file), + petab_problem.parameter_df, + reparametrize, + ) + ): + problem.to_files(output_dir) + if not validate: return @@ -519,13 +587,32 @@ def extract_placeholders(row: pd.Series, type_: str) -> str: return df +def _has_objective_priors(parameter_df: pd.DataFrame) -> bool: + """Whether any estimated parameter of a v1 parameter table has an + objective prior.""" + if v1.C.OBJECTIVE_PRIOR_TYPE not in parameter_df: + return False + return bool( + ( + parameter_df[v1.C.OBJECTIVE_PRIOR_TYPE].notna() + & (parameter_df[v1.C.ESTIMATE].astype(int) == 1) + ).any() + ) + + def v1v2_parameter_df( parameter_df: pd.DataFrame, + reparametrize: Reparametrize = "none", ) -> pd.DataFrame: """Convert parameter table from petab v1 to v2. Do all the necessary conversions to the parameter table that can be done with the parameter table alone. + + :param parameter_df: The PEtab v1 parameter table. + :param reparametrize: + The reparametrization option passed to :func:`petab1to2`. + Only used for warnings, the reparametrization itself is done later. """ df = parameter_df.copy().reset_index() @@ -534,10 +621,16 @@ def v1v2_parameter_df( lambda x: str(bool(int(x))).lower() ) - def update_prior(row) -> tuple[str, str | None]: + has_priors = _has_objective_priors(parameter_df) + + def update_prior(row) -> tuple[str | None, str | None]: """Convert prior type and prior parameters to v2 format.""" prior_type = row.get(v1.C.OBJECTIVE_PRIOR_TYPE) if pd.isna(prior_type): + if not has_priors or row[v2.C.ESTIMATE] != "true": + # no prior + return None, None + # the v2 default if other parameters have priors prior_type = v1.C.UNIFORM prior_pars = row.get(v1.C.OBJECTIVE_PRIOR_PARAMETERS) @@ -628,10 +721,35 @@ def update_prior(row) -> tuple[str, str | None]: stacklevel=9, ) if not (df[v1.C.PARAMETER_SCALE] == v1.C.LIN).all(): - warnings.warn( - "Parameter scales are not supported in PEtab v2.", - stacklevel=9, + msg = "Parameter scales are not supported in PEtab v2." + reparametrizable = ( + parameter_df[v1.C.PARAMETER_SCALE].isin([v1.C.LOG, v1.C.LOG10]) + & (parameter_df[v1.C.ESTIMATE].astype(int) == 1) + & (parameter_df[v1.C.LOWER_BOUND] >= 0) ) + if reparametrize == "none" and reparametrizable.any(): + prior_types = parameter_df.get(v1.C.OBJECTIVE_PRIOR_TYPE) + scale_priors = ( + [] + if prior_types is None + else parameter_df.index[ + reparametrizable + & prior_types.isin(v1.C.PARAMETER_SCALE_PRIOR_TYPES) + ] + ) + msg += " Log-scaled parameters are estimated on linear scale." + if len(scale_priors): + msg += ( + f" Parameter-scale priors of {', '.join(scale_priors)} " + "were converted to the corresponding priors of the " + "unscaled parameters, which changes the maximum a " + "posteriori estimate." + ) + msg += ( + " Use `reparametrize='priors'` or `'all'` to estimate them " + "on log scale." + ) + warnings.warn(msg, stacklevel=9) df = df.drop( columns=[ v1.C.INITIALIZATION_PRIOR_TYPE, @@ -659,3 +777,177 @@ def update_prior_pars(row): df[v2.C.PRIOR_PARAMETERS] = df.apply(update_prior_pars, axis=1) return df + + +def _reparametrize_log_scaled( + problem: v2.Problem, + v1_parameter_df: pd.DataFrame, + reparametrize: Reparametrize, +) -> v2.Problem | None: + """Reparametrize log-scaled parameters of a converted problem. + + See :func:`petab1to2`. + + :param problem: The converted problem. + :param v1_parameter_df: The original PEtab v1 parameter table. + :param reparametrize: The parameters to reparametrize. + :return: + The reparametrized problem, or ``None`` if no parameter was + reparametrized. + """ + df = v1_parameter_df + has_priors = _has_objective_priors(df) + # for each base: candidates and their priors over the scaled parameter + candidates = {10: {}, "e": {}} + # parameters with parameter-scale priors + required = {10: set(), "e": set()} + selected = {10: set(), "e": set()} + # selected parameters that cannot be reparametrized due to their bounds + # and priors + zero_bound = [] + # non-estimated parameters whose priors cannot be converted + no_prior = set() + for par_id, row in df.iterrows(): + pscale = row[v1.C.PARAMETER_SCALE] + if pscale not in (v1.C.LOG, v1.C.LOG10): + continue + base = 10 if pscale == v1.C.LOG10 else "e" + log_function = _log10 if base == 10 else _log + prior_type = row.get(v1.C.OBJECTIVE_PRIOR_TYPE) + if pd.isna(prior_type): + # v1 default + prior_type = None + elif prior_type not in v1.C.PARAMETER_SCALE_PRIOR_TYPES: + # A prior over the unscaled parameter. Its conversion already + # preserves the objective function. + continue + estimate = int(row[v1.C.ESTIMATE]) + is_selected = estimate and ( + reparametrize == "all" or prior_type is not None + ) + if estimate and ( + row[v1.C.LOWER_BOUND] < 0 or row[v1.C.UPPER_BOUND] <= 0 + ): + # No valid bounds on log scale. A negative lower bound is invalid + # in v1; an upper bound of 0 implies `lb == ub == 0`. + continue + + # The prior over the scaled parameter. Derived from the v1 table + # rather than from the converted v2 prior to avoid rounding errors. + # Without a prior type, but with other estimated parameters having + # priors, the v1 default is made explicit, as in + # `v1v2_parameter_df`. + prior = None + if prior_type is not None or (estimate and has_priors): + lb, ub = row[v1.C.LOWER_BOUND], row[v1.C.UPPER_BOUND] + if lb >= 0 and ub > 0: + # v1 allows a lower bound of 0 for log-scaled parameters + lb, ub = log_function(lb), log_function(ub) + else: + # missing or invalid bounds of a non-estimated parameter + lb, ub = -math.inf, math.inf + pars = row.get(v1.C.OBJECTIVE_PRIOR_PARAMETERS) + if prior_type is None or pd.isna(pars): + pars = [lb, ub] + else: + pars = [ + float(x) for x in str(pars).split(v1.C.PARAMETER_SEPARATOR) + ] + if prior_type in (None, v1.C.PARAMETER_SCALE_UNIFORM): + prior = (v2.C.UNIFORM, [max(pars[0], lb), min(pars[1], ub)]) + else: + name = prior_type.removeprefix("parameterScale").lower() + prior = (name, pars) + if prior[0] == v2.C.UNIFORM and not all( + map(math.isfinite, prior[1]) + ): + if estimate: + # improper prior for a lower bound of 0 + if is_selected: + zero_bound.append(par_id) + continue + prior = None + no_prior.add(par_id) + + candidates[base][par_id] = prior + if estimate and prior_type is not None: + required[base].add(par_id) + if is_selected: + selected[base].add(par_id) + + if zero_bound: + warnings.warn( + f"Log-scaled parameters {', '.join(zero_bound)} were not " + "reparametrized, because their lower bound of 0 would result in " + "an improper prior on log scale. " + "They are estimated on linear scale, and any parameter-scale " + "priors were converted to the corresponding priors of the " + "unscaled parameters, which changes the maximum a posteriori " + "estimate.", + stacklevel=4, + ) + + # Placeholders can only be transformed if all their overrides are + # reparametrized or non-negative numbers. Add log-scaled parameters that + # share a placeholder with selected parameters. If that is not possible, + # drop selected parameters that don't have parameter-scale priors. + # Otherwise, `reparametrize_log` will raise. + slots = _get_override_slots(problem).values() + excluded = set() + for base, base_selected in selected.items(): + changed = True + while changed: + changed = False + for values in slots: + ids = {str(v) for v in values if isinstance(v, sp.Symbol)} + if not ids & base_selected: + continue + other = ids - base_selected + if other and other <= candidates[base].keys() - excluded: + base_selected |= other + changed = True + elif other or any(v.is_number and v < 0 for v in values): + optional = (ids & base_selected) - required[base] + if optional: + base_selected -= optional + excluded |= optional + changed = True + if excluded: + warnings.warn( + "Log-scaled parameters " + f"{', '.join(sorted(excluded))} were not reparametrized, " + "because they are used as overrides for placeholders that " + "cannot be log-transformed.", + stacklevel=4, + ) + + if not any(selected.values()): + return None + + if no_prior := no_prior & set().union(*selected.values()): + warnings.warn( + "The priors of the non-estimated parameters " + f"{', '.join(sorted(no_prior))} were dropped, because they " + "cannot be converted to log scale without valid bounds.", + stacklevel=4, + ) + + for base, base_selected in selected.items(): + if not base_selected: + continue + try: + problem, _ = reparametrize_log( + problem, + # keep the parameter table order + [p for p in candidates[base] if p in base_selected], + base=base, + priors={p: candidates[base][p] for p in base_selected}, + ) + except NotImplementedError as e: + raise NotImplementedError( + "Log-scaled parameters could not be reparametrized. " + "Use `reparametrize='none'` to estimate them on linear " + "scale.\n" + str(e) + ) from e + + return problem diff --git a/petab/v2/reparametrize.py b/petab/v2/reparametrize.py index 5dfbc130..7490e0c6 100644 --- a/petab/v2/reparametrize.py +++ b/petab/v2/reparametrize.py @@ -84,8 +84,10 @@ def reparametrize_log( * Measurement table: Overrides ``p`` are replaced by ``q``. The respective placeholder ``ph`` is substituted by ``base^ph`` in the observable (and noise) formula. This requires that all overrides for that - placeholder are reparametrized parameters or positive numbers. - Numbers are log-transformed. + placeholder are reparametrized parameters or non-negative numbers. + Numbers are log-transformed. Since overrides must be finite, ``0`` is + replaced by a new non-estimated parameter ``negative_infinity`` with + nominal value ``-inf``. The new ID is ``log10_

`` for base 10 and ``log_

`` for base e, with a numeric suffix if necessary to make it unique. @@ -356,8 +358,10 @@ def _reparametrize_overrides( A placeholder ``ph`` is substituted by ``b^ph`` if any of its overrides is a reparametrized parameter. This requires all its overrides (across all measurements of the observable) to be reparametrized - parameters or positive numbers. Otherwise, the observable would have to - be split, which is not supported. + parameters or non-negative numbers. Otherwise, the observable would have + to be split, which is not supported. + Numeric overrides are log-transformed. Since overrides must be finite, + ``0`` is replaced by a new fixed parameter with nominal value ``-inf``. """ observables = {o.id: o for o in problem.observables} @@ -380,16 +384,14 @@ def quote(values) -> str: ) if not reparametrized: continue - non_positive = sorted( - float(v) for v in values if v.is_number and not v > 0 - ) + negative = sorted(float(v) for v in values if v.is_number and v < 0) others = sorted( str(v) for v in values if not v.is_number and not (isinstance(v, sp.Symbol) and v.name in id_map) ) - if not non_positive and not others: + if not negative and not others: slots_to_transform.append(slot) continue @@ -404,9 +406,9 @@ def quote(values) -> str: f"\n - parameters that are not reparametrized: " f"{quote(others)}" ) - if non_positive: - error += "\n - non-positive numbers: " + ", ".join( - f"{v:g}" for v in non_positive + if negative: + error += "\n - negative numbers: " + ", ".join( + f"{v:g}" for v in negative ) errors.append(error) @@ -415,7 +417,7 @@ def quote(values) -> str: "Cannot log-transform the following placeholders, which are " "overridden by parameters to be reparametrized. This requires " "all overrides of a placeholder to be reparametrized parameters " - "or positive numbers.\n" + "or non-negative numbers.\n" + "\n".join(errors) + "\nSplitting observables is not supported. Reparametrize " "all or none of the parameters overriding a placeholder, " @@ -439,9 +441,21 @@ def quote(values) -> str: ) # update override values + neg_inf_id = None + def transform(value: sp.Basic) -> sp.Basic: + nonlocal neg_inf_id if isinstance(value, sp.Symbol): return petab_symbol(id_map[value.name]) + if value.is_zero: + if neg_inf_id is None: + neg_inf_id = problem.get_unused_id( + "negative_infinity", exclude=id_map.values() + ) + problem.add_parameter( + neg_inf_id, estimate=False, nominal_value=-math.inf + ) + return petab_symbol(neg_inf_id) return sp.Float(log_function(float(value))) slots_to_transform = set(slots_to_transform) diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 1a2a528f..bb6c3628 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -12,6 +12,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.math import petab_symbol from petab.v2.petab1to2 import petab1to2, v1v2_observable_df @@ -84,12 +85,30 @@ def test_v1v2_observable_df_log10(): def _write_v1_problem( path: Path, *, - observables: list[dict], - measurements: list[dict], + observables: list[dict] | None = None, + measurements: list[dict] | None = None, parameters: list[dict], ) -> str: """Write a PEtab v1 problem for the conversion A -> B with the given - tables (as lists of rows) to ``path`` and return the YAML file name.""" + tables (as lists of rows) to ``path`` and return the YAML file name. + + By default, there is one measurement of the observable ``obs_B = B``. + """ + observables = observables or [ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.OBSERVABLE_FORMULA: "B", + v1.C.NOISE_FORMULA: "0.1", + } + ] + measurements = measurements or [ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.SIMULATION_CONDITION_ID: "c0", + v1.C.TIME: 1.0, + v1.C.MEASUREMENT: 0.3, + } + ] problem = v1.Problem( model=SbmlModel.from_antimony( "r1: A -> B; k1 * A; A = 1; B = 0; k1 = 0.5" @@ -260,6 +279,556 @@ def test_petab1to2_priors( assert v2_prior.pdf(x) == pytest.approx(v1_prior.pdf(x), rel=1e-12) +@pytest.mark.parametrize( + "estimate_k1, sigma_prior", [(0, False), (1, False), (1, True)] +) +def test_petab1to2_default_priors(tmp_path, estimate_k1, sigma_prior): + """The default prior is only made explicit for estimated parameters, + and only if any estimated parameter has a prior. Priors of + non-estimated parameters are kept.""" + fixed = { + v1.C.PARAMETER_SCALE: v1.C.LIN, + v1.C.LOWER_BOUND: np.nan, + v1.C.UPPER_BOUND: np.nan, + v1.C.NOMINAL_VALUE: 0.05, + v1.C.ESTIMATE: 0, + } + sigma = _K1 | {v1.C.PARAMETER_ID: "sigma"} + if sigma_prior: + sigma |= { + v1.C.OBJECTIVE_PRIOR_TYPE: v1.C.NORMAL, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "0.1;1", + } + yaml_file = _write_v1_problem( + tmp_path, + observables=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.OBSERVABLE_FORMULA: "B + offset1 + offset2", + v1.C.NOISE_FORMULA: "sigma", + } + ], + parameters=[ + _K1 | {v1.C.ESTIMATE: estimate_k1}, + sigma, + # (the prior is ignored for non-estimated parameters) + fixed + | { + v1.C.PARAMETER_ID: "offset1", + v1.C.OBJECTIVE_PRIOR_TYPE: v1.C.NORMAL, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "0;1", + }, + fixed | {v1.C.PARAMETER_ID: "offset2"}, + ], + ) + problem = petab1to2(yaml_file) + + if estimate_k1 and sigma_prior: + assert problem["k1"].prior_distribution == v2.C.UNIFORM + assert problem["sigma"].prior_distribution == v2.C.NORMAL + else: + assert problem["k1"].prior_distribution is None + assert problem["sigma"].prior_distribution is None + assert problem["offset1"].prior_distribution == v2.C.NORMAL + assert problem["offset1"].prior_parameters == [0, 1] + assert problem["offset2"].prior_distribution is None + + +@pytest.mark.filterwarnings( + "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning" +) +@pytest.mark.parametrize( + "parameter_scale, prior_type, prior_parameters, reparametrize", + [ + (v1.C.LOG10, v1.C.PARAMETER_SCALE_NORMAL, "-0.5;0.7", "priors"), + (v1.C.LOG10, v1.C.PARAMETER_SCALE_LAPLACE, "-0.5;0.7", "priors"), + (v1.C.LOG, v1.C.PARAMETER_SCALE_NORMAL, "-0.5;0.7", "priors"), + (v1.C.LOG, v1.C.PARAMETER_SCALE_LAPLACE, "-0.5;0.7", "all"), + (v1.C.LOG10, v1.C.PARAMETER_SCALE_UNIFORM, "-0.5;0.7", "priors"), + (v1.C.LOG10, v1.C.PARAMETER_SCALE_UNIFORM, "-400;400", "priors"), + (v1.C.LOG10, v1.C.PARAMETER_SCALE_UNIFORM, np.nan, "all"), + (v1.C.LOG, v1.C.PARAMETER_SCALE_UNIFORM, "-1000;1000", "priors"), + (v1.C.LOG10, np.nan, np.nan, "priors"), + (v1.C.LOG10, np.nan, np.nan, "all"), + ], +) +def test_petab1to2_reparametrize_priors( + tmp_path, parameter_scale, prior_type, prior_parameters, reparametrize +): + """Test that reparametrization preserves parameter-scale priors, + including their normalization constant up to truncation.""" + k1 = _K1 | { + v1.C.PARAMETER_SCALE: parameter_scale, + v1.C.OBJECTIVE_PRIOR_TYPE: prior_type, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: prior_parameters, + } + yaml_file = _write_v1_problem( + tmp_path, + parameters=[k1], + ) + problem = petab1to2(yaml_file, reparametrize=reparametrize) + + if pd.isna(prior_type) and reparametrize == "priors": + # the default prior is only reparametrized with "all" + assert [p.id for p in problem.parameters] == ["k1"] + return + + new_id = "log10_k1" if parameter_scale == v1.C.LOG10 else "log_k1" + assert [p.id for p in problem.parameters] == [new_id] + par = problem[new_id] + lb, ub = ( + v1.parameters.scale(x, parameter_scale) + for x in (_K1[v1.C.LOWER_BOUND], _K1[v1.C.UPPER_BOUND]) + ) + assert (par.lb, par.ub) == pytest.approx((lb, ub)) + assert par.nominal_value == pytest.approx( + v1.parameters.scale(_K1[v1.C.NOMINAL_VALUE], parameter_scale) + ) + + # k1 = b^q in the model + sbml_model = problem.model.sbml_model + assert sbml_model.getInitialAssignmentBySymbol("k1") is not None + assert sbml_model.getParameter(new_id).getConstant() + + if pd.isna(prior_type): + # no priors at all + assert par.prior_distribution is None + return + + # the v2 prior of the new parameter is the v1 prior density over the + # scaled parameter (both truncated to the bounds) + v1_prior = Prior.from_par_dict(k1, type_="objective") + if prior_type == v1.C.PARAMETER_SCALE_UNIFORM and not pd.isna( + prior_parameters + ): + # evaluate within the support + low, high = map(float, prior_parameters.split(";")) + lb, ub = max(lb, low), min(ub, high) + theta = np.linspace(lb, ub, 13)[1:-1] + v1_pdf = v1_prior.pdf(theta, x_scaled=True, rescale=True) + v2_pdf = par.prior_dist.pdf(theta) + if prior_type in ( + v1.C.PARAMETER_SCALE_NORMAL, + v1.C.PARAMETER_SCALE_LAPLACE, + ): + assert np.ptp(np.log(v1_pdf)) > 1 + else: + # uniform within the bounds + assert np.ptp(np.log(v1_pdf)) < 1e-12 + assert np.all(v1_pdf > 0) + np.testing.assert_allclose(v2_pdf, v1_pdf, rtol=1e-12) + + +@pytest.mark.filterwarnings( + "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning" +) +def test_petab1to2_reparametrize_llh(tmp_path): + """Test that reparametrizing overrides and parameters used in formulas + preserves the likelihood, for both bases.""" + + def par(par_id, scale, nominal_value): + return { + v1.C.PARAMETER_ID: par_id, + v1.C.PARAMETER_SCALE: scale, + v1.C.LOWER_BOUND: 1e-3, + v1.C.UPPER_BOUND: 1e3, + v1.C.NOMINAL_VALUE: nominal_value, + v1.C.ESTIMATE: 1, + v1.C.OBJECTIVE_PRIOR_TYPE: v1.C.PARAMETER_SCALE_NORMAL, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "0;1", + } + + yaml_file = _write_v1_problem( + tmp_path, + observables=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.OBSERVABLE_FORMULA: "observableParameter1_obs_B * B " + "+ offset", + v1.C.NOISE_FORMULA: "noiseParameter1_obs_B", + } + ], + measurements=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.SIMULATION_CONDITION_ID: "c0", + v1.C.TIME: float(t), + v1.C.MEASUREMENT: 0.3 + t, + v1.C.OBSERVABLE_PARAMETERS: obs_par, + v1.C.NOISE_PARAMETERS: noise_par, + } + for t, (obs_par, noise_par) in enumerate( + [("s1", "sd"), ("s2", "0.2"), ("4", "sd")] + ) + ], + parameters=[ + _K1, + par("s1", v1.C.LOG10, 3.0), + par("s2", v1.C.LOG10, 0.7), + par("sd", v1.C.LOG, 0.25), + par("offset", v1.C.LOG10, 0.13), + ], + ) + v1_problem = v1.Problem.from_yaml(yaml_file) + v2_problem = petab1to2(yaml_file, reparametrize="priors") + assert [p.id for p in v2_problem.parameters] == [ + "k1", + "log10_s1", + "log10_s2", + "log_sd", + "log10_offset", + ] + + def simulate(measurement_df): + return measurement_df.rename( + columns={v1.C.MEASUREMENT: v1.C.SIMULATION} + ).assign(**{v1.C.SIMULATION: [0.5, 0.6, 0.7]}) + + llh_v1 = calculate_llh_v1( + v1_problem.measurement_df, + simulate(v1_problem.measurement_df), + v1_problem.observable_df, + v1_problem.parameter_df, + ) + llh_v2 = calculate_llh_v2( + v2_problem.measurement_df, + simulate(v2_problem.measurement_df), + v2_problem.observable_df, + v2_problem.parameter_df, + ) + assert np.isfinite(llh_v1) + assert llh_v2 == pytest.approx(llh_v1, rel=1e-12) + + +def _write_v1_problem_zero_lower_bound( + tmp_path, prior_type=None, prior_parameters=None, sigma_prior=False +) -> str: + """Write a v1 problem with a log10-scaled `k1` with lower bound 0. + + If `sigma_prior`, the noise parameter `sigma` has a prior. + """ + # v1 allows a lower bound of 0 for log-scaled parameters + # if there is an initialization prior + k1 = _K1 | { + v1.C.PARAMETER_SCALE: v1.C.LOG10, + v1.C.LOWER_BOUND: 0.0, + v1.C.INITIALIZATION_PRIOR_TYPE: v1.C.PARAMETER_SCALE_UNIFORM, + v1.C.INITIALIZATION_PRIOR_PARAMETERS: "-3;2", + } + if prior_type is not None: + k1[v1.C.OBJECTIVE_PRIOR_TYPE] = prior_type + k1[v1.C.OBJECTIVE_PRIOR_PARAMETERS] = prior_parameters + parameters = [k1] + observables = None + if sigma_prior: + parameters.append( + _K1 + | { + v1.C.PARAMETER_ID: "sigma", + v1.C.OBJECTIVE_PRIOR_TYPE: v1.C.NORMAL, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "0.1;1", + } + ) + observables = [ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.OBSERVABLE_FORMULA: "B", + v1.C.NOISE_FORMULA: "sigma", + } + ] + yaml_file = _write_v1_problem( + tmp_path, observables=observables, parameters=parameters + ) + assert not v1.lint_problem(v1.Problem.from_yaml(yaml_file)) + return yaml_file + + +@pytest.mark.filterwarnings( + "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning", + "ignore:Initialisation priors.*:UserWarning", +) +@pytest.mark.parametrize( + "prior_type,prior_parameters,reparametrize,expected_prior", + [ + ( + v1.C.PARAMETER_SCALE_NORMAL, + "-1;0.5", + "priors", + (v2.C.NORMAL, [-1, 0.5]), + ), + ( + v1.C.PARAMETER_SCALE_NORMAL, + "-1;0.5", + "all", + (v2.C.NORMAL, [-1, 0.5]), + ), + ( + v1.C.PARAMETER_SCALE_UNIFORM, + "-3;1", + "priors", + (v2.C.UNIFORM, [-3, 1]), + ), + # no priors at all + (None, None, "all", (None, [])), + # no priors, but an empty prior column + (np.nan, np.nan, "all", (None, [])), + ], +) +def test_petab1to2_reparametrize_zero_lower_bound( + tmp_path, prior_type, prior_parameters, reparametrize, expected_prior +): + """Parameters with a lower bound of 0 and a proper prior on log scale + are reparametrized, with a lower bound of -inf.""" + yaml_file = _write_v1_problem_zero_lower_bound( + tmp_path, prior_type, prior_parameters + ) + + problem = petab1to2(yaml_file, reparametrize=reparametrize) + + assert [p.id for p in problem.parameters] == ["log10_k1"] + par = problem["log10_k1"] + assert (par.lb, par.ub) == (-np.inf, 2) + assert (par.prior_distribution, par.prior_parameters) == expected_prior + ia = problem.model.sbml_model.getInitialAssignmentBySymbol("k1") + assert ia is not None + + +@pytest.mark.filterwarnings( + "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning", + "ignore:Initialisation priors.*:UserWarning", +) +def test_petab1to2_reparametrize_zero_lower_bound_improper(tmp_path): + """Parameters with a lower bound of 0 and an improper prior on log scale + are skipped.""" + # the v1 default prior (uniform on parameter scale) is made explicit + # if other estimated parameters have priors + yaml_file = _write_v1_problem_zero_lower_bound( + tmp_path, np.nan, np.nan, sigma_prior=True + ) + + with pytest.warns(UserWarning, match=r"k1 were not reparametrized"): + problem = petab1to2(yaml_file, reparametrize="all") + + # as with reparametrize="none" + assert [p.id for p in problem.parameters] == ["k1", "sigma"] + par = problem["k1"] + assert (par.lb, par.ub) == (0.0, _K1[v1.C.UPPER_BOUND]) + assert par.prior_distribution == v2.C.UNIFORM + assert par.prior_parameters == [0.0, _K1[v1.C.UPPER_BOUND]] + sbml_model = problem.model.sbml_model + assert sbml_model.getInitialAssignmentBySymbol("k1") is None + assert sbml_model.getAssignmentRuleByVariable("k1") is None + + +def _write_v1_problem_with_overrides( + tmp_path, parameters: list[dict], overrides: list[str] +) -> str: + return _write_v1_problem( + tmp_path, + observables=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.OBSERVABLE_FORMULA: "observableParameter1_obs_B * B", + v1.C.NOISE_FORMULA: "0.1", + } + ], + measurements=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.SIMULATION_CONDITION_ID: "c0", + v1.C.TIME: float(i), + v1.C.MEASUREMENT: 0.3, + v1.C.OBSERVABLE_PARAMETERS: override, + } + for i, override in enumerate(overrides) + ], + parameters=[_K1, *parameters], + ) + + +_LOG10_PAR = { + v1.C.PARAMETER_SCALE: v1.C.LOG10, + v1.C.LOWER_BOUND: 1e-3, + v1.C.UPPER_BOUND: 1e3, + v1.C.NOMINAL_VALUE: 2, + v1.C.ESTIMATE: 1, +} + + +@pytest.mark.filterwarnings( + "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning" +) +def test_petab1to2_reparametrize_overrides(tmp_path): + """Parameters that share a placeholder with selected parameters are + reparametrized as well.""" + yaml_file = _write_v1_problem_with_overrides( + tmp_path, + parameters=[ + _LOG10_PAR + | { + v1.C.PARAMETER_ID: "s1", + v1.C.OBJECTIVE_PRIOR_TYPE: v1.C.PARAMETER_SCALE_NORMAL, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "0;1", + }, + _LOG10_PAR | {v1.C.PARAMETER_ID: "s2"}, + # non-estimated, fixed at 0, with a prior + _LOG10_PAR + | { + v1.C.PARAMETER_ID: "s3", + v1.C.ESTIMATE: 0, + v1.C.NOMINAL_VALUE: 0, + v1.C.OBJECTIVE_PRIOR_TYPE: v1.C.PARAMETER_SCALE_NORMAL, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "0;1", + }, + # non-estimated, with a prior that depends on missing bounds + _LOG10_PAR + | { + v1.C.PARAMETER_ID: "s4", + v1.C.ESTIMATE: 0, + v1.C.LOWER_BOUND: np.nan, + v1.C.UPPER_BOUND: np.nan, + v1.C.OBJECTIVE_PRIOR_TYPE: v1.C.PARAMETER_SCALE_UNIFORM, + }, + ], + overrides=["s1", "s2", "100", "s3", "s4", "0"], + ) + assert not v1.lint_problem(v1.Problem.from_yaml(yaml_file)) + # s2-s4 aren't selected, but share the placeholder with s1 + with pytest.warns(UserWarning, match="non-estimated parameters s4 were"): + problem = petab1to2(yaml_file, reparametrize="priors") + assert [p.id for p in problem.parameters] == [ + "k1", + "log10_s1", + "log10_s2", + "log10_s3", + "log10_s4", + # replaces the override `0` + "negative_infinity", + ] + assert problem["negative_infinity"].nominal_value == -np.inf + assert problem.measurements[-1].observable_parameters == [ + petab_symbol("negative_infinity") + ] + assert problem["log10_s1"].prior_distribution == v2.C.NORMAL + assert problem["log10_s1"].prior_parameters == [0, 1] + # explicit default prior, as other parameters have priors + assert problem["log10_s2"].prior_distribution == v2.C.UNIFORM + assert problem["log10_s2"].prior_parameters == pytest.approx([-3, 3]) + s3 = problem["log10_s3"] + assert not s3.estimate + assert s3.nominal_value == -np.inf + # the parameter-scale prior is kept + assert s3.prior_distribution == v2.C.NORMAL + assert s3.prior_parameters == [0, 1] + s4 = problem["log10_s4"] + assert s4.nominal_value == pytest.approx(np.log10(2)) + assert s4.prior_distribution is None + + +@pytest.mark.filterwarnings( + "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning" +) +def test_petab1to2_reparametrize_all_skips_unresolvable(tmp_path): + """With "all", parameters without parameter-scale priors are skipped if + their placeholder can't be log-transformed.""" + yaml_file = _write_v1_problem_with_overrides( + tmp_path, + parameters=[ + _LOG10_PAR | {v1.C.PARAMETER_ID: "s1"}, + _LOG10_PAR | {v1.C.PARAMETER_ID: "s2"}, + ], + overrides=["s1", "s2", "-1"], + ) + with pytest.warns(UserWarning, match="s1, s2 were not reparametrized"): + problem = petab1to2(yaml_file, reparametrize="all") + assert [p.id for p in problem.parameters] == ["k1", "s1", "s2"] + + yaml_file = _write_v1_problem_with_overrides( + tmp_path, + parameters=[ + _LOG10_PAR + | { + v1.C.PARAMETER_ID: "s1", + v1.C.OBJECTIVE_PRIOR_TYPE: v1.C.PARAMETER_SCALE_NORMAL, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "0;1", + }, + _LOG10_PAR | {v1.C.PARAMETER_ID: "s2"}, + ], + overrides=["s1", "s2", "-1"], + ) + # s1 has a prior, so it must be reparametrized + with ( + pytest.raises(NotImplementedError) as excinfo, + pytest.warns(UserWarning, match="s2 were not reparametrized"), + ): + petab1to2(yaml_file, reparametrize="all") + msg = str(excinfo.value) + assert "Use `reparametrize='none'`" in msg + assert "(overridden by `s1')" in msg + assert "not reparametrized: `s2'" in msg + assert "negative numbers: -1" in msg + + +@pytest.mark.filterwarnings( + "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning" +) +def test_petab1to2_reparametrize_output_parameter(tmp_path): + """Test reparametrization of parameters that are not model entities, + and that x-space priors are not reparametrized.""" + yaml_file = _write_v1_problem( + tmp_path, + observables=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.OBSERVABLE_FORMULA: "scaling * B", + v1.C.NOISE_FORMULA: "sigma", + } + ], + parameters=[ + _K1, + _LOG10_PAR | {v1.C.PARAMETER_ID: "scaling"}, + _LOG10_PAR + | { + v1.C.PARAMETER_ID: "sigma", + v1.C.OBJECTIVE_PRIOR_TYPE: v1.C.LOG_NORMAL, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "0;1", + }, + ], + ) + problem = petab1to2(yaml_file, reparametrize="all") + assert [p.id for p in problem.parameters] == [ + "k1", + "log10_scaling", + "sigma", + ] + assert problem["log10_scaling"].prior_distribution == v2.C.UNIFORM + + +@pytest.mark.parametrize( + "prior_type", + [ + v1.C.PARAMETER_SCALE_NORMAL, + v1.C.PARAMETER_SCALE_LAPLACE, + v1.C.PARAMETER_SCALE_UNIFORM, + ], +) +def test_petab1to2_reparametrize_warnings_and_errors(tmp_path, prior_type): + k1 = _K1 | { + v1.C.PARAMETER_SCALE: v1.C.LOG10, + v1.C.OBJECTIVE_PRIOR_TYPE: prior_type, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "-1;1", + } + yaml_file = _write_v1_problem( + tmp_path, + parameters=[k1], + ) + with pytest.warns(UserWarning, match=r"of k1 were converted.*'priors'"): + petab1to2(yaml_file) + + with pytest.raises(ValueError, match="Invalid value"): + petab1to2(yaml_file, reparametrize="some") + + def test_petab1to2_placeholder_order(tmp_path): """Test that placeholders are ordered by index, not as text, so that the positional overrides still match for 10+ placeholders.""" @@ -698,8 +1267,15 @@ def test_petab1to2_existing_experiment_id_column(tmp_path, clash): @pytest.mark.filterwarnings( "ignore:.*which is reserved in PEtab v2.*:UserWarning" ) +@pytest.mark.filterwarnings( + "ignore:Log-scaled parameters .* were not reparametrized.*:UserWarning" +) +@pytest.mark.filterwarnings( + "ignore:Parameter .* has no non-negative nominal value.*:UserWarning" +) +@pytest.mark.parametrize("reparametrize", ["none", "all"]) @parametrize_or_skip -def test_benchmark_collection(problem_id): +def test_benchmark_collection(problem_id, reparametrize): """Test that we can upgrade all benchmark collection models.""" logging.basicConfig(level=logging.DEBUG) @@ -709,8 +1285,10 @@ def test_benchmark_collection(problem_id): yaml_path = benchmark_models_petab.get_problem_yaml_path(problem_id) try: - problem = petab1to2(yaml_path) + problem = petab1to2(yaml_path, reparametrize=reparametrize) except NotImplementedError as e: + if "could not be reparametrized" in str(e): + raise pytest.skip(str(e)) assert isinstance(problem, Problem) assert len(problem.measurements) diff --git a/tests/v2/test_reparametrize.py b/tests/v2/test_reparametrize.py index 8be6c879..5d0f9307 100644 --- a/tests/v2/test_reparametrize.py +++ b/tests/v2/test_reparametrize.py @@ -383,24 +383,60 @@ def test_reparametrize_log_unresolvable_override(): "(overridden by `offset1') also receives:" in msg ) assert "not reparametrized: `offset2'\n" in msg - assert "non-positive" not in msg + assert "- negative numbers" not in msg - # non-positive numeric overrides, in several placeholders, reported at - # once + # negative numeric overrides, in several placeholders, reported at once problem = _problem() problem.measurements[0].noise_parameters = ["offset1"] - problem.measurements[1].observable_parameters = [0] + problem.measurements[1].observable_parameters = [-2] problem.measurements[1].noise_parameters = [-1] with pytest.raises(NotImplementedError) as excinfo: reparametrize_log(problem, ["offset1"]) msg = str(excinfo.value) assert "Observable placeholder `offset_A'" in msg assert "Noise placeholder `noiseParameter1_obs_A'" in msg - assert "non-positive numbers: 0\n" in msg - assert "non-positive numbers: -1\n" in msg + assert "negative numbers: -2\n" in msg + assert "negative numbers: -1\n" in msg assert "not reparametrized" not in msg +def test_reparametrize_log_zero_override(): + """Zero overrides are replaced by a fixed parameter with value -inf.""" + problem = _problem() + problem.measurements[0].noise_parameters = ["offset1"] + problem.measurements[1].observable_parameters = [0] + problem.measurements[1].noise_parameters = [0] + problem.parameter_tables[0].elements = [ + p for p in problem.parameters if p.id != "offset2" + ] + # the preferred ID is taken + sbml_par = problem.model.sbml_model.createParameter() + sbml_par.setId("negative_infinity") + sbml_par.setConstant(True) + assert not problem.validate().has_errors() + + problem, id_map = reparametrize_log(problem, ["offset1"]) + assert id_map == {"offset1": "log10_offset1"} + # a single parameter for all zero overrides + neg_inf = petab_symbol("negative_infinity_1") + assert [ + (m.observable_parameters, m.noise_parameters) + for m in problem.measurements + ] == [ + ([petab_symbol("log10_offset1")], [petab_symbol("log10_offset1")]), + ([neg_inf], [neg_inf]), + ] + par = problem["negative_infinity_1"] + assert not par.estimate + assert par.nominal_value == -math.inf + assert "negative_infinity_1" in [p.id for p in problem.parameters] + assert problem.observables[0].formula == ( + petab_symbol("scaling") * petab_symbol("A") + + 10 ** petab_symbol("offset_A") + ) + assert not problem.validate().has_errors() + + def test_reparametrize_log_noise_override(): """Observable placeholders in noise formulas and noise overrides.""" problem = _problem()