From dd0cea086e556524b39cdfe60ee1013f4e618b83 Mon Sep 17 00:00:00 2001 From: Bill Hlavacek Date: Wed, 4 Nov 2026 07:36:08 -0700 Subject: [PATCH 1/2] fix(petab): undo petab 0.10.0's log(10) factor on a converted log10 observable's sigma petab 0.10.0 converts a v1 log10 observable to v2's natural-log family and writes its noise formula as log(10) * (), which is exact for that family. petab1to2_preserve_scale states the scale differently (ADR-0073): it re-injects observableTransformation = log10 and resets noiseDistribution to the v1 linear base, and under that column the v1 sigma is the right one. Together the two gave a log10 residual with sigma ln 10 times the stated one. A numeric sigma imported as "sigma = formula log(10) * (0.1)" with no error (petab 0.9.0 imports "fix_at 0.1"); a placeholder sigma failed at import, which is the failure main's CI showed. restore_v1_noise_formulas puts each log10 observable's v1 noiseFormula back. It removes exactly the log(10) * (...) wrapper, around the v1 formula as petab's own v1 reader prints it, leaves an unchanged formula alone, and refuses any other rewrite with both formulas named. petab 0.10.0's new warning for this case joins the silenced patterns, since the converter repairs what it says. The output is now the same under petab 0.9 and 0.10. petab 0.10.0 also converts the v1 priors 0.9.x refused (logNormal, logLaplace, log10 parameterScaleUniform / parameterScaleLaplace). The converter already rewrites every declared prior from v1, so nothing changes there; on petab >= 0.10 those cases now run through the three existing prior oracles (hand-written v1 and v2 densities, PyBNF's imported prior, libpetab's own reading), and the refusal test runs only where petab still refuses. Tests: tests/test_petab_convert.py passes on petab 0.10.0 (156 passed, 8 skipped) and 0.9.0 (140 passed). With the restore disabled, both log10 full-conversion tests fail on 0.10.0. Signed-off-by: Bill Hlavacek --- .../+petab-log10-noise-formula.fixed.md | 8 ++ ...imported-as-the-native-lognormal-family.md | 26 +++++ pybnf/petab/convert.py | 96 +++++++++++++++++-- tests/test_petab_convert.py | 91 ++++++++++++++++-- 4 files changed, 205 insertions(+), 16 deletions(-) create mode 100644 changelog.d/+petab-log10-noise-formula.fixed.md diff --git a/changelog.d/+petab-log10-noise-formula.fixed.md b/changelog.d/+petab-log10-noise-formula.fixed.md new file mode 100644 index 00000000..7fa7f575 --- /dev/null +++ b/changelog.d/+petab-log10-noise-formula.fixed.md @@ -0,0 +1,8 @@ +- **Converting a PEtab v1 problem with petab 0.10.0 no longer inflates a log10 observable's + sigma by ln 10.** petab 0.10.0 writes a v1 `log10` observable's noise formula as + `log(10) * ()`, which is right for the natural-log noise family petab writes but + wrong under the `observableTransformation = log10` column that + `petab1to2_preserve_scale` restores. A log10 observable with a numeric sigma then imported + with a sigma 2.3 times the stated one, with no error. The converter now puts the v1 noise + formula back, and refuses any other rewrite of it. It gives the same result under petab + 0.9 and 0.10. diff --git a/docs/adr/0073-petab-v1-log10-observabletransformation-is-re-injected-as-a-preserved-column-and-imported-as-the-native-lognormal-family.md b/docs/adr/0073-petab-v1-log10-observabletransformation-is-re-injected-as-a-preserved-column-and-imported-as-the-native-lognormal-family.md index 83928b66..fef9b223 100644 --- a/docs/adr/0073-petab-v1-log10-observabletransformation-is-re-injected-as-a-preserved-column-and-imported-as-the-native-lognormal-family.md +++ b/docs/adr/0073-petab-v1-log10-observabletransformation-is-re-injected-as-a-preserved-column-and-imported-as-the-native-lognormal-family.md @@ -229,3 +229,29 @@ PyBNF searches that parameter linearly. The objective is unchanged. This is the declared-prior counterpart of the bare-scale loss #548 fixed, and it is not addressed here; the "Parameter scales are not supported" warning that would have covered it is still silenced. None of the benchmark problems above has such a row. + +## Addendum (2026-11-04): petab 0.10.0 rescales the noise formula, and the converter undoes it + +petab 0.10.0 converts a v1 `log10` observable faithfully in v2's own terms: `log10-normal` +becomes `log-normal` and the noise formula becomes `log(10) * ()`, since a normal +`log10(y)` with sigma is a normal `ln(y)` with `ln(10)` sigma. This decision states the scale +in the re-injected `observableTransformation` column over the linear base family instead, and +under that column the v1 sigma is the right one. With petab 0.10.0 the converter therefore +produced `normal` + `log10` + `log(10) * (sigma)`: a log10 residual whose sigma was `ln 10` +times the stated one. A numeric sigma imported as `sigma = formula log(10) * (0.1)` with no +error, where petab 0.9.0 imports `fix_at 0.1`; a placeholder sigma failed at import. + +The converter now puts each `log10` observable's v1 `noiseFormula` back +(`restore_v1_noise_formulas`). It removes exactly the `log(10) * (...)` wrapper petab 0.10.0 +adds, around the v1 formula as petab's own v1 reader prints it, and leaves an unchanged +formula alone, so the output is the same under petab 0.9 and 0.10. Any other rewrite is +refused, naming both formulas, rather than carried over: the cell sets sigma. The new +wording of petab's warning for this case joins the silenced patterns, because the converter +repairs all of what it describes. + +The representation chosen above is kept rather than switching to petab 0.10.0's: it is the +same on both petab releases pyproject allows, and the importer reads it today. petab 0.10.0 +also converts the v1 priors it used to refuse (`logNormal`, `logLaplace`, a `log10` +`parameterScaleUniform` / `parameterScaleLaplace`). Nothing changes for them here, since the +converter rewrites every declared prior from v1; on petab 0.10.0 they are checked by the same +three oracles as every other prior. diff --git a/pybnf/petab/convert.py b/pybnf/petab/convert.py index 5f94c47f..13572593 100644 --- a/pybnf/petab/convert.py +++ b/pybnf/petab/convert.py @@ -33,9 +33,10 @@ ``log-uniform(e^a;e^b)`` (no warning). So petab1to2's prior is not kept for any row whose v1 author declared one: :func:`v2_prior_from_v1` rewrites each from the v1 row, mapping every v1 prior type on every scale to the v2 prior with the same distribution over the parameter. -The cases petab1to2 refuses outright (``logNormal``, ``logLaplace``, a ``log10`` -``parameterScaleUniform`` / ``parameterScaleLaplace``) still stop the conversion there; -ADR-0073's 2026-09-25 addendum tabulates every case. +petab < 0.10.0 refuses some cases outright (``logNormal``, ``logLaplace``, a ``log10`` +``parameterScaleUniform`` / ``parameterScaleLaplace``), and the conversion stops there; +petab 0.10.0 converts them, and the converter rewrites them from v1 like every other declared +prior. ADR-0073's 2026-09-25 addendum tabulates every case. **observableTransformation (the residual scale; issue #499, ADR-0073).** A v1 observable with ``observableTransformation = log10`` fits the residual on the log10 scale (with the @@ -45,15 +46,19 @@ log10 form**; ``petab1to2`` therefore cannot carry a v1 ``log10`` transformation faithfully. petab < 0.9.0 dropped it entirely (a missing ``return`` left ``noiseDistribution`` blank, so the observable imported as a linear Gaussian and the fit optimised the wrong objective); -petab >= 0.9.0 substitutes the natural-log family (``log10-normal`` -> ``log-normal``, with a -warning), which silently rescales sigma by ``ln 10``. Since v2 has no faithful -representation, we re-inject ``observableTransformation`` as a **preserved extra column** on +petab 0.9.x substitutes the natural-log family (``log10-normal`` -> ``log-normal``, with a +warning), which silently rescales sigma by ``ln 10``; petab 0.10.0 makes that substitution +exact by also writing the noise formula as ``log(10) * ()``. Rather than depend on +which petab ran, we re-inject ``observableTransformation`` as a **preserved extra column** on the v2 observables table and reset ``noiseDistribution`` to the v1 **linear base family** (``normal`` / ``laplace``), so the residual scale is stated in exactly one place: PyBNF's importer reads the column to select the noise family's additive scale (``lin`` / ``log10`` / ``log``; :mod:`pybnf.petab.observables`, :mod:`pybnf.petab.import_`) and refuses a log transformation over an already-log distribution as a contradiction (issue #679), while other -PEtab v2 tools ignore the unknown column (it passes v2 lint). This is directly parallel to +PEtab v2 tools ignore the unknown column (it passes v2 lint). A log10 observable's +``noiseFormula`` is put back to the v1 one as well, undoing petab 0.10.0's ``log(10)`` +factor, since under the re-injected column its sigma is on the log10 scale, as v1 stated it +(:func:`restore_v1_noise_formulas`; any other rewrite is refused). This is directly parallel to the parameterScale re-injection above -- the same "re-add the scale petab1to2 dropped" migration, on the observable axis. @@ -109,6 +114,10 @@ # log10 + normal observable -> log-normal. Every log observable's noiseDistribution is # reset to its v1 base and its observableTransformation re-injected (#499, #679). r"Noise distribution `log10-normal' for observable `[^']*' is not supported in PEtab v2\.", + # The same, as petab >= 0.10.0 words it: log10- -> log- with the noise + # formula multiplied by log(10). The distribution and transformation are repaired as + # above, and the multiplication is undone by restore_v1_noise_formulas. + r'PEtab v2 does not support log10-transformed observables\.', # log10 parameterScaleNormal -> log-normal with the numbers unchanged. Every declared # prior is rewritten from its v1 row by v2_prior_from_v1 (#893). r"Prior distribution `log10-normal' for parameter `[^']*' is not supported in PEtab v2\.", @@ -212,14 +221,21 @@ def petab1to2_preserve_scale(v1_yaml_path, out_dir): # family (petab >= 0.9.0; blank before), so that column is reset to the v1 linear base # (normal / laplace) -- the importer refuses a log transformation over an already-log # distribution as a contradiction (issue #679). + # petab >= 0.10.0 also rewrites a log10 observable's noiseFormula as + # ``log(10) * ()``, which is right for the natural-log family it writes and + # wrong under the log10 transformation re-injected here: the sigma would be ln 10 times + # the one v1 stated, with no error. So that formula is put back (restore_v1_noise_formulas). scales = _v1_observable_scales(v1_yaml_path, v1_spec) if scales: transformations = {oid: transformation for oid, (transformation, _) in scales.items()} distributions = {oid: distribution for oid, (_, distribution) in scales.items()} + log10_formulas = _v1_noise_formulas(v1_yaml_path, v1_spec, { + oid for oid, transformation in transformations.items() if transformation == 'log10'}) for obs_file in v2_spec.get('observable_files', []): v2_obs_path = out_dir / obs_file v2_odf = pd.read_csv(v2_obs_path, sep='\t') inject_observable_transformations(v2_odf, transformations, distributions) + restore_v1_noise_formulas(v2_odf, log10_formulas) v2_odf.to_csv(v2_obs_path, sep='\t', index=False) return v2_yaml @@ -499,6 +515,72 @@ def inject_observable_transformations(v2_odf, transformations, distributions=Non return v2_odf +#: How petab >= 0.10.0 rewrites a log10 observable's noise formula (``v1v2_observable_df``): +#: ``log(10) * ()``, the factor that turns a log10-scale sigma into a natural-log one. +_LN10_NOISE_PREFIX = 'log(10) * (' +_LN10_NOISE_SUFFIX = ')' + + +def restore_v1_noise_formulas(v2_odf, v1_formulas): + """Put back the v1 ``noiseFormula`` of each log10 observable petab1to2 rescaled. + + ``v1_formulas`` is ``{observableId: v1 noiseFormula text}`` for the observables whose v1 + ``observableTransformation`` is ``log10``. PyBNF states that scale in the re-injected + ``observableTransformation`` column over the linear base family + (:func:`inject_observable_transformations`), so the noise formula has to be the v1 one, + whose sigma is on the log10 scale. petab < 0.10.0 leaves the formula as it was; petab + >= 0.10.0 writes ``log(10) * ()`` for its natural-log family, and that wrapper is + removed here. Any other rewrite raises :class:`PybnfError` rather than being guessed at: + the cell decides the sigma, and a wrong sigma is a wrong objective, not a failed fit. + Mutates and returns ``v2_odf``. + """ + col = 'noiseFormula' + if not v1_formulas or col not in v2_odf.columns: + return v2_odf + v2_odf[col] = v2_odf[col].astype('object') + for i, row in v2_odf.iterrows(): + oid = str(row['observableId']) + if oid not in v1_formulas: + continue + v1 = v1_formulas[oid] + v2 = '' if _cell(row[col]) is None else str(row[col]) + if v2.strip() == v1.strip(): + continue + if v2.startswith(_LN10_NOISE_PREFIX) and v2.endswith(_LN10_NOISE_SUFFIX): + inner = v2[len(_LN10_NOISE_PREFIX):-len(_LN10_NOISE_SUFFIX)] + if inner.strip() == v1.strip(): + v2_odf.at[i, col] = inner + continue + raise PybnfError( + f"Converting PEtab v1 observable '{oid}' (observableTransformation log10): " + f"petab1to2 rewrote its noiseFormula from {v1!r} to {v2!r}, which is not the " + f"'log(10) * (...)' rescaling this converter knows how to undo.", + hint='The noise formula sets sigma, so it is not carried over on a guess. Check ' + 'what this petab release does to log10 observables, or convert the problem ' + 'with petab 0.9 or 0.10.') + return v2_odf + + +def _v1_noise_formulas(v1_yaml_path, v1_spec, observable_ids): + """``{observableId: noiseFormula}`` from the v1 observable file(s), as text, for each id + in ``observable_ids``. Read with petab's own v1 reader, the one petab1to2 reads with, so a + numeric cell prints here exactly as petab1to2 printed it inside ``log(10) * (...)``.""" + import petab.v1 as petab_v1 + import petab.v1.C as C1 + + formulas = {} + if not observable_ids: + return formulas + for problem in v1_spec.get('problems', []): + for obs_file in problem.get('observable_files', []): + odf = petab_v1.get_observable_df(str(Path(v1_yaml_path).parent / obs_file)) + for oid, row in odf.iterrows(): + if str(oid) in observable_ids: + value = row.get(C1.NOISE_FORMULA) + formulas[str(oid)] = '' if _cell(value) is None else str(value) + return formulas + + def _v1_observable_scales(v1_yaml_path, v1_spec): """``{observableId: ('log' | 'log10', 'normal' | 'laplace')}`` for every v1 observable with a log residual scale: its transformation and the linear base family of its v1 diff --git a/tests/test_petab_convert.py b/tests/test_petab_convert.py index a9e18a39..95fb9a38 100644 --- a/tests/test_petab_convert.py +++ b/tests/test_petab_convert.py @@ -32,6 +32,7 @@ inject_log_uniform_priors, inject_observable_transformations, petab1to2_preserve_scale, + restore_v1_noise_formulas, v2_prior_from_v1, ) from pybnf.printing import PybnfError @@ -140,6 +141,27 @@ def test_converted_log10_problem_imports_as_lognormal(self, tmp_path): assert 'objective = chi_sq' not in conf # the linear (wrong) import is gone assert 'loguniform_var = k' in conf # parameterScale=log10 preserved too + def test_a_log10_observables_sigma_imports_as_v1_stated_it(self, tmp_path): + # v1: log10(y) ~ N(log10(h), 0.1). petab 0.10.0 writes the faithful v2 natural-log + # form, log(y) ~ N(log(h), log(10) * 0.1), but the converter states the scale as + # observableTransformation = log10 over `normal`, where the v1 sigma is the right one. + # Left alone, the import read sigma = formula log(10) * (0.1): a log10 residual with a + # sigma 2.30x too large, imported with no error (petab 0.9.0 imported fix_at 0.1). + from pybnf.petab import import_job + yaml = _write_v1_problem(tmp_path / 'v1', 'log10') + obs = tmp_path / 'v1' / 'observables.tsv' + obs.write_text(obs.read_text().replace('\tnoiseParameter1_obs_V\n', '\t0.1\n')) + meas = tmp_path / 'v1' / 'measurements.tsv' + meas.write_text(meas.read_text().replace('\t0.1\n', '\t\n')) + v2_yaml = petab1to2_preserve_scale(str(yaml), str(tmp_path / 'v2')) + row = pd.read_csv(v2_yaml.parent / 'observables.tsv', sep='\t').set_index( + 'observableId').loc['obs_V'] + assert (row['observableTransformation'], row['noiseDistribution']) == ('log10', 'normal') + assert float(row['noiseFormula']) == 0.1 + conf = (import_job(v2_yaml, tmp_path / 'imported') / 'imported.conf').read_text() + assert 'noise_model = lognormal, sigma = fix_at 0.1' in conf + assert 'log(10)' not in conf + def test_converted_log_problem_imports_as_lnnormal(self, tmp_path): # v1 ``log`` means natural log. It reaches Gaussian(LN) through the explicit # ``lnnormal`` token, distinct from PyBNF's log10 ``lognormal``. @@ -265,10 +287,11 @@ def test_converted_problem_imports_as_a_preequilibrated_scan(self, tmp_path, con ('lin', '', '0.01;100'), ('log', '', '0.5;2'), ('log10', '', '-2;2'), ] -# The v1 priors petab1to2 (petab 0.9.0) refuses outright. Each is exactly representable in v2 -# (v2_prior_from_v1 maps it; TestV2PriorFromV1 pins those values), but petab1to2 raises before -# the converter sees its output. If a petab release starts accepting one, move it into -# _ACCEPTED_V1_PRIORS: the converter rewrites it from v1 whatever petab1to2 writes. +# The v1 priors petab1to2 refused outright through petab 0.9.x. Each is exactly representable in +# v2 (v2_prior_from_v1 maps it; TestV2PriorFromV1 pins those values), but petab 0.9.x raises +# before the converter sees its output. petab 0.10.0 converts every one of them, so on such a +# release they join _ACCEPTED_V1_PRIORS' oracles: the converter rewrites each from v1 whatever +# petab1to2 writes. pyproject allows both releases, so each is tested where it applies. _REFUSED_UPSTREAM_V1_PRIORS = [ ('lin', 'logNormal', '-1;0.5'), ('log', 'logNormal', '-1;0.5'), ('log10', 'logNormal', '-1;0.5'), @@ -278,6 +301,16 @@ def test_converted_problem_imports_as_a_preequilibrated_scan(self, tmp_path, con ] +def _petab_release(): + import petab + return tuple(int(part) for part in petab.__version__.split('.')[:2]) + + +_PETAB1TO2_REFUSES_SOME_V1_PRIORS = _petab_release() < (0, 10) +_CONVERTED_V1_PRIORS = _ACCEPTED_V1_PRIORS + ( + [] if _PETAB1TO2_REFUSES_SOME_V1_PRIORS else _REFUSED_UPSTREAM_V1_PRIORS) + + def _v1_logpdf(prior_type, parameters, scale, theta): """Log density over theta of an (untruncated) PEtab v1 objective prior, written out by hand from the v1 specification: uniform/normal/laplace act on theta, logNormal/logLaplace on @@ -400,7 +433,7 @@ def test_blank_prior_type_under_parameters_is_parameter_scale_uniform(self, tmp_ assert fp.type == 'loguniform_var' assert (fp.p1, fp.p2) == (pytest.approx(0.01, rel=1e-15), pytest.approx(100.0, rel=1e-15)) - @pytest.mark.parametrize('scale,prior_type,parameters', _ACCEPTED_V1_PRIORS) + @pytest.mark.parametrize('scale,prior_type,parameters', _CONVERTED_V1_PRIORS) def test_converted_prior_has_the_v1_distribution(self, tmp_path, scale, prior_type, parameters): # Oracle 1: the v2 row, read by hand from the v2 specification, has exactly the v1 @@ -423,7 +456,7 @@ def test_converted_prior_has_the_v1_distribution(self, tmp_path, scale, prior_ty offset = pybnf[finite] - v1[finite] np.testing.assert_allclose(offset, offset[0], rtol=0, atol=1e-7) - @pytest.mark.parametrize('scale,prior_type,parameters', _ACCEPTED_V1_PRIORS) + @pytest.mark.parametrize('scale,prior_type,parameters', _CONVERTED_V1_PRIORS) def test_libpetab_reads_the_v1_and_v2_rows_alike(self, tmp_path, scale, prior_type, parameters): # Oracle 3: libpetab's own reading of the v1 row (petab.v1.priors.Prior) and of the @@ -438,10 +471,12 @@ def test_libpetab_reads_the_v1_and_v2_rows_alike(self, tmp_path, scale, prior_ty (param,) = petab_v2.Problem.from_yaml(v2_yaml).parameters np.testing.assert_allclose(param.prior_dist.pdf(_THETAS), v1, rtol=1e-9, atol=0) + @pytest.mark.skipif(not _PETAB1TO2_REFUSES_SOME_V1_PRIORS, + reason='petab >= 0.10.0 converts these; _CONVERTED_V1_PRIORS checks them') @pytest.mark.parametrize('scale,prior_type,parameters', _REFUSED_UPSTREAM_V1_PRIORS) def test_priors_petab1to2_refuses_still_stop_the_conversion(self, tmp_path, scale, prior_type, parameters): - # Loud, not silent: petab1to2 raises NotImplementedError (log10-uniform/-laplace) or a + # Loud, not silent: petab 0.9.x raises NotImplementedError (log10-uniform/-laplace) or a # pydantic ValidationError, a ValueError (logNormal/logLaplace, not v2 spellings). with pytest.raises((NotImplementedError, ValueError)): _convert_one(tmp_path, scale, prior_type, parameters) @@ -480,7 +515,7 @@ def test_an_exact_petab1to2_table_is_left_byte_for_byte(self, tmp_path): == (tmp_path / 'plain' / 'parameters.tsv').read_bytes()) -# A narrow box that every prior in _ACCEPTED_V1_PRIORS crosses on at least one side, so the +# A narrow box that every prior in _CONVERTED_V1_PRIORS crosses on at least one side, so the # bounds truncate each one; nominalValue 3 lies inside every truncated support. The points # straddle both walls and the interior support edges (e^0.5 = 1.649 for the natural-log # uniform cases). @@ -522,7 +557,7 @@ class TestTruncatedV1PriorConversion: point for point, not just up to a constant: the same support, the same walls, the same mass. A second oracle is libpetab's own truncated reading of the v1 row.""" - @pytest.mark.parametrize('scale,prior_type,parameters', _ACCEPTED_V1_PRIORS) + @pytest.mark.parametrize('scale,prior_type,parameters', _CONVERTED_V1_PRIORS) def test_imported_prior_is_the_truncated_v1_prior(self, tmp_path, scale, prior_type, parameters): from petab.v1.priors import Prior @@ -890,6 +925,44 @@ def test_without_distributions_the_column_is_untouched(self): assert df.set_index('observableId').loc['obs_V', 'noiseDistribution'] == 'log-normal' +class TestRestoreV1NoiseFormulas: + """petab 0.10.0 writes a log10 observable's noiseFormula as ``log(10) * ()``; under + the re-injected log10 transformation the v1 formula is the right sigma, so it goes back.""" + + def _odf(self, formulas): + return pd.DataFrame({'observableId': list(formulas), 'noiseFormula': list(formulas.values())}) + + def test_petab_0_10s_wrapper_is_removed(self): + odf = self._odf({'obs_V': 'log(10) * (noiseParameter1_obs_V)', 'obs_c': 'log(10) * (0.1)'}) + restore_v1_noise_formulas(odf, {'obs_V': 'noiseParameter1_obs_V', 'obs_c': '0.1'}) + assert list(odf['noiseFormula']) == ['noiseParameter1_obs_V', '0.1'] + + def test_petab_0_9s_untouched_formula_is_left_alone(self): + odf = self._odf({'obs_V': 'sigma_V'}) + restore_v1_noise_formulas(odf, {'obs_V': 'sigma_V'}) + assert list(odf['noiseFormula']) == ['sigma_V'] + + def test_a_v1_formula_that_is_itself_scaled_by_log_10_is_unwrapped_once(self): + # v1 wrote the factor; petab 0.9.x keeps it, petab 0.10.0 adds a second one. + v1 = 'log(10) * (s)' + for written in (v1, 'log(10) * (log(10) * (s))'): + odf = self._odf({'obs_V': written}) + restore_v1_noise_formulas(odf, {'obs_V': v1}) + assert list(odf['noiseFormula']) == [v1], written + + def test_observables_that_are_not_log10_are_not_touched(self): + # A natural-log observable keeps its formula under every petab release, and a linear + # one is never in the map; neither is rewritten even if its text looks wrapped. + odf = self._odf({'obs_ln': 'log(10) * (s)', 'obs_lin': 'log(10) * (s)'}) + restore_v1_noise_formulas(odf, {}) + assert list(odf['noiseFormula']) == ['log(10) * (s)', 'log(10) * (s)'] + + def test_any_other_rewrite_is_refused_naming_both_formulas(self): + odf = self._odf({'obs_V': '2.302585 * (sigma_V)'}) + with pytest.raises(PybnfError, match=r"obs_V.*'sigma_V'.*'2\.302585 \* \(sigma_V\)'"): + restore_v1_noise_formulas(odf, {'obs_V': 'sigma_V'}) + + class TestHelpers: @pytest.mark.parametrize('value,expected', [ From 1478403d02853252748fe20dfb901a8b7f96f43a Mon Sep 17 00:00:00 2001 From: Bill Hlavacek Date: Wed, 4 Nov 2026 08:12:56 -0700 Subject: [PATCH 2/2] fix(petab): read the converted observables table as text Independent review of the previous commit found a loud regression, reproduced by running it: on petab 0.9 with pandas < 3.1, a log10 observable whose numeric sigma has 17 significant digits was refused. The converter re-read the v2 observables table with pandas' default float parse, which moves about 40% of 17-digit values by an ulp, so the formula no longer matched the v1 text petab's own reader printed, and restore_v1_noise_formulas refused it with a message blaming petab1to2. The table is now read with dtype=str and keep_default_na=False, so nothing in it is reparsed; the hint no longer tells the user to switch petab release. Also from the review: the claim that the output is the same under petab 0.9 and 0.10 is narrowed to the noise formula (log10 + laplace, and some pre-equilibration shapes, still differ between the releases, loudly or as described below), and the test reads petab's version from package metadata, since petab.__version__ is deprecated. The reviewer's tests are added (tests/test_petab_convert_review.py). They check the objective PyBNF builds from a converted, imported problem against libpetab's own v1 likelihood (calculate_llh) across numeric, placeholder, parameter and prediction-dependent sigmas, several observables and two files, and the refusal of an unknown rewrite through the whole converter. One strict xfail remains, on petab 0.9 only: petab 0.9's petab1to2 drops a pre-equilibration whose condition changes nothing, with no warning. That is upstream and pre-existing, and is not addressed here. Tests: test_petab_convert.py + the review file pass on petab 0.10.0 / pandas 3.1.0 (173 passed, 8 skipped) and petab 0.9.0 / pandas 3.0.5 (156 passed, 1 xfailed). Signed-off-by: Bill Hlavacek --- .../+petab-log10-noise-formula.fixed.md | 4 +- pybnf/petab/convert.py | 10 +- tests/test_petab_convert.py | 4 +- tests/test_petab_convert_review.py | 267 ++++++++++++++++++ 4 files changed, 277 insertions(+), 8 deletions(-) create mode 100644 tests/test_petab_convert_review.py diff --git a/changelog.d/+petab-log10-noise-formula.fixed.md b/changelog.d/+petab-log10-noise-formula.fixed.md index 7fa7f575..36793701 100644 --- a/changelog.d/+petab-log10-noise-formula.fixed.md +++ b/changelog.d/+petab-log10-noise-formula.fixed.md @@ -4,5 +4,5 @@ wrong under the `observableTransformation = log10` column that `petab1to2_preserve_scale` restores. A log10 observable with a numeric sigma then imported with a sigma 2.3 times the stated one, with no error. The converter now puts the v1 noise - formula back, and refuses any other rewrite of it. It gives the same result under petab - 0.9 and 0.10. + formula back, and refuses any other rewrite of it, so the noise formula is the same under + petab 0.9 and 0.10. diff --git a/pybnf/petab/convert.py b/pybnf/petab/convert.py index 13572593..6d70c8df 100644 --- a/pybnf/petab/convert.py +++ b/pybnf/petab/convert.py @@ -233,7 +233,9 @@ def petab1to2_preserve_scale(v1_yaml_path, out_dir): oid for oid, transformation in transformations.items() if transformation == 'log10'}) for obs_file in v2_spec.get('observable_files', []): v2_obs_path = out_dir / obs_file - v2_odf = pd.read_csv(v2_obs_path, sep='\t') + # As text: a float parse can move a 17-digit noise formula by an ulp (pandas + # < 3.1), which then no longer matches its v1 text in the restore below. + v2_odf = pd.read_csv(v2_obs_path, sep='\t', dtype=str, keep_default_na=False) inject_observable_transformations(v2_odf, transformations, distributions) restore_v1_noise_formulas(v2_odf, log10_formulas) v2_odf.to_csv(v2_obs_path, sep='\t', index=False) @@ -555,9 +557,9 @@ def restore_v1_noise_formulas(v2_odf, v1_formulas): f"Converting PEtab v1 observable '{oid}' (observableTransformation log10): " f"petab1to2 rewrote its noiseFormula from {v1!r} to {v2!r}, which is not the " f"'log(10) * (...)' rescaling this converter knows how to undo.", - hint='The noise formula sets sigma, so it is not carried over on a guess. Check ' - 'what this petab release does to log10 observables, or convert the problem ' - 'with petab 0.9 or 0.10.') + hint='The noise formula sets sigma, so it is not carried over on a guess. This ' + 'converter knows what petab 0.9 and 0.10 write for a log10 observable; ' + 'check what the installed petab release writes.') return v2_odf diff --git a/tests/test_petab_convert.py b/tests/test_petab_convert.py index 95fb9a38..a3715d88 100644 --- a/tests/test_petab_convert.py +++ b/tests/test_petab_convert.py @@ -302,8 +302,8 @@ def test_converted_problem_imports_as_a_preequilibrated_scan(self, tmp_path, con def _petab_release(): - import petab - return tuple(int(part) for part in petab.__version__.split('.')[:2]) + from importlib.metadata import version + return tuple(int(part) for part in version('petab').split('.')[:2]) _PETAB1TO2_REFUSES_SOME_V1_PRIORS = _petab_release() < (0, 10) diff --git a/tests/test_petab_convert_review.py b/tests/test_petab_convert_review.py new file mode 100644 index 00000000..836064f1 --- /dev/null +++ b/tests/test_petab_convert_review.py @@ -0,0 +1,267 @@ +"""Independent review of the petab 0.10 noise-formula restore in ``petab1to2_preserve_scale``. + +Written by a reviewer who did not write the change. The author's tests check the converted +``noiseFormula`` cell and the imported ``.conf`` text. These tests go one step further and +score the objective PyBNF actually builds from the converted problem, against an oracle that +never touches PyBNF's code: + +* **libpetab's own v1 likelihood** (``petab.v1.calculate.calculate_llh``) on the *original v1* + tables, for every noise source libpetab can evaluate (number, placeholder, fixed parameter), + on the ``log10`` observable the change repairs and on its ``log`` / ``lin`` siblings, one + observable table and two. +* **the v1 Gaussian-on-log10 density written out by hand**, for a prediction-dependent sigma + (``a + b * V``), which libpetab's calculator cannot evaluate. + +PyBNF drops the data-only constants of the likelihood, so the oracle is that +``PyBNF objective - v1 negative log-likelihood`` is the same number for every simulation. A +sigma that is ``ln 10`` too large (the bug) changes the quadratic term by a factor 5.3 and +breaks that by orders of magnitude more than the tolerance. + +It also drives the refusal end to end (a petab1to2 rewrite the converter does not recognise +must stop the conversion), and pins two places where petab 0.9 and 0.10 still differ. +""" +import math +import re +import warnings +from types import SimpleNamespace + +import numpy as np +import pytest + +pd = pytest.importorskip('pandas') +petab = pytest.importorskip('petab') + +from pybnf.printing import PybnfError + +_SBML = """\ + + + + + + + {extra} + + + + kV + + + + +""" + +_OBS_HEADER = ('observableId\tobservableFormula\tobservableTransformation\tnoiseDistribution' + '\tnoiseFormula') +_TIMES = (0.0, 1.0, 2.0) +_DATA = {'V': (5.0, 3.0, 2.0), 'W': (4.0, 2.5, 1.5)} + + +def _petab_minor(): + from importlib.metadata import version + return tuple(int(p) for p in version('petab').split('.')[:2]) + + +def _pandas_minor(): + return tuple(int(p) for p in pd.__version__.split('.')[:2]) + + +def _write_v1(root, obs_files, noise_parameters=None, fixed=(), conditions='conditionId\nc0\n', + preequilibration=None, estimated=()): + """A v1 problem. ``obs_files`` maps a file name to ``[(id, species, transformation, + noiseFormula)]``; every observable is the bare species it names, measured at _TIMES. + ``noise_parameters`` maps an observable id to its measurement ``noiseParameters`` cell; + ``fixed`` lists ``(parameterId, value)`` non-estimated parameters, ``estimated`` the + ``(parameterId, nominalValue)`` of estimated linear ones.""" + root.mkdir(parents=True, exist_ok=True) + extra = ''.join(f'' + for pid, _ in (*fixed, *estimated)) + (root / 'model.xml').write_text(_SBML.format(extra=extra)) + for name, rows in obs_files.items(): + (root / name).write_text('\n'.join( + [_OBS_HEADER] + [f'{oid}\t{sp}\t{tr}\tnormal\t{nf}' for oid, sp, tr, nf in rows]) + + '\n') + noise_parameters = noise_parameters or {} + header = ['observableId', 'simulationConditionId', 'measurement', 'time', 'noiseParameters'] + if preequilibration: + header.insert(1, 'preequilibrationConditionId') + lines = ['\t'.join(header)] + for rows in obs_files.values(): + for oid, species, _, _ in rows: + for t, y in zip(_TIMES, _DATA[species]): + cells = [oid, 'c0', str(y), str(t), noise_parameters.get(oid, '')] + if preequilibration: + cells.insert(1, preequilibration) + lines.append('\t'.join(cells)) + (root / 'measurements.tsv').write_text('\n'.join(lines) + '\n') + (root / 'conditions.tsv').write_text(conditions) + (root / 'parameters.tsv').write_text( + 'parameterId\tparameterScale\tlowerBound\tupperBound\tnominalValue\testimate\n' + 'k\tlog10\t1e-3\t1e3\t0.1\t1\n' + + ''.join(f'{pid}\tlin\t1e-3\t10\t{value}\t0\n' for pid, value in fixed) + + ''.join(f'{pid}\tlin\t1e-3\t10\t{value}\t1\n' for pid, value in estimated)) + (root / 'problem.yaml').write_text( + 'format_version: 1\nparameter_file: parameters.tsv\nproblems:\n' + ' - sbml_files: [model.xml]\n condition_files: [conditions.tsv]\n' + f' measurement_files: [measurements.tsv]\n observable_files: [{", ".join(obs_files)}]\n') + return root / 'problem.yaml' + + +def _convert_and_import(tmp_path, yaml): + from pybnf.petab import import_job + from pybnf.petab.convert import petab1to2_preserve_scale + v2_yaml = petab1to2_preserve_scale(str(yaml), str(tmp_path / 'v2')) + return import_job(v2_yaml, tmp_path / 'imported') + + +def _pybnf_objective(out, monkeypatch): + """``score(sim)``: the fitter's own objective, as built from imported.conf, on a fixed + simulation ``{species: values at _TIMES}`` (through ``evaluate_multiple``, the entry + point the fitter calls, so a prediction-dependent sigma sees the parameter values).""" + from pybnf.parse import load_config + monkeypatch.chdir(out) + cfg = load_config('imported.conf') + (model,) = cfg.exp_data + (exp_name,) = cfg.exp_data[model] + + def score(sim, parameters=None): + pset = [SimpleNamespace(name=n, value=v) for n, v in (parameters or {}).items()] + return cfg.obj.evaluate_multiple({model: {exp_name: _sim(sim)}}, cfg.exp_data, pset) + return score + + +def _sim(values): + from pybnf.data import Data + names = sorted(values) + arr = np.column_stack([_TIMES] + [values[n] for n in names]) + return Data.from_columns(arr, ['time'] + names, indvar='time') + + +def _random_simulations(n=4, seed=0): + rng = np.random.default_rng(seed) + return [{sp: rng.uniform(0.5, 8.0, size=len(_TIMES)) for sp in _DATA} for _ in range(n)] + + +def _libpetab_v1_nllh(v1_root, obs_files, sim): + """-log L of the ORIGINAL v1 problem, by libpetab's v1 calculator (independent of PyBNF).""" + import petab.v1 as pv1 + from petab.v1.calculate import calculate_llh + mdf = pv1.get_measurement_df(str(v1_root / 'measurements.tsv')) + odf = pd.concat([pv1.get_observable_df(str(v1_root / f)) for f in obs_files]) + pdf = pv1.get_parameter_df(str(v1_root / 'parameters.tsv')) + species = {oid: sp for rows in obs_files.values() for oid, sp, _, _ in rows} + sdf = mdf.rename(columns={'measurement': 'simulation'}) + sdf['simulation'] = [sim[species[o]][_TIMES.index(float(t))] + for o, t in zip(sdf['observableId'], sdf['time'])] + return -calculate_llh(mdf, sdf, odf, pdf) + + +# (case id, observable files, measurement noiseParameters, fixed parameters) +_LIBPETAB_CASES = [ + ('log10-number', {'o.tsv': [('obs_V', 'V', 'log10', '0.1')]}, {}, ()), + ('log10-number-1e-1', {'o.tsv': [('obs_V', 'V', 'log10', '1e-1')]}, {}, ()), + ('log10-number-0.10', {'o.tsv': [('obs_V', 'V', 'log10', '0.10')]}, {}, ()), + ('log10-placeholder', {'o.tsv': [('obs_V', 'V', 'log10', 'noiseParameter1_obs_V')]}, + {'obs_V': '0.1'}, ()), + ('log10-fixed-parameter', {'o.tsv': [('obs_V', 'V', 'log10', 'sigma_V')]}, {}, + (('sigma_V', 0.1),)), + ('log10-numeric-expression', {'o.tsv': [('obs_V', 'V', 'log10', '0.05 * 2')]}, {}, ()), + ('log-number', {'o.tsv': [('obs_V', 'V', 'log', '0.1')]}, {}, ()), + ('lin-number', {'o.tsv': [('obs_V', 'V', 'lin', '0.5')]}, {}, ()), + # Several observables, a mixed numeric/text noiseFormula column, two observable files. + ('two-files-mixed', {'a.tsv': [('obs_V', 'V', 'log10', '0.1')], + 'b.tsv': [('obs_W', 'W', 'log10', 'noiseParameter1_obs_W')]}, + {'obs_W': '0.2'}, ()), + ('one-table-log10-and-log', {'o.tsv': [('obs_V', 'V', 'log10', '0.1'), + ('obs_W', 'W', 'log', '0.3')]}, {}, ()), +] + + +@pytest.mark.parametrize('case,obs_files,noise_parameters,fixed', _LIBPETAB_CASES, + ids=[c[0] for c in _LIBPETAB_CASES]) +def test_objective_of_the_converted_problem_is_libpetabs_v1_likelihood( + tmp_path, monkeypatch, case, obs_files, noise_parameters, fixed): + yaml = _write_v1(tmp_path / 'v1', obs_files, noise_parameters, fixed) + score = _pybnf_objective(_convert_and_import(tmp_path, yaml), monkeypatch) + offsets = [score(sim) - _libpetab_v1_nllh(tmp_path / 'v1', obs_files, sim) + for sim in _random_simulations()] + # Data-only constants differ; anything that depends on the simulation must not. + np.testing.assert_allclose(offsets, offsets[0], rtol=0, atol=1e-9) + + +@pytest.mark.parametrize('noise,noise_parameters,estimated', [ + ('noiseParameter1_obs_V + noiseParameter2_obs_V * V', {'obs_V': '0.1;0.2'}, ()), + ('a + b * V', {}, (('a', 0.1), ('b', 0.2))), +], ids=['placeholders', 'estimated-coefficients']) +def test_prediction_dependent_log10_sigma_is_the_v1_one(tmp_path, monkeypatch, noise, + noise_parameters, estimated): + # v1: log10(y) ~ N(log10(h), sigma), sigma = 0.1 + 0.2 h, written out by hand. petab 0.10 + # would make it ln(10) * (0.1 + 0.2 h) on a natural-log residual; under the restored log10 + # column it must be the v1 sigma, including in the log(sigma) normalizer. + obs_files = {'o.tsv': [('obs_V', 'V', 'log10', noise)]} + yaml = _write_v1(tmp_path / 'v1', obs_files, noise_parameters, estimated=estimated) + score = _pybnf_objective(_convert_and_import(tmp_path, yaml), monkeypatch) + y = np.array(_DATA['V']) + offsets = [] + for sim in _random_simulations(): + h = sim['V'] + sigma = 0.1 + 0.2 * h + v1_nllh = np.sum(0.5 * np.log(2 * np.pi * sigma ** 2 * y ** 2 * math.log(10) ** 2) + + 0.5 * ((np.log10(y) - np.log10(h)) / sigma) ** 2) + offsets.append(score(sim, dict(estimated)) - v1_nllh) + np.testing.assert_allclose(offsets, offsets[0], rtol=0, atol=1e-9) + + +@pytest.mark.parametrize('rewrite', [ + 'log(10)*(0.1)', # the same factor, printed without spaces + '2.302585092994046 * (0.1)', # the factor as a number + '0.1 * log(10)', # the factor on the other side +]) +def test_an_unrecognised_petab1to2_rewrite_stops_the_conversion(tmp_path, monkeypatch, rewrite): + # End to end through petab1to2_preserve_scale, not just the helper: if a petab release + # rescales a log10 observable's sigma in any form the converter does not undo, the + # conversion must stop, naming the observable and both formulas, not carry the rescaled + # sigma into the fit. + import petab.v2.petab1to2 as p12 + real = p12.petab1to2 + + def petab1to2_that_rewrites(yaml_config, output_dir=None, **kwargs): + result = real(yaml_config, output_dir, **kwargs) + path = f'{output_dir}/o.tsv' + odf = pd.read_csv(path, sep='\t', dtype=str, keep_default_na=False) + odf.loc[odf['observableId'] == 'obs_V', 'noiseFormula'] = rewrite + odf.to_csv(path, sep='\t', index=False) + return result + + monkeypatch.setattr(p12, 'petab1to2', petab1to2_that_rewrites) + yaml = _write_v1(tmp_path / 'v1', {'o.tsv': [('obs_V', 'V', 'log10', '0.1')]}) + from pybnf.petab.convert import petab1to2_preserve_scale + with pytest.raises(PybnfError, match=rf"obs_V.*'0\.1'.*{re.escape(repr(rewrite))}"): + petab1to2_preserve_scale(str(yaml), str(tmp_path / 'v2')) + + +# Was a strict xfail on petab 0.9 with pandas < 3.1 (review finding): the v2 observables +# table was re-read with a float parse, so a 17-digit sigma no longer matched its v1 text and +# the conversion was refused. The converter now reads that table as text. +def test_a_full_precision_numeric_log10_sigma_converts_exactly(tmp_path): + # A v1 sigma written with 17 significant digits (what any float printed by Python looks + # like) must come through with no change of value and no refusal, on every petab release. + sigma = '0.0009486494471372439' + yaml = _write_v1(tmp_path / 'v1', {'o.tsv': [('obs_V', 'V', 'log10', sigma)]}) + conf = (_convert_and_import(tmp_path, yaml) / 'imported.conf').read_text() + assert f'noise_model = lognormal, sigma = fix_at {sigma}' in conf + + +@pytest.mark.xfail(_petab_minor() < (0, 10), strict=True, + reason='pre-existing, upstream: petab 0.9 petab1to2 drops a ' + 'pre-equilibration whose condition changes nothing, silently') +def test_pre_equilibration_under_an_empty_condition_survives_conversion(tmp_path): + # v1: pre-equilibrate to steady state (no changes), then simulate. The conditions table + # has a column, but c0 leaves it blank. The imported experiment must still pre-equilibrate. + yaml = _write_v1(tmp_path / 'v1', {'o.tsv': [('obs_V', 'V', 'lin', '0.5')]}, + conditions='conditionId\tW\nc0\t\nc1\t2\n', preequilibration='c0') + with warnings.catch_warnings(): + warnings.simplefilter('ignore') + conf = (_convert_and_import(tmp_path, yaml) / 'imported.conf').read_text() + (experiment,) = [ln for ln in conf.splitlines() if ln.startswith('experiment:')] + assert 'preequilibrate:' in experiment