Add bounded exhaustive gap-block top-K and threshold counts

CHANGELOG.md

@@ -7,6 +7,15 @@ versioning: breaking changes need a **major** bump.
 
 ### Added
 
+- **1.0.2 candidate: bounded exhaustive gapblock reductions.** `topk_batch` keeps
+  stable `(score, ref_id)` top-k rows without a dense score matrix; `count_batch`
+  counts all ties at arbitrary per-query thresholds. Linked two-view variants use
+  the maximum lane penalty and preserve both lane scores. Exact identity exclusion
+  precedes selection and removes only both-view identity for linked inputs.
+  `Index.ref_seqs()` retrieves stored reference strings in ID order in one native
+  call. All scoring reuses the existing single-gap-block cell; no calibration,
+  domain payload or runtime dependency is introduced.
+
 - **`TextIndex` searches with insertions and deletions.** `search_batch` takes `max_indels`
   alongside `max_subs`, capped independently — `max_subs=2, max_indels=1` accepts two
   substitutions *and* one gap, not three edits of any kind.

docs/gapblock.rst

@@ -314,3 +314,44 @@ See also
 * :doc:`api` — full signatures for :mod:`seqtree.gapblock`, :mod:`seqtree.seeds` and
   :mod:`seqtree.layout`.
 * ``bench/bench_gapblock.py`` — reproduces the gap-freedom ladder table above.
+
+Bounded exhaustive ranking and counts
+-------------------------------------
+
+``topk_batch`` evaluates every reference with the same single-gap-block scorer as
+``score_matrix`` and returns at most k ``Hit`` rows per query, ordered by
+``(score, ref_id)``. Reference IDs retain input row order; repeated sequences remain
+separate rows. Hit edit fields are unset, rather than an alignment decomposition.
+``exclude_exact=True`` excludes identical encoded sequences before selection;
+nonidentical sequences with a zero substitution penalty remain eligible.
+
+``count_batch`` takes one list of integer thresholds per query and returns the
+number of reference scores at or below each threshold, including all ties.
+Threshold rows may be empty, unsorted or contain duplicates; output follows the
+original threshold order. Negative thresholds count zero.
+
+The paired variants score linked two-view query/reference rows and compare the
+maximum of their two lane penalties. ``paired_topk_batch`` returns
+``(ref_id, score_a, score_b)`` tuples; ``paired_count_batch`` counts all linked rows
+at each requested maximum-score threshold. Paired exact exclusion removes only
+rows identical in both views. No cross-row combinations are created.
+
+.. code-block:: python
+
+   from seqtree.gapblock import topk_batch, count_batch
+
+   queries, refs = ["ACDE", "ACD"], ["ACDE", "ACD", "ACDF"]
+   hits = topk_batch(queries, refs, k=2, threads=1)
+   thresholds = [[h.score for h in row] for row in hits]
+   counts = count_batch(queries, refs, thresholds, threads=1)
+
+Each operation runs as one native batch with the GIL released. Top-k results use
+O(Q*k) storage and each worker uses O(k + maximum sequence length) scratch.
+Count output uses O(total requested thresholds), with worker scratch proportional
+to maximum sequence length and the largest threshold row. Encoded inputs and the
+existing optional prior cube are shared. Neither path allocates a dense Q*N
+score matrix; both remain exhaustive scoring operations. They introduce no
+statistical model or calibration.
+
+``Index.ref_seqs()`` retrieves all stored strings in reference-ID order through
+one native call, preserving duplicate rows.

include/seqtree/seqtree.hpp

@@ -5,6 +5,7 @@
 #include <memory>
 #include <string>
 #include <string_view>
+#include <tuple>
 #include <vector>
 
 namespace seqtree {
@@ -221,6 +222,29 @@ std::vector<int32_t> gapblock_matrix(const std::vector<std::string>& queries,
                                      int32_t gap_extend, const std::vector<int32_t>& prior,
                                      uint32_t prior_width, int threads = 0);
 
+// 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.
+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,
+        uint32_t width, int threads, bool exclude_exact);
+std::vector<std::vector<std::tuple<uint32_t, int32_t, int32_t>>> gapblock_paired_topk_batch(
+        const std::vector<std::string>& qa, const std::vector<std::string>& qb,
+        const std::vector<std::string>& ra, const std::vector<std::string>& rb,
+        uint32_t k, 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);
+std::vector<std::vector<uint64_t>> gapblock_count_batch(const std::vector<std::string>& queries,
+        const std::vector<std::string>& refs, const std::vector<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);
+std::vector<std::vector<uint64_t>> gapblock_paired_count_batch(
+        const std::vector<std::string>& qa, const std::vector<std::string>& qb,
+        const std::vector<std::string>& ra, const std::vector<std::string>& rb,
+        const std::vector<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);
+
 // ---------------------------------------------------------------------------------------
 // Plain string edit distances on raw characters, unit costs. Deliberately alphabet-agnostic:
 // no Codec, no substitution matrix, no gap model -- the distance you want when you just need a

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "seqtree"
-version = "1.0.1"
+version = "1.0.2"
 description = "Fast fuzzy search over biological sequences (C++ core, Python bindings)"
 readme = "README.md"
 requires-python = ">=3.10"

python/seqtree/gapblock.py

@@ -42,9 +42,13 @@
 
 from ._core import Index, ScoreMatrix, SearchParams, SubstitutionMatrix
 from ._core import gapblock_matrix as _gapblock_matrix
+from ._core import gapblock_topk_batch as _gapblock_topk_batch
+from ._core import gapblock_count_batch as _gapblock_count_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
 
 __all__ = [
-    "gapblock_score", "score_matrix", "deletion_variants", "central_prior", "profile_prior",
+    "gapblock_score", "score_matrix", "topk_batch", "count_batch", "paired_topk_batch", "paired_count_batch", "deletion_variants", "central_prior", "profile_prior",
     "frame_prior", "positions_prior", "embed_in_frame", "gap_cost", "GapBlockIndex", "ScoreMatrix",
     "IslandProfile",
 ]
@@ -564,6 +568,111 @@ def score_matrix(
     return _gapblock_matrix(q, r, alphabet, matrix, gap_open, gap_extend, cube, width, threads)
 
 
+def _batch_integer(value, name, lower=0, upper=(1 << 31) - 1):
+    import operator
+    if isinstance(value, bool):
+        raise ValueError(f"{name} must be an integer")
+    try:
+        result = operator.index(value)
+    except TypeError:
+        raise ValueError(f"{name} must be an integer") from None
+    if not lower <= result <= upper:
+        raise ValueError(f"{name} must lie in [{lower}, {upper}]")
+    return result
+
+
+def _batch_options(axes, matrix, gap_open, gap_extend, gap_prior, threads):
+    go = _batch_integer(_default_gap_open(matrix) if gap_open is None else gap_open, "gap_open")
+    ge = _batch_integer(gap_extend, "gap_extend")
+    threads = _batch_integer(threads, "threads")
+    width = max((len(s) for axis in axes for s in axis), default=0)
+    cube = _prior_cube(gap_prior, width) if gap_prior is not None else []
+    return go, ge, cube, width, threads
+
+
+def _threshold_rows(thresholds, n_queries):
+    rows = [list(row) for row in thresholds]
+    if len(rows) != n_queries:
+        raise ValueError("thresholds must have one row per query")
+    return [[_batch_integer(value, "threshold", -(1 << 31)) for value in row] for row in rows]
+
+
+def topk_batch(queries, refs, k=10, matrix=None, gap_open=None, gap_extend=1,
+               gap_prior=None, alphabet="aa", threads=0, exclude_exact=False):
+    """Exhaustive single-gap-block top-k penalties, one native batch, GIL released.
+
+    Return one list of ``Hit(ref_id, score)`` per query, ordered by
+    ``(score, ref_id)``. At most k reference ROWS survive; repeated reference
+    sequences remain distinct rows. Hit edit fields are unset (zero), not an
+    alignment decomposition. Empty references yield empty query rows.
+
+    Arguments match :func:`score_matrix`, with positive integer k added. Exact
+    exclusion compares case-insensitive alphabet-encoded sequences BEFORE
+    selection; a nonidentity zero penalty remains eligible.
+
+    Shared storage holds encoded inputs and the existing optional prior cube.
+    Each worker holds O(k + max_sequence_length) scratch; results take O(Q*k),
+    never a dense Q*N matrix. This remains exhaustive O(Q*N*sequence_length)
+    scoring; no radius predicate, approximation or statistical model is added.
+    """
+    q, r = list(queries), list(refs)
+    k = _batch_integer(k, "k", 1, (1 << 32) - 1)
+    go, ge, cube, width, threads = _batch_options((q, r), matrix, gap_open, gap_extend, gap_prior, threads)
+    return _gapblock_topk_batch(q, r, k, alphabet, matrix, go, ge, cube, width, threads, exclude_exact)
+
+
+def count_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1,
+                gap_prior=None, alphabet="aa", threads=0, exclude_exact=False):
+    """Count all reference scores <= each query's requested integer thresholds.
+
+    ``thresholds`` has one row per query; rows may be empty, unsorted or contain
+    duplicates. Return integer counts in their original order, including ALL
+    score ties. Negative cutoffs count zero. Same scorer/exclusion as
+    :func:`topk_batch`; one native batch with O(Q*S) output, no dense score matrix.
+    Worker scratch is O(max_sequence_length + maximum_threshold_row_length).
+    """
+    q, r = list(queries), list(refs)
+    thresholds = _threshold_rows(thresholds, len(q))
+    go, ge, cube, width, threads = _batch_options((q, r), matrix, gap_open, gap_extend, gap_prior, threads)
+    return _gapblock_count_batch(q, r, thresholds, alphabet, matrix, go, ge, cube, width, threads, exclude_exact)
+
+
+def paired_topk_batch(queries_a, queries_b, refs_a, refs_b, k=10, matrix=None,
+                      gap_open=None, gap_extend=1, gap_prior=None, alphabet="aa",
+                      threads=0, exclude_exact=False):
+    """Top-k linked two-view rows, ordered by (max(score_a,score_b), ref_id).
+
+    Return one list of ``(ref_id, score_a, score_b)`` tuples per linked query.
+    Both query axes and both reference axes must have matching lengths.
+    Exclusion removes only rows identical in BOTH views; one-view identities
+    remain eligible. Scoring arguments and memory bounds match :func:`topk_batch`.
+    No cross-row combinations, payloads or statistical assumptions are introduced.
+    """
+    qa, qb, ra, rb = map(list, (queries_a, queries_b, refs_a, refs_b))
+    if len(qa) != len(qb) or len(ra) != len(rb):
+        raise ValueError("paired query/reference axes must agree")
+    k = _batch_integer(k, "k", 1, (1 << 32) - 1)
+    go, ge, cube, width, threads = _batch_options((qa, qb, ra, rb), matrix, gap_open, gap_extend, gap_prior, threads)
+    return _gapblock_paired_topk_batch(qa, qb, ra, rb, k, alphabet, matrix, go, ge, cube, width, threads, exclude_exact)
+
+
+def paired_count_batch(queries_a, queries_b, refs_a, refs_b, thresholds, matrix=None,
+                       gap_open=None, gap_extend=1, gap_prior=None, alphabet="aa",
+                       threads=0, exclude_exact=False):
+    """Count all linked two-view rows with max(lane penalties) <= each cutoff.
+
+    Shape, threshold order, ties and negative-cutoff semantics match
+    :func:`count_batch`; exact exclusion is the both-view predicate described in
+    :func:`paired_topk_batch`. One native batch; no dense score matrix.
+    """
+    qa, qb, ra, rb = map(list, (queries_a, queries_b, refs_a, refs_b))
+    if len(qa) != len(qb) or len(ra) != len(rb):
+        raise ValueError("paired query/reference axes must agree")
+    thresholds = _threshold_rows(thresholds, len(qa))
+    go, ge, cube, width, threads = _batch_options((qa, qb, ra, rb), matrix, gap_open, gap_extend, gap_prior, threads)
+    return _gapblock_paired_count_batch(qa, qb, ra, rb, thresholds, alphabet, matrix, go, ge, cube, width, threads, exclude_exact)
+
+
 class GapBlockIndex:
     """Search a reference set under the single-gap-block model, reusing the Hamming engine.
 

skills/seqtree/SKILL.md

@@ -29,6 +29,7 @@ only, and nothing in the result says why. Passing a matrix to `seqtrie` without
 
 ```python
 Index.build(refs: list[str], alphabet: str = "aa") -> Index      # 'aa' | 'nt' | 'iupac'
+Index.ref_seqs() -> list[str]                                  # one native call, ref_id order
 Index.search(query, params) -> list[Hit]                          # Hit(ref_id, score, n_subs, n_ins, n_dels)
 Index.search_top(query, params, k=1) -> list[Hit]
 Index.search_batch(queries, params, threads=0) -> list[list[Hit]] # releases the GIL
@@ -288,6 +289,27 @@ The result *is* the distance: the Gram transform is applied per residue when the
 so there is no `d = s(a,a) + s(b,b) - 2·s(a,b)` step, and non-negativity, symmetry and a zero
 diagonal hold by construction. Budget `4 * n * K` bytes and chunk the queries.
 
+### Bounded exhaustive reductions
+
+`gapblock.topk_batch(queries, refs, k=10, matrix=None, gap_open=None, gap_extend=1,
+gap_prior=None, alphabet="aa", threads=0, exclude_exact=False)` returns nested Hit
+rows sorted by `(score, ref_id)`, scoring every reference with the same native
+single-gap-block cell. Only score/ref_id are provided; edit fields remain unset.
+Results use O(Q*k), worker scratch O(k+maxlen); encoded inputs/prior cube shared.
+Exact encoded identity is excluded before selection; zero-cost nonidentity remains.
+
+`gapblock.count_batch(queries, refs, thresholds, **same_scoring_args)` returns
+nested integer counts for scores<=each threshold, including ALL ties. Thresholds
+have one row/query, arbitrary order/duplicates/empty rows allowed, output order
+preserved, negative cutoffs count zero. No dense Q*N allocation.
+
+`gapblock.paired_topk_batch(queries_a, queries_b, refs_a, refs_b, k=10,
+**same_scoring_args)` ranks linked rows by `(max(score_a,score_b),ref_id)` and
+returns nested `(ref_id,score_a,score_b)` tuples. Matching lane axis lengths required.
+`paired_count_batch(..., thresholds, **same_scoring_args)` counts all linked rows
+at those maximum-score cutoffs. Exclusion removes only both-view identity.
+These are generic exhaustive score reductions, not statistical calibration.
+
 ### Performance
 
 Profiled over the bundled 250k control: 91% of query time is the query-deletion-variant branch,