Stream component and grouped score audits without hit materialization

docs/gapblock.rst

@@ -414,3 +414,19 @@ Supply both axes and do not combine them with ``position_weights_by_length``.
 ``paired_sum_count_batch`` accepts separate ``_alpha`` and ``_beta`` versions of
 these arguments. Each lane retains its own coordinates; distances sum against the
 same linked reference row. These are geometry options, not confidence estimates.
+
+Streaming component and group audit
+-----------------------------------
+
+``group_audit_batch`` accepts the same geometry as ``ball_batch``, with one
+threshold per query. It returns three lists: the component count before group
+offsets, the full-score count, and the minimum full score within the component
+ball. An empty component ball has minimum ``-1``; duplicate references count
+separately and encoded exact exclusion is preserved.
+
+If the component penalty is c and the nonnegative group offset is g, every
+full-score hit c+g <= t is also a component hit c <= t. The reduction visits
+references once with shared encoded inputs and per-worker alignment scratch;
+output occupies O(Q) storage without sparse hit objects or a dense matrix.
+The bounded minimum is not an unbounded nearest result. No statistical model
+or cutoff recommendation is introduced.

include/seqtree/seqtree.hpp

@@ -255,6 +255,20 @@ std::vector<std::vector<Hit>> gapblock_ball_batch(const std::vector<std::string>
 // Exhaustive bounded reductions using the same gapblock cell. Top-k ties use
 // (penalty, ref_id); paired penalty is max(lane scores). Exact exclusion uses
 // encoded sequence identity (both lanes for pairs), never zero penalty alone.
+// Count the component ball before/after group offsets; minimum full score within it.
+// One threshold per query; nearest=-1 denotes an empty component ball. O(Q) output.
+std::tuple<std::vector<uint64_t>, std::vector<uint64_t>, std::vector<int32_t>>
+gapblock_group_audit_batch(const std::vector<std::string>& queries,
+        const std::vector<std::string>& refs, const std::vector<int32_t>& thresholds,
+        Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge,
+        const std::vector<int32_t>& prior, uint32_t width, int threads, bool exclude_exact,
+        const std::vector<std::vector<int32_t>>& weights = {},
+        const std::vector<int32_t>& query_groups = {},
+        const std::vector<int32_t>& reference_groups = {},
+        const std::vector<std::vector<int32_t>>& group_distances = {},
+        const std::vector<std::vector<int32_t>>& query_weights = {},
+        const std::vector<std::vector<int32_t>>& reference_weights = {});
+
 std::vector<std::vector<Hit>> gapblock_topk_batch(const std::vector<std::string>& queries,
         const std::vector<std::string>& refs, uint32_t k, Alphabet alphabet,
         const SubstitutionMatrix* matrix, int32_t go, int32_t ge, const std::vector<int32_t>& prior,

python/seqtree/gapblock.py

@@ -46,12 +46,13 @@
 from ._core import gapblock_count_batch as _gapblock_count_batch
 from ._core import gapblock_count_mass_batch as _gapblock_count_mass_batch
 from ._core import gapblock_ball_batch as _gapblock_ball_batch
+from ._core import gapblock_group_audit_batch as _gapblock_group_audit_batch
 from ._core import gapblock_paired_topk_batch as _gapblock_paired_topk_batch
 from ._core import gapblock_paired_count_batch as _gapblock_paired_count_batch
 from ._core import gapblock_paired_sum_count_batch as _gapblock_paired_sum_count_batch
 
 __all__ = [
-    "gapblock_score", "score_matrix", "ball_batch", "topk_batch", "count_batch", "paired_topk_batch", "paired_count_batch", "paired_sum_count_batch", "deletion_variants", "central_prior", "profile_prior",
+    "gapblock_score", "score_matrix", "ball_batch", "group_audit_batch", "topk_batch", "count_batch", "paired_topk_batch", "paired_count_batch", "paired_sum_count_batch", "deletion_variants", "central_prior", "profile_prior",
     "frame_prior", "positions_prior", "embed_in_frame", "gap_cost", "GapBlockIndex", "ScoreMatrix",
     "IslandProfile",
 ]
@@ -730,6 +731,30 @@ def count_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_exten
                             *_sequence_weights(query_position_weights, reference_position_weights, position_weights_by_length, q, r))
 
 
+def group_audit_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1,
+               gap_prior=None, alphabet="aa", threads=0, exclude_exact=False,
+               position_weights_by_length=None,
+               query_group_ids=None, reference_group_ids=None, group_distances=None,
+               query_position_weights=None, reference_position_weights=None):
+    """Stream component-ball count, full-score count and bounded nearest in one batch.
+
+    Geometry/offset options match ball_batch. Component excludes the nonnegative
+    group offset; full score includes it. For each query threshold, return three
+    lists: component count, full count and minimum full score among component hits.
+    Empty component balls have nearest=-1. Duplicate references count separately;
+    full encoded exact puncturing is identical. O(Q) output; no sparse hit objects.
+    """
+    q, r = list(queries), list(refs)
+    thresholds = [_batch_integer(v, "threshold", -(1 << 31)) for v in thresholds]
+    if len(thresholds) != len(q):
+        raise ValueError("group audit thresholds must have one value per query")
+    go, ge, cube, width, threads = _batch_options((q, r), matrix, gap_open, gap_extend, gap_prior, threads)
+    return _gapblock_group_audit_batch(q, r, thresholds, alphabet, matrix, go, ge, cube, width,
+                               threads, exclude_exact, _position_weights(position_weights_by_length, q, r),
+                            *_group_options(query_group_ids, reference_group_ids, group_distances, q, r),
+                            *_sequence_weights(query_position_weights, reference_position_weights, position_weights_by_length, q, r))
+
+
 def ball_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1,
                gap_prior=None, alphabet="aa", threads=0, exclude_exact=False,
                position_weights_by_length=None,

src/_bindings.cpp

@@ -824,6 +824,12 @@ NB_MODULE(_core, m) {
         return gapblock_count_mass_batch(queries, refs, thresholds, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights);
     }, nb::arg("queries"), nb::arg("refs"), nb::arg("thresholds"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector<int32_t>{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, nb::arg("position_weights_by_length")=std::vector<std::vector<int32_t>>{}, nb::arg("query_group_ids")=std::vector<int32_t>{}, nb::arg("reference_group_ids")=std::vector<int32_t>{}, nb::arg("group_distances")=std::vector<std::vector<int32_t>>{}, nb::arg("query_position_weights")=std::vector<std::vector<int32_t>>{}, nb::arg("reference_position_weights")=std::vector<std::vector<int32_t>>{}, "Counts and integer linear mass in one bounded traversal, GIL released.");
 
+    m.def("gapblock_group_audit_batch", [](const std::vector<std::string>& queries, const std::vector<std::string>& refs, const std::vector<int32_t>& thresholds, const std::string& alphabet, const std::optional<SubstitutionMatrix>& matrix, int32_t go, int32_t ge, const std::vector<int32_t>& prior, uint32_t width, int threads, bool exclude_exact, const std::vector<std::vector<int32_t>>& weights, const std::vector<int32_t>& query_groups, const std::vector<int32_t>& reference_groups, const std::vector<std::vector<int32_t>>& group_distances, const std::vector<std::vector<int32_t>>& query_weights, const std::vector<std::vector<int32_t>>& reference_weights) {
+        const auto alph = parse_alphabet(alphabet);
+        nb::gil_scoped_release release;
+        return gapblock_group_audit_batch(queries, refs, thresholds, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights);
+    }, nb::arg("queries"), nb::arg("refs"), nb::arg("thresholds"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector<int32_t>{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, nb::arg("position_weights_by_length")=std::vector<std::vector<int32_t>>{}, nb::arg("query_group_ids")=std::vector<int32_t>{}, nb::arg("reference_group_ids")=std::vector<int32_t>{}, nb::arg("group_distances")=std::vector<std::vector<int32_t>>{}, nb::arg("query_position_weights")=std::vector<std::vector<int32_t>>{}, nb::arg("reference_position_weights")=std::vector<std::vector<int32_t>>{}, "Streaming component/full ball counts and minimum full score within the component ball, GIL released.");
+
     m.def("gapblock_ball_batch", [](const std::vector<std::string>& queries, const std::vector<std::string>& refs, const std::vector<int32_t>& thresholds, const std::string& alphabet, const std::optional<SubstitutionMatrix>& matrix, int32_t go, int32_t ge, const std::vector<int32_t>& prior, uint32_t width, int threads, bool exclude_exact, const std::vector<std::vector<int32_t>>& weights, const std::vector<int32_t>& query_groups, const std::vector<int32_t>& reference_groups, const std::vector<std::vector<int32_t>>& group_distances, const std::vector<std::vector<int32_t>>& query_weights, const std::vector<std::vector<int32_t>>& reference_weights) {
         const auto alph = parse_alphabet(alphabet);
         nb::gil_scoped_release release;

src/gapblock.cpp

@@ -371,6 +371,43 @@ std::vector<std::vector<Hit>> gapblock_ball_batch(const std::vector<std::string>
     return out;
 }
 
+std::tuple<std::vector<uint64_t>, std::vector<uint64_t>, std::vector<int32_t>>
+gapblock_group_audit_batch(const std::vector<std::string>& queries,
+        const std::vector<std::string>& refs, const std::vector<int32_t>& thresholds,
+        Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge,
+        const std::vector<int32_t>& prior, uint32_t width, int threads, bool exclude_exact,
+        const std::vector<std::vector<int32_t>>& weights,
+        const std::vector<int32_t>& query_groups, const std::vector<int32_t>& reference_groups,
+        const std::vector<std::vector<int32_t>>& group_distances,
+        const std::vector<std::vector<int32_t>>& query_weights,
+        const std::vector<std::vector<int32_t>>& reference_weights) {
+    if (thresholds.size() != queries.size())
+        throw std::invalid_argument("group audit thresholds must have one value per query");
+    const auto in = prepare(queries, refs, alphabet, matrix, go, ge, prior, width, threads,
+                           weights, query_groups, reference_groups, group_distances,
+                           query_weights, reference_weights);
+    std::tuple<std::vector<uint64_t>, std::vector<uint64_t>, std::vector<int32_t>> out{
+        std::vector<uint64_t>(queries.size()), std::vector<uint64_t>(queries.size()),
+        std::vector<int32_t>(queries.size(), -1)};
+    auto& [before, after, nearest] = out;
+    parallel_for(queries.size(), threads,
+        [&] { return std::vector<int32_t>(size_t(in.longest) + 1); },
+        [&](size_t i, std::vector<int32_t>& scratch) {
+            if (thresholds[i] < 0) return;
+            int32_t minimum = std::numeric_limits<int32_t>::max();
+            for (size_t j = 0; j < refs.size(); ++j) {
+                if (exclude_exact && in.q[i] == in.r[j]) continue;
+                const auto score = score_cell(in, i, j, go, ge, prior, scratch.data());
+                if (score - group_offset(in, i, j) > thresholds[i]) continue;
+                ++before[i];
+                after[i] += score <= thresholds[i];
+                minimum = std::min(minimum, score);
+            }
+            if (before[i]) nearest[i] = minimum;
+        });
+    return out;
+}
+
 std::vector<std::vector<Hit>> gapblock_topk_batch(const std::vector<std::string>& queries,
         const std::vector<std::string>& refs, uint32_t k, Alphabet alphabet,
         const SubstitutionMatrix* matrix, int32_t go, int32_t ge, const std::vector<int32_t>& prior,

tests/python/test_gapblock_group_audit.py

@@ -0,0 +1,21 @@
+"""Streaming group audit preserves exact component/full balls without hit materialization."""
+import numpy as np
+import pytest
+from seqtree import gapblock
+
+def test_group_audit_matches_dense_components_and_threads():
+    q=['AA','CA','aa',''];r=['AA','CC','CAA','AA','']
+    opts=dict(gap_open=2,gap_extend=1,position_weights_by_length={0:[],2:[0,3],3:[2,1,0]})
+    group=dict(query_group_ids=[0,1,0,1],reference_group_ids=[0,1,0,1,1],group_distances=[[0,5],[3,0]])
+    base=np.asarray(gapblock.score_matrix(q,r,**opts));full=np.asarray(gapblock.score_matrix(q,r,**opts,**group))
+    limits=[0,3,10,-1]
+    for threads in [1,3]:
+        for excluded in [False,True]:
+            before,after,nearest=gapblock.group_audit_batch(q,r,limits,threads=threads,exclude_exact=excluded,**opts,**group)
+            for i,t in enumerate(limits):
+                valid=np.array([not excluded or q[i].upper()!=s.upper() for s in r]);selected=valid&(base[i]<=t)
+                assert before[i]==int(selected.sum())
+                assert after[i]==int((valid&(full[i]<=t)).sum())
+                assert nearest[i]==(int(full[i,selected].min()) if selected.any() else -1)
+    assert gapblock.group_audit_batch(['AA'],[],[0])==([0],[0],[-1])
+    with pytest.raises(ValueError):gapblock.group_audit_batch(q,r,[])