From 15e0afacbcf11e9ed676141c17cef7a1dc8ad180 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Francisco=20Miguel=20P=C3=A9rez=20Canales?= Date: Thu, 25 Jun 2026 06:18:24 +0200 Subject: [PATCH] fix(graph): sparse prop='fill' bit-parity with dense + upstream (changes metrics, needs review) (#20) ## Summary (CORRECTNESS-CRITICAL: changes metric numbers , please review) The `CAFAEVAL_SPARSE` env switched not only the confusion-matrix kernel but also the **propagation algorithm**. Under `prop="fill"` (PROTEA's production default), the sparse path (`_propagate_sparse_pushup`) and the dense serial sweep produced **different propagated prediction matrices** on real DAG inputs, shifting Fmax: | Setting | sparse (before) | dense (before) | |, -|, -|, -| | NK Fmax | 0.58859 | 0.58914 | | LK Fmax | 0.61612 | 0.61391 | | PK Fmax | 0.40195 | 0.40195 | ## Which impl was canonical The **dense serial sweep is canonical**: it is bit-identical to upstream cafaeval `prop="fill"` (BioComputingUP/CAFA-evaluator and claradepaolis/CAFA-evaluator-PK, the fork's ancestors). Verified on 5 hand-worked DAGs plus a 200-trial randomized stress; the dense path matched the upstream oracle in 100% of cases. The **sparse pushup was the buggy one**. `fill` is a stepwise leaves-to-roots recurrence where a non-zero intermediate term *blocks* deeper descendants (it keeps its own value, and that value, not the descendant's, flows up to the parent). The pushup scattered every input non-zero straight into every ancestor and group-maxed, ignoring those blockers and overshooting. Worked counter-example (chain leaf->mid->root, leaf=1.0, mid=0.5, root=0): - upstream / dense fill: root = **0.5** (mid blocks the leaf) - sparse pushup: root = **1.0** (leaf scattered straight to root), WRONG This matches the observed "sparse=1.0 where dense=0.86" divergence, concentrated at high tau on biological_process multi-term shared-ancestor rows. ## The fix (which impl was fixed) The **sparse** path is fixed. `mode="fill"` is now routed to a new `_propagate_sparse_fill` that reproduces the recurrence sparsely (walk `ont.order`, per-row max over each term's children's *current* values, never overwrite originally-non-zero cells). The fast scatter pushup is retained for `mode="max"`, where it is proven correct. Minimal, localized change to the propagation path only. ## After: bit-identical - **Propagated prediction matrices: bit-identical** between sparse and dense, `max|diff| = 0.000e+00` across BPO/CCO/MFO (including the BPO rows where the bug lived). - NK/LK/PK Fmax now identical between both paths and equal to the dense (canonical) values: | Setting | sparse (after) | dense (after) | |, -|, -|, -| | NK Fmax | 0.5891445 | 0.5891445 | | LK Fmax | 0.6139093 | 0.6139093 | | PK Fmax | 0.4019451 | 0.4019451 | (Residual ~1e-16 / 1 ULP in some weighted count columns is float reduction-order noise in the counting kernel, downstream of the now-identical propagated matrices, below any metric significance.) - Sparse stays ~5x faster than dense on the fill path (no perf regression). ## Test added `tests/test_propagation_fill_parity.py`: the worked blocked-chain counter-example, shared-ancestor and high-tau BPO-shaped rows, plus a 200-trial randomized DAG stress, all asserting `sparse == dense == upstream-oracle` bit-for-bit. Closes the coverage gap left by the existing self-parity test, which only exercises `prop="max"`. ## Local CI ruff clean; mypy clean (py3.12 target, matching the CI runner); full `tests/` + `tests/diff` (Phase B) green. ## Deploy note PROTEA pins `cafaeval @ git+main`. A redeploy will pick this up once it lands on `main`. I did NOT bump the pin. --- src/cafaeval/graph.py | 108 ++++++++++-- tests/test_propagation_fill_parity.py | 241 ++++++++++++++++++++++++++ 2 files changed, 338 insertions(+), 11 deletions(-) create mode 100644 tests/test_propagation_fill_parity.py diff --git a/src/cafaeval/graph.py b/src/cafaeval/graph.py index 754e2cb..fd7d364 100644 --- a/src/cafaeval/graph.py +++ b/src/cafaeval/graph.py @@ -1,7 +1,7 @@ import numpy as np import logging import os -from collections import deque +from collections import defaultdict, deque logger = logging.getLogger(__name__) logger.addHandler(logging.NullHandler()) @@ -275,6 +275,92 @@ def _ancestors_csr(ont): return indptr, indices +def _propagate_sparse_fill(matrix, ont, triples=None): + """Sparse, in-place ``mode='fill'`` propagation, bit-identical to the dense + serial sweep (and to upstream cafaeval ``prop='fill'``). + + ``fill`` is **not** a scatter-to-all-ancestors group-max. It is a stepwise + recurrence over the DAG (leaves -> roots): + + v[t] = orig[t] if orig[t] != 0 (blocker) + v[t] = max over children c of v[c] if orig[t] == 0 + + where ``v[c]`` is each child's *already-finalised* value. A non-zero + intermediate node "absorbs" deeper descendants: it keeps its own value and + that value (not the descendant's) is what flows up. The plain pushup, which + scatters every input non-zero directly into every ancestor, ignores those + intermediate blockers and overshoots (e.g. it would set a parent to a + deep leaf's score even though a lower-scored non-zero node sits between + them). See ``tests/test_propagation_fill_parity.py`` for the worked + counter-example and the sparse/dense/upstream parity gate. + + We reproduce the recurrence sparsely by walking ``ont.order`` (the + topological leaves->roots order) and, for each term that has children, + taking the per-row max over its children's current values. Only terms that + actually have active descendants in a given row contribute work, so the + cost tracks the propagated non-zero count rather than ``n_prot * n_terms``. + + If ``triples`` is ``(nz_rows, nz_cols, nz_scores)`` the caller already knows + the input non-zero positions and we skip the dense ``np.nonzero`` scan. + """ + n_prot, n_terms = matrix.shape + if n_prot == 0 or n_terms == 0: + return + + if triples is not None: + nz_rows, nz_cols, nz_scores = triples + nz_rows = np.asarray(nz_rows, dtype=np.int64) + nz_cols = np.asarray(nz_cols, dtype=np.int64) + nz_scores = np.asarray(nz_scores) + else: + nz_rows, nz_cols = np.nonzero(matrix) + if nz_rows.size == 0: + return + nz_scores = matrix[nz_rows, nz_cols] + + if nz_rows.size == 0: + return + + chi_indptr = ont._chi_indptr + chi_idx = ont._chi_idx + + # ``current[t]`` holds the current (post-fill) ``{row: value}`` map for + # term ``t``; ``orig_rows[t]`` is the set of rows that were originally + # non-zero at ``t`` (blockers that must never be overwritten). Both are + # seeded from the input non-zeros and grown as we sweep upward. + current = defaultdict(dict) + orig_rows = defaultdict(set) + for r, c, s in zip(nz_rows.tolist(), nz_cols.tolist(), nz_scores.tolist()): + current[c][r] = s + orig_rows[c].add(r) + + for t in ont.order: + c0 = chi_indptr[t] + c1 = chi_indptr[t + 1] + if c1 == c0: + continue + # Per-row max over this term's children's current values. + row_max = {} + for child in chi_idx[c0:c1]: + child_vals = current.get(int(child)) + if not child_vals: + continue + for r, v in child_vals.items(): + if v > row_max.get(r, -np.inf): + row_max[r] = v + if not row_max: + continue + blocked = orig_rows.get(t, ()) + cur_t = current[t] + for r, v in row_max.items(): + if r in blocked: + # Originally non-zero -> keep the input value (it already feeds + # parents via ``cur_t``); descendants do not overwrite it. + continue + cur_t[r] = v + matrix[r, t] = v + + def _propagate_sparse_pushup(matrix, ont, mode, triples=None): """Sparse alternative to ``_propagate_serial``. @@ -292,10 +378,15 @@ def _propagate_sparse_pushup(matrix, ont, mode, triples=None): large win for ground-truth matrices where the non-zero density is ~1e-4 and ``np.nonzero`` would otherwise touch every cell. - For ``mode='fill'`` the update semantics are "only overwrite cells that - were originally zero"; this is restored by snapshotting the input - non-zero positions and writing their original values back at the end. + ``mode='fill'`` is delegated to ``_propagate_sparse_fill`` because the + zero-only-overwrite semantics are a stepwise recurrence, not a + scatter-to-ancestors group-max; the pushup below handles ``mode='max'`` + only. """ + if mode == "fill": + _propagate_sparse_fill(matrix, ont, triples=triples) + return + n_prot, n_terms = matrix.shape if n_prot == 0 or n_terms == 0: return @@ -351,17 +442,12 @@ def _propagate_sparse_pushup(matrix, ont, mode, triples=None): out_rows = unique_flat // n_terms out_cols = unique_flat % n_terms - # Max against what is currently in the matrix. For mode='max' this is the - # final answer; for mode='fill' we restore original non-zeros below. + # Max against what is currently in the matrix; for mode='max' this is the + # final answer (mode='fill' is handled by _propagate_sparse_fill above). current = matrix[out_rows, out_cols] np.maximum(current, group_max, out=current) matrix[out_rows, out_cols] = current - if mode == "fill": - # Only zero cells may be filled by propagation; originally non-zero - # cells keep their input value regardless of what descendants say. - matrix[nz_rows, nz_cols] = nz_scores - def propagate_to_coo(triples, ont, mode="max"): """Sparse-native propagation that never materialises a dense matrix. diff --git a/tests/test_propagation_fill_parity.py b/tests/test_propagation_fill_parity.py new file mode 100644 index 0000000..810c238 --- /dev/null +++ b/tests/test_propagation_fill_parity.py @@ -0,0 +1,241 @@ +"""Sparse vs dense parity for ``mode='fill'`` propagation. + +``prop='fill'`` is the production default in the PROTEA stack, and the metric +numbers (f_micro_w, Fmax) depend on it. The sparse path (``CAFAEVAL_SPARSE=1``) +and the dense serial sweep (``CAFAEVAL_SPARSE=0``) must produce a BIT-IDENTICAL +propagated matrix, both matching the canonical upstream cafaeval ``fill``. + +The original sparse pushup scattered every input non-zero straight into every +ancestor and group-maxed. That is correct for ``mode='max'`` but WRONG for +``fill``: ``fill`` is a stepwise recurrence where a non-zero intermediate node +"blocks" deeper descendants (it keeps its own value and that value, not the +descendant's, flows up to the parent). The pushup ignored those blockers and +overshot, e.g. setting a parent to a deep leaf's score even though a +lower-scored non-zero node sat between them. ``_propagate_sparse_fill`` walks +the topological order instead and is bit-identical to the dense sweep. +""" +from __future__ import annotations + +import os + +import numpy as np + +from cafaeval.graph import Graph, propagate + + +def _make_graph(term_ids, parent_edges): + """Build a minimal single-namespace Graph. + + ``term_ids`` fixes the term -> index order; ``parent_edges`` is a list of + ``(child_id, parent_id)`` (child ``is_a`` parent). + """ + children_parents = {t: [] for t in term_ids} + for child, parent in parent_edges: + children_parents[child].append(parent) + terms_dict = { + t: {"name": t, "namespace": "test", "def": "", + "alt_id": [], "rel": children_parents[t]} + for t in term_ids + } + return Graph("test", terms_dict) + + +def _upstream_fill(matrix, ont): + """Verbatim canonical upstream cafaeval ``fill`` (BioComputingUP / + claradepaolis CAFA-evaluator-PK), kept here as an independent oracle so the + test does not merely compare the fork against itself. + """ + m = matrix.copy() + order = np.asarray(ont.order) + n = ont.idxs + # dense parent adjacency: dag[i, j] == 1 iff j is a parent of i + dag = np.zeros((n, n), dtype=bool) + for t in range(n): + for p in ont._par_idx[ont._par_indptr[t]:ont._par_indptr[t + 1]]: + dag[t, p] = True + deepest = np.where(np.sum(m[:, order], axis=0) > 0)[0][0] + order_ = np.delete(order, list(range(0, deepest))) + for i in order_: + children = np.where(dag[:, i] != 0)[0] + if children.size > 0: + cols = np.concatenate((children, [i])) + rows = np.where(m[:, i] == 0)[0] + if rows.size: + idx = np.ix_(rows, cols) + m[rows, i] = m[idx].max(axis=1) + return m + + +def _propagate_fill(matrix, ont, sparse): + os.environ["CAFAEVAL_SPARSE"] = "1" if sparse else "0" + for attr in ("_ancestors_csr", "_children_by_term"): + if hasattr(ont, attr): + delattr(ont, attr) + out = matrix.copy() + propagate(out, ont, ont.order, mode="fill", parallel=0) + return out + + +def _run_all(ont, base): + dense = _propagate_fill(base, ont, sparse=False) + sparse = _propagate_fill(base, ont, sparse=True) + upstream = _upstream_fill(base, ont) + return dense, sparse, upstream + + +def test_fill_blocked_intermediate(): + """The exact counter-example the buggy pushup got wrong. + + Chain leaf -> mid -> root with mid originally non-zero and *lower* than the + leaf. ``fill`` must give root = mid's value (0.5), not the leaf's (1.0): + the non-zero mid blocks the leaf. The pushup produced 1.0. + """ + ont = _make_graph(["leaf", "mid", "root"], + [("leaf", "mid"), ("mid", "root")]) + idx = {t: ont.terms_dict[t]["index"] for t in ("leaf", "mid", "root")} + base = np.zeros((1, ont.idxs)) + base[0, idx["leaf"]] = 1.0 + base[0, idx["mid"]] = 0.5 + + dense, sparse, upstream = _run_all(ont, base) + + expected = base.copy() + expected[0, idx["root"]] = 0.5 # mid (0.5) flows up, leaf (1.0) is blocked + + np.testing.assert_array_equal(dense, expected) + np.testing.assert_array_equal(sparse, expected) + np.testing.assert_array_equal(upstream, expected) + np.testing.assert_array_equal(sparse, dense) + + +def test_fill_shared_ancestor_multi_term(): + """Multi-term row with a shared ancestor and a non-zero intermediate. + + Two leaves feed a mid (originally 0.6); mid feeds the root. The root must + take the mid's blocked value (0.6), not the higher leaf (0.9). + """ + ont = _make_graph(["l1", "l2", "mid", "root"], + [("l1", "mid"), ("l2", "mid"), ("mid", "root")]) + idx = {t: ont.terms_dict[t]["index"] for t in ("l1", "l2", "mid", "root")} + base = np.zeros((1, ont.idxs)) + base[0, idx["l1"]] = 0.9 + base[0, idx["l2"]] = 0.3 + base[0, idx["mid"]] = 0.6 + + dense, sparse, upstream = _run_all(ont, base) + + expected = base.copy() + expected[0, idx["root"]] = 0.6 + + np.testing.assert_array_equal(sparse, expected) + np.testing.assert_array_equal(dense, expected) + np.testing.assert_array_equal(upstream, expected) + + +def test_fill_pure_chain_no_block(): + """No intermediate blocker: a single leaf value fills the whole chain.""" + ont = _make_graph(["leaf", "mid", "root"], + [("leaf", "mid"), ("mid", "root")]) + idx = {t: ont.terms_dict[t]["index"] for t in ("leaf", "mid", "root")} + base = np.zeros((1, ont.idxs)) + base[0, idx["leaf"]] = 0.8 + + dense, sparse, upstream = _run_all(ont, base) + + expected = base.copy() + expected[0, idx["mid"]] = 0.8 + expected[0, idx["root"]] = 0.8 + + np.testing.assert_array_equal(sparse, expected) + np.testing.assert_array_equal(dense, expected) + np.testing.assert_array_equal(upstream, expected) + + +def test_fill_two_level_block(): + """Blocker two levels below the root still caps everything above it.""" + ont = _make_graph(["leaf", "m1", "m2", "root"], + [("leaf", "m1"), ("m1", "m2"), ("m2", "root")]) + idx = {t: ont.terms_dict[t]["index"] for t in ("leaf", "m1", "m2", "root")} + base = np.zeros((1, ont.idxs)) + base[0, idx["leaf"]] = 1.0 + base[0, idx["m1"]] = 0.5 + + dense, sparse, upstream = _run_all(ont, base) + + expected = base.copy() + expected[0, idx["m2"]] = 0.5 + expected[0, idx["root"]] = 0.5 + + np.testing.assert_array_equal(sparse, expected) + np.testing.assert_array_equal(dense, expected) + np.testing.assert_array_equal(upstream, expected) + + +def test_fill_high_tau_bpo_shaped(): + """A wide, deep, multi-term row shaped like the high-tau biological_process + rows where the divergence was originally observed: many leaves, several + shared ancestors, a few non-zero intermediate blockers at high score. + """ + # 3 leaves -> 2 mids -> 1 hub -> root. Cross edges make mids share leaves. + term_ids = ["la", "lb", "lc", "ma", "mb", "hub", "root"] + edges = [ + ("la", "ma"), ("lb", "ma"), ("lb", "mb"), ("lc", "mb"), + ("ma", "hub"), ("mb", "hub"), ("hub", "root"), + ] + ont = _make_graph(term_ids, edges) + idx = {t: ont.terms_dict[t]["index"] for t in term_ids} + base = np.zeros((1, ont.idxs)) + # leaves at high tau-ish scores + base[0, idx["la"]] = 0.99 + base[0, idx["lb"]] = 0.80 + base[0, idx["lc"]] = 0.75 + # mb is a non-zero blocker LOWER than its descendant lb (0.80) + base[0, idx["mb"]] = 0.50 + + dense, sparse, upstream = _run_all(ont, base) + + # ma: zero -> max(la=0.99, lb=0.80) = 0.99 + # mb: blocked at 0.50 (lb/lc do not overwrite it) + # hub: zero -> max(ma=0.99, mb=0.50) = 0.99 + # root: zero -> hub = 0.99 + expected = base.copy() + expected[0, idx["ma"]] = 0.99 + expected[0, idx["hub"]] = 0.99 + expected[0, idx["root"]] = 0.99 + + np.testing.assert_array_equal(sparse, expected) + np.testing.assert_array_equal(dense, expected) + np.testing.assert_array_equal(upstream, expected) + # The headline invariant: sparse and dense agree bit-for-bit. + np.testing.assert_array_equal(sparse, dense) + + +def test_fill_randomized_dag_stress(): + """Random DAGs + random multi-row sparse inputs: sparse, dense and the + upstream oracle must agree bit-for-bit on every trial. + """ + rng = np.random.default_rng(20240625) + for _ in range(200): + n = int(rng.integers(5, 16)) + term_ids = [f"T{i}" for i in range(n)] + # i < j edge guarantees acyclicity; index order respects topology + edges = [] + for i in range(n): + candidates = list(range(i + 1, n)) + rng.shuffle(candidates) + for j in candidates[: int(rng.integers(0, 3))]: + edges.append((term_ids[i], term_ids[j])) + ont = _make_graph(term_ids, edges) + + n_rows = int(rng.integers(1, 6)) + base = np.zeros((n_rows, ont.idxs)) + for r in range(n_rows): + n_cells = int(rng.integers(1, n)) + for c in rng.choice(n, size=n_cells, replace=False): + base[r, c] = round(float(rng.uniform(0.01, 1.0)), 4) + if base.sum() == 0: # propagate raises on an all-zero matrix + continue + + dense, sparse, upstream = _run_all(ont, base) + np.testing.assert_array_equal(sparse, dense) + np.testing.assert_array_equal(sparse, upstream)