From a7b9c8f27ce022da56d2e0e415a655fa96fcb568 Mon Sep 17 00:00:00 2001 From: Bill Hlavacek Date: Wed, 16 Sep 2026 17:27:40 -0600 Subject: [PATCH] fix(objective): ave_norm_sos's normalizer averages the measured points only (#707) AveNormSumOfSquaresObjective.evaluate precomputed each column's normalizer as np.average(exp_data[name]) -- over the raw column, NaNs included. A NaN in an experimental column is missing data and a first-class supported value: the exp file reader parses it natively, docs/config.rst tells users to write it, SummationObjective.evaluate skips those rows, and the #479 follow-up already made every Data normalization reduce past them, because a sparse multi-observable column carries NaN wherever that observable was not measured. But this mean is a divisor. Taken over the raw column it was NaN, so ((sim - exp) / ybar)**2 returned NaN for every present point of the column too, and the whole objective went NaN. Nothing downstream recovered it. evaluate gates the prediction for NaN and inf but not the accumulated func_value, so it returned a float rather than None -- a NaN is not a failed simulation, so it never took the failed-simulation path and was handed to the optimizer as a real score. Every NaN comparison is False, so new_score < best_score was never true: every parameter set was silently rejected, and the fit ran to completion reporting no improvement, with no error and nothing in the log to say why. Trajectory.add turns the NaN into inf to keep the heap ordered, which hides it further. Four sites computed this same column mean and all four had it, so the fix is one shared pybnf.data.observed_mean that averages the non-NaN entries: - the legacy objfunc = ave_norm_sos (pybnf/objective.py); - the modern objective = ave_norm_sos, which desugars to ColumnMeanSigma, whose value() was the same np.average -- its comment claimed the result is identical to the legacy path, and it was, poisoning included. Both now read Data.column_mean, so that claim is structural rather than a parallel implementation; - the PEtab export, which writes the mean as an observable's noiseFormula constant. A sparse problem exported the literal 'nan', which is not a usable PEtab problem -- and PEtab's own round trip writes 'nan' for an unmeasured point, so sparse data is exactly what arrives here; - the PEtab import, which compares a constant sigma against the column mean to tell ave_norm_sos from sos. NaN compares equal to nothing, so a round-tripped ave_norm_sos came back silently downgraded to sos. On a dense column observed_mean is the plain mean, so every existing fit and every existing golden value is byte-identical. A column with no measured point at all has no mean and gives NaN; that is harmless, because every one of its rows is skipped in scoring and the normalizer is never read. It is written out rather than using np.nanmean so that case returns quietly instead of warning "Mean of empty slice" once per unscored column on every evaluate() of the fit. Tests: the NaN-in-experiment case was untested for this objfunc -- the only ave_norm_sos NaN test put the NaN in the simulation. Added the sparse-column value for the legacy objfunc (test_objective_funcs), for the desugared modern surface (test_objective_surface), for the PEtab export's noiseFormula (test_petab_export), and the wholly-unmeasured-column case, which also asserts no warning is raised. All four fail against the old code. Two unit tests for observed_mean itself sit in test_data_class beside the #479 nan-aware normalization test. The independent-variable column and the _SD columns get a mean computed too and never read; that is unchanged. Docs: algorithms.rst described ybar as "the average of the entire data column", which is now wrong, and the objective-function preamble did not mention that missing points are skipped at all. Both say what the code does now, as do the two column_mean noise-source descriptions in config_keys.rst and noise_models.rst. Signed-off-by: Bill Hlavacek --- docs/algorithms.rst | 4 ++-- docs/config_keys.rst | 5 ++-- docs/noise_models.rst | 5 ++-- pybnf/data.py | 36 +++++++++++++++++++++++++++++ pybnf/noise/source.py | 5 ++-- pybnf/objective.py | 5 +++- pybnf/petab/export.py | 7 ++++-- pybnf/petab/import_.py | 7 +++++- tests/test_data_class.py | 30 ++++++++++++++++++++++++ tests/test_objective_funcs.py | 38 ++++++++++++++++++++++++++++++ tests/test_objective_surface.py | 18 +++++++++++++++ tests/test_petab_export.py | 41 +++++++++++++++++++++++++++++++++ 12 files changed, 189 insertions(+), 12 deletions(-) diff --git a/docs/algorithms.rst b/docs/algorithms.rst index d42bbf8e2..7287e2e22 100644 --- a/docs/algorithms.rst +++ b/docs/algorithms.rst @@ -69,13 +69,13 @@ The ``latin_hypercube`` option for initialization is enabled by default. This op Objective functions ^^^^^^^^^^^^^^^^^^^ -All algorithms use an objective function to evaluate the quality of fit for each parameter set. The objective function is set with the ``objfunc`` key. The following options are available. Note that :math:`y_i` are the experimental data points and :math:`a_i` are the simulated data points. The summation is over all experimental data points. +All algorithms use an objective function to evaluate the quality of fit for each parameter set. The objective function is set with the ``objfunc`` key. The following options are available. Note that :math:`y_i` are the experimental data points and :math:`a_i` are the simulated data points. The summation is over all experimental data points; a point recorded as ``nan`` is missing data, not an observation, and is skipped (see the :ref:`exp file `). * Chi squared (``obj_func = chi_sq``): :math:`f(y, a) = \sum_i \frac{(y_i - a_i)^2}{2 \sigma_i^2}` , where :math:`\sigma_i` is the standard deviation of point :math:`y_i`, which must be specified in the :ref:`exp file `. * Sum of squares (``obj_func = sos``): :math:`f(y, a) = \sum_i (y_i - a_i)^2` * Sum of differences (``obj_func = sod``): :math:`f(y, a) = \sum_i |y_i - a_i|` * Normalized sum of squares (``obj_func = norm_sos``): :math:`f(y, a) = \sum_i \frac{(y_i - a_i)^2}{y_i^2}` - * Average-normalized sum of squares (``obj_func = ave_norm_sos``): :math:`f(y, a) = \sum_i \frac{(y_i - a_i)^2}{\bar{y}^2}`, where :math:`\bar{y}` is the average of the entire data column :math:`y`. + * Average-normalized sum of squares (``obj_func = ave_norm_sos``): :math:`f(y, a) = \sum_i \frac{(y_i - a_i)^2}{\bar{y}^2}`, where :math:`\bar{y}` is the average of the measured values in data column :math:`y` -- the ``nan`` (missing) entries of a sparse column are left out of it, as they are of the summation. If you include any :ref:`constraints ` in your fit, the constraints add extra terms to the objective function. diff --git a/docs/config_keys.rst b/docs/config_keys.rst index 3c242b055..1ff38880b 100644 --- a/docs/config_keys.rst +++ b/docs/config_keys.rst @@ -602,8 +602,9 @@ Required Keys - ``relative []`` - constant coefficient of variation: ``sigma = cv * |value|``, so the noise scales with the measurement (``cv`` defaults to 1). This is the heteroscedastic model the legacy ``norm_sos`` fits. - - ``column_mean`` - ``sigma`` is the observable's experimental column mean (one - scale per column). This is the model the legacy ``ave_norm_sos`` fits. + - ``column_mean`` - ``sigma`` is the mean of the observable's measured experimental + values (one scale per column; ``nan`` entries are missing data and do not enter the + mean). This is the model the legacy ``ave_norm_sos`` fits. - ``formula `` - an arithmetic expression over free parameters (and constants), evaluated per point against the current fit; the PEtab ``noiseFormula`` source. - ``prediction_formula `` - an expression whose ``sigma`` scales with the diff --git a/docs/noise_models.rst b/docs/noise_models.rst index 4c7a60cff..95811df9a 100644 --- a/docs/noise_models.rst +++ b/docs/noise_models.rst @@ -74,8 +74,9 @@ The sources, in the vocabulary of the :ref:`noise_model ` key: - **relative** (``relative []``) — a constant coefficient of variation, :math:`\sigma = \mathrm{cv}\cdot|\mathrm{value}|`; the heteroscedastic model the legacy ``norm_sos`` fits. -- **column mean** (``column_mean``) — one scale per column, the observable's - experimental column mean; the model the legacy ``ave_norm_sos`` fits. +- **column mean** (``column_mean``) — one scale per column, the mean of the + observable's measured experimental values (``nan`` entries are missing data and + do not enter it); the model the legacy ``ave_norm_sos`` fits. - **formula** (``formula ``) — an expression over free parameters (and, row by row, PEtab noise placeholders); the PEtab ``noiseFormula`` source. - **prediction formula** (``prediction_formula ``) — an expression whose diff --git a/pybnf/data.py b/pybnf/data.py index eec162c8d..f85c016aa 100644 --- a/pybnf/data.py +++ b/pybnf/data.py @@ -57,6 +57,31 @@ def slice_for(self, selector, axis='parameter'): return tensor[:, col, :] +def observed_mean(values): + """The mean of a data column's *observed* values -- the column's scale, which + ``ave_norm_sos`` and the ``column_mean`` sigma source normalize by. + + NaN is missing data, not an observation: scoring skips those rows + (``SummationObjective.evaluate``) and :class:`Data`'s normalizations reduce past + them (#479). A plain ``np.average`` over the raw column returns NaN for any sparse + multi-observable column, and because this mean is a *divisor*, that NaN then + poisons every **present** point of the column too -- the whole objective goes NaN, + which is not a failed simulation, so it reaches the optimizer as a real score and + silently rejects every parameter set (#707). + + Values with no observed entry have no mean and give NaN. Harmless in scoring: none + of such a column's rows are scored, so the normalizer is never read. Written out + rather than ``np.nanmean`` so that case returns quietly instead of warning + ("Mean of empty slice") on every evaluate() of the fit. + + :param values: 1D array of column values, possibly containing NaN + :return: Mean over the non-NaN entries, or NaN if there are none + """ + values = np.asarray(values, dtype=float) + observed = values[~np.isnan(values)] + return np.mean(observed) if observed.size else np.nan + + def stack_scan_sensitivities(per_point): """Stack per-dose-point forward-sensitivity tensors into one scan :class:`OutputSensitivities`. @@ -402,6 +427,17 @@ def __setitem__(self, key, value): idx = self.cols[key] self.data[:, idx] = value + def column_mean(self, col_header): + """ + The mean of a column's observed (non-NaN) values -- see :func:`observed_mean` + for why the NaNs must be excluded (#707). + + :param col_header: Data column name + :type col_header: str + :return: Mean over the non-NaN entries, or NaN if there are none + """ + return observed_mean(self[col_header]) + def get_row(self, col_header, value): """ Returns the (first) data row in which field col_header is equal to value. diff --git a/pybnf/noise/source.py b/pybnf/noise/source.py index 83604a955..2aaae5f8b 100644 --- a/pybnf/noise/source.py +++ b/pybnf/noise/source.py @@ -644,5 +644,6 @@ def value(self, owner, sim_data, sim_row, exp_data, exp_row, col_name): # The mean is a per-column constant; recomputing it per point is O(1) amortized # over the small data PyBNF fits and keeps the source stateless (no per-exp_data # cache to invalidate across models/suffixes). The legacy ave_norm_sos - # precomputes it once in evaluate(); the result is identical. - return np.average(exp_data[col_name]) + # precomputes it once in evaluate(); the result is identical -- structurally so, + # since both paths go through Data.column_mean (which is nan-aware, #707). + return exp_data.column_mean(col_name) diff --git a/pybnf/objective.py b/pybnf/objective.py index 5b499929c..1fc31fbaf 100644 --- a/pybnf/objective.py +++ b/pybnf/objective.py @@ -2590,7 +2590,10 @@ class AveNormSumOfSquaresObjective(SummationObjective): def evaluate(self, sim_data, exp_data, show_warnings=True, data_key=None): # Precalculate the average of each exp column to use for all points in this call. - self.aves = {name: np.average(exp_data[name]) for name in exp_data.cols} + # Over the OBSERVED values only (#707): the mean is a divisor, so a NaN in a sparse + # multi-observable column would otherwise poison every present point of that column + # too -- see Data.column_mean. + self.aves = {name: exp_data.column_mean(name) for name in exp_data.cols} return super().evaluate(sim_data, exp_data, show_warnings, data_key=data_key) def eval_point(self, sim_data, exp_data, sim_row, exp_row, col_name): diff --git a/pybnf/petab/export.py b/pybnf/petab/export.py index b24d22864..4310e89d9 100644 --- a/pybnf/petab/export.py +++ b/pybnf/petab/export.py @@ -69,7 +69,7 @@ import numpy as np from .. import edition -from ..data import Data +from ..data import Data, observed_mean from ..objective import _OBJECTIVE_DESUGAR from ..parse import ploop from ..printing import PybnfError @@ -1202,7 +1202,10 @@ def _noise_source_for_column(verb, arg, col, datas): if verb == 'fix_at': return ('constant', float(arg)) if verb == 'column_mean': - return ('constant', float(np.average(np.concatenate([d[col] for d in holders])))) + # Over the OBSERVED values only (#707), matching what ColumnMeanSigma computes at + # fit time: a sparse multi-observable column carries NaN in its unmeasured rows, and + # a plain average would export 'nan' as this observable's noiseFormula constant. + return ('constant', float(observed_mean(np.concatenate([d[col] for d in holders])))) if verb == 'fit': # A free-parameter (estimated) sigma -> a bare-id noiseFormula naming the noise # parameter (declared estimated in parameters.tsv; admitted as an observation-layer diff --git a/pybnf/petab/import_.py b/pybnf/petab/import_.py index 6c36f4b3f..9d9101d73 100644 --- a/pybnf/petab/import_.py +++ b/pybnf/petab/import_.py @@ -82,6 +82,7 @@ import numpy as np +from ..data import observed_mean from ..printing import PybnfError from ..priors import PRIOR_KEYWORD_MAP from .conditions import ( @@ -713,7 +714,11 @@ def column_mean_of(observable_id): col = observable_id_to_column[observable_id] values = [data[col] for group in datas.values() for data in group if col in data.cols] - return float(np.average(np.concatenate(values))) + # Observed values only (#707) -- the same mean the export wrote. A plain average + # over a sparse column is NaN, which compares equal to nothing, so the sigma + # constant would fail to match and a round-tripped ave_norm_sos would silently + # come back as sos. + return float(observed_mean(np.concatenate(values))) return column_mean_of diff --git a/tests/test_data_class.py b/tests/test_data_class.py index 75978cbe2..5f96af412 100644 --- a/tests/test_data_class.py +++ b/tests/test_data_class.py @@ -1,4 +1,6 @@ import math +import warnings + import numpy as np import numpy.testing as npt import pytest @@ -282,6 +284,34 @@ def mk(): npt.assert_allclose(du.data[[0, 2], 1], np.array([0.0, 1.0])) assert np.isnan(du.data[1, 1]) + def test_column_mean_is_nan_aware(self): + # #707, the same sparse-column hazard one step later: column_mean is the DIVISOR + # ave_norm_sos / the column_mean sigma source normalize by, so a plain np.average + # returning NaN poisons every *present* point of the column too. Average the + # measured points only. obs1 = [3, NaN, 4] -> 3.5, not NaN. + d = data.Data() + d.data = d._read_file_lines(['# x obs1 obs2 obs3\n', + ' 0 3 4 5\n', + ' 1 nan 3 6\n', + ' 2 4 nan 10\n'], r'\s+') + assert d.column_mean('obs1') == 3.5 + assert d.column_mean('obs2') == 3.5 + npt.assert_allclose(d.column_mean('obs3'), 7.0) # dense column: plain mean + npt.assert_allclose(d.column_mean('x'), 1.0) + + def test_column_mean_of_an_unmeasured_column_is_nan_and_quiet(self): + # A column with no observed value has no mean. It returns NaN rather than raising, + # and must do so WITHOUT a warning -- np.nanmean would emit "Mean of empty slice" + # on every evaluate() call, once per unscored column, for the whole fit. The NaN is + # harmless: every row of such a column is skipped in scoring, so it is never read. + d = data.Data() + d.data = d._read_file_lines(['# x obs1 obs2\n', + ' 0 1 nan\n', + ' 1 2 nan\n'], r'\s+') + with warnings.catch_warnings(): + warnings.simplefilter('error') + assert np.isnan(d.column_mean('obs2')) + def test_zero_normalization(self): d0 = copy.deepcopy(self.d0) d0.normalize_to_zero() diff --git a/tests/test_objective_funcs.py b/tests/test_objective_funcs.py index d98b0d4fa..c2706440a 100644 --- a/tests/test_objective_funcs.py +++ b/tests/test_objective_funcs.py @@ -1,3 +1,5 @@ +import warnings + from .context import data, noise, objective, printing, raises import numpy as np import numpy.testing as npt @@ -152,6 +154,42 @@ def test_obj_nan(self): assert self.norm_sos.evaluate(self.d1s_nan, self.d1e) is None assert self.ave_norm_sos.evaluate(self.d1s_nan, self.d1e) is None + def test_ave_norm_sos_normalizer_skips_missing_experimental_points(self): + """#707: the ``ybar`` normalizer is averaged over the OBSERVED points only. + + A NaN in an experimental column is missing data -- ``evaluate`` already skips + those rows. But ``ybar`` is a divisor, so averaging it over the raw column made + it NaN and turned every *present* point of that column NaN too, and hence the + whole objective. That NaN is not a failed simulation (``evaluate`` returns a + float, not None), so it reached the optimizer as a real score; since every NaN + comparison is False, no parameter set ever improved on another and the fit + silently reported nothing. Sparse multi-observable exp files -- one observable + per measurement time -- are the ordinary way to write this data. + + obs1 = [1, NaN, 3] -> ybar = 2, and only rows 0 and 2 are scored. + """ + exp = _mkdata(['# x obs1\n', ' 0 1\n', ' 1 nan\n', ' 2 3\n']) + sim = _mkdata(['# x obs1\n', ' 0 1.1\n', ' 1 2.0\n', ' 2 3.1\n']) + obj = objective.AveNormSumOfSquaresObjective() + npt.assert_almost_equal(obj.evaluate(sim, exp), + ((1.1 - 1) / 2.) ** 2 + ((3.1 - 3) / 2.) ** 2) + assert obj.aves['obs1'] == 2.0 + + def test_ave_norm_sos_tolerates_a_wholly_unmeasured_column(self): + """An exp column with no observed point has no mean. Its NaN normalizer is never + read (every one of its rows is skipped), so the scored columns are unaffected -- + and computing it must not warn, since it happens on every evaluate() of the fit.""" + exp = _mkdata(['# x obs1 obs3\n', + ' 0 1 nan\n', ' 1 nan nan\n', ' 2 3 nan\n']) + sim = _mkdata(['# x obs1 obs3\n', + ' 0 1.1 9\n', ' 1 2.0 9\n', ' 2 3.1 9\n']) + obj = objective.AveNormSumOfSquaresObjective() + with warnings.catch_warnings(): + warnings.simplefilter('error') + value = obj.evaluate(sim, exp) + npt.assert_almost_equal(value, ((1.1 - 1) / 2.) ** 2 + ((3.1 - 3) / 2.) ** 2) + assert np.isnan(obj.aves['obs3']) + def test_obj_inf(self): assert self.chi_sq.evaluate(self.d1s_inf, self.d1e_sd) is None assert self.sos.evaluate(self.d1s_inf, self.d1e) is None diff --git a/tests/test_objective_surface.py b/tests/test_objective_surface.py index f9935bb74..ba76d845d 100644 --- a/tests/test_objective_surface.py +++ b/tests/test_objective_surface.py @@ -88,6 +88,24 @@ def test_objective_ave_norm_sos_uses_column_mean_source(): assert isinstance(src, noise.ColumnMeanSigma) +def test_column_mean_source_skips_missing_experimental_points(): + """#707: the ``column_mean`` sigma is the mean of the OBSERVED points. + + The desugared path shares the legacy defect because it shares the quantity -- a + sigma averaged over the raw column is NaN on any sparse multi-observable column, + and a NaN sigma makes every present point's density NaN. Both now read + ``Data.column_mean``, so the "the result is identical" claim in ``ColumnMeanSigma`` + holds structurally. obs1 = [1, NaN, 3] -> sigma = 2; value is the legacy score + times the proper 1/2. + """ + exp = _mkdata(['# x obs1\n', ' 0 1\n', ' 1 nan\n', ' 2 3\n']) + sim = _mkdata(['# x obs1\n', ' 0 1.1\n', ' 1 2.0\n', ' 2 3.1\n']) + legacy = objective.AveNormSumOfSquaresObjective().evaluate(sim, exp) + modern = _modern({'objective': 'ave_norm_sos'}).evaluate(sim, exp) + assert modern == pytest.approx(0.5 * legacy) + assert modern == pytest.approx(0.5 * (((1.1 - 1) / 2.) ** 2 + ((3.1 - 3) / 2.) ** 2)) + + def test_chi_sq_desugar_is_value_identical_to_legacy(): """chi_sq already carried the 1/2, so the desugared form is value-identical, not merely argmin-identical.""" diff --git a/tests/test_petab_export.py b/tests/test_petab_export.py index cf4ce1f69..a84624a97 100644 --- a/tests/test_petab_export.py +++ b/tests/test_petab_export.py @@ -407,6 +407,47 @@ def test_ave_norm_sos_is_the_column_mean(self, tmp_path_factory): assert r['noisePlaceholders'] == '' assert _petab_validation_errors(out / 'problem.yaml') == [] + def test_column_mean_noise_formula_skips_unmeasured_points(self, tmp_path_factory): + """#707: the exported column-mean sigma averages the OBSERVED points only. + + PEtab's own round-trip writes 'nan' for a point an observable does not have + (see the exporter's measurement rows), so a sparse multi-observable data file is + exactly what comes back through this path. A plain average over the raw column + put 'nan' in noiseFormula -- an unusable PEtab problem, and on reimport a sigma + that matches no objective.""" + import shutil + src = tmp_path_factory.mktemp('colmean_sparse_src') + shutil.copy(DEMO_DIR / DEMO_MODEL, src / DEMO_MODEL) + # par1.exp with x unmeasured on the odd rows -- y keeps its full time course. + rows = [r for r in (DEMO_DIR / 'par1.exp').read_text().splitlines() if r.strip()] + sparse = [rows[0]] + for i, row in enumerate(rows[1:]): + fields = row.split() + if i % 2: + fields[1] = 'nan' # x unmeasured at this time point + sparse.append('\t'.join(fields)) + (src / 'par1.exp').write_text('\n'.join(sparse) + '\n') + (src / 'job.conf').write_text( + f'edition = 2\njob_type = de\nobjective = ave_norm_sos\n' + f'model: {DEMO_MODEL}\n' + 'experiment: par1, data: par1.exp\n' + 'uniform_var = v1 0 10\nuniform_var = v2 0 10\n' + 'uniform_var = v3 0 10\n') + out = tmp_path_factory.mktemp('colmean_sparse_out') + export_job(src / 'job.conf', out) + + data = Data(file_name=str(src / 'par1.exp')) + by_id = {r['observableId']: r for r in _tsv_rows(out / 'observables.tsv')} + for oid, col in (('obs_x', 'x'), ('func_y', 'y')): + column = data[col] + expected = np.average(column[~np.isnan(column)]) + written = float(by_id[oid]['noiseFormula']) + assert not np.isnan(written) # the regression: 'nan' was written here + assert written == pytest.approx(expected) + # The sparse column really is sparse -- otherwise this test proves nothing. + assert np.isnan(data['x']).any() + assert _petab_validation_errors(out / 'problem.yaml') == [] + def test_chi_sq_still_uses_the_sd_placeholder(self, tmp_path_factory): # The existing _SD path is unchanged: per-point placeholder + noiseParameters. out = self._export(tmp_path_factory, 'chi_sq')