Merge pull request 'fix: resolve test race conditions, add logging, and fix CI deadlock' (#66) from push-kywuzlnvrqyx into main

Reviewed-on: #66

.gitea/workflows/ci.yml

@@ -1,4 +1,4 @@
-name: CI
+pname: CI
 
 on:
   pull_request:
@@ -25,8 +25,8 @@ jobs:
             ~/.cargo/registry
             ~/.cargo/git
             src/target
-          key: ${{ runner.os }}-cargo-${{ hashFiles('src/Cargo.lock') }}
-          restore-keys: ${{ runner.os }}-cargo-
+          key: ${{ runner.os }}-cargo-v2-${{ hashFiles('src/Cargo.lock') }}
+          restore-keys: ${{ runner.os }}-cargo-v2-
 
       - name: Build
         run: cargo build --release

.gitignore

@@ -9,6 +9,7 @@ data-stress
 ./**/*.json
 *.bin
 *.log
+*.csv
 Betula_exilis--IGA-24-33
 benchmark/genomes
 benchmark/simulated_data
@@ -23,3 +24,5 @@ benchmark/reference_dist
 benchmark/obikmer_dist
 benchmark/specific_index_count
 benchmark/specific_index_presence
+TNT
+phyg

docmd/implementation/filtering.md

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

docmd/implementation/obicompactvec.md

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

docmd/index.md

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

docmd/references.bib

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

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

@@ -1701,17 +1701,22 @@ dependencies = [
  "obikpartitionner",
  "obikseq",
  "obilayeredmap",
+ "obipipeline",
+ "obiread",
+ "obiskbuilder",
  "obiskio",
  "obisys",
  "rayon",
  "serde",
  "serde_json",
+ "tempfile",
  "tracing",
+ "tracing-subscriber",
 ]
 
 [[package]]
 name = "obikmer"
-version = "1.1.39"
+version = "1.1.44"
 dependencies = [
  "clap",
  "csv",

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/obidebruinj/src/tests/debruijn.rs

@@ -1,6 +1,17 @@
 use super::*;
 use obikseq::{k, set_k, unitig::Unitig, Kmer};
 
+// `obikseq::params` is process-wide (see obikseq/src/params.rs): tests in this
+// file don't all use the same `k` (`push_palindrome_single_node` needs an
+// even k=4 — no odd-length self-revcomp palindrome exists — while the rest
+// use k=5), and they run concurrently by default. Serialize the
+// set_k-through-use critical section across this file's tests so one test's
+// `k` can never be overwritten mid-flight by another.
+static K_LOCK: std::sync::Mutex<()> = std::sync::Mutex::new(());
+fn lock_k() -> std::sync::MutexGuard<'static, ()> {
+    K_LOCK.lock().unwrap_or_else(|e| e.into_inner())
+}
+
 // Build a graph from an ASCII sequence, inserting all canonical k-mers.
 fn graph_from_ascii(seq: &[u8]) -> GraphDeBruijn {
     let mut g = GraphDeBruijn::new();
@@ -37,6 +48,7 @@ fn collect_unitigs(g: &GraphDeBruijn) -> Vec<Unitig> {
 #[test]
 fn push_deduplicates_revcomp() {
     let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let kmer = Kmer::from_ascii(b"ACGTA").unwrap();
     let mut g = GraphDeBruijn::new();
@@ -49,6 +61,7 @@ fn push_deduplicates_revcomp() {
 fn push_palindrome_single_node() {
     // ACGT is its own revcomp
     let k = 4;
+    let _guard = lock_k();
     set_k(k);
     let kmer = Kmer::from_ascii(b"ACGT").unwrap();
     assert_eq!(kmer, kmer.revcomp(), "test requires a palindrome");
@@ -71,6 +84,7 @@ fn linear_chain_graph() -> (GraphDeBruijn, Vec<CanonicalKmer>) {
 #[test]
 fn degrees_linear_chain_node_count() {
     let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let (g, kmers) = linear_chain_graph();
     assert_eq!(g.len(), kmers.len());
@@ -82,6 +96,7 @@ fn degrees_linear_chain_extensions() {
     // Note: start_iter must not be consumed standalone — its second pass only
     // finds true cycle nodes when interleaved with chain traversal (iter_unitig).
     let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let seq = b"AAAAGGGG";
     let g = graph_from_ascii(seq);
@@ -118,6 +133,7 @@ fn kmers_from_unitigs(unitigs: &[Unitig]) -> Vec<CanonicalKmer> {
 fn unitig_roundtrip_linear() {
     // Non-repetitive sequence: all k-mers must be recovered across unitigs.
     let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let seq = b"ACCTGGCTA";
     let g = graph_from_ascii(seq);
@@ -136,6 +152,7 @@ fn unitig_roundtrip_longer_sequence() {
     // Longer non-repetitive sequence with no repeated k-mer of length k.
     // ACGTGGCTATCGAC with k=5 → 10 distinct k-mers, one linear chain.
     let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let seq = b"ACGTGGCTATCGAC";
     let g = graph_from_ascii(seq);
@@ -152,6 +169,7 @@ fn unitig_roundtrip_longer_sequence() {
 fn unitig_isolated_node() {
     // Single k-mer with no neighbours
     let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let kmer = Kmer::from_ascii(b"ACGTA").unwrap();
     let mut g = GraphDeBruijn::new();
@@ -165,6 +183,7 @@ fn unitig_isolated_node() {
 #[test]
 fn unitig_two_isolated_nodes() {
     let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let mut g = GraphDeBruijn::new();
     // Two k-mers that share no (k-1)-overlap
@@ -177,6 +196,7 @@ fn unitig_two_isolated_nodes() {
 #[test]
 fn unitig_two_truly_distinct_isolated_nodes() {
     let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let mut g = GraphDeBruijn::new();
     g.push(Kmer::from_ascii(b"AAAAC").unwrap().canonical());
@@ -192,7 +212,8 @@ fn unitig_two_truly_distinct_isolated_nodes() {
 
 #[test]
 fn no_kmer_lost_or_duplicated() {
-    let k = 7;
+    let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let seq = b"ACGTACGTACGTTTTTACGTACGT";
     let g = graph_from_ascii(seq);
@@ -218,6 +239,7 @@ fn cycle_kmers_not_lost() {
     // start_iter first pass yields nothing (all nodes internal); second pass
     // picks up cycle entries. All 4 k-mers must appear in the unitigs.
     let k = 5;
+    let _guard = lock_k();
     set_k(k);
     let seq = b"ACGTACGT";
     let g = graph_from_ascii(seq);
@@ -240,6 +262,7 @@ fn branching_graph_no_kmer_lost_or_duplicated() {
     // Each "node" is a distinct 5-mer; edges share a 4-mer suffix/prefix.
     // We use long non-repetitive sequences and extract only the required kmers.
     let k: usize = 5;
+    let _guard = lock_k();
     set_k(k);
     let mut g = GraphDeBruijn::new();
 

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,11 @@ 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"
+tracing-subscriber = { version = "0.3", features = ["fmt", "env-filter"] }
+
 [features]
 default = ["numa"]
 numa    = ["hwlocality"]

src/obikindex/src/distance.rs

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

src/obikindex/src/lib.rs

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

src/obikindex/src/numa.rs

@@ -79,9 +79,7 @@ pub fn build() -> NumaSetup {
     }
 
     // UMA fallback: single synthetic node, all cores, no pool, no pinning.
-    let n_cores = std::thread::available_parallelism()
-        .map(|n| n.get())
-        .unwrap_or(1);
+    let n_cores = obisys::effective_parallelism();
     debug!("UMA: single synthetic node, {} core(s)", n_cores);
     NumaSetup {
         pools: vec![None],
@@ -91,9 +89,7 @@ pub fn build() -> NumaSetup {
 
 #[cfg(not(feature = "numa"))]
 pub fn build() -> NumaSetup {
-    let n_cores = std::thread::available_parallelism()
-        .map(|n| n.get())
-        .unwrap_or(1);
+    let n_cores = obisys::effective_parallelism();
     debug!("UMA: single synthetic node, {} core(s)", n_cores);
     NumaSetup {
         pools: vec![None],
@@ -299,20 +295,27 @@ impl PartitionRunner {
                     let pool = node.pool.clone();
 
                     s.spawn(move || {
+                        let tid = std::thread::current().id();
+                        debug!(?tid, "PartitionRunner worker: waiting on activation");
                         if arx.recv().is_err() {
+                            debug!(?tid, "PartitionRunner worker: activation channel closed, exiting");
                             return;
                         }
+                        debug!(?tid, "PartitionRunner worker: activated");
                         if !cpu_ids.is_empty() {
                             pin_current_thread(cpu_ids);
                         }
                         for i in &prx {
+                            debug!(?tid, partition = i, "PartitionRunner worker: picked partition");
                             let t = Instant::now();
                             let r = match &pool {
                                 Some(p) => p.install(|| f(i)),
                                 None => f(i),
                             };
+                            debug!(?tid, partition = i, "PartitionRunner worker: partition done");
                             etx.send(WorkerEvent::Completed(i, r, t.elapsed())).ok();
                         }
+                        debug!(?tid, "PartitionRunner worker: no more partitions, exiting");
                     });
                 }
             }
@@ -323,13 +326,18 @@ impl PartitionRunner {
             // ── Controller ────────────────────────────────────────────────────
             let mut activation = NodeActivation::new(&activate_txs, &node_caps, max_workers);
             activation.activate_initial(INITIAL_DIVISOR, n_total);
+            debug!(n_total, activated = activation.total(), "PartitionRunner controller: initial activation");
 
             let mut cpu_sample = CpuSample::now();
             let mut io_sample = IoSample::now();
             let mut completed = 0usize;
 
             while completed < n_total {
-                let Ok(event) = event_rx.recv() else { break };
+                debug!(completed, n_total, "PartitionRunner controller: waiting for an event");
+                let Ok(event) = event_rx.recv() else {
+                    debug!("PartitionRunner controller: event channel closed, stopping");
+                    break;
+                };
                 match event {
                     WorkerEvent::Completed(i, r, dur) => {
                         match r {

src/obikmer/Cargo.toml

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

src/obikmer/src/cli.rs

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

src/obikmer/src/cmd/distance.rs

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

src/obikmer/src/cmd/predicate.rs

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

src/obikmer/src/cmd/query.rs

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

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/params.rs

@@ -7,12 +7,28 @@
 //! different value panics.  This prevents silent divergence between the global
 //! parameter and the values used to build data structures.
 //!
-//! In test builds (`#[cfg(test)]`) the same public API is backed by
-//! `thread_local!` [`Cell`]s instead.  Each test thread gets its own
-//! independent copies of `K` and `M`, so tests can use arbitrary values
-//! without coordinating with one another and without any reset mechanism.
-//! The `OnceLock` constraint is deliberately absent: test isolation is
-//! provided by thread locality, not by write-once semantics.
+//! In test builds (`#[cfg(test)]`) the same public API is backed by plain
+//! process-wide atomics instead, freely overwritable (no write-once
+//! constraint) so tests don't need a reset mechanism between runs.
+//!
+//! An earlier version of this module used `thread_local!` `Cell`s here,
+//! reasoning that "each test thread gets its own copy" gives isolation
+//! between tests using different `k`/`m` values. That assumption broke as
+//! soon as any code under test fanned work out to *other* threads it
+//! doesn't control — `PartitionRunner`'s pre-spawned workers, or a bare
+//! `rayon::par_iter()` — since a freshly spawned thread never inherits the
+//! calling test thread's thread-local state, silently reading back `k=0`
+//! there instead (surfaced as a `bitvec`/slice-indexing panic deep inside
+//! whatever used the bogus length). Process-wide atomics make `k()`/`m()`
+//! correct on *any* thread without every call site having to know to
+//! re-propagate them. Unset still silently reads back as `0` (same as the
+//! old `Cell` default) rather than panicking: several existing tests read
+//! `m()` without ever calling `set_m` themselves, relying on that default.
+//! The trade-off: tests that genuinely need different `k`/`m` values from
+//! other tests must not run concurrently with them in the same process (in
+//! practice: every test file in this workspace already uses one fixed
+//! `k`/`m` pair for all its own tests, so this doesn't currently cost
+//! anything).
 
 // ── Production implementation ─────────────────────────────────────────────────
 
@@ -44,22 +60,22 @@ mod state {
 
 // ── Test implementation ───────────────────────────────────────────────────────
 //
-// Each test thread owns its private K and M via thread_local!, so tests may
-// call set_k / set_m with any value without affecting other tests.
+// Process-wide, freely overwritable (no write-once constraint), visible from
+// any thread — including threads a test doesn't spawn itself (rayon workers,
+// PartitionRunner workers, ...). `0` (never explicitly set) is returned as-is,
+// same default as the old thread-local `Cell`.
 
 #[cfg(any(test, feature = "test-utils"))]
 mod state {
-    use std::cell::Cell;
+    use std::sync::atomic::{AtomicUsize, Ordering};
 
-    thread_local! {
-        static K: Cell<usize> = Cell::new(0);
-        static M: Cell<usize> = Cell::new(0);
-    }
+    static K: AtomicUsize = AtomicUsize::new(0);
+    static M: AtomicUsize = AtomicUsize::new(0);
 
-    pub fn set_k(k: usize) { K.with(|c| c.set(k)); }
-    pub fn k()    -> usize { K.with(|c| c.get()) }
-    pub fn set_m(m: usize) { M.with(|c| c.set(m)); }
-    pub fn m()    -> usize { M.with(|c| c.get()) }
+    pub fn set_k(k: usize) { K.store(k, Ordering::SeqCst); }
+    pub fn k() -> usize { K.load(Ordering::SeqCst) }
+    pub fn set_m(m: usize) { M.store(m, Ordering::SeqCst); }
+    pub fn m() -> usize { M.load(Ordering::SeqCst) }
 }
 
 // ── Public API (identical signature in both configurations) ───────────────────

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/src/lib.rs

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

src/obiskbuilder/src/rolling_stat.rs

@@ -196,13 +196,6 @@ 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;
-        }
-        Some(self.entropy.entropy(order))
-    }
-
     pub fn normalized_entropy(&self) -> Option<f64> {
         if !self.ready() {
             return None;

src/obiskio/src/tests/reader.rs

@@ -102,9 +102,9 @@ fn roundtrip_single() {
 
 #[test]
 fn roundtrip_all_lengths() {
-    obikseq::params::set_k(11);
+    setup();
     let bases: Vec<u8> = (0..300).map(|i| b"ACGT"[i % 4]).collect();
-    for len in (11..=19).chain([255, 256, 257]) {
+    for len in (TEST_K..=19).chain([255, 256, 257]) {
         let sk = make_sk(&bases[..len]);
         let mut buf = Vec::new();
         sk.write_to_binary(&mut buf).unwrap();

src/obisys/src/lib.rs

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