Add reparametrize_log and petab1to2(reparametrize=...)

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
 

doc/modules.rst

@@ -39,3 +39,4 @@ API Reference
    petab.v2.lint
    petab.v2.models
    petab.v2.petab1to2
+   petab.v2.reparametrize

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

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)

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_<id>`` (``log_<id>``) instead of
+        ``<id>``, 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())