petab1to2: rescale noise and priors when converting log10 to log

petab/v2/petab1to2.py

@@ -2,6 +2,7 @@
 
 from __future__ import annotations
 
+import math
 import shutil
 import warnings
 from contextlib import suppress
@@ -405,20 +406,12 @@ def get_noise_dist(row):
 
             if trans == v1.C.LIN:
                 new_dist = dist
+            elif trans == v1.C.LOG10:
+                # see `get_noise_formula`
+                new_dist = f"{v1.C.LOG}-{dist}"
             else:
                 new_dist = f"{trans}-{dist}"
 
-            if new_dist == "log10-normal":
-                warnings.warn(
-                    f"Noise distribution `{new_dist}' for "
-                    f"observable `{row[v1.C.OBSERVABLE_ID]}'"
-                    f" is not supported in PEtab v2. "
-                    "Using `log-normal` instead.",
-                    # call to `petab1to2`
-                    stacklevel=9,
-                )
-                new_dist = v2.C.LOG_NORMAL
-
             if new_dist not in v2.C.NOISE_DISTRIBUTIONS:
                 raise NotImplementedError(
                     f"Noise distribution `{new_dist}' for "
@@ -428,7 +421,36 @@ def get_noise_dist(row):
 
             return new_dist
 
+        # PEtab v2 has no log10-based noise distributions. Since
+        #  ln(y) = ln(10) * log10(y), a normal (Laplace) distribution of
+        #  log10(y) with standard deviation (scale) sigma is a normal
+        #  (Laplace) distribution of ln(y) with ln(10) * sigma.
+        #  Therefore, log10-<dist> becomes log-<dist> with the noise
+        #  formula multiplied by ln(10).
+        def get_noise_formula(row):
+            formula = row[v1.C.NOISE_FORMULA]
+            if row[v1.C.OBSERVABLE_TRANSFORMATION] != v1.C.LOG10:
+                return formula
+            return f"log(10) * ({formula})"
+
         df[v2.C.NOISE_DISTRIBUTION] = df.apply(get_noise_dist, axis=1)
+        df[v1.C.NOISE_FORMULA] = df.apply(get_noise_formula, axis=1)
+
+        is_log10 = df[v1.C.OBSERVABLE_TRANSFORMATION] == v1.C.LOG10
+        if is_log10.any():
+            log10_ids = ", ".join(df.loc[is_log10, v1.C.OBSERVABLE_ID])
+            warnings.warn(
+                "PEtab v2 does not support log10-transformed observables. "
+                f"The noise distributions of observables {log10_ids} were "
+                "converted to log-normal or log-laplace, and their noise "
+                "formulas were multiplied by `log(10)`. Parameters in these "
+                "noise formulas (e.g., estimated noise parameters) keep their "
+                "PEtab v1 meaning, i.e., they refer to the log10-transformed "
+                "observable, not to its natural logarithm.",
+                # call to `petab1to2`
+                stacklevel=4,
+            )
+
         df.drop(columns=[v1.C.OBSERVABLE_TRANSFORMATION], inplace=True)
 
     def extract_placeholders(row: pd.Series, type_: str) -> str:
@@ -482,44 +504,52 @@ def v1v2_parameter_df(
         lambda x: str(bool(int(x))).lower()
     )
 
-    def update_prior(row):
-        """Convert prior to v2 format."""
+    def update_prior(row) -> tuple[str, str | None]:
+        """Convert prior type and prior parameters to v2 format."""
         prior_type = row.get(v1.C.OBJECTIVE_PRIOR_TYPE)
         if pd.isna(prior_type):
             prior_type = v1.C.UNIFORM
+        prior_pars = row.get(v1.C.OBJECTIVE_PRIOR_PARAMETERS)
 
         pscale = row.get(v1.C.PARAMETER_SCALE)
         if pd.isna(pscale):
             pscale = v1.C.LIN
 
         if prior_type not in v1.C.PARAMETER_SCALE_PRIOR_TYPES:
-            return prior_type
+            return prior_type, prior_pars
+
+        if pscale == v1.C.LOG10 and prior_type in (
+            v1.C.PARAMETER_SCALE_NORMAL,
+            v1.C.PARAMETER_SCALE_LAPLACE,
+        ):
+            # PEtab v2 has no log10-based priors. Since
+            #  ln(x) = ln(10) * log10(x), a normal (Laplace) distribution of
+            #  log10(x) with mean (location) mu and standard deviation
+            #  (scale) sigma is a normal (Laplace) distribution of ln(x)
+            #  with ln(10) * mu and ln(10) * sigma.
+            pscale = v1.C.LOG
+            prior_pars = v2.C.PARAMETER_SEPARATOR.join(
+                str(float(p) * math.log(10))
+                for p in str(prior_pars).split(v1.C.PARAMETER_SEPARATOR)
+            )
 
         new_prior_type = prior_type.removeprefix("parameterScale").lower()
         if pscale != v1.C.LIN:
             new_prior_type = f"{pscale}-{new_prior_type}"
 
-        if new_prior_type == "log10-normal":
-            warnings.warn(
-                f"Prior distribution `{new_prior_type}' for parameter "
-                f"`{row[v1.C.PARAMETER_ID]}' is not supported in PEtab v2. "
-                "Using `log-normal` instead.",
-                # call to `petab1to2`
-                stacklevel=9,
-            )
-            new_prior_type = v2.C.LOG_NORMAL
-
         if new_prior_type not in v2.C.PRIOR_DISTRIBUTIONS:
             raise NotImplementedError(
                 f"PEtab v2 does not support prior type `{new_prior_type}' "
                 f"required for parameter `{row[v1.C.PARAMETER_ID]}'."
             )
 
-        return new_prior_type
+        return new_prior_type, prior_pars
 
     # update parameterScale*-priors
     if v1.C.OBJECTIVE_PRIOR_TYPE in df.columns:
-        df[v1.C.OBJECTIVE_PRIOR_TYPE] = df.apply(update_prior, axis=1)
+        df[[v1.C.OBJECTIVE_PRIOR_TYPE, v1.C.OBJECTIVE_PRIOR_PARAMETERS]] = (
+            df.apply(update_prior, axis=1, result_type="expand")
+        )
 
     # rename objectivePrior* to prior*
     df.rename(

tests/v2/test_conversion.py

@@ -1,12 +1,16 @@
 import logging
 from pathlib import Path
 
+import numpy as np
 import pandas as pd
 import pytest
 
 from petab import v1, v2
+from petab.v1.calculate import calculate_llh as calculate_llh_v1
 from petab.v1.models.sbml_model import SbmlModel
+from petab.v1.priors import Prior
 from petab.v2 import Problem
+from petab.v2.calculate import calculate_llh as calculate_llh_v2
 from petab.v2.petab1to2 import petab1to2, v1v2_observable_df
 
 
@@ -31,6 +35,213 @@ def test_v1v2_observable_df_noise_distribution():
     ]
 
 
+def test_v1v2_observable_df_log10():
+    """Test that log10-transformed observables are converted to
+    log-transformed ones with the noise formula scaled by ln(10)."""
+    observable_df = pd.DataFrame(
+        data={
+            v1.C.OBSERVABLE_ID: ["obs1", "obs2", "obs3"],
+            v1.C.OBSERVABLE_FORMULA: ["a", "b", "c"],
+            v1.C.NOISE_FORMULA: [
+                "noiseParameter1_obs1 + 0.1",
+                "noiseParameter1_obs2 + 0.1",
+                "noiseParameter1_obs3 + 0.1",
+            ],
+            v1.C.OBSERVABLE_TRANSFORMATION: [
+                v1.C.LOG,
+                v1.C.LOG10,
+                v1.C.LOG10,
+            ],
+            v1.C.NOISE_DISTRIBUTION: [
+                v1.C.NORMAL,
+                v1.C.NORMAL,
+                v1.C.LAPLACE,
+            ],
+        }
+    ).set_index(v1.C.OBSERVABLE_ID)
+
+    with pytest.warns(UserWarning, match=r"observables obs2, obs3 were"):
+        new_df = v1v2_observable_df(observable_df)
+
+    assert list(new_df[v2.C.NOISE_DISTRIBUTION]) == [
+        v2.C.LOG_NORMAL,
+        v2.C.LOG_NORMAL,
+        v2.C.LOG_LAPLACE,
+    ]
+    assert list(new_df[v2.C.NOISE_FORMULA]) == [
+        "noiseParameter1_obs1 + 0.1",
+        "log(10) * (noiseParameter1_obs2 + 0.1)",
+        "log(10) * (noiseParameter1_obs3 + 0.1)",
+    ]
+    assert list(new_df[v2.C.NOISE_PLACEHOLDERS]) == [
+        "noiseParameter1_obs1",
+        "noiseParameter1_obs2",
+        "noiseParameter1_obs3",
+    ]
+
+
+def _write_v1_problem(
+    path: Path,
+    *,
+    observables: list[dict],
+    measurements: list[dict],
+    parameters: list[dict],
+) -> str:
+    """Write a PEtab v1 problem for the conversion A -> B with the given
+    tables (as lists of rows) to ``path`` and return the YAML file name."""
+    problem = v1.Problem(
+        model=SbmlModel.from_antimony(
+            "r1: A -> B; k1 * A; A = 1; B = 0; k1 = 0.5"
+        ),
+        condition_df=v1.get_condition_df(
+            pd.DataFrame([{v1.C.CONDITION_ID: "c0", "A": 1.0}])
+        ),
+        observable_df=v1.get_observable_df(pd.DataFrame(observables)),
+        measurement_df=v1.get_measurement_df(pd.DataFrame(measurements)),
+        parameter_df=v1.get_parameter_df(pd.DataFrame(parameters)),
+    )
+    return problem.to_files_generic(path)
+
+
+_K1 = {
+    v1.C.PARAMETER_ID: "k1",
+    v1.C.PARAMETER_SCALE: v1.C.LIN,
+    v1.C.LOWER_BOUND: 1e-3,
+    v1.C.UPPER_BOUND: 1e2,
+    v1.C.NOMINAL_VALUE: 0.5,
+    v1.C.ESTIMATE: 1,
+}
+
+
+@pytest.mark.parametrize("noise_distribution", [v1.C.NORMAL, v1.C.LAPLACE])
+@pytest.mark.parametrize("noise_placeholder", [False, True])
+def test_petab1to2_log10_noise_llh(
+    tmp_path, noise_distribution, noise_placeholder
+):
+    """Test that the conversion of log10-transformed observables preserves
+    the likelihood."""
+    if noise_placeholder:
+        noise_formula = "noiseParameter1_obs_B + offset"
+        noise_parameters = "sigma_B"
+        noise_parameter_rows = [
+            {
+                v1.C.PARAMETER_ID: "sigma_B",
+                v1.C.PARAMETER_SCALE: v1.C.LIN,
+                v1.C.LOWER_BOUND: 1e-3,
+                v1.C.UPPER_BOUND: 1e2,
+                v1.C.NOMINAL_VALUE: 0.1,
+                v1.C.ESTIMATE: 1,
+            },
+            {
+                v1.C.PARAMETER_ID: "offset",
+                v1.C.PARAMETER_SCALE: v1.C.LIN,
+                v1.C.LOWER_BOUND: 0,
+                v1.C.UPPER_BOUND: 1,
+                v1.C.NOMINAL_VALUE: 0.05,
+                v1.C.ESTIMATE: 0,
+            },
+        ]
+    else:
+        noise_formula = "0.1"
+        noise_parameters = ""
+        noise_parameter_rows = []
+
+    yaml_file = _write_v1_problem(
+        tmp_path,
+        observables=[
+            {
+                v1.C.OBSERVABLE_ID: "obs_B",
+                v1.C.OBSERVABLE_FORMULA: "B",
+                v1.C.OBSERVABLE_TRANSFORMATION: v1.C.LOG10,
+                v1.C.NOISE_FORMULA: noise_formula,
+                v1.C.NOISE_DISTRIBUTION: noise_distribution,
+            }
+        ],
+        measurements=[
+            {
+                v1.C.OBSERVABLE_ID: "obs_B",
+                v1.C.SIMULATION_CONDITION_ID: "c0",
+                v1.C.TIME: t,
+                v1.C.MEASUREMENT: y,
+                v1.C.NOISE_PARAMETERS: noise_parameters,
+            }
+            for t, y in ((1.0, 0.3), (2.0, 0.5))
+        ],
+        parameters=[_K1, *noise_parameter_rows],
+    )
+    v1_problem = v1.Problem.from_yaml(yaml_file)
+    with pytest.warns(UserWarning, match=r"observables obs_B were"):
+        v2_problem = petab1to2(yaml_file)
+
+    def simulate(measurement_df):
+        return measurement_df.rename(
+            columns={v1.C.MEASUREMENT: v1.C.SIMULATION}
+        ).assign(**{v1.C.SIMULATION: [0.35, 0.45]})
+
+    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 llh_v2 == pytest.approx(llh_v1, rel=1e-12)
+
+
+@pytest.mark.filterwarnings(
+    "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning"
+)
+@pytest.mark.parametrize(
+    "parameter_scale, prior_type",
+    [
+        (v1.C.LOG10, v1.C.PARAMETER_SCALE_NORMAL),
+        (v1.C.LOG10, v1.C.PARAMETER_SCALE_LAPLACE),
+        (v1.C.LOG, v1.C.PARAMETER_SCALE_NORMAL),
+        (v1.C.LOG, v1.C.PARAMETER_SCALE_LAPLACE),
+        (v1.C.LIN, v1.C.PARAMETER_SCALE_NORMAL),
+    ],
+)
+def test_petab1to2_parameter_scale_priors(
+    tmp_path, parameter_scale, prior_type
+):
+    """Test that parameterScale{Normal,Laplace} priors are converted to
+    priors with the same density."""
+    k1 = _K1 | {
+        v1.C.PARAMETER_SCALE: parameter_scale,
+        v1.C.OBJECTIVE_PRIOR_TYPE: prior_type,
+        v1.C.OBJECTIVE_PRIOR_PARAMETERS: "-0.5;0.7",
+    }
+    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],
+    )
+    v1_prior = Prior.from_par_dict(k1, type_="objective")
+    v2_prior = petab1to2(yaml_file)["k1"].prior_dist
+
+    x = np.logspace(-3, 2, 21)
+    assert v2_prior.pdf(x) == pytest.approx(v1_prior.pdf(x), rel=1e-12)
+
+
 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."""
@@ -357,16 +568,16 @@ def test_petab1to2_validate(tmp_path):
     )
 
 
-@pytest.mark.filterwarnings(
-    "ignore:.*Using `log-normal` instead.*:UserWarning"
-)
 @pytest.mark.filterwarnings(
     "ignore:.*Initialisation priors in parameter table are not supported.*:"
     "UserWarning"
 )
 @pytest.mark.filterwarnings(
     "ignore:.*Parameter scales are not supported in PEtab v2.*:UserWarning"
 )
+@pytest.mark.filterwarnings(
+    "ignore:.*does not support log10-transformed observables.*:UserWarning"
+)
 @parametrize_or_skip
 def test_benchmark_collection(problem_id):
     """Test that we can upgrade all benchmark collection models."""