Prepare seqtree 1.0.2 ranked gap-block APIs

CHANGELOG.md

@@ -3,6 +3,75 @@
 All notable changes to `seqtree`. Dates are release dates. From 1.0.0 the project follows semantic
 versioning: breaking changes need a **major** bump.
 
+## [Unreleased]
+
+### 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.
+
+  ```python
+  res = ix.search_batch(peptides, max_subs=1, max_indels=1, threads=0)
+  h = res[0][0]
+  ix.ref_seq(h.ref_id)[h.offset : h.offset + h.length]   # the matched text
+  h.n_subs, h.n_ins, h.n_dels
+  ```
+
+  **The index, the build and the probe enumeration are unchanged.** With
+  `b = max_subs + max_indels + 1` disjoint blocks the total edit count cannot reach the block
+  count, so the same pigeonhole argument still leaves a **zero-error** block — and a zero-error
+  block matches the text exactly whether or not the occurrence carries gaps, so the existing
+  exact seed lookup finds it. Verified before any of this was written: 2,800 planted
+  occurrences with random insertions and deletions, **recall 1.000** from the unmodified exact
+  lookup.
+
+  What changes is verification. A seed no longer pins the start — the prefix before the block
+  may have gained or lost up to `max_indels` residues — so every start in
+  `[pos - base - max_indels, pos - base + max_indels]` is aligned by a banded DP over the states
+  `(query position, net offset, indels spent)`. An alignment must begin and end on an aligned
+  pair: a gap at either edge would stretch the reported interval over a residue that matches
+  nothing, and a leading gap is the same occurrence starting one residue over, so allowing it
+  would report one match several times.
+
+  `TextHit` gains `n_ins`, `n_dels` and **`length`** — a gapped match is no longer `len(query)`
+  residues wide, and `length` is `len(query) + n_dels - n_ins`. `TextResult.arrays()` gains the
+  three matching columns.
+
+  Two limits, both refusals rather than partial answers. Every seed block must be exact, so a
+  query must be at least `(max_subs + max_indels + 1) * k` long — at `k = 4` with
+  `max_subs=1, max_indels=1` that is 12 residues — and a shorter one raises with the bound in
+  the message. And `mismatches` is empty on this path, with `matrix=` ignored: recovering which
+  columns were substituted needs an alignment traceback the verifier does not keep.
+
+  Measured on a synthetic 5,000,000-residue corpus at `k = 4`, one thread, against the same
+  queries with `max_indels = 0`:
+
+  | query length | `max_subs` | ms/query, no indels | ms/query, `max_indels=1` |
+  |--:|--:|--:|--:|
+  | 16 | 1 | 0.001 | 0.028 |
+  | 16 | 2 | 0.001 | 0.043 |
+  | 20 | 3 | 0.002 | 0.065 |
+  | 24 | 3 | 0.002 | 0.067 |
+
+  25–37× an ungapped search at the same length and `max_subs`, and still under 0.1 ms/query.
+  An early exit when every banded state has already breached a cap is worth 2.5–3× of that; the
+  DP was otherwise walking all L rows on candidates the Hamming path rejects after two residues.
+
+  Completeness is held to the same standard as the substitution path: **set equality against an
+  independent brute force**, over `k` ∈ {3, 4, 5} × `max_subs` 0–3 × `max_indels` 1–2,
+  **6,732 occurrences, zero missing and zero extra**. `max_indels=0` is byte-identical to the
+  previous substitution path, hit for hit, and pays none of the cost.
+
 ## [1.0.1] — 2026-10-09
 
 ### Fixed

README.md

@@ -103,7 +103,10 @@ Beyond search, seqtree ships:
   completeness pinned by brute-force **set equality** over L 6–30 × `max_subs` 0–3 ×
   k ∈ {3,4,5}, and by the answer being identical across `k`. On the human proteome
   (69,578,135 residues) one thread answers a 9-mer within 2 substitutions in **1.3 ms** and a
-  15-mer within 3 in **1.5 ms**. Results come back as flat CSR arrays with zero-copy
+  15-mer within 3 in **1.5 ms**. `max_indels` adds insertions and deletions on the same index
+  and the same scheme — with `max_subs + max_indels + 1` blocks the edit count is still below
+  the block count, so pigeonhole leaves a zero-error block and the *unchanged* exact seed lookup
+  finds gapped occurrences too; only verification becomes a banded DP. Results come back as flat CSR arrays with zero-copy
   `numpy` views, mismatches as `(pos, query_aa, text_aa)` pairs, an optional fold onto
   caller-supplied group ids that makes a tie explicit, and a cap that is always reported.
 - **E-values / significance** — calibrate hit counts against a background control repertoire

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.

docs/text-index.rst

@@ -211,6 +211,81 @@ A ``.sti`` is treated as untrusted input: the header's record, bucket and postin
 cross-checked against each other and against the file length before any of them is used to size
 an allocation or to index into the mapping.
 
+Insertions and deletions
+------------------------
+
+``max_indels`` allows gaps. The two caps are independent, so ``max_subs=2, max_indels=1``
+accepts two substitutions **and** one gap — not three edits of any kind:
+
+.. code-block:: python
+
+   res = ix.search_batch(peptides, max_subs=1, max_indels=1, threads=0)
+   for h in res[0]:
+       print(h.ref_id, h.offset, h.length, h.n_subs, h.n_ins, h.n_dels)
+       print(ix.ref_seq(h.ref_id)[h.offset : h.offset + h.length])
+
+**The seed table does not change, and neither does the scheme.** With
+``b = max_subs + max_indels + 1`` disjoint blocks the total edit count cannot reach the block
+count, so the same pigeonhole argument still puts a **zero-error** block somewhere — and a
+zero-error block matches the text exactly, gaps or not, so the existing exact lookup finds it.
+Nothing about the index, the build, or the probe enumeration is different.
+
+What changes is verification. A seed no longer pins the start: the prefix before the block may
+have gained or lost up to ``max_indels`` residues, so every start in
+``[pos - base - max_indels, pos - base + max_indels]`` is aligned, by a banded DP over the
+states ``(query position, net offset, indels spent)``.
+
+Two consequences worth knowing before you call it:
+
+* **A match is no longer** ``len(query)`` **residues wide.** Read
+  :attr:`~seqtree.TextHit.length`; it is ``len(query) + n_dels - n_ins``.
+* **Every seed block must be exact**, so a query must be at least
+  ``(max_subs + max_indels + 1) * k`` long. At ``k = 4``, ``max_subs=1, max_indels=1`` needs 12
+  residues. A shorter query raises, naming the bound — it is not answered partially. Build a
+  second index at a smaller ``k`` if your queries are shorter than that.
+
+An alignment must **begin and end on an aligned pair**. A gap at either edge would stretch the
+reported interval over a residue that matches nothing, and a leading gap is just the same
+occurrence starting one residue over — allowing it would report one match two or three times.
+
+``mismatches`` is empty on this path: recovering which columns were substituted needs an
+alignment traceback the verifier does not keep. The counts are exact; the per-column detail is
+only available at ``max_indels = 0``. ``matrix=`` is likewise ignored.
+
+Cost, measured on a synthetic 5,000,000-residue corpus, ``k = 4``, one thread, against the same
+query set with ``max_indels = 0``:
+
+.. list-table::
+   :header-rows: 1
+
+   * - query length
+     - ``max_subs``
+     - ms/query, no indels
+     - ms/query, ``max_indels=1``
+   * - 16
+     - 1
+     - 0.001
+     - 0.028
+   * - 16
+     - 2
+     - 0.001
+     - 0.043
+   * - 20
+     - 3
+     - 0.002
+     - 0.065
+   * - 24
+     - 3
+     - 0.002
+     - 0.067
+
+Gapped search costs 25–37× an ungapped one at the same length and ``max_subs``: the banded DP
+replaces a byte comparison that early-exits after a residue or two, and it runs once per start
+in the window. It is still well under 0.1 ms/query. Completeness is pinned the same way the
+substitution path is — **set equality against an independent brute force**, over
+``k`` ∈ {3, 4, 5} × ``max_subs`` 0–3 × ``max_indels`` 1–2, 6,732 occurrences, zero missing and
+zero extra.
+
 Limits
 ------
 

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

include/seqtree/text_index.hpp

@@ -13,12 +13,24 @@ struct TextStore;  // internal, src/text_index.cpp
 // One row of the optional group fold: hits collapsed onto the build-time group ids.
 struct TextGroup {
     uint32_t group_id = 0;
-    uint16_t min_subs = 0;  // best (smallest) Hamming distance reaching this group
+    uint16_t min_subs = 0;  // best (smallest) substitution count reaching this group
     uint32_t n_hits   = 0;  // how many hits fell in it, at any distance
 };
 
 struct TextQueryOpts {
     uint16_t max_subs      = 0;
+    // Insertions + deletions allowed, capped independently of max_subs. 0 keeps the pure
+    // Hamming predicate and its verification path untouched.
+    //
+    // The seed table does not change: with b = max_subs + max_indels + 1 disjoint blocks the
+    // total edit count is below the block count, so by pigeonhole SOME block still matches
+    // EXACTLY -- the existing exact lookup finds every occurrence, indels included. What that
+    // costs is verification, which becomes a banded alignment over the 2*max_indels+1 possible
+    // starts rather than one Hamming scan at a fixed start.
+    //
+    // Because every block must be exact, this needs L >= (max_subs + max_indels + 1) * k;
+    // search_batch refuses a shorter query rather than answering incompletely.
+    uint16_t max_indels    = 0;
     bool     exclude_exact = false;  // drop 0-mismatch hits (a candidate that IS the text)
     bool     best_only     = false;  // keep only the minimum-n_subs shell (still ALL of it)
     bool     group_by      = false;  // also fold hits onto the build-time group ids
@@ -40,6 +52,11 @@ struct TextResult {
     std::vector<uint32_t> query_begin;  // n_queries + 1, CSR into the hit arrays
     std::vector<uint32_t> ref_id, offset;
     std::vector<uint16_t> n_subs;
+    // Indels in the alignment, and the matched substring's length in the text. With
+    // max_indels = 0 these are 0, 0 and the query length for every hit. With indels the match
+    // is no longer len(query) residues wide, so `length` is what bounds it:
+    // ref_seq(ref_id)[offset : offset + length].
+    std::vector<uint16_t> n_ins, n_dels, length;
     std::vector<int32_t>  score;
     std::vector<uint32_t> mm_begin;     // n_hits + 1, CSR into the mismatch arrays
     std::vector<uint16_t> mm_pos;

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"