Merge pull request 'ci: disable numa feature, bump obikmer, and document Sankoff costs' (#65) from push-oruynkvporsn into main

Reviewed-on: #65

.gitea/workflows/ci.yml

@@ -1,4 +1,4 @@
-name: CI
+pname: CI
 
 on:
   pull_request:
@@ -25,11 +25,19 @@ jobs:
             ~/.cargo/registry
             ~/.cargo/git
             src/target
-          key: ${{ runner.os }}-cargo-${{ hashFiles('src/Cargo.lock') }}
-          restore-keys: ${{ runner.os }}-cargo-
+          key: ${{ runner.os }}-cargo-v2-${{ hashFiles('src/Cargo.lock') }}
+          restore-keys: ${{ runner.os }}-cargo-v2-
 
+      # Both `obikmer` and `obikindex` default to the `numa` feature
+      # (hwloc-based topology detection + CPU pinning), which is only useful
+      # on bare-metal multi-socket indexing hosts. Under this runner's
+      # container/cgroup setup it deadlocks at startup — confirmed live
+      # (2026-08-11): the same test binary hangs indefinitely with `numa` on
+      # and passes instantly, repeatedly, with it off, on the same
+      # container. Disable it for CI; it has nothing to do with test
+      # correctness.
       - name: Build
-        run: cargo build --release
+        run: cargo build --release --no-default-features
 
       - name: Test
-        run: cargo test --release
+        run: cargo test --release --no-default-features

.gitignore

@@ -9,6 +9,7 @@ data-stress
 ./**/*.json
 *.bin
 *.log
+*.csv
 Betula_exilis--IGA-24-33
 benchmark/genomes
 benchmark/simulated_data

docmd/implementation/filtering.md

@@ -92,18 +92,48 @@ For each genome:
 
 | Flag | Applies to | Meaning |
 |------|-----------|---------|
-| `--min-count N` | ingroup | k-mer present in at least N ingroup genomes |
-| `--max-count N` | ingroup | k-mer present in at most N ingroup genomes |
+| `--min-count N` | ingroup | k-mer present in at least N ingroup genomes (N may be negative, see below) |
+| `--max-count N` | ingroup | k-mer present in at most N ingroup genomes (N may be negative, see below) |
 | `--min-frac F` | ingroup | k-mer present in at least fraction F of ingroup genomes |
 | `--max-frac F` | ingroup | k-mer present in at most fraction F of ingroup genomes |
-| `--min-outgroup-count N` | outgroup | k-mer present in at least N outgroup genomes |
-| `--max-outgroup-count N` | outgroup | k-mer present in at most N outgroup genomes |
+| `--min-outgroup-count N` | outgroup | k-mer present in at least N outgroup genomes (N may be negative, see below) |
+| `--max-outgroup-count N` | outgroup | k-mer present in at most N outgroup genomes (N may be negative, see below) |
 | `--min-outgroup-frac F` | outgroup | k-mer present in at least fraction F of outgroup genomes |
 | `--max-outgroup-frac F` | outgroup | k-mer present in at most fraction F of outgroup genomes |
 | `--min-total-count N` | all genomes | sum of per-genome counts ≥ N (`filter` only) |
 | `--max-total-count N` | all genomes | sum of per-genome counts ≤ N (`filter` only) |
 | `--presence-threshold N` | all | per-genome count > N to be considered "present" (default 0) |
 
+### Negative counts — offset from group size
+
+The four integer count flags (`--min-count`, `--max-count`, `--min-outgroup-count`,
+`--max-outgroup-count`) accept **negative** values, interpreted as an offset counted
+down from the group size `n`, resolved at run time once `n` is known:
+
+| Value | Effective threshold |
+|-------|---------------------|
+| `N ≥ 0` | literal absolute count `N` |
+| `-x` (x > 0) | `max(1, n − x)` — "all but x" |
+
+`-1` literally means *all but one*, `-2` *all but two*, and so on. This expresses
+a quorum relative to the group size that a plain fraction cannot state exactly
+(e.g. "present in every genome except at most one" is `n−1`, which is `0.9` for
+`n = 10` but `0.857…` for `n = 7`).
+
+The threshold is **floored at 1**, never 0: the negative form always keeps
+constraining the group. Without the floor, `--min-count -1` on a singleton
+ingroup (`n = 1`) would resolve to `0` ("at least 0") and silently drop the
+constraint; the floor makes it `1` ("present in that one genome") instead.
+
+To express a count of `0` (e.g. "absent from the ingroup"), use the literal `0`,
+not a negative — `0` and `-0` are indistinguishable, so the offset form starts at
+`-1`.
+
+> **Edge case** — on an *empty* group (`n = 0`, e.g. a predicate matching no
+> genome), a negative count still resolves to `1`, an impossible constraint that
+> rejects every k-mer. This is consistent with an empty group letting nothing
+> through, but differs from the "no constraint" behaviour of the fraction flags.
+
 **Conditional defaults** — the defaults for `--min-frac` and `--max-outgroup-count` depend on two conditions:
 whether the corresponding group was declared, **and** whether any quorum flag for that group was explicitly set.
 
@@ -215,6 +245,17 @@ obikmer filter src --output dst \
   --max-outgroup-count 0
 ```
 
+Noise-tolerant core — keep k-mers present in *all but one* ingroup genome
+(`-1` = `n−1`) and absent from *all but one* of the outgroup:
+
+```sh
+obikmer filter src --output dst \
+  --ingroup  "genus=Betula" \
+  --outgroup "*" \
+  --min-count -1 \
+  --max-outgroup-count -1
+```
+
 To dump only k-mers specific to *Betula nana*:
 
 ```sh

docmd/implementation/obicompactvec.md

@@ -347,11 +347,24 @@ Provided finalisations:
 | `relfreq_euclidean_dist_matrix()` | `√partial_relfreq_euclidean[i,j]` |
 | `hellinger_dist_matrix()` | `√partial_hellinger[i,j] / √2` |
 | `hellinger_euclidean_dist_matrix()` | `√partial_hellinger[i,j]` |
+| `threshold_mash_dist_matrix(k, t)` | Mash distance, derived from `threshold_jaccard_dist_matrix(t)` — no separate partial |
 
 ### BitPartials
 
 Required: `partial_jaccard() -> (Array2<u64>, Array2<u64>)`, `partial_hamming() -> Array2<u64>`. Both additive across layers and partitions.
 
+Provided finalisations also include `jaccard_dist_matrix()`, `hamming_dist_matrix()`, and `mash_dist_matrix(k)`.
+
+### Mash distance
+
+`mash_dist_matrix`/`threshold_mash_dist_matrix` add no new additive primitive: both are a pointwise transform of the existing Jaccard distance matrix, per the Mash mutation-rate estimator [@Mash-distances-doc; @Fan2015-mash-formula]:
+
+```
+D = -1/k · ln(2J / (1+J)),  J = 1 - d_jaccard
+```
+
+`J ≤ 0` (i.e. `d_jaccard ≥ 1`, no shared k-mers) maps to `D = 1` (maximal distance) rather than the `ln` singularity at `J = 0`.
+
 ---
 
 ## Temp-file-backed types

docmd/index.md

@@ -13,7 +13,7 @@
 | `query`     | Query an index with sequences and annotate matches |
 | `dump`      | Dump all indexed k-mers as CSV (kmer + per-genome counts or presence); supports the shared [kmer filtering](implementation/filtering.md) system; `--head N` limits output to the first N k-mers |
 | `annotate`  | Add or update genome metadata from a CSV file; or dump metadata as CSV |
-| `distance`  | Compute pairwise distance matrix between genomes; optionally build NJ/UPGMA trees; `--presence-threshold N` sets the minimum count to consider a k-mer present when computing Jaccard on count indexes (default 1) |
+| `distance`  | Compute pairwise distance matrix between genomes (`--metric jaccard\|mash\|hamming\|bray-curtis\|relfreq-bray-curtis\|euclidean\|relfreq-euclidean\|hellinger\|hellinger-euclidean`); optionally build NJ/UPGMA trees; `--presence-threshold N` sets the minimum count to consider a k-mer present when computing Jaccard/Mash on count indexes (default 1) |
 | `unitig`    | Build a global de Bruijn graph across all partitions and enumerate its unitigs as FASTA; supports the shared [kmer filtering](implementation/filtering.md) system |
 | `select`    | Project and/or aggregate genome columns into a new or in-place index; the column-axis counterpart of `filter` (see [select](implementation/select.md)) |
 | `estimate`  | Estimate approximate-index parameters (z, evidence bits, FP rates) before indexing |

docmd/references.bib

@@ -241,3 +241,21 @@
 	volume = 33,
 	year = 2017,
 	bdsk-url-1 = {http://dx.doi.org/10.1093/bioinformatics/btw832}}
+
+@misc{Mash-distances-doc,
+	author = {{Marbl Lab}},
+	howpublished = {Mash documentation},
+	title = {Mash Distance},
+	url = {https://mash.readthedocs.io/en/latest/distances.html},
+	urldate = {2026-07-09},
+	year = 2026}
+
+@article{Fan2015-mash-formula,
+	author = {Fan, Huan and Ives, Anthony R and Surget-Groba, Yann and Cannon, Charles H},
+	doi = {10.1186/s12864-015-1647-5},
+	journal = {BMC Genomics},
+	number = 1,
+	title = {An assembly and alignment-free method of phylogeny reconstruction from next-generation sequencing data},
+	url = {https://doi.org/10.1186/s12864-015-1647-5},
+	volume = 16,
+	year = 2015}

docmd/theory/entropy.md

@@ -1,6 +1,6 @@
 # Kmer entropy filter
 
-Low-complexity kmers (polyA, polyT, tandem repeats) are detected and excluded during phase 1. The filter computes a **normalized Shannon entropy** over sub-words of multiple sizes, corrected for two sources of bias: the small number of observations within a single kmer, and the unequal sizes of circular equivalence classes.
+Low-complexity kmers (polyA, polyT, tandem repeats) are detected and excluded during phase 1. The filter computes a **normalized Shannon entropy** over sub-words of multiple sizes, corrected for one source of bias: the small number of observations within a single kmer relative to the number of possible sub-words.
 
 ## Sub-word frequencies
 
@@ -8,17 +8,15 @@ For a kmer of length k and a sub-word size ws (1 ≤ ws ≤ ws_max, typically ws
 
 $$w_i = \text{kmer}[i \mathinner{..} i+ws-1], \quad i = 0, \ldots, n_{\text{words}}-1$$
 
-Each sub-word is mapped to its **circular canonical form**: the lexicographic minimum among all cyclic rotations of the word **and all cyclic rotations of its reverse complement**. This extended equivalence relation ensures that entropy(K) = entropy(revcomp(K)) — the filter is strand-symmetric. Let $s_j$ be the size of equivalence class $j$ (number of distinct raw words mapping to canonical form $j$), and $f_j$ the count of canonical form $j$ among the $n_{\text{words}}$ sub-words ($\sum_j f_j = n_{\text{words}}$).
+Each sub-word is tallied under its own raw 2-bit-packed value — **no canonicalization**. Let $f_j$ be the count of raw word $j$ among the $n_{\text{words}}$ sub-words ($\sum_j f_j = n_{\text{words}}$), over the $4^{ws}$ possible raw words.
+
+An earlier version of this filter first folded each sub-word into a circular+reverse-complement equivalence class, then "unfolded" the observed class frequency back onto its members to correct for unequal class sizes. That machinery bought nothing it was claimed for — see *Why no equivalence classes* below — while measurably weakening detection of the very sequences the filter exists to catch, so it was removed.
 
 ## Corrected Shannon entropy
 
-The circular equivalence classes have unequal sizes: under a uniform distribution over all $4^{ws}$ raw words, class $j$ is visited with probability $s_j / 4^{ws}$, not $1/n_a$. Computing entropy directly over canonical classes therefore underestimates the entropy of a random sequence.
+$$H_{\text{corr}} = \log(n_{\text{words}}) - \frac{1}{n_{\text{words}}} \sum_j f_j \log f_j$$
 
-The correction "unfolds" each canonical class back to its member raw words, redistributing each observation of class $j$ equally among its $s_j$ members:
-
-$$H_{\text{corr}} = \log(n_{\text{words}}) - \frac{1}{n_{\text{words}}} \sum_j f_j \log f_j + \frac{1}{n_{\text{words}}} \sum_j f_j \log s_j$$
-
-The last term is the correction for unequal class sizes. For a uniformly random sequence ($f_j \approx n_{\text{words}} \cdot s_j / 4^{ws}$), this gives $H_{\text{corr}} \approx \log(4^{ws}) = 2 \cdot ws \cdot \log 2$, the maximum entropy over raw words.
+This is a plain Shannon entropy over the observed raw-word frequencies.
 
 ## Maximum entropy correction for small samples
 
@@ -42,27 +40,45 @@ $$\text{entropy}(kmer) = \min_{ws=1}^{ws_{\max}} \hat{H}(ws)$$
 
 A value near 0 indicates low complexity (e.g. AAAA…); near 1 indicates high complexity. A kmer is rejected if $\text{entropy}(kmer) < \theta$, where $\theta$ is a collection parameter (default 0.7). The minimum across word sizes ensures that any scale of repetition is detected independently: polyA is caught at ws=1, dinucleotide repeats at ws=2, etc.
 
+## Why no equivalence classes
+
+A prior design folded each sub-word into the canonical form of its circular-rotation + reverse-complement equivalence class before tallying, on the reasoning that (a) it guarantees $\text{entropy}(K) = \text{entropy}(\text{revcomp}(K))$, and (b) collapsing phase-shifted repeats (e.g. `ATG` ≡ `TGA` ≡ `GAT`) into one class better reflects that they are "the same" low-complexity pattern.
+
+Both properties already hold for the raw, unfolded entropy above, without any class machinery:
+
+- **Reverse complement**: for any K of length n, window $j$ of $\text{revcomp}(K)$ equals $\text{revcomp}$ of window $(n{-}ws{-}j)$ of K. This is a bijection between the window sets under which each window maps to its own revcomp — and revcomp is itself a bijection (involution) on the space of raw ws-mers. So the multiset of raw-word frequencies for $\text{revcomp}(K)$ is exactly a relabeling of the multiset for K, and Shannon entropy — a function of the frequency multiset alone — is exactly invariant. No folding required, for any K.
+- **Tandem repeats**: a period-p repeat sampled by a stride-1 sliding window naturally cycles through its own rotations as raw tokens (e.g. `ATGATGATG…` yields the raw words `ATG`, `TGA`, `GAT` in rotation as the window slides). The low diversity this represents (few distinct raw words out of $4^{ws}$ possible) is already visible in the raw frequency distribution — no folding needed to detect it.
+
+What the fold-then-unfold step actually did was credit each observed class with the frequency of equivalence-class members that were **never observed on the read strand**, inflating $H_{\text{corr}}$ for genuine repeats. Worked example: k=31, ws=3, kmer = `ATG` repeated ($n_{\text{words}}=29$, all 29 windows fall into one class of size 6 under the old scheme — 3 rotations × forward/revcomp):
+
+| | $H_{\text{corr}}$ | normalized |
+|---|---|---|
+| old (folded, class size 6) | $\log 6 \approx 1.79$ | $\approx 0.53$ |
+| current (raw, unfolded) | $\log 3 \approx 1.10$ | $\approx 0.33$ |
+
+The gap is not a rounding artifact: per sub-word order, the folded score for this same repeat swings from 0.53 (ws=3, aligned with the period) up to **1.03** (ws=5, misaligned with the period) — i.e. a period-3 repeat could score *above* the theoretical maximum for a random sequence, depending on which ws happens to divide the repeat's period. The raw formula stays flat at ≈0.33–0.40 across ws=2..6 regardless of alignment, which is the robustness the "minimum across ws" design was meant to provide in the first place.
+
 ## Interpretation as an effective number of classes
 
-$H_{\text{corr}}$ is a standard Shannon entropy over raw words (after unfolding the equivalence classes), so the classical perplexity interpretation holds directly: $N_{\text{eff}} = e^{H_{\text{corr}}}$ is the number of equiprobable classes that would yield the same entropy.
+$H_{\text{corr}}$ is a standard Shannon entropy over raw words, so the classical perplexity interpretation holds directly: $N_{\text{eff}} = e^{H_{\text{corr}}}$ is the number of equiprobable raw words that would yield the same entropy.
 
-For the normalised score $\hat{H}$, dividing by $H_{\text{max}}$ changes the logarithm base:
+For the normalised score $\hat{H}$, dividing by $H_{\max}$ changes the logarithm base:
 
-$$\hat{H} = \frac{\log N_{\text{eff}}}{\log N_{\text{max}}} = \log_{N_{\text{max}}} N_{\text{eff}} \quad \Longleftrightarrow \quad N_{\text{eff}} = N_{\text{max}}^{\,\hat{H}}$$
+$$\hat{H} = \frac{\log N_{\text{eff}}}{\log N_{\max}} = \log_{N_{\max}} N_{\text{eff}} \quad \Longleftrightarrow \quad N_{\text{eff}} = N_{\max}^{\,\hat{H}}$$
 
-The property is preserved: $\hat{H}$ is the logarithm (in base $N_{\text{max}}$) of the effective number of equi-represented classes.
+The property is preserved: $\hat{H}$ is the logarithm (in base $N_{\max}$) of the effective number of equi-represented raw words.
 
-In the large-sample limit ($n_{\text{words}} \gg 4^{ws}$), $N_{\text{max}} \approx 4^{ws}$, giving:
+In the large-sample limit ($n_{\text{words}} \gg 4^{ws}$), $N_{\max} \approx 4^{ws}$, giving:
 
 $$N_{\text{eff}} \approx 4^{ws \cdot \hat{H}}$$
 
-This has a clean interpretation: $ws \cdot \hat{H}$ is the **effective word length** (in bases) of a perfectly uniform distribution that would produce the same entropy. At $\hat{H} = 1$ the full space of $4^{ws}$ words is used; at $\hat{H} = 0.5$ with ws=2, only $4^1 = 4$ effective classes out of 16 are occupied.
+This has a clean interpretation: $ws \cdot \hat{H}$ is the **effective word length** (in bases) of a perfectly uniform distribution that would produce the same entropy. At $\hat{H} = 1$ the full space of $4^{ws}$ words is used; at $\hat{H} = 0.5$ with ws=2, only $4^1 = 4$ effective words out of 16 are occupied.
 
-In our actual regime, $n_{\text{words}}$ is small and $4^{ws}$ can exceed $n_{\text{words}}$, so $H_{\text{max}} < \log(4^{ws})$ due to the small-sample correction. The exact effective count is $N_{\text{max}}^{\hat{H}}$, not $4^{ws \cdot \hat{H}}$.
+In our actual regime, $n_{\text{words}}$ is small and $4^{ws}$ can exceed $n_{\text{words}}$, so $H_{\max} < \log(4^{ws})$ due to the small-sample correction. The exact effective count is $N_{\max}^{\hat{H}}$, not $4^{ws \cdot \hat{H}}$.
 
 ## Properties
 
 The entropy score is a function of the kmer sequence alone — it does not depend on the surrounding context or on the position within any genome. Two consequences:
 
-- **Orientation invariance**: $\text{entropy}(K) = \text{entropy}(\text{revcomp}(K))$, guaranteed by the strand-symmetric canonical form.
+- **Orientation invariance**: $\text{entropy}(K) = \text{entropy}(\text{revcomp}(K))$ — see *Why no equivalence classes* above for why this holds without any explicit strand-folding step.
 - **Context independence**: the same kmer is always rejected or always kept, regardless of which genome it occurs in, where in that genome it appears, or which strand is considered. The filter defines a fixed partition of the kmer space into low-complexity and valid kmers.

docmd/theory/entropy.refs.md

@@ -3,10 +3,14 @@
 
 ## Code couvert
 
-- `obiskbuilder/src/entropy_table.rs` — filtre Shannon sur les kmers à basse complexité
-- `obiskbuilder/src/lib.rs` — application du filtre lors du scatter (phase 1)
+- `obikentropy/src/table.rs`, `obikentropy/src/tracker.rs` — formule d'entropie et tables de correction petits effectifs
+- `obikentropy/src/kmer_entropy.rs` — entropie d'un kmer isolé (`KmerEntropy`)
+- `obiskbuilder/src/rolling_stat.rs` — composition de `obikentropy::EntropyTracker` dans le suivi streaming (sélection de minimiseur + entropie)
+- `obiskbuilder/src/iter.rs`, `obiskbuilder/src/stream_iter.rs` — application du filtre lors du scatter (phase 1)
 
 ## Notes
 
-Document théorique stable. Vérifier que les paramètres `theta` et `level_max` dans le CLI
+Le repli en classes d'équivalence circulaires + brin inverse (décrit dans une version antérieure de ce document) a été supprimé : voir la section « Why no equivalence classes » de `entropy.md` pour la justification théorique et numérique.
+
+Vérifier que les paramètres `theta` et `level_max` dans le CLI
 (`obikmer/src/cli.rs` → `CommonArgs`) correspondent bien à ce qui est décrit.

mkdocs.yml

@@ -36,6 +36,7 @@ nav:
       - Entropy filter: theory/entropy.md
       - Minimizer selection: theory/minimizer.md
       - Partitioning architecture: theory/indexing.md
+      - Central-position SNP distance (discussion): theory/evolutionary_distances.md
   - Implementation:
       - SuperKmer: implementation/superkmer.md
       - Kmer: implementation/kmer.md

src/Cargo.lock

@@ -1682,6 +1682,13 @@ dependencies = [
  "xxhash-rust",
 ]
 
+[[package]]
+name = "obikentropy"
+version = "0.1.0"
+dependencies = [
+ "obikseq",
+]
+
 [[package]]
 name = "obikindex"
 version = "0.1.0"
@@ -1694,17 +1701,21 @@ dependencies = [
  "obikpartitionner",
  "obikseq",
  "obilayeredmap",
+ "obipipeline",
+ "obiread",
+ "obiskbuilder",
  "obiskio",
  "obisys",
  "rayon",
  "serde",
  "serde_json",
+ "tempfile",
  "tracing",
 ]
 
 [[package]]
 name = "obikmer"
-version = "1.1.37"
+version = "1.1.43"
 dependencies = [
  "clap",
  "csv",
@@ -1742,6 +1753,7 @@ dependencies = [
  "niffler 3.0.0",
  "obicompactvec",
  "obidebruinj",
+ "obikentropy",
  "obikrope",
  "obikseq",
  "obilayeredmap",
@@ -1824,7 +1836,9 @@ dependencies = [
 name = "obiskbuilder"
 version = "0.1.0"
 dependencies = [
+ "criterion2",
  "lazy_static",
+ "obikentropy",
  "obikrope",
  "obikseq",
  "obiread",

src/Cargo.toml

@@ -1,5 +1,5 @@
 [workspace]
 resolver = "3"
-members = ["obikseq", "obiread", "obiskbuilder", "obifastwrite", "obikmer","obikrope","obipipeline", "obikpartitionner","obiskio","obidebruinj","obilayeredmap", "obicompactvec", "obisys", "obikindex", "obitaxonomy"]
+members = ["obikseq", "obiread", "obiskbuilder", "obifastwrite", "obikmer","obikrope","obipipeline", "obikpartitionner","obiskio","obidebruinj","obilayeredmap", "obicompactvec", "obisys", "obikindex", "obitaxonomy", "obikentropy"]
 [profile.release]
 debug = 1

src/obicompactvec/src/lib.rs

@@ -7,6 +7,7 @@ mod intmatrix;
 mod layer_meta;
 mod meta;
 mod reader;
+mod siblingannex;
 mod tempbitvec;
 mod tempintvec;
 mod views;
@@ -18,6 +19,7 @@ pub use builder::PersistentCompactIntVecBuilder;
 pub use colgroup::{ColGroup, FilterMask, MatrixGroupOps, eval_filter_mask};
 pub use intmatrix::{PersistentCompactIntMatrix, PersistentCompactIntMatrixBuilder, pack_compact_int_matrix};
 pub use layer_meta::LayerMeta;
+pub use siblingannex::{FamilyMask, SiblingAnnex, SiblingAnnexBuilder};
 pub use reader::{PersistentCompactIntVec, Iter as CompactIntVecIter};
 pub use tempbitvec::{TempBitVec, TempBitVecBuilder};
 pub use tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder};

src/obicompactvec/src/siblingannex.rs

@@ -0,0 +1,245 @@
+//! Family presence-mask annex: a compact, read-only-after-build, per-slot
+//! derived value used by the central-position SNP distance estimator (see
+//! `docmd/theory/evolutionary_distances.md`, "Step 2b" and "Definitions:
+//! family, and the canonical form of a family").
+//!
+//! One byte is stored per MPHF slot of a partition/layer, its low 4 bits
+//! encoding a **presence mask** for the slot's k-mer's "family" (the up to 4
+//! k-mers sharing the same flanks, differing only at the central base):
+//! bit `b` (`b` = 0..3, in the fixed A/C/G/T = 0/1/2/3 encoding already used
+//! for a single nucleotide) is set iff the family member whose *own* central
+//! base — in its own canonical orientation — is `b`, is observed anywhere in
+//! the current multi-genome index. This is a property of the whole index,
+//! not of any one genome.
+//!
+//! Both facts the earlier (superseded) 3-bit design stored explicitly are
+//! derived from the mask instead, not stored:
+//! - sibling count = `popcount(mask) - 1`;
+//! - minorant = regenerate the family's 4 canonical forms from the slot's
+//!   own k-mer (`CanonicalKmerOf::central_canonical_neighbors`, cheap, no
+//!   lookup), compare the raw encodings of whichever are set in the mask,
+//!   take the smallest — see `obikindex::siblings`.
+//!
+//! Mask value 0 is logically unreachable as a real result (a slot's own base
+//! is always present in its own family) and is reused as the "not yet
+//! computed" sentinel: annex files are pre-initialised to all-zero, and a
+//! real value is only ever written once, by the computation pass.
+//!
+//! Deliberately simpler than a true 4-bit pack (1 byte/slot instead of 4
+//! bits/slot): correctness and simplicity first, for a first implementation.
+//! Packing to 4 bits/slot is a pure storage-density follow-up, not a
+//! behavioural change, left for later.
+
+use std::fs::{File, OpenOptions};
+use std::io;
+use std::path::{Path, PathBuf};
+
+use memmap2::{Mmap, MmapMut};
+
+const MAGIC: [u8; 4] = *b"PSIB";
+
+// Header: magic(4) + _pad(4) + n(8) = 16 bytes. Data (1 byte/slot) follows.
+const HEADER_SIZE: usize = 16;
+
+/// A family presence mask: bit `b` set iff the member whose own canonical
+/// central base is `b` (0=A, 1=C, 2=G, 3=T) is observed in the index.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub struct FamilyMask(u8);
+
+impl FamilyMask {
+    /// The empty mask — never a valid *computed* result (a slot's own base
+    /// is always present in its own family) — used only to build up a mask
+    /// via repeated [`with`](Self::with) calls before storing it.
+    pub const EMPTY: FamilyMask = FamilyMask(0);
+
+    /// Set bit `base` (0=A, 1=C, 2=G, 3=T).
+    #[inline]
+    pub fn with(self, base: u8) -> Self {
+        debug_assert!(base < 4, "base out of range: {base}");
+        FamilyMask(self.0 | (1 << base))
+    }
+
+    /// Is the member with central base `base` (0..3) present?
+    #[inline]
+    pub fn has(self, base: u8) -> bool {
+        debug_assert!(base < 4, "base out of range: {base}");
+        self.0 & (1 << base) != 0
+    }
+
+    /// Number of family members observed anywhere in the index (1..=4).
+    #[inline]
+    pub fn family_size(self) -> u32 {
+        self.0.count_ones()
+    }
+
+    /// Number of *other* members observed (0..=3) — `family_size() - 1`.
+    #[inline]
+    pub fn siblings(self) -> u32 {
+        self.family_size() - 1
+    }
+
+    /// Raw bitmask (bit `b` = base `b` present) — for callers that build up
+    /// a mask via their own bit operations (e.g. concurrently, via an
+    /// `AtomicU8`) and only need the `FamilyMask` wrapper at the end.
+    #[inline]
+    pub fn bits(self) -> u8 {
+        self.0
+    }
+
+    /// Construct from a raw bitmask (only the low 4 bits are kept).
+    #[inline]
+    pub fn from_bits(bits: u8) -> Self {
+        FamilyMask(bits & 0b1111)
+    }
+
+    #[inline]
+    fn encode(self) -> u8 {
+        self.0
+    }
+
+    #[inline]
+    fn decode(byte: u8) -> Option<Self> {
+        if byte == 0 {
+            // Unreachable for a real result — reserved as the "not yet
+            // computed" sentinel.
+            return None;
+        }
+        Some(FamilyMask(byte & 0b1111))
+    }
+}
+
+// ── SiblingAnnex (reader) ───────────────────────────────────────────────────
+
+pub struct SiblingAnnex {
+    mmap: Mmap,
+    n: usize,
+    path: PathBuf,
+}
+
+impl SiblingAnnex {
+    pub fn open(path: &Path) -> io::Result<Self> {
+        let mmap = unsafe { Mmap::map(&File::open(path)?)? };
+        if mmap.len() < HEADER_SIZE {
+            return Err(io::Error::new(io::ErrorKind::InvalidData, "PSIB file too short"));
+        }
+        if mmap[0..4] != MAGIC {
+            return Err(io::Error::new(io::ErrorKind::InvalidData, "bad PSIB magic"));
+        }
+        let n = u64::from_le_bytes(mmap[8..16].try_into().unwrap()) as usize;
+        if mmap.len() < HEADER_SIZE + n {
+            return Err(io::Error::new(io::ErrorKind::InvalidData, "PSIB file truncated"));
+        }
+        Ok(Self { mmap, n, path: path.to_path_buf() })
+    }
+
+    pub fn path(&self) -> &Path { &self.path }
+    pub fn len(&self) -> usize { self.n }
+    pub fn is_empty(&self) -> bool { self.n == 0 }
+
+    /// `None` means the slot has not (yet) been computed — see module docs.
+    pub fn get(&self, slot: usize) -> Option<FamilyMask> {
+        FamilyMask::decode(self.mmap[HEADER_SIZE + slot])
+    }
+}
+
+// ── SiblingAnnexBuilder (writer) ────────────────────────────────────────────
+
+pub struct SiblingAnnexBuilder {
+    mmap: MmapMut,
+    n: usize,
+    path: PathBuf,
+}
+
+impl SiblingAnnexBuilder {
+    /// Create a new annex of `n` slots at `path`, pre-initialised to the
+    /// "not yet computed" sentinel (all-zero).
+    pub fn new(n: usize, path: &Path) -> io::Result<Self> {
+        let file_size = HEADER_SIZE + n;
+        let file = OpenOptions::new()
+            .read(true).write(true).create(true).truncate(true)
+            .open(path)?;
+        file.set_len(file_size as u64)?;
+        let mut mmap = unsafe { MmapMut::map_mut(&file)? };
+        mmap[0..4].copy_from_slice(&MAGIC);
+        mmap[4..8].copy_from_slice(&[0u8; 4]);
+        mmap[8..16].copy_from_slice(&(n as u64).to_le_bytes());
+        // Data region left at 0 by `set_len`/mmap — the sentinel value.
+        Ok(Self { mmap, n, path: path.to_path_buf() })
+    }
+
+    pub fn len(&self) -> usize { self.n }
+    pub fn is_empty(&self) -> bool { self.n == 0 }
+
+    pub fn get(&self, slot: usize) -> Option<FamilyMask> {
+        FamilyMask::decode(self.mmap[HEADER_SIZE + slot])
+    }
+
+    pub fn set(&mut self, slot: usize, mask: FamilyMask) {
+        // Redundant concurrent writes from independent recomputation paths
+        // converge to the same encoded byte for a given slot, so a plain
+        // store here is safe even without external synchronisation, as long
+        // as the byte write itself is atomic (true for a single aligned
+        // byte on every platform this project targets).
+        self.mmap[HEADER_SIZE + slot] = mask.encode();
+    }
+
+    pub fn close(self) -> io::Result<()> { self.mmap.flush() }
+
+    pub fn finish(self) -> io::Result<SiblingAnnex> {
+        let path = self.path.clone();
+        self.close()?;
+        SiblingAnnex::open(&path)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use tempfile::tempdir;
+
+    #[test]
+    fn sentinel_is_zero_and_unset_slots_read_as_uncomputed() {
+        let dir = tempdir().unwrap();
+        let path = dir.path().join("test.psib");
+        let builder = SiblingAnnexBuilder::new(4, &path).unwrap();
+        for slot in 0..4 {
+            assert_eq!(builder.get(slot), None);
+        }
+        builder.close().unwrap();
+    }
+
+    #[test]
+    fn roundtrip_all_valid_masks() {
+        let dir = tempdir().unwrap();
+        let path = dir.path().join("test.psib");
+        let mut builder = SiblingAnnexBuilder::new(4, &path).unwrap();
+
+        let masks = [
+            FamilyMask::EMPTY.with(0),                   // just A: family size 1
+            FamilyMask::EMPTY.with(0).with(3),            // A + T: size 2
+            FamilyMask::EMPTY.with(1).with(2).with(3),    // C+G+T: size 3
+            FamilyMask::EMPTY.with(0).with(1).with(2).with(3), // all 4
+        ];
+        for (slot, mask) in masks.iter().enumerate() {
+            builder.set(slot, *mask);
+        }
+        let annex = builder.finish().unwrap();
+        for (slot, mask) in masks.iter().enumerate() {
+            assert_eq!(annex.get(slot), Some(*mask));
+        }
+        assert_eq!(annex.get(0).unwrap().siblings(), 0);
+        assert_eq!(annex.get(1).unwrap().siblings(), 1);
+        assert_eq!(annex.get(2).unwrap().siblings(), 2);
+        assert_eq!(annex.get(3).unwrap().siblings(), 3);
+        assert_eq!(annex.get(3).unwrap().family_size(), 4);
+    }
+
+    #[test]
+    fn has_reflects_individual_bits() {
+        let mask = FamilyMask::EMPTY.with(0).with(2);
+        assert!(mask.has(0));
+        assert!(!mask.has(1));
+        assert!(mask.has(2));
+        assert!(!mask.has(3));
+    }
+}

src/obicompactvec/src/traits.rs

@@ -1,5 +1,16 @@
 use ndarray::{Array1, Array2};
 
+/// Convert a Jaccard distance matrix (`1 - J`) into a Mash distance matrix, per
+/// https://mash.readthedocs.io/en/latest/distances.html:
+/// `D = -1/k * ln(2J / (1+J))`.
+fn jaccard_to_mash(d_jaccard: &Array2<f64>, k: usize) -> Array2<f64> {
+    d_jaccard.mapv(|d| {
+        let j = 1.0 - d;
+        if j <= 0.0 { 1.0 }
+        else { -1.0 / k as f64 * (2.0 * j / (1.0 + j)).ln() }
+    })
+}
+
 // ── Column-level weight statistic — total count or presence count per column.
 /// Additive across layers and partitions; used as denominator in normalised distances.
 ///
@@ -74,6 +85,12 @@ pub trait CountPartials: ColumnWeights {
         m
     }
 
+    /// Mash distance (https://mash.readthedocs.io/en/latest/distances.html), derived
+    /// from the presence-threshold Jaccard distance.
+    fn threshold_mash_dist_matrix(&self, k: usize, threshold: u32) -> Array2<f64> {
+        jaccard_to_mash(&self.threshold_jaccard_dist_matrix(threshold), k)
+    }
+
     fn relfreq_bray_dist_matrix(&self) -> Array2<f64> {
         let global = self.col_weights();
         let mut m = self.partial_relfreq_bray(&global).mapv(|v| 1.0 - v);
@@ -126,6 +143,12 @@ pub trait BitPartials: ColumnWeights {
         m
     }
 
+    /// Mash distance (https://mash.readthedocs.io/en/latest/distances.html), derived
+    /// from the Jaccard distance.
+    fn mash_dist_matrix(&self, k: usize) -> Array2<f64> {
+        jaccard_to_mash(&self.jaccard_dist_matrix(), k)
+    }
+
     fn hamming_dist_matrix(&self) -> Array2<u64> {
         self.partial_hamming()
     }

src/obikentropy/Cargo.toml

@@ -0,0 +1,10 @@
+[package]
+name = "obikentropy"
+version = "0.1.0"
+edition = "2024"
+
+[dependencies]
+obikseq = { path = "../obikseq" }
+
+[dev-dependencies]
+obikseq = { path = "../obikseq", features = ["test-utils"] }

src/obikentropy/build.rs

@@ -4,57 +4,6 @@ use std::path::PathBuf;
 const K_MAX: usize = 32;
 const WS_MAX: usize = 6;
 
-fn normalize_circular(kmer: u64, ws: usize) -> u64 {
-    let mask = (1u64 << (ws * 2)) - 1;
-    let mut canonical = kmer & mask;
-    let mut current = canonical;
-    for _ in 0..ws - 1 {
-        let top = (current >> ((ws - 1) * 2)) & 3;
-        current = ((current << 2) | top) & mask;
-        if current < canonical {
-            canonical = current;
-        }
-    }
-    canonical
-}
-
-fn revcomp_raw(x: u64, k: usize) -> u64 {
-    let x = !x;
-    let x = x.swap_bytes();
-    let x = ((x >> 4) & 0x0F0F0F0F0F0F0F0F) | ((x & 0x0F0F0F0F0F0F0F0F) << 4);
-    let x = ((x >> 2) & 0x3333333333333333) | ((x & 0x3333333333333333) << 2);
-    x << (64 - 2 * k)
-}
-
-fn build_normalized_kmer(k: usize) -> Vec<u64> {
-    let n = 1usize << (k * 2);
-    let shift = 64 - k * 2;
-    let mut result = vec![0u64; n];
-    for i in 0..n {
-        let la = (i as u64) << shift;
-        let ra = i as u64;
-        let rc_ra = revcomp_raw(la, k) >> shift;
-        let circ = normalize_circular(ra, k);
-        let circ_rc = normalize_circular(rc_ra, k);
-        result[i] = if circ < circ_rc { circ } else { circ_rc };
-    }
-    result
-}
-
-fn build_ln_class(norm: &[u64]) -> Vec<f64> {
-    let n = norm.len();
-    let mut sizes = vec![0u32; n];
-    for &c in norm {
-        sizes[c as usize] += 1;
-    }
-    norm.iter()
-        .map(|&c| {
-            let s = sizes[c as usize];
-            if s > 0 { (s as f64).ln() } else { 0.0 }
-        })
-        .collect()
-}
-
 fn build_n_log_n() -> [f64; K_MAX + 1] {
     let mut t = [0.0f64; K_MAX + 1];
     for n in 1..=K_MAX {
@@ -63,6 +12,9 @@ fn build_n_log_n() -> [f64; K_MAX + 1] {
     t
 }
 
+/// Max achievable entropy over `4^ws` raw sub-words given only `nwords`
+/// observations (most-uniform integer partition), per
+/// `docmd/theory/entropy.md`.
 fn build_emax() -> [[f64; WS_MAX + 1]; K_MAX + 1] {
     let mut t = [[0.0f64; WS_MAX + 1]; K_MAX + 1];
     for k in 2..=K_MAX {
@@ -125,13 +77,6 @@ fn main() {
     let out_dir = PathBuf::from(std::env::var("OUT_DIR").unwrap());
     let mut out = String::new();
 
-    for k in 1..=6usize {
-        let n = 1usize << (k * 2);
-        let norm = build_normalized_kmer(k);
-        let ln_class = build_ln_class(&norm);
-        emit_f64_1d(&mut out, &format!("LN_CLASS{k}"), n, &ln_class);
-    }
-
     let n_log_n = build_n_log_n();
     emit_f64_1d(&mut out, "N_LOG_N", K_MAX + 1, &n_log_n);
 
@@ -141,5 +86,5 @@ fn main() {
     let log_nwords = build_log_nwords();
     emit_f64_2d(&mut out, "LOG_NWORDS", K_MAX + 1, WS_MAX + 1, &log_nwords);
 
-    fs::write(out_dir.join("ln_class_tables.rs"), out).unwrap();
+    fs::write(out_dir.join("entropy_tables.rs"), out).unwrap();
 }

src/obikentropy/src/kmer_entropy.rs

@@ -0,0 +1,41 @@
+//! Normalized entropy of an isolated, already-built k-mer (e.g. one
+//! reconstructed from an index's `unitigs.bin`, with no surrounding
+//! sequence) — drives the window through [`EntropyTracker`] one base at a
+//! time, exactly like the streaming path, so a `theta` threshold means the
+//! same thing whether applied during index construction or after the fact
+//! (e.g. `obikmer filter`).
+
+use obikseq::CanonicalKmer;
+
+use crate::tracker::EntropyTracker;
+
+/// Extension trait: compute the normalized entropy of a single canonical
+/// k-mer, independent of any surrounding sequence.
+pub trait KmerEntropy {
+    /// Normalized entropy across sub-word orders `1..=level_max` (the
+    /// minimum is taken across orders). Lower means less complex; `theta`
+    /// in `index`/`filter` rejects k-mers with a score `< theta`.
+    fn entropy(&self, level_max: usize) -> f64;
+}
+
+impl KmerEntropy for CanonicalKmer {
+    fn entropy(&self, level_max: usize) -> f64 {
+        let raw = self.raw(); // left-aligned, 2 bits/base, MSB-first
+        let k = obikseq::params::k();
+        let mask = (!0u64) >> (64 - k * 2);
+
+        let mut tracker = EntropyTracker::new(k);
+        let mut rolling: u64 = 0;
+        for i in 0..k {
+            let shift = 64 - 2 * (i + 1);
+            let base = (raw >> shift) & 3;
+            rolling = ((rolling << 2) | base) & mask;
+            tracker.push(i + 1, rolling);
+        }
+        tracker.normalized_entropy(level_max)
+    }
+}
+
+#[cfg(test)]
+#[path = "tests/kmer_entropy.rs"]
+mod tests;

src/obikentropy/src/lib.rs

@@ -0,0 +1,17 @@
+//! Normalized k-mer entropy: formulas, tables, and a streaming tracker.
+//!
+//! This crate holds every piece of the entropy computation described in
+//! `docmd/theory/entropy.md`: the compile-time tables ([`table`], private),
+//! the incremental accumulator ([`EntropyTracker`]) that callers compose
+//! into their own streaming state, and the [`KmerEntropy`] convenience trait
+//! for scoring a single, already-built k-mer.
+
+#![deny(missing_docs)]
+
+mod kmer_entropy;
+mod ring;
+mod table;
+mod tracker;
+
+pub use kmer_entropy::KmerEntropy;
+pub use tracker::EntropyTracker;

src/obikentropy/src/ring.rs

@@ -0,0 +1,40 @@
+//! Stack-allocated ring buffer backing the sliding sub-word windows.
+
+/// Fixed-capacity ring buffer backed by a stack array.
+/// N must be a power of two; operations are branchless via `% N`.
+pub(crate) struct Ring<T: Copy + Default, const N: usize> {
+    buf: [T; N],
+    head: usize,
+    len: usize,
+}
+
+impl<T: Copy + Default, const N: usize> Ring<T, N> {
+    #[inline]
+    pub(crate) fn new() -> Self {
+        Self {
+            buf: [T::default(); N],
+            head: 0,
+            len: 0,
+        }
+    }
+
+    #[inline]
+    pub(crate) fn clear(&mut self) {
+        self.len = 0;
+        self.head = 0;
+    }
+
+    #[inline]
+    pub(crate) fn push_back(&mut self, val: T) {
+        self.buf[(self.head + self.len) % N] = val;
+        self.len += 1;
+    }
+
+    #[inline]
+    pub(crate) fn pop_front(&mut self) -> T {
+        let val = self.buf[self.head];
+        self.head = (self.head + 1) % N;
+        self.len -= 1;
+        val
+    }
+}

src/obikentropy/src/table.rs

@@ -0,0 +1,30 @@
+//! Compile-time tables backing the normalized k-mer entropy formula: the
+//! max-entropy correction for small samples. See `docmd/theory/entropy.md`.
+//!
+//! Entropy is computed directly on raw (non-canonicalized) sub-words — no
+//! equivalence-class folding. Empirically (see the discussion that produced
+//! this crate's history), folding sub-words into circular/revcomp classes
+//! before unfolding them back buys nothing for the invariances it was meant
+//! to guarantee (both hold for raw sub-word entropy already, by a direct
+//! bijection argument for revcomp and by the sliding window's own dynamics
+//! for tandem repeats), while it measurably *weakens* detection of the
+//! low-complexity sequences the filter exists to catch.
+
+include!(concat!(env!("OUT_DIR"), "/entropy_tables.rs"));
+
+pub(crate) const WS_MAX: usize = 6;
+
+#[inline(always)]
+pub(crate) const fn n_log_n(n: usize) -> f64 {
+    N_LOG_N[n]
+}
+
+#[inline(always)]
+pub(crate) const fn emax(k: usize, ws: usize) -> f64 {
+    EMAX[k][ws]
+}
+
+#[inline(always)]
+pub(crate) const fn log_nwords(k: usize, ws: usize) -> f64 {
+    LOG_NWORDS[k][ws]
+}

src/obikentropy/src/tests/kmer_entropy.rs

@@ -0,0 +1,52 @@
+use super::*;
+use obikseq::Sequence;
+use obikseq::kmer::Kmer;
+
+const K: usize = 21;
+const LEVEL_MAX: usize = 6;
+
+fn kmer_from_ascii(seq: &[u8]) -> CanonicalKmer {
+    obikseq::set_k(K);
+    Kmer::from_ascii(seq).expect("valid k-mer sequence").canonical()
+}
+
+#[test]
+fn homopolymer_scores_lower_than_diverse_sequence() {
+    let homopolymer = kmer_from_ascii(b"AAAAAAAAAAAAAAAAAAAAA"); // 21 bases
+    let diverse = kmer_from_ascii(b"CATTAGCGTACCTGATCAGGT"); // 21 bases, same as used elsewhere in this workspace's tests
+
+    let e_homopolymer = homopolymer.entropy(LEVEL_MAX);
+    let e_diverse = diverse.entropy(LEVEL_MAX);
+
+    assert!(
+        e_homopolymer < e_diverse,
+        "homopolymer ({e_homopolymer}) should score lower than a diverse sequence ({e_diverse})"
+    );
+    // A pure homopolymer is the most degenerate case representable — its
+    // score should sit near the bottom of the range, not just "somewhat lower".
+    assert!(e_homopolymer < 0.3, "homopolymer entropy unexpectedly high: {e_homopolymer}");
+}
+
+#[test]
+fn entropy_is_deterministic_for_the_same_kmer() {
+    let a = kmer_from_ascii(b"CATTAGCGTACCTGATCAGGT");
+    let b = kmer_from_ascii(b"CATTAGCGTACCTGATCAGGT");
+    assert_eq!(a.entropy(LEVEL_MAX), b.entropy(LEVEL_MAX));
+}
+
+#[test]
+fn entropy_is_within_zero_one_range() {
+    let mut repeat = "AT".repeat(K / 2 + 1);
+    repeat.truncate(K);
+
+    for seq in [
+        "AAAAAAAAAAAAAAAAAAAAA".to_string(),
+        repeat,
+        "CATTAGCGTACCTGATCAGGT".to_string(),
+    ] {
+        assert_eq!(seq.len(), K, "test sequence must be exactly K bases: {seq:?}");
+        let kmer = kmer_from_ascii(seq.as_bytes());
+        let e = kmer.entropy(LEVEL_MAX);
+        assert!((0.0..=1.0).contains(&e), "entropy {e} out of [0,1] for {seq:?}");
+    }
+}

src/obikentropy/src/tracker.rs

@@ -0,0 +1,255 @@
+//! Incremental (streaming) normalized k-mer entropy.
+//!
+//! [`EntropyTracker`] maintains, over a sliding window of the last `k` bases,
+//! the per-sub-word-size raw-word frequency statistics needed to evaluate
+//! the corrected Shannon entropy described in `docmd/theory/entropy.md`,
+//! updated in O(1) per base rather than recomputed from scratch. No
+//! canonicalization is applied — each sub-word is tallied under its own raw
+//! 2-bit-packed value; only the small-sample max-entropy correction departs
+//! from a textbook Shannon entropy.
+//!
+//! It carries no notion of minimizers or superkmer segmentation — callers
+//! that need both (e.g. `obiskbuilder::RollingStat`) compose an
+//! `EntropyTracker` as a plain field alongside their own state, so the two
+//! concerns update in the same streaming pass without being conflated in one
+//! struct.
+
+use crate::ring::Ring;
+use crate::table::{WS_MAX, emax, log_nwords, n_log_n};
+
+/// Incremental normalized-entropy accumulator over a sliding window of `k`
+/// bases. Composed as a plain field by callers that also need other
+/// per-base state (e.g. minimizer selection) in the same streaming pass.
+pub struct EntropyTracker {
+    k: usize,
+    steady: bool,
+
+    // Sliding-window queues over the last `k` raw sub-words, one per word
+    // size — stack-allocated, capacity ≤ k ≤ 31.
+    k1q: Ring<u64, 32>,
+    k2q: Ring<u64, 32>,
+    k3q: Ring<u64, 32>,
+    k4q: Ring<u64, 32>,
+    k5q: Ring<u64, 32>,
+    k6q: Ring<u64, 32>,
+
+    // Frequency count arrays, indexed by the raw sub-word value (2 bits per
+    // base). Max count per cell ≤ k ≤ 31 → u8 is sufficient.
+    k1c: [u8; 4],
+    k2c: [u8; 16],
+    k3c: [u8; 64],
+    k4c: [u8; 256],
+    k5c: [u8; 1024],
+    k6c: [u8; 4096],
+
+    sum_f_log_f: [f64; WS_MAX + 1],
+}
+
+impl EntropyTracker {
+    /// New tracker for a window of `k` bases (1..=31).
+    pub fn new(k: usize) -> Self {
+        Self {
+            k,
+            steady: false,
+            k1q: Ring::new(),
+            k2q: Ring::new(),
+            k3q: Ring::new(),
+            k4q: Ring::new(),
+            k5q: Ring::new(),
+            k6q: Ring::new(),
+            k1c: [0; 4],
+            k2c: [0; 16],
+            k3c: [0; 64],
+            k4c: [0; 256],
+            k5c: [0; 1024],
+            k6c: [0; 4096],
+            sum_f_log_f: [0.0; WS_MAX + 1],
+        }
+    }
+
+    /// Clear all accumulated state, ready to track a new window from
+    /// scratch (`k` is unchanged).
+    pub fn reset(&mut self) {
+        self.steady = false;
+
+        self.k1c.fill(0);
+        self.k2c.fill(0);
+        self.k3c.fill(0);
+        self.k4c.fill(0);
+        self.k5c.fill(0);
+        self.k6c.fill(0);
+
+        self.k1q.clear();
+        self.k2q.clear();
+        self.k3q.clear();
+        self.k4q.clear();
+        self.k5q.clear();
+        self.k6q.clear();
+
+        self.sum_f_log_f = [0.0; WS_MAX + 1];
+    }
+
+    #[inline]
+    fn update_sums_decrement<const K: usize>(sum_f_log_f: &mut [f64; WS_MAX + 1], f: usize) {
+        sum_f_log_f[K] += n_log_n(f - 1) - n_log_n(f);
+    }
+
+    #[inline]
+    fn update_sums_increment<const K: usize>(sum_f_log_f: &mut [f64; WS_MAX + 1], g: usize) {
+        sum_f_log_f[K] += n_log_n(g + 1) - n_log_n(g);
+    }
+
+    /// Advance the window by one base. `received` is the caller's running
+    /// count of bases pushed so far (1-based, i.e. after this base);
+    /// `rolling_kmer` is the current right-aligned, 2-bit-packed k-mer
+    /// window (same convention as `obiskbuilder::RollingStat::rolling_k`).
+    pub fn push(&mut self, received: usize, rolling_kmer: u64) {
+        let raw1 = rolling_kmer & 3;
+        let raw2 = rolling_kmer & 15;
+        let raw3 = rolling_kmer & 63;
+        let raw4 = rolling_kmer & 255;
+        let raw5 = rolling_kmer & 1023;
+        let raw6 = rolling_kmer & 4095;
+
+        if received > self.k {
+            let old1 = self.k1q.pop_front();
+            let f1 = self.k1c[old1 as usize] as usize;
+            Self::update_sums_decrement::<1>(&mut self.sum_f_log_f, f1);
+            self.k1c[old1 as usize] -= 1;
+
+            let old2 = self.k2q.pop_front();
+            let f2 = self.k2c[old2 as usize] as usize;
+            Self::update_sums_decrement::<2>(&mut self.sum_f_log_f, f2);
+            self.k2c[old2 as usize] -= 1;
+
+            let old3 = self.k3q.pop_front();
+            let f3 = self.k3c[old3 as usize] as usize;
+            Self::update_sums_decrement::<3>(&mut self.sum_f_log_f, f3);
+            self.k3c[old3 as usize] -= 1;
+
+            let old4 = self.k4q.pop_front();
+            let f4 = self.k4c[old4 as usize] as usize;
+            Self::update_sums_decrement::<4>(&mut self.sum_f_log_f, f4);
+            self.k4c[old4 as usize] -= 1;
+
+            let old5 = self.k5q.pop_front();
+            let f5 = self.k5c[old5 as usize] as usize;
+            Self::update_sums_decrement::<5>(&mut self.sum_f_log_f, f5);
+            self.k5c[old5 as usize] -= 1;
+
+            let old6 = self.k6q.pop_front();
+            let f6 = self.k6c[old6 as usize] as usize;
+            Self::update_sums_decrement::<6>(&mut self.sum_f_log_f, f6);
+            self.k6c[old6 as usize] -= 1;
+        }
+
+        if self.steady {
+            let g1 = self.k1c[raw1 as usize] as usize;
+            Self::update_sums_increment::<1>(&mut self.sum_f_log_f, g1);
+            self.k1c[raw1 as usize] += 1;
+            self.k1q.push_back(raw1);
+
+            let g2 = self.k2c[raw2 as usize] as usize;
+            Self::update_sums_increment::<2>(&mut self.sum_f_log_f, g2);
+            self.k2c[raw2 as usize] += 1;
+            self.k2q.push_back(raw2);
+
+            let g3 = self.k3c[raw3 as usize] as usize;
+            Self::update_sums_increment::<3>(&mut self.sum_f_log_f, g3);
+            self.k3c[raw3 as usize] += 1;
+            self.k3q.push_back(raw3);
+
+            let g4 = self.k4c[raw4 as usize] as usize;
+            Self::update_sums_increment::<4>(&mut self.sum_f_log_f, g4);
+            self.k4c[raw4 as usize] += 1;
+            self.k4q.push_back(raw4);
+
+            let g5 = self.k5c[raw5 as usize] as usize;
+            Self::update_sums_increment::<5>(&mut self.sum_f_log_f, g5);
+            self.k5c[raw5 as usize] += 1;
+            self.k5q.push_back(raw5);
+
+            let g6 = self.k6c[raw6 as usize] as usize;
+            Self::update_sums_increment::<6>(&mut self.sum_f_log_f, g6);
+            self.k6c[raw6 as usize] += 1;
+            self.k6q.push_back(raw6);
+        } else {
+            self.push_warmup_increments(received, raw1, raw2, raw3, raw4, raw5, raw6);
+        }
+    }
+
+    #[cold]
+    #[inline(never)]
+    fn push_warmup_increments(
+        &mut self,
+        received: usize,
+        raw1: u64, raw2: u64, raw3: u64,
+        raw4: u64, raw5: u64, raw6: u64,
+    ) {
+        let g1 = self.k1c[raw1 as usize] as usize;
+        Self::update_sums_increment::<1>(&mut self.sum_f_log_f, g1);
+        self.k1c[raw1 as usize] += 1;
+        self.k1q.push_back(raw1);
+
+        if received >= 2 {
+            let g2 = self.k2c[raw2 as usize] as usize;
+            Self::update_sums_increment::<2>(&mut self.sum_f_log_f, g2);
+            self.k2c[raw2 as usize] += 1;
+            self.k2q.push_back(raw2);
+
+            if received >= 3 {
+                let g3 = self.k3c[raw3 as usize] as usize;
+                Self::update_sums_increment::<3>(&mut self.sum_f_log_f, g3);
+                self.k3c[raw3 as usize] += 1;
+                self.k3q.push_back(raw3);
+
+                if received >= 4 {
+                    let g4 = self.k4c[raw4 as usize] as usize;
+                    Self::update_sums_increment::<4>(&mut self.sum_f_log_f, g4);
+                    self.k4c[raw4 as usize] += 1;
+                    self.k4q.push_back(raw4);
+
+                    if received >= 5 {
+                        let g5 = self.k5c[raw5 as usize] as usize;
+                        Self::update_sums_increment::<5>(&mut self.sum_f_log_f, g5);
+                        self.k5c[raw5 as usize] += 1;
+                        self.k5q.push_back(raw5);
+
+                        if received >= 6 {
+                            let g6 = self.k6c[raw6 as usize] as usize;
+                            Self::update_sums_increment::<6>(&mut self.sum_f_log_f, g6);
+                            self.k6c[raw6 as usize] += 1;
+                            self.k6q.push_back(raw6);
+                            self.steady = true;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    /// Normalized entropy at sub-word size `order` (1..=6). The caller is
+    /// responsible for not calling this before the window is full (`k`
+    /// bases pushed) — an empty/partial window yields a meaningless value.
+    pub fn entropy(&self, order: usize) -> f64 {
+        let k = self.k;
+        let em = emax(k, order);
+        if em <= 0.0 {
+            return 1.0;
+        }
+        let nwords = k - order + 1;
+        let log_nw = log_nwords(k, order);
+        let nw_f = nwords as f64;
+        let h_corr = log_nw - self.sum_f_log_f[order] / nw_f;
+        (h_corr / em).max(0.0)
+    }
+
+    /// Minimum of [`Self::entropy`] over sub-word sizes `1..=order_max`, same
+    /// caller responsibility re: window readiness as `entropy`.
+    pub fn normalized_entropy(&self, order_max: usize) -> f64 {
+        let min_e = (1..=order_max)
+            .map(|ws| self.entropy(ws))
+            .fold(f64::MAX, f64::min);
+        if min_e == f64::MAX { 1.0 } else { min_e }
+    }
+}

src/obikindex/Cargo.toml

@@ -10,6 +10,8 @@ obiskio          = { path = "../obiskio" }
 obisys           = { path = "../obisys" }
 obicompactvec    = { path = "../obicompactvec" }
 obilayeredmap    = { path = "../obilayeredmap" }
+obiskbuilder     = { path = "../obiskbuilder" }
+obipipeline      = { path = "../obipipeline" }
 ndarray          = "0.16"
 rayon             = "1"
 crossbeam-channel = "0.5"
@@ -19,6 +21,10 @@ indicatif        = "0.17"
 tracing          = "0.1.44"
 hwlocality       = { version = "1.0.0-alpha.11", features = ["vendored"], optional = true }
 
+[dev-dependencies]
+obiread  = { path = "../obiread" }
+tempfile = "3"
+
 [features]
 default = ["numa"]
 numa    = ["hwlocality"]

src/obikindex/src/distance.rs

@@ -14,6 +14,8 @@ pub enum DistanceMetric {
     Jaccard,
     /// Hamming distance (number of differing kmer positions) on presence/absence data.
     Hamming,
+    /// Mash distance on presence/absence data (Jaccard-derived mutation-rate estimate).
+    Mash,
     /// Bray-Curtis dissimilarity on raw counts.
     BrayCurtis,
     /// Bray-Curtis dissimilarity normalised by per-genome total counts.
@@ -84,6 +86,7 @@ impl KmerIndex {
                 DistanceMetric::Hellinger         => CountPartials::hellinger_dist_matrix(&global),
                 DistanceMetric::HellingerEuclidean => CountPartials::hellinger_euclidean_dist_matrix(&global),
                 DistanceMetric::Jaccard           => CountPartials::threshold_jaccard_dist_matrix(&global, presence_threshold),
+                DistanceMetric::Mash              => CountPartials::threshold_mash_dist_matrix(&global, self.kmer_size(), presence_threshold),
                 DistanceMetric::Hamming           => {
                     return Err(OKIError::InvalidInput(
                         "Hamming is only available for presence/absence indexes".into(),
@@ -108,6 +111,7 @@ impl KmerIndex {
 
             let matrix = match metric {
                 DistanceMetric::Jaccard => BitPartials::jaccard_dist_matrix(&global),
+                DistanceMetric::Mash    => BitPartials::mash_dist_matrix(&global, self.kmer_size()),
                 DistanceMetric::Hamming => {
                     BitPartials::hamming_dist_matrix(&global).mapv(|v| v as f64)
                 }

src/obikindex/src/lib.rs

@@ -9,6 +9,7 @@ mod numa;
 mod rebuild;
 mod reindex;
 mod select;
+mod siblings;
 mod stats;
 
 pub use error::{OKIError, OKIResult};
@@ -18,3 +19,4 @@ pub use merge::MergeMode;
 pub use meta::{validate_label, GenomeInfo, IndexConfig, IndexMeta, META_FILENAME};
 pub use state::{IndexState, SENTINEL_COUNTED, SENTINEL_INDEXED, SENTINEL_SCATTERED};
 pub use stats::IndexBitsPerKmer;
+pub use siblings::{RawSnpDistanceOutput, SiblingAnnexStats, SnpAlignment};

src/obikindex/src/numa.rs

@@ -79,9 +79,7 @@ pub fn build() -> NumaSetup {
     }
 
     // UMA fallback: single synthetic node, all cores, no pool, no pinning.
-    let n_cores = std::thread::available_parallelism()
-        .map(|n| n.get())
-        .unwrap_or(1);
+    let n_cores = obisys::effective_parallelism();
     debug!("UMA: single synthetic node, {} core(s)", n_cores);
     NumaSetup {
         pools: vec![None],
@@ -91,9 +89,7 @@ pub fn build() -> NumaSetup {
 
 #[cfg(not(feature = "numa"))]
 pub fn build() -> NumaSetup {
-    let n_cores = std::thread::available_parallelism()
-        .map(|n| n.get())
-        .unwrap_or(1);
+    let n_cores = obisys::effective_parallelism();
     debug!("UMA: single synthetic node, {} core(s)", n_cores);
     NumaSetup {
         pools: vec![None],

src/obikmer/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "obikmer"
-version = "1.1.37"
+version = "1.1.43"
 edition = "2024"
 
 [[bin]]

src/obikmer/src/cli.rs

@@ -38,9 +38,7 @@ pub struct CommonArgs {
     #[arg(
         short = 'T',
         long,
-        default_value_t = std::thread::available_parallelism()
-            .map(|n| n.get())
-            .unwrap_or(1)
+        default_value_t = obisys::effective_parallelism()
     )]
     pub threads: usize,
 

src/obikmer/src/cmd/distance.rs

@@ -3,13 +3,15 @@ use std::path::PathBuf;
 
 use clap::Args;
 use kodama::{Method, linkage};
-use obikindex::{DistanceMetric, KmerIndex};
+use obifastwrite::{JsonVal, write_record};
+use obikindex::{DistanceMetric, KmerIndex, RawSnpDistanceOutput, SiblingAnnexStats, SnpAlignment};
 use speedytree::{DistanceMatrix, Hybrid, NeighborJoiningSolver, to_newick};
 use tracing::info;
 
 #[derive(clap::ValueEnum, Clone, Copy, Debug)]
 pub enum MetricArg {
     Jaccard,
+    Mash,
     Hamming,
     BrayCurtis,
     #[value(name = "relfreq-bray-curtis")]
@@ -26,6 +28,7 @@ impl From<MetricArg> for DistanceMetric {
     fn from(m: MetricArg) -> Self {
         match m {
             MetricArg::Jaccard            => DistanceMetric::Jaccard,
+            MetricArg::Mash               => DistanceMetric::Mash,
             MetricArg::Hamming            => DistanceMetric::Hamming,
             MetricArg::BrayCurtis         => DistanceMetric::BrayCurtis,
             MetricArg::RelfreqBrayCurtis  => DistanceMetric::RelfreqBrayCurtis,
@@ -62,7 +65,37 @@ pub struct DistanceArgs {
     #[arg(long)]
     pub upgma: bool,
 
+    /// Build the sibling-count/minorant annex on this (multi-genome) index
+    /// — see `docmd/theory/evolutionary_distances.md`, Step 2b. Construction
+    /// only; does not by itself compute or write any statistics.
+    #[arg(long)]
+    pub sibling_annex: bool,
+
+    /// Tally the sibling-count distribution (CSV) of an already-built annex
+    /// (run with `--sibling-annex` first, in this invocation or an earlier
+    /// one). A separate, occasional diagnostic pass — not run every time the
+    /// annex itself is (re)built.
+    #[arg(long)]
+    pub sibling_stats: bool,
+
+    /// Compute the raw p-distance restricted to loci that are single-copy
+    /// in both genomes of each pair (an already-built sibling annex is
+    /// required — run with `--sibling-annex` first, in this invocation or
+    /// an earlier one). A quick way to test the central-position SNP
+    /// estimator against a real index; not the full `SnpTally` design.
+    #[arg(long)]
+    pub raw_snp_distance: bool,
+
+    /// Write a SNP-only pseudo-alignment (FASTA, IUPAC-coded) from an
+    /// already-built sibling annex — one row per genome, one column per
+    /// variable family (monomorphic families skipped), no flanking
+    /// sequence. See `docmd/theory/evolutionary_distances.md`,
+    /// "Multi-genome framing: family as pseudo-alignment column".
+    #[arg(long)]
+    pub snp: bool,
+
     /// Output prefix: <prefix>_dist.csv, <prefix>_shared.csv,
+    /// <prefix>_siblings.csv, <prefix>_rawsnp.csv, <prefix>_snp.fasta,
     /// <prefix>_nj.nwk, <prefix>_upgma.nwk.
     /// If omitted, the distance matrix is written to stdout.
     #[arg(short, long)]
@@ -78,6 +111,51 @@ pub fn run(args: DistanceArgs) {
 
     let labels: Vec<String> = idx.meta().genomes.iter().map(|g| g.label.clone()).collect();
     let n = labels.len();
+
+    // ── Sibling-count/minorant annex (independent of the distance metric) ──
+    // Construction (`--sibling-annex`) and stats (`--sibling-stats`) are
+    // deliberately decoupled: the annex is meant to be (re)built routinely,
+    // the distribution only occasionally, on demand.
+    if args.sibling_annex {
+        info!("building sibling-count/minorant annex");
+        idx.build_sibling_annex().unwrap_or_else(|e| {
+            eprintln!("error building sibling annex: {e}");
+            std::process::exit(1);
+        });
+    }
+    if args.sibling_stats {
+        let stats = idx.sibling_annex_stats().unwrap_or_else(|e| {
+            eprintln!("error computing sibling-annex stats: {e}");
+            std::process::exit(1);
+        });
+        write_sibling_stats_csv(&stats, &labels, &args.output);
+    }
+    if args.raw_snp_distance {
+        let result = idx.raw_snp_distance().unwrap_or_else(|e| {
+            eprintln!("error computing raw SNP distance: {e}");
+            std::process::exit(1);
+        });
+        write_raw_snp_distance_csv(&result, &labels, &args.output);
+    }
+    if args.snp {
+        let alignment = idx.snp_pseudo_alignment().unwrap_or_else(|e| {
+            eprintln!("error computing SNP pseudo-alignment: {e}");
+            std::process::exit(1);
+        });
+        write_snp_fasta(&alignment, &labels, &args.output);
+    }
+
+    // `--sibling-annex`/`--sibling-stats`/`--raw-snp-distance`/`--snp` are
+    // their own operation, not a modifier on top of a distance-metric
+    // computation — a metric was never requested by asking for any of them,
+    // so there is nothing for the rest of this function to compute. Not a
+    // historical accident to keep: stop here rather than always also
+    // running a Jaccard (or whichever `--metric` defaults to) pass and
+    // printing an unrequested matrix.
+    if args.sibling_annex || args.sibling_stats || args.raw_snp_distance || args.snp {
+        return;
+    }
+
     info!(
         "computing {:?} distances for {} genome(s)",
         args.metric, n
@@ -189,6 +267,103 @@ pub fn run(args: DistanceArgs) {
     }
 }
 
+// ── Family-size distribution → CSV ──────────────────────────────────────────
+//
+// Each row is a family (the up-to-4 k-mers sharing flanks, differing only at
+// the centre), counted once — at its minorant — regardless of how many of
+// its members are observed. Family size 1..4 (not "sibling count" 0..3):
+// see `docmd/theory/evolutionary_distances.md`, "Definitions".
+
+fn write_sibling_stats_csv(stats: &SiblingAnnexStats, labels: &[String], output: &Option<PathBuf>) {
+    // One row per genome (4 columns, family size 1-4: number of families of
+    // that size for which the genome carries at least one member), plus a
+    // `global` row — the actual deduplicated family-size histogram
+    // (`stats.counts`), NOT a sum of the per-genome columns (a family shared
+    // by several genomes would otherwise be counted once per genome it
+    // appears in, inflating the total beyond the real family count).
+    let path = output.as_ref()
+        .map(|p| format!("{}_siblings.csv", p.display()))
+        .unwrap_or_else(|| "siblings.csv".into());
+    let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| {
+        eprintln!("error creating {path}: {e}");
+        std::process::exit(1);
+    }));
+    writeln!(f, "genome,1,2,3,4").unwrap();
+    for (label, counts) in labels.iter().zip(stats.per_genome.iter()) {
+        writeln!(f, "{label},{},{},{},{}", counts[0], counts[1], counts[2], counts[3]).unwrap();
+    }
+    writeln!(
+        f, "global,{},{},{},{}",
+        stats.counts[0], stats.counts[1], stats.counts[2], stats.counts[3],
+    ).unwrap();
+    let total: u64 = stats.counts.iter().sum();
+    info!("family-size distribution → {path} (total {total} famil{})",
+        if total == 1 { "y" } else { "ies" });
+}
+
+// ── Raw single-copy SNP distance → CSV ──────────────────────────────────────
+//
+// p_hat[i,j] = snp[i,j] / (snp[i,j] + shared[i,j]) over loci single-copy in
+// both i and j — see `RawSnpDistanceOutput` / `KmerIndex::raw_snp_distance`.
+// A single file: the distance matrix, with an eligible-loci count alongside
+// each value so a 0/0 pair (no eligible locus at all) is distinguishable
+// from a genuinely identical pair.
+
+fn write_raw_snp_distance_csv(result: &RawSnpDistanceOutput, labels: &[String], output: &Option<PathBuf>) {
+    let path = output.as_ref()
+        .map(|p| format!("{}_rawsnp.csv", p.display()))
+        .unwrap_or_else(|| "rawsnp.csv".into());
+    let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| {
+        eprintln!("error creating {path}: {e}");
+        std::process::exit(1);
+    }));
+    let n = labels.len();
+    write!(f, "genome").unwrap();
+    for g in labels { write!(f, ",{g}").unwrap(); }
+    writeln!(f).unwrap();
+    for (i, g) in labels.iter().enumerate() {
+        write!(f, "{g}").unwrap();
+        for j in 0..n {
+            let snp = result.snp[[i, j]];
+            let shared = result.shared[[i, j]];
+            let eligible = snp + shared;
+            if eligible == 0 {
+                write!(f, ",NA").unwrap();
+            } else {
+                write!(f, ",{:.6}", snp as f64 / eligible as f64).unwrap();
+            }
+        }
+        writeln!(f).unwrap();
+    }
+    info!("raw single-copy SNP distance matrix → {path}");
+}
+
+// ── SNP-only pseudo-alignment → FASTA ───────────────────────────────────────
+//
+// One record per genome, IUPAC-coded, no flanking sequence — see
+// `SnpAlignment` / `KmerIndex::snp_pseudo_alignment`. Uses the project's
+// existing FASTA writer (`obifastwrite::write_record`) rather than
+// hand-rolling one.
+
+fn write_snp_fasta(alignment: &SnpAlignment, labels: &[String], output: &Option<PathBuf>) {
+    let path = output.as_ref()
+        .map(|p| format!("{}_snp.fasta", p.display()))
+        .unwrap_or_else(|| "snp.fasta".into());
+    let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| {
+        eprintln!("error creating {path}: {e}");
+        std::process::exit(1);
+    }));
+    let n_sites = alignment.sequences.first().map(|s| s.len()).unwrap_or(0);
+    for (label, seq) in labels.iter().zip(alignment.sequences.iter()) {
+        write_record(seq, label, &[("n_sites", JsonVal::Num(n_sites as u64))], &mut f).unwrap_or_else(|e| {
+            eprintln!("error writing {path}: {e}");
+            std::process::exit(1);
+        });
+    }
+    info!("SNP pseudo-alignment → {path} ({n_sites} site{})",
+        if n_sites == 1 { "" } else { "s" });
+}
+
 // ── UPGMA Newick from kodama dendrogram ───────────────────────────────────────
 
 fn upgma_to_newick(dendro: &kodama::Dendrogram<f64>, names: &[String]) -> String {

src/obikmer/src/cmd/filter.rs

@@ -2,14 +2,14 @@ use std::path::PathBuf;
 
 use clap::Args;
 use obikindex::{KmerIndex, MergeMode};
-use obikpartitionner::filter::{MaxTotalCount, MinTotalCount};
+use obikpartitionner::filter::{MaxTotalCount, MinComplexity, MinTotalCount};
 use obisys::Reporter;
 use tracing::info;
 
 use super::predicate::FilterArgs as KmerFilterArgs;
 
 #[derive(Args)]
-pub struct FilterArgs {
+pub struct FilterCmdArgs {
     /// Source index directory
     pub source: PathBuf,
 
@@ -28,6 +28,18 @@ pub struct FilterArgs {
     #[arg(long)]
     pub max_total_count: Option<u32>,
 
+    /// Minimum normalized entropy (complexity) to keep a k-mer — same metric
+    /// as `obikmer index`'s --theta, applied here to k-mers already committed
+    /// to the source index (reconstructed from unitigs.bin). K-mers scoring
+    /// below this are removed.
+    #[arg(long)]
+    pub min_complexity: Option<f64>,
+
+    /// Maximum sub-word size for the complexity computation (see `obikmer
+    /// index`'s --level-max). Only used when --min-complexity is set.
+    #[arg(long, default_value_t = 6)]
+    pub complexity_level_max: usize,
+
     /// Output as presence/absence instead of counts
     #[arg(long)]
     pub presence: bool,
@@ -37,7 +49,7 @@ pub struct FilterArgs {
     pub force: bool,
 }
 
-pub fn run(args: FilterArgs) {
+pub fn run(args: FilterCmdArgs) {
     let src = KmerIndex::open(&args.source).unwrap_or_else(|e| {
         eprintln!("error opening source index: {e}");
         std::process::exit(1);
@@ -62,6 +74,9 @@ pub fn run(args: FilterArgs) {
     if let Some(v) = args.max_total_count {
         filters.push(Box::new(MaxTotalCount { total: v }));
     }
+    if let Some(theta) = args.min_complexity {
+        filters.push(Box::new(MinComplexity { level_max: args.complexity_level_max, theta }));
+    }
 
     let mut rep = Reporter::new();
     KmerIndex::rebuild(&args.output, &src, &filters, mode, args.force, &mut rep)

src/obikmer/src/cmd/predicate.rs

@@ -151,12 +151,14 @@ pub struct FilterArgs {
     pub outgroup: Vec<String>,
 
     /// Minimum number of ingroup genomes containing the k-mer
-    #[arg(long)]
-    pub min_count: Option<usize>,
+    /// (negative: offset from group size, e.g. -1 = all but one)
+    #[arg(long, allow_hyphen_values = true)]
+    pub min_count: Option<isize>,
 
     /// Maximum number of ingroup genomes containing the k-mer
-    #[arg(long)]
-    pub max_count: Option<usize>,
+    /// (negative: offset from group size, e.g. -1 = all but one)
+    #[arg(long, allow_hyphen_values = true)]
+    pub max_count: Option<isize>,
 
     /// Minimum fraction of ingroup genomes containing the k-mer [0.0–1.0]
     /// (default 1.0 when --ingroup is set, 0.0 otherwise)
@@ -168,13 +170,15 @@ pub struct FilterArgs {
     pub max_frac: Option<f64>,
 
     /// Minimum number of outgroup genomes containing the k-mer
-    #[arg(long)]
-    pub min_outgroup_count: Option<usize>,
+    /// (negative: offset from outgroup size, e.g. -1 = all but one)
+    #[arg(long, allow_hyphen_values = true)]
+    pub min_outgroup_count: Option<isize>,
 
     /// Maximum number of outgroup genomes containing the k-mer
-    /// (default 0 when --outgroup is set, no constraint otherwise)
-    #[arg(long)]
-    pub max_outgroup_count: Option<usize>,
+    /// (default 0 when --outgroup is set, no constraint otherwise;
+    /// negative: offset from outgroup size, e.g. -1 = all but one)
+    #[arg(long, allow_hyphen_values = true)]
+    pub max_outgroup_count: Option<isize>,
 
     /// Minimum fraction of outgroup genomes containing the k-mer [0.0–1.0]
     #[arg(long)]
@@ -239,12 +243,12 @@ pub fn matching_genome_indices(pred_str: &str, genomes: &[GenomeInfo]) -> Result
 
 pub struct GroupFilterParams {
     pub threshold:          u32,
-    pub min_count:          Option<usize>,
-    pub max_count:          Option<usize>,
+    pub min_count:          Option<isize>,
+    pub max_count:          Option<isize>,
     pub min_frac:           Option<f64>,
     pub max_frac:           Option<f64>,
-    pub min_outgroup_count: Option<usize>,
-    pub max_outgroup_count: Option<usize>,
+    pub min_outgroup_count: Option<isize>,
+    pub max_outgroup_count: Option<isize>,
     pub min_outgroup_frac:  Option<f64>,
     pub max_outgroup_frac:  Option<f64>,
 }
@@ -279,12 +283,20 @@ pub fn build_group_filter(
     let default_min_frac           = if !ingroup_preds.is_empty()  && !ingroup_quorum_explicit  { 1.0      } else { 0.0      };
     let default_max_outgroup_count = if !outgroup_preds.is_empty() && !outgroup_quorum_explicit { 0        } else { out_size };
 
-    let min_count          = p.min_count.unwrap_or(0);
-    let max_count          = p.max_count.unwrap_or(in_size);
+    // Resolve a signed count: negative means an offset from the group size
+    // (e.g. -1 = all but one), floored at 1 so the negative form always keeps
+    // constraining the group — even a singleton group, where n-1 would be 0
+    // and would otherwise drop the constraint entirely.
+    let resolve = |v: isize, size: usize| -> usize {
+        if v < 0 { (size as isize + v).max(1) as usize } else { v as usize }
+    };
+
+    let min_count          = p.min_count.map(|v| resolve(v, in_size)).unwrap_or(0);
+    let max_count          = p.max_count.map(|v| resolve(v, in_size)).unwrap_or(in_size);
     let min_frac           = p.min_frac.unwrap_or(default_min_frac);
     let max_frac           = p.max_frac.unwrap_or(1.0);
-    let min_outgroup_count = p.min_outgroup_count.unwrap_or(0);
-    let max_outgroup_count = p.max_outgroup_count.unwrap_or(default_max_outgroup_count);
+    let min_outgroup_count = p.min_outgroup_count.map(|v| resolve(v, out_size)).unwrap_or(0);
+    let max_outgroup_count = p.max_outgroup_count.map(|v| resolve(v, out_size)).unwrap_or(default_max_outgroup_count);
     let min_outgroup_frac  = p.min_outgroup_frac.unwrap_or(0.0);
     let max_outgroup_frac  = p.max_outgroup_frac.unwrap_or(1.0);
 

src/obikmer/src/cmd/query.rs

@@ -70,9 +70,7 @@ pub struct QueryArgs {
     #[arg(
         short = 'T',
         long,
-        default_value_t = std::thread::available_parallelism()
-            .map(|n| n.get())
-            .unwrap_or(1)
+        default_value_t = obisys::effective_parallelism()
     )]
     pub threads: usize,
 
@@ -325,6 +323,42 @@ fn process_chunk(
     let batch = QueryBatch::from_records(records, k, 6, 0.7, n_partitions);
     let n_seqs = batch.ids.len();
 
+    // Estimate QueryBatch::by_partition's actual memory footprint: the
+    // k-mer-level dedup map (roadmap point 5) — one HashMap<CanonicalKmer,
+    // Vec<KmerDesc>> per partition, sized by *unique* k-mers, not shrunk by
+    // dedup. On real workloads with a low intra-chunk duplication rate this
+    // can dwarf every other per-chunk structure, including the sparse
+    // Findere ones logged further down — unlike those, chunk_bytes's formula
+    // (run()) does not account for this at all today. Measured by allocated
+    // capacity, not logical length, to reflect real memory pressure
+    // (HashMap/Vec growth slack) — `by_partition` is alive for the entire
+    // process_chunk call (never drained, only iterated by reference), so
+    // this is its footprint for the whole chunk lifetime, not a transient.
+    let hashmap_slot_bytes = (std::mem::size_of::<CanonicalKmer>()
+        + std::mem::size_of::<Vec<KmerDesc>>()
+        + 1) as u64; // +1 ≈ hashbrown control byte per slot
+    let by_partition_map_bytes: u64 = batch
+        .by_partition
+        .iter()
+        .map(|m| m.capacity() as u64 * hashmap_slot_bytes)
+        .sum();
+    let by_partition_desc_bytes: u64 = batch
+        .by_partition
+        .iter()
+        .flat_map(|m| m.values())
+        .map(|v| v.capacity() as u64 * std::mem::size_of::<KmerDesc>() as u64)
+        .sum();
+    let by_partition_bytes = by_partition_map_bytes + by_partition_desc_bytes;
+
+    debug!(
+        n_unique_kmers_total = batch.by_partition.iter().map(|m| m.len() as u64).sum::<u64>(),
+        by_partition_map_bytes,
+        by_partition_desc_bytes,
+        by_partition_bytes,
+        chunk_bytes,
+        "by_partition memory retained"
+    );
+
     // Sparse bookkeeping for the whole chunk:
     // - smer_index: O(total_smers) — is this s-mer in the index at all.
     // - by_genome[g]: raw (seq_idx, pos_smer, value) hits for genome g, only

src/obikmer/src/main.rs

@@ -21,7 +21,7 @@ enum Commands {
     /// Merge multiple built indexes into one
     Merge(cmd::merge::MergeArgs),
     /// Apply row-level selection (σ) to an index: retain only k-mers matching the predicates
-    Filter(cmd::filter::FilterArgs),
+    Filter(cmd::filter::FilterCmdArgs),
     /// Project and/or aggregate genome columns into a new or in-place index
     Select(cmd::select::SelectArgs),
     /// Query an index with sequences and annotate matches

src/obikpartitionner/Cargo.toml

@@ -6,7 +6,6 @@ edition = "2024"
 [dev-dependencies]
 tempfile = "3"
 obikseq   = { path = "../obikseq",   features = ["test-utils"] }
-obiskbuilder = { path = "../obiskbuilder" }
 obiread   = { path = "../obiread" }
 obikrope  = { path = "../obikrope" }
 
@@ -14,6 +13,8 @@ obikrope  = { path = "../obikrope" }
 niffler = "3.0.0"
 remove_dir_all = "0.8"
 obikseq      = { path = "../obikseq" }
+obikentropy  = { path = "../obikentropy" }
+obiskbuilder = { path = "../obiskbuilder" }
 obiskio      = { path = "../obiskio" }
 obidebruinj  = { path = "../obidebruinj" }
 obilayeredmap = { path = "../obilayeredmap" }

src/obikpartitionner/src/dump_layer.rs

@@ -62,7 +62,7 @@ impl KmerPartition {
                 for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
                     if let Some(slot) = mphf.find(kmer) {
                         let row = mat.row(slot);
-                        if passes_all(filters, &row, n_genomes) {
+                        if passes_all(filters, kmer, &row, n_genomes) {
                             cont = cb(kmer, row);
                             if !cont { break; }
                         }
@@ -75,7 +75,7 @@ impl KmerPartition {
                 for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
                     if let Some(slot) = mphf.find(kmer) {
                         let row: Box<[u32]> = mat.row(slot).iter().map(|&b| b as u32).collect();
-                        if passes_all(filters, &row, n_genomes) {
+                        if passes_all(filters, kmer, &row, n_genomes) {
                             cont = cb(kmer, row);
                             if !cont { break; }
                         }
@@ -83,16 +83,17 @@ impl KmerPartition {
                 }
                 cont
             } else {
-                // No data matrix: implicit presence — all values are 1.
-                // The filter result is identical for every kmer, so evaluate once.
+                // No data matrix: implicit presence — all values are 1. `row`
+                // is identical for every kmer, but a filter can still depend
+                // on the kmer's own sequence (e.g. MinComplexity), so this
+                // cannot be evaluated once for the whole layer — filters must
+                // still be tested per kmer.
                 let all_present: Box<[u32]> = vec![1u32; n_genomes].into();
                 let mut cont = true;
-                if passes_all(filters, &all_present, n_genomes) {
-                    for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
-                        if mphf.find(kmer).is_some() {
-                            cont = cb(kmer, all_present.clone());
-                            if !cont { break; }
-                        }
+                for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
+                    if mphf.find(kmer).is_some() && passes_all(filters, kmer, &all_present, n_genomes) {
+                        cont = cb(kmer, all_present.clone());
+                        if !cont { break; }
                     }
                 }
                 cont
@@ -140,7 +141,7 @@ impl KmerPartition {
                 for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
                     if let Some(slot) = mphf.find(kmer) {
                         let row = mat.row(slot);
-                        if passes_all(filters, &row, n_genomes) {
+                        if passes_all(filters, kmer, &row, n_genomes) {
                             cont = cb(part, layer, kmer, row);
                             if !cont { break; }
                         }
@@ -153,7 +154,7 @@ impl KmerPartition {
                 for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
                     if let Some(slot) = mphf.find(kmer) {
                         let row: Box<[u32]> = mat.row(slot).iter().map(|&b| b as u32).collect();
-                        if passes_all(filters, &row, n_genomes) {
+                        if passes_all(filters, kmer, &row, n_genomes) {
                             cont = cb(part, layer, kmer, row);
                             if !cont { break; }
                         }
@@ -161,14 +162,15 @@ impl KmerPartition {
                 }
                 cont
             } else {
+                // Same as iter_partition_kmers: row is constant but a filter
+                // may still depend on the kmer's own sequence, so this must
+                // be tested per kmer, not once for the whole layer.
                 let all_present: Box<[u32]> = vec![1u32; n_genomes].into();
                 let mut cont = true;
-                if passes_all(filters, &all_present, n_genomes) {
-                    for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
-                        if mphf.find(kmer).is_some() {
-                            cont = cb(part, layer, kmer, all_present.clone());
-                            if !cont { break; }
-                        }
+                for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
+                    if mphf.find(kmer).is_some() && passes_all(filters, kmer, &all_present, n_genomes) {
+                        cont = cb(part, layer, kmer, all_present.clone());
+                        if !cont { break; }
                     }
                 }
                 cont

src/obikpartitionner/src/filter.rs

@@ -1,17 +1,24 @@
 use obicompactvec::FilterMask;
+use obikseq::CanonicalKmer;
 
-/// Trait for kmer row filters.
+/// Trait for kmer filters.
 ///
+/// `kmer` is the k-mer's own canonical sequence, reconstructed from the
+/// source index's `unitigs.bin` (always present — see `rebuild_layer.rs`);
 /// `row` contains raw per-genome counts (or 0/1 for presence/absence data).
-/// `n_genomes` equals `row.len()`.
+/// `n_genomes` equals `row.len()`. Most filters only need `row` — `kmer` is
+/// there for filters that reason about the k-mer's sequence itself (e.g.
+/// [`MinComplexity`]).
 pub trait KmerFilter: Send + Sync {
-    fn passes(&self, row: &[u32], n_genomes: usize) -> bool;
+    fn passes(&self, kmer: CanonicalKmer, row: &[u32], n_genomes: usize) -> bool;
 
     /// Express this filter as a [`FilterMask`] column-operation expression.
     ///
     /// Returns `Some(expr)` if the filter can be evaluated solely from matrix
     /// column aggregates (no per-kmer row scan needed).  Returns `None` if the
-    /// filter requires row-level inspection.
+    /// filter requires row-level inspection — always the case for a filter
+    /// that needs the k-mer's sequence, since a `FilterMask` only expresses
+    /// per-genome column aggregates, never per-slot sequence data.
     ///
     /// `threshold` semantics in the returned mask use `>= threshold`, matching
     /// [`obicompactvec::MatrixGroupOps`].  Implementations must add 1 to any
@@ -23,8 +30,13 @@ pub trait KmerFilter: Send + Sync {
 
 /// True when `row` passes every filter in `filters`.
 /// Returns `true` if `filters` is empty.
-pub fn passes_all(filters: &[Box<dyn KmerFilter>], row: &[u32], n_genomes: usize) -> bool {
-    filters.iter().all(|f| f.passes(row, n_genomes))
+pub fn passes_all(
+    filters: &[Box<dyn KmerFilter>],
+    kmer: CanonicalKmer,
+    row: &[u32],
+    n_genomes: usize,
+) -> bool {
+    filters.iter().all(|f| f.passes(kmer, row, n_genomes))
 }
 
 // ── Quorum filters ─────────────────────────────────────────────────────────────
@@ -40,7 +52,7 @@ pub struct MinGenomeFraction {
 }
 
 impl KmerFilter for MinGenomeFraction {
-    fn passes(&self, row: &[u32], n_genomes: usize) -> bool {
+    fn passes(&self, _kmer: CanonicalKmer, row: &[u32], n_genomes: usize) -> bool {
         let p = present_count(row, self.threshold);
         p as f64 / n_genomes as f64 >= self.frac
     }
@@ -63,7 +75,7 @@ pub struct MaxGenomeFraction {
 }
 
 impl KmerFilter for MaxGenomeFraction {
-    fn passes(&self, row: &[u32], n_genomes: usize) -> bool {
+    fn passes(&self, _kmer: CanonicalKmer, row: &[u32], n_genomes: usize) -> bool {
         let p = present_count(row, self.threshold);
         p as f64 / n_genomes as f64 <= self.frac
     }
@@ -86,7 +98,7 @@ pub struct MinGenomeCount {
 }
 
 impl KmerFilter for MinGenomeCount {
-    fn passes(&self, row: &[u32], _n_genomes: usize) -> bool {
+    fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool {
         present_count(row, self.threshold) >= self.count
     }
 
@@ -107,7 +119,7 @@ pub struct MaxGenomeCount {
 }
 
 impl KmerFilter for MaxGenomeCount {
-    fn passes(&self, row: &[u32], _n_genomes: usize) -> bool {
+    fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool {
         present_count(row, self.threshold) <= self.count
     }
 
@@ -129,7 +141,7 @@ pub struct MinTotalCount {
 }
 
 impl KmerFilter for MinTotalCount {
-    fn passes(&self, row: &[u32], _n_genomes: usize) -> bool {
+    fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool {
         row.iter().sum::<u32>() >= self.total
     }
 
@@ -147,7 +159,7 @@ pub struct MaxTotalCount {
 }
 
 impl KmerFilter for MaxTotalCount {
-    fn passes(&self, row: &[u32], _n_genomes: usize) -> bool {
+    fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool {
         row.iter().sum::<u32>() <= self.total
     }
 
@@ -212,7 +224,7 @@ impl GroupQuorumFilter {
 }
 
 impl KmerFilter for GroupQuorumFilter {
-    fn passes(&self, row: &[u32], _n_genomes: usize) -> bool {
+    fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool {
         if !self.ingroup_idx.is_empty() {
             let n = self.ingroup_idx.iter()
                 .filter(|&&i| row.get(i).copied().unwrap_or(0) > self.threshold)
@@ -260,3 +272,27 @@ impl KmerFilter for GroupQuorumFilter {
         Some(FilterMask::And(parts))
     }
 }
+
+// ── Complexity filter (post-hoc, sequence-based) ──────────────────────────────
+
+/// Reject k-mers with normalized entropy below `theta` — the same complexity
+/// metric `obikmer index`'s `--theta`/`--level-max` apply *during* superkmer
+/// construction (see [`obikentropy::KmerEntropy`]), applied here after the
+/// fact, to k-mers already committed to a built index.
+///
+/// Unlike every other filter in this module, this one needs the k-mer's own
+/// sequence, not its per-genome row — `column_mask_expr` is never overridden
+/// (stays `None`), so this filter always forces the row-level scan path in
+/// `rebuild_layer.rs` (which reconstructs the sequence from `unitigs.bin`
+/// regardless, so no extra I/O beyond what filtering already requires).
+pub struct MinComplexity {
+    pub level_max: usize,
+    pub theta: f64,
+}
+
+impl KmerFilter for MinComplexity {
+    fn passes(&self, kmer: CanonicalKmer, _row: &[u32], _n_genomes: usize) -> bool {
+        use obikentropy::KmerEntropy;
+        kmer.entropy(self.level_max) >= self.theta
+    }
+}

src/obikpartitionner/src/rebuild_layer.rs

@@ -126,7 +126,7 @@ fn iter_src_kmers_masked(
                 Some(m) => m.get(slot),
                 None => {
                     let row = src_data.fill_row_by_slot(slot, n_genomes);
-                    filters.iter().all(|f| f.passes(&row, n_genomes))
+                    filters.iter().all(|f| f.passes(kmer, &row, n_genomes))
                 }
             };
             if passes { cb(kmer); }
@@ -165,7 +165,7 @@ fn iter_src_layers(
                 cb(kmer, row.into_boxed_slice());
             } else {
                 let row = src_data.fill_row_by_slot(slot, n_genomes);
-                if filters.iter().all(|f| f.passes(&row, n_genomes)) {
+                if filters.iter().all(|f| f.passes(kmer, &row, n_genomes)) {
                     cb(kmer, row.into_boxed_slice());
                 }
             }

src/obikseq/src/kmer.rs

@@ -341,6 +341,27 @@ impl<L: KmerLength> CanonicalKmerOf<L> {
         ]
     }
 
+    /// Return the four central canonical neighbours (each already canonical),
+    /// substituting the base at the middle position `m = (L::len()-1)/2`
+    /// (well-defined for odd `L::len()`). Each of the 4 substitutions is
+    /// canonicalised independently — this correctly handles the case where a
+    /// substitution flips the canonical orientation, unlike inferring the
+    /// variant from a fixed-orientation flank key. One of the 4 equals
+    /// `self`'s own canonical form (the identity substitution); callers that
+    /// only want the 3 genuine variants should skip it.
+    pub fn central_canonical_neighbors(&self) -> [CanonicalKmerOf<L>; 4] {
+        let k = L::len();
+        let m = (k - 1) / 2;
+        let shift = KMER_BITS - 2 - 2 * m;
+        let cleared = self.0 & !((0b11 as RawKmer) << shift);
+        [
+            KmerOf::<L>(cleared | ((0 as RawKmer) << shift), PhantomData).canonical(),
+            KmerOf::<L>(cleared | ((1 as RawKmer) << shift), PhantomData).canonical(),
+            KmerOf::<L>(cleared | ((2 as RawKmer) << shift), PhantomData).canonical(),
+            KmerOf::<L>(cleared | ((3 as RawKmer) << shift), PhantomData).canonical(),
+        ]
+    }
+
     /// Return the inner value as a raw [`KmerOf<L>`].
     #[inline]
     pub fn into_kmer(self) -> KmerOf<L> {

src/obikseq/src/tests/kmer.rs

@@ -210,4 +210,46 @@ mod tests {
         check!(31);
         check!(32);
     }
+
+    // ── central_canonical_neighbors ─────────────────────────────────────────
+
+    #[test]
+    fn central_canonical_neighbors_hand_checked_k3() {
+        // k=3, centre = index 1. For "ACG", every one of the 4 central
+        // substitutions ("AAG","ACG","AGG","ATG") happens to stay in forward
+        // orientation when canonicalised (verified by hand: each is already
+        // lexicographically <= its own reverse complement), so this case
+        // exercises the substitution logic without the RC-flip edge case.
+        let ck = KmerOf::<ConstLen<3>>::from_ascii(b"ACG").unwrap().canonical();
+        let neighbours = ck.central_canonical_neighbors();
+        let ascii: Vec<Vec<u8>> = neighbours.iter().map(|n| n.to_ascii()).collect();
+        assert_eq!(ascii, vec![b"AAG".to_vec(), b"ACG".to_vec(), b"AGG".to_vec(), b"ATG".to_vec()]);
+        // The identity substitution (centre unchanged) must reproduce `ck`.
+        assert!(neighbours.contains(&ck));
+    }
+
+    #[test]
+    fn central_canonical_neighbors_identity_present_for_various_k() {
+        macro_rules! check {
+            ($n:expr) => {{
+                let ck = KmerOf::<ConstLen<$n>>::from_ascii(&make_seq::<$n>())
+                    .unwrap()
+                    .canonical();
+                let neighbours = ck.central_canonical_neighbors();
+                assert!(
+                    neighbours.contains(&ck),
+                    "identity substitution missing from central_canonical_neighbors for k={}",
+                    $n
+                );
+                // Every returned neighbour must itself already be canonical.
+                for n in &neighbours {
+                    assert_eq!(n.into_kmer().canonical(), *n, "neighbour not canonical for k={}", $n);
+                }
+            }};
+        }
+        check!(1);
+        check!(3);
+        check!(5);
+        check!(31);
+    }
 }

src/obiskbuilder/Cargo.toml

@@ -7,7 +7,13 @@ edition = "2024"
 obikseq = { path = "../obikseq" }
 obikrope = { path = "../obikrope" }
 obiread = { path = "../obiread" }
+obikentropy = { path = "../obikentropy" }
 lazy_static = "1.5.0"
 
 [dev-dependencies]
 obikseq = { path = "../obikseq", features = ["test-utils"] }
+criterion2 = { version = "3", features = ["cargo_bench_support"] }
+
+[[bench]]
+name = "superkmer_stream"
+harness = false

src/obiskbuilder/benches/superkmer_stream.rs

@@ -0,0 +1,58 @@
+//! Throughput of the streaming superkmer pipeline (`RollingStat`'s hot path:
+//! minimizer selection + entropy tracking fused in a single pass).
+//!
+//! Reference point for the `obikentropy` extraction: the entropy bookkeeping
+//! that used to live inline in `RollingStat` was pulled out into a composed
+//! `EntropyTracker`. This benchmark is run before and after that change to
+//! confirm no regression.
+
+use criterion::{Criterion, Throughput, criterion_group, criterion_main};
+use obikrope::Rope;
+use obiskbuilder::SuperKmerIter;
+
+const K: usize = 21;
+const M: usize = 9;
+const LEVEL_MAX: usize = 6;
+const THETA: f64 = 0.7;
+const SEQ_LEN: usize = 200_000;
+
+/// Deterministic pseudo-random ACGT sequence — high enough complexity that
+/// the entropy filter rarely rejects, so the bench stays on the steady-state
+/// path rather than repeatedly resetting.
+fn make_sequence(len: usize) -> Vec<u8> {
+    let mut state: u64 = 0x9E3779B97F4A7C15;
+    (0..len)
+        .map(|_| {
+            state ^= state << 13;
+            state ^= state >> 7;
+            state ^= state << 17;
+            b"ACGT"[(state % 4) as usize]
+        })
+        .collect()
+}
+
+fn make_rope(seq: &[u8]) -> Rope {
+    let mut rope = Rope::new(None);
+    rope.push(seq.to_vec());
+    rope
+}
+
+fn bench_build_superkmers(c: &mut Criterion) {
+    obikseq::set_k(K);
+    obikseq::set_m(M);
+
+    let seq = make_sequence(SEQ_LEN);
+    let rope = make_rope(&seq);
+
+    let mut group = c.benchmark_group("build_superkmers");
+    group.throughput(Throughput::Bytes(SEQ_LEN as u64));
+    group.bench_function("stream", |b| {
+        b.iter(|| {
+            SuperKmerIter::new(std::hint::black_box(&rope), K, LEVEL_MAX, THETA).count()
+        });
+    });
+    group.finish();
+}
+
+criterion_group!(benches, bench_build_superkmers);
+criterion_main!(benches);

src/obiskbuilder/src/entropy_table.rs

@@ -1,108 +0,0 @@
-pub(crate) const NORMK1: [u64; 4] = build_normalized_kmer::<4>();
-pub(crate) const NORMK2: [u64; 16] = build_normalized_kmer::<16>();
-pub(crate) const NORMK3: [u64; 64] = build_normalized_kmer::<64>();
-pub(crate) const NORMK4: [u64; 256] = build_normalized_kmer::<256>();
-pub(crate) const NORMK5: [u64; 1024] = build_normalized_kmer::<1024>();
-pub(crate) const NORMK6: [u64; 4096] = build_normalized_kmer::<4096>();
-
-include!(concat!(env!("OUT_DIR"), "/ln_class_tables.rs"));
-
-const fn normalize_circular(kmer: u64, ws: usize) -> u64 {
-    let mask = (1u64 << (ws * 2)) - 1;
-    let mut canonical = kmer & mask;
-    let mut current = canonical;
-    let mut i = 0;
-    while i < (ws - 1) {
-        let top = (current >> ((ws - 1) * 2)) & 3;
-        current = ((current << 2) | top) & mask;
-        if current < canonical {
-            canonical = current;
-        }
-        i += 1;
-    }
-    canonical
-}
-
-const fn build_normalized_kmer<const N: usize>() -> [u64; N] {
-    let mut result = [0u64; N];
-    let k = k_from_n::<N>();
-    let shift = 64 - k * 2;
-    let mut i = 0;
-    while i < N {
-        let la = (i as u64) << shift;
-        let ra = i as u64;
-        let rc_ra = revcomp_raw(la, k) >> shift;
-        let circ = normalize_circular(ra, k);
-        let circ_rc = normalize_circular(rc_ra, k);
-        result[i] = if circ < circ_rc { circ } else { circ_rc };
-        i += 1;
-    }
-    result
-}
-
-const fn revcomp_raw(x: u64, k: usize) -> u64 {
-    let x = !x;
-    let x = x.swap_bytes();
-    let x = ((x >> 4) & 0x0F0F0F0F0F0F0F0F) | ((x & 0x0F0F0F0F0F0F0F0F) << 4);
-    let x = ((x >> 2) & 0x3333333333333333) | ((x & 0x3333333333333333) << 2);
-    x << (64 - 2 * k)
-}
-
-const fn k_from_n<const N: usize>() -> usize {
-    match N {
-        4 => 1,
-        16 => 2,
-        64 => 3,
-        256 => 4,
-        1024 => 5,
-        4096 => 6,
-        _ => panic!("N must be a power of 4"),
-    }
-}
-
-pub(crate) const WS_MAX: usize = 6;
-
-#[inline(always)]
-pub(crate) const fn n_log_n(n: usize) -> f64 {
-    N_LOG_N[n]
-}
-
-#[inline(always)]
-pub(crate) const fn emax(k: usize, ws: usize) -> f64 {
-    EMAX[k][ws]
-}
-
-#[inline(always)]
-pub(crate) const fn log_nwords(k: usize, ws: usize) -> f64 {
-    LOG_NWORDS[k][ws]
-}
-
-#[inline(always)]
-pub(crate) const fn entropy_norm_kmer<const LEFT: bool, const K: usize>(kmer: u64) -> u64 {
-    const SHIFT: [usize; 7] = [0, 62, 60, 58, 56, 54, 52];
-    const NORM: [&[u64]; 7] = [&[], &NORMK1, &NORMK2, &NORMK3, &NORMK4, &NORMK5, &NORMK6];
-
-    let shift = SHIFT[K];
-    let ra = if LEFT { kmer >> shift } else { kmer };
-    let canonical_ra = NORM[K][ra as usize];
-    if LEFT {
-        canonical_ra << shift
-    } else {
-        canonical_ra
-    }
-}
-
-#[inline(always)]
-pub(crate) const fn ln_class_size<const LEFT: bool, const K: usize>(kmer: u64) -> f64 {
-    const SHIFT: [usize; 7] = [0, 62, 60, 58, 56, 54, 52];
-    let ra = if LEFT { kmer >> SHIFT[K] } else { kmer };
-    match K {
-        1 => LN_CLASS1[ra as usize],
-        2 => LN_CLASS2[ra as usize],
-        3 => LN_CLASS3[ra as usize],
-        4 => LN_CLASS4[ra as usize],
-        5 => LN_CLASS5[ra as usize],
-        6 => LN_CLASS6[ra as usize],
-        _ => panic!("k must be 1..=6"),
-    }
-}

src/obiskbuilder/src/iter.rs

@@ -149,157 +149,5 @@ impl Iterator for SuperKmerIter<'_> {
 // ── tests ─────────────────────────────────────────────────────────────────────
 
 #[cfg(test)]
-mod tests {
-    use super::*;
-    use obikrope::Rope;
-
-    fn setup() {
-        obikseq::params::set_k(K);
-        obikseq::params::set_m(5);
-    }
-
-    fn make_rope(data: &[u8]) -> Rope {
-        let mut r = Rope::new(None);
-        r.push(data.to_vec());
-        r
-    }
-
-    fn run_nofilter(data: &[u8], k: usize) -> Vec<Vec<u8>> {
-        let rope = make_rope(data);
-        SuperKmerIter::new(&rope, k, 1, 0.0)
-            .map(|rsk| rsk.superkmer().to_ascii())
-            .collect()
-    }
-
-    // k=11, m=5 — valeurs minimales du projet (k ∈ [11,31])
-    const K: usize = 11;
-
-    /// Collect the set of canonical k-mers from a raw ASCII sequence (no NUL).
-    fn direct_canonical_kmers(seq: &[u8]) -> std::collections::HashSet<Vec<u8>> {
-        (0..seq.len().saturating_sub(K - 1))
-            .map(|i| obikseq::SuperKmer::from_ascii(&seq[i..i + K]).to_ascii())
-            .collect()
-    }
-
-    /// Collect the set of canonical k-mers emitted by SuperKmerIter over a rope.
-    fn iter_canonical_kmers(rope: &Rope) -> std::collections::HashSet<Vec<u8>> {
-        SuperKmerIter::new(rope, K, 1, 0.0)
-            .flat_map(|rsk| {
-                rsk.superkmer()
-                    .iter_canonical_kmers()
-                    .map(|km| km.to_ascii())
-                    .collect::<Vec<_>>()
-            })
-            .collect()
-    }
-
-    #[test]
-    fn coverage_single_segment() {
-        setup();
-        let seq = b"ACGTACGTACGTACGTACGT";
-        let rope = make_rope(&[seq.as_ref(), b"\x00"].concat());
-        let direct = direct_canonical_kmers(seq);
-        let from_iter = iter_canonical_kmers(&rope);
-        let missing: Vec<_> = direct.difference(&from_iter).collect();
-        assert!(
-            missing.is_empty(),
-            "k-mers perdus dans segment unique : {missing:?}"
-        );
-    }
-
-    #[test]
-    fn coverage_two_segments() {
-        setup();
-        let seg1 = b"ACGTACGTACGTACGTACGT";
-        let seg2 = b"TGCATGCATGCATGCATGCA";
-        let rope = make_rope(&[seg1.as_ref(), b"\x00", seg2.as_ref(), b"\x00"].concat());
-        let mut direct = direct_canonical_kmers(seg1);
-        direct.extend(direct_canonical_kmers(seg2));
-        let from_iter = iter_canonical_kmers(&rope);
-        let missing: Vec<_> = direct.difference(&from_iter).collect();
-        assert!(
-            missing.is_empty(),
-            "k-mers perdus dans deux segments : {missing:?}"
-        );
-    }
-
-    #[test]
-    fn coverage_minimizer_boundary() {
-        setup();
-        // sequence assez longue pour forcer plusieurs changements de minimiseur
-        let seq: Vec<u8> = (0..80).map(|i| b"ACGT"[i % 4]).collect();
-        let rope = make_rope(&[seq.as_slice(), b"\x00"].concat());
-        let direct = direct_canonical_kmers(&seq);
-        let from_iter = iter_canonical_kmers(&rope);
-        let missing: Vec<_> = direct.difference(&from_iter).collect();
-        assert!(
-            missing.is_empty(),
-            "k-mers perdus à la frontière de minimiseur : {missing:?}"
-        );
-    }
-
-    #[test]
-    fn single_segment_one_superkmer() {
-        setup();
-        let out = run_nofilter(b"ACGTACGTACGTACGTACGT\x00", K);
-        assert!(!out.is_empty());
-        let total: Vec<u8> = out.into_iter().flatten().collect();
-        assert!(total.len() >= K);
-    }
-
-    #[test]
-    fn segment_shorter_than_k_emits_nothing() {
-        setup();
-        let out = run_nofilter(b"ACGTACGT\x00", K);
-        assert_eq!(out, Vec::<Vec<u8>>::new());
-    }
-
-    #[test]
-    fn empty_input_emits_nothing() {
-        setup();
-        let out = run_nofilter(b"", K);
-        assert_eq!(out, Vec::<Vec<u8>>::new());
-    }
-
-    #[test]
-    fn two_segments_both_emitted() {
-        setup();
-        let out = run_nofilter(b"ACGTACGTACGTACGT\x00TGCATGCATGCATGCA\x00", K);
-        assert!(!out.is_empty());
-    }
-
-    #[test]
-    fn low_complexity_kmer_is_rejected() {
-        setup();
-        let out_pass = run_nofilter(b"AAAAAAAAAAAACGTACGTACGT\x00", K);
-        assert!(!out_pass.is_empty());
-
-        let rope = make_rope(b"AAAAAAAAAAAAAAAAAAAA\x00");
-        let out_reject: Vec<Vec<u8>> = SuperKmerIter::new(&rope, K, 6, 0.9)
-            .map(|rsk| rsk.superkmer().to_ascii())
-            .collect();
-        assert!(out_reject.is_empty());
-    }
-
-    #[test]
-    fn multi_slice_rope() {
-        setup();
-        let data = b"ACGTACGTACGTACGTACGT\x00";
-        let mid = data.len() / 2;
-        let mut rope = Rope::new(None);
-        rope.push(data[..mid].to_vec());
-        rope.push(data[mid..].to_vec());
-        let out: Vec<Vec<u8>> = SuperKmerIter::new(&rope, K, 1, 0.0)
-            .map(|rsk| rsk.superkmer().to_ascii())
-            .collect();
-        assert!(!out.is_empty());
-    }
-
-    #[test]
-    fn yields_minimizer_value() {
-        setup();
-        let rope = make_rope(b"ACGTACGTACGTACGTACGT\x00");
-        let results: Vec<RoutableSuperKmer> = SuperKmerIter::new(&rope, K, 1, 0.0).collect();
-        assert!(!results.is_empty());
-    }
-}
+#[path = "tests/iter.rs"]
+mod tests;

src/obiskbuilder/src/lib.rs

@@ -10,8 +10,8 @@ pub mod stream_iter;
 mod scratch;
 
 pub(crate) mod encoding;
-pub(crate) mod entropy_table;
-pub(crate) mod rolling_stat;
+#[allow(missing_docs)]
+pub mod rolling_stat;
 
 pub use iter::SuperKmerIter;
 pub use scratch::SuperKmerScratch;

src/obiskbuilder/src/rolling_stat.rs

@@ -1,8 +1,8 @@
+use obikentropy::EntropyTracker;
 use obikseq::kmer::{Minimizer, hash_kmer};
 use obikseq::params;
 
 use crate::encoding::encode_nuc;
-use crate::entropy_table::{WS_MAX, emax, entropy_norm_kmer, ln_class_size, log_nwords, n_log_n};
 
 // ── Stack-allocated ring buffer ───────────────────────────────────────────────
 
@@ -83,33 +83,19 @@ pub struct RollingStat {
     entropy_max_k: usize,
     k: usize,
     m: usize,
-    steady: bool,
     rolling_k: u64,
     rolling_rck: u64,
     k_mask: u64,
     m_mask: u64,
     received: usize,
 
-    // Sliding-window queues — stack-allocated, capacity ≤ k ≤ 31.
-    k1q: Ring<u64, 32>,
-    k2q: Ring<u64, 32>,
-    k3q: Ring<u64, 32>,
-    k4q: Ring<u64, 32>,
-    k5q: Ring<u64, 32>,
-    k6q: Ring<u64, 32>,
+    // Minimizer selection state.
     minimier: Ring<MmerItem, 32>,
 
-    // Frequency count arrays.
-    // Max count per cell ≤ k ≤ 31 → u8 is sufficient.
-    k1c: [u8; 4],
-    k2c: [u8; 16],
-    k3c: [u8; 64],
-    k4c: [u8; 256],
-    k5c: [u8; 1024],
-    k6c: [u8; 4096],
-
-    sum_f_log_f: [f64; WS_MAX + 1],
-    sum_f_log_s: [f64; WS_MAX + 1],
+    // Entropy tracking, composed as a plain inline field so both concerns
+    // update in the same streaming pass without being conflated in one
+    // struct — see `obikentropy::EntropyTracker`.
+    entropy: EntropyTracker,
 }
 
 impl RollingStat {
@@ -120,27 +106,13 @@ impl RollingStat {
             entropy_max_k,
             k,
             m,
-            steady: false,
             rolling_k: 0,
             rolling_rck: 0,
             k_mask: (!0u64) >> (64 - k * 2),
             m_mask: (!0u64) >> (64 - m * 2),
             received: 0,
-            k1q: Ring::new(),
-            k2q: Ring::new(),
-            k3q: Ring::new(),
-            k4q: Ring::new(),
-            k5q: Ring::new(),
-            k6q: Ring::new(),
             minimier: Ring::new(),
-            k1c: [0; 4],
-            k2c: [0; 16],
-            k3c: [0; 64],
-            k4c: [0; 256],
-            k5c: [0; 1024],
-            k6c: [0; 4096],
-            sum_f_log_f: [0.0; WS_MAX + 1],
-            sum_f_log_s: [0.0; WS_MAX + 1],
+            entropy: EntropyTracker::new(k),
         }
     }
 
@@ -148,54 +120,9 @@ impl RollingStat {
         self.rolling_k = 0;
         self.rolling_rck = 0;
         self.received = 0;
-        self.steady = false;
 
-        //        for i in self.k1q.iter() { self.k1c[i as usize] = 0; }
-        //        for i in self.k2q.iter() { self.k2c[i as usize] = 0; }
-        //        for i in self.k3q.iter() { self.k3c[i as usize] = 0; }
-        //        for i in self.k4q.iter() { self.k4c[i as usize] = 0; }
-        //        for i in self.k5q.iter() { self.k5c[i as usize] = 0; }
-        //        for i in self.k6q.iter() { self.k6c[i as usize] = 0; }
-
-        self.k1c.fill(0);
-        self.k2c.fill(0);
-        self.k3c.fill(0);
-        self.k4c.fill(0);
-        self.k5c.fill(0);
-        self.k6c.fill(0);
-
-        self.k1q.clear();
-        self.k2q.clear();
-        self.k3q.clear();
-        self.k4q.clear();
-        self.k5q.clear();
-        self.k6q.clear();
         self.minimier.clear();
-
-        self.sum_f_log_f = [0.0; WS_MAX + 1];
-        self.sum_f_log_s = [0.0; WS_MAX + 1];
-    }
-
-    #[inline]
-    fn update_sums_decrement<const K: usize>(
-        sum_f_log_f: &mut [f64; WS_MAX + 1],
-        sum_f_log_s: &mut [f64; WS_MAX + 1],
-        canonical: u64,
-        f: usize,
-    ) {
-        sum_f_log_f[K] += n_log_n(f - 1) - n_log_n(f);
-        sum_f_log_s[K] -= ln_class_size::<false, K>(canonical);
-    }
-
-    #[inline]
-    fn update_sums_increment<const K: usize>(
-        sum_f_log_f: &mut [f64; WS_MAX + 1],
-        sum_f_log_s: &mut [f64; WS_MAX + 1],
-        canonical: u64,
-        g: usize,
-    ) {
-        sum_f_log_f[K] += n_log_n(g + 1) - n_log_n(g);
-        sum_f_log_s[K] += ln_class_size::<false, K>(canonical);
+        self.entropy.reset();
     }
 
     pub fn push(&mut self, nuc: u8) {
@@ -209,13 +136,6 @@ impl RollingStat {
         self.rolling_rck =
             ((self.rolling_rck >> 2) | ((cnuc as u64) << ((k - 1) * 2))) & self.k_mask;
 
-        let canonical_k1 = entropy_norm_kmer::<false, 1>(self.rolling_k & 3);
-        let canonical_k2 = entropy_norm_kmer::<false, 2>(self.rolling_k & 15);
-        let canonical_k3 = entropy_norm_kmer::<false, 3>(self.rolling_k & 63);
-        let canonical_k4 = entropy_norm_kmer::<false, 4>(self.rolling_k & 255);
-        let canonical_k5 = entropy_norm_kmer::<false, 5>(self.rolling_k & 1023);
-        let canonical_k6 = entropy_norm_kmer::<false, 6>(self.rolling_k & 4095);
-
         self.received += 1;
 
         if self.received >= m {
@@ -248,153 +168,7 @@ impl RollingStat {
             }
         }
 
-        if self.received > k {
-            let old1 = self.k1q.pop_front();
-            let f1 = self.k1c[old1 as usize] as usize;
-            Self::update_sums_decrement::<1>(
-                &mut self.sum_f_log_f,
-                &mut self.sum_f_log_s,
-                old1,
-                f1,
-            );
-            self.k1c[old1 as usize] -= 1;
-
-            let old2 = self.k2q.pop_front();
-            let f2 = self.k2c[old2 as usize] as usize;
-            Self::update_sums_decrement::<2>(
-                &mut self.sum_f_log_f,
-                &mut self.sum_f_log_s,
-                old2,
-                f2,
-            );
-            self.k2c[old2 as usize] -= 1;
-
-            let old3 = self.k3q.pop_front();
-            let f3 = self.k3c[old3 as usize] as usize;
-            Self::update_sums_decrement::<3>(
-                &mut self.sum_f_log_f,
-                &mut self.sum_f_log_s,
-                old3,
-                f3,
-            );
-            self.k3c[old3 as usize] -= 1;
-
-            let old4 = self.k4q.pop_front();
-            let f4 = self.k4c[old4 as usize] as usize;
-            Self::update_sums_decrement::<4>(
-                &mut self.sum_f_log_f,
-                &mut self.sum_f_log_s,
-                old4,
-                f4,
-            );
-            self.k4c[old4 as usize] -= 1;
-
-            let old5 = self.k5q.pop_front();
-            let f5 = self.k5c[old5 as usize] as usize;
-            Self::update_sums_decrement::<5>(
-                &mut self.sum_f_log_f,
-                &mut self.sum_f_log_s,
-                old5,
-                f5,
-            );
-            self.k5c[old5 as usize] -= 1;
-
-            let old6 = self.k6q.pop_front();
-            let f6 = self.k6c[old6 as usize] as usize;
-            Self::update_sums_decrement::<6>(
-                &mut self.sum_f_log_f,
-                &mut self.sum_f_log_s,
-                old6,
-                f6,
-            );
-            self.k6c[old6 as usize] -= 1;
-        }
-
-        if self.steady {
-            let g1 = self.k1c[canonical_k1 as usize] as usize;
-            Self::update_sums_increment::<1>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k1, g1);
-            self.k1c[canonical_k1 as usize] += 1;
-            self.k1q.push_back(canonical_k1);
-
-            let g2 = self.k2c[canonical_k2 as usize] as usize;
-            Self::update_sums_increment::<2>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k2, g2);
-            self.k2c[canonical_k2 as usize] += 1;
-            self.k2q.push_back(canonical_k2);
-
-            let g3 = self.k3c[canonical_k3 as usize] as usize;
-            Self::update_sums_increment::<3>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k3, g3);
-            self.k3c[canonical_k3 as usize] += 1;
-            self.k3q.push_back(canonical_k3);
-
-            let g4 = self.k4c[canonical_k4 as usize] as usize;
-            Self::update_sums_increment::<4>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k4, g4);
-            self.k4c[canonical_k4 as usize] += 1;
-            self.k4q.push_back(canonical_k4);
-
-            let g5 = self.k5c[canonical_k5 as usize] as usize;
-            Self::update_sums_increment::<5>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k5, g5);
-            self.k5c[canonical_k5 as usize] += 1;
-            self.k5q.push_back(canonical_k5);
-
-            let g6 = self.k6c[canonical_k6 as usize] as usize;
-            Self::update_sums_increment::<6>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k6, g6);
-            self.k6c[canonical_k6 as usize] += 1;
-            self.k6q.push_back(canonical_k6);
-        } else {
-            self.push_warmup_increments(
-                canonical_k1, canonical_k2, canonical_k3,
-                canonical_k4, canonical_k5, canonical_k6,
-            );
-        }
-    }
-
-    #[cold]
-    #[inline(never)]
-    fn push_warmup_increments(
-        &mut self,
-        canonical_k1: u64, canonical_k2: u64, canonical_k3: u64,
-        canonical_k4: u64, canonical_k5: u64, canonical_k6: u64,
-    ) {
-        let g1 = self.k1c[canonical_k1 as usize] as usize;
-        Self::update_sums_increment::<1>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k1, g1);
-        self.k1c[canonical_k1 as usize] += 1;
-        self.k1q.push_back(canonical_k1);
-
-        if self.received >= 2 {
-            let g2 = self.k2c[canonical_k2 as usize] as usize;
-            Self::update_sums_increment::<2>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k2, g2);
-            self.k2c[canonical_k2 as usize] += 1;
-            self.k2q.push_back(canonical_k2);
-
-            if self.received >= 3 {
-                let g3 = self.k3c[canonical_k3 as usize] as usize;
-                Self::update_sums_increment::<3>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k3, g3);
-                self.k3c[canonical_k3 as usize] += 1;
-                self.k3q.push_back(canonical_k3);
-
-                if self.received >= 4 {
-                    let g4 = self.k4c[canonical_k4 as usize] as usize;
-                    Self::update_sums_increment::<4>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k4, g4);
-                    self.k4c[canonical_k4 as usize] += 1;
-                    self.k4q.push_back(canonical_k4);
-
-                    if self.received >= 5 {
-                        let g5 = self.k5c[canonical_k5 as usize] as usize;
-                        Self::update_sums_increment::<5>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k5, g5);
-                        self.k5c[canonical_k5 as usize] += 1;
-                        self.k5q.push_back(canonical_k5);
-
-                        if self.received >= 6 {
-                            let g6 = self.k6c[canonical_k6 as usize] as usize;
-                            Self::update_sums_increment::<6>(&mut self.sum_f_log_f, &mut self.sum_f_log_s, canonical_k6, g6);
-                            self.k6c[canonical_k6 as usize] += 1;
-                            self.k6q.push_back(canonical_k6);
-                            self.steady = true;
-                        }
-                    }
-                }
-            }
-        }
+        self.entropy.push(self.received, self.rolling_k);
     }
 
     pub fn ready(&self) -> bool {
@@ -422,29 +196,10 @@ impl RollingStat {
             .map(|raw| Minimizer::from_raw_unchecked(raw << (64 - self.m * 2)))
     }
 
-    pub fn entropy(&self, order: usize) -> Option<f64> {
-        if !self.ready() {
-            return None;
-        }
-        let k = self.k;
-        let em = emax(k, order);
-        if em <= 0.0 {
-            return Some(1.0);
-        }
-        let nwords = k - order + 1;
-        let log_nw = log_nwords(k, order);
-        let nw_f = nwords as f64;
-        let h_corr = log_nw + (self.sum_f_log_s[order] - self.sum_f_log_f[order]) / nw_f;
-        Some((h_corr / em).max(0.0))
-    }
-
     pub fn normalized_entropy(&self) -> Option<f64> {
         if !self.ready() {
             return None;
         }
-        let min_e = (1..=self.entropy_max_k)
-            .filter_map(|ws| self.entropy(ws))
-            .fold(f64::MAX, f64::min);
-        Some(if min_e == f64::MAX { 1.0 } else { min_e })
+        Some(self.entropy.normalized_entropy(self.entropy_max_k))
     }
 }

src/obiskbuilder/src/tests/iter.rs

@@ -0,0 +1,152 @@
+use super::*;
+use obikrope::Rope;
+
+fn setup() {
+    obikseq::params::set_k(K);
+    obikseq::params::set_m(5);
+}
+
+fn make_rope(data: &[u8]) -> Rope {
+    let mut r = Rope::new(None);
+    r.push(data.to_vec());
+    r
+}
+
+fn run_nofilter(data: &[u8], k: usize) -> Vec<Vec<u8>> {
+    let rope = make_rope(data);
+    SuperKmerIter::new(&rope, k, 1, 0.0)
+        .map(|rsk| rsk.superkmer().to_ascii())
+        .collect()
+}
+
+// k=11, m=5 — valeurs minimales du projet (k ∈ [11,31])
+const K: usize = 11;
+
+/// Collect the set of canonical k-mers from a raw ASCII sequence (no NUL).
+fn direct_canonical_kmers(seq: &[u8]) -> std::collections::HashSet<Vec<u8>> {
+    (0..seq.len().saturating_sub(K - 1))
+        .map(|i| obikseq::SuperKmer::from_ascii(&seq[i..i + K]).to_ascii())
+        .collect()
+}
+
+/// Collect the set of canonical k-mers emitted by SuperKmerIter over a rope.
+fn iter_canonical_kmers(rope: &Rope) -> std::collections::HashSet<Vec<u8>> {
+    SuperKmerIter::new(rope, K, 1, 0.0)
+        .flat_map(|rsk| {
+            rsk.superkmer()
+                .iter_canonical_kmers()
+                .map(|km| km.to_ascii())
+                .collect::<Vec<_>>()
+        })
+        .collect()
+}
+
+#[test]
+fn coverage_single_segment() {
+    setup();
+    let seq = b"ACGTACGTACGTACGTACGT";
+    let rope = make_rope(&[seq.as_ref(), b"\x00"].concat());
+    let direct = direct_canonical_kmers(seq);
+    let from_iter = iter_canonical_kmers(&rope);
+    let missing: Vec<_> = direct.difference(&from_iter).collect();
+    assert!(
+        missing.is_empty(),
+        "k-mers perdus dans segment unique : {missing:?}"
+    );
+}
+
+#[test]
+fn coverage_two_segments() {
+    setup();
+    let seg1 = b"ACGTACGTACGTACGTACGT";
+    let seg2 = b"TGCATGCATGCATGCATGCA";
+    let rope = make_rope(&[seg1.as_ref(), b"\x00", seg2.as_ref(), b"\x00"].concat());
+    let mut direct = direct_canonical_kmers(seg1);
+    direct.extend(direct_canonical_kmers(seg2));
+    let from_iter = iter_canonical_kmers(&rope);
+    let missing: Vec<_> = direct.difference(&from_iter).collect();
+    assert!(
+        missing.is_empty(),
+        "k-mers perdus dans deux segments : {missing:?}"
+    );
+}
+
+#[test]
+fn coverage_minimizer_boundary() {
+    setup();
+    // sequence assez longue pour forcer plusieurs changements de minimiseur
+    let seq: Vec<u8> = (0..80).map(|i| b"ACGT"[i % 4]).collect();
+    let rope = make_rope(&[seq.as_slice(), b"\x00"].concat());
+    let direct = direct_canonical_kmers(&seq);
+    let from_iter = iter_canonical_kmers(&rope);
+    let missing: Vec<_> = direct.difference(&from_iter).collect();
+    assert!(
+        missing.is_empty(),
+        "k-mers perdus à la frontière de minimiseur : {missing:?}"
+    );
+}
+
+#[test]
+fn single_segment_one_superkmer() {
+    setup();
+    let out = run_nofilter(b"ACGTACGTACGTACGTACGT\x00", K);
+    assert!(!out.is_empty());
+    let total: Vec<u8> = out.into_iter().flatten().collect();
+    assert!(total.len() >= K);
+}
+
+#[test]
+fn segment_shorter_than_k_emits_nothing() {
+    setup();
+    let out = run_nofilter(b"ACGTACGT\x00", K);
+    assert_eq!(out, Vec::<Vec<u8>>::new());
+}
+
+#[test]
+fn empty_input_emits_nothing() {
+    setup();
+    let out = run_nofilter(b"", K);
+    assert_eq!(out, Vec::<Vec<u8>>::new());
+}
+
+#[test]
+fn two_segments_both_emitted() {
+    setup();
+    let out = run_nofilter(b"ACGTACGTACGTACGT\x00TGCATGCATGCATGCA\x00", K);
+    assert!(!out.is_empty());
+}
+
+#[test]
+fn low_complexity_kmer_is_rejected() {
+    setup();
+    let out_pass = run_nofilter(b"AAAAAAAAAAAACGTACGTACGT\x00", K);
+    assert!(!out_pass.is_empty());
+
+    let rope = make_rope(b"AAAAAAAAAAAAAAAAAAAA\x00");
+    let out_reject: Vec<Vec<u8>> = SuperKmerIter::new(&rope, K, 6, 0.9)
+        .map(|rsk| rsk.superkmer().to_ascii())
+        .collect();
+    assert!(out_reject.is_empty());
+}
+
+#[test]
+fn multi_slice_rope() {
+    setup();
+    let data = b"ACGTACGTACGTACGTACGT\x00";
+    let mid = data.len() / 2;
+    let mut rope = Rope::new(None);
+    rope.push(data[..mid].to_vec());
+    rope.push(data[mid..].to_vec());
+    let out: Vec<Vec<u8>> = SuperKmerIter::new(&rope, K, 1, 0.0)
+        .map(|rsk| rsk.superkmer().to_ascii())
+        .collect();
+    assert!(!out.is_empty());
+}
+
+#[test]
+fn yields_minimizer_value() {
+    setup();
+    let rope = make_rope(b"ACGTACGTACGTACGTACGT\x00");
+    let results: Vec<RoutableSuperKmer> = SuperKmerIter::new(&rope, K, 1, 0.0).collect();
+    assert!(!results.is_empty());
+}

src/obisys/src/lib.rs

@@ -202,6 +202,94 @@ fn cgroup_v1_available() -> Option<u64> {
     Some(limit.saturating_sub(used))
 }
 
+// ── CPU parallelism query ────────────────────────────────────────────────────
+
+/// Returns the number of cores this process can actually use concurrently.
+///
+/// `std::thread::available_parallelism()` reads CPU affinity
+/// (`sched_getaffinity`), not the container's CPU quota — a Docker/cgroup
+/// container commonly reports the *host's* full core count this way while
+/// actually being throttled (via `cpu.max`/`cpu.cfs_quota_us`) to a fraction
+/// of a core. Sizing a thread/worker pool off the unthrottled count causes
+/// severe oversubscription: dozens of threads contending for a sliver of
+/// real CPU time, which can look indistinguishable from a hang for minutes
+/// or hours (observed in CI). On Linux, this reads the cgroup CPU quota
+/// first and returns `min(cgroup_quota, host_parallelism)` when a finite
+/// quota is found; falls back to `available_parallelism()` otherwise (same
+/// convention as [`available_memory_bytes`]).
+pub fn effective_parallelism() -> usize {
+    let host = std::thread::available_parallelism().map(|n| n.get()).unwrap_or(1);
+    #[cfg(target_os = "linux")]
+    {
+        if let Some(quota) = cgroup_v2_cpu_quota() {
+            let effective = quota.clamp(1, host);
+            tracing::debug!(host, quota, effective, source = "cgroup v2", "effective_parallelism");
+            return effective;
+        }
+        if let Some(quota) = cgroup_v1_cpu_quota() {
+            let effective = quota.clamp(1, host);
+            tracing::debug!(host, quota, effective, source = "cgroup v1", "effective_parallelism");
+            return effective;
+        }
+    }
+    tracing::debug!(host, effective = host, source = "available_parallelism (no cgroup quota found)", "effective_parallelism");
+    host
+}
+
+/// cgroup v2 (unified hierarchy): reads `cpu.max` ("<quota> <period>", or
+/// "max <period>" when unlimited) for the current process's cgroup, rounded
+/// up to whole cores. Returns `None` if unlimited or on any parse error.
+#[cfg(target_os = "linux")]
+fn cgroup_v2_cpu_quota() -> Option<usize> {
+    let cgroup = std::fs::read_to_string("/proc/self/cgroup").ok()?;
+    let rel = cgroup
+        .lines()
+        .find(|l| l.starts_with("0::"))?
+        .strip_prefix("0::")?
+        .trim();
+    let base = format!("/sys/fs/cgroup{rel}");
+    let raw = std::fs::read_to_string(format!("{base}/cpu.max")).ok()?;
+    let mut parts = raw.split_whitespace();
+    let quota_str = parts.next()?;
+    let period: f64 = parts.next()?.parse().ok()?;
+    if quota_str == "max" {
+        return None; // unlimited
+    }
+    let quota: f64 = quota_str.parse().ok()?;
+    Some((quota / period).ceil().max(1.0) as usize)
+}
+
+/// cgroup v1 (cpu subsystem): reads `cpu.cfs_quota_us`/`cpu.cfs_period_us`,
+/// rounded up to whole cores. Returns `None` if unlimited (quota <= 0) or on
+/// any parse error.
+#[cfg(target_os = "linux")]
+fn cgroup_v1_cpu_quota() -> Option<usize> {
+    let cgroup = std::fs::read_to_string("/proc/self/cgroup").ok()?;
+    let path = cgroup
+        .lines()
+        .find(|l| l.contains(":cpu:") || l.contains(":cpu,cpuacct:"))?
+        .split(':')
+        .nth(2)?;
+    let base = format!("/sys/fs/cgroup/cpu{path}");
+    let quota: i64 = std::fs::read_to_string(format!("{base}/cpu.cfs_quota_us"))
+        .ok()?
+        .trim()
+        .parse()
+        .ok()?;
+    if quota <= 0 {
+        return None; // unlimited
+    }
+    let period: i64 = std::fs::read_to_string(format!("{base}/cpu.cfs_period_us"))
+        .ok()?
+        .trim()
+        .parse()
+        .ok()?;
+    if period <= 0 {
+        return None;
+    }
+    Some(((quota as f64) / (period as f64)).ceil().max(1.0) as usize)
+}
+
 // ── raw helpers ───────────────────────────────────────────────────────────────
 
 fn get_rusage() -> rusage {
@@ -654,9 +742,7 @@ impl fmt::Display for Reporter {
             return Ok(());
         }
 
-        let n_cores = std::thread::available_parallelism()
-            .map(|n| n.get())
-            .unwrap_or(1);
+        let n_cores = effective_parallelism();
 
         // column widths
         let nw = self