fix(graph): propagation fill bit-parity to main (PROTEA pin)

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.

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)