From bdad27e8782e967dd745465777e0a4c17f7604ec Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 13 Oct 2026 22:12:37 +0200 Subject: [PATCH 1/6] Add `reparametrize_log` and `petab1to2(reparametrize=...)` PEtab v2 has no `parameterScale`. Converting a v1 `parameterScaleNormal` prior on a log-scaled parameter to `log-normal` keeps the prior distribution, but adds a Jacobian term to the MAP objective and moves the optimum. No v2 prior over the unscaled parameter can preserve both. Add `petab.v2.reparametrize_log`, which estimates `log10_` (or `log_`) instead of `` and defines ` = 10^log10_`: * parameter table: log-transformed bounds and nominal value; priors are taken from the `priors` argument or derived from `log-normal`, `log-laplace` and `log-uniform` priors * SBML model: new parameter and assignment rule for model parameters * observable/noise formulas and condition target values: substitution * measurement overrides: replaced by the new IDs, numbers are log-transformed, and the placeholder is exponentiated in the formula `petab1to2` and `petab_files_1to2` get a `reparametrize` argument: `"none"` (default, current behaviour), `"priors"` (log-scaled parameters with `parameterScale*` priors) or `"all"` (all estimated log-scaled parameters without priors on linear scale). Log-scaled parameters that share a placeholder with selected parameters are added. With `"all"`, parameters without parameter-scale priors are skipped if their placeholder can't be log-transformed. With `"none"`, a warning points to the new option. References https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/issues/333 Generated with AI assistance. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 6 + doc/modules.rst | 1 + petab/v2/__init__.py | 1 + petab/v2/petab1to2.py | 229 +++++++++++++++++- petab/v2/reparametrize.py | 428 +++++++++++++++++++++++++++++++++ tests/v2/test_conversion.py | 326 +++++++++++++++++++++++++ tests/v2/test_reparametrize.py | 293 ++++++++++++++++++++++ 7 files changed, 1280 insertions(+), 4 deletions(-) create mode 100644 petab/v2/reparametrize.py create mode 100644 tests/v2/test_reparametrize.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 87ad83c9..75e4986c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -63,6 +63,12 @@ `Problem.to_files(columns=...)`, or globally via `petab.v2.set_options` and `petab.v2.option_context` (by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/543) +* `petab.v2.reparametrize_log` estimates parameters on log scale by replacing + them with new log-transformed parameters. `petab1to2` and + `petab_files_1to2` can use it via `reparametrize="priors"` or + `reparametrize="all"` to preserve the objective function of log-scaled + parameters with `parameterScale*` priors + (by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/539) ## 0.9 series diff --git a/doc/modules.rst b/doc/modules.rst index 66bdcc37..1658aa1d 100644 --- a/doc/modules.rst +++ b/doc/modules.rst @@ -40,3 +40,4 @@ API Reference petab.v2.models petab.v2.options petab.v2.petab1to2 + petab.v2.reparametrize diff --git a/petab/v2/__init__.py b/petab/v2/__init__.py index f6c36ca3..93c07518 100644 --- a/petab/v2/__init__.py +++ b/petab/v2/__init__.py @@ -38,3 +38,4 @@ from .lint import lint_problem # noqa: F401 from .models import MODEL_TYPE_PYSB, MODEL_TYPE_SBML, Model # noqa: F401 from .options import * +from .reparametrize import reparametrize_log # noqa: F401 diff --git a/petab/v2/petab1to2.py b/petab/v2/petab1to2.py index e3cb19c3..ebce0256 100644 --- a/petab/v2/petab1to2.py +++ b/petab/v2/petab1to2.py @@ -8,10 +8,12 @@ from contextlib import suppress from pathlib import Path from tempfile import TemporaryDirectory +from typing import Literal, get_args from urllib.parse import urlparse from uuid import uuid4 import pandas as pd +import sympy as sp from pandas.io.common import get_handle, is_url from .. import v1, v2 @@ -20,15 +22,20 @@ 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, 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. @@ -56,16 +63,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 reparametrize. + PEtab v2 has no ``parameterScale``. Parameters are estimated on linear + scale, and priors are densities over the unscaled parameter. + Therefore, a v1 ``parameterScale*`` prior on a log-scaled parameter + cannot be converted such that both the prior distribution and the + maximum a posteriori estimate are preserved. + Reparametrization estimates ``log10_`` (``log_``) instead of + ````, which preserves both + (see :func:`petab.v2.reparametrize_log`): + + * ``"none"`` (default): don't reparametrize. Priors are converted to + the corresponding densities over the unscaled parameters + (e.g., ``parameterScaleNormal`` to ``log-normal``). + This preserves the prior distribution, but not the objective + function. + * ``"priors"``: reparametrize all estimated log-scaled parameters with + a ``parameterScale*`` prior. + * ``"all"``: reparametrize all estimated log-scaled parameters, + except for those with priors on linear scale (``normal``, ...). + This also preserves the default prior + (uniform on parameter scale). + + Further log-scaled parameters are reparametrized if they share + a placeholder with reparametrized parameters 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. """ 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)) @@ -74,6 +116,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. @@ -85,6 +128,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 @@ -102,6 +148,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: @@ -121,10 +173,24 @@ def petab_files_1to2( new_yaml_config = _update_yaml(yaml_config) new_yaml_config = v2.ProblemConfig(**new_yaml_config) + if reparametrize != "none": + # The model will be modified. Don't overwrite the original one. + for file in ( + model.location for model in new_yaml_config.model_files.values() + ): + with suppress(FileNotFoundError, TypeError): + if Path(get_dest_path(file)).samefile(get_src_path(file)): + raise ValueError( + "Reparametrization requires the output directory to " + "differ from the input directory." + ) + # 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]) ) @@ -266,6 +332,13 @@ 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": + problem = v2.Problem.from_yaml(new_yaml_file) + if _reparametrize_log_scaled( + problem, petab_problem.parameter_df, reparametrize + ): + problem.to_files(output_dir) + if not validate: return @@ -498,11 +571,17 @@ def extract_placeholders(row: pd.Series, type_: str) -> str: 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() @@ -606,10 +685,33 @@ def update_prior(row) -> tuple[str, str | None]: stacklevel=9, ) if not (df[v1.C.PARAMETER_SCALE] == v1.C.LIN).all(): + if reparametrize == "none": + hint = ( + " Log-scaled parameters are estimated on linear scale. " + "Use `reparametrize='all'` to estimate them on log scale." + ) + else: + hint = "" warnings.warn( - "Parameter scales are not supported in PEtab v2.", + f"Parameter scales are not supported in PEtab v2.{hint}", stacklevel=9, ) + if reparametrize == "none" and v1.C.OBJECTIVE_PRIOR_TYPE in parameter_df: + has_scale_prior = parameter_df[v1.C.PARAMETER_SCALE].isin( + [v1.C.LOG, v1.C.LOG10] + ) & parameter_df[v1.C.OBJECTIVE_PRIOR_TYPE].isin( + [v1.C.PARAMETER_SCALE_NORMAL, v1.C.PARAMETER_SCALE_LAPLACE] + ) + if has_scale_prior.any(): + warnings.warn( + "Parameter-scale priors of log-scaled parameters " + f"{', '.join(parameter_df.index[has_scale_prior])} were " + "converted to the corresponding priors of the unscaled " + "parameters. This preserves the prior distribution, but " + "changes the maximum a posteriori estimate. " + "Use `reparametrize='priors'` to preserve both.", + stacklevel=9, + ) df.drop( columns=[ v1.C.INITIALIZATION_PRIOR_TYPE, @@ -638,3 +740,122 @@ 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, +) -> bool: + """Reparametrize log-scaled parameters of a converted problem. + + See :func:`petab1to2`. + + :param problem: The converted problem, modified in place. + :param v1_parameter_df: The original PEtab v1 parameter table. + :param reparametrize: The parameters to reparametrize. + :return: Whether any parameter was reparametrized. + """ + df = v1_parameter_df + has_priors = v1.C.OBJECTIVE_PRIOR_TYPE in 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()} + for par_id, row in df.iterrows(): + pscale = row.get(v1.C.PARAMETER_SCALE, v1.C.LIN) + if pscale not in (v1.C.LOG, v1.C.LOG10): + continue + base = 10 if pscale == v1.C.LOG10 else "e" + 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 + + prior = None + if not int(row[v1.C.ESTIMATE]): + pass + elif prior_type in (None, v1.C.PARAMETER_SCALE_UNIFORM): + lb = scale(row[v1.C.LOWER_BOUND], pscale) + ub = scale(row[v1.C.UPPER_BOUND], pscale) + prior_pars = row.get(v1.C.OBJECTIVE_PRIOR_PARAMETERS) + if prior_type is not None and not pd.isna(prior_pars): + low, high = map( + float, str(prior_pars).split(v1.C.PARAMETER_SEPARATOR) + ) + prior = (v2.C.UNIFORM, [max(low, lb), min(high, ub)]) + elif has_priors: + # The v2 default prior is uniform over the parameter, but + # other parameters have explicit priors. + # Be explicit for consistency with `v1v2_parameter_df`. + prior = (v2.C.UNIFORM, [lb, ub]) + else: + prior = ( + prior_type.removeprefix("parameterScale").lower(), + list( + map( + float, + str(row[v1.C.OBJECTIVE_PRIOR_PARAMETERS]).split( + v1.C.PARAMETER_SEPARATOR + ), + ) + ), + ) + + candidates[base][par_id] = prior + if int(row[v1.C.ESTIMATE]) and prior_type is not None: + required[base].add(par_id) + if int(row[v1.C.ESTIMATE]) and ( + reparametrize == "all" or prior_type is not None + ): + selected[base].add(par_id) + + # Placeholders can only be transformed if all their overrides are + # reparametrized or positive 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() + for base, base_selected in selected.items(): + excluded = set() + 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 not 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, + ) + + for base, base_selected in selected.items(): + 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}, + ) + + return any(selected.values()) diff --git a/petab/v2/reparametrize.py b/petab/v2/reparametrize.py new file mode 100644 index 00000000..8dc62560 --- /dev/null +++ b/petab/v2/reparametrize.py @@ -0,0 +1,428 @@ +"""Reparametrization of PEtab problems.""" + +from __future__ import annotations + +import math +import warnings +from collections.abc import Iterable +from itertools import chain +from typing import Literal + +import libsbml +import sympy as sp + +from . import C +from .core import Parameter, PriorDistribution, Problem +from .lint import get_placeholders +from .models._sbml_utils import add_sbml_parameter, check +from .models.sbml_model import SbmlModel + +__all__ = ["reparametrize_log"] + +#: x-space priors and the corresponding priors of log_b(x), given as +#: (prior type, function(prior parameters, ln(b)) -> prior parameters) +_LOG_PRIOR_TO_THETA = { + C.LOG_NORMAL: (C.NORMAL, lambda pars, ln_b: [p / ln_b for p in pars]), + C.LOG_LAPLACE: (C.LAPLACE, lambda pars, ln_b: [p / ln_b for p in pars]), + C.LOG_UNIFORM: ( + C.UNIFORM, + lambda pars, ln_b: [math.log(p) / ln_b for p in pars], + ), +} + + +def reparametrize_log( + problem: Problem, + parameter_ids: Iterable[str], + *, + base: Literal[10, "e"] = 10, + priors: dict[str, tuple[str, list[float]] | None] | None = None, +) -> dict[str, str]: + """Estimate parameters on log scale. + + For each parameter ``p``, replace ``p`` in the parameter table by a new + parameter ``q = log_base(p)`` and define ``p = base^q``. + This is the PEtab v2 way of estimating parameters on log scale + (and of specifying priors for the log-transformed parameter). + + The problem is modified in place: + + * Parameter table: ``p`` is replaced by ``q`` with log-transformed bounds + and nominal value. + Parameters without a positive nominal value get no nominal value for + ``q``. + * Model: If ``p`` is a model parameter, ``q`` is added to the model, + and ``p`` becomes the target of the assignment rule ``p = base^q``. + Only SBML models are supported. + * Observable and noise formulas, and condition table target values: + ``p`` is substituted by ``base^q``. + * 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. + + The new ID is ``log10_

`` for base 10 and ``log_

`` for base e, + with a numeric suffix if necessary to make it unique. + + :param problem: + The PEtab problem to modify. + :param parameter_ids: + IDs of the parameters to reparametrize. + :param base: + The base of the logarithm. + :param priors: + Prior distributions of the new parameters as + ``{p: (prior_distribution, prior_parameters)}``, or ``{p: None}`` + for no prior. + For parameters not included here, the prior is derived from the prior + of ``p``, where possible: + ``log-normal(m, s)`` -> ``normal(m / ln(b), s / ln(b))``, + ``log-laplace(m, s)`` -> ``laplace(m / ln(b), s / ln(b))``, + ``log-uniform(a, b)`` -> ``uniform(log_b(a), log_b(b))``, + no prior -> no prior (i.e., uniform on log scale). + For non-estimated parameters, no prior is derived. + :return: + A dictionary mapping old to new parameter IDs. + :raises ValueError: + If a parameter cannot be reparametrized, e.g., because its lower bound + is not positive, or if its prior cannot be expressed in terms of the + new parameter. + :raises NotImplementedError: + If the problem has features that are not supported yet, + e.g., multiple or non-SBML models, mapping table entries for the + parameters, or placeholders with overrides that cannot be + log-transformed. + """ + parameter_ids = list(dict.fromkeys(parameter_ids)) + priors = priors or {} + if base == 10: + ln_b = math.log(10) + b_expr = sp.Integer(10) + prefix = "log10_" + elif base == "e": + ln_b = 1.0 + b_expr = sp.E + prefix = "log_" + else: + raise ValueError(f"Unsupported base: {base}") + + if unknown := set(priors) - set(parameter_ids): + raise ValueError(f"Priors given for unknown parameters: {unknown}") + + if not parameter_ids: + return {} + + if problem.extensions.sciml is not None: + raise NotImplementedError( + "Reparametrization of PEtab SciML problems is not supported." + ) + + if len(problem.models) > 1: + raise NotImplementedError( + "Reparametrization of multi-model problems is not supported." + ) + model = problem.model + + def log_b(x: float) -> float: + return math.log(x) / ln_b + + new_parameters = {} + model_entities = [] + for p in parameter_ids: + par = _get_parameter(problem, p) + if any(p in (m.petab_id, m.model_id) for m in problem.mappings): + raise NotImplementedError( + f"Parameter `{p}' occurs in the mapping table." + ) + if model is not None and model.has_entity_with_id(p): + if not isinstance(model, SbmlModel): + raise NotImplementedError( + "Reparametrization of model parameters is only " + "supported for SBML models." + ) + if model.sbml_model.getParameter(p) is None: + raise NotImplementedError( + f"Model entity `{p}' is not an SBML parameter." + ) + model_entities.append(p) + if any(bound is not None and bound <= 0 for bound in (par.lb, par.ub)): + # (bounds are required for estimated parameters) + raise ValueError( + f"Bounds of parameter `{p}' must be positive, " + f"but got [{par.lb}, {par.ub}]." + ) + + if ( + par.nominal_value is None + or par.nominal_value == "array" + or par.nominal_value <= 0 + ): + if not par.estimate: + raise ValueError( + f"Non-estimated parameter `{p}' must have a positive " + f"nominal value, but got {par.nominal_value}." + ) + warnings.warn( + f"Parameter `{p}' has no positive nominal value " + f"({par.nominal_value}). The nominal value of its " + "reparametrization is left empty.", + stacklevel=2, + ) + nominal_value = None + else: + nominal_value = log_b(par.nominal_value) + + if p in priors and priors[p] is not None: + prior_distribution, prior_parameters = priors[p] + prior_distribution = PriorDistribution(prior_distribution) + elif p in priors or par.prior_distribution is None or not par.estimate: + # priors of non-estimated parameters are irrelevant + prior_distribution, prior_parameters = None, [] + elif par.prior_distribution in _LOG_PRIOR_TO_THETA: + prior_distribution, f = _LOG_PRIOR_TO_THETA[par.prior_distribution] + prior_parameters = f(par.prior_parameters, ln_b) + else: + raise ValueError( + f"Prior `{par.prior_distribution}' of parameter `{p}' " + "cannot be expressed in terms of the log-transformed " + "parameter." + ) + + new_parameters[p] = Parameter( + id=p, + name=par.name and f"{prefix.removesuffix('_')}({par.name})", + lb=None if par.lb is None else log_b(par.lb), + ub=None if par.ub is None else log_b(par.ub), + nominal_value=nominal_value, + estimate=par.estimate, + prior_distribution=prior_distribution, + prior_parameters=prior_parameters, + **(par.model_extra or {}), + ) + + # assign new IDs + used_ids = _get_used_ids(problem) + id_map = {} + for p in parameter_ids: + q = f"{prefix}{p}" + i = 0 + while q in used_ids or ( + model is not None and model.has_entity_with_id(q) + ): + i += 1 + q = f"{prefix}{p}_{i}" + used_ids.add(q) + id_map[p] = q + new_parameters[p].id = q + + _reparametrize_overrides(problem, id_map, b_expr, log_b) + + # substitute p -> b^q in all expressions + def subs(expr: sp.Basic | None) -> sp.Basic | None: + if expr is None: + return expr + return expr.subs( + { + s: b_expr ** sp.Symbol(id_map[s.name], real=True) + for s in expr.free_symbols + if isinstance(s, sp.Symbol) and s.name in id_map + } + ) + + for observable in problem.observables: + observable.formula = subs(observable.formula) + observable.noise_formula = subs(observable.noise_formula) + for condition in problem.conditions: + for change in condition.changes: + change.target_value = subs(change.target_value) + + # update parameter tables + for table in problem.parameter_tables: + table.elements = [ + new_parameters.get(par.id, par) for par in table.elements + ] + + # update the model + for p in model_entities: + _add_sbml_log_rule( + model.sbml_model, + p, + id_map[p], + b_expr, + new_parameters[p].nominal_value, + ) + + return id_map + + +def _get_parameter(problem: Problem, parameter_id: str) -> Parameter: + """Get the parameter table entry for the given ID.""" + for par in problem.parameters: + if par.id == parameter_id: + return par + raise ValueError( + f"Parameter `{parameter_id}' is not in the parameter table." + ) + + +def _get_used_ids(problem: Problem) -> set[str]: + """Get all IDs and symbols used in the PEtab tables.""" + ids = set(get_placeholders(problem)) + ids |= {p.id for p in problem.parameters} + ids |= {o.id for o in problem.observables} + ids |= {c.id for c in problem.conditions} + ids |= {e.id for e in problem.experiments} + ids |= {m.petab_id for m in problem.mappings} + ids |= {m.model_id for m in problem.mappings if m.model_id} + ids |= { + change.target_id + for cond in problem.conditions + for change in cond.changes + } + expressions = chain( + (o.formula for o in problem.observables), + (o.noise_formula for o in problem.observables), + ( + change.target_value + for cond in problem.conditions + for change in cond.changes + ), + chain.from_iterable( + chain(m.observable_parameters, m.noise_parameters) + for m in problem.measurements + ), + ) + ids |= { + str(s) + for expr in expressions + if expr is not None + for s in expr.free_symbols + } + return ids + + +def _get_override_slots( + problem: Problem, +) -> dict[tuple[str, str, int], set[sp.Basic]]: + """Get the values of all placeholder overrides. + + :return: + Dictionary mapping ``(observable ID, placeholder type, index)`` + to the set of override values. The placeholder type is + ``"observable"`` or ``"noise"``. + """ + slots = {} + for m in problem.measurements: + for type_, overrides in ( + ("observable", m.observable_parameters), + ("noise", m.noise_parameters), + ): + for i, value in enumerate(overrides): + slots.setdefault((m.observable_id, type_, i), set()).add(value) + return slots + + +def _reparametrize_overrides( + problem: Problem, + id_map: dict[str, str], + b_expr: sp.Expr, + log_b, +) -> None: + """Log-transform placeholders overridden by reparametrized parameters.""" + slots_to_transform = [] + for slot, values in _get_override_slots(problem).items(): + if not any( + isinstance(v, sp.Symbol) and v.name in id_map for v in values + ): + continue + invalid = { + v + for v in values + if not ( + (isinstance(v, sp.Symbol) and v.name in id_map) + or (v.is_number and v > 0) + ) + } + if invalid: + obs_id, type_, i = slot + raise NotImplementedError( + f"Cannot reparametrize {type_} placeholder {i + 1} of " + f"observable `{obs_id}': all its overrides must be " + "reparametrized parameters or positive numbers, but got " + f"{sorted(map(str, invalid))}." + ) + slots_to_transform.append(slot) + + if not slots_to_transform: + return + + # placeholder -> b^placeholder in formulas + observables = {o.id: o for o in problem.observables} + for obs_id, type_, i in slots_to_transform: + observable = observables[obs_id] + if type_ == "observable": + placeholder = observable.observable_placeholders[i] + observable.formula = observable.formula.subs( + placeholder, b_expr**placeholder + ) + else: + placeholder = observable.noise_placeholders[i] + # observable placeholders may also occur in the noise formula + observable.noise_formula = observable.noise_formula.subs( + placeholder, b_expr**placeholder + ) + + # update override values + def transform(value: sp.Basic) -> sp.Basic: + if isinstance(value, sp.Symbol): + return sp.Symbol(id_map[value.name], real=True) + return sp.Float(log_b(float(value))) + + slots_to_transform = set(slots_to_transform) + for m in problem.measurements: + for type_, attr in ( + ("observable", "observable_parameters"), + ("noise", "noise_parameters"), + ): + overrides = getattr(m, attr) + if any( + (m.observable_id, type_, i) in slots_to_transform + for i in range(len(overrides)) + ): + setattr( + m, + attr, + [ + transform(v) + if (m.observable_id, type_, i) in slots_to_transform + else v + for i, v in enumerate(overrides) + ], + ) + + +def _add_sbml_log_rule( + model: libsbml.Model, + parameter_id: str, + new_id: str, + b_expr: sp.Expr, + value: float | None, +) -> None: + """Add parameter ``new_id`` and the rule ``parameter_id = b^new_id``.""" + add_sbml_parameter(model, new_id, value=value, constant=True) + + par = model.getParameter(parameter_id) + check(par.setConstant(False)) + # an assignment rule and an initial assignment must not have the same + # target + if model.getInitialAssignmentBySymbol(parameter_id) is not None: + model.removeInitialAssignment(parameter_id) + + rule = model.createAssignmentRule() + check(rule.setVariable(parameter_id)) + formula = f"exp({new_id})" if b_expr == sp.E else f"10^({new_id})" + ast = libsbml.parseL3Formula(formula) + if ast is None: + raise RuntimeError(f"Failed to parse `{formula}'.") + check(rule.setMath(ast)) diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 99625b94..7afd9562 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -4,9 +4,11 @@ import numpy as np import pandas as pd import pytest +import sympy as sp from petab import v1, v2 from petab.v1.calculate import calculate_llh as calculate_llh_v1 +from petab.v1.math import sympify_petab from petab.v1.models.sbml_model import SbmlModel from petab.v1.priors import Prior from petab.v2 import Problem @@ -200,6 +202,9 @@ def simulate(measurement_df): @pytest.mark.filterwarnings( "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning" ) +@pytest.mark.filterwarnings( + "ignore:.*Parameter-scale priors of log-scaled parameters.*:UserWarning" +) @pytest.mark.parametrize( "parameter_scale, prior_type, prior_parameters", [ @@ -263,6 +268,324 @@ def test_petab1to2_priors( assert v2_prior.pdf(x) == pytest.approx(v1_prior.pdf(x), rel=1e-12) +@pytest.mark.filterwarnings( + "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning" +) +@pytest.mark.parametrize("reparametrize", ["priors", "all"]) +@pytest.mark.parametrize( + "parameter_scale, prior_type, prior_parameters", + [ + (v1.C.LOG10, v1.C.PARAMETER_SCALE_NORMAL, "-0.5;0.7"), + (v1.C.LOG10, v1.C.PARAMETER_SCALE_LAPLACE, "-0.5;0.7"), + (v1.C.LOG, v1.C.PARAMETER_SCALE_NORMAL, "-0.5;0.7"), + (v1.C.LOG, v1.C.PARAMETER_SCALE_LAPLACE, "-0.5;0.7"), + (v1.C.LOG10, v1.C.PARAMETER_SCALE_UNIFORM, "-0.5;0.7"), + (v1.C.LOG10, v1.C.PARAMETER_SCALE_UNIFORM, "-400;400"), + (v1.C.LOG10, v1.C.PARAMETER_SCALE_UNIFORM, np.nan), + (v1.C.LOG, v1.C.PARAMETER_SCALE_UNIFORM, "-1000;1000"), + (v1.C.LOG10, np.nan, np.nan), + ], +) +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, + observables=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.OBSERVABLE_FORMULA: "B", + v1.C.NOISE_FORMULA: "0.1", + } + ], + measurements=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.SIMULATION_CONDITION_ID: "c0", + v1.C.TIME: 1.0, + v1.C.MEASUREMENT: 0.3, + } + ], + 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.getAssignmentRuleByVariable("k1") is not None + assert sbml_model.getParameter(new_id).getConstant() + + # 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) + + +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): + """Test reparametrization of placeholder overrides, including + parameters that are added because they share a placeholder.""" + 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"}, + _LOG10_PAR | {v1.C.PARAMETER_ID: "s3"}, + ], + overrides=["s1", "s2", "100", "s3"], + ) + # s2 and s3 have no prior, but share the placeholder with s1 + problem = petab1to2(yaml_file, reparametrize="priors") + assert [p.id for p in problem.parameters] == [ + "k1", + "log10_s1", + "log10_s2", + "log10_s3", + ] + 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]) + + observable = problem["obs_B"] + placeholder = observable.observable_placeholders[0] + assert observable.formula == sympify_petab(f"10^{placeholder} * B") + assert [m.observable_parameters[0] for m in problem.measurements] == [ + sp.Symbol("log10_s1", real=True), + sp.Symbol("log10_s2", real=True), + pytest.approx(2.0), + sp.Symbol("log10_s3", real=True), + ] + + +@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", "0"], + ) + 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", "0"], + ) + # s1 has a prior, so it must be reparametrized + with ( + pytest.raises(NotImplementedError, match=r"\['0(\.0)?', 's2'\]"), + pytest.warns(UserWarning, match="s2 were not reparametrized"), + ): + petab1to2(yaml_file, reparametrize="all") + + +@pytest.mark.filterwarnings( + "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning" +) +def test_petab1to2_reparametrize_overrides_unresolvable(tmp_path): + 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", v1.C.PARAMETER_SCALE: v1.C.LIN}, + ], + overrides=["s1", "s2"], + ) + with pytest.raises(NotImplementedError, match=r"\['s2'\]"): + petab1to2(yaml_file, reparametrize="priors") + + +@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", + } + ], + measurements=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.SIMULATION_CONDITION_ID: "c0", + v1.C.TIME: 1.0, + v1.C.MEASUREMENT: 0.3, + } + ], + 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 + assert problem["obs_B"].formula == sympify_petab("10^log10_scaling * B") + + +def test_petab1to2_reparametrize_warnings_and_errors(tmp_path): + k1 = _K1 | { + v1.C.PARAMETER_SCALE: v1.C.LOG10, + 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: "B", + v1.C.NOISE_FORMULA: "0.1", + } + ], + measurements=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.SIMULATION_CONDITION_ID: "c0", + v1.C.TIME: 1.0, + v1.C.MEASUREMENT: 0.3, + } + ], + parameters=[k1], + ) + with ( + pytest.warns(UserWarning, match="reparametrize='all'"), + pytest.warns(UserWarning, match=r"k1 were converted.*'priors'"), + ): + petab1to2(yaml_file) + + with pytest.raises(ValueError, match="Invalid value"): + petab1to2(yaml_file, reparametrize="some") + + with pytest.raises(ValueError, match="output directory"): + petab1to2(yaml_file, tmp_path, reparametrize="all") + + 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.""" @@ -638,6 +961,9 @@ def test_petab1to2_validate(tmp_path): @pytest.mark.filterwarnings( "ignore:.*does not support log10-transformed observables.*:UserWarning" ) +@pytest.mark.filterwarnings( + "ignore:.*Parameter-scale priors of log-scaled parameters.*:UserWarning" +) @parametrize_or_skip def test_benchmark_collection(problem_id): """Test that we can upgrade all benchmark collection models.""" diff --git a/tests/v2/test_reparametrize.py b/tests/v2/test_reparametrize.py new file mode 100644 index 00000000..9a8994b6 --- /dev/null +++ b/tests/v2/test_reparametrize.py @@ -0,0 +1,293 @@ +import math + +import libsbml +import numpy as np +import pytest +import sympy as sp + +from petab.v1.math import sympify_petab +from petab.v2 import C, Problem, reparametrize_log +from petab.v2.models.sbml_model import SbmlModel + + +def _problem() -> Problem: + """A problem with a model parameter, an output parameter, and a + placeholder override.""" + problem = Problem() + problem.model = SbmlModel.from_antimony( + """ + k1 = 1 + k2 = 1 + A' = -k1 * A + A = 1 + """ + ) + problem.add_observable( + "obs_A", + "scaling * A + offset_A", + noise_formula="noiseParameter1_obs_A", + observable_placeholders=["offset_A"], + noise_placeholders=["noiseParameter1_obs_A"], + ) + problem.add_condition("c1", k2="k2_c1") + problem.add_experiment("e1", 0, "c1") + for offset in ("offset1", "offset2"): + problem.add_measurement( + "obs_A", + experiment_id="e1", + time=1, + measurement=1, + observable_parameters=[offset], + noise_parameters=[0.5], + ) + problem.add_parameter( + "k1", + nominal_value=2, + lb=1e-3, + ub=1e3, + prior_dist=C.LOG_NORMAL, + prior_pars=[1, 2], + parameterName="k one", + ) + problem.add_parameter( + "k2_c1", + nominal_value=3, + lb=1e-2, + ub=1e2, + prior_dist=C.LOG_UNIFORM, + prior_pars=[1e-1, 1e1], + ) + problem.add_parameter("scaling", nominal_value=4, lb=1e-3, ub=1e3) + problem.add_parameter("offset1", nominal_value=5, lb=1e-3, ub=1e3) + problem.add_parameter("offset2", estimate=False, nominal_value=6) + assert not problem.validate().has_errors() + return problem + + +def _sym(x: str) -> sp.Symbol: + return sp.Symbol(x, real=True) + + +@pytest.mark.parametrize("base", [10, "e"]) +def test_reparametrize_log_parameter_table(base): + problem = _problem() + ln_b = math.log(10) if base == 10 else 1.0 + prefix = "log10_" if base == 10 else "log_" + + id_map = reparametrize_log( + problem, ["k1", "k2_c1", "scaling", "offset1", "offset2"], base=base + ) + assert id_map == { + "k1": f"{prefix}k1", + "k2_c1": f"{prefix}k2_c1", + "scaling": f"{prefix}scaling", + "offset1": f"{prefix}offset1", + "offset2": f"{prefix}offset2", + } + assert [p.id for p in problem.parameters] == list(id_map.values()) + offset2 = problem[f"{prefix}offset2"] + assert offset2.estimate is False + assert offset2.nominal_value == pytest.approx(math.log(6) / ln_b) + + k1 = problem[f"{prefix}k1"] + assert k1.lb == pytest.approx(math.log(1e-3) / ln_b) + assert k1.ub == pytest.approx(math.log(1e3) / ln_b) + assert k1.nominal_value == pytest.approx(math.log(2) / ln_b) + assert k1.estimate is True + assert k1.prior_distribution == C.NORMAL + assert k1.prior_parameters == pytest.approx([1 / ln_b, 2 / ln_b]) + assert k1.name == (f"{prefix.removesuffix('_')}(k one)") + + k2 = problem[f"{prefix}k2_c1"] + assert k2.prior_distribution == C.UNIFORM + assert k2.prior_parameters == pytest.approx( + [math.log(1e-1) / ln_b, math.log(1e1) / ln_b] + ) + + scaling = problem[f"{prefix}scaling"] + assert scaling.prior_distribution is None + assert scaling.name is None + + assert not problem.validate().has_errors() + + +def test_reparametrize_log_expressions_and_model(): + problem = _problem() + reparametrize_log( + problem, ["k1", "k2_c1", "scaling", "offset1", "offset2"] + ) + + # observable formula + obs = problem["obs_A"] + assert obs.formula == ( + 10 ** _sym("log10_scaling") * _sym("A") + 10 ** _sym("offset_A") + ) + assert obs.noise_formula == sympify_petab("noiseParameter1_obs_A") + + # condition table + assert problem["c1"].changes[0].target_value == 10 ** _sym("log10_k2_c1") + + # overrides + assert [m.observable_parameters for m in problem.measurements] == [ + [_sym("log10_offset1")], + [_sym("log10_offset2")], + ] + assert [m.noise_parameters for m in problem.measurements] == [ + [0.5], + [0.5], + ] + + # model + sbml_model = problem.model.sbml_model + assert sbml_model.getParameter("k1").getConstant() is False + q = sbml_model.getParameter("log10_k1") + assert q.getConstant() is True + assert q.getValue() == pytest.approx(math.log10(2)) + rule = sbml_model.getAssignmentRuleByVariable("k1") + assert rule is not None + assert sympify_petab(rule.getFormula()) == sympify_petab("10^log10_k1") + assert problem.model.is_valid() + # k2 is not reparametrized, it's a condition target + assert sbml_model.getAssignmentRuleByVariable("k2") is None + + assert not problem.validate().has_errors() + + +def test_reparametrize_log_initial_assignment_removed(): + problem = Problem() + problem.model = SbmlModel.from_antimony("k1 = 2; A' = -k1 * A; A = 1") + sbml_model = problem.model.sbml_model + ia = sbml_model.createInitialAssignment() + ia.setSymbol("k1") + ia.setMath(libsbml.parseL3Formula("3")) + problem.add_parameter("k1", nominal_value=2, lb=1, ub=10) + reparametrize_log(problem, ["k1"]) + assert sbml_model.getInitialAssignmentBySymbol("k1") is None + assert sbml_model.getAssignmentRuleByVariable("k1") is not None + + +def test_reparametrize_log_id_collision(): + problem = _problem() + problem.add_parameter("log10_scaling", nominal_value=1, lb=1, ub=2) + problem.observables[0].formula += _sym("log10_scaling") + problem.model.sbml_model.createParameter().setId("log10_scaling_1") + id_map = reparametrize_log(problem, ["scaling"]) + assert id_map == {"scaling": "log10_scaling_2"} + + +def test_reparametrize_log_nominal(): + problem = _problem() + problem["scaling"].nominal_value = None + problem["offset1"].nominal_value = 0 + with pytest.warns(UserWarning, match="no positive nominal value"): + reparametrize_log(problem, ["scaling", "offset1", "offset2"]) + assert problem["log10_scaling"].nominal_value is None + assert problem["log10_offset1"].nominal_value is None + + +def test_reparametrize_log_explicit_priors(): + problem = _problem() + reparametrize_log( + problem, + ["k1", "scaling"], + priors={"k1": (C.LAPLACE, [0, 1]), "scaling": (C.NORMAL, [1, 2])}, + ) + assert problem["log10_k1"].prior_distribution == C.LAPLACE + assert problem["log10_k1"].prior_parameters == [0, 1] + assert problem["log10_scaling"].prior_distribution == C.NORMAL + assert problem["log10_scaling"].prior_parameters == [1, 2] + + with pytest.raises(ValueError, match="unknown parameters"): + reparametrize_log(_problem(), ["k1"], priors={"k2": (C.NORMAL, [])}) + + +@pytest.mark.parametrize("prior", [C.UNIFORM, C.NORMAL, C.GAMMA]) +def test_reparametrize_log_inexpressible_prior(prior): + problem = _problem() + problem["scaling"].prior_distribution = prior + problem["scaling"].prior_parameters = [1, 2] + with pytest.raises(ValueError, match="cannot be expressed"): + reparametrize_log(problem, ["scaling"]) + + +def test_reparametrize_log_invalid(): + problem = _problem() + problem["scaling"].lb = 0 + with pytest.raises(ValueError, match="must be positive"): + reparametrize_log(problem, ["scaling"]) + + with pytest.raises(ValueError, match="not in the parameter table"): + reparametrize_log(_problem(), ["nonexistent"]) + + problem = _problem() + problem.add_mapping("scaling", "scaling") + with pytest.raises(NotImplementedError, match="mapping table"): + reparametrize_log(problem, ["scaling"]) + + +def test_reparametrize_log_unresolvable_override(): + # offset1 is reparametrized, offset2 is not + problem = _problem() + problem["offset2"].lb = 1 + problem["offset2"].ub = 10 + problem["offset2"].estimate = True + with pytest.raises( + NotImplementedError, match=r"observable placeholder 1 of .*offset2" + ): + reparametrize_log(problem, ["offset1"]) + + # non-positive numeric override + problem = _problem() + problem.measurements[1].observable_parameters = [0] + with pytest.raises(NotImplementedError, match=r"\['0'\]"): + reparametrize_log(problem, ["offset1"]) + + +def test_reparametrize_log_noise_override(): + """Observable placeholders in noise formulas and noise overrides.""" + problem = _problem() + problem.observables[0].noise_formula = sympify_petab( + "noiseParameter1_obs_A * offset_A" + ) + for m in problem.measurements: + m.noise_parameters = ["sigma"] + problem.add_parameter("sigma", nominal_value=0.1, lb=1e-3, ub=1) + assert not problem.validate().has_errors() + + reparametrize_log(problem, ["sigma", "offset1", "offset2"]) + assert problem.observables[0].noise_formula == ( + 10 ** _sym("noiseParameter1_obs_A") * 10 ** _sym("offset_A") + ) + assert [m.noise_parameters for m in problem.measurements] == [ + [_sym("log10_sigma")], + [_sym("log10_sigma")], + ] + assert not problem.validate().has_errors() + + +@pytest.mark.parametrize("base", [10, "e"]) +@pytest.mark.parametrize( + "prior", + [ + (C.LOG_NORMAL, [0.5, 1.5]), + (C.LOG_LAPLACE, [0.5, 1.5]), + (C.LOG_UNIFORM, [1e-2, 1e2]), + ], +) +def test_reparametrize_log_prior_densities(base, prior): + """The derived prior over q is the pushforward of the prior over p.""" + problem = _problem() + ( + problem["scaling"].prior_distribution, + problem["scaling"].prior_parameters, + ) = prior + old_dist = problem["scaling"].prior_dist + id_map = reparametrize_log(problem, ["scaling"], base=base) + new_dist = problem[id_map["scaling"]].prior_dist + + b = 10 if base == 10 else math.e + thetas = np.linspace(-1.5, 1.5, 7) + xs = b**thetas + # density transformation: pi_q(theta) = pi_p(b^theta) * b^theta * ln(b) + expected = old_dist.pdf(xs) * xs * math.log(b) + np.testing.assert_allclose(new_dist.pdf(thetas), expected, rtol=1e-10) From 0afcdffe6377016af47d5154ab0e0701f8a000bd Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 15 Oct 2026 14:47:46 +0200 Subject: [PATCH 2/6] Improve error for placeholders that can't be log-transformed Report all affected placeholders at once, by ID, together with the reparametrized parameters overriding them. Separate overrides by parameters that are not reparametrized (which can be fixed by including them in `parameter_ids`) from non-positive numbers. Co-Authored-By: Claude Opus 5.5 --- petab/v2/reparametrize.py | 78 +++++++++++++++++++++++++--------- tests/v2/test_conversion.py | 10 ++++- tests/v2/test_reparametrize.py | 24 ++++++++--- 3 files changed, 84 insertions(+), 28 deletions(-) diff --git a/petab/v2/reparametrize.py b/petab/v2/reparametrize.py index 8dc62560..d2734766 100644 --- a/petab/v2/reparametrize.py +++ b/petab/v2/reparametrize.py @@ -330,44 +330,80 @@ def _reparametrize_overrides( log_b, ) -> None: """Log-transform placeholders overridden by reparametrized parameters.""" + observables = {o.id: o for o in problem.observables} + + def get_placeholder(obs_id: str, type_: str, i: int) -> sp.Symbol: + observable = observables[obs_id] + if type_ == "observable": + return observable.observable_placeholders[i] + return observable.noise_placeholders[i] + + def quote(values) -> str: + return ", ".join(f"`{v}'" for v in values) + slots_to_transform = [] + errors = [] for slot, values in _get_override_slots(problem).items(): - if not any( - isinstance(v, sp.Symbol) and v.name in id_map for v in values - ): + reparametrized = sorted( + v.name + for v in values + if isinstance(v, sp.Symbol) and v.name in id_map + ) + if not reparametrized: continue - invalid = { - v + non_positive = sorted( + float(v) for v in values if v.is_number and not v > 0 + ) + others = sorted( + str(v) for v in values - if not ( - (isinstance(v, sp.Symbol) and v.name in id_map) - or (v.is_number and v > 0) + if not v.is_number + and not (isinstance(v, sp.Symbol) and v.name in id_map) + ) + if not non_positive and not others: + slots_to_transform.append(slot) + continue + + obs_id, type_, i = slot + error = ( + f"* {type_.capitalize()} placeholder `{get_placeholder(*slot)}' " + f"of observable `{obs_id}' (overridden by " + f"{quote(reparametrized)}) also receives:" + ) + if others: + error += ( + f"\n - parameters that are not reparametrized: " + f"{quote(others)} (consider including them in " + "`parameter_ids')" ) - } - if invalid: - obs_id, type_, i = slot - raise NotImplementedError( - f"Cannot reparametrize {type_} placeholder {i + 1} of " - f"observable `{obs_id}': all its overrides must be " - "reparametrized parameters or positive numbers, but got " - f"{sorted(map(str, invalid))}." + if non_positive: + error += "\n - non-positive numbers: " + ", ".join( + f"{v:g}" for v in non_positive ) - slots_to_transform.append(slot) + errors.append(error) + + if errors: + raise NotImplementedError( + "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" + + "\n".join(errors) + + "\nSplitting observables is not supported. Exclude the " + "affected parameters or split the observables manually." + ) if not slots_to_transform: return # placeholder -> b^placeholder in formulas - observables = {o.id: o for o in problem.observables} for obs_id, type_, i in slots_to_transform: observable = observables[obs_id] + placeholder = get_placeholder(obs_id, type_, i) if type_ == "observable": - placeholder = observable.observable_placeholders[i] observable.formula = observable.formula.subs( placeholder, b_expr**placeholder ) - else: - placeholder = observable.noise_placeholders[i] # observable placeholders may also occur in the noise formula observable.noise_formula = observable.noise_formula.subs( placeholder, b_expr**placeholder diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 7afd9562..8fc8ca2c 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -476,10 +476,14 @@ def test_petab1to2_reparametrize_all_skips_unresolvable(tmp_path): ) # s1 has a prior, so it must be reparametrized with ( - pytest.raises(NotImplementedError, match=r"\['0(\.0)?', 's2'\]"), + pytest.raises(NotImplementedError) as excinfo, pytest.warns(UserWarning, match="s2 were not reparametrized"), ): petab1to2(yaml_file, reparametrize="all") + msg = str(excinfo.value) + assert "(overridden by `s1')" in msg + assert "not reparametrized: `s2'" in msg + assert "non-positive numbers: 0" in msg @pytest.mark.filterwarnings( @@ -500,7 +504,9 @@ def test_petab1to2_reparametrize_overrides_unresolvable(tmp_path): ], overrides=["s1", "s2"], ) - with pytest.raises(NotImplementedError, match=r"\['s2'\]"): + with pytest.raises( + NotImplementedError, match=r"not reparametrized: `s2' \(" + ): petab1to2(yaml_file, reparametrize="priors") diff --git a/tests/v2/test_reparametrize.py b/tests/v2/test_reparametrize.py index 9a8994b6..163b35e2 100644 --- a/tests/v2/test_reparametrize.py +++ b/tests/v2/test_reparametrize.py @@ -231,16 +231,30 @@ def test_reparametrize_log_unresolvable_override(): problem["offset2"].lb = 1 problem["offset2"].ub = 10 problem["offset2"].estimate = True - with pytest.raises( - NotImplementedError, match=r"observable placeholder 1 of .*offset2" - ): + with pytest.raises(NotImplementedError) as excinfo: reparametrize_log(problem, ["offset1"]) + msg = str(excinfo.value) + assert ( + "Observable placeholder `offset_A' of observable `obs_A' " + "(overridden by `offset1') also receives:" in msg + ) + assert "not reparametrized: `offset2' (consider including" in msg + assert "non-positive" not in msg - # non-positive numeric override + # non-positive numeric overrides, in several placeholders, reported at + # once problem = _problem() + problem.measurements[0].noise_parameters = ["offset1"] problem.measurements[1].observable_parameters = [0] - with pytest.raises(NotImplementedError, match=r"\['0'\]"): + 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 "not reparametrized" not in msg def test_reparametrize_log_noise_override(): From 724c3d05c89bb34256343a5950ec4e850e443099 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 15 Oct 2026 16:20:33 +0200 Subject: [PATCH 3/6] Use `math.log10` in `reparametrize_log` for base 10 `math.log(x) / math.log(10)` is inexact for powers of 10 (e.g. `2.9999999999999996` for 1000), which gave bounds and log-uniform prior parameters that differ from those computed with `log10` elsewhere. Co-Authored-By: Claude Opus 5.5 --- petab/v2/reparametrize.py | 22 ++++++++++++++-------- tests/v2/test_reparametrize.py | 5 +++++ 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/petab/v2/reparametrize.py b/petab/v2/reparametrize.py index d2734766..43e33e07 100644 --- a/petab/v2/reparametrize.py +++ b/petab/v2/reparametrize.py @@ -20,13 +20,19 @@ __all__ = ["reparametrize_log"] #: x-space priors and the corresponding priors of log_b(x), given as -#: (prior type, function(prior parameters, ln(b)) -> prior parameters) +#: (prior type, function(prior parameters, log_b, ln(b)) -> prior parameters) _LOG_PRIOR_TO_THETA = { - C.LOG_NORMAL: (C.NORMAL, lambda pars, ln_b: [p / ln_b for p in pars]), - C.LOG_LAPLACE: (C.LAPLACE, lambda pars, ln_b: [p / ln_b for p in pars]), + C.LOG_NORMAL: ( + C.NORMAL, + lambda pars, log_b, ln_b: [p / ln_b for p in pars], + ), + C.LOG_LAPLACE: ( + C.LAPLACE, + lambda pars, log_b, ln_b: [p / ln_b for p in pars], + ), C.LOG_UNIFORM: ( C.UNIFORM, - lambda pars, ln_b: [math.log(p) / ln_b for p in pars], + lambda pars, log_b, ln_b: [log_b(p) for p in pars], ), } @@ -98,10 +104,13 @@ def reparametrize_log( priors = priors or {} if base == 10: ln_b = math.log(10) + # not `math.log(x) / ln_b`, which is inexact for powers of 10 + log_b = math.log10 b_expr = sp.Integer(10) prefix = "log10_" elif base == "e": ln_b = 1.0 + log_b = math.log b_expr = sp.E prefix = "log_" else: @@ -124,9 +133,6 @@ def reparametrize_log( ) model = problem.model - def log_b(x: float) -> float: - return math.log(x) / ln_b - new_parameters = {} model_entities = [] for p in parameter_ids: @@ -181,7 +187,7 @@ def log_b(x: float) -> float: prior_distribution, prior_parameters = None, [] elif par.prior_distribution in _LOG_PRIOR_TO_THETA: prior_distribution, f = _LOG_PRIOR_TO_THETA[par.prior_distribution] - prior_parameters = f(par.prior_parameters, ln_b) + prior_parameters = f(par.prior_parameters, log_b, ln_b) else: raise ValueError( f"Prior `{par.prior_distribution}' of parameter `{p}' " diff --git a/tests/v2/test_reparametrize.py b/tests/v2/test_reparametrize.py index 163b35e2..36c0233f 100644 --- a/tests/v2/test_reparametrize.py +++ b/tests/v2/test_reparametrize.py @@ -108,6 +108,11 @@ def test_reparametrize_log_parameter_table(base): assert scaling.prior_distribution is None assert scaling.name is None + if base == 10: + # powers of 10 are mapped exactly + assert [k1.lb, k1.ub] == [-3.0, 3.0] + assert k2.prior_parameters == [-1.0, 1.0] + assert not problem.validate().has_errors() From 83d34029f86ff3f3153f0ea15a20f6e388383ce9 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 15 Oct 2026 21:16:43 +0200 Subject: [PATCH 4/6] Define reparametrized model parameters by initial assignments In `reparametrize_log`, define `p = base^q` by an initial assignment instead of an assignment rule. `q` is a parameter table parameter and therefore constant during the simulation, so `p` keeps its `constant=true` attribute and stays a constant entity in the sense of the PEtab v2 spec (e.g., it can still be a condition target). An existing initial assignment for `p` is replaced. Co-Authored-By: Claude Opus 5.5 --- petab/v2/reparametrize.py | 22 +++++++++++----------- tests/v2/test_conversion.py | 2 +- tests/v2/test_reparametrize.py | 21 +++++++++++++-------- 3 files changed, 25 insertions(+), 20 deletions(-) diff --git a/petab/v2/reparametrize.py b/petab/v2/reparametrize.py index 43e33e07..e1a397b1 100644 --- a/petab/v2/reparametrize.py +++ b/petab/v2/reparametrize.py @@ -58,7 +58,8 @@ def reparametrize_log( Parameters without a positive nominal value get no nominal value for ``q``. * Model: If ``p`` is a model parameter, ``q`` is added to the model, - and ``p`` becomes the target of the assignment rule ``p = base^q``. + and ``p`` is defined by the initial assignment ``p = base^q``, + replacing any existing initial assignment for ``p``. Only SBML models are supported. * Observable and noise formulas, and condition table target values: ``p`` is substituted by ``base^q``. @@ -251,7 +252,7 @@ def subs(expr: sp.Basic | None) -> sp.Basic | None: # update the model for p in model_entities: - _add_sbml_log_rule( + _add_sbml_log_assignment( model.sbml_model, p, id_map[p], @@ -444,27 +445,26 @@ def transform(value: sp.Basic) -> sp.Basic: ) -def _add_sbml_log_rule( +def _add_sbml_log_assignment( model: libsbml.Model, parameter_id: str, new_id: str, b_expr: sp.Expr, value: float | None, ) -> None: - """Add parameter ``new_id`` and the rule ``parameter_id = b^new_id``.""" + """Add parameter ``new_id`` and the initial assignment + ``parameter_id = b^new_id``.""" add_sbml_parameter(model, new_id, value=value, constant=True) - par = model.getParameter(parameter_id) - check(par.setConstant(False)) - # an assignment rule and an initial assignment must not have the same - # target + # any existing initial assignment is overridden by the parameter table + # value anyway if model.getInitialAssignmentBySymbol(parameter_id) is not None: model.removeInitialAssignment(parameter_id) - rule = model.createAssignmentRule() - check(rule.setVariable(parameter_id)) + ia = model.createInitialAssignment() + check(ia.setSymbol(parameter_id)) formula = f"exp({new_id})" if b_expr == sp.E else f"10^({new_id})" ast = libsbml.parseL3Formula(formula) if ast is None: raise RuntimeError(f"Failed to parse `{formula}'.") - check(rule.setMath(ast)) + check(ia.setMath(ast)) diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 8fc8ca2c..350b2bba 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -336,7 +336,7 @@ def test_petab1to2_reparametrize_priors( # k1 = b^q in the model sbml_model = problem.model.sbml_model - assert sbml_model.getAssignmentRuleByVariable("k1") is not None + assert sbml_model.getInitialAssignmentBySymbol("k1") is not None assert sbml_model.getParameter(new_id).getConstant() # the v2 prior of the new parameter is the v1 prior density over the diff --git a/tests/v2/test_reparametrize.py b/tests/v2/test_reparametrize.py index 36c0233f..9fa319df 100644 --- a/tests/v2/test_reparametrize.py +++ b/tests/v2/test_reparametrize.py @@ -144,21 +144,24 @@ def test_reparametrize_log_expressions_and_model(): # model sbml_model = problem.model.sbml_model - assert sbml_model.getParameter("k1").getConstant() is False + assert sbml_model.getParameter("k1").getConstant() is True q = sbml_model.getParameter("log10_k1") assert q.getConstant() is True assert q.getValue() == pytest.approx(math.log10(2)) - rule = sbml_model.getAssignmentRuleByVariable("k1") - assert rule is not None - assert sympify_petab(rule.getFormula()) == sympify_petab("10^log10_k1") + assert sbml_model.getAssignmentRuleByVariable("k1") is None + ia = sbml_model.getInitialAssignmentBySymbol("k1") + assert ia is not None + assert sympify_petab(libsbml.formulaToL3String(ia.getMath())) == ( + sympify_petab("10^log10_k1") + ) assert problem.model.is_valid() # k2 is not reparametrized, it's a condition target - assert sbml_model.getAssignmentRuleByVariable("k2") is None + assert sbml_model.getInitialAssignmentBySymbol("k2") is None assert not problem.validate().has_errors() -def test_reparametrize_log_initial_assignment_removed(): +def test_reparametrize_log_initial_assignment_replaced(): problem = Problem() problem.model = SbmlModel.from_antimony("k1 = 2; A' = -k1 * A; A = 1") sbml_model = problem.model.sbml_model @@ -167,8 +170,10 @@ def test_reparametrize_log_initial_assignment_removed(): ia.setMath(libsbml.parseL3Formula("3")) problem.add_parameter("k1", nominal_value=2, lb=1, ub=10) reparametrize_log(problem, ["k1"]) - assert sbml_model.getInitialAssignmentBySymbol("k1") is None - assert sbml_model.getAssignmentRuleByVariable("k1") is not None + assert sbml_model.getNumInitialAssignments() == 1 + ia = sbml_model.getInitialAssignmentBySymbol("k1") + assert libsbml.formulaToL3String(ia.getMath()) == "10^log10_k1" + assert problem.model.is_valid() def test_reparametrize_log_id_collision(): From fd3334da7b8d1dfc5e2fde8aec06901f528199f4 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 15 Oct 2026 21:56:52 +0200 Subject: [PATCH 5/6] Skip log-scaled parameters with a lower bound of 0 in `petab1to2` PEtab v1 allows a lower bound of 0 for log-scaled parameters if an initialization prior is given. On log scale, this would be `-inf`, for which `reparametrize_log` raises, and neither the v1 default prior nor `parameterScaleUniform` would be proper. Instead of failing the whole conversion, keep such parameters on linear scale and warn. Generated with AI assistance. Co-Authored-By: Claude Opus 5.5 --- petab/v2/petab1to2.py | 26 +++++++++++++++ tests/v2/test_conversion.py | 64 +++++++++++++++++++++++++++++++++++++ 2 files changed, 90 insertions(+) diff --git a/petab/v2/petab1to2.py b/petab/v2/petab1to2.py index ebce0256..0c97b9b7 100644 --- a/petab/v2/petab1to2.py +++ b/petab/v2/petab1to2.py @@ -88,6 +88,7 @@ def petab1to2( Further log-scaled parameters are reparametrized if they share a placeholder with reparametrized parameters in the measurement table. + Parameters with a lower bound of ``0`` are not reparametrized. :raises ValueError: If the input is invalid or, if ``validate`` is ``True``, if the input @@ -763,6 +764,8 @@ def _reparametrize_log_scaled( # 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 + zero_bound = [] for par_id, row in df.iterrows(): pscale = row.get(v1.C.PARAMETER_SCALE, v1.C.LIN) if pscale not in (v1.C.LOG, v1.C.LOG10): @@ -776,6 +779,18 @@ def _reparametrize_log_scaled( # A prior over the unscaled parameter. Its conversion already # preserves the objective function. continue + if ( + row.get(v1.C.LOWER_BOUND, 1) <= 0 + or row.get(v1.C.UPPER_BOUND, 1) <= 0 + ): + # v1 allows a lower bound of 0 for log-scaled parameters. + # This would be -inf on log scale, where neither the v1 default + # prior nor a parameterScaleUniform prior are proper. + if int(row[v1.C.ESTIMATE]) and ( + reparametrize == "all" or prior_type is not None + ): + zero_bound.append(par_id) + continue prior = None if not int(row[v1.C.ESTIMATE]): @@ -815,6 +830,17 @@ def _reparametrize_log_scaled( ): selected[base].add(par_id) + if zero_bound: + warnings.warn( + f"Log-scaled parameters {', '.join(zero_bound)} were not " + "reparametrized, because their bounds are not positive. " + "They are estimated on linear scale, and their 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 positive numbers. Add log-scaled parameters that # share a placeholder with selected parameters. If that is not possible, diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 350b2bba..e310d0be 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -363,6 +363,70 @@ def test_petab1to2_reparametrize_priors( np.testing.assert_allclose(v2_pdf, v1_pdf, rtol=1e-12) +@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", + [ + (v1.C.PARAMETER_SCALE_NORMAL, "-1;0.5", "priors"), + (v1.C.PARAMETER_SCALE_NORMAL, "-1;0.5", "all"), + (np.nan, np.nan, "all"), + ], +) +def test_petab1to2_reparametrize_zero_lower_bound( + tmp_path, prior_type, prior_parameters, reparametrize +): + """Parameters with a lower bound of 0 are skipped.""" + # 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", + v1.C.OBJECTIVE_PRIOR_TYPE: prior_type, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: prior_parameters, + } + yaml_file = _write_v1_problem( + tmp_path, + observables=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.OBSERVABLE_FORMULA: "B", + v1.C.NOISE_FORMULA: "0.1", + } + ], + measurements=[ + { + v1.C.OBSERVABLE_ID: "obs_B", + v1.C.SIMULATION_CONDITION_ID: "c0", + v1.C.TIME: 1.0, + v1.C.MEASUREMENT: 0.3, + } + ], + parameters=[k1], + ) + assert not v1.lint_problem(v1.Problem.from_yaml(yaml_file)) + + with pytest.warns(UserWarning, match=r"k1 were not reparametrized"): + problem = petab1to2(yaml_file, reparametrize=reparametrize) + + assert [p.id for p in problem.parameters] == ["k1"] + par = problem["k1"] + assert (par.lb, par.ub) == (0.0, _K1[v1.C.UPPER_BOUND]) + if pd.isna(prior_type): + # as with reparametrize="none" + assert par.prior_distribution == v2.C.UNIFORM + assert par.prior_parameters == [0.0, _K1[v1.C.UPPER_BOUND]] + else: + assert par.prior_distribution == v2.C.LOG_NORMAL + 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: From 4be4d4afe461ace065e8e5cc0ea60a8afdaae6f6 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 15 Oct 2026 22:35:20 +0200 Subject: [PATCH 6/6] Address review comments on log reparametrization * `petab1to2`: also warn about converted `parameterScaleUniform` priors with `reparametrize="none"`, which change the maximum a posteriori estimate as well * `petab1to2`: point to `reparametrize="none"` if reparametrization fails, and drop the `parameter_ids` hint from the `reparametrize_log` error, which doesn't apply there * `petab1to2`: warn once about parameters skipped because of placeholders, not once per base * `reparametrize_log`: ignore non-positive bounds of non-estimated parameters, and raise `NotImplementedError` for SBML models before L2V2, which have no initial assignments * tests: likelihood after reparametrizing overrides with both bases, and the benchmark collection with `reparametrize="all"` Generated with AI assistance. Co-Authored-By: Claude Opus 5.5 --- petab/v2/petab1to2.py | 56 +++++++++------- petab/v2/reparametrize.py | 27 +++++--- tests/v2/test_conversion.py | 115 ++++++++++++++++++++++++++++++--- tests/v2/test_reparametrize.py | 20 +++++- 4 files changed, 175 insertions(+), 43 deletions(-) diff --git a/petab/v2/petab1to2.py b/petab/v2/petab1to2.py index 0c97b9b7..66d8d01d 100644 --- a/petab/v2/petab1to2.py +++ b/petab/v2/petab1to2.py @@ -701,7 +701,11 @@ def update_prior(row) -> tuple[str, str | None]: has_scale_prior = parameter_df[v1.C.PARAMETER_SCALE].isin( [v1.C.LOG, v1.C.LOG10] ) & parameter_df[v1.C.OBJECTIVE_PRIOR_TYPE].isin( - [v1.C.PARAMETER_SCALE_NORMAL, v1.C.PARAMETER_SCALE_LAPLACE] + [ + v1.C.PARAMETER_SCALE_NORMAL, + v1.C.PARAMETER_SCALE_LAPLACE, + v1.C.PARAMETER_SCALE_UNIFORM, + ] ) if has_scale_prior.any(): warnings.warn( @@ -779,16 +783,13 @@ def _reparametrize_log_scaled( # A prior over the unscaled parameter. Its conversion already # preserves the objective function. continue - if ( - row.get(v1.C.LOWER_BOUND, 1) <= 0 - or row.get(v1.C.UPPER_BOUND, 1) <= 0 + if int(row[v1.C.ESTIMATE]) and ( + row[v1.C.LOWER_BOUND] <= 0 or row[v1.C.UPPER_BOUND] <= 0 ): # v1 allows a lower bound of 0 for log-scaled parameters. # This would be -inf on log scale, where neither the v1 default # prior nor a parameterScaleUniform prior are proper. - if int(row[v1.C.ESTIMATE]) and ( - reparametrize == "all" or prior_type is not None - ): + if reparametrize == "all" or prior_type is not None: zero_bound.append(par_id) continue @@ -834,7 +835,7 @@ def _reparametrize_log_scaled( warnings.warn( f"Log-scaled parameters {', '.join(zero_bound)} were not " "reparametrized, because their bounds are not positive. " - "They are estimated on linear scale, and their parameter-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.", @@ -847,8 +848,8 @@ def _reparametrize_log_scaled( # 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(): - excluded = set() changed = True while changed: changed = False @@ -866,22 +867,29 @@ def _reparametrize_log_scaled( 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 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, + ) for base, base_selected in selected.items(): - 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}, - ) + try: + 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 any(selected.values()) diff --git a/petab/v2/reparametrize.py b/petab/v2/reparametrize.py index e1a397b1..3f339518 100644 --- a/petab/v2/reparametrize.py +++ b/petab/v2/reparametrize.py @@ -92,9 +92,9 @@ def reparametrize_log( :return: A dictionary mapping old to new parameter IDs. :raises ValueError: - If a parameter cannot be reparametrized, e.g., because its lower bound - is not positive, or if its prior cannot be expressed in terms of the - new parameter. + If a parameter cannot be reparametrized, e.g., because it is estimated + and its lower bound is not positive, or if its prior cannot be + expressed in terms of the new parameter. :raises NotImplementedError: If the problem has features that are not supported yet, e.g., multiple or non-SBML models, mapping table entries for the @@ -152,8 +152,14 @@ def reparametrize_log( raise NotImplementedError( f"Model entity `{p}' is not an SBML parameter." ) + sbml_model = model.sbml_model + if (sbml_model.getLevel(), sbml_model.getVersion()) < (2, 2): + raise NotImplementedError( + "Reparametrization of model parameters requires " + "initial assignments (SBML L2V2 or later)." + ) model_entities.append(p) - if any(bound is not None and bound <= 0 for bound in (par.lb, par.ub)): + if par.estimate and (par.lb <= 0 or par.ub <= 0): # (bounds are required for estimated parameters) raise ValueError( f"Bounds of parameter `{p}' must be positive, " @@ -199,8 +205,9 @@ def reparametrize_log( new_parameters[p] = Parameter( id=p, name=par.name and f"{prefix.removesuffix('_')}({par.name})", - lb=None if par.lb is None else log_b(par.lb), - ub=None if par.ub is None else log_b(par.ub), + # bounds of non-estimated parameters are irrelevant + lb=log_b(par.lb) if par.lb is not None and par.lb > 0 else None, + ub=log_b(par.ub) if par.ub is not None and par.ub > 0 else None, nominal_value=nominal_value, estimate=par.estimate, prior_distribution=prior_distribution, @@ -380,8 +387,7 @@ def quote(values) -> str: if others: error += ( f"\n - parameters that are not reparametrized: " - f"{quote(others)} (consider including them in " - "`parameter_ids')" + f"{quote(others)}" ) if non_positive: error += "\n - non-positive numbers: " + ", ".join( @@ -396,8 +402,9 @@ def quote(values) -> str: "all overrides of a placeholder to be reparametrized parameters " "or positive numbers.\n" + "\n".join(errors) - + "\nSplitting observables is not supported. Exclude the " - "affected parameters or split the observables manually." + + "\nSplitting observables is not supported. Reparametrize " + "all or none of the parameters overriding a placeholder, " + "or split the observables manually." ) if not slots_to_transform: diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index e310d0be..ef1607d8 100644 --- a/tests/v2/test_conversion.py +++ b/tests/v2/test_conversion.py @@ -363,6 +363,87 @@ def test_petab1to2_reparametrize_priors( 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) + + @pytest.mark.filterwarnings( "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning", "ignore:Initialisation priors.*:UserWarning", @@ -568,10 +649,11 @@ def test_petab1to2_reparametrize_overrides_unresolvable(tmp_path): ], overrides=["s1", "s2"], ) - with pytest.raises( - NotImplementedError, match=r"not reparametrized: `s2' \(" - ): + with pytest.raises(NotImplementedError) as excinfo: petab1to2(yaml_file, reparametrize="priors") + msg = str(excinfo.value) + assert "Use `reparametrize='none'`" in msg + assert "not reparametrized: `s2'\n" in msg @pytest.mark.filterwarnings( @@ -618,11 +700,19 @@ def test_petab1to2_reparametrize_output_parameter(tmp_path): assert problem["obs_B"].formula == sympify_petab("10^log10_scaling * B") -def test_petab1to2_reparametrize_warnings_and_errors(tmp_path): +@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: v1.C.PARAMETER_SCALE_NORMAL, - v1.C.OBJECTIVE_PRIOR_PARAMETERS: "0;1", + v1.C.OBJECTIVE_PRIOR_TYPE: prior_type, + v1.C.OBJECTIVE_PRIOR_PARAMETERS: "-1;1", } yaml_file = _write_v1_problem( tmp_path, @@ -1034,8 +1124,15 @@ def test_petab1to2_validate(tmp_path): @pytest.mark.filterwarnings( "ignore:.*Parameter-scale priors of log-scaled parameters.*:UserWarning" ) +@pytest.mark.filterwarnings( + "ignore:Log-scaled parameters .* were not reparametrized.*:UserWarning" +) +@pytest.mark.filterwarnings( + "ignore:Parameter .* has no positive 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) @@ -1045,8 +1142,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 9fa319df..d717a57b 100644 --- a/tests/v2/test_reparametrize.py +++ b/tests/v2/test_reparametrize.py @@ -234,6 +234,24 @@ def test_reparametrize_log_invalid(): with pytest.raises(NotImplementedError, match="mapping table"): reparametrize_log(problem, ["scaling"]) + problem = _problem() + doc = problem.model.sbml_model.getSBMLDocument() + assert doc.setLevelAndVersion(2, 1) + with pytest.raises(NotImplementedError, match="L2V2"): + reparametrize_log(problem, ["k1"]) + + +def test_reparametrize_log_fixed_nonpositive_bounds(): + """Bounds of non-estimated parameters are irrelevant.""" + problem = _problem() + problem["offset2"].lb = 0 + problem["offset2"].ub = 10 + reparametrize_log(problem, ["offset1", "offset2"]) + par = problem["log10_offset2"] + assert par.lb is None + assert par.ub == 1 + assert par.nominal_value == pytest.approx(math.log10(6)) + def test_reparametrize_log_unresolvable_override(): # offset1 is reparametrized, offset2 is not @@ -248,7 +266,7 @@ def test_reparametrize_log_unresolvable_override(): "Observable placeholder `offset_A' of observable `obs_A' " "(overridden by `offset1') also receives:" in msg ) - assert "not reparametrized: `offset2' (consider including" in msg + assert "not reparametrized: `offset2'\n" in msg assert "non-positive" not in msg # non-positive numeric overrides, in several placeholders, reported at