diff --git a/docs/algorithms.rst b/docs/algorithms.rst index d42bbf8e..7287e2e2 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 3c242b05..1ff38880 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 4c7a60cf..95811df9 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 eec162c8..f85c016a 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 83604a95..2aaae5f8 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 5b499929..1fc31fba 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 b24d2286..4310e89d 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 6c36f4b3..9d9101d7 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 75978cbe..5f96af41 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 d98b0d4f..c2706440 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 f9935bb7..ba76d845 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 cf4ce1f6..a84624a9 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')