From 21a4b0bc33bfd7e80392b65346f1a8d0d8691073 Mon Sep 17 00:00:00 2001 From: Bill Hlavacek Date: Wed, 16 Sep 2026 18:26:41 -0600 Subject: [PATCH] fix(data): every normalization reduces over a column's measured rows (#726) A sparse column now normalizes to exactly what the dense column of its measured values would: each reduction (max, min, mean, std) and each reference row (argmax, argmin, the init/unit baseline) skips the NaNs, so a NaN row stays a NaN row and moves nothing else. #479 established that for peak and floor; zero, init and unit were missed, and neither issue covered a column with no measured value at all. A simulated column can legitimately carry NaN at some output rows -- one failed integration step, or an observable that is 0/0 at t=0 -- while the exp file measures a handful of the simulated times, so that row is frequently one nothing scores. normalize_to_zero reduces the whole column, so a single NaN made the mean NaN, hence every centered value NaN, hence the std NaN. _subtract_baseline read row 0, so a NaN there did the same to init and unit. The column then scored as a failed simulation and the parameter set was thrown away, with nothing logged to say that a normalizer rather than the model had made the call. Under normalization = peak the identical simulation scored a real number, which is what showed this was a defect and not a policy. init and unit now take the first MEASURED row as their baseline, which is what "the initial value" meant; _subtract_baseline returns the row it used so normalize_to_unit_scale can record the real baseline_row, since the gradient reads the baseline back by index. A column with no measured value at all has no peak, no min, no mean and no baseline. np.nanargmax / np.nanargmin raise ValueError('All-NaN slice encountered') rather than returning a sentinel, so peak, floor and unit crashed outright and zero and init reached the same state by poisoning. For that case #479's nan-awareness had converted silent corruption into an unexplained traceback. Such a column is now left untouched, with an identity record. Normalization is a transform; whether an all-NaN simulated column is a failure is scoring's call, and scoring already makes it. Both fit scoring paths absorbed the exception in #388's blanket handler and penalized the evaluation with inf, so a fit never died, but two callers do not have that handler: - model_check.run_check calls normalize outside the try below it (which covers only postprocess_data), and neither it nor its caller in pybnf.py catches anything. So --check-simulation on a model producing an all-NaN observable died with a numpy traceback, bypassing the 'Simulation contained NaN or Inf values' diagnostic three lines later that the command exists to print. - the floor is applied to the EXPERIMENTAL data at config load, so an exp file carrying a wholly unmeasured observable column -- the shape #707 supports -- crashed before the fit began. Neither caller grows a handler; the transform stops raising instead. gradient/assembly.py's z-score rule changes with it. Its two reductions (s_bar and the dsigma sum) ran over every row and divided by nrows - ddof. Once a partially-NaN column scores it also reaches the gradient, where those full-column reductions return NaN and poison the whole Jacobian, so they now mask to the measured rows -- read off ``normed``, which marks them -- and count the same points the value-side std did. That keeps record.ddof at its plain user-facing meaning rather than smuggling the NaN count into it. peak, init and unit need no gradient change: their rules read two or three specific rows rather than reducing the column, and nanargmax never returns a NaN row. Tests. test_data_class gains the invariant as a property over random draws -- a sparse column against the dense column of its measured values, per method -- which is a stronger oracle than the hand-computed cases beside it and states the rule in one line; a hand-computed sparse case pinning the recorded scale, ref_row and baseline_row; and the untouched all-NaN column. test_gradient_assembly extends its finite-difference oracle to a column with a NaN at an unscored row, on the baseline row and mid-column. test_model_check pins that an all-NaN column reaches the NaN message rather than a traceback, driving the real Result.normalize rather than the file's recording fake. Against the old code 3 of the 5 property cases fail (peak and floor pass, the controls), 4 of the 8 finite-difference cases, and all 3 model-check cases. On a dense column every reduction is the plain one it always was, so existing fits are byte-identical. Two degenerate cases are unchanged and stay as they were: a z-score over a single measured value has an undefined std, and unit on a constant column divides by 0. Both behave exactly as the equivalent dense column always has -- verified against the old code -- so this change makes them reachable from sparse data rather than introducing them. Signed-off-by: Bill Hlavacek --- docs/config_keys.rst | 9 +++ pybnf/data.py | 120 +++++++++++++++++++++++++++++--- pybnf/gradient/assembly.py | 18 +++-- tests/test_data_class.py | 99 ++++++++++++++++++++++++++ tests/test_gradient_assembly.py | 57 +++++++++++++++ tests/test_model_check.py | 41 ++++++++++- 6 files changed, 328 insertions(+), 16 deletions(-) diff --git a/docs/config_keys.rst b/docs/config_keys.rst index 1ff38880b..143b9e004 100644 --- a/docs/config_keys.rst +++ b/docs/config_keys.rst @@ -1405,6 +1405,15 @@ Algorithm Options ``peak`` / ``init`` / ``zero`` / ``unit`` rescale the **simulated** column only (the data is assumed pre-normalized by the user); ``floor`` and ``scale`` are applied symmetrically to the model and the data (a floor or an analytic scale is only meaningful applied to both). + + Every transform reads only a column's **measured** rows. A ``nan`` is missing data, not a + value: it stays ``nan`` and contributes to no maximum, minimum, mean or standard deviation, + and ``init`` / ``unit`` take their baseline from the first row that actually holds a value + rather than from row 0. So a simulation that produces ``nan`` at one output time -- a failed + integration step, or an observable that is 0/0 at ``t = 0`` -- keeps the rest of its column, + and is discarded only if a ``nan`` falls on a time the data actually measures. A column with + no measured value at all is left untouched; that simulation is then rejected by scoring, as + it would be without any normalization. Together with ``objective = lognormal`` the chain ``floor 0.03, scale`` spells the "sum of squared log-differences of geometric-mean-normalized trajectories" objective common to arbitrary-unit fluorescence / blot fits. diff --git a/pybnf/data.py b/pybnf/data.py index f85c016aa..e0a6973e7 100644 --- a/pybnf/data.py +++ b/pybnf/data.py @@ -82,6 +82,49 @@ def observed_mean(values): return np.mean(observed) if observed.size else np.nan +def has_no_observed_value(column): + """True when a column holds no measured entry at all -- every row NaN. + + Such a column has no peak, no minimum, no mean and no baseline, so there is nothing + for a normalization to read: ``np.nanargmax`` / ``np.nanargmin`` raise + ``ValueError: All-NaN slice encountered`` rather than returning a sentinel (#726). + The normalizations treat it as a no-op instead of raising. Normalization is a + transform; whether an all-NaN simulated column is a *failure* is scoring's call, and + scoring already makes it -- ``SummationObjective.evaluate`` returns ``None`` for a NaN + prediction, the failed-simulation path. Raising here instead pre-empts that decision + in the two callers that do not wrap ``normalize`` in #388's handler. + """ + return not np.any(~np.isnan(np.asarray(column, dtype=float))) + + +def _no_op_record(method, ddof=0): + """The :class:`NormalizationRecord` for a transform that was skipped because its column + holds no measured value (#726). + + The column is left exactly as it stands, so this records the identity: divide by 1, add + ``rho = 0``, read row 0. Recording rather than skipping keeps "every transform applied + leaves one record" true, which is the invariant the chain folding relies on (#539, + ADR-0102). The gradient never reads it -- an all-NaN column cannot score, so no gradient + is taken through it -- but a chain that puts a real transform after this one still finds + a well-formed stage beneath it rather than an empty chain. + """ + return NormalizationRecord(method, 1.0, ref_row=0, + baseline_row=0 if method == 'unit' else None, + rho=0.0, ddof=ddof) + + +def first_observed_row(column): + """The index of a column's first measured entry, or ``None`` if it has none. + + The baseline row of ``init`` / ``unit``. Row 0 is the intended baseline only because + it is normally the first measured row; when it is NaN -- an observable that is 0/0 at + its initial condition, say -- subtracting or dividing by it turns the whole column NaN + (#726). The first row that actually holds a value is what "initial value" meant. + """ + observed = np.flatnonzero(~np.isnan(np.asarray(column, dtype=float))) + return int(observed[0]) if observed.size else None + + def stack_scan_sensitivities(per_point): """Stack per-dose-point forward-sensitivity tensors into one scan :class:`OutputSensitivities`. @@ -568,6 +611,12 @@ def normalize_to_peak(self, idx=0, cols='all'): cols.remove(idx) for c in cols: column = self.data[:, c] + # A column with no measured point has no peak to divide by (#726): leave it + # as it is -- already all-NaN, so any divisor would leave it all-NaN anyway -- + # and let scoring decide it is a failed simulation. np.nanargmax would raise. + if has_no_observed_value(column): + self._record_normalization(c, _no_op_record('peak'), column) + continue # Record N = peak and its row before the in-place divide overwrites them # (#453): the gradient threads d(raw/N)/d theta. Additive, value-preserving. # nan-aware (#479 follow-up): a sparse multi-observable column carries NaN in the @@ -593,11 +642,21 @@ def normalize_to_init(self, idx=0, cols='all'): cols = list(range(self.data.shape[1])) cols.remove(idx) for c in cols: - # Record N = initial value before the in-place divide overwrites row 0 - # (#453): ref_row 0 is the divisor's source row for the gradient chain rule. + column = self.data[:, c] + # The divisor is the first MEASURED value, not row 0 (#726). Row 0 is the + # intended baseline only because it is normally the first measured row; when + # it is NaN (an observable that is 0/0 at its initial condition, say) dividing + # by it turns every row NaN, and the column is then scored as a failed + # simulation even though the NaN row may be one no exp point ever reads. + base = first_observed_row(column) + if base is None: + self._record_normalization(c, _no_op_record('init'), column) + continue + # Record N = initial value before the in-place divide overwrites its row + # (#453): ref_row is the divisor's source row for the gradient chain rule. self._record_normalization(c, NormalizationRecord( - 'init', float(self.data[0, c]), ref_row=0), self.data[:, c]) - self.data[:, c] = self.data[:, c] / self.data[0, c] + 'init', float(self.data[base, c]), ref_row=base), column) + self.data[:, c] = self.data[:, c] / self.data[base, c] def normalize_to_zero(self, idx=0, bc=True, cols='all'): """ @@ -622,22 +681,46 @@ def normalize_to_zero(self, idx=0, bc=True, cols='all'): # previous stage, whose rule reads them -- ADR-0102). The arithmetic is the # same subtract-then-divide as the in-place form it replaces. col = self.data[:, c] - centered = col - np.mean(col) - std = np.std(centered, ddof=ddof) + if has_no_observed_value(col): + self._record_normalization(c, _no_op_record('zero', ddof=ddof), col) + continue + # Over the MEASURED points only (#726). The z-score is the one method that + # reduces the whole column, so a single NaN row -- one failed integration step, + # or a 0/0 observable at t=0 -- made the mean NaN, hence every centered value + # NaN, hence the std NaN: the entire column was destroyed and the parameter set + # discarded as a failed simulation, even when no exp point reads that row. + # nanstd == std on a dense column, so a NaN-free fit is byte-identical. + centered = col - observed_mean(col) + std = np.nanstd(centered, ddof=ddof) # Record the z-score scale (std; 0 means the column was left un-divided) and # the K - ddof denominator the gradient's d std/d theta uses (#453). z-score # couples every row, so only these scalars are recorded -- the per-row mean of - # the sensitivities is recomputed from the tensor at gradient time. + # the sensitivities is recomputed from the tensor at gradient time, over the + # same measured rows this std used (gradient/assembly.py masks on ``normed``). self._record_normalization(c, NormalizationRecord('zero', float(std), ddof=ddof), col) self.data[:, c] = centered / std if std != 0 else centered def _subtract_baseline(self, idx=0, cols='all'): + """Shift each column so its baseline sits at 0, and return ``{col_index: baseline_row}``. + + The baseline is the first **measured** row rather than row 0 (#726): row 0 is the + intended baseline only because it is normally the first measured row, and subtracting + a NaN row 0 turns the whole column NaN. The row is returned because the caller records + it as ``NormalizationRecord.baseline_row``, which the gradient reads back + (``tensor_sens(col_name, record.baseline_row)``), so the recorded row has to be the one + actually subtracted. A column with no measured row is left untouched and maps to + ``None``. + """ if cols == 'all': cols = list(range(self.data.shape[1])) cols.remove(idx) + baselines = {} for c in cols: col = self.data[:, c] - self.data[:, c] = col - self.data[0, c] + baselines[c] = base = first_observed_row(col) + if base is not None: + self.data[:, c] = col - self.data[base, c] + return baselines def normalize_to_unit_scale(self, idx=0, cols='all'): """ @@ -659,8 +742,15 @@ def normalize_to_unit_scale(self, idx=0, cols='all'): # this transform consumed, not the baseline-subtracted ones (ADR-0102). A no-op (None) # unless this column is already normalized, so a single unit-scaling copies nothing. consumed = {c: self._chain_stage_input(c) for c in cols} - self._subtract_baseline(idx, cols) + # The baseline is the first measured row, and which row that was is recorded below -- + # the gradient reads the baseline back by row index (#726). + baselines = self._subtract_baseline(idx, cols) for c in cols: + # A column with no measured point has no baseline and no max (#726): nothing was + # subtracted, nothing is divided, and np.nanargmax/argmin would raise on it. + if baselines[c] is None: + self._record_normalization(c, _no_op_record('unit'), consumed[c]) + continue # nan-aware (#479 follow-up): skip NaN rows (a sparse column's unmeasured points) so a # multi-observable target is not poisoned by np.max/np.min seeing a NaN. cmax = np.nanmax(self.data[:, c]) @@ -670,14 +760,15 @@ def normalize_to_unit_scale(self, idx=0, cols='all'): # sensitivity enters with a flipped sign (#453); the baseline is still row 0. self._record_normalization(c, NormalizationRecord( 'unit', float(np.abs(np.nanmin(self.data[:, c]))), - ref_row=int(np.nanargmin(self.data[:, c])), baseline_row=0, sign=-1.0), + ref_row=int(np.nanargmin(self.data[:, c])), + baseline_row=baselines[c], sign=-1.0), consumed[c]) self.data[:, c] = self.data[:, c] / np.abs(np.nanmin(self.data[:, c])) else: # N = the max-after-baseline; ref_row is its argmax, baseline is row 0 (#453). self._record_normalization(c, NormalizationRecord( 'unit', float(cmax), ref_row=int(np.nanargmax(self.data[:, c])), - baseline_row=0, sign=1.0), consumed[c]) + baseline_row=baselines[c], sign=1.0), consumed[c]) self.data[:, c] = self.data[:, c] / np.nanmax(self.data[:, c]) def normalize_to_floor(self, rho, idx=0, cols='all'): @@ -709,6 +800,13 @@ def normalize_to_floor(self, rho, idx=0, cols='all'): # whole column (every point -> NaN -> silently skipped in scoring -> objective 0.0), # so take the max/argmax over the measured (non-NaN) points only. On a dense column # (no NaN) nanmax == max, so this is byte-identical for the common case. + # No measured point means no max to take a fraction of (#726). The floor is the one + # transform applied to the EXPERIMENTAL data too (config.py, ADR-0066), so this is + # reached at config load by an exp file carrying a wholly unmeasured observable + # column -- a shape #707 supports. np.nanargmax would raise before the fit starts. + if has_no_observed_value(column): + self._record_normalization(c, _no_op_record('floor'), column) + continue cmax = float(np.nanmax(column)) # Record the added amount (rho) and the max its argmax row before the offset -- the # gradient's ∂(x + rho*max)/∂θ = s_i + rho*s_argmax reads them (#533). diff --git a/pybnf/gradient/assembly.py b/pybnf/gradient/assembly.py index 3fe9cedc8..d69258b50 100644 --- a/pybnf/gradient/assembly.py +++ b/pybnf/gradient/assembly.py @@ -914,13 +914,23 @@ def _zscore_sensitivity(record, col_name, row, tensor_sens, normed, s_i): (``n_k = (raw_k - μ)/σ`` is the recorded normalized value, so ``(raw_k - μ) = σ·n_k`` cancels the σ in ``∂σ/∂θ``). A σ of 0 means ``Data`` left the column un-divided - (``n_i = raw_i - μ``), so ``∂n_i/∂θ = s_i - s_bar``.""" - nrows = len(normed) - all_s = np.array([tensor_sens(col_name, k) for k in range(nrows)]) # (K, n_param) + (``n_i = raw_i - μ``), so ``∂n_i/∂θ = s_i - s_bar``. + + Both reductions run over the column's **measured** rows -- the same points + ``normalize_to_zero`` centred and scaled over (#726) -- so ``K`` is the measured count and + ``s_bar`` their mean. A NaN row carries no value for ``n_k`` to weight, and including one + would make ``s_bar``, ``∂σ/∂θ`` and hence the whole gradient NaN. A simulated column can + reach here with a NaN in a row no exp point scores, which is precisely the case #726 stopped + poisoning the rest of the column; ``normed`` marks those rows, so the mask is read off it. + On a dense column every row is measured and this is the plain reduction it always was.""" + normed = np.asarray(normed, dtype=float) + measured = np.flatnonzero(~np.isnan(normed)) + all_s = np.array([tensor_sens(col_name, k) for k in measured]) # (K, n_param) s_bar = all_s.mean(axis=0) if record.scale == 0.0: return s_i - s_bar - dsigma = (normed[:, np.newaxis] * (all_s - s_bar)).sum(axis=0) / (nrows - record.ddof) + dsigma = ((normed[measured, np.newaxis] * (all_s - s_bar)).sum(axis=0) + / (len(measured) - record.ddof)) return (s_i - s_bar) / record.scale - normed[row] * dsigma / record.scale diff --git a/tests/test_data_class.py b/tests/test_data_class.py index 5f96af412..0eb0f421e 100644 --- a/tests/test_data_class.py +++ b/tests/test_data_class.py @@ -284,6 +284,105 @@ 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_zero_init_unit_are_nan_aware_on_sparse_columns(self): + # #726, completing the #479 sweep above: a simulated column can carry NaN at a row no + # exp point scores (one failed integration step, or a 0/0 observable at t=0). zero + # reduces the WHOLE column, so one NaN made the mean NaN -> every centered value NaN -> + # the std NaN; init and unit read row 0 specifically, so a NaN there did the same. The + # column was then scored as a failed simulation. Reduce over the measured points, and + # take the first MEASURED row as the baseline. obs1 = [NaN, 3, 5, 7]. + lines = ['# x obs1\n', ' 0 nan\n', ' 1 3\n', ' 2 5\n', ' 3 7\n'] + + def mk(): + d = data.Data() + d.data = d._read_file_lines(lines, r'\s+') + return d + + # zero: centred on the measured mean (5) and scaled by their std; NaN row stays NaN. + dz = mk(); dz.normalize_to_zero(bc=False) + npt.assert_allclose(dz.data[1:, 1], np.array([-2., 0., 2.]) / np.std([3., 5., 7.])) + assert np.isnan(dz.data[0, 1]) + npt.assert_allclose(dz.normalization['obs1'][0].scale, np.std([3., 5., 7.])) + + # init: divides by the first MEASURED value (3, row 1), not the NaN row 0. + di = mk(); di.normalize_to_init() + npt.assert_allclose(di.data[1:, 1], np.array([1., 5. / 3., 7. / 3.])) + assert np.isnan(di.data[0, 1]) + assert di.normalization['obs1'][0].scale == 3.0 and di.normalization['obs1'][0].ref_row == 1 + + # unit: baseline is row 1 (the first measured), then /nanmax-after-baseline(=4). + du = mk(); du.normalize_to_unit_scale() + npt.assert_allclose(du.data[1:, 1], np.array([0., 0.5, 1.])) + assert np.isnan(du.data[0, 1]) + rec = du.normalization['obs1'][0] + # The recorded baseline_row must be the row actually subtracted -- the gradient reads + # the baseline back by index (gradient/assembly.py), so 0 here would be a wrong row. + assert rec.baseline_row == 1 and rec.ref_row == 3 and rec.scale == 4.0 + + @pytest.mark.parametrize('method', ['peak', 'init', 'zero', 'unit', ('floor', 0.03)], + ids=['peak', 'init', 'zero', 'unit', 'floor']) + def test_a_sparse_column_normalizes_like_the_dense_column_of_its_measured_values(self, method): + # The invariant #479 and #726 together are reaching for, stated once and checked over + # random draws rather than hand-computed cases: a NaN is missing data, so normalizing a + # column with NaN rows must leave its measured rows exactly where normalizing a dense + # column of just those values would put them. Every reduction (max, min, mean, std) and + # every reference row (argmax, argmin, the init/unit baseline) has to skip the NaNs for + # this to hold -- one that does not shifts or NaNs the whole column and the comparison + # fails. This also fixes the degenerate cases in place: a sparse column with a single + # measured value behaves like a one-row dense column (it does not acquire some separate + # sparse-only behaviour), which is what makes the rule easy to state to a user. + rng = np.random.default_rng(20260916) + compared = 0 + for _ in range(40): + n = int(rng.integers(3, 9)) + values = rng.normal(5.0, 3.0, n) + missing = rng.random(n) < 0.35 + if missing.all() or (~missing).sum() < 2: + continue + + def mk(vals): + d = data.Data() + d.data = d._read_file_lines( + ['# x obs1\n'] + [' %g %r\n' % (i, float(v)) for i, v in enumerate(vals)], + r'\s+') + return d + + sparse = values.copy() + sparse[missing] = np.nan + with warnings.catch_warnings(): + warnings.simplefilter('ignore') # degenerate draws divide by 0 as ever + ds, dd = mk(sparse), mk(values[~missing]) + ds.normalize(method) + dd.normalize(method) + npt.assert_allclose(ds.data[~missing, 1], dd.data[:, 1], rtol=1e-12, atol=1e-12) + assert np.isnan(ds.data[missing, 1]).all() # and the missing rows stay missing + compared += 1 + assert compared >= 20 # the draws actually exercised the comparison + + def test_an_unmeasured_column_is_left_alone_instead_of_raising(self): + # #726: a column with no measured point has no peak, no min, no mean and no baseline. + # np.nanargmax/nanargmin raise ValueError('All-NaN slice encountered') rather than + # returning a sentinel, so peak/floor/unit crashed outright and zero/init reached the + # same state by poisoning. Normalization is a transform; whether an all-NaN simulated + # column is a FAILURE is scoring's call, and scoring already makes it (evaluate returns + # None). Two callers do not wrap normalize in #388's handler -- model_check.run_check, + # which has its own 'Simulation contained NaN or Inf values' message three lines later, + # and the config-load floor on EXPERIMENTAL data -- so raising here escaped as a numpy + # traceback. Leave the column untouched and record the identity. + for method in ('peak', 'init', 'zero', 'unit', ('floor', 0.03)): + d = data.Data() + d.data = d._read_file_lines(['# x obs1 obs2\n', + ' 0 1 nan\n', + ' 1 2 nan\n'], r'\s+') + d.normalize(method) + assert np.isnan(d.data[:, 2]).all(), method # untouched, still all-NaN + npt.assert_allclose(d.data[:, 0], np.array([0., 1.])) + record = d.normalization['obs2'][0] + assert record.scale == 1.0 and record.rho == 0.0, method # the identity + # The measured sibling column is normalized as usual -- one dead column does not + # disturb the rest of the file. + assert not np.isnan(d.data[:, 1]).any(), method + 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 diff --git a/tests/test_gradient_assembly.py b/tests/test_gradient_assembly.py index 98ac9f2ff..bdeeb4fb6 100644 --- a/tests/test_gradient_assembly.py +++ b/tests/test_gradient_assembly.py @@ -1076,6 +1076,63 @@ def normed_of(column): np.testing.assert_allclose(res.jacobian[:, 0], fd / sigma, rtol=1e-5, atol=1e-7) +@pytest.mark.parametrize('nan_row', [0, 2], ids=['nan-at-baseline-row', 'nan-mid-column']) +@pytest.mark.parametrize('chain', [['peak'], ['init'], ['zero'], ['unit']], + ids=['peak', 'init', 'zero', 'unit']) +def test_normalization_chain_rule_matches_finite_difference_with_an_unscored_nan(chain, nan_row): + """#726: the threaded derivative still matches the FD oracle when the simulated column + carries a NaN at a row no experimental point scores. + + A simulation can legitimately produce NaN at some output rows -- one failed integration + step, or an observable that is 0/0 at t=0 -- while the exp file measures only a few of the + simulated times, so that row is often never read. ``normalize_to_zero`` used to make the + whole column NaN from it (and ``init``/``unit`` did when the NaN sat on the baseline row), + so the parameter set was discarded as a failed simulation; now the measured rows normalize + and the column scores, which means it reaches the gradient. The z-score's two reductions + (``s_bar`` and ``∂σ/∂θ``) therefore have to run over the measured rows only -- over all + rows they return NaN and poison the whole Jacobian. + + The NaN is placed on the **baseline** row as well as mid-column, because ``init`` and + ``unit`` read that row specifically and now take the first *measured* one instead. + """ + raw = np.array([2.0, 9.0, 5.0, 3.0]) + dk = np.array([0.5, -2.0, 1.3, -0.7]) + sigma = 1.0 + raw_nan = raw.copy() + raw_nan[nan_row] = np.nan + method = chain[0] + + sim = _sim_with_sensitivities(raw_nan.copy(), d_param=dk) + sim.normalize(method) + # The NaN stays on its own row; every other row carries a real normalized value. + normed_col = sim.data[:, 1] + assert np.isnan(normed_col[nan_row]) + assert not np.isnan(np.delete(normed_col, nan_row)).any() + + # That row is unmeasured in the exp file, so scoring skips it (objective.py's row skip). + obs = np.zeros(4) + obs[nan_row] = np.nan + exp = _exp(obs, sigma) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) + free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) + + res = assemble_gaussian_gradient(ChiSquareObjective(), [(sim, exp, routing)], free) + + def normed_of(column): + d = Data.from_columns(np.column_stack([TIMES, column]), ['time', 'Stot']) + d.normalize(method) + return d.data[:, 1] + + h = 1e-6 + fd = (normed_of(raw_nan + h * dk) - normed_of(raw_nan - h * dk)) / (2.0 * h) + # The Jacobian carries one row per SCORED observation, so the skipped row is absent. + scored = [r for r in range(4) if r != nan_row] + assert res.jacobian.shape[0] == len(scored) + assert not np.isnan(res.jacobian[:, 0]).any() # the regression: NaN Jacobian + np.testing.assert_allclose(res.jacobian[:, 0], fd[scored] / sigma, + rtol=1e-5, atol=1e-7) + + def test_floor_normalization_closed_form(): """A floor (ADR-0066, #533) is additive and separable: ``x' = x + rho*max(x)``, so its threaded derivative is ``∂x'_i/∂θ = s_i + rho*s_argmax`` -- every row picks up the *same* diff --git a/tests/test_model_check.py b/tests/test_model_check.py index 3c422493f..b3180614b 100644 --- a/tests/test_model_check.py +++ b/tests/test_model_check.py @@ -39,7 +39,7 @@ import pytest -from .context import algorithms, printing, pset +from .context import algorithms, data, printing, pset import pybnf.pybnf as pybnf_main @@ -394,6 +394,45 @@ def test_none_score_prints_nan_message_and_skips_constraints(monkeypatch, capsys assert 'inst' not in holder # constraint block skipped +# 'peak'/'unit' are the whole-fit string form; the floor is spelled as the per-data-key dict +# Configuration emits ({data_key: [(transform, [cols])]}), which is the only shape +# Result.normalize forwards a parameterized transform through -- a bare ('floor', 0.03) tuple +# matches neither of its branches and would silently normalize nothing. +@pytest.mark.parametrize('method', ['peak', 'unit', {'s': [(('floor', 0.03), ['obs1'])]}], + ids=['peak', 'unit', 'floor']) +def test_all_nan_column_reaches_the_nan_message_instead_of_a_traceback(monkeypatch, capsys, + method): + """#726: an all-NaN simulated column must not raise out of ``normalize``. + + ``run_check`` calls ``result.normalize(...)`` *outside* the ``try`` that follows it (that + one covers only ``postprocess_data``), and neither ``run_check`` nor its caller in + ``pybnf.py`` has a handler -- unlike the two fit scoring paths, which wrap normalize in + #388's blanket except. A column with no measured point has no peak, min or baseline, so + ``np.nanargmax`` / ``np.nanargmin`` raised ``ValueError('All-NaN slice encountered')`` and + ``--check-simulation`` on a broken model died with a numpy traceback -- bypassing the + diagnostic three lines further down that the command exists to print. + + This uses the **real** ``Result.normalize`` over a real ``Data`` (not ``_SpyResult``'s + recorder), because the defect is in the normalization itself. + """ + sim = data.Data() + sim.data = sim._read_file_lines(['# t obs1\n', ' 0 nan\n', ' 1 nan\n'], r'\s+') + result = algorithms.core.Result(pset.PSet([]), {'m': {'s': sim}}, 'check') + _patch_job(monkeypatch) + _patch_run_job(monkeypatch, result) + holder = _patch_counter(monkeypatch, fail_count=0) + # score=None is what evaluate_multiple returns for a NaN prediction -- the failed-simulation + # path run_check already handles. Reaching it is the whole point. + mc = _make_check(objective=_SpyObjective(score=None), constraints=(_FakeCset(2),), + normalization=method) + + mc.run_check() # must not raise + + out = capsys.readouterr().out + assert 'NaN or Inf' in out + assert 'inst' not in holder # returned before the constraint block, as designed + + # =========================================================================== # # run_check: constraint reporting block # =========================================================================== #