Merge pull request 'refactor: centralize CPU core detection using cgroup-aware utility' (#63) from push-lqzukpulzykz into main

Reviewed-on: #63

.gitea/workflows/ci.yml

@@ -1,4 +1,4 @@
-name: CI
+pname: CI
 
 on:
   pull_request:

.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 |
+| `--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 |
+| `--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 |
+| `--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 |
+| `--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:
+This is a plain Shannon entropy over the observed raw-word frequencies.
-
-$$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.
 
 ## 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é
+- `obikentropy/src/table.rs`, `obikentropy/src/tracker.rs` — formule d'entropie et tables de correction petits effectifs
-- `obiskbuilder/src/lib.rs` — application du filtre lors du scatter (phase 1)
+- `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.

docmd/theory/evolutionary_distances.md

@@ -0,0 +1,820 @@
+# Central-position SNP distance (discussion)
+
+Not implemented. Design discussion for a substitution-rate estimator that
+observes SNPs directly from paired-genome k-mer comparison, as an alternative
+to Mash's Poisson-Jaccard inference (see [obicompactvec](../implementation/obicompactvec.md)
+for the implemented Jaccard/Mash distances).
+
+## Motivation
+
+**Primary intent: restrict the comparison to what is actually comparable.**
+Mash's Jaccard is computed over the **union** of both genomes' k-mer content:
+anything not identically shared is folded into a single undifferentiated
+mass, whether the cause is a point substitution, a genuinely absent
+homologous region (lineage-specific content, gene-family expansion, HGT,
+genome-size asymmetry), or a diverged paralogous copy. The model then
+back-infers a single mutation rate from that mass, silently attributing
+non-homology to mutation. The central-SNP approach instead conditions every
+comparison on local, positive evidence of homology: a locus only enters the
+statistic if its `2m` flanking bases (`m = (k-1)/2`) are found intact in
+*both* genomes — genuinely absent or non-homologous content is excluded from
+the comparison entirely (neither numerator nor denominator), rather than
+silently counted as divergence. This is a conditioning on comparability, not
+just a richer summary statistic — see "Statistic and correspondence with
+`shared`" below for how it plays out against genome-size asymmetry and
+diverged gene families, and "Heterozygosity, ploidy, and consensus-assembly
+inputs" for the corresponding paralogy/heterozygosity filter.
+
+**Secondary benefit: access to the substitution's nature.** Because the
+central base of an odd-k window is directly observable once the flanks are
+confirmed conserved, this also yields more than a rate — the
+transition/transversion split — enabling classical corrected distances
+(Jukes-Cantor, Kimura 2-parameter, LogDet) that a single Jaccard scalar
+cannot support.
+
+## Statistic and correspondence with `shared`
+
+A genomic position `p` is covered by `k` overlapping k-mer windows. Requiring
+the substitution to sit at the window's **center** makes exactly one window
+per SNP eligible — a 1:1 correspondence between SNP and center-neighbor k-mer
+pair, avoiding the ~k-fold overcount of an any-position neighbor search.
+
+A locus with a fully conserved `k`-window (flanks **and** center) is an
+exact-shared k-mer at that locus; a locus with conserved flanks but a
+substituted center is a "central SNP". Both count each locus exactly once, in
+matching units:
+
+```
+p_hat[i,j] = SNP[i,j] / (SNP[i,j] + shared[i,j])
+```
+
+`p_hat` is `P(center substituted | 2m flanks conserved)`. `shared[i,j]` here
+is **not** the general-purpose `shared_kmers` matrix used by Jaccard/Mash
+(`--shared-kmers`, `BitPartials::partial_jaccard` /
+`CountPartials::partial_threshold_jaccard`) — that matrix counts raw k-mer
+identity with no per-genome copy-number constraint, whereas `p_hat`'s
+denominator applies the eligibility rule defined below (raw or
+paralogy-filtered). Both `SNP` and `shared` are accumulated by the same
+sweep, from the same per-locus candidate set (source k-mer + 3 variants),
+under the same eligibility rule — see "Locus eligibility" below, and
+"Heterozygosity, ploidy, and consensus-assembly inputs" for why the
+copy-number constraint matters and what it costs.
+
+**Canonical invariance**: for odd k, the central position maps to itself under
+reverse-complement (`m -> k-1-m = m`, base complemented). A transition maps to
+a transition, a transversion to a transversion — the transition/transversion
+split is well-defined in canonical space.
+
+### Definitions: family, and the canonical form of a family
+
+**Family.** The family of a k-mer `x` is the set of (up to) 4 k-mers sharing
+`x`'s `2m` flanking bases, differing only at the central base `m`. Membership
+is a property of the flank pattern, not of `x` itself: any of the 4 possible
+central substitutions belongs to the same family.
+
+**`central_canonical_neighbors()`** (`obikseq`, `CanonicalKmerOf::central_canonical_neighbors`)
+generates all 4 members from any one of them (observed or not), each
+independently canonicalised (`.canonical()`, i.e. `min(kmer, revcomp(kmer))`).
+This independent canonicalisation is necessary because a central substitution
+can flip which orientation is lexicographically smaller — two members of the
+same family can end up canonicalised in *different* orientations. Despite
+that, the **set** of 4 resulting canonical k-mers is invariant: calling
+`central_canonical_neighbors()` on any member of a family — present in the
+index or not — yields the same 4 values. This is relied upon throughout the
+rest of this document.
+
+**Canonical form of a family.** Because orientation can differ member to
+member, "which of the 4 is the reference" cannot be defined relative to
+*whichever member happened to be visited first*, nor relative to the
+minorant (see below) — both are data-dependent (they depend on what is
+actually observed), so using either as the reference would make the
+reference itself vary depending on what happens to be present in a given
+index. Instead: **the canonical form of a family is, by definition, the
+member whose own central base — read in its own already-canonical
+orientation — is `A`.** This is well-defined for every family, computed
+purely from the flank pattern, whether or not that specific member (or any
+member at all) is actually observed anywhere in the index. Concretely: call
+`central_canonical_neighbors()` on any member (observed or not) to get the
+family's 4 canonical forms; the one among them whose own centre nucleotide is
+`A` is the family's canonical form. The other 3 (`C`, `G`, `T`) are labelled
+relative to *that* fixed reference, not relative to the calling member's own
+orientation.
+
+**Consequence for the minorant.** With this fixed A-referenced labelling,
+`minorant` (the smallest raw encoding among the family's *observed* members,
+introduced further below) becomes directly computable rather than needing to
+be tracked as extra state: regenerate the family's 4 canonical forms from
+any member's own k-mer (cheap, no lookup), compare the raw encodings of
+whichever are marked present, and take the smallest. No separate stored bit
+is required — see Step 2b below, where this replaces the earlier
+minorant-bit design.
+
+## Locus eligibility: raw definition vs. paralogy filter
+
+For each k-mer `x` observed in genome A (source, one MPHF slot; the 3
+central-position variants generated as in the sweep below): check whether
+A's locus (flanks fixed) is resolvable in genome B under one of the 4
+central forms.
+
+**Raw / no model.** The locus counts in the denominator iff at least one of
+the 4 forms is present in B; it counts in the numerator iff the form found in
+B differs from A's own. No constraint on A's or B's own copy number at this
+locus. Open question, not resolved: what if **more than one** of the 4 forms
+is present in B simultaneously (ambiguous target — count once arbitrarily,
+count all, or drop)? The stringent filter below sidesteps the question by
+construction rather than answering it.
+
+**Stringent / paralogy-aware.** The locus counts only if exactly one of the
+4 forms is present in A **and** exactly one is present in B (`count == 1` at
+that slot too, when a count index is available, to also exclude same-allele
+duplicates that presence alone cannot see). This drops the raw definition's
+ambiguous-B case automatically, at the cost of also dropping heterozygous
+sites indiscriminately alongside true duplications (see "Heterozygosity,
+ploidy, and consensus-assembly inputs" below).
+
+**Rejected: parsimony-based multiset pairing for multiplicity > 1.** Rather
+than dropping ambiguous loci, pair identical alleles between A and B first
+(0-mutation explanation preferred), then take `min(unmatched_A, unmatched_B)`
+as inferred SNP pairs. Rejected on two grounds: (1) circularity — selecting
+pairs by minimal apparent divergence, then measuring divergence on those same
+pairs, deflates the estimate by construction, not a neutral heuristic; (2)
+the discriminating signal is a single base among 4 possible values, and the
+flanks are *already* guaranteed identical for every candidate by
+construction (that is how the locus was selected) — no information remains
+in a k-mer window to tell which copy in B truly corresponds to which copy in
+A once multiplicity > 1 on either side. Any pairing rule invents a
+correspondence the data cannot support. Multiplicity > 1 is treated as
+non-identifiable, not as a puzzle to solve with a heuristic.
+
+## Multi-genome framing: family as pseudo-alignment column
+
+**Idea.** Instead of resolving locus eligibility and correspondence one
+genome pair at a time, treat a family as a column of a pseudo multiple
+alignment across *all* genomes simultaneously: for each family, each genome
+has either a net single-copy state (`A`/`C`/`G`/`T`, when the genome carries
+exactly one of the 4 forms) or "missing" (`?`, multi-copy or absent). Flank
+conservation (the `2m` bases fixed by construction) supplies positional
+homology for free — the same role a real MSA would play, without alignment
+software, gap penalties, or progressive-alignment approximations. Stacking
+one such column per family, genomes as rows, produces a genuine SNP
+pseudo-alignment matrix, not just a bag of pairwise distances.
+
+**Precedent.** This is the same principle behind reference-free
+k-mer-based phylogenomics tools — SKA (Split K-mer Analysis, Harris 2018) and
+kSNP: split the k-mer around a variable center, use flank identity to call
+homologous columns across arbitrarily many genomes with no reference and no
+MSA step, then feed the resulting pseudo-alignment to standard phylogenetic
+tools. Landing on the same design independently is a good sign, not a
+coincidence.
+
+**Resolves the pairwise-correspondence problem, properly.** The "Rejected:
+parsimony-based multiset pairing" case above failed because, with only two
+genomes' cardinalities to look at, there is no external constraint to justify
+picking one correspondence between leftover alleles over another — `min(a,b)`
+is a lower bound dressed up as a point estimate (see the follow-up discussion
+on Felsenstein-style parsimony inconsistency: minimum-event explanations are
+systematically biased low whenever homoplasy/multiplicity is real, not
+noise-cancelling). With `N` genomes and many families jointly, the same
+question can be answered the way real phylogenetics answers it: ancestral
+state reconstruction / ML mapping over a tree estimated from the whole
+column set. The tree supplies the missing constraint that two isolated
+columns cannot — this is the principled way out, not a heuristic replacement
+for one.
+
+**Relation to what's already implemented.** `KmerIndex::raw_snp_distance`
+already computes, internally, per family, exactly this row — `single_form:
+Vec<Option<u8>>`, one entry per genome, `None` where ambiguous/absent —
+before immediately collapsing it into pairwise `snp[i,j]`/`shared[i,j]`
+tallies. The pivot this section proposes is small at the implementation
+level: stop collapsing early, and surface the per-family row as a first-class
+artifact (a `families x genomes` matrix). Pairwise raw p-distance becomes one
+projection of that matrix (what's computed today), not the primary object;
+downstream, the matrix itself could feed real phylogenetic tools (parsimony/
+ML, e.g. RAxML/IQ-TREE-style) instead of only NJ/UPGMA on a homemade
+pairwise-distance matrix.
+
+**Caveat: column completeness shrinks with `N`.** The probability that a
+family's flanks stay intact simultaneously across all `N` genomes decays with
+`N` (same ascertainment-bias mechanism as Bias 1 above, compounded over more
+genomes) — fully-resolved columns (no `?` anywhere) become rare as more
+genomes are added. Same missing-data situation any real multi-species
+alignment faces, and phylogenetic tools already handle it well; the practical
+implication is that columns should be allowed partial coverage (>=2 resolved
+genomes, not unanimous) rather than requiring every genome to be net
+single-copy at that locus.
+
+## Heterozygosity, ploidy, and consensus-assembly inputs
+
+A within-genome multiplicity signal (more than one of the 4 central forms
+present at a locus) is produced identically by two distinct causes:
+paralogous duplication and diploid/polyploid heterozygosity. K-mer data alone
+cannot distinguish them. The one real discriminator is sequencing depth
+(heterozygous site: total depth of the present forms ~= the genome's
+single-copy average; duplication: ~2x or more) — but that signal only exists
+if genome "counts" are raw-read depth (FASTQ input), not occurrence counts in
+an assembled FASTA, where per-locus depth is not preserved.
+
+**Magnitude is taxon- and mating-system-dependent, not universal.**
+Heterozygosity density: mammals ~1 site / 1-1.5 kb (~0.1%); highly
+outcrossing plants (maize, poplar) reported an order of magnitude higher
+(~1%); self-fertilising plants (*Arabidopsis thaliana*) near zero — but with
+a documented failure mode where segmental duplication masquerades as
+"pseudo-heterozygosity"; fungi split between haploid vegetative stages
+(non-issue) and dikaryotic Basidiomycetes, where two long-diverged haploid
+nuclei coexist without fusing. The estimator's target use case (closely
+related genomes, k=31) is exactly where the stringent filter above costs the
+least for low-heterozygosity taxa and the most for outcrossing/dikaryotic
+ones — no universal threshold; this is a scope caveat to document, not a
+problem to solve generically.
+
+**Why assembled-consensus inputs don't make measured distances wrong.**
+Phylogenetic inputs are near-universally assemblies, not raw reads, and
+assemblers collapse heterozygous sites to one consensus allele per
+position — effectively an arbitrary, largely uncorrelated-between-assemblies
+choice at each het site. This does not inject unbounded noise: standard
+population genetics gives `d_xy = d_a + (pi_A + pi_B)/2` — the expected
+pairwise difference between a random allele of population A and a random
+allele of population B equals the net (fixed) divergence `d_a` plus the
+average of the two populations' own within-population diversity `pi`.
+Consensus flattening realises exactly this random-allele draw, so the
+measured genome-to-genome distance is a `d_xy`-like quantity, not `d_a` —
+inflated by heterozygosity by a well-characterised additive term, not
+distorted unpredictably. The term is negligible when `pi << d_xy` (the common
+case for cross-species comparisons), and becomes material precisely in the
+two cases already flagged above: very closely related genomes (this
+estimator's explicit target) and highly heterozygous outcrossing organisms,
+where `pi` and `d_xy` are the same order of magnitude.
+
+Caveat: this assumes the flattening is uncorrelated with the phylogenetic
+signal — plausible for de novo assembly, not guaranteed for reference-guided
+assembly biased toward one allele (e.g. the reference's) at each het site,
+which would turn the noise term into a systematic bias toward the reference
+lineage. Not evaluated here.
+
+**Forward-looking implication, not part of the current design.** The
+multiplicity > 1 signal discarded by the stringent filter is a crude
+per-genome proxy for `pi` (under low background paralogy). If a `pi_hat` per
+genome were tallied alongside `SnpTally`, a `d_a` correction
+(`p_hat - mean(pi_hat_i, pi_hat_j)/2`, roughly) could recover an estimate
+closer to net divergence instead of `d_xy` — a possible extension, not
+scoped here.
+
+## Sufficient statistic: 4x4 base-pair tally
+
+Tabulating the joint distribution of `(center_i, center_j)` over conserved-flank
+loci, per genome pair, is sufficient for every downstream correction:
+
+| Estimator | Input | Formula |
+|---|---|---|
+| Raw p-distance | total off-diagonal / total | `p = SNP / (SNP + shared)` |
+| Jukes-Cantor | p | `d = -3/4 * ln(1 - 4p/3)` |
+| Kimura 2-parameter | transition rate P, transversion rate Q | `d = 1/2 ln(1/(1-2P-Q)) + 1/4 ln(1/(1-2Q))` |
+| LogDet/paralinear | full 4x4 + base-composition margins | `d ~= -1/4 ln det(F)`, robust to non-stationary base composition |
+
+JC/K2P need only the total and the transition/transversion split (the
+diagonal collapses to a single "shared" total). LogDet needs the full 4x4,
+already populated at no extra cost (see Step 1/2 below).
+
+Memory for the 4x4 tally: `n^2 * 16` counters. Trivial for the project's
+genome-scale use case (tens to hundreds of genomes); ~13 GB at n=10^4 — outside
+scope but worth flagging if n grows.
+
+## Biases (properties of the estimator, not defects)
+
+1. **Conserved-flank ascertainment bias.** Only SNPs with intact `2m`-base
+   flanks are visible; window-intact probability decays as `(1-p)^{2m}`. For
+   k=31 (2m=30): 0.74 at p=1%, 0.21 at p=5%, 0.04 at p=10%. This estimator
+   targets **closely related genomes**. Under rate heterogeneity across sites
+   (universal in practice), conserved flanks correlate with slow centers, so
+   `p_hat` underestimates the genome-wide average rate — it specifically
+   estimates the substitution rate of **conserved regions**.
+   Two distinct factors are at play here, not one: `P(centre of a given
+   window is a SNP) = p` exactly, **independent of k** — a direct restatement
+   of the raw per-site rate via the bijective window<->centre-position
+   correspondence (Statistic section above), not a k-dependent quantity.
+   `(1-p)^{2m}` is the *separate*, genuinely k-dependent ascertainment factor
+   (are the flanks also intact). The two multiply:
+   `P(usable window showing a central SNP) = p * (1-p)^{2m}` — e.g. at
+   p=1/31 (~3.2%), k=31: `p * (1-p)^30 ~= 0.0323 * 0.374 ~= 1.2%`, i.e. about
+   1 window in 83, not 1 in 31 (which is only the centre-mutated fraction,
+   before requiring intact flanks).
+2. **Bias toward isolated SNPs.** Two SNPs within k of each other disqualify
+   each other's flanks. Hypervariable regions are invisible by construction.
+3. **Indels are invisible.** A frameshift destroys k-mer matches in a block;
+   this channel captures substitutions only. Indel divergence shows up as lost
+   shared k-mers (lower Jaccard/Mash), not as SNP signal.
+4. **k-dependent specificity.** "A k-mer match implies common ancestry" is
+   quantitative. For a 3 Gbp genome, expected random flank-30 collisions
+   (k=31): `(3e9)^2 / 4^30 ~= 8` — negligible. At k=21: `(3e9)^2 / 4^20 ~= 2e6`
+   — no longer negligible. k=31 is safe; k<=21 is marginal to unreliable for
+   large genomes. The large k that guarantees homology is the same k that
+   shrinks the detectable-divergence window — an inherent tension.
+
+## Implementation: avoid materializing a de Bruijn graph
+
+A central-SNP pair is topologically a simple bubble in the colored de Bruijn
+graph (source/sink k-mer shared, two length-k branches differing only at the
+midpoint). Classical bubble-calling (Cortex/discoSNP-style) finds these, but
+requires the graph — nodes plus adjacency for ~10^9 colored k-mers — resident
+in memory. **Rejected**: prohibitive RAM for this project's scale.
+
+A naive per-pair generalisation of variant lookup across n genomes (query each
+non-shared k-mer's 3 central variants against every counterpart genome's
+index) costs `O(n^2 . N . 3)` random lookups, with the same k-mer's 3 variants
+regenerated and requeried once per counterpart genome — pure redundant work.
+**Rejected** as the basis for an n-genome design.
+
+## Implementation: sequential per-partition sweep (no scratch, no graph)
+
+`KmerIndex::distance()` already opens every partition's `presence_store`/
+`count_store` simultaneously, memory-mapped, into one `LayeredStore`
+(`distance.rs:73-77`). "Querying another partition" is therefore not a new
+I/O pattern to design — it is the same O(1) MPHF+evidence lookup the `query`
+command already performs at scale. This lets the SNP tally be computed with
+**no scratch files and no auxiliary graph**, by sweeping partitions once each
+as a source:
+
+1. For source partition `p`, enumerate its **distinct** k-mers (one per MPHF
+   slot; each already carries its full multi-genome presence/count vector —
+   no need to explode per (k-mer, genome) occurrence).
+2. For each, generate the 3 central-substitution variants and **canonicalise
+   each independently** (`min(kmer, revcomp)`, exactly as any normal query) —
+   this avoids the orientation edge case a masked-flank grouping would have
+   (a substitution that flips canonical orientation is handled correctly
+   because each variant is canonicalised on its own, not inferred from a
+   fixed-orientation flank key).
+3. Compute each variant's target partition `q` via its minimizer; batch/sort
+   the partition's outgoing variant queries by `q` for locality.
+4. Look up each variant in `q`'s already-mmap'd MPHF+evidence; on a hit,
+   combine the source's presence vector (base `a`) with the variant's
+   presence vector (base `b`): for every `i` carrying `a` and `j` carrying
+   `b`, `tally[i,j][a,b] += 1`.
+
+**Deduplication needs no persisted state.** Sweeping partitions in a fixed
+order `p = 0, 1, ..., P-1` and only acting on a variant when its target
+partition `q >= p` guarantees each unordered SNP pair is counted exactly
+once: a pair with `q < p` was already resolved earlier, when `q` was itself
+the source partition and `p` (being `>= q`) was a valid forward target. No
+cross-partition flag array is needed — the sweep order *is* the
+deduplication rule. Within the same partition (`q == p`), a lightweight
+transient tie-break suffices: either a `#slots(p)`-bit scratch flag reset per
+partition, or simply comparing the two k-mers' raw `u64` encodings and only
+counting when `kmer_source < kmer_variant` — no storage at all.
+
+This is a distinct computation stage, not a `partial_*` in the existing
+additive-by-partition sense: step 3-4 read across partition boundaries by
+construction, unlike the row-local `partial_jaccard`/`partial_threshold_jaccard`
+primitives. But it requires no new index files, permanent or scratch:
+`unitigs.bin`, `mphf.bin`, `evidence.bin`, and the presence/count columns are
+read as-is, and the only extra memory is the current partition's small
+outgoing-query batch (`#kmers(p) * 3`, released once `p` is done) plus the
+persistent `tally` accumulator (`n^2 * 16` counters, see above).
+
+**Outer loop (over source partitions `p`) must stay sequential.** Two
+independent reasons, not just one: (a) memory — the bounded-footprint claim
+above only holds with one partition's outgoing-query batch in flight; running
+`T` source partitions concurrently multiplies that batch by `T`, exactly the
+blowup the design avoids; (b) correctness — the `q >= p` deduplication rule
+requires partitions to be claimed as sources in a fixed order; running `p1 <
+p2` concurrently gives no guarantee `p1` has finished claiming its `q >= p1`
+targets before `p2` starts claiming its own, breaking the "counted exactly
+once" property.
+
+**Inner loop (target-partition lookups for a fixed `p`) parallelises safely.**
+Each lookup is an O(1) read against an already-mmap'd structure, independent
+of the others, with no growing allocation — no memory blowup, no ordering
+dependency between different `q`. The only shared mutable state is `tally`;
+give each worker thread a **thread-local partial tally** (fixed `n^2 * 16`
+size, independent of partition size) and merge into the global `tally` once
+`p`'s inner loop completes — the same reduce-then-merge pattern Rayon already
+uses elsewhere in this codebase to open partitions in parallel. Extra memory:
+`#threads * n^2 * 16`, negligible (~512 MB at n~1000, 32 threads) and
+unrelated to partition size.
+
+**Cost**: `3 * N_distinct` MPHF lookups total across the whole index (each
+partition swept once as source) — the same order of magnitude and the same
+operation as running `query` over the index's entire k-mer content against
+itself, three times. This is the tool's already-optimized regime, not a new
+I/O profile to validate.
+
+## Cheaper: subsampling
+
+Since the target is a ratio, restricting the source-partition sweep to a
+bottom-`s` hash sketch (only enumerate k-mers with `hash < threshold` as
+sources) divides the lookup count by the sampling factor without biasing
+`p_hat`. Mash-like tradeoff: rate estimated from a sample, not the full
+k-mer set.
+
+## Recommendation
+
+Sequential per-partition sweep (Route D): reuse the already-mmap'd
+per-partition MPHF/evidence/presence structures for O(1) variant lookups,
+dedup via the fixed sweep-order rule (`q >= p`, plus an in-partition
+tie-break), no scratch files, no graph materialisation. Both the SNP
+(off-diagonal) and shared (diagonal) counts are accumulated by this same
+sweep, under the locus-eligibility rule chosen (raw or paralogy-filtered) —
+not reused from the general-purpose `shared_kmers` matrix, whose raw-identity
+definition does not apply the same copy-number constraint. Distances (p, JC,
+K2P, LogDet) as finalisations of the resulting 4x4 tally, mirroring the
+`partial_* -> *_dist_matrix` pattern used for Jaccard/Mash/Bray-Curtis/etc.
+
+## Detailed implementation plan
+
+Grounded in the current codebase. File/type references are anchors, not
+prescriptions; adjust to reality when implementing.
+
+### Step 0 — new low-level primitives (`obikseq`)
+
+Two helpers do not yet exist and are prerequisites:
+
+1. **Central neighbours.** `CanonicalKmerOf<L>` already exposes
+   `left_canonical_neighbors()` / `right_canonical_neighbors()`
+   (`obikseq/src/kmer.rs`), each returning the 4 canonicalised neighbours at
+   an end position. Add `central_canonical_neighbors()` returning the 4
+   variants at position `m = (k-1)/2` (each independently canonicalised via
+   `.canonical()`). The 3 that differ from the source are the query variants;
+   skip the identity. Building on `nucleotide(i)` / the raw 2-bit layout keeps
+   it O(1).
+2. **Lone-k-mer minimiser.** Routing a *synthetic* variant to its partition
+   needs its minimiser, but `RollingStat` (`obiskbuilder/src/rolling_stat.rs`)
+   only computes minimisers incrementally along a sequence. Add a standalone
+   `minimizer(kmer) -> Minimizer` that scans the `k-m+1` m-mer windows
+   (`PackedSeq::mmer`, `obikseq/src/packed_seq.rs`), canonicalises each, and
+   takes the min by `seq_hash()` — the same selection `RollingStat` performs,
+   evaluated once. Partition index is then
+   `(minimizer.seq_hash() & (n_partitions - 1)) as usize`, exactly as
+   `QueryBatch::from_records` (`obikmer/src/cmd/query.rs:142`); `n_partitions`
+   is a power of two so the mask is valid.
+
+### Step 1 — the tally accumulator (`obikindex`)
+
+A `SnpTally` holding, per genome pair, the 4x4 joint count of central bases:
+`n * n * 4 * 4` `u64` (or a packed lower-triangular form since it is
+symmetric). Provide `merge(&mut self, other: &SnpTally)` for the thread-local
+reduce, and accessors yielding, per pair `(i,j)`: total off-diagonal (SNP),
+diagonal (shared, i.e. `p_hat`'s denominator minus SNP), transition count
+`P`, transversion count `Q`. The diagonal is always populated — it is not an
+optional LogDet-only extra, since `p_hat`'s denominator is no longer sourced
+from the external `shared_kmers` matrix (see "Locus eligibility" and
+"Statistic and correspondence with `shared`" above): the source k-mer's own
+presence/count vector, already in hand when it is enumerated, supplies the
+diagonal entry directly, at no extra lookup cost.
+
+### Step 2 — the sweep (`obikindex`, new `snp.rs`)
+
+Mirror `distance.rs`: open the presence or count store per partition. But
+instead of a per-partition `partial_*`, run the sequential source sweep:
+
+```text
+for p in 0..n_partitions:                       # OUTER — sequential
+    open source partition p's layers (QueryLayer-style, obikpartitionner)
+    enumerate distinct canonical k-mers of p (one per MPHF slot) with their
+        presence/count vectors                   # column-major, as query stage 2
+    par_iter over these source k-mers:           # INNER — rayon, thread-local tally
+        apply eligibility rule to the source's own vector (raw: none;         # diagonal
+            stringent: exactly one of the 4 forms present in each genome)     # gate
+        for i in genomes eligible with source base a:
+            for j in genomes eligible with source base a:
+                thread_tally[i,j][a,a] += 1        # diagonal — no extra lookup
+        for each of the 3 central variants:
+            q = partition_of(variant)
+            if q < p: continue                    # dedup: forward targets only
+            if q == p and variant <= source.raw(): continue   # in-partition tie-break
+            slot = layers[q].find_slot(variant)   # MphfLayer::find, mmap'd
+            if hit:
+                vb = variant presence/count vector
+                apply eligibility rule to vb (as above)
+                for i in eligible genomes with source base a:
+                    for j in eligible genomes with variant base b:
+                        thread_tally[i,j][a,b] += 1
+    merge thread-local tallies into global SnpTally
+```
+
+The inner lookup is precisely `QueryLayer::find_slot` +
+`col_value(g, slot)` (`obikpartitionner/src/query_layer.rs`) — reuse or factor
+out that path rather than reimplementing MPHF access. Enumerating "all distinct
+k-mers of a partition with their vectors" is the `dump`/`query` stage-2
+column-major scan already implemented in `dump_layer.rs` /
+`query_partition_with`; factor a reusable iterator if none fits.
+
+`presence_threshold` applies exactly as elsewhere: a genome "carries base b"
+iff its count at that slot is `>= presence_threshold` (trivially `>= 1` for
+presence indexes).
+
+### Open problem (unresolved, session end — not yet fully convinced)
+
+The `q >= p` / tie-break dedup rule in Step 2's pseudocode above is **flawed**
+for the stringent (paralogy-filtered) eligibility rule: it only ever brings
+two family members into view at once (the source and one looked-up variant),
+never all four simultaneously, and which subset gets compared depends on
+partition sweep order. "Exactly one of the 4 forms present in genome A" is a
+whole-family property and cannot be decided correctly from a sequence of
+pairwise, order-dependent glimpses — the pseudocode above needs revision, not
+just the eligibility gate bolted onto it as written.
+
+Direction discussed, **not yet settled**:
+
+1. **Every distinct source k-mer looks up all 3 variants unconditionally**
+   (drop the `q < p` skip entirely) so that every observed family member
+   independently gathers all 4 vectors (its own + whichever of the 3
+   variants exist) at once — a whole-family, order-independent view, computed
+   redundantly once per observed member. Same total lookup order of
+   magnitude as already budgeted (`3 * N_distinct`), just organised
+   differently (no lookup actually skipped, versus the original rule which
+   skipped roughly half).
+2. **Tie-break after gathering, not before**: only the member whose own
+   canonical encoding is the smallest *among the members actually observed*
+   (now known, since all were just looked up) writes to `SnpTally`; the
+   others silently discard their redundant computation. Deterministic,
+   order-independent — as a side effect this also removes the "outer loop
+   must stay sequential" constraint from the cost/parallelism discussion
+   above, since no step depends on partition processing order any more.
+3. **Proposed optimisation**: precompute, once at index build time, a
+   compact global (not per-genome) annex per MPHF slot — the count of
+   *other* family members observed anywhere in the dataset (0-3). Slots with
+   count 0 (majority under low divergence and few genomes, but see the
+   scaling caveat below) need no cross-lookup at all: eligibility reduces to
+   a local `count == 1` check at that single slot, and only slots with count
+   >= 1 enter the 3-lookup sweep machinery above. Revised (see Step 2b
+   below): minorant status *is* stored alongside the count after all, on 3
+   bits rather than 2 — it comes for free from the same lookups needed to
+   count siblings, and storing it lets the sweep discard non-minorant slots
+   without re-fetching anything.
+
+**Minorant/sibling-count relationship, worked out precisely.** "Minorant" is
+a one-way implication from sibling count, not an equivalence: `0 siblings
+=> minorant` (trivially — with no other observed member, the k-mer is by
+definition the smallest of the observed set, itself alone), and its
+contrapositive `not minorant => >= 1 sibling`. The converse does not hold:
+being the minorant says nothing about sibling count — a minorant can have 0,
+1, 2 or 3 siblings, all with larger encodings than itself. Consequence: this
+confirms, as a logical necessity rather than a heuristic, that a 0-sibling
+slot can always write its diagonal contribution with zero ambiguity and no
+lookup (it is unconditionally its own minorant) — but it gives no shortcut
+for the >= 1-sibling case, where minorant status still requires the actual
+comparison of gathered encodings; sibling count alone never determines it.
+
+**When to compute the annex, and cache invalidation.** Sibling count is a
+property of the whole set of columns (genomes/groups) currently in the
+index, not of any single genome — it cannot be computed correctly at
+mono-genome build time (a family may gain siblings, or its minorant may
+change, once more genomes are merged in later). Computing it eagerly at
+every `merge` would also waste work on intermediate merged states nobody
+ever queries. Instead: compute it lazily, on first `distance` call against a
+given index, and persist the result alongside that index for subsequent
+calls — the same lazy-derived-cache pattern `PersistentBitMatrix` already
+uses for `Columnar` -> `Packed`. This requires no explicit invalidation for
+`merge` or `filter` (`obikindex/src/merge.rs`, `obikmer/src/cmd/filter.rs`):
+both only ever write to a fresh `--output` directory, never mutate an input
+index in place, so a re-merged/re-filtered index is simply a new state with
+no annex yet. `select --in-place` (`select_layer.rs:139-235`) is the
+exception: it aggregates genome columns into groups (Any/All/None/Sum/Min/
+Max) by mutating the existing index's files without changing its location.
+It does not remove k-mer rows, but it can still change eligibility and
+sibling counts derived from those rows (e.g. a `Sum` over several
+single-copy genomes can read as multi-copy at the group level). Because it
+mutates in place, **`select --in-place` must explicitly invalidate (delete
+or mark stale) any cached sibling-count annex for that index** — the one
+operation in the current pipeline where this doesn't happen for free.
+
+Not yet convinced this is the right shape, and Step 2's pseudocode above has
+not been rewritten to match — flagged for the next pass rather than resolved
+here.
+
+### Step 2b — sibling-count / minorant annex (consolidated plan)
+
+Scope: only the precursor annex — not the SNP tally itself, whose Step 2
+sweep remains unresolved above. This piece is simpler than the sweep,
+because it writes to an independent per-slot value, not a shared
+cross-k-mer accumulator, so it needs no dedup/ownership logic at all at this
+stage.
+
+**Revised annex encoding — 4-bit presence mask, not 3-bit (minorant +
+count).** Superseded after settling the "canonical form of a family"
+definition above. The 3-bit design (1 minorant bit + 2-bit sibling count,
+§ below, kept for the historical record) had two problems: it discards
+*which* variants are present (only how many), so any future consumer
+(the SNP sweep, or a stats pass — see below) that needs to know which bases
+exist still has to regenerate and blindly re-query all 3 candidates; and
+the minorant bit's meaning was tied to whichever member was visited, not to
+a fixed reference. Storing instead a **4-bit mask** — one bit per base
+(A/C/G/T), set iff that member of the family (labelled relative to the
+family's fixed canonical form, i.e. the member with `A` at the centre — see
+above) is observed anywhere in the index — fixes both:
+- **Sibling count is derived, not stored**: `siblings = popcount(mask) - 1`.
+- **Minorant is derived, not stored**: regenerate the family's 4 canonical
+  forms from the slot's own k-mer (cheap, no lookup — see above), compare
+  the raw encodings of whichever bits are set in the mask, take the
+  smallest.
+- **A future consumer knows exactly which variants to (re-)query** —
+  `popcount(mask) - 1` lookups instead of always 3, and it knows *which*
+  3 (or fewer) to issue, not just how many hits to expect.
+- The all-zero value (no base present at all) is still logically
+  unreachable as a real result — the slot's *own* base is always present in
+  its own family — so it remains available as a free "not yet computed"
+  sentinel, exactly as before.
+
+1. **Primitive.** Reuse `central_canonical_neighbors()` from Step 0
+   unchanged — the 3 canonicalised central-substitution variants of a k-mer
+   (plus the identity, i.e. all 4 members of the family — see "Definitions"
+   above).
+2. **New annex type** (`obicompactvec`, alongside `bitmatrix.rs`): a 4-bit-
+   per-slot packed array (the presence mask above), one per partition — same
+   on-disk shape family as `PersistentBitMatrix`'s `Packed` variant, but
+   simpler (no per-genome columns, a single derived read-only value per
+   slot).
+   <details><summary>Superseded 3-bit design (historical)</summary>
+   3 bits, storing minorant status alongside sibling count directly, since
+   it came for free from the same lookups (point 3 below) — 5 real states
+   (not-minorant; minorant with 0/1/2/3 siblings) fit in 3 bits (8 states,
+   3 unused). This let the future SNP sweep discard a non-minorant slot
+   instantly, with no lookup at all. The otherwise-unreachable combination
+   "not-minorant + 0 siblings" (0 siblings always implies minorant) doubled
+   as the "not yet computed" sentinel. Replaced by the 4-bit mask above,
+   which subsumes this benefit (minorant still derivable, now for free at
+   read time rather than stored) while also fixing the "which variant"
+   blindness.
+   </details>
+3. **Computation pass** (`obikindex`, new `siblings.rs`): **one
+   `obipipeline` run per layer, iterated sequentially over the index's
+   layers** — settled after two false starts, worth recording both.
+   - *False start 1*: "fully parallel over every partition/slot at once,
+     no ordering at all". Correctness is fine with this (sibling count and
+     minorant are order-independent, unlike the old `q >= p` dedup they
+     replace), but it reintroduces, at a larger scale, exactly the
+     memory-blowup the original Step 2 sweep's sequential-outer-loop
+     constraint existed to prevent: scattering every source partition at
+     once multiplies the in-flight outgoing-query volume by the number of
+     partitions.
+   - *False start 2*: push the layer loop itself into the pipeline (source
+     = the index's layers, a first `Flat` stage expands each layer into
+     its k-mers). `obipipeline`'s scheduler already bounds memory on its
+     own — it dispatches every item through a **shared** worker pool at
+     each stage boundary (`scheduler.rs:217-372`, `dispatch()` into a
+     common `worker_tx` queue, any free worker picks up any pending item;
+     not "one worker owns a chunk end to end"), with a biased `Select`
+     that prioritises draining items already advanced in the chain over
+     admitting new source items (`scheduler.rs:271-282`: stage results
+     outrank the source, "vider le pipeline en priorité" / "dernier
+     recours" for new data) — so bounded channel `capacity` plus this
+     drain-first bias already caps in-flight work without any external
+     sequential discipline. Correct, but it means k-mers from several
+     layers can be completing concurrently, so the sink would need to
+     track several open per-layer annex-file writers at once — real,
+     avoidable complexity.
+   - **Settled design**: keep the layer loop external and sequential —
+     not for memory (the pipeline's own `capacity`/priority mechanism
+     already provides that, for free, regardless), but so each pipeline
+     run's sink targets exactly one layer's annex file, no concurrent
+     multi-writer bookkeeping. Per layer: source = that layer's distinct
+     k-mers; a `Flat` (1->N) stage generates the 3 central variants of a
+     k-mer, each tagged with its origin (local slot); a transform stage
+     routes each variant to its target partition (unchanged per-k-mer
+     minimiser); a transform stage performs the lookup (existence-only —
+     `find_slot` hit/miss, cheaper than the SNP sweep's full column
+     fetch); a final stage/sink folds each answer into its origin's
+     running state (below) and, once a layer's k-mers are all resolved,
+     flushes the completed array to that layer's annex file. Many small,
+     single-purpose stages on purpose, to let the scheduler interleave
+     them finely across many in-flight items — this deliberately does
+     **not** mirror how `obipipeline` is used elsewhere today: `query.rs`'s
+     `process_chunk` lumps parse+route+query+serialise into one closure
+     (`query.rs:325,743-758`), and `scatter.rs` only pipelines file-
+     reading/superkmer construction, routing partitions afterwards in a
+     plain sequential loop (`KmerPartition::write_batch`,
+     `partition.rs:140`) — both under-use the fine-grained scheduling the
+     mechanism offers, so they are not precedents to copy, only existing
+     (and arguably improvable, out of scope here) usages. Cross-partition
+     lookups (querying another layer's MPHF for a variant) remain
+     necessary as before — only the *output* side is kept single-layer.
+   - **Reconciliation**: processed at the granularity of one *answer batch
+     per destination partition*, not one source k-mer at a time — this is a
+     proper shuffle, not a per-k-mer wait. Each source partition `p` holds a
+     small array of running states `(minorant = true, siblings = 0)`, one
+     per local slot, initialised at scatter time and **persisting across
+     however many destination-partition batches answer it** (up to 3, one
+     per variant, not necessarily all from the same `q`). Every scattered
+     query carries an origin tag (source partition + local slot) so its
+     answer can be routed back. When target partition `q` returns its batch
+     (all answers for every query that named `q`, regardless of which source
+     k-mer or which source partition they came from), that batch is walked
+     once, locally, and each answer updates — via its origin tag — the
+     matching entry in *its* source partition's array: a miss changes
+     nothing; a hit does `siblings += 1`, and if the found sibling's own
+     encoding is smaller than the source's, `minorant = false`. A given
+     source k-mer's state is final only once every destination batch
+     concerning it has been folded in; its partition's array is flushed to
+     the persistent annex once complete. Commutative per entry, so the order
+     in which destination batches arrive and get folded in doesn't matter.
+
+   **Open optimisation, not adopted yet — real tradeoff, not a free win.**
+   Since looking up sibling `y` from `x`'s visit already yields everything
+   needed to fill `y`'s own annex entry too, one visit per *family* could in
+   principle replace one visit per *observed family member* — cutting this
+   pass's cost roughly by the average family size instead of paying
+   `3 * N_distinct` regardless. But it means threads processing different
+   source k-mers can end up writing the *same* sibling's slot concurrently —
+   the fully independent, ownership-free parallelism of the plan above is
+   deliberately traded away for this gain. It stays safe only because the
+   computed value for a given slot is deterministic regardless of who
+   computes it, so redundant concurrent writes converge to the same
+   value — correct as long as each write is atomic, no locking needed — but
+   it is a real design complexity increase over "every member redoes its
+   own 3 lookups independently," not a strict improvement to adopt by
+   default.
+4. **Trigger and caching** (`obikindex::KmerIndex`/`distance.rs`): compute
+   lazily on first `distance` call for an SNP-family metric against a given
+   index; check for an existing annex file first (mirrors
+   `PersistentBitMatrix::open()`'s auto-detect-and-fall-back,
+   `bitmatrix.rs:264-287`); if absent, run step 3 and persist; if present,
+   mmap and reuse.
+5. **Invalidation.** `merge` and `filter` always write to a fresh `--output`
+   directory (`obikindex/src/merge.rs`, `obikmer/src/cmd/filter.rs`) so a
+   re-merged/re-filtered index simply has no annex yet — nothing to
+   invalidate. `select --in-place` (`select_layer.rs:139-235`) mutates
+   columns of an existing index without changing its location, which can
+   change sibling counts without removing rows — it must explicitly delete
+   any cached annex for that index as part of its in-place rewrite.
+6. **Testing**: hand-built tiny indexes with known sibling counts (0-3);
+   order-independence (recompute twice on a static index, identical
+   result, given the fully-parallel no-ownership design); invalidation
+   (annex absent/correctly recomputed after `select --in-place`); once
+   Step 2's sweep is fixed, a regression check that sibling_count == 0
+   slots are never looked up cross-partition during the sweep.
+
+Cost: `3 * N_distinct` existence-only lookups, computed once per index
+state and amortised over every subsequent `distance` call that reuses the
+cached annex — cheaper per-lookup than the sweep itself (hit/miss only, no
+column fetch).
+
+### Step 3 — finalisation (`obikindex`)
+
+From the global `SnpTally` alone (diagonal and off-diagonal both populated by
+the sweep, see Step 1/2 — no dependency on the external `shared_kmers`
+matrix), derive n x n distance matrices, each a pure function of the
+accumulated counts (same shape as `jaccard_to_mash`):
+
+- `p_hat[i,j] = SNP / (SNP + shared)`
+- Jukes-Cantor, Kimura-2P (from `P`, `Q`), optionally LogDet (needs the
+  diagonal + base-composition margins).
+
+Guard the singularities (`p >= 3/4` for JC, `1-2P-Q <= 0` or `1-2Q <= 0` for
+K2P) by clamping to a max distance, as `jaccard_to_mash` clamps `J <= 0`.
+
+### Step 4 — surfacing (`obikindex` + `obikmer` CLI)
+
+These metrics do not fit `DistanceMetric`'s current `LayeredStore`-partial
+dispatch (they need the cross-partition sweep and produce a different
+intermediate). Two options, to decide:
+
+- **(a)** New `DistanceMetric` variants (`Pdistance`, `JukesCantor`,
+  `Kimura2P`, `LogDet`) whose `KmerIndex::distance` arm calls the sweep
+  (`snp.rs`) instead of the partial path, still returning `DistanceOutput`.
+  Keeps one CLI surface (`--metric jukes-cantor`), at the cost of a branch in
+  `distance()` that ignores the `LayeredStore` it built.
+- **(b)** A dedicated pathway (`KmerIndex::snp_distance`) and a distinct CLI
+  entry, if mixing a cross-partition sweep into the partition-local `distance`
+  command is judged architecturally muddy.
+
+Recommendation: (a) for user ergonomics (all pairwise distances under
+`distance`, all feeding NJ/UPGMA/`--shared-kmers` unchanged), but compute the
+sweep lazily only when an SNP-family metric is requested, so the existing
+metrics keep their partition-local fast path untouched.
+
+### Step 5 — subsampling flag
+
+Add `--snp-sample <fraction>` (or a bottom-`s` hash threshold): restrict the
+source-k-mer enumeration in Step 2 to `seq_hash(kmer) < threshold`. Divides
+lookups proportionally; `p_hat` is unbiased. Off by default (exact).
+
+### Testing
+
+- **Primitive unit tests**: `central_canonical_neighbors` on hand-checked
+  k-mers incl. palindrome-boundary cases; lone-k-mer `minimizer` against
+  `RollingStat`'s incremental result on the same k-mer.
+- **End-to-end tiny index**: two 1-genome indexes differing by a handful of
+  known isolated SNPs (transitions and transversions placed by hand), assert
+  exact `SNP`, `P`, `Q` counts and the resulting JC/K2P values.
+- **Dedup invariant**: assert the tally is identical regardless of genome/
+  partition order and that no pair is double-counted (compare against a
+  brute-force all-pairs reference on a small index).
+- **Subsampling**: `p_hat` within sampling error of the exact run.
+
+### Suggested phasing
+
+1. Step 0 primitives + their unit tests (self-contained, no distance wiring).
+   This also unblocks the long-declared-but-unimplemented `query --mismatch`
+   (`obikmer/src/cmd/query.rs:676`, currently a warning), which needs the same
+   neighbour + routing machinery.
+2. `SnpTally` + finalisation math with a brute-force (non-swept) reference
+   backend, validated on a tiny index.
+3. The real per-partition sweep (Step 2) behind the same finalisation; assert
+   it matches the brute-force backend.
+4. CLI surfacing (Step 4a) and NJ/UPGMA integration (already generic over the
+   matrix).
+5. Subsampling (Step 5).
+
+## References
+
+The Mash mutation-rate model this discussion contrasts with:
+[@Mash-distances-doc; @Fan2015-mash-formula].

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.38"
+version = "1.1.41"
 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()
+    let n_cores = obisys::effective_parallelism();
-        .map(|n| n.get())
-        .unwrap_or(1);
     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()
+    let n_cores = obisys::effective_parallelism();
-        .map(|n| n.get())
-        .unwrap_or(1);
     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.38"
+version = "1.1.41"
 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()
+        default_value_t = obisys::effective_parallelism()
-            .map(|n| n.get())
-            .unwrap_or(1)
     )]
     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)]
+    /// (negative: offset from group size, e.g. -1 = all but one)
-    pub min_count: Option<usize>,
+    #[arg(long, allow_hyphen_values = true)]
+    pub min_count: Option<isize>,
 
     /// Maximum number of ingroup genomes containing the k-mer
-    #[arg(long)]
+    /// (negative: offset from group size, e.g. -1 = all but one)
-    pub max_count: Option<usize>,
+    #[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)]
+    /// (negative: offset from outgroup size, e.g. -1 = all but one)
-    pub min_outgroup_count: Option<usize>,
+    #[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)
+    /// (default 0 when --outgroup is set, no constraint otherwise;
-    #[arg(long)]
+    /// negative: offset from outgroup size, e.g. -1 = all but one)
-    pub max_outgroup_count: Option<usize>,
+    #[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 min_count:          Option<isize>,
-    pub max_count:          Option<usize>,
+    pub max_count:          Option<isize>,
     pub min_frac:           Option<f64>,
     pub max_frac:           Option<f64>,
-    pub min_outgroup_count: Option<usize>,
+    pub min_outgroup_count: Option<isize>,
-    pub max_outgroup_count: Option<usize>,
+    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);
+    // Resolve a signed count: negative means an offset from the group size
-    let max_count          = p.max_count.unwrap_or(in_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 min_outgroup_count = p.min_outgroup_count.map(|v| resolve(v, out_size)).unwrap_or(0);
-    let max_outgroup_count = p.max_outgroup_count.unwrap_or(default_max_outgroup_count);
+    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()
+        default_value_t = obisys::effective_parallelism()
-            .map(|n| n.get())
-            .unwrap_or(1)
     )]
     pub threads: usize,
 

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.
+                // No data matrix: implicit presence — all values are 1. `row`
-                // The filter result is identical for every kmer, so evaluate once.
+                // 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() {
-                    for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
+                    if mphf.find(kmer).is_some() && passes_all(filters, kmer, &all_present, n_genomes) {
-                        if mphf.find(kmer).is_some() {
+                        cont = cb(kmer, all_present.clone());
-                            cont = cb(kmer, all_present.clone());
+                        if !cont { break; }
-                            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() {
-                    for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
+                    if mphf.find(kmer).is_some() && passes_all(filters, kmer, &all_present, n_genomes) {
-                        if mphf.find(kmer).is_some() {
+                        cont = cb(part, layer, kmer, all_present.clone());
-                            cont = cb(part, layer, kmer, all_present.clone());
+                        if !cont { break; }
-                            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 {
+pub fn passes_all(
-    filters.iter().all(|f| f.passes(row, n_genomes))
+    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 {
+#[path = "tests/iter.rs"]
-    use super::*;
+mod tests;
-    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/obiskbuilder/src/lib.rs

@@ -10,8 +10,8 @@ pub mod stream_iter;
 mod scratch;
 
 pub(crate) mod encoding;
-pub(crate) mod entropy_table;
+#[allow(missing_docs)]
-pub(crate) mod rolling_stat;
+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.
+    // Minimizer selection state.
-    k1q: Ring<u64, 32>,
-    k2q: Ring<u64, 32>,
-    k3q: Ring<u64, 32>,
-    k4q: Ring<u64, 32>,
-    k5q: Ring<u64, 32>,
-    k6q: Ring<u64, 32>,
     minimier: Ring<MmerItem, 32>,
 
-    // Frequency count arrays.
+    // Entropy tracking, composed as a plain inline field so both concerns
-    // Max count per cell ≤ k ≤ 31 → u8 is sufficient.
+    // update in the same streaming pass without being conflated in one
-    k1c: [u8; 4],
+    // struct — see `obikentropy::EntropyTracker`.
-    k2c: [u8; 16],
+    entropy: EntropyTracker,
-    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],
 }
 
 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],
+            entropy: EntropyTracker::new(k),
-            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],
         }
     }
 
@@ -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.entropy.reset();
-        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);
     }
 
     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 {
+        self.entropy.push(self.received, self.rolling_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;
-                        }
-                    }
-                }
-            }
-        }
     }
 
     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)
+        Some(self.entropy.normalized_entropy(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 })
     }
 }

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()
+        let n_cores = effective_parallelism();
-            .map(|n| n.get())
-            .unwrap_or(1);
 
         // column widths
         let nw = self