diff --git a/CHANGELOG.md b/CHANGELOG.md index 44569474..7ef13cb3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,19 @@ All notable changes to PyBNF are documented below. This project adheres to ## [Unreleased] ### Added +- **Scale-preserving PEtab v1→v2 conversion: `pybnf.petab.petab1to2_preserve_scale`.** The + official `petab.v2.petab1to2` **drops** the v1 `parameterScale` column (PEtab v2 removed it) and + only *warns* — so a `parameterScale = log10` estimated parameter carrying no objective prior (the + common case for a multi-decade kinetic parameter) converts to a *linear* `uniform_var` over the raw + bounds: the same argmin, but a far harder, worse-conditioned optimization than the log10 search the + modeler specified. This wrapper runs the standard converter and re-injects the dropped estimation + scale in the **v2-native** form — `priorDistribution = log-uniform` over each such parameter's + bounds — which PyBNF imports as a `loguniform_var` on the Log10 scale. Because the optimizer + objective excludes the prior, this sets only the search scale and initial sampling, not the + objective, so the fit stays the pure-MLE problem v1 specified; parameters petab1to2 already folded + into a prior (`parameterScale*Normal` → `log-normal`, …) are left untouched, as are linear ones. + `import_job` stays a pure v2 importer — the conversion is an explicit, named step, not a reach-back + to v1 in the read path. Intended as the opt-in migration `petab1to2` itself should offer. - **General-objective trust-region optimizer: `fit_type = gntr` (ADR-0068, #481).** Fills the last empty cell of the gradient-fitting (objective × curvature-model) matrix. `trf` gives a trust-region step with a `JᵀJ` (Gauss-Newton / empirical-Fisher) Hessian but only for an exact diff --git a/pybnf/petab/__init__.py b/pybnf/petab/__init__.py index a81edf46..18275323 100644 --- a/pybnf/petab/__init__.py +++ b/pybnf/petab/__init__.py @@ -25,6 +25,7 @@ read_condition_table, read_experiment_table, ) +from .convert import petab1to2_preserve_scale from .export import clean_model_for_petab, export_job, write_problem_yaml from .import_ import import_job, read_problem_yaml from .measurements import ( @@ -90,4 +91,6 @@ # job importer 'import_job', 'read_problem_yaml', + # scale-preserving v1 -> v2 converter + 'petab1to2_preserve_scale', ] diff --git a/pybnf/petab/convert.py b/pybnf/petab/convert.py new file mode 100644 index 00000000..c3a858aa --- /dev/null +++ b/pybnf/petab/convert.py @@ -0,0 +1,136 @@ +"""Scale-preserving PEtab v1 -> v2 conversion (wraps ``petab.v2.petab1to2``). + +``petab.v2.petab1to2`` drops the v1 ``parameterScale`` column -- PEtab v2 removed it, on +the view that estimation scale is an *optimiser* concern, not part of the problem spec -- +and only *warns* ("Parameter scales are not supported in PEtab v2"). For a parameter that +was ``parameterScale = log10`` **without** an objective prior (the common case: a +maximum-likelihood fit over a multi-decade kinetic parameter), the converted v2 problem is +read as *linear* ``uniform_var`` over the raw bounds. That is the same argmin, but a +vastly harder, worse-conditioned optimisation than the log10 search the modeller specified +-- and running the problem *as specified* is the whole point of a benchmark. + +petab1to2 already preserves scale where it is attached to an *objective prior* +(``parameterScale{Normal,Uniform,Laplace}`` -> ``log-normal`` / ``log-laplace`` / ...); it +is only the **bare** estimation scale it drops. This module wraps the standard converter +and re-injects that scale in the **v2-native** form: each bare log/log10 estimated +parameter gets ``priorDistribution = log-uniform`` over its bounds. PyBNF -- and any PEtab +v2 tool -- then reads a log search (a ``log-uniform`` prior maps to a ``loguniform_var`` on +the Log10 scale; :mod:`pybnf.petab.parameters`). PyBNF's optimiser objective *excludes* the +prior, so the log-uniform prior sets only the search scale and initial sampling, not the +objective -- the fit stays the pure-MLE problem v1 specified. + +This is the migration ``petab1to2`` should offer as an opt-in; it lives here as an +explicit, named converter so :func:`pybnf.petab.import_job` stays a pure v2 importer with +no reach-back to v1 in the read path. +""" + +import warnings +from pathlib import Path + +from ..printing import PybnfError + +#: v1 ``parameterScale`` values that mean "estimate in log space" (base-independent for a +#: *uniform* prior: uniform-in-ln and uniform-in-log10 are the same distribution over the +#: same bounds, and PyBNF searches log-uniform on its Log10 scale either way). +_LOG_SCALES = frozenset({'log', 'log10'}) + + +def petab1to2_preserve_scale(v1_yaml_path, out_dir): + """Convert a PEtab **v1** problem to **v2**, preserving log estimation scales. + + Runs :func:`petab.v2.petab1to2`, then rewrites the converted v2 parameter table so each + estimated v1 ``parameterScale`` in {``log``, ``log10``} that carries no prior gains a + v2-native ``priorDistribution = log-uniform`` over its bounds. Returns the ``Path`` to + the converted v2 ``problem.yaml``. + + Parameters petab1to2 already gave a prior (a v1 ``parameterScale*`` objective prior) are + left untouched; ``lin``-scale parameters stay plain linear ``uniform_var``. + """ + try: + import pandas as pd + import petab.v1 as petab_v1 + import petab.v1.C as C1 + from petab.v2.petab1to2 import petab1to2 + except ImportError as e: + raise PybnfError( + 'Scale-preserving PEtab v1->v2 conversion needs the petab library.', + "Install the PEtab extra: pip install 'pybnf[petab]'.") from e + + v1_yaml_path = Path(v1_yaml_path) + out_dir = Path(out_dir) + + # 1. Standard conversion. Silence the "Parameter scales are not supported" warning -- + # re-adding that scale is exactly this function's job. + with warnings.catch_warnings(): + warnings.simplefilter('ignore') + petab1to2(str(v1_yaml_path), str(out_dir)) + v2_yaml = _sole_yaml(out_dir) + + # 2. The v1 estimation scales, per parameterId (estimated log/log10 only). + v1_spec = petab_v1.yaml.load_yaml(str(v1_yaml_path)) + v1_pdf = petab_v1.get_parameter_df(str(v1_yaml_path.parent / v1_spec['parameter_file'])) + log_estimated = { + str(pid) + for pid, row in v1_pdf.iterrows() + if str(row.get(C1.PARAMETER_SCALE, C1.LIN)) in _LOG_SCALES + and _is_estimated(row.get(C1.ESTIMATE, 1)) + } + if not log_estimated: + return v2_yaml + + # 3. Re-inject the dropped scale as a v2-native log-uniform prior over the bounds. + v2_spec = petab_v1.yaml.load_yaml(str(v2_yaml)) + v2_param_path = out_dir / v2_spec['parameter_files'][0] + v2_pdf = pd.read_csv(v2_param_path, sep='\t') + inject_log_uniform_priors(v2_pdf, log_estimated) + v2_pdf.to_csv(v2_param_path, sep='\t', index=False) + return v2_yaml + + +def inject_log_uniform_priors(v2_pdf, log_estimated_ids): + """Give each v2 parameter row in ``log_estimated_ids`` a ``log-uniform`` prior in place. + + For every row whose ``parameterId`` is in ``log_estimated_ids`` **and** that carries no + prior yet, sets ``priorDistribution = log-uniform`` and ``priorParameters`` to its + ``[lowerBound, upperBound]``. Rows with an existing prior (a scale petab1to2 already + folded into one) and rows not in the set are left untouched. Mutates and returns + ``v2_pdf`` (a v2 parameter :class:`pandas.DataFrame`). + """ + import petab.v2.C as C2 + + for col in (C2.PRIOR_DISTRIBUTION, C2.PRIOR_PARAMETERS): + if col not in v2_pdf.columns: + v2_pdf[col] = '' + # petab1to2 emits an all-empty priorParameters as float64 (NaN); coerce to object + # so the string cells below don't raise a dtype error. NaNs still write as blank. + v2_pdf[col] = v2_pdf[col].astype('object') + for i, row in v2_pdf.iterrows(): + if str(row[C2.PARAMETER_ID]) not in log_estimated_ids: + continue + if _has_prior(row.get(C2.PRIOR_DISTRIBUTION)): + continue # petab1to2 already carried this scale into a prior -- don't clobber. + lb, ub = row[C2.LOWER_BOUND], row[C2.UPPER_BOUND] + v2_pdf.at[i, C2.PRIOR_DISTRIBUTION] = C2.LOG_UNIFORM + v2_pdf.at[i, C2.PRIOR_PARAMETERS] = f'{lb}{C2.PARAMETER_SEPARATOR}{ub}' + return v2_pdf + + +def _sole_yaml(out_dir): + """The single ``problem.yaml`` petab1to2 wrote into ``out_dir``.""" + yamls = sorted(Path(out_dir).glob('*.yaml')) + if not yamls: + raise PybnfError(f'petab1to2 wrote no problem.yaml into {out_dir}.') + return yamls[0] + + +def _is_estimated(value): + """v1 ``estimate`` cell (``1``/``0`` or truthy) -> whether the parameter is fit.""" + try: + return int(value) == 1 + except (TypeError, ValueError): + return str(value).strip().lower() in ('1', 'true') + + +def _has_prior(value): + """Whether a v2 ``priorDistribution`` cell already names a prior (non-empty, non-NaN).""" + return value is not None and str(value).strip().lower() not in ('', 'nan') diff --git a/tests/test_petab_convert.py b/tests/test_petab_convert.py new file mode 100644 index 00000000..217eece9 --- /dev/null +++ b/tests/test_petab_convert.py @@ -0,0 +1,74 @@ +"""Scale-preserving PEtab v1->v2 conversion (:mod:`pybnf.petab.convert`). + +Unit coverage for the log-scale re-injection that `petab1to2_preserve_scale` layers on top +of the standard `petab.v2.petab1to2` (which drops the v1 `parameterScale` column). The +end-to-end conversion is exercised against the real benchmark problems in the fitting +harness; here we pin the transform that could regress: a bare log/log10 estimated parameter +gains a v2 `log-uniform` prior over its bounds, a linear one does not, and an existing prior +is never clobbered. +""" +import pytest + +pd = pytest.importorskip('pandas') +pytest.importorskip('petab') + +from pybnf.petab.convert import _has_prior, _is_estimated, inject_log_uniform_priors + + +class TestInjectLogUniformPriors: + + def _petab1to2_shape(self): + # Exactly what petab1to2 emits: an all-empty priorParameters as float64, and no + # priorDistribution column at all. + return pd.DataFrame({ + 'parameterId': ['klog', 'klin'], + 'lowerBound': [1e-5, 0.0], + 'upperBound': [1e5, 5.0], + 'nominalValue': [0.02, 0.5], + 'estimate': [True, False], + 'priorParameters': [float('nan'), float('nan')], + }) + + def test_log_param_gets_log_uniform_over_its_bounds(self): + df = self._petab1to2_shape() + inject_log_uniform_priors(df, {'klog'}) + r = df.set_index('parameterId').loc['klog'] + assert r['priorDistribution'] == 'log-uniform' + lo, hi = (float(x) for x in r['priorParameters'].split(';')) + assert (lo, hi) == (1e-5, 1e5) + + def test_param_not_in_log_set_is_untouched(self): + df = self._petab1to2_shape() + inject_log_uniform_priors(df, {'klog'}) # klin is linear -> not in the set + r = df.set_index('parameterId').loc['klin'] + assert not _has_prior(r['priorDistribution']) + + def test_float64_priorparameters_column_is_coerced_not_raised(self): + # Regression: the string cell must not raise on petab1to2's float64 NaN column. + inject_log_uniform_priors(self._petab1to2_shape(), {'klog'}) # must not raise + + def test_existing_prior_is_not_clobbered(self): + df = self._petab1to2_shape() + df['priorDistribution'] = [None, 'normal'] # klin already carries a prior + df['priorParameters'] = df['priorParameters'].astype('object') + df.loc[df.parameterId == 'klin', 'priorParameters'] = '0;1' + inject_log_uniform_priors(df, {'klog', 'klin'}) # klin in set, but already priored + r = df.set_index('parameterId').loc['klin'] + assert r['priorDistribution'] == 'normal' and r['priorParameters'] == '0;1' + + +class TestHelpers: + + @pytest.mark.parametrize('value,expected', [ + (1, True), ('1', True), (True, True), ('true', True), + (0, False), ('0', False), (float('nan'), False), + ]) + def test_is_estimated(self, value, expected): + assert _is_estimated(value) is expected + + @pytest.mark.parametrize('value,expected', [ + ('normal', True), ('log-uniform', True), + (None, False), ('', False), ('nan', False), (float('nan'), False), + ]) + def test_has_prior(self, value, expected): + assert _has_prior(value) is expected