Prepare seqtree 1.1.1 grouped native reductions

docs/gapblock.rst

@@ -432,3 +432,113 @@ 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.
+
+Factored paired-feature offsets
+-------------------------------
+
+``PairedFeatureOffsets(query_codes, reference_codes, costs, weight=1.0, scale=1)``
+represents an additive offset without a query-by-reference penalty table. Each
+code row contains ``[head_0, tail_0, ..., head_(D-1), tail_(D-1)]`` for a positive,
+shared depth ``D``. Integer codes index the query rows/reference columns of the
+shared nonnegative rectangular ``costs`` table. The offset is::
+
+    scale * round_even(weight * sum(max(costs[q_head_i][r_head_i],
+                                       costs[q_tail_i][r_tail_i])) / D)
+
+``round_even`` rounds to nearest integer with ties to even. Binary64 multiplication
+by ``weight`` precedes division by ``D``; ``scale`` multiplies after rounding.
+Weight must be finite and nonnegative, scale a nonnegative integer. Code rows,
+axes, costs and combined int32 score bounds are checked before worker execution.
+Both empty axes are accepted with empty code rows; nonempty axes require positive
+depth. Zero weight or scale supplies a zero offset with the same validation.
+The descriptor owns its inputs and exposes no mutable Python fields; one batch
+copy remains shared and immutable while the GIL is released.
+
+``score_matrix``, ``count_batch`` (including ``linear_mass=True``), ``ball_batch``
+and ``group_audit_batch`` accept ``feature_offsets=None``. The offset adds to any
+existing ``group_distances`` penalty. The audit component excludes both additive
+offsets. ``paired_sum_count_batch`` accepts ``feature_offsets_alpha=None`` and
+``feature_offsets_beta=None`` independently. Exact exclusion still uses full
+encoded sequence identity, including both lanes for linked rows. Omission
+preserves the existing scorer. Storage is O((Q+R)*D+C), where ``C`` is the number
+of code-pair costs; no Q-by-R offset allocation is made.
+
+Native capability can be checked by importing ``PairedFeatureOffsets`` from
+``seqtree._core`` and checking the native reducer docstring for ``feature_offsets``.
+The Python wrapper alone does not prove that an installed extension supports it.
+
+Shared batch cutoffs
+--------------------
+
+``count_batch`` (including ``linear_mass=True``) and ``paired_sum_count_batch``
+accept ``shared_thresholds=False``. With ``True``, ``thresholds`` is one flat
+integer cutoff list for the whole query batch. Python validates that list once;
+the native reducer owns and sorts one immutable cutoff vector shared by all
+workers. Outputs retain their usual query-by-cutoff shape, cutoff order,
+duplicate cutoffs, boundary ties, negative-cutoff behavior and exact exclusion.
+There is no cross-call cache. An empty cutoff list produces an empty row per
+query; empty query axes produce empty outputs. Default ``False`` preserves the
+existing per-query row API. The native entry point represents shared cutoffs as
+exactly one row, including for empty query axes, and checks that shape before
+workers. Int32 cutoff bounds and uint64 mass overflow validation remain intact.
+
+Native capability is exposed by ``shared_thresholds`` in the reducer docstrings;
+check both the Python signature and the installed native extension.
+
+Grouped label populations
+-------------------------
+
+``grouped_count_batch`` and ``grouped_paired_sum_count_batch`` prepare one immutable
+query/reference scoring population and return a native iterator over contiguous
+query chunks. The geometry is exactly the existing count/summed-pair geometry,
+including positional weights, additive group distances and factored feature
+offsets. Each query/reference key is scored once, then its single cutoff bin
+contributes to every attached target label. No statistical reducer or domain data
+is added. Default count APIs are unchanged.
+
+Required label metadata is generic integer coding:
+
+* ``reference_target_offsets`` and ``reference_target_ids`` form CSR rows for
+  the reference axis; each row contains sorted unique target IDs.
+* ``target_peptide_ids`` assigns each target to a parent label (peptide downstream).
+* ``reference_exclusive_peptide_ids`` is the parent label a reference belongs to
+  exclusively, or -1 for shared/unselected labels. Compute it from all original
+  labels, before restricting requested targets. Empty target memberships are valid;
+  an exclusive ID may be provided even when that reference has no selected target
+  membership. Included memberships must agree with any supplied exclusive ID.
+* ``query_gene_ids`` and ``reference_gene_ids`` code a separate canonical identity
+  axis for same-identity target subsets. This equality does not alter allele/group
+  costs. For linked pairs, encode a tuple of the two identities into one common ID.
+
+Thresholds are one strictly increasing int32 list, shared across the population;
+negative or empty cutoffs are supported. Every chunk owns read-only uint32 buffers,
+with no NumPy dependency or nested Python count lists. ``numpy.asarray(buffer)``
+and ``memoryview(buffer)`` expose the same allocation without copying. For ``C``
+chunk rows, ``T`` targets, ``P`` parent labels and ``K`` cutoffs:
+
+* ``chunk.start`` / ``chunk.stop`` identify its original query interval.
+* ``chunk.counts`` has shape ``(C, 1+2*T+P, K)``. Channels are whole population,
+  targets, same-identity targets, exclusive parent populations, in that order.
+* ``chunk.sizes`` has shape ``(C, 1+2*T+P)`` in the same channel order.
+* ``chunk.exact_counts`` has shape ``(C, T)`` and always reports raw full-key exact
+  multiplicity, including when exact references are excluded from other outputs.
+
+With ``exclude_exact=True``, every curve and corresponding population size is
+punctured by the full encoded query key, irrespective of additive groups/identities.
+For paired sum geometry both full sequences must match. Subtracting an exclusive
+parent channel from the whole channel therefore yields the exact complement,
+including references carrying shared labels. Count overflow is excluded explicitly
+by requiring at most uint32-max reference rows; addressable buffer sizes and all
+metadata axes/ranges are validated. Empty query axes yield no chunks; empty
+reference axes yield zero curves and sizes with ordinary shapes.
+
+``chunk_size=64`` controls output memory: O(C*(1+2*T+P)*(K+1)+C*T), plus O(Q+R)
+encoded inputs, CSR memberships and geometry tables. Worker scratch is only
+O(maximum sequence length); histogram bins live directly in disjoint chunk rows.
+Keep only the current chunk for bounded memory. A retained buffer owns its chunk
+and safely outlives the iterator. Workers use the existing CPU-budget scaffold;
+preparation and each native chunk release the GIL. No persisted computed cache,
+per-query Python search loop, global Q-by-label allocation or new distance formula
+is introduced. Native capability is ``gapblock_grouped_count_batch`` plus
+``GroupedCountBatch``/``GroupedCountChunk`` in ``seqtree._core``; the constructor
+and chunk methods document the buffer/feature support.

include/seqtree/seqtree.hpp

@@ -3,6 +3,7 @@
 
 #include <cstdint>
 #include <memory>
+#include <optional>
 #include <string>
 #include <string_view>
 #include <tuple>
@@ -228,6 +229,20 @@ std::vector<int32_t> align_dist_matrix(const std::vector<std::string>& queries,
 // Both axes are required; these rows cannot combine with length profiles.
 // Optional group IDs index a nonnegative rectangular penalty table added per pair.
 // Exact exclusion remains encoded junction identity, independent of these groups.
+// Factored nonnegative offsets: paired feature codes, one shared code cost table.
+// Round weight * sum(max(head cost, tail cost)) / depth to nearest, ties even,
+// then multiply by scale. Input descriptors are copied into immutable batch state.
+struct PairedFeatureOffsets {
+    std::vector<std::vector<int32_t>> query_codes, reference_codes, costs;
+    double weight;
+    int32_t scale;
+    PairedFeatureOffsets(std::vector<std::vector<int32_t>> q,
+                         std::vector<std::vector<int32_t>> r,
+                         std::vector<std::vector<int32_t>> c, double w = 1.0, int32_t s = 1)
+        : query_codes(std::move(q)), reference_codes(std::move(r)), costs(std::move(c)),
+          weight(w), scale(s) {}
+};
+
 std::vector<int32_t> gapblock_matrix(const std::vector<std::string>& queries,
                                      const std::vector<std::string>& refs, Alphabet,
                                      const SubstitutionMatrix* matrix, int32_t gap_open,
@@ -238,7 +253,8 @@ std::vector<int32_t> gapblock_matrix(const std::vector<std::string>& queries,
                                      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 std::vector<std::vector<int32_t>>& reference_weights = {},
+        const std::optional<PairedFeatureOffsets>& feature_offsets = std::nullopt);
 
 // Exhaustive ball rows in reference order; hit edit fields are unset (zero).
 std::vector<std::vector<Hit>> gapblock_ball_batch(const std::vector<std::string>& queries,
@@ -250,7 +266,8 @@ std::vector<std::vector<Hit>> gapblock_ball_batch(const std::vector<std::string>
         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 std::vector<std::vector<int32_t>>& reference_weights = {},
+        const std::optional<PairedFeatureOffsets>& feature_offsets = std::nullopt);
 
 // Exhaustive bounded reductions using the same gapblock cell. Top-k ties use
 // (penalty, ref_id); paired penalty is max(lane scores). Exact exclusion uses
@@ -267,7 +284,8 @@ gapblock_group_audit_batch(const std::vector<std::string>& queries,
         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 std::vector<std::vector<int32_t>>& reference_weights = {},
+        const std::optional<PairedFeatureOffsets>& feature_offsets = std::nullopt);
 
 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,
@@ -287,7 +305,8 @@ std::vector<std::vector<uint64_t>> gapblock_count_batch(const std::vector<std::s
         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 std::vector<std::vector<int32_t>>& reference_weights = {},
+        const std::optional<PairedFeatureOffsets>& feature_offsets = std::nullopt, bool shared_thresholds = false);
 std::pair<std::vector<std::vector<uint64_t>>, std::vector<std::vector<uint64_t>>> gapblock_count_mass_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,
@@ -297,7 +316,8 @@ std::pair<std::vector<std::vector<uint64_t>>, std::vector<std::vector<uint64_t>>
         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 std::vector<std::vector<int32_t>>& reference_weights = {},
+        const std::optional<PairedFeatureOffsets>& feature_offsets = std::nullopt, bool shared_thresholds = false);
 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,
@@ -319,7 +339,9 @@ std::pair<std::vector<std::vector<uint64_t>>, std::vector<std::vector<uint64_t>>
         const std::vector<std::vector<int32_t>>& query_weights_alpha = {},
         const std::vector<std::vector<int32_t>>& query_weights_beta = {},
         const std::vector<std::vector<int32_t>>& reference_weights_alpha = {},
-        const std::vector<std::vector<int32_t>>& reference_weights_beta = {});
+        const std::vector<std::vector<int32_t>>& reference_weights_beta = {},
+        const std::optional<PairedFeatureOffsets>& feature_offsets_alpha = std::nullopt,
+        const std::optional<PairedFeatureOffsets>& feature_offsets_beta = std::nullopt, bool shared_thresholds = false);
 
 // ---------------------------------------------------------------------------------------
 // Plain string edit distances on raw characters, unit costs. Deliberately alphabet-agnostic:
@@ -341,4 +363,48 @@ std::vector<int32_t> hamming_matrix(const std::vector<std::string>& a,
 std::vector<int32_t> levenshtein_matrix(const std::vector<std::string>& a,
                                         const std::vector<std::string>& b, int threads = 0);
 
+// Generic label histograms over one prepared native query/reference population.
+struct GroupedGeometry {
+    std::vector<PairedFeatureOffsets> additional_features;
+    std::optional<std::vector<int32_t>> gap_prior;
+    std::string alignment_mode = "gapblock";
+    int32_t alignment_scale = 1;
+    uint32_t trim_left = 0, trim_right = 0;
+    int32_t score_quantum = 1;
+    std::vector<std::vector<int32_t>> weights, group_distances, query_weights, reference_weights;
+    std::vector<int32_t> query_groups, reference_groups;
+    std::optional<PairedFeatureOffsets> features;
+};
+struct CountBuffer {
+    std::vector<uint32_t> data;
+    std::vector<size_t> shape;
+};
+struct GroupedCountChunk {
+    size_t start = 0, stop = 0;
+    CountBuffer counts, sizes, exact_counts;
+};
+class GroupedCountBatch {
+public:
+    GroupedCountBatch(const std::vector<std::string>& qa, const std::vector<std::string>& ra,
+        std::vector<int32_t> thresholds, std::vector<uint32_t> target_offsets,
+        std::vector<uint32_t> target_ids, std::vector<int32_t> target_peptides,
+        std::vector<int32_t> exclusive_peptides, std::vector<int32_t> query_genes,
+        std::vector<int32_t> reference_genes, size_t chunk_size, Alphabet alphabet,
+        const SubstitutionMatrix* matrix, int32_t go, int32_t ge, std::vector<int32_t> prior,
+        uint32_t width, int threads, bool exclude_exact, const GroupedGeometry& alpha,
+        const std::optional<GroupedGeometry>& beta = std::nullopt,
+        const std::vector<std::string>& qb = {}, const std::vector<std::string>& rb = {},
+        std::optional<std::vector<int32_t>> query_withheld_peptides = std::nullopt,
+        std::vector<uint32_t> reference_peptide_offsets = {},
+        std::vector<int32_t> reference_peptide_ids = {});
+    std::optional<GroupedCountChunk> next();
+    size_t query_count() const;
+    size_t reference_count() const;
+    size_t target_count() const;
+    size_t peptide_count() const;
+private:
+    struct Impl;
+    std::shared_ptr<Impl> impl_;
+};
+
 }  // namespace seqtree

pyproject.toml

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