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..36793701 --- /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, so the noise formula is the same 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..6d70c8df 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,23 @@ 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') + # 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) return v2_yaml @@ -499,6 +517,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. 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 + + +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..a3715d88 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(): + 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) +_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', [ 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