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

@@ -24,3 +24,5 @@ benchmark/reference_dist
 benchmark/obikmer_dist
 benchmark/specific_index_count
 benchmark/specific_index_presence
+TNT
+phyg

docmd/theory/evolutionary_distances.md

@@ -203,6 +203,349 @@ 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.
 
+## Context, detectability, and a 3-way ordinal distance per pair
+
+Empirical follow-up to the pseudo-alignment idea above: `obikmer distance
+--snp` was run on a real 20-genome benchmark index and the resulting FASTA
+fed to `raxml-ng`. Two problems surfaced, both traced back to conflating
+distinct notions under one symbol.
+
+**"Context", precisely.** Sharing a central base between two genomes is not
+just sharing a nucleotide — it is sharing a **context**: the `2m` flanking
+bases, identical, which is a homology claim about that flanked window
+(guaranteed non-coincidental by k-specificity, Bias 4 above), *not* a claim
+about orthology or paralogy of the copy each genome carries. `A` opposite `C`
+= same context, divergent centre. `A` opposite nothing = **this context is
+not observed in one of the two genomes** — informative, not neutral.
+
+**Why the IUPAC/DNA encoding used for the first `--snp` test was wrong.**
+Feeding IUPAC-coded ambiguity into a standard DNA model (`raxml-ng --model
+GTR+G`) is a semantic mismatch: Felsenstein-pruning ML treats an ambiguous
+tip as "exactly one true state, unknown which" (a uniform partial-likelihood
+vector over compatible bases), not "these states are simultaneously
+present". The two encodings look identical (same IUPAC letters) but the
+software reads them backwards from what was intended — this invalidates the
+literal branch lengths from that first experiment (topology-level groupings
+by genus were still informative, see the worked example further down).
+
+**Why `-` (absence) must not be scored as similarity, but also must not be
+scored as a shared character between two absences.** Two genomes both
+lacking a context are not observed to resemble each other at that locus —
+neither is observed to differ from the other either. It is a symmetric
+non-observation, uninformative for that pair, and should contribute nothing
+(not a small positive nor a small negative signal) to their distance. A
+genome carrying a state (`A`) against one carrying none is a different case
+entirely: informative, and should not be scored as neutral "missing data"
+the way a generic DNA/ML pipeline would.
+
+**Detectability vs existence — a deliberate simplification, accepted.**
+"Context not observed" conflates two different biological events: (1) true
+loss of the locus, (2) the locus still exists but a mutation/indel *outside*
+the centre, anywhere in the `2m` flanks, broke k-mer recognition. The design
+adopts a rigorist stance on purpose: any flank-breaking mutation counts as
+"this context no longer exists", full stop — because both causes (1) and (2)
+independently require *at least* one more mutational event than a lone
+central substitution would. This licenses treating "context absent in one of
+the two genomes" as a **lower bound** on distance strictly greater than a
+plain central SNP, without needing to know which of the two causes applies.
+Coarser than a true event count, and accepted as such (fine substitution-type
+modelling, e.g. transition/transversion weighting, is a secondary
+refinement, not required for this to be useful).
+
+**Resulting ordinal distance between two genomes at one family/context:**
+
+| Comparison | Distance | Meaning |
+|---|---|---|
+| same centre (`A`/`A`) | `0` | identical |
+| different centre, both single-copy (`A`/`C`) | `1` | plain central SNP |
+| one genome has a state, the other has none | `>1` (lower bound) | context undetectable in one genome — at least one extra mutational event, of unknown type |
+| neither genome has any state (`∅`/`∅`) | excluded | symmetric non-observation, not comparable, contributes nothing |
+
+This is a direct extension of `KmerIndex::raw_snp_distance` (`obikindex/src/siblings.rs`),
+which today only implements the `0`/`1` rows and silently drops everything
+else (including the informative `>1` row) rather than scoring it.
+
+**Open, not yet resolved:**
+- Calibrating `>1` to a real number for tools expecting continuous distances
+  (NJ/UPGMA, ML branch lengths), rather than an arbitrary placeholder.
+  Natural route: estimate `p_hat` from the resolved (`0`/`1`) sites first,
+  then use the already-derived ascertainment formula (`P(usable window
+  showing a central SNP) = p * (1-p)^(2m)`, Bias 1 above) to derive a
+  model-consistent value for the `>1` bucket instead of guessing a constant.
+- Where multi-copy/ambiguous states (the IUPAC case: a genome carrying more
+  than one form) fit into this ordinal scheme — plausibly also `>1` by the
+  same "at least one extra event" argument (a second form appearing is a
+  gain, itself an event), but not yet worked out.
+
+**Practical alternative validated for the pseudo-alignment output itself**
+(orthogonal to the ordinal-distance question above, useful regardless of how
+`>1` ends up calibrated): re-encode each family as 4 independent binary
+presence/absence characters (`A`,`C`,`G`,`T` columns) instead of one IUPAC
+column, feed to a `BIN`-type model instead of `DNA`. `∅` becomes an explicit
+`0000` state (identity with another `0000`, not missing data) rather than a
+gap — removes the semantic mismatch above by construction. Known cost,
+accepted for now: a plain central substitution (`A` -> `C`) becomes 2 binary
+flips (`1000` -> `0100`), overweighting substitutions relative to true
+gain/loss events, and the 4 sub-characters of one family are not
+statistically independent the way a generic `BIN` model assumes. A proper
+fix (single 16-state alphabet, i.e. the powerset of `{A,C,G,T}`, with a
+substitution-rate structure that respects the subset lattice rather than a
+fully general 16x16 GTR-analogue) is very likely not expressible in
+`raxml-ng`'s `MULTI` datatype as-is (Mk or fully-general rates only) and a
+fully general 16-state rate matrix is almost certainly unidentifiable here
+(states of cardinality >=3 are ~2% of sites in the benchmark run). Treated as
+a longer-term research question, not a near-term implementation target.
+
+## Sankoff parsimony as the resolution of the 16-state model problem
+
+The "longer-term research question" just above (a 16-state alphabet — the
+powerset of `{A,C,G,T}` — with a substitution structure that respects the
+subset lattice) turns out to have a near-term answer, once the *unification*
+question below is worked through.
+
+**Distance methods (NJ/UPGMA/ME) vs. character methods (parsimony/ML): not a
+deep philosophical divide, but a real practical distinction for this
+project.** Historically "phenetic" (characters -> distances -> tree) and
+"cladistic" (characters -> tree directly) approaches were presented as
+opposed schools; the modern view is a mathematical continuity, not a
+dichotomy — Minimum Evolution (ME: find the tree minimizing total branch
+length from a distance matrix) and Maximum Parsimony (MP: find the tree
+minimizing total character-state changes) are both instances of "minimise a
+global explanatory cost", and coincide under simple encodings (see Farris
+1983, "The Logical Basis of Phylogenetic Analysis"; the MP/ME connection is
+developed in the Minimum Evolution / Balanced Minimum Evolution literature,
+e.g. Nei and colleagues — citations not independently re-verified here, flag
+before quoting further). NJ's own agglomeration step already uses the whole
+distance matrix jointly (the Q-matrix), not just the pair being merged — an
+earlier claim in this discussion that distance methods are "blind" to
+cross-taxon structure at every stage was too strong.
+
+What *does* remain a real, structural distinction for this project: in a
+character method, a given character's cost is **re-evaluated per candidate
+topology** during tree search (the same family can cost 1 change under one
+topology, 2 under another). In a pairwise-distance pipeline
+(`raw_snp_distance` as it exists today), each family's contribution to
+`d(i,j)` is computed **once**, independent of any candidate topology, before
+NJ/UPGMA ever runs — so a question like "does this shared `∅` look like a
+synapomorphy under topology T" can never be posed in that pipeline, for any
+T. That question is only answerable by a method that tests candidate
+topologies and re-scores characters under each — i.e. a character method.
+
+**Sankoff parsimony directly resolves the `∅`/gain-loss/substitution
+question, without the identifiability problem of a fitted 16-state model.**
+Sankoff's algorithm generalises Fitch parsimony to an arbitrary
+user-supplied cost matrix between states (`obikseq`/`obikindex` would treat
+each family as a `2^{4}`-state character, state = subset of `{A,C,G,T}`
+observed in that genome, `∅` included as a real state, not a gap). The
+previous 16-state idea failed specifically because *fitting* a full 16x16
+rate matrix by ML is unidentifiable at this data volume; Sankoff sidesteps
+that because the cost matrix is **fixed a priori from domain knowledge**, not
+estimated — e.g. `c({A},{C}) = 1` (a substitution), `c({A},{A,C}) = 1` (a
+gain), `c({A,C},{A}) = 1` (a loss), `c({A,C},{G,T}) = 2` (two changes) — no
+estimation, no overparameterisation. This reframes gain/loss and central
+substitution as two cost categories with independently chosen weights,
+exactly the "two families of parameters" (`mu_substitution`, `mu_gain/loss`)
+floated earlier in this discussion, now with an actual algorithmic home.
+
+**Caveat, not blocking for this project's scope.** Sankoff is still
+parsimony: in principle exposed to Felsenstein's statistical-inconsistency
+result under long-branch attraction (already invoked earlier against
+`D_F = min(a,b)`) — parsimony and ML only provably coincide in the
+short-branch regime. This is not a practical concern here because it is
+exactly this estimator's declared target (closely related genomes, short
+branches) — the regime where parsimony's known failure mode does not apply —
+but worth stating explicitly as a scope guard rather than leaving it
+implicit.
+
+**Cheapest next experiment: don't write a Sankoff tree-search engine, use
+one that exists.** The hard part of a from-scratch implementation is not the
+Sankoff DP itself (a straightforward dynamic program over a *fixed* tree) but
+the topology search (SPR/NNI with incremental re-scoring) that comes for
+free with `raxml-ng` on the ML side. **TNT** (Tree analysis using New
+Technology, free, standard in morphological cladistics) already implements
+Sankoff parsimony with a custom cost matrix plus topology search — the
+family-state matrix (already close to what `--snp` produces, minus the
+IUPAC/DNA-model mismatch) could be fed there directly, no new code required,
+before considering a bespoke engine.
+
+**The concrete comparison this unlocks:** run both pipelines on the same
+family data —
+`k-mer families -> D_ij -> NJ/ME/UPGMA` (phenetic, what exists today) vs.
+`k-mer families -> characters -> argmin_T Sankoff-cost(T)` (cladistic, via
+TNT) — and compare the resulting topologies. Agreement would validate that
+the pairwise-distance projection preserves the phylogenetic signal;
+disagreement would pinpoint exactly what the projection to a single number
+per pair loses. Not yet run.
+
+### A concrete Sankoff cost matrix for the 16-state alphabet
+
+**Set-edit-distance formula.** For two states `X, Y ⊆ {A,C,G,T}`, split into
+`seulement_X = X \ Y` (size `a`) and `seulement_Y = Y \ X` (size `b`).
+Elements present in both cost nothing. Pair up to `min(a,b)` of the
+remaining elements as **substitutions** (cheaper than treating them as an
+unrelated loss + gain whenever `c_sub < 2*c_gl`, which any sane parameter
+choice satisfies); whatever is left over after pairing is a pure
+**gain/loss**:
+
+```
+cost(X,Y) = min(a,b)*c_sub + |a-b|*c_gl
+```
+
+Worked examples (`c_sub = c_gl = 1`): `{A}->{C}` = 1 (one substitution);
+`{A}->{A,C}` = 1 (one gain, no substitution pair available since nothing is
+only-in-Y that matches an only-in-X element after the shared `A` is
+excluded); `{A,C}->{G,T}` = 2 (two substitution pairs, `A/C` vs `G/T`, both
+same-size sets share nothing); `{A,C,G}->{A,C,T}` = 1 (`G`/`T` is the only
+mismatched pair, `A,C` shared).
+
+**`∅` is a flat-cost special case, not a instance of the general formula.**
+Applying the formula naively to e.g. `{A,C,G} -> ∅` would charge `3*c_gl`
+(three independent losses). Reject that: total context disappearance
+(flank-breaking, or true structural loss) is plausibly **one** event, not
+`|X|` of them, so:
+
+```
+cost(X, ∅) = cost(∅, X) = c_ctx   (constant, independent of |X|)
+cost(∅, ∅) = 0
+```
+
+`c_ctx` should not be guessed — it is exactly the value already derived for
+the `>1` bucket in "Context, detectability, and a 3-way ordinal distance per
+pair" above (`(2m*p_hat)/(1-(1-p_hat)^(2m)) + p_hat`), reused rather than
+invented.
+
+**Better construction method: shortest path in a small state graph, not the
+closed-form formula directly.** Build a graph on the 16 states with two edge
+types — substitution edges between same-cardinality sets differing by one
+element (weight `c_sub`, or `c_ts`/`c_tv` if split further below), and
+gain/loss edges between sets whose cardinality differs by one (weight
+`c_gl`) — then define `cost(X,Y)` as shortest-path distance in that graph,
+precomputed once (16 nodes, trivial) into a dense 16x16 matrix before
+feeding it to Sankoff/TNT. Verified equivalent to the closed-form formula
+above in the uniform-cost case (checked by hand on `{A}->{C,G,T}`: both give
+`c_sub + 2*c_gl`). The graph construction is not just a reformulation for
+its own sake: it is the version that generalises correctly once substitution
+costs stop being uniform (next point) — the closed-form's `min(a,b)`
+counting silently assumes *any* pairing costs the same, which breaks the
+moment `c_sub` depends on which two bases are involved.
+
+**Transition/transversion refinement.** Split `c_sub` into `c_ts` (A<->G or
+C<->T) and `c_tv` (the other four pairs) — already well-defined in canonical
+space (see "Canonical invariance" above: a transition maps to a transition,
+a transversion to a transversion, regardless of orientation). This turns
+"pick `min(a,b)` substitution pairs" into a genuine (tiny, <=4 elements per
+side, trivially enumerable) minimum-cost bipartite matching problem instead
+of a plain count — the state-graph shortest-path construction handles this
+automatically, no separate logic needed. Calibration is not a guess either:
+the "Sufficient statistic: 4x4 base-pair tally" section below already plans
+to collect the joint `(centre_i, centre_j)` distribution over resolved (`0`
+or `1`) sites — that tally directly gives the empirical Ts/Tv ratio, from
+which `c_ts`/`c_tv` follow (e.g. `cost ∝ -log(observed rate)`, the standard
+generalised-parsimony step-weighting heuristic), reusing a statistic already
+planned rather than adding a new one.
+
+**`gamma`/`mu` (i.e. `c_gl`/`c_sub`) sensitivity sweep before calibration.**
+No strong prior on whether gain/loss events should cost more or less than a
+point substitution — duplication/deletion rates are not a priori equal to
+point-mutation rates, but the direction isn't obvious, and setting it too
+high risks the parsimony search *eliminating* exactly the
+heterozygosity/paralogy signal the design is meant to tolerate (see
+"Heterozygosity, ploidy..." below). Cheap first step: run the topology at a
+handful of ratios (`0.5, 1, 2, 5`) and check whether it's stable — a robust
+topology across that range is far more trustworthy than one built on a
+single, unvalidated guess. Empirical calibration of `c_gl` from the
+family-size distribution already available (`sibling_annex_stats`) is a
+natural follow-up once the sensitivity sweep shows the topology is worth
+refining further.
+
+**Caveat carried over from the `D_F = min(a,b)` rejection earlier:** this
+cost matrix is only valid as **Sankoff step-cost input**, re-evaluated for
+every branch of every candidate topology during tree search. Reusing
+`cost(leaf_A, leaf_B)` directly as a standalone pairwise distance (bypassing
+the tree) would reintroduce the exact circularity already rejected — the
+`min(a,b)` pairing here is a locally-defined edit distance between two
+states, not a claim about the true evolutionary history between two
+specific genomes.
+
+**Feasibility confirmed:** TNT's `costs` command accepts custom step
+matrices for multistate characters, so this whole construction (16x16
+matrix derived from the state graph, `c_ts`/`c_tv`/`c_gl`/`c_ctx` as
+tunable parameters) is directly usable there — no new tooling required
+before testing it.
+
+### Experiment: TNT run on real data (2026-08-11)
+
+**Status:** validated at genus/family/order scale within Bacteria; not
+informative across domains with this character type. Exploratory only —
+run entirely outside the repo (`/tmp/tnt_run`, TNT installed locally under
+`TNT/`), no new Rust code. Kept here as the record of what was learned.
+
+**Pipeline.** `obikmer distance --snp` emits one IUPAC-coded pseudo-alignment
+row per genome (`snp_pseudo_alignment`, one column per family with
+`family_size() >= 2`). A small external Python script decodes IUPAC back to
+the 16-state bitmask, applies the set-edit-distance formula above, and
+emits a complete TNT script (`xread` matrix + `smatrix` step-matrix + `hold`/
+`mult` search). No new Rust code was needed for this pass.
+
+**Calibration ("quick option").** Rather than modifying
+`write_raw_snp_distance_csv` to emit raw counts, `c_ctx` was approximated as
+the unweighted mean of the pairwise `snp/(snp+shared)` ratios already
+present in `rawsnp.csv`, restricted to the relevant taxon subset. Used
+values: `mu = gamma = 10`, `c_ctx = 18` (ratio `c_ctx/mu = 1.8`, consistent
+with the observed pairwise ratios for genomes within Enterobacteriaceae/
+Eubacteria, which cluster around 1.5-2 and barely move as the taxon set
+widens — the context signal is stable, not sensitive to which subset is
+chosen). The more principled route (raw counts, weighted `p_hat`, Ts/Tv
+split) remains a follow-up, not yet done.
+
+**Run 1 — Enterobacteriaceae (11 taxa: 4 *E. coli*, 4 *Salmonella enterica*,
+3 *Klebsiella pneumoniae*).** All three genera recovered as monophyletic.
+Initial read of the exported (unrooted, TNT/Nexus `[&U]`) tree as showing a
+genus-arrangement disagreement with known systematics (*Escherichieae*:
+*Escherichia*+*Salmonella* sister vs. more distant *Klebsielleae*) was
+**wrong** — diagnosed via `force = (taxa);` monophyly-constraint test
+(identical score constrained vs. free ⟹ no real disagreement). Root cause of
+the misreading: with exactly 3 clades and no outgroup, an unrooted tree has
+only **one possible topology** (a single trifurcation) — there is no
+internal arrangement to get right or wrong. This run cannot test
+inter-genus relationships at all; it can only test intra-genus monophyly
+(which held).
+
+**Run 2 — Eubacteria (18 taxa: run 1 + *Acidobacterium capsulatum*,
+*Opitutus terrae*, *Bacillus subtilis*, *Shouchella clausii*, *Wolbachia*
+endosymbiont, *Proteus mirabilis*, *Yersinia ruckeri*).** The
+Enterobacteriaceae substructure from run 1 is reproduced identically, now
+correctly rooted by real outgroups, resolving the tribal arrangement left
+undetermined in run 1: *Escherichia*+*Salmonella* sister, *Klebsiella* more
+distant — matching known systematics. Outgroup placement: `(Bacillus,
+Shouchella)` sister pair (Firmicutes/*Bacillales*) splits from all
+Proteobacteria — a correct phylum-level split; `Wolbachia`
+(Alphaproteobacteria) splits from the Gammaproteobacteria block
+(`Proteus`, `Yersinia`, Enterobacteriaceae) — a correct class-level split.
+Lower confidence: the fine nested order `(Proteus, (Yersinia,
+Enterobacteriaceae))` and the relative position of `Acidobacterium` vs.
+`Opitutus` — plausible, not independently verified against current
+Enterobacterales family-level literature.
+
+**Run 3 — full domain set (20 taxa: run 2 + *Candidozyma auris* [yeast] and
+*Saccharolobus islandicus* [archaeon]).** The bacterial clade from run 2 is
+reproduced **unchanged and intact** — a real robustness signal, the method
+does not fragment the ingroup when unrelated deep taxa are added. But the
+result carries **no information on Bacteria/Archaea/Eukarya relationships**:
+with exactly one eukaryote, one archaeon and one bacterial clade, the
+unrooted tree is again forced into the single 3-clade trifurcation from run
+1's caveat — there is no second representative of either outgroup domain to
+resolve internal arrangement, so nothing about their relative position can
+be read from the topology (a "ladder" ordering in the exported tree is
+serialization, not signal). Independently, `rawsnp.csv` shows *why* this
+character type cannot reach further: pairwise ratios involving
+*Candidozyma*/*Saccharolobus* are almost all `NA` (no central-position
+family shared at all) or saturated at `1.0` (every shared family differs) —
+central-position families require literal 31 bp context conservation, which
+simply does not survive domain-level divergence. Cross-domain placement
+would need conserved-marker characters (rRNA, ribosomal proteins), not this
+estimator.
+
 ## Heterozygosity, ploidy, and consensus-assembly inputs
 
 A within-genome multiplicity signal (more than one of the 4 central forms

src/Cargo.lock

@@ -1711,11 +1711,12 @@ dependencies = [
  "serde_json",
  "tempfile",
  "tracing",
+ "tracing-subscriber",
 ]
 
 [[package]]
 name = "obikmer"
-version = "1.1.40"
+version = "1.1.44"
 dependencies = [
  "clap",
  "csv",

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

@@ -24,6 +24,7 @@ hwlocality       = { version = "1.0.0-alpha.11", features = ["vendored"], option
 [dev-dependencies]
 obiread  = { path = "../obiread" }
 tempfile = "3"
+tracing-subscriber = { version = "0.3", features = ["fmt", "env-filter"] }
 
 [features]
 default = ["numa"]

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

@@ -350,7 +350,7 @@ impl KmerIndex {
         // every single variant — fine at the scale of a handful of test
         // k-mers, but with billions of lookups against a real index this
         // manifested as ~90% system time, observed in practice. ───────────
-        let n_workers = std::thread::available_parallelism().map(|n| n.get()).unwrap_or(4);
+        let n_workers = obisys::effective_parallelism();
         let capacity = 256;
 
         // Throttling is not optional once a `Flat` stage is in the pipeline
@@ -976,7 +976,23 @@ mod tests {
     /// the same primitives `obikmer`'s `scatter` step uses (minus the
     /// multi-file `obipipeline` wrapper — a single sequence needs none of
     /// that): normalise -> build superkmers -> route -> write.
+    /// `cargo test` doesn't install a `tracing` subscriber the way `obikmer`'s
+    /// CLI does, so `debug!`/etc. are silent no-ops by default — including the
+    /// `PartitionRunner` instrumentation that would matter most for
+    /// re-diagnosing a hang here. `try_init` is idempotent across concurrently
+    /// running tests (later calls just find a subscriber already installed).
+    fn init_tracing() {
+        let _ = tracing_subscriber::fmt()
+            .with_env_filter(
+                tracing_subscriber::EnvFilter::try_from_default_env()
+                    .unwrap_or_else(|_| tracing_subscriber::EnvFilter::new("info")),
+            )
+            .with_writer(std::io::stderr)
+            .try_init();
+    }
+
     fn build_single_genome_index(dir: &Path, label: &str, seq: &[u8]) -> KmerIndex {
+        init_tracing();
         let fasta_path = dir.join(format!("{label}.fasta"));
         let mut f = std::fs::File::create(&fasta_path).unwrap();
         writeln!(f, ">{label}").unwrap();

src/obikmer/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "obikmer"
-version = "1.1.40"
+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/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/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/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