diff --git a/CHANGELOG.md b/CHANGELOG.md index fd39f282..a0c59acd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,10 +29,17 @@ 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) +* `petab.v2.ConditionTable.to_df` writes non-numeric target values in PEtab + math syntax instead of sympy syntax (by @dweindl) **Features** * `petab1to2` and `petab_files_1to2` can skip validation via `validate=False` (by @FFroehlich in https://github.com/PEtab-dev/libpetab-python/pull/534) +* `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) ## 0.9 series diff --git a/doc/modules.rst b/doc/modules.rst index 08e150a3..4f79e8f4 100644 --- a/doc/modules.rst +++ b/doc/modules.rst @@ -39,3 +39,4 @@ API Reference petab.v2.lint petab.v2.models petab.v2.petab1to2 + petab.v2.reparametrize diff --git a/petab/v2/__init__.py b/petab/v2/__init__.py index e67f2dea..c4c810e9 100644 --- a/petab/v2/__init__.py +++ b/petab/v2/__init__.py @@ -37,3 +37,4 @@ ) from .lint import lint_problem # noqa: F401 from .models import MODEL_TYPE_PYSB, MODEL_TYPE_SBML, Model # noqa: F401 +from .reparametrize import reparametrize_log # noqa: F401 diff --git a/petab/v2/core.py b/petab/v2/core.py index 2efba572..8399ce06 100644 --- a/petab/v2/core.py +++ b/petab/v2/core.py @@ -553,7 +553,7 @@ def to_df(self) -> pd.DataFrame: record[C.TARGET_VALUE] = ( float(record[C.TARGET_VALUE]) if record[C.TARGET_VALUE].is_number - else str(record[C.TARGET_VALUE]) + else petab_math_str(record[C.TARGET_VALUE]) ) return ( pd.DataFrame(records) diff --git a/petab/v2/petab1to2.py b/petab/v2/petab1to2.py index 20263434..e214ec53 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]) ) @@ -260,6 +326,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 @@ -492,11 +565,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() @@ -600,10 +679,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, @@ -632,3 +734,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..5d198076 --- /dev/null +++ b/petab/v2/reparametrize.py @@ -0,0 +1,431 @@ +"""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." + ) + + extra = dict(par.model_extra or {}) + if name := extra.get(C.PARAMETER_NAME): + extra[C.PARAMETER_NAME] = f"{prefix.removesuffix('_')}({name})" + + new_parameters[p] = Parameter( + id=p, + 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, + **extra, + ) + + # 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 68685c37..bfa056c2 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 @@ -196,6 +198,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", [ @@ -259,6 +264,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.""" @@ -595,6 +918,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_core.py b/tests/v2/test_core.py index e33378e9..bdb3c804 100644 --- a/tests/v2/test_core.py +++ b/tests/v2/test_core.py @@ -64,6 +64,24 @@ def test_condition_table_round_trip(): assert conditions == conditions2 +def test_condition_table_round_trip_expression(tmp_path): + """Expressions are written in PEtab math syntax.""" + conditions = ConditionTable( + [ + Condition( + id="c1", changes=[Change(target_id="k", target_value="2^p")] + ) + ] + ) + assert conditions.to_df()[C.TARGET_VALUE].tolist() == ["2.0 ^ p"] + conditions.to_tsv(tmp_path / "conditions.tsv") + conditions2 = ConditionTable.from_tsv(tmp_path / "conditions.tsv") + assert ( + conditions2.conditions[0].changes[0].target_value + == conditions.conditions[0].changes[0].target_value + ) + + @pytest.mark.filterwarnings( "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning" ) diff --git a/tests/v2/test_reparametrize.py b/tests/v2/test_reparametrize.py new file mode 100644 index 00000000..02b84016 --- /dev/null +++ b/tests/v2/test_reparametrize.py @@ -0,0 +1,295 @@ +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.model_extra[C.PARAMETER_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 C.PARAMETER_NAME not in scaling.model_extra + + 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)