diff --git a/petab/v2/petab1to2.py b/petab/v2/petab1to2.py index cd978e58..bd731691 100644 --- a/petab/v2/petab1to2.py +++ b/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- becomes log- 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( diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index d2d33bb1..a09acb22 100644 --- a/tests/v2/test_conversion.py +++ b/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,9 +568,6 @@ 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" @@ -367,6 +575,9 @@ def test_petab1to2_validate(tmp_path): @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."""