From a4864518b37a45abfa3d7fb5b187c57b4cbe3f86 Mon Sep 17 00:00:00 2001 From: mikessh Date: Sun, 11 Oct 2026 01:34:11 +0300 Subject: [PATCH] Stream component and grouped score audit in one native batch Co-Authored-By: Codex --- docs/gapblock.rst | 16 ++++++++++ include/seqtree/seqtree.hpp | 14 +++++++++ python/seqtree/gapblock.py | 27 ++++++++++++++++- src/_bindings.cpp | 6 ++++ src/gapblock.cpp | 37 +++++++++++++++++++++++ tests/python/test_gapblock_group_audit.py | 21 +++++++++++++ 6 files changed, 120 insertions(+), 1 deletion(-) create mode 100644 tests/python/test_gapblock_group_audit.py diff --git a/docs/gapblock.rst b/docs/gapblock.rst index 3630a60..7fdd412 100644 --- a/docs/gapblock.rst +++ b/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. diff --git a/include/seqtree/seqtree.hpp b/include/seqtree/seqtree.hpp index dd7d6c4..16dbeb4 100644 --- a/include/seqtree/seqtree.hpp +++ b/include/seqtree/seqtree.hpp @@ -255,6 +255,20 @@ std::vector> gapblock_ball_batch(const std::vector // 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, std::vector> +gapblock_group_audit_batch(const std::vector& queries, + const std::vector& refs, const std::vector& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights = {}, + const std::vector& query_groups = {}, + const std::vector& reference_groups = {}, + const std::vector>& group_distances = {}, + const std::vector>& query_weights = {}, + const std::vector>& reference_weights = {}); + std::vector> gapblock_topk_batch(const std::vector& queries, const std::vector& refs, uint32_t k, Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, const std::vector& prior, diff --git a/python/seqtree/gapblock.py b/python/seqtree/gapblock.py index 35ccf73..ab2e7bf 100644 --- a/python/seqtree/gapblock.py +++ b/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, diff --git a/src/_bindings.cpp b/src/_bindings.cpp index d91575c..e7b4a76 100644 --- a/src/_bindings.cpp +++ b/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{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, nb::arg("position_weights_by_length")=std::vector>{}, nb::arg("query_group_ids")=std::vector{}, nb::arg("reference_group_ids")=std::vector{}, nb::arg("group_distances")=std::vector>{}, nb::arg("query_position_weights")=std::vector>{}, nb::arg("reference_position_weights")=std::vector>{}, "Counts and integer linear mass in one bounded traversal, GIL released."); + m.def("gapblock_group_audit_batch", [](const std::vector& queries, const std::vector& refs, const std::vector& thresholds, const std::string& alphabet, const std::optional& matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact, const std::vector>& weights, const std::vector& query_groups, const std::vector& reference_groups, const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& 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{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, nb::arg("position_weights_by_length")=std::vector>{}, nb::arg("query_group_ids")=std::vector{}, nb::arg("reference_group_ids")=std::vector{}, nb::arg("group_distances")=std::vector>{}, nb::arg("query_position_weights")=std::vector>{}, nb::arg("reference_position_weights")=std::vector>{}, "Streaming component/full ball counts and minimum full score within the component ball, GIL released."); + m.def("gapblock_ball_batch", [](const std::vector& queries, const std::vector& refs, const std::vector& thresholds, const std::string& alphabet, const std::optional& matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact, const std::vector>& weights, const std::vector& query_groups, const std::vector& reference_groups, const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { const auto alph = parse_alphabet(alphabet); nb::gil_scoped_release release; diff --git a/src/gapblock.cpp b/src/gapblock.cpp index 1e67733..d59fd04 100644 --- a/src/gapblock.cpp +++ b/src/gapblock.cpp @@ -371,6 +371,43 @@ std::vector> gapblock_ball_batch(const std::vector return out; } +std::tuple, std::vector, std::vector> +gapblock_group_audit_batch(const std::vector& queries, + const std::vector& refs, const std::vector& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights, + const std::vector& query_groups, const std::vector& reference_groups, + const std::vector>& group_distances, + const std::vector>& query_weights, + const std::vector>& 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, std::vector> out{ + std::vector(queries.size()), std::vector(queries.size()), + std::vector(queries.size(), -1)}; + auto& [before, after, nearest] = out; + parallel_for(queries.size(), threads, + [&] { return std::vector(size_t(in.longest) + 1); }, + [&](size_t i, std::vector& scratch) { + if (thresholds[i] < 0) return; + int32_t minimum = std::numeric_limits::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> gapblock_topk_batch(const std::vector& queries, const std::vector& refs, uint32_t k, Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, const std::vector& prior, diff --git a/tests/python/test_gapblock_group_audit.py b/tests/python/test_gapblock_group_audit.py new file mode 100644 index 0000000..a505942 --- /dev/null +++ b/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,[])