fix(objective): ave_norm_sos's normalizer averages the measured points only (#707)

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 <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 <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 <con-file>` in your fit, the constraints add extra terms to the objective function. 
 

docs/config_keys.rst

@@ -602,8 +602,9 @@ Required Keys
    - ``relative [<cv>]`` - 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 <expr>`` - an arithmetic expression over free parameters (and constants),
      evaluated per point against the current fit; the PEtab ``noiseFormula`` source.
    - ``prediction_formula <expr>`` - an expression whose ``sigma`` scales with the

docs/noise_models.rst

@@ -74,8 +74,9 @@ The sources, in the vocabulary of the :ref:`noise_model <noise_model_key>` key:
 - **relative** (``relative [<cv>]``) — 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 <expr>``) — an expression over free parameters (and,
   row by row, PEtab noise placeholders); the PEtab ``noiseFormula`` source.
 - **prediction formula** (``prediction_formula <expr>``) — an expression whose

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.

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)

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):

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

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
 
 

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()

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

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."""

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')