From 4f34a646c5e8b44486e9a247c62e2ecdcea771e1 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Fri, 28 Aug 2026 09:24:35 +0200 Subject: [PATCH] Implement SNP distance models with gamma correction and PHYLIP output Renames the CLI flag from --metric to --distance and introduces eight closed-form SNP distance models with optional Jin-Nei gamma correction. Integrates the ndarray crate for matrix operations and adds relaxed PHYLIP output formatting. Updates architecture and theory documentation to cover the new sparse matrix variants, algorithmic fixes, and distance metric implementations. --- DevDoc/architecture/siblings/index.html | 114 +++ .../implementation/obicompactvec/index.html | 11 +- .../partition_layer_cache/index.html | 882 +++++++++++++++++- DevDoc/sitemap.xml.gz | Bin 127 -> 127 bytes .../theory/evolutionary_distances/index.html | 343 +++++++ DevDocMD/theory/evolutionary_distances.md | 144 ++- src/Cargo.lock | 1 + src/obikmer2/Cargo.toml | 1 + src/obikmer2/src/cmd/phylo/args.rs | 135 ++- src/obikmer2/src/cmd/phylo/mod.rs | 103 +- src/obikmer2/src/cmd/phylo/phylip.rs | 22 + src/obikphylo/src/lib.rs | 1 + src/obikphylo/src/siblings/algorithms/mod.rs | 3 + .../src/siblings/algorithms/pairwise.rs | 60 ++ .../src/siblings/algorithms/snp_distance.rs | 398 ++++++++ .../src/siblings/extensions/sibling_ext.rs | 50 +- src/obikphylo/src/siblings/mod.rs | 4 +- 17 files changed, 2199 insertions(+), 73 deletions(-) create mode 100644 src/obikmer2/src/cmd/phylo/phylip.rs create mode 100644 src/obikphylo/src/siblings/algorithms/snp_distance.rs diff --git a/DevDoc/architecture/siblings/index.html b/DevDoc/architecture/siblings/index.html index a010b3cc..8be09b91 100644 --- a/DevDoc/architecture/siblings/index.html +++ b/DevDoc/architecture/siblings/index.html @@ -1367,6 +1367,17 @@ + + +
  • + + + + PersistentCompactIntMatrix::Sparse — implemented (2026-08-26) + + + +
    • @@ -1652,6 +1663,17 @@ + + +
  • + + + + PersistentCompactIntMatrix::Sparse — implemented (2026-08-26) + + + +
    • @@ -2582,6 +2604,98 @@ mismatches. The dense/sparse performance gap is gone — previously sparse systematic gap. pack --sparse's claimed query win isn't confirmed outright by this (sparse should arguably now beat dense on truly sparse real data, not just tie), but the pathological regression is fixed.

    +

    PersistentCompactIntMatrix::Sparse — implemented (2026-08-26)

    +

    Closes the gap flagged throughout this document ("no sparse count format +exists yet", traits.rs:9-12's "Explicitly deferred"): obicompactvec +already had PersistentSparseCompactIntMatrix (row-major, built on top of +PersistentSparseBitMatrix as its "which columns are non-zero" support, +values not deduplicated — see that struct's own doc comment), but it was +never wired into PersistentCompactIntMatrix, the dense-dispatching enum +every real consumer (TypedLayer<PersistentCompactIntMatrix>, +KmerLayer::Count) actually holds. Concretely: kmer_index.rs:: +pack_matrices(sparse=true) already called pack_sparse_compact_int_matrix +on every layer's counts/ — but PersistentCompactIntMatrix::open had no +code path back to what that just wrote, so a Count layer became +unreadable ("no count matrix found ... run 'obikmer upgrade'") the moment +anyone ran pack --sparse on an index with count layers. Root cause, not a +workaround: add the missing Sparse variant.

    + +

    Neither surfaced through jaccard_dist_matrix/hamming_dist_matrix + (both explicitly zero their own diagonal at finalisation, and the + off-diagonal partial_hamming bug had gone untested against dense + entirely) — only visible to a caller of the raw partial_* methods + directly, which is exactly what partial_threshold_jaccard(1)'s new + shortcut became. Fixed at the source, not patched around at the call + site; regression test added: + tests::sparse::partial_jaccard_and_hamming_match_dense_including_diagonal.

    + diff --git a/DevDoc/implementation/obicompactvec/index.html b/DevDoc/implementation/obicompactvec/index.html index 1f73a5a6..9f6c2fff 100644 --- a/DevDoc/implementation/obicompactvec/index.html +++ b/DevDoc/implementation/obicompactvec/index.html @@ -1168,7 +1168,7 @@ Pass 1 — byte max, SIMD-vectorizable, O(n)

    Matrix types

    -

    Both matrix types are enums behind a transparent API — the caller never matches on the variant. PersistentCompactIntMatrix has two variants (Columnar, Packed). PersistentBitMatrix has four:

    +

    Both matrix types are enums behind a transparent API — the caller never matches on the variant. PersistentCompactIntMatrix has three variants (Columnar, Packed, Sparse). PersistentBitMatrix has four:

    @@ -1185,12 +1185,12 @@ Pass 1 — byte max, SIMD-vectorizable, O(n) - + - - + + @@ -1200,7 +1200,8 @@ Pass 1 — byte max, SIMD-vectorizable, O(n)
    Packedsingle matrix.pbmx mmap filesingle matrix.pbmx/matrix.pcmx mmap file query-optimised, produced by pack_bit_matrix/pack_compact_int_matrix
    Sparse (bit only)sparse_meta.json + PFIV/Elias-Fano component files, row-majorSparsebit: sparse_meta.json + PFIV/Elias-Fano component files, row-major. Int: same support files (built on PersistentSparseBitMatrix internally) plus singleton_values.pciv/multi_values.pciv/multi_offsets for the per-row, non-deduplicated values pack --sparse; see siblings.md for the sparse-vs-dense access-pattern trade-off
    -

    PersistentBitMatrix::open(layer_dir) auto-detects the variant, in order: matrix.pbmx → Packed, presence/meta.json → Columnar, presence/sparse_meta.json → Sparse, layer_meta.json (no presence dir at all) → Implicit. col_view/col/sub_matrix panic on Sparse/Implicit where the operation has no direct-slice equivalent (Sparse is k-mer-major, not column-major; Implicit has no backing storage) — callers needing per-column data on those variants go through row/fill_row.

    +

    PersistentBitMatrix::open(layer_dir) auto-detects the variant, in order: matrix.pbmx → Packed, presence/meta.json → Columnar, presence/sparse_meta.json → Sparse, layer_meta.json (no presence dir at all) → Implicit. PersistentCompactIntMatrix::open(layer_dir) mirrors the same priority order minus Implicit (there's no implicit count matrix — counts always have at least one on-disk column): matrix.pcmx → Packed, counts/meta.json → Columnar, counts/singleton_values.pciv → Sparse. col_view/col/sub_matrix panic on Sparse/Implicit where the operation has no direct-slice equivalent (Sparse is k-mer-major, not column-major; Implicit has no backing storage) — callers needing per-column data on those variants go through row/fill_row.

    +

    Unlike the bit side, PersistentSparseCompactIntMatrix's values are not deduplicated across rows — two rows can share the same non-zero column set (same dict_id in the shared support) while carrying different counts — so its CountPartials impl can't reuse the support's dict-multiplicity shortcut the way BitPartials for PersistentSparseBitMatrix does. It still avoids the naive O(n_cols² × n) column-pair scan via a single row-major pass (row_major_pairwise in sparse_intmatrix.rs), reconstructing the squared-difference formulas (euclidean/relfreq_euclidean/hellinger) from per-column marginals via Σ(a-b)² = Σa²+Σb²-2Σab — see siblings.md's "PersistentCompactIntMatrix::Sparse — implemented" entry for the full derivation.

    col_view(c) returns the appropriate view directly:

    // PersistentBitMatrix
 pub fn col_view(&self, c: usize) -> BitSliceView<'_>
diff --git a/DevDoc/implementation/partition_layer_cache/index.html b/DevDoc/implementation/partition_layer_cache/index.html
index f3f36012..1cce3095 100644
--- a/DevDoc/implementation/partition_layer_cache/index.html
+++ b/DevDoc/implementation/partition_layer_cache/index.html
@@ -1213,6 +1213,122 @@
       
     
   
+
+      
+        
  • 
+  
+    
+      
+        (6) done (2026-08-21): Counter — third algorithm, extracted the same way as Dereplicator
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (7) done (2026-08-21): LayerBuilder — fourth and last pipeline algorithm
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (8) design agreed, not yet implemented (2026-08-21): obikalgorithm::Algorithm trait + obikindexer::extensions — private/public extension-trait split, KmerLayer rename
+      
+    
+  
+  
+    
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (9) done (2026-08-21): obikindexer::extensions::PrivateBuilder — item 1 above, implemented
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (10) done (2026-08-21): obikindex::IndexBuilder — the public counterpart, same session
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (11) done (2026-08-21): KmerIndex/IndexMeta made fully stateless, IndexState moved off sentinel files
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (12) done (2026-08-21): obikalgorithm::Algorithm — the shared trait, resumed and closed in one session
+      
+    
+  
+  
 
    • 
       
         
  • 
@@ -1719,6 +1835,122 @@
       
     
   
+
    • 
+      
+        
  • 
+  
+    
+      
+        (6) done (2026-08-21): Counter — third algorithm, extracted the same way as Dereplicator
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (7) done (2026-08-21): LayerBuilder — fourth and last pipeline algorithm
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (8) design agreed, not yet implemented (2026-08-21): obikalgorithm::Algorithm trait + obikindexer::extensions — private/public extension-trait split, KmerLayer rename
+      
+    
+  
+  
+    
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (9) done (2026-08-21): obikindexer::extensions::PrivateBuilder — item 1 above, implemented
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (10) done (2026-08-21): obikindex::IndexBuilder — the public counterpart, same session
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (11) done (2026-08-21): KmerIndex/IndexMeta made fully stateless, IndexState moved off sentinel files
+      
+    
+  
+  
+
    • 
+      
+        
  • 
+  
+    
+      
+        (12) done (2026-08-21): obikalgorithm::Algorithm — the shared trait, resumed and closed in one session
+      
+    
+  
+  
 
    • 
       
         
  • 
@@ -1853,6 +2085,21 @@ The underlying module boundary and its rationale (Layer tier / Partition
 tier / Index tier, each depending only downward) are unchanged; only the
 crate-vs-module packaging changed. See obikindex::layer
 for the current module doc.
    
+
    Superseded, second event, same day (2026-08-21): obikpartitionner
+and obikderep, the two algorithm crates, are also gone — but unlike
+obikpartition/obilayeredmap above, they were not folded into
+obikindex. They were first (mistakenly) merged into obikindex as an
+algorithms submodule, then corrected into a new sibling crate,
+obikindexer, holding obikindexer::algorithms::{partitionner,
+dereplicator} and depending on obikindex — never the reverse, same
+dependency direction obikpartitionner/obikderep already had. Read
+obikpartitionner::X as obikindexer::algorithms::partitionner::X and
+obikderep::X as obikindexer::algorithms::dereplicator::X throughout
+what follows. The distinction the mistake surfaced, worth keeping: data
+crates (obikindex, holding the index/partition/layer model) merge
+naturally into one crate as submodules; algorithm crates that operate on
+that model from outside stay separate, so the dependency only ever runs
+one way.
    
 
    Status (2026-08-20, latest pass): (1) done — obilayeredmap::Layer
 exists, Mat is gone. (1b) done — Layer::Empty, the first non-ready
 state, added (panics on every read method). (2a) done — the
@@ -1872,7 +2119,24 @@ turned out to be bigger than KmerPartition alone: Layerread it before touching KmerPartition/Layer
-signatures, the shape is fully specified. Earlier mix-up, for
+signatures, the shape is fully specified. (6) done — Counter, a third
+algorithm, extracted from PartitionRouter the same way Dereplicator
+was in (4). (7) done — LayerBuilder, the fourth and last pipeline
+algorithm; the indexing pipeline is now fully decomposed into
+obikindexer::algorithms::{partitionner, dereplicator, counter,
+layer_builder}. (8) design agreed, item 1 done in (9) —
+obikindexer::extensions::PrivateBuilder, private, six construction-only
+KmerIndex methods moved out. (10) done, same session — obikindex::
+IndexBuilder, public, the four maintenance methods
+(clear_output_for_create/create_skeleton/finalize_indexed/state)
+shared with merge/select/rebuild/reindex. Item 2 from (8)
+(obikalgorithm::Algorithm) done in (12) — new crate, type Output +
+fn run(&mut self) -> SKResult<Self::Output>, on_progress moved off
+run()'s signature entirely into a per-algorithm .on_progress(...)
+setter. Note: Layer renamed
+KmerLayer (2026-08-21, outside this conversation). (5) itself still not
+implemented, still first on the "order of remaining work" list. Earlier
+mix-up, for
 context: an earlier
 version of this doc used the name KmerPartition (singular) for what was
 actually the collection type (later renamed KmerPartitions, later
@@ -2515,6 +2779,622 @@ itself or a new crate.
    
    examples were judged not enough to be sure of the shape (Fn+Sync vs
    FnMut callback bound already diverged between the two that exist).
    • 
 
+
    (6) done (2026-08-21): Counter — third algorithm, extracted the same way as Dereplicator
    
+
    Between (5) and this, the user did a session of their own crate
+restructuring (see the two "Superseded" notes at the top of this file):
+obikpartition/obilayeredmap folded into obikindex as submodules
+(obikindex::partition, obikindex::layer), and obikpartitionner/
+obikderep merged into one sibling crate, obikindexer, holding
+obikindexer::algorithms::{partitionner, dereplicator}. (5)'s design
+(Layer/KmerPartition self-naming by number, PartitionRouter's
+&mut → & fix) was not part of that — pure crate/module packaging,
+confirmed by reading the actual code (Layer::open/create still take an
+external dir: &Path, KmerPartition still eagerly opens all layers,
+PartitionRouter still holds &mut KmerIndex). (5) remains exactly as
+specified, not yet implemented.
    
+
    This step: count_kmer (still living on PartitionRouter, per (4)'s own
+"still not done" note) extracted into obikindexer::algorithms::counter::
+Counter, mirroring Dereplicator exactly — third data point for the
+eventual obikalgorithm trait, still not extracted (still only 3 examples
+with 2 different callback bounds; holding off per (5)'s "order of
+remaining work").
    
+
    pub struct Counter<'a> {
+    index: &'a KmerIndex,
+    n_partitions: usize,
+    keep_partial: bool,
+}
+
+impl<'a> Counter<'a> {
+    pub fn new(index: &'a KmerIndex) -> Self;
+    pub fn keep_partial(mut self, v: bool) -> Self;   // setter, mirrors PartitionRouter's style; defaults to false
+    pub fn run(&self, on_progress: Option<impl Fn(Progress) + Sync>) -> SKResult<KmerSpectrum>;
+}
+
    
+
    Same shape as Dereplicator throughout: Fn(Progress) + Sync (not
+FnMut) since counting is also a parallel par_iter over partitions, an
+AtomicU64 position counter incremented from inside the parallel closure
+so progress reports arrive in real time rather than bursting at the end
+once .collect() finishes, total: Some(n_partitions) (known up front).
+KmerSpectrum (the {f0, f1, counts} aggregate) moved from
+partitionner::router to counter, since it's Counter::run's return
+value now, not PartitionRouter's. count.rs/kmer_sort.rs moved
+verbatim from partitionner/ to counter/ (unchanged bodies — only
+count_kmer itself, KmerSpectrum, and the imports they pulled in were
+removed from router.rs).
    
+
    One divergence from Dereplicator: a keep_partial setter exists (no
+equivalent on Dereplicator, which has no setters at all) — a real,
+already-present parameter (keep_intermediate at the CLI), not a
+speculative addition.
    
+
    count_kmer's three former callers (obikmer::cmd::index, obikphylo's
+test harness, obikindexer::algorithms::partitionner's own
+pipeline_counts test helper) all updated to Counter::new(&idx).
+run(...) — the last one simplified further: it used to read back
+kmer_spectrum_raw.json from disk after calling count_partition
+directly (white-box), now it just uses the KmerSpectrum Counter::run
+already returns.
    
+
    Full workspace suite green (cargo check --workspace --all-targets +
+cargo test --workspace, exit code 0), plus an end-to-end CLI smoke test
+against real FASTA data (scatter → dereplicate → count → index-build →
+query) — required every time per (3)'s lesson, and it earned its keep
+again: the very first smoke-test query returned zero matches, which
+looked like a regression until traced to the query sequence itself being
+low-complexity ("GGCCCCCCACG", six same-base runs) and rejected by
+query's own default entropy threshold — nothing to do with this change.
+Re-tested with a different substring, confirmed working (kmer found,
+count matched the index).
    
+
    Still not done: (5) (Layer/KmerPartition redesign, PartitionRouter's
+&mut→&), the future cache crate, build_layers (still a KmerIndex
+inherent method, not an algorithm), and obikalgorithm itself.
    
+
    (7) done (2026-08-21): LayerBuilder — fourth and last pipeline algorithm
    
+
    Closes out the indexing pipeline: build_layers/build_index_layer
+(the last stage still living as KmerIndex inherent methods, flagged as
+inconsistent since (6)) extracted into obikindexer::algorithms::
+layer_builder::LayerBuilder, same two-phase shape as the other three.
    
+
    pub struct LayerBuilder<'a> {
+    index: &'a KmerIndex,
+    n_partitions: usize,
+    min_abundance: u32,
+    max_abundance: Option<u32>,
+    keep_intermediate: bool,
+}
+
+impl<'a> LayerBuilder<'a> {
+    pub fn new(index: &'a KmerIndex) -> Self;
+    pub fn min_abundance(mut self, v: u32) -> Self;
+    pub fn max_abundance(mut self, v: Option<u32>) -> Self;
+    pub fn keep_intermediate(mut self, v: bool) -> Self;
+    pub fn run(&self, on_progress: Option<impl FnMut(Progress) + Send>) -> SKResult<usize>;  // returns total kmers built
+}
+
    
+
    Different from all three prior extractions in one respect, deliberately
+not "fixed" to match them: the actual per-partition construction logic
+(De Bruijn graph from dereplicated superkmers + provisional counts →
+unitigs → MPHF → matrix) stayed put as KmerIndex::build_index_layer/
+remove_build_artifacts (both already pub) — not moved into
+obikindexer. Checked first: unlike dereplicate_partition/
+count_partition (which only ever had one caller), build_index_layer
+depends on several obikindex-internal helpers (graph_pipeline::
+{write_graph_as_unitigs, materialize_layer}, common::olm_to_sk) that
+are pub(crate) and shared with merge/select/rebuild's own
+layer-construction paths — moving build_index_layer out would have
+meant either exporting that internal surface just for this one algorithm
+or duplicating it. Neither was needed: build_index_layer/
+remove_build_artifacts were already public KmerIndex methods, so
+LayerBuilder's job is purely the orchestration around them (scheduling,
+config, progress) — the exact same "algorithm calls already-public
+KmerIndex primitives" shape PartitionRouter/Dereplicator/Counter
+already have, just at a coarser grain for this one stage. This is the
+"is the producer's API actually deficient?" check from
+[[feedback_no_spaghetti_petits_pois]] applied and answered "no" — not
+skipped.
    
+
    Two more real divergences, both forced by PartitionRunner, not
+arbitrary:
+- Uses obikindex::PartitionRunner (NUMA-aware scheduler, already
+  pub used from obikindex) instead of plain rayon::into_par_iter
+  like Dereplicator/Counter — matches what build_layers already used
+  before extraction; this stage is more CPU/memory-intensive per partition
+  (graph construction) than scatter/dereplicate/count.
+- Callback bound is FnMut(Progress) + Send — a third variant, not
+  matching either prior shape. PartitionRunner::run's on_done is
+  invoked from its own single controller thread (never concurrently, so
+  no Sync needed, unlike Dereplicator/Counter's Fn + Sync), but
+  that controller thread is itself std::thread::scope-spawned, so the
+  closure still has to be Send to cross into it — caught immediately by
+  the compiler (cannot be sent between threads safely) when Send was
+  first omitted, not a design guess. obikalgorithm's eventual shared
+  trait now has three real callback-bound data points to reconcile
+  (FnMut alone for PartitionRouter::run's sequential loop, FnMut +
+  Send here, Fn + Sync for Dereplicator/Counter's rayon
+  par_iter), not two.
    
+
    KmerIndex::build_layers deleted outright (KmerIndex stays a pure data
+structure — no compute orchestration methods, consistent with dereplicate/
+count_kmer's removal in (4)/(6)). New KmerIndex::mark_indexed() added,
+symmetric to mark_scattered/mark_counted, replacing the inline
+touch(SENTINEL_INDEXED) that used to live inside build_layers.
+Stage::start("index")/rep.push(...) and the progress_bar/
+"{n} total kmers indexed" log line both moved to cmd/index/mod.rs,
+same pattern as (3)/(4)/(6) — LayerBuilder renders nothing itself, just
+reports Progress.
    
+
    All callers updated: cmd/index/mod.rs (Stage 3), obikphylo's test
+harness (also gained a mark_indexed() call it was missing — harmless
+before since nothing checked IndexState::Indexed in that test, but now
+correct).
    
+
    Full workspace suite green (cargo check --workspace --all-targets +
+cargo test --workspace, exit code 0), plus the end-to-end CLI smoke test
+(scripts/smoke_test_index.sh, built earlier specifically so this
+verification step is a one-liner from now on) — 870 kmers indexed, query
+round-trip confirmed, same numbers as (6).
    
+
    The indexing pipeline is now fully decomposed: obikindexer::
+algorithms::{partitionner, dereplicator, counter, layer_builder}, each a
+new/(setters)/run algorithm operating on a &KmerIndex (or &mut for
+PartitionRouter, not yet fixed — see (5)), KmerIndex itself holding no
+pipeline-orchestration logic anymore. Still not done: (5), the future
+cache crate, obikalgorithm (now unblocked — three real callback-bound
+variants observed, worth revisiting whether a single trait can express
+all three or whether that's itself the answer: it can't, and the trait
+should not force it).
    
+
    (8) design agreed, not yet implemented (2026-08-21): obikalgorithm::Algorithm trait + obikindexer::extensions — private/public extension-trait split, KmerLayer rename
    
+
    Session note: Layer was renamed KmerLayer (user, outside this
+conversation, alongside other naming homogenisation with KmerIndex/
+KmerPartition) — every reference to Layer in this doc from before
+2026-08-21 means today's obikindex::layer::KmerLayer.
    
+
    Why this came up
    
+
    Verifying "does cmd/index now rest entirely on the algorithm structs"
+(it doesn't quite — see below) led to sorting KmerIndex's own methods by
+a criterion the user was explicit is semantic, not mechanical: "les
+méthodes qui, sémantiquement, n'ont pas d'intérêt hors de la construction
+de l'index" (methods that have no semantic interest outside index
+construction) — not "methods only called from cmd/index today," which
+a grep could answer but would miss methods construction-adjacent code
+elsewhere (merge/select/rebuild/reindex) also depends on for the
+same reason.
    
+
    Checked, not assumed (grepped every call site before classifying):
    
+
      
+
    • Construction-only, real candidates for a private extension trait:
+  KmerIndex::{mark_scattered, mark_counted, mark_indexed, write_spectrum,
+  build_index_layer, remove_build_artifacts, clear_output_for_create,
+  create_skeleton, finalize_indexed, state}. The last four are called
+  from merge.rs/select.rs/rebuild.rs/reindex.rs too (as
+  precondition checks — "is my source Indexed?" — or shared
+  skeleton/finalize machinery), not just from the 4-stage pipeline — so
+  this extension trait's scope is "construction of any kind," not
+  narrowly "the initial build pipeline."
      • 
+
    • Looked construction-only by name, checked, and kept on KmerIndex:
+  layer_unitigs_path (unitigs are the only way to recover a built
+  index's kmer sequences — read by rebuild_layer.rs and others, well
+  beyond construction — see [[project_unitigs_always_kept]]),
+  pack_matrices (re-runnable maintenance on an already-finished index
+  via obikmer pack, not just a pipeline step), upgrade_layer_meta
+  (migration, runnable on any existing index at any time).
      • 
+
    
+
    The general pattern (not obikindexer-specific)
    
+
    KmerIndex/KmerPartition/KmerLayer stay generic, in obikindex —
+every domain-specific consumer crate gets to attach its own extension
+trait(s), of two kinds:
    
+
      
+
    • Private (pub(crate), invisible outside the defining crate) — for
+  plumbing only that crate's own algorithms need. obikindexer gets
+  exactly one of these (see below); no public counterpart makes sense for
+  it — "l'index est tellement central que le second trait n'a pas
+  vraiment d'intérêt" for construction specifically: nothing external
+  should ever want to call mark_scattered or build_index_layer.
      • 
+
    • Public — for a genuinely reusable domain extension. The user's own
+  example, found while discussing this, not hypothetical: obikindex/src/
+  index/distance.rs (phylogenetic distance metrics) is currently an
+  impl KmerIndex block inside obikindex itself — under this
+  principle it should be a public extension trait owned by obikphylo
+  instead (distance metrics are a phylo concept, obikindex has no more
+  business defining them than obikindex::layer has defining
+  "family"/"minorant", the reasoning SiblingLayerExt already followed
+  for KmerLayer — see obikphylo/src/siblings/iter.rs). Explicitly
+  deferred — noted here so it isn't lost, not part of this round.
      • 
+
    • The future cache-manager crate (still blocked on (5), see above) will
+  add its own public extension trait mirroring part of KmerIndex's/
+  KmerPartition's own read API in cached form (e.g. a cached
+  .partition(i) that doesn't re-touch disk) — same pattern, third data
+  point once built.
      • 
+
    
+
    Concretely, next to implement (two items, in order)
    
+
      
+
    1. obikindexer::extensions — a private (pub(crate)) extension
+   trait, most likely named something like IndexBuildExt (final name
+   not yet chosen), implemented for KmerIndex, carrying the ten methods
+   listed above, moved out of obikindex::index::{kmer_index,
+   index_layer}. Every algorithm in obikindexer::algorithms::* that
+   currently calls idx.mark_scattered()/etc. keeps the same call syntax
+   (extension trait methods are called the same way as inherent ones,
+   just need the trait in scope) — cmd/index/mod.rs itself would need
+   use obikindexer::extensions::IndexBuildExt; (or the module re-exports
+   it) to keep compiling, since it's the one place outside obikindexer's
+   own algorithms that currently calls mark_scattered/write_spectrum/
+   mark_counted/mark_indexed directly. Not yet decided: exact
+   trait name, whether it's one trait or split further (e.g. sentinel
+   marking vs. skeleton/finalize machinery), and whether merge/select/
+   rebuild/reindex (not yet extracted into algorithms themselves) move
+   onto it now too or keep calling the soon-to-be-inherent-no-longer
+   methods some other way in the meantime — ask before implementing,
+   this changes the blast radius significantly (4 more obikindex
+   internal files depend on clear_output_for_create/create_skeleton/
+   finalize_indexed/state).
      2. 
+
    2. obikalgorithm::Algorithm trait — new crate, the shared trait
+   obikpartitionner→obikindexer merge (session start of 2026-08-21)
+   and (6)/(7) were deliberately building toward, now with four real
+   new/(setters)/run examples and three distinct callback-bound
+   shapes to reconcile (plain FnMut for PartitionRouter, FnMut +
+   Send for LayerBuilder, Fn + Sync for Dereplicator/Counter —
+   see (7)). Exact shape not yet drafted in this doc — do that as its own
+   design pass before coding, same discipline as everything above.
      3. 
+
    
+
    Both items: design only, nothing implemented yet — this section is
+the record to resume from, not a plan already executed.
    
+
    (9) done (2026-08-21): obikindexer::extensions::PrivateBuilder — item 1 above, implemented
    
+
    Renamed from IndexBuilder to PrivateBuilder immediately after (same
+session), freeing the name IndexBuilder for (10)'s public trait — read
+IndexBuilder below as PrivateBuilder throughout this section.
    
+
    Scoped down from (8)'s six-method list to the concrete set that's
+genuinely movable without further ripple — checked, not assumed, before
+writing anything:
    
+
    pub(crate) trait PrivateBuilder {
+    fn mark_scattered(&mut self) -> OKIResult<()>;
+    fn mark_counted(&self) -> OKIResult<()>;
+    fn mark_indexed(&self) -> OKIResult<()>;
+    fn write_spectrum(&self, f0: u64, f1: u64, counts: &BTreeMap<u32, u64>) -> OKIResult<()>;
+    fn build_index_layer(&self, i: usize, min_ab: u32, max_ab: Option<u32>, with_counts: bool, mode: &IndexMode, block_bits: u8) -> Result<usize, SKError>;
+    fn remove_build_artifacts(&self, i: usize);
+}
+impl PrivateBuilder for KmerIndex { ... }
+
    
+
    All six moved bodily out of obikindex::index::{kmer_index, index_layer}
+into obikindexer::extensions (new module, pub(crate)) — index_layer.rs
+is now empty and deleted outright.
+clear_output_for_create/create_skeleton/finalize_indexed/state
+stayed inherent on KmerIndex, per (8)'s reasoning: merge/select/
+rebuild/reindex — living inside obikindex itself — call them too,
+and obikindex can never depend on obikindexer to reach a trait defined
+there. Moving those four is real future work (extract
+merge/select/rebuild/reindex into algorithms first), not part of this
+step.
    
+
    One new, small, deliberate API widening in obikindex: build_index_layer
+depends on three helpers that were pub(crate) to obikindex
+(graph_pipeline::{write_graph_as_unitigs, materialize_layer},
+common::olm_to_sk) — widened to pub (re-exported from obikindex's
+crate root) so obikindexer could reach them. This is exactly the
+"enrich shared/lower-level APIs instead of ad hoc local code" call the
+project's own rules ask for, made explicitly rather than routed around:
+three functions, already generically written (no rewrite needed), now
+serve a second caller instead of being duplicated.
    
+
    Why the trait had to be defined in obikindexer, not obikindex:
+Rust's orphan rule — implementing a trait for a foreign type requires
+either the trait or the type to be local to the current crate. KmerIndex
+is foreign to obikindexer, so the trait must be the local half; if it
+were defined in obikindex instead, pub(crate) there would make it
+invisible to obikindexer too (crate-private means private to that
+crate, not "private except to one named dependent") — the opposite of
+what was wanted.
    
+
    A real design decision made while wiring callers up, not a mechanical
+rename: PrivateBuilder being genuinely pub(crate) to obikindexer
+means obikmer::cmd::index (a different crate) can no longer call
+mark_scattered/mark_counted/mark_indexed/write_spectrum directly —
+it never could have, once privacy was real rather than aspirational. Each
+algorithm now marks its own completion as part of run()/close()
+instead of leaving it to the caller:
+- PartitionRouter::close() (not run()) calls mark_scattered() —
+  close(), not run(), is the actual shared completion point between
+  the file-driven run() path and the manual write/write_batch+
+  close() path low-level callers (tests) use; putting it in run()
+  alone would have silently skipped marking for every caller that never
+  calls run(). run() already calls self.close() at its own end, so
+  this covers both paths through one line, not two.
+- Counter::run calls write_spectrum then mark_counted before
+  returning.
+- LayerBuilder::run calls mark_indexed before returning.
    
+
    cmd/index/mod.rs lost all four direct calls (mark_scattered/
+write_spectrum/mark_counted/mark_indexed) — each stage's if
+idx.state() < IndexState::X { ... } block is now purely "run the
+algorithm," no separate bookkeeping call after it. Confirms, precisely
+this time (checked by re-reading the whole file, not assumed): cmd/index
+now rests on the four algorithms for every read/write of pipeline state
+except KmerIndex::{exists, create, state, n_partitions}, which are
+genuinely index-identity concerns, not construction bookkeeping — the
+original question this whole design pass started from.
    
+
    Same fix applied to obikphylo's test harness (its four explicit
+mark_*/write_spectrum calls removed, relying on the algorithms now
+doing it themselves) — obikindexer::algorithms::partitionner's own
+pipeline_counts test needed no change (never called mark_* directly).
    
+
    Full workspace suite green (cargo check --workspace --all-targets +
+cargo test --workspace, exit code 0), plus the CLI smoke test — 870
+kmers, same as (6)/(7).
    
+
    Still not done at the time of writing: item 2 from (8) (obikalgorithm::
+Algorithm), (5), the future cache crate, the distance.rs →
+obikphylo relocation (noted in (8), explicitly deferred), and
+extracting merge/select/rebuild/reindex into algorithms.
    
+
    (10) done (2026-08-21): obikindex::IndexBuilder — the public counterpart, same session
    
+
    Immediate correction to (9): the private trait built there was renamed
+PrivateBuilder (freeing the name), and the four methods (9) had left
+inherent on KmerIndex — clear_output_for_create/create_skeleton/
+finalize_indexed/state — got their own trait after all: IndexBuilder,
+public, defined in obikindex itself (not obikindexer):
    
+
    pub trait IndexBuilder: Sized {
+    fn clear_output_for_create<P: AsRef<Path>>(output: P, force: bool) -> OKIResult<()>;
+    fn create_skeleton<P: AsRef<Path>>(output: P, meta: &IndexMeta) -> OKIResult<Self>;
+    fn finalize_indexed<P: AsRef<Path>>(output: P, rep: &mut Reporter) -> OKIResult<Self>;
+    fn state(&self) -> IndexState;
+}
+impl IndexBuilder for KmerIndex { ... }
+
    
+
    User's framing: these four are "maintenance", not "scientific computation
+on an index" — a different kind of non-generic-ness than (9)'s six
+(mark_*/write_spectrum/build_index_layer/remove_build_artifacts,
+exclusive to the 4-stage pipeline). Maintenance is used more broadly
+(merge/select/rebuild/reindex), so it gets a real, public trait —
+not folded back into KmerIndex's inherent surface, and not private
+either.
    
+
    Where it lives, and why that's not arbitrary: (9) needed the orphan
+rule to force its trait into obikindexer, to achieve genuine
+crate-private visibility. Here the requirement is the opposite:
+merge.rs/select.rs/rebuild.rs/reindex.rs — the trait's own
+heaviest users — live inside obikindex. A trait they need to reach
+must be local to obikindex (or a crate obikindex itself depends on,
+which doesn't exist for this). So IndexBuilder lives in a new
+obikindex/src/index/builder.rs, pub trait (no orphan-rule tension at
+all here — both trait and type are local to the same crate), re-exported
+from obikindex's crate root alongside PrivateBuilder's sibling
+obikindexer::extensions::PrivateBuilder staying where it is. Two
+traits, two crates, two different reasons, not a contradiction.
    
+
    Blast radius, all inside obikindex plus one external crate: every
+internal caller of these four methods needs the trait imported now that
+they're no longer inherent — merge.rs, select.rs, rebuild.rs,
+reindex.rs (use crate::index::builder::IndexBuilder;) and, externally,
+obikmer::cmd::index::mod (use obikindex::IndexBuilder;, for the three
+idx.state() < IndexState::X resumability checks). Call syntax at every
+site is unchanged (KmerIndex::create_skeleton(...),
+self.state()) — only trait-in-scope requirements are new, which is
+exactly the point: same ergonomics, less surface baked into KmerIndex
+itself.
    
+
    Verification went one step further than (9): beyond
+cargo check --workspace --all-targets + cargo test --workspace +
+scripts/smoke_test_index.sh (all green, 870 kmers again), ran
+obikmer merge end to end on two freshly built indexes (exercises
+clear_output_for_create/finalize_indexed directly, the two methods
+scripts/smoke_test_index.sh itself never touches) — exit 0, pack
+stage completed. Test suite alone would not have caught a regression
+here: no existing test builds two real indexes and merges them through
+the CLI.
    
+
    KmerIndex itself now carries only: identity/config accessors
+(root_path/meta/kmer_size/n_bits/evidence_mode/genomes/...),
+path resolution (partition_dir/index_dir/layer_dir/
+partition_meta/n_layers), and a few index-maintenance operations not
+yet sorted into either trait (layer_unitigs_path, pack_matrices,
+upgrade_layer_meta — see (8)'s "tested and discarded" list; still
+correctly inherent, not construction-only by the semantic criterion) —
+create/open/exists (identity, can't be trait methods needing Self
+before one exists) round that out.
    
+
    (11) done (2026-08-21): KmerIndex/IndexMeta made fully stateless, IndexState moved off sentinel files
    
+
    Triggered mid-discussion of obikalgorithm::Algorithm (still not started —
+see "Still not done" below): user asked why PartitionRouter::new still
+took &mut KmerIndex at all, and questioned whether mark_scattered
+belonged in the algorithm or in cmd/index. Investigation found &mut
+had become newly necessary since (9) — mark_scattered was mutating
+self.meta.genomes in memory so Counter's later write_spectrum call
+(same idx instance) would see the derived label. User's resolution: the
+"disk is truth, stateless" principle already agreed for KmerPartition/
+Layer in (5) (still unimplemented for those two) should extend to
+KmerIndex itself — move IndexState (Empty/Scattered/Counted/
+Indexed) off the three sentinel files (scatter.done/count.done/
+index.done, detected by existence) and into a field of index.meta's
+own JSON, so the mark_* calls become plain disk writes an algorithm can
+legitimately make on &self — no in-memory mutation left to protect.
    
+
    Shape of IndexMeta, per the user's explicit spec: one JSON file per
+index (index.meta), one IndexMeta instance per index, held and
+returned as Arc<IndexMeta> (not &IndexMeta) by KmerIndex::meta().
+config (kmer_size/minimizer_size/n_bits/with_counts/evidence/
+block_bits) is fixed at construction, cached as a pub field (getter
+kept alongside, for symmetry) — "les champs constants restent des champs
+de la structure", read once, never re-read from disk. genomes and
+state are the opposite: no in-memory cache at all, every accessor
+(genomes(), state()) re-reads index.meta from disk, every mutator
+(push_genome/rename_genome/set_genomes/set_state/mark_scattered/
+mark_counted/mark_indexed) does a full read-modify-write of the same
+file. An internal std::sync::RwLock<()> is held across each
+read-modify-write sequence (not just the write) so two callers sharing
+the same Arc<IndexMeta> can't lose an update to each other — this is
+not a cross-process lock (that's obisys::DirLock, already held by
+cmd/index for the whole build); it only serialises access through one
+shared in-process instance.
    
+
    Construction, chain-of-responsibility style, matching (5)'s pattern:
+IndexMeta::create(&KmerIndex, config, genomes) / IndexMeta::open(&KmerIndex)
+ask the index for its own root path rather than taking one directly. Since
+KmerIndex::create doesn't have a complete KmerIndex yet to hand in
+(it's what's being built), added lower-level pub(crate) path-based
+primitives create_at(&Path, ...) / open_at(&Path) that KmerIndex::create/
+open and builder.rs's create_skeleton call directly, bypassing the
+convenience wrappers for that one bootstrap case.
    
+
    The one deliberate exception: select_in_place and reindex
+genuinely rewrite config after an index already exists (output
+type/evidence mode changes in place) — contradicting "config never
+changes" for the general case. Resolved with a separate, explicitly
+rare-labelled IndexMeta::rewrite_config(config, genomes) (preserves
+state, overwrites everything else); callers refresh their own cached
+Arc<IndexMeta> afterward (self.meta = Arc::new(IndexMeta::open(self)?))
+since IndexMeta has no way to reach back into whichever KmerIndex
+holds it.
    
+
    Consequence confirmed, not just hoped for: with mark_scattered no
+longer touching anything in memory, PartitionRouter genuinely never
+needs &mut KmerIndex — PartitionRouter<'a> { index: &'a KmerIndex },
+new(&'a KmerIndex). This is effectively the PartitionRouter half of
+(5)'s "order of remaining work" item done as a side effect; KmerPartition/
+Layer themselves are still unimplemented for (5).
    
+
    Blast radius — much larger than (9)/(10), touched nearly every crate:
+every .meta().genomes/.meta.genomes field access became a fallible
+.genomes()? method call (genomes reads io::Result<Vec<GenomeInfo>>
+now, not a field), and .meta_mut() was removed outright (no more direct
+field mutation from outside IndexMeta). Fixed across:
+- obikindex internals: meta.rs/state.rs/kmer_index.rs/builder.rs
+  (full rewrites), reindex.rs/select.rs (switched to rewrite_config),
+  merge.rs (heaviest single file — genome counts precomputed once per
+  source into a Vec<Vec<GenomeInfo>> up front rather than re-reading
+  index.meta from disk repeatedly through the function, sentinel write
+  replaced with dst2.meta.mark_indexed()), stats.rs, distance.rs,
+  dump.rs, predicate.rs (its IndexMeta-inherent matching_genome_indices/
+  build_group_filter now read genomes fresh internally), mod.rs/lib.rs
+  (sentinel constant re-exports removed — IndexState no longer has
+  SENTINEL_*/detect() at all).
+- obikindexer::extensions::PrivateBuilder: mark_scattered signature
+  dropped &mut self → &self; the four mark_*/write_spectrum bodies
+  became one-line delegations to self.meta().mark_*().
+- obikphylo::siblings: alignment.rs/cardinality.rs/distance.rs/
+  entropy.rs/sankoff_bundle.rs/stats.rs/tests.rs — all had
+  self.meta().genomes.len()-shaped reads, mechanically fixed to
+  .genomes().map_err(OKIError::Io)?.len() (tests: .unwrap()).
+- obikmer::cmd::*: annotate (rewrote its rename path to load genomes
+  once, mutate the in-memory Vec, then idx.meta().set_genomes(...)
+  instead of meta_mut()), filter/pack/dump/unitig/merge/select/
+  phylo (fetch-once-then-use pattern for genome counts/labels),
+  utils/maintenance.rs (run_rename now calls the pre-existing
+  IndexMeta::rename_genome, dropping its own hand-rolled field mutation
+  entirely), index/mod.rs (three idx.state() < IndexState::X
+  resumability checks needed a fallible read — factored into a small
+  current_state(&KmerIndex) -> IndexState helper rather than repeating
+  the same unwrap_or_else three times), query/* (emit_batch's
+  signature changed from &IndexMeta to &[GenomeInfo], and genomes is
+  now fetched once in run() and threaded down through process_chunk
+  as Arc<Vec<GenomeInfo>> rather than re-reading index.meta from disk
+  on every chunk — a deliberate deviation from the "always re-read"
+  default, justified because this is a genuine per-chunk hot path, unlike
+  every other call site touched in this pass).
+- One &IndexMeta-vs-Arc<IndexMeta> argument-type mismatch pattern
+  recurred at several CLI call sites (build_filters/build_specs/
+  emit_batch's original signature) — resolved via Arc's deref
+  coercion (&idx.meta() coerces to &IndexMeta) rather than changing
+  every downstream signature to accept Arc<IndexMeta>.
    
+
    Verification: cargo check --workspace --all-targets and
+cargo test --workspace both green (0 failures) after the full
+propagation, scripts/smoke_test_index.sh green (870 kmers, same as every
+prior round), plus a manual CLI run of index (×2) → merge → select
+→ reindex → utils --new-label (rename) → utils --stats, all exit 0,
+confirming the four most-affected commands (the ones (10)'s verification
+already flagged as under-covered by the automated test suite) still work
+end to end against the new Arc<IndexMeta>/on-disk-IndexState shape.
    
+
    Still not done: (5)'s KmerPartition/Layer self-naming redesign itself
+(only the PartitionRouter-&mut-removal piece landed, as a side
+effect); the distance.rs → obikphylo relocation ((9), explicitly
+deferred); the future cache-manager crate. (8)'s obikalgorithm::
+Algorithm trait, resumed and closed in (12) below.
    
+
    (12) done (2026-08-21): obikalgorithm::Algorithm — the shared trait, resumed and closed in one session
    
+
    Resumed (8)'s point 2 through a point-by-point discussion of what's
+actually common across the four pipeline algorithms, now that (11) made
+KmerIndex itself immutable everywhere. Four sub-points, each closed
+before moving to the next:
    
+
    1. Receiver (&self vs &mut self) — investigated whether (11)'s
+removal of &mut KmerIndex also removed the need for PartitionRouter::
+run to take &mut self. It didn't: PartitionRouter holds real
+per-run state of its own (writers: Vec<Option<SKFileWriter>>, open file
+handles, purely in-process RAM — confirmed by checking where writers is
+stored, nothing to do with KmerIndex/disk truth), unrelated to the
+index. First proposal (wrap writers in RefCell so all four could
+share a uniform &self) was retracted on pushback: manufacturing
+interior mutability with runtime borrow checks to satisfy a cosmetic
+uniformity that Rust doesn't even require is over-engineering — a trait
+method's receiver must match exactly across implementors, but nothing
+stops that shared receiver from being &mut self with three of the four
+implementations simply not using the mutability. Settled: trait declares
+&mut self; Dereplicator/Counter/LayerBuilder (previously &self)
+now also take &mut self, unused.
    
+
    2. path_source as a PartitionRouter setter, not a run() param —
+added a files: Option<Box<dyn Iterator<Item = PathBuf> + Send>> field +
+.files(impl Iterator<Item = PathBuf> + Send + 'static) -> Self setter
+(boxed rather than a generic type parameter on PartitionRouter<'a>:
+negligible cost — one PathBuf per input file, not per k-mer — for a
+much more usable type when passing the builder around). run now does
+self.files.take().ok_or_else(...), erroring if .files(...) was never
+called, instead of taking path_source as a parameter.
    
+
    3. Unifying the three progress-callback bound shapes — reopened, then
+resolved differently than (8) originally framed it. First proposal
+(force everything to FnMut(Progress) + Send) was rejected on the same
+principle as point 1: Dereplicator/Counter's Fn(Progress) + Sync
+isn't arbitrary — their callback is invoked concurrently from multiple
+rayon worker threads, and FnMut requires exclusive access, so forcing
+it would mean wrapping the callback in a Mutex for zero benefit at the
+one real call site (pb.inc(1), already thread-safe). The actual
+resolution: move on_progress off run()'s signature entirely, onto a
+per-algorithm .on_progress(...) setter — same treatment as point 2's
+path_source — so each algorithm keeps its own bound (PartitionRouter:
+FnMut(Progress) + 'a, sequential, no Send needed; LayerBuilder:
+FnMut(Progress) + Send + 'a, crosses into PartitionRunner's
+thread::scope-spawned controller thread once; Dereplicator/Counter:
+Fn(Progress) + Sync + 'a, invoked concurrently from rayon workers).
+This also dissolves the original problem run() had: once the
+callback isn't part of run's signature at all, there's nothing left to
+unify there, and point 4 (below) becomes trivial.
    
+
    4. Output as an associated type, Error fixed — trivial once (3)
+moved the callback out: Error was already uniform (all four return
+obiskio::SKResult<T> = Result<T, SKError>), only Output varied
+(()/()/KmerSpectrum/usize). First cut reused obiskio::SKResult
+directly as the trait's return type — caught and corrected the same
+session: SKError enumerates I/O-specific cases (BadMagic/
+Truncated/Compression/...), meaningless at the level of a generic
+"algorithm" abstraction, and borrowing it made obikalgorithm — meant to
+be minimal and neutral — depend on a low-level I/O crate purely to reuse
+its error type. Textbook instance of the "petits pois" failure mode
+(patch around a convenient existing type instead of asking what this
+crate should actually own). Fixed to a genuinely generic, boxed error
+type owned by obikalgorithm itself:
    
+
    // obikalgorithm — no dependency on obiskio or any other crate
+pub type Error = Box<dyn std::error::Error + Send + Sync>;
+pub type Result<T> = std::result::Result<T, Error>;
+
+pub trait Algorithm {
+    type Output;
+    fn run(&mut self) -> Result<Self::Output>;
+}
+
    
+
    Any concrete error (SKError, std::io::Error, ...) converts
+automatically via ?, through std's own blanket impl<E: Error + Send
++ Sync> From<E> for Box<dyn Error + Send + Sync> — no custom From impl
+needed, no dependency on the crate that defines the concrete error type.
+The four algorithms' run bodies needed no change beyond the signature's
+return type (every existing ? on an SKError-returning subcall keeps
+compiling, converting through the same blanket impl at the boundary).
    
+
    PartitionRouter/Dereplicator/Counter/LayerBuilder each impl
+Algorithm for X<'_> { type Output = ...; fn run(&mut self) -> obikalgorithm::Result<...> { ... } }
+— the old inherent run methods were removed outright (not kept as
+duplicates), so callers now use obikalgorithm::Algorithm; to call
+.run(). Every field-lifetime-bound boxed callback (Box<dyn
+FnMut(Progress) + 'a> etc.) is tied to the algorithm's own 'a (the
+&'a KmerIndex lifetime already on the struct), not 'static — avoids
+forcing callers' progress closures to move-capture (and therefore clone
+or Arc-wrap) local state like TracedBar/EMA-rate accumulators that
+they'd otherwise want to keep using by reference after run() returns.
    
+
    Why a new crate, not a submodule of obikindexer: obikmer::cmd::
+index::mod and obikphylo's own test helpers both need to call .run()
+on these algorithms, so the trait has to be reachable from outside
+obikindexer — putting it in obikindexer itself would work file-wise
+but conflates "the trait every algorithm implements" with "one crate's
+particular four implementations of it", the same reasoning that already
+separated obikindex (data model) from obikindexer (algorithms
+operating on it). obikalgorithm has no dependencies at all (see
+above); obikindexer, obikmer, and obikphylo (dev-dependency, for
+its test helper) all depend on it.
    
+
    Blast radius: obikindexer's four algorithm modules (struct field +
+setter + trait impl each); obikmer::cmd::index::mod (three call sites:
+.on_progress(cb) before .run(), unqualified now that the trait is in
+scope); obikindexer::algorithms::partitionner::tests and obikphylo::
+siblings::tests (both had direct .run(None::<fn(Progress)>)-shaped
+calls needing the same treatment). New obikalgorithm crate registered
+in the workspace Cargo.toml, depended on by obikindexer/obikmer
+(regular) and obikphylo (dev).
    
+
    Verification: cargo check --workspace --all-targets and cargo test
+--workspace both green (0 failures), scripts/smoke_test_index.sh green
+(870 kmers, same as every prior round) — this round didn't repeat the
+manual merge/select/reindex CLI exercise from (11), since nothing
+in this pass touched those commands' code paths (only the four pipeline
+algorithms and cmd::index, already covered by the smoke test). Reverified
+after the obiskio-dependency fix above (same three checks, still green,
+obikalgorithm/Cargo.toml now has zero [dependencies]).
    
+
    Still not done: (5)'s KmerPartition/Layer self-naming redesign; the
+distance.rs → obikphylo relocation ((9), explicitly deferred); the
+future cache-manager crate (mentioned in (8) as a later, mirrored
+extension-trait exercise, not started).
    
 
    The problem
    
 
    Reading a layer's data (MPHF + matrix) is not free: MphfLayer::open mmaps
 mphf.bin plus (evidence.bin/fingerprint.bin + unitigs.bin), and the
diff --git a/DevDoc/sitemap.xml.gz b/DevDoc/sitemap.xml.gz
index f3981d08318db42eeb6724d9e49e46634451d7f7..03a78e939757e3344acaf2aa3115fcca5d967fa9 100644
GIT binary patch
delta 13
Ucmb=gXP58h;Al8MVIq4403Ox^uK)l5

delta 13
Ucmb=gXP58h;9!{AK9Riw02+`3>;M1&

diff --git a/DevDoc/theory/evolutionary_distances/index.html b/DevDoc/theory/evolutionary_distances/index.html
index 542fbc37..77c3d6cf 100644
--- a/DevDoc/theory/evolutionary_distances/index.html
+++ b/DevDoc/theory/evolutionary_distances/index.html
@@ -909,6 +909,56 @@
     
   
   
+
+      
+        
  • 
+  
+    
+      
+        --distance unification: SNP corrections as first-class metrics (implemented, 2026-08-28)
+      
+    
+  
+  
+    
+  
 
    • 
       
         
  • 
@@ -2195,6 +2245,56 @@
     
   
   
+
    • 
+      
+        
  • 
+  
+    
+      
+        --distance unification: SNP corrections as first-class metrics (implemented, 2026-08-28)
+      
+    
+  
+  
+    
+  
 
    • 
       
         
  • 
@@ -4261,6 +4361,249 @@ among the survivors) — a single extra pass is sufficient.
    
 M call at ~1/62 frequency, --iqtree-min-freq 0.05; asserts M absent
 from the written _iqtree_states.csv and A/C still present). Full
 workspace cargo test green.
    
+
    --distance unification: SNP corrections as first-class metrics (implemented, 2026-08-28)
    
+
    Implemented. --metric (renamed --distance — several of
+its existing values, e.g. Bray-Curtis, aren't metrics in the strict sense,
+--metric was a misnomer) gains a family of snp-* values computed from the
+central-position SNP pipeline, routed internally to the sibling-annex
+machinery (PairwiseTally, obikphylo::siblings::algorithms::pairwise)
+instead of cache.distance(...)'s existing per-layer traversal — a different
+code path behind the same CLI surface, not just another branch of one
+formula function.
    
+
    Why unify at the CLI level despite the implementation split: phylogenetically
+a SNP-corrected distance is a distance like any other — NJ/UPGMA are agnostic
+to how the matrix was produced, so exposing it as a special-cased subcommand
+instead of a --distance value would misrepresent its role. The
+implementation divergence (sibling-annex-based vs. plain index scan) is real
+but belongs at the routing layer, invisible to the CLI's own vocabulary.
    
+
    --subsample becomes optional for snp-* distances (it stays mandatory
+for --sankoff/--pseudo-alignment, unrelated commands): absent means
+exhaustive, achieved for free by reusing sample_index's existing
+proportional-per-layer-quota mechanism with n set to the index-wide total
+non-monomorphic-minorant count (already available from the sibling-annex
+stats) — every layer's quota then equals its own full count, giving Bernoulli
+p = 1 everywhere, i.e. every eligible family is drawn. No second,
+exhaustive-only driver needed. Present means sampled, exactly as --sankoff
+already behaves.
    
+
    One shared tally, many derived formulas. PairwiseTally's subst[4][4]
+per-pair substitution counts (plus marginal base frequencies derived from it)
+are the sufficient statistic for every closed-form correction below — each
+is a small pure function PairwiseTally -> Array2<f64>, at the same level as
+the already-implemented raw_snp_distance/base_pair_tally/
+cardinality_tally. No new full scan per formula, whether the tally itself
+was built exhaustively or from a subsample.
    
+
    --raw-snp-counts stays a separate, unrelated flag — same underlying
+tally, but a diagnostic (n_snp/n_shared/n_eligible per genome pair, one
+row per pair) rather than a distance value, and its long-table shape doesn't
+fold into a single N×N matrix the way a distance does. No change to its
+existing CSV format.
    
+
    snp-* distance catalog
    
+
    All closed-form (method-of-moments / direct formula), none requiring
+per-pair or per-tree maximum-likelihood fitting — that excludes HKY85's
+tree-ML usage but not its pairwise estimator, which is closed-form like
+F84/TN93 and is included below. snp- prefix on every CLI value.
    
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
    valuecorrects forinputs beyond raw counts
    snp-rawnothing (uncorrected p-distance)—
    snp-jc (Jukes-Cantor, JC69)multiple substitutions per site—
    snp-k2p (Kimura 2-parameter, K80)+ transition/transversion rate biasts/tv split
    snp-k81 (Kimura 3-parameter, K3ST)+ splits transversions into 2 categoriests/tv split, by category
    snp-f81 (Felsenstein 81)+ unequal base frequencies (no ts/tv split)empirical base freqs
    snp-tajima-nei (Tajima-Nei 1984)same goal as F81 (equal-input model), different formula, better small-sample behaviorempirical base freqs
    snp-t92 (Tamura 3-parameter)K2P + GC-content biasts/tv split, GC content
    snp-f84 (Felsenstein 84)full empirical base freqs + single ts/tv rateempirical base freqs, ts/tv split
    snp-hky85 (Hasegawa-Kishino-Yano, pairwise estimator)same inputs as F84, different formulaempirical base freqs, ts/tv split
    snp-tn93 (Tamura-Nei)full empirical base freqs + separate purine/pyrimidine transition rates + transversion rateempirical base freqs, purine-ts/pyrimidine-ts/tv split
    snp-logdet (LogDet / paralinear)no shared-model or stationarity assumption at all — general divergence-matrix determinantfull empirical 4×4 divergence matrix (already subst[4][4])
    snp-tv (transversions-only p-distance)diagnostic/deep-divergence variant — drops transitions entirely (they saturate first)tv-only counts
    
+
    +Γ rate-heterogeneity modifier, applicable to snp-jc, snp-k2p,
+snp-k81, snp-t92, snp-f84, snp-hky85, snp-tn93 (not snp-raw,
+nothing to correct; not snp-logdet, no standard gamma formulation) — same
+formula as the base correction, weighted by a shape parameter α supplied
+by the user (--gamma-shape <alpha>), not estimated by ML. A modifier on
+existing values, not a separate enum arm per distance.
    
+
    Implemented now: snp-raw, snp-jc, snp-k2p, snp-k81, snp-f81,
+snp-t92, snp-tn93, snp-tv, all with +Γ except raw/tv — see
+"Exact formulas" below. snp-tajima-nei, snp-f84, snp-hky85,
+snp-logdet are catalogued above but not implemented: snp-logdet
+needs the true directional per-pair base co-occurrence matrix
+(PairwiseTally only keeps the symmetrised substitution counts
+BasePairTally itself wants — see snp_distance.rs's own module docs for
+why that loses exactly the compositional-asymmetry information LogDet
+exists to detect), snp-tajima-nei needs each genome's own base
+composition (not the pair-pooled estimate the formulas below use), and
+snp-f84/snp-hky85 had no formula independently verified against a
+primary source at implementation time (unlike every formula below, checked
+line-by-line against ape's own
+src/dist_dna.c, not re-derived from memory). Adding any of these later is
+a new function in obikphylo::siblings::algorithms::snp_distance, plus for
+snp-logdet/snp-tajima-nei a new field on PairStats/a per-genome
+accumulator — not an architecture change.
    
+
    Exact formulas (implemented, 2026-08-28)
    
+
    Sufficient statistic, per genome pair (i, j), from
+PairwiseTally::categories/PairwiseTally::base_freq (base order always
+0=A, 1=C, 2=G, 3=T, matching FamilyMask/STATE_SYMBOL):
    
+
      
+
    • \(n_{ts1}\): A↔G substitutions (purine transitions), \(n_{ts2}\): C↔T
+  (pyrimidine transitions)
      • 
+
    • \(n_{tv1}\): A↔C and G↔T substitutions, \(n_{tv2}\): A↔T and C↔G
+  (Kimura's two transversion categories)
      • 
+
    • \(n_{shared}\): loci where both genomes agree
      • 
+
    • \(L = n_{ts1} + n_{ts2} + n_{tv1} + n_{tv2} + n_{shared}\) (total eligible
+  loci for the pair)
      • 
+
    • \(\pi_A, \pi_C, \pi_G, \pi_T\): pair-pooled base frequencies,
+  \(\pi_a = \dfrac{2 \cdot (\text{agreements on } a) + \sum_b n_{a \leftrightarrow b}}{2L}\)
+  (both genomes' calls at this pair's eligible loci, pooled — Nei & Kumar's
+  standard pairwise estimator, not a whole-index average)
      • 
+
    
+
    Derived proportions used below:
    
+
    \[
+p = \frac{n_{ts1}+n_{ts2}+n_{tv1}+n_{tv2}}{L}, \quad
+P = \frac{n_{ts1}+n_{ts2}}{L}, \quad
+Q = \frac{n_{tv1}+n_{tv2}}{L}, \quad
+Q_1 = \frac{n_{tv1}}{L}, \quad
+Q_2 = \frac{n_{tv2}}{L}, \quad
+P_1 = \frac{n_{ts1}}{L}, \quad
+P_2 = \frac{n_{ts2}}{L}
+\]
    
+
    Every formula below was checked term-by-term against ape's own
+src/dist_dna.c (not re-derived from memory) before being ported to
+obikphylo::siblings::algorithms::snp_distance.
    
+
    snp-raw — uncorrected p-distance:
    
+
    \[
+d_{raw} = p
+\]
    
+
    snp-tv — transversions-only p-distance (deliberately uncorrected —
+dropping transitions, which saturate first, is the correction):
    
+
    \[
+d_{tv} = Q
+\]
    
+
    snp-jc (Jukes-Cantor, JC69):
    
+
    \[
+d_{JC} = -\frac{3}{4} \ln\!\left(1 - \frac{4p}{3}\right)
+\]
    
+
    snp-k2p (Kimura 2-parameter, K80), with \(a_1 = 1-2P-Q\), \(a_2 = 1-2Q\):
    
+
    \[
+d_{K2P} = -\frac{1}{2}\ln a_1 - \frac{1}{4}\ln a_2
+\]
    
+
    snp-k81 (Kimura 3-parameter, K3ST), with \(a_1 = 1-2P-2Q_1\),
+\(a_2 = 1-2P-2Q_2\), \(a_3 = 1-2Q_1-2Q_2\):
    
+
    \[
+d_{K81} = -\frac{1}{4}\left(\ln a_1 + \ln a_2 + \ln a_3\right)
+\]
    
+
    snp-f81 (Felsenstein 81), with \(E = 1 - \left(\pi_A^2+\pi_C^2+\pi_G^2+\pi_T^2\right)\):
    
+
    \[
+d_{F81} = -E \ln\!\left(1 - \frac{p}{E}\right)
+\]
    
+
    snp-t92 (Tamura 3-parameter), with GC content
+\(g = \pi_C+\pi_G\), \(w = 2g(1-g)\), \(a_1 = 1 - \dfrac{P}{w} - Q\),
+\(a_2 = 1-2Q\):
    
+
    \[
+d_{T92} = -w \ln a_1 - \frac{1}{2}(1-w)\ln a_2
+\]
    
+
    snp-tn93 (Tamura-Nei), with purine/pyrimidine pooled frequencies
+\(g_R = \pi_A+\pi_G\), \(g_Y = \pi_C+\pi_T\), and
    
+
    \[
+k_1 = \frac{2\pi_A\pi_G}{g_R}, \quad
+k_2 = \frac{2\pi_C\pi_T}{g_Y}, \quad
+k_3 = 2\left(g_R g_Y - \frac{\pi_A\pi_G\, g_Y}{g_R} - \frac{\pi_C\pi_T\, g_R}{g_Y}\right)
+\]
    
+
    \[
+w_1 = 1 - \frac{P_1}{k_1} - \frac{Q}{2g_R}, \quad
+w_2 = 1 - \frac{P_2}{k_2} - \frac{Q}{2g_Y}, \quad
+w_3 = 1 - \frac{Q}{2g_R g_Y}
+\]
    
+
    \[
+d_{TN93} = -k_1 \ln w_1 - k_2 \ln w_2 - k_3 \ln w_3
+\]
    
+
    +Γ gamma correction (Jin & Nei 1990): every formula above is a
+weighted sum of \(-\ln(x)\) terms; the gamma-corrected version replaces
+each such term with the same weight applied to
+\(\alpha\left(x^{-1/\alpha} - 1\right)\) instead — the standard mechanical
+substitution (as \(\alpha \to \infty\), this expression → \(-\ln(x)\),
+recovering the uncorrected formula exactly). E.g. for JC:
    
+
    \[
+d_{JC,\Gamma} = \frac{3}{4}\,\alpha\left[\left(1-\frac{4p}{3}\right)^{-1/\alpha} - 1\right]
+\]
    
+
    Verified term-by-term against ape's own gamma branches for JC69/K80/F81
+(including K80's two-term form — algebraically identical to the generic
+substitution applied to snp-k2p's own \(a_1\)/\(a_2\) terms above, checked
+both symbolically and numerically before simplifying the implementation to
+share one corrected_log helper across every model rather than
+special-casing K80). K81/T92/TN93's gamma branches follow the same
+mechanical substitution but weren't independently checked against an
+ape-equivalent reference for those three specifically — flagged here, not
+silently assumed correct.
    
+
    Output format: PHYLIP-relaxed by default for the distance matrix
    
+
    Implemented. The primary distance-matrix output
+(_dist.csv today) gains multiple formats: PHYLIP-relaxed becomes the
+default (widely read by external NJ tools — PHYLIP neighbor, FastME,
+T-REX, SplitsTree — relaxed rather than strict to avoid the 10-character
+label truncation, since genome labels here routinely exceed it), a --csv
+flag opts back into the current CSV format, PHYLIP-strict is a possible
+future addition (not now). This changes the default output of every
+existing --distance value (jaccard, hamming, bray-curtis, ...), not just
+the new snp-* ones — accepted explicitly (pre-release, single developer
+user, no external consumers to break). Scoped to the distance matrix only:
+--shared-kmers and --raw-snp-counts are counts, not distances, and keep
+their existing CSV-only format.
    
 
    References
    
 
    The Mash mutation-rate model this discussion contrasts with:
 (Fan et al. 2015; Marbl Lab 2026)1 2.
    
diff --git a/DevDocMD/theory/evolutionary_distances.md b/DevDocMD/theory/evolutionary_distances.md
index bdc6e9e5..7b1a3506 100644
--- a/DevDocMD/theory/evolutionary_distances.md
+++ b/DevDocMD/theory/evolutionary_distances.md
@@ -2182,9 +2182,9 @@ Covered by `iqtree::tests::iqtree_min_freq_folds_rare_states_into_missing`
 from the written `_iqtree_states.csv` and `A`/`C` still present). Full
 workspace `cargo test` green.
 
-## `--distance` unification: SNP corrections as first-class metrics (discussion, 2026-08-28)
+## `--distance` unification: SNP corrections as first-class metrics (implemented, 2026-08-28)
 
-**Decided, not yet implemented.** `--metric` (renamed `--distance` — several of
+**Implemented.** `--metric` (renamed `--distance` — several of
 its existing values, e.g. Bray-Curtis, aren't metrics in the strict sense,
 `--metric` was a misnomer) gains a family of `snp-*` values computed from the
 central-position SNP pipeline, routed internally to the sibling-annex
@@ -2253,9 +2253,147 @@ formula as the base correction, weighted by a shape parameter `α` supplied
 by the user (`--gamma-shape `), not estimated by ML. A modifier on
 existing values, not a separate enum arm per distance.
 
+**Implemented now: `snp-raw`, `snp-jc`, `snp-k2p`, `snp-k81`, `snp-f81`,
+`snp-t92`, `snp-tn93`, `snp-tv`, all with `+Γ` except `raw`/`tv`** — see
+"Exact formulas" below. `snp-tajima-nei`, `snp-f84`, `snp-hky85`,
+`snp-logdet` are catalogued above but **not implemented**: `snp-logdet`
+needs the true *directional* per-pair base co-occurrence matrix
+(`PairwiseTally` only keeps the symmetrised substitution counts
+`BasePairTally` itself wants — see `snp_distance.rs`'s own module docs for
+why that loses exactly the compositional-asymmetry information LogDet
+exists to detect), `snp-tajima-nei` needs each genome's *own* base
+composition (not the pair-pooled estimate the formulas below use), and
+`snp-f84`/`snp-hky85` had no formula independently verified against a
+primary source at implementation time (unlike every formula below, checked
+line-by-line against [ape](https://github.com/emmanuelparadis/ape)'s own
+`src/dist_dna.c`, not re-derived from memory). Adding any of these later is
+a new function in `obikphylo::siblings::algorithms::snp_distance`, plus for
+`snp-logdet`/`snp-tajima-nei` a new field on `PairStats`/a per-genome
+accumulator — not an architecture change.
+
+### Exact formulas (implemented, 2026-08-28)
+
+Sufficient statistic, per genome pair `(i, j)`, from
+`PairwiseTally::categories`/`PairwiseTally::base_freq` (base order always
+`0=A, 1=C, 2=G, 3=T`, matching `FamilyMask`/`STATE_SYMBOL`):
+
+- \(n_{ts1}\): A↔G substitutions (purine transitions), \(n_{ts2}\): C↔T
+  (pyrimidine transitions)
+- \(n_{tv1}\): A↔C and G↔T substitutions, \(n_{tv2}\): A↔T and C↔G
+  (Kimura's two transversion categories)
+- \(n_{shared}\): loci where both genomes agree
+- \(L = n_{ts1} + n_{ts2} + n_{tv1} + n_{tv2} + n_{shared}\) (total eligible
+  loci for the pair)
+- \(\pi_A, \pi_C, \pi_G, \pi_T\): pair-pooled base frequencies,
+  \(\pi_a = \dfrac{2 \cdot (\text{agreements on } a) + \sum_b n_{a \leftrightarrow b}}{2L}\)
+  (both genomes' calls at this pair's eligible loci, pooled — Nei & Kumar's
+  standard pairwise estimator, not a whole-index average)
+
+Derived proportions used below:
+
+\[
+p = \frac{n_{ts1}+n_{ts2}+n_{tv1}+n_{tv2}}{L}, \quad
+P = \frac{n_{ts1}+n_{ts2}}{L}, \quad
+Q = \frac{n_{tv1}+n_{tv2}}{L}, \quad
+Q_1 = \frac{n_{tv1}}{L}, \quad
+Q_2 = \frac{n_{tv2}}{L}, \quad
+P_1 = \frac{n_{ts1}}{L}, \quad
+P_2 = \frac{n_{ts2}}{L}
+\]
+
+Every formula below was checked term-by-term against `ape`'s own
+`src/dist_dna.c` (not re-derived from memory) before being ported to
+`obikphylo::siblings::algorithms::snp_distance`.
+
+**`snp-raw`** — uncorrected p-distance:
+
+\[
+d_{raw} = p
+\]
+
+**`snp-tv`** — transversions-only p-distance (deliberately uncorrected —
+dropping transitions, which saturate first, *is* the correction):
+
+\[
+d_{tv} = Q
+\]
+
+**`snp-jc`** (Jukes-Cantor, JC69):
+
+\[
+d_{JC} = -\frac{3}{4} \ln\!\left(1 - \frac{4p}{3}\right)
+\]
+
+**`snp-k2p`** (Kimura 2-parameter, K80), with \(a_1 = 1-2P-Q\), \(a_2 = 1-2Q\):
+
+\[
+d_{K2P} = -\frac{1}{2}\ln a_1 - \frac{1}{4}\ln a_2
+\]
+
+**`snp-k81`** (Kimura 3-parameter, K3ST), with \(a_1 = 1-2P-2Q_1\),
+\(a_2 = 1-2P-2Q_2\), \(a_3 = 1-2Q_1-2Q_2\):
+
+\[
+d_{K81} = -\frac{1}{4}\left(\ln a_1 + \ln a_2 + \ln a_3\right)
+\]
+
+**`snp-f81`** (Felsenstein 81), with \(E = 1 - \left(\pi_A^2+\pi_C^2+\pi_G^2+\pi_T^2\right)\):
+
+\[
+d_{F81} = -E \ln\!\left(1 - \frac{p}{E}\right)
+\]
+
+**`snp-t92`** (Tamura 3-parameter), with GC content
+\(g = \pi_C+\pi_G\), \(w = 2g(1-g)\), \(a_1 = 1 - \dfrac{P}{w} - Q\),
+\(a_2 = 1-2Q\):
+
+\[
+d_{T92} = -w \ln a_1 - \frac{1}{2}(1-w)\ln a_2
+\]
+
+**`snp-tn93`** (Tamura-Nei), with purine/pyrimidine pooled frequencies
+\(g_R = \pi_A+\pi_G\), \(g_Y = \pi_C+\pi_T\), and
+
+\[
+k_1 = \frac{2\pi_A\pi_G}{g_R}, \quad
+k_2 = \frac{2\pi_C\pi_T}{g_Y}, \quad
+k_3 = 2\left(g_R g_Y - \frac{\pi_A\pi_G\, g_Y}{g_R} - \frac{\pi_C\pi_T\, g_R}{g_Y}\right)
+\]
+
+\[
+w_1 = 1 - \frac{P_1}{k_1} - \frac{Q}{2g_R}, \quad
+w_2 = 1 - \frac{P_2}{k_2} - \frac{Q}{2g_Y}, \quad
+w_3 = 1 - \frac{Q}{2g_R g_Y}
+\]
+
+\[
+d_{TN93} = -k_1 \ln w_1 - k_2 \ln w_2 - k_3 \ln w_3
+\]
+
+**`+Γ` gamma correction** (Jin & Nei 1990): every formula above is a
+weighted sum of \(-\ln(x)\) terms; the gamma-corrected version replaces
+each such term with the same weight applied to
+\(\alpha\left(x^{-1/\alpha} - 1\right)\) instead — the standard mechanical
+substitution (as \(\alpha \to \infty\), this expression → \(-\ln(x)\),
+recovering the uncorrected formula exactly). E.g. for JC:
+
+\[
+d_{JC,\Gamma} = \frac{3}{4}\,\alpha\left[\left(1-\frac{4p}{3}\right)^{-1/\alpha} - 1\right]
+\]
+
+Verified term-by-term against `ape`'s own gamma branches for JC69/K80/F81
+(including K80's two-term form — algebraically identical to the generic
+substitution applied to `snp-k2p`'s own \(a_1\)/\(a_2\) terms above, checked
+both symbolically and numerically before simplifying the implementation to
+share one `corrected_log` helper across every model rather than
+special-casing K80). K81/T92/TN93's gamma branches follow the same
+mechanical substitution but weren't independently checked against an
+`ape`-equivalent reference for those three specifically — flagged here, not
+silently assumed correct.
+
 ### Output format: PHYLIP-relaxed by default for the distance matrix
 
-**Decided, not yet implemented.** The primary distance-matrix output
+**Implemented.** The primary distance-matrix output
 (`_dist.csv` today) gains multiple formats: **PHYLIP-relaxed becomes the
 default** (widely read by external NJ tools — PHYLIP `neighbor`, FastME,
 T-REX, SplitsTree — relaxed rather than strict to avoid the 10-character
diff --git a/src/Cargo.lock b/src/Cargo.lock
index f4877475..cffa8612 100644
--- a/src/Cargo.lock
+++ b/src/Cargo.lock
@@ -1670,6 +1670,7 @@ version = "1.2.2"
 dependencies = [
  "clap",
  "csv",
+ "ndarray",
  "obifastwrite",
  "obikalgorithm",
  "obikdump",
diff --git a/src/obikmer2/Cargo.toml b/src/obikmer2/Cargo.toml
index c8cdd7cb..c5387907 100644
--- a/src/obikmer2/Cargo.toml
+++ b/src/obikmer2/Cargo.toml
@@ -29,6 +29,7 @@ obifastwrite  = { path = "../obifastwrite" }
 obiskbuilder  = { path = "../obiskbuilder" }
 clap          = { version = "4", features = ["derive"] }
 csv           = "1"
+ndarray       = "0.17"
 serde         = { version = "1", features = ["derive"] }
 serde_yaml    = "0.9"
 tracing       = "0.1.44"
diff --git a/src/obikmer2/src/cmd/phylo/args.rs b/src/obikmer2/src/cmd/phylo/args.rs
index 32ae6c2f..5ae930d5 100644
--- a/src/obikmer2/src/cmd/phylo/args.rs
+++ b/src/obikmer2/src/cmd/phylo/args.rs
@@ -1,10 +1,16 @@
 use std::path::PathBuf;
 
 use clap::Args;
-use obikphylo::DistanceMetric;
+use obikphylo::{DistanceMetric, SnpDistanceKind};
 
+/// `--distance` value — either one of `obikphylo::DistanceMetric`'s
+/// whole-index metrics (routed to `IndexCache::distance`) or one of
+/// `obikphylo::SnpDistanceKind`'s `snp-*` corrections (routed to
+/// `SiblingExt::snp_distance`, the sibling-annex pipeline) — two genuinely
+/// different code paths behind one CLI vocabulary, see
+/// `DevDocMD/theory/evolutionary_distances.md`, "`--distance` unification".
 #[derive(clap::ValueEnum, Clone, Copy, Debug)]
-pub enum MetricArg {
+pub enum DistanceArg {
     Jaccard,
     Mash,
     Hamming,
@@ -17,35 +23,69 @@ pub enum MetricArg {
     Hellinger,
     #[value(name = "hellinger-euclidean")]
     HellingerEuclidean,
+    #[value(name = "snp-raw")]
+    SnpRaw,
+    #[value(name = "snp-jc")]
+    SnpJc,
+    #[value(name = "snp-k2p")]
+    SnpK2p,
+    #[value(name = "snp-k81")]
+    SnpK81,
+    #[value(name = "snp-f81")]
+    SnpF81,
+    #[value(name = "snp-t92")]
+    SnpT92,
+    #[value(name = "snp-tn93")]
+    SnpTn93,
+    #[value(name = "snp-tv")]
+    SnpTv,
 }
 
-impl From 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,
-            MetricArg::Euclidean          => DistanceMetric::Euclidean,
-            MetricArg::RelfreqEuclidean   => DistanceMetric::RelfreqEuclidean,
-            MetricArg::Hellinger          => DistanceMetric::Hellinger,
-            MetricArg::HellingerEuclidean => DistanceMetric::HellingerEuclidean,
-        }
+impl DistanceArg {
+    /// `Some` for the whole-index metrics, `None` for `snp-*` values.
+    pub fn as_classic(self) -> Option {
+        Some(match self {
+            DistanceArg::Jaccard            => DistanceMetric::Jaccard,
+            DistanceArg::Mash               => DistanceMetric::Mash,
+            DistanceArg::Hamming            => DistanceMetric::Hamming,
+            DistanceArg::BrayCurtis         => DistanceMetric::BrayCurtis,
+            DistanceArg::RelfreqBrayCurtis  => DistanceMetric::RelfreqBrayCurtis,
+            DistanceArg::Euclidean          => DistanceMetric::Euclidean,
+            DistanceArg::RelfreqEuclidean   => DistanceMetric::RelfreqEuclidean,
+            DistanceArg::Hellinger          => DistanceMetric::Hellinger,
+            DistanceArg::HellingerEuclidean => DistanceMetric::HellingerEuclidean,
+            _ => return None,
+        })
+    }
+
+    /// `Some` for the `snp-*` values, `None` for the whole-index metrics.
+    pub fn as_snp(self) -> Option {
+        Some(match self {
+            DistanceArg::SnpRaw  => SnpDistanceKind::Raw,
+            DistanceArg::SnpJc   => SnpDistanceKind::Jc,
+            DistanceArg::SnpK2p  => SnpDistanceKind::K2p,
+            DistanceArg::SnpK81  => SnpDistanceKind::K81,
+            DistanceArg::SnpF81  => SnpDistanceKind::F81,
+            DistanceArg::SnpT92  => SnpDistanceKind::T92,
+            DistanceArg::SnpTn93 => SnpDistanceKind::Tn93,
+            DistanceArg::SnpTv   => SnpDistanceKind::Tv,
+            _ => return None,
+        })
     }
 }
 
-/// Partial transfer of `obikmer`'s `phylo` command: the plain distance-metric
-/// path (`--metric`/NJ/UPGMA), annex construction (`--sibling-annex`),
-/// annex diagnostics (`--sibling-stats`, `--sibling-hist`), entropy
-/// reporting (`--shannon`), SNP pseudo-alignment sampling
-/// (`--pseudo-alignment`, `--subsample`, `--free-loss`, `--no-ambiguity`,
-/// `--entropy`/`--entropy-sd`), Sankoff cost-matrix calibration
-/// (`--sankoff`, `--sankoff-ratio-ceiling`) and its TNT/PhyG/IQ-TREE exports
-/// (`--tnt`, `--phyg`, `--iqtree`/`--iqtree-min-freq`,
-/// `--sankoff-cost-scale`) — everything else sibling-annex-based (raw SNP
-/// distance, family overlap, ...) stays in `obikmer` until the rest of
-/// `obikphylo::siblings` is reconnected (see the project memory on this).
+/// Partial transfer of `obikmer`'s `phylo` command: the whole-index
+/// `--distance` path (classic metrics + `snp-*` corrections/NJ/UPGMA),
+/// annex construction (`--sibling-annex`), annex diagnostics
+/// (`--sibling-stats`, `--sibling-hist`), entropy reporting (`--shannon`),
+/// SNP pseudo-alignment sampling (`--pseudo-alignment`, `--subsample`,
+/// `--free-loss`, `--no-ambiguity`, `--entropy`/`--entropy-sd`), Sankoff
+/// cost-matrix calibration (`--sankoff`, `--sankoff-ratio-ceiling`) and its
+/// TNT/PhyG/IQ-TREE exports (`--tnt`, `--phyg`, `--iqtree`/
+/// `--iqtree-min-freq`, `--sankoff-cost-scale`) — everything else
+/// sibling-annex-based (family overlap, ...) stays in `obikmer` until the
+/// rest of `obikphylo::siblings` is reconnected (see the project memory on
+/// this).
 #[derive(Args)]
 pub struct PhyloArgs {
     /// Index directory
@@ -97,9 +137,13 @@ pub struct PhyloArgs {
     pub pseudo_alignment: bool,
 
     /// Target number of variable sites to sample index-wide for
-    /// `--pseudo-alignment` — a target, not a guarantee (proportional
-    /// per-layer sampling; see `obikphylo::siblings::SiblingExt::snp_pseudo_alignment`'s
-    /// own docs).
+    /// `--pseudo-alignment`/`--sankoff` (mandatory for both) and for a
+    /// `snp-*` `--distance` value (optional there: omitted means exhaustive
+    /// — every non-monomorphic minorant of the whole index, not an
+    /// approximation, see `obikphylo::siblings::SiblingExt::snp_distance`'s
+    /// own docs) — a target, not a guarantee when given (proportional
+    /// per-layer sampling; see
+    /// `obikphylo::siblings::SiblingExt::snp_pseudo_alignment`'s own docs).
     #[arg(long)]
     pub subsample: Option,
 
@@ -205,14 +249,43 @@ pub struct PhyloArgs {
     #[arg(long, default_value = "100")]
     pub sankoff_cost_scale: f64,
 
-    /// Distance metric to compute
+    /// Distance to compute — either a whole-index metric (`jaccard`,
+    /// `mash`, `hamming`, `bray-curtis`, ...) or a `snp-*` correction over
+    /// the central-position SNP substitution spectrum (`snp-raw`, `snp-jc`,
+    /// `snp-k2p`, `snp-k81`, `snp-f81`, `snp-t92`, `snp-tn93`, `snp-tv`) —
+    /// the latter route to a different computation entirely
+    /// (`SiblingExt::snp_distance`, requires `--sibling-annex` first; see
+    /// `DevDocMD/theory/evolutionary_distances.md`, "`--distance`
+    /// unification" for the full catalog and why LogDet/Tajima-Nei/F84/
+    /// HKY85 aren't offered yet).
     #[arg(long, value_enum, default_value = "jaccard")]
-    pub metric: MetricArg,
+    pub distance: DistanceArg,
+
+    /// Rate-heterogeneity correction (Jin-Nei gamma shape parameter `α`)
+    /// for `snp-*` `--distance` values that support it
+    /// (`obikphylo::SnpDistanceKind::supports_gamma`: every one except
+    /// `snp-raw`/`snp-tv`, which have nothing to correct/are deliberately
+    /// uncorrected). Has no effect on the whole-index metrics. Rejected at
+    /// runtime if given alongside an unsupported `--distance` value.
+    #[arg(long, value_name = "ALPHA")]
+    pub gamma_shape: Option,
 
     /// Minimum count to consider a kmer present when computing Jaccard on count indexes
     #[arg(long, default_value = "1")]
     pub presence_threshold: u32,
 
+    /// Write the primary distance matrix as plain CSV instead of the
+    /// default relaxed-PHYLIP format (`n` on the first line, then one
+    /// `labelvalue...` row per genome — no 10-character label
+    /// truncation, unlike strict PHYLIP, not yet offered here). PHYLIP is
+    /// the default because it's what external NJ tools (PHYLIP `neighbor`,
+    /// FastME, T-REX, SplitsTree) actually read; CSV stays available for
+    /// scripting/inspection. Only affects the primary distance matrix —
+    /// `--shared-kmers` keeps its own CSV-only format regardless of this
+    /// flag.
+    #[arg(long)]
+    pub csv: bool,
+
     /// Also output the shared-kmer count matrix (CSV)
     #[arg(long)]
     pub shared_kmers: bool,
diff --git a/src/obikmer2/src/cmd/phylo/mod.rs b/src/obikmer2/src/cmd/phylo/mod.rs
index a272a987..fc2d62f8 100644
--- a/src/obikmer2/src/cmd/phylo/mod.rs
+++ b/src/obikmer2/src/cmd/phylo/mod.rs
@@ -1,6 +1,7 @@
 mod args;
 mod iqtree;
 mod phyg;
+mod phylip;
 mod sankoff;
 mod tnt;
 
@@ -19,6 +20,7 @@ use tracing::info;
 
 use iqtree::write_iqtree;
 use phyg::write_sankoff_phyg;
+use phylip::write_phylip_relaxed;
 use sankoff::{write_sankoff_alignment_fasta, write_sankoff_matrix_csv, write_sankoff_params};
 use tnt::write_sankoff_tnt;
 
@@ -271,57 +273,100 @@ pub fn run(args: PhyloArgs) {
         }
     }
 
-    info!("computing {:?} distances for {} genome(s)", args.metric, n);
+    // ── Distance computation: classic whole-index metric vs. `snp-*` ───────────
+    // Two genuinely different code paths behind one `--distance` value — see
+    // `args::DistanceArg`'s own docs.
+    let (matrix, shared_kmers) = match args.distance.as_classic() {
+        Some(metric) => {
+            info!("computing {metric:?} distances for {n} genome(s)");
+            let need_shared = args.shared_kmers || args.nj || args.upgma;
+            let t = Stage::start("distance");
+            let result = cache
+                .distance(metric, need_shared, args.presence_threshold)
+                .unwrap_or_else(|e| {
+                    eprintln!("error computing distances: {e}");
+                    std::process::exit(1);
+                });
+            rep.push(t.stop());
+            (result.matrix, result.shared_kmers)
+        }
+        None => {
+            if args.shared_kmers {
+                eprintln!("error: --shared-kmers has no meaning for a snp-* --distance value");
+                std::process::exit(1);
+            }
+            let kind = args.distance.as_snp().expect("DistanceArg is always classic or snp");
+            info!(
+                "computing {kind:?} SNP distance for {n} genome(s){}",
+                match args.subsample {
+                    Some(n) => format!(" (subsampled, target {n} site(s))"),
+                    None => " (exhaustive)".into(),
+                }
+            );
+            let t = Stage::start("snp_distance");
+            let matrix = cache
+                .snp_distance(
+                    kind,
+                    args.subsample,
+                    args.free_loss,
+                    args.no_ambiguity,
+                    &exclude_mask,
+                    entropy_bias,
+                    args.gamma_shape,
+                )
+                .unwrap_or_else(|e| {
+                    eprintln!("error computing SNP distance: {e}");
+                    std::process::exit(1);
+                });
+            rep.push(t.stop());
+            (matrix, None)
+        }
+    };
 
-    let need_shared = args.shared_kmers || args.nj || args.upgma;
-    let t = Stage::start("distance");
-    let result = cache
-        .distance(args.metric.into(), need_shared, args.presence_threshold)
-        .unwrap_or_else(|e| {
-            eprintln!("error computing distances: {e}");
-            std::process::exit(1);
-        });
-    rep.push(t.stop());
-
-    // Rows/columns kept in every matrix CSV below — `cache.distance(...)`
-    // above is computed over every genome regardless; only the writers skip
+    // Rows/columns kept in every matrix output below — the computation
+    // above runs over every genome regardless; only the writers skip
     // excluded ones.
     let kept: Vec = (0..n).filter(|&i| !exclude_mask[i]).collect();
 
-    // ── Distance matrix → CSV ─────────────────────────────────────────────────
-    let write_dist_csv = |w: &mut dyn Write| {
-        write!(w, "genome").unwrap();
-        for &j in &kept { write!(w, ",{}", labels[j]).unwrap(); }
-        writeln!(w).unwrap();
-        for &i in &kept {
-            write!(w, "{}", labels[i]).unwrap();
-            for &j in &kept {
-                write!(w, ",{:.6}", result.matrix[[i, j]]).unwrap();
-            }
+    // ── Distance matrix → relaxed PHYLIP (default) or CSV (`--csv`) ────────────
+    let write_dist = |w: &mut dyn Write| {
+        if args.csv {
+            write!(w, "genome").unwrap();
+            for &j in &kept { write!(w, ",{}", labels[j]).unwrap(); }
             writeln!(w).unwrap();
+            for &i in &kept {
+                write!(w, "{}", labels[i]).unwrap();
+                for &j in &kept {
+                    write!(w, ",{:.6}", matrix[[i, j]]).unwrap();
+                }
+                writeln!(w).unwrap();
+            }
+        } else {
+            write_phylip_relaxed(w, &labels, &kept, &matrix);
         }
     };
 
     match &args.output {
         Some(prefix) => {
-            let path = format!("{}_dist.csv", prefix.display());
+            let suffix = if args.csv { "_dist.csv" } else { "_dist.phy" };
+            let path = format!("{}{suffix}", prefix.display());
             let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| {
                 eprintln!("error creating {path}: {e}");
                 std::process::exit(1);
             }));
-            write_dist_csv(&mut f);
+            write_dist(&mut f);
             info!("distance matrix → {path}");
         }
         None => {
             let stdout = io::stdout();
             let mut out = BufWriter::new(stdout.lock());
-            write_dist_csv(&mut out);
+            write_dist(&mut out);
         }
     }
 
     // ── Shared-kmer matrix → CSV ──────────────────────────────────────────────
     if args.shared_kmers {
-        if let Some(shared) = &result.shared_kmers {
+        if let Some(shared) = &shared_kmers {
             let path = args.output.as_ref()
                 .map(|p| format!("{}_shared.csv", p.display()))
                 .unwrap_or_else(|| "shared.csv".into());
@@ -343,7 +388,7 @@ pub fn run(args: PhyloArgs) {
 
     // ── NJ tree ────────────────────────────────────────────────────────────────
     if args.nj {
-        let tree = neighbor_joining(&result.matrix, &labels).unwrap_or_else(|e| {
+        let tree = neighbor_joining(&matrix, &labels).unwrap_or_else(|e| {
             eprintln!("error computing NJ tree: {e}");
             std::process::exit(1);
         });
@@ -360,7 +405,7 @@ pub fn run(args: PhyloArgs) {
 
     // ── UPGMA tree ───────────────────────────────────────────────────────────────
     if args.upgma {
-        let newick = upgma(&result.matrix, &labels).to_newick();
+        let newick = upgma(&matrix, &labels).to_newick();
         let path = args.output.as_ref()
             .map(|p| format!("{}_upgma.nwk", p.display()))
             .unwrap_or_else(|| "upgma.nwk".into());
diff --git a/src/obikmer2/src/cmd/phylo/phylip.rs b/src/obikmer2/src/cmd/phylo/phylip.rs
new file mode 100644
index 00000000..4c616260
--- /dev/null
+++ b/src/obikmer2/src/cmd/phylo/phylip.rs
@@ -0,0 +1,22 @@
+//! Relaxed-PHYLIP distance-matrix writer — the default output format for
+//! `--distance` (see `args::PhyloArgs::csv`'s own docs for why): `n` on the
+//! first line, then one `labelvalue...` row per genome. "Relaxed"
+//! (unlike strict PHYLIP) means no 10-character label truncation/padding —
+//! just whitespace-separated fields, which every external NJ tool this
+//! targets (PHYLIP `neighbor`, FastME, T-REX, SplitsTree) already reads,
+//! and genome labels here routinely exceed strict PHYLIP's 10 characters.
+
+use std::io::Write;
+
+use ndarray::Array2;
+
+pub(super) fn write_phylip_relaxed(w: &mut dyn Write, labels: &[String], indices: &[usize], matrix: &Array2) {
+    writeln!(w, "{}", indices.len()).unwrap();
+    for &i in indices {
+        write!(w, "{}", labels[i]).unwrap();
+        for &j in indices {
+            write!(w, "\t{:.6}", matrix[[i, j]]).unwrap();
+        }
+        writeln!(w).unwrap();
+    }
+}
diff --git a/src/obikphylo/src/lib.rs b/src/obikphylo/src/lib.rs
index b422a7ab..842313bb 100644
--- a/src/obikphylo/src/lib.rs
+++ b/src/obikphylo/src/lib.rs
@@ -19,4 +19,5 @@ mod tree;
 pub mod siblings;
 
 pub use distance::{DistanceMetric, DistanceOutput, Metrics};
+pub use siblings::SnpDistanceKind;
 pub use tree::{Tree, neighbor_joining, upgma};
diff --git a/src/obikphylo/src/siblings/algorithms/mod.rs b/src/obikphylo/src/siblings/algorithms/mod.rs
index 9ec79f52..71890175 100644
--- a/src/obikphylo/src/siblings/algorithms/mod.rs
+++ b/src/obikphylo/src/siblings/algorithms/mod.rs
@@ -22,6 +22,7 @@ mod masking;
 mod minorant_selection;
 mod pairwise;
 mod sankoff;
+mod snp_distance;
 mod stats;
 mod subsample;
 
@@ -33,6 +34,7 @@ pub use cardcomp::{
 };
 pub use pairwise::{BasePairTally, CardinalityTally, RawSnpDistanceOutput};
 pub use sankoff::SankoffBundle;
+pub use snp_distance::SnpDistanceKind;
 pub use stats::SiblingAnnexStats;
 pub use subsample::{EntropyBias, SurvivingFamily};
 
@@ -41,6 +43,7 @@ pub(crate) use annex::build_layer_sibling_annex;
 pub(crate) use entropy::{ensure_layer_entropy_annex, family_entropy, family_entropy_4};
 pub(crate) use family_scan::{Selection, scan_layer_families};
 pub(crate) use sankoff::sankoff_bundle;
+pub(crate) use snp_distance::snp_distance;
 pub(crate) use stats::{sibling_annex_stats, sibling_family_size_histogram};
 
 /// Whether every layer number in `cache` fits in a `FamilyMask` field
diff --git a/src/obikphylo/src/siblings/algorithms/pairwise.rs b/src/obikphylo/src/siblings/algorithms/pairwise.rs
index d29dccaf..07e5cff1 100644
--- a/src/obikphylo/src/siblings/algorithms/pairwise.rs
+++ b/src/obikphylo/src/siblings/algorithms/pairwise.rs
@@ -175,6 +175,66 @@ impl PairwiseTally {
         }
         CardinalityTally { counts }
     }
+
+    /// Categorised substitution counts for pair `(i, j)`, plus its eligible
+    /// ("shared") locus count — the sufficient statistic every closed-form
+    /// `algorithms::snp_distance` correction is built from. Base order
+    /// `0=A,1=C,2=G,3=T` (same convention as `FamilyMask`/`STATE_SYMBOL`).
+    /// `ts1`/`ts2` split transitions by purine (A↔G) vs pyrimidine (C↔T)
+    /// pair — required for TN93, collapse to a single `ts1+ts2` for
+    /// K80/K81/T92, which don't distinguish them. `tv1`/`tv2` split
+    /// transversions the way Kimura's 3-parameter model does (A↔C/G↔T vs
+    /// A↔T/C↔G) — collapse to `tv1+tv2` for every model that doesn't need
+    /// the distinction (K80, F81, T92, TN93).
+    pub(crate) fn categories(&self, i: usize, j: usize) -> PairCategories {
+        let stats = self.pair(i, j);
+        PairCategories {
+            ts1: stats.subst[0][2],
+            ts2: stats.subst[1][3],
+            tv1: stats.subst[0][1] + stats.subst[2][3],
+            tv2: stats.subst[0][3] + stats.subst[1][2],
+            shared: stats.same_all.iter().sum(),
+        }
+    }
+
+    /// Pooled base composition for pair `(i, j)`, estimated from both
+    /// genomes' observed calls at their `shared`/differing eligible loci —
+    /// `2 * same_all[a]` (each agreement locus contributes `a` to both
+    /// genomes) plus `Σ_b subst[a][b]` (each differing locus contributes
+    /// `a` to whichever of the two genomes carried it — `subst` is
+    /// symmetric by construction, so summing one row already counts each
+    /// such event exactly once, see `reduce_pairwise`'s own docs), over
+    /// `2 * n_eligible_loci` total base observations. Feeds F81/T92/TN93,
+    /// which all correct for base-composition bias using exactly this kind
+    /// of pooled empirical frequency (Nei & Kumar's standard estimator for
+    /// a pairwise comparison, not a whole-index average).
+    pub(crate) fn base_freq(&self, i: usize, j: usize) -> [f64; 4] {
+        let stats = self.pair(i, j);
+        let mut counts = [0u64; 4];
+        for a in 0..4 {
+            counts[a] = 2 * stats.same_all[a] + (0..4).map(|b| stats.subst[a][b]).sum::();
+        }
+        let total: u64 = counts.iter().sum();
+        if total == 0 {
+            return [0.25; 4];
+        }
+        counts.map(|c| c as f64 / total as f64)
+    }
+}
+
+/// See [`PairwiseTally::categories`]'s own docs.
+pub(crate) struct PairCategories {
+    pub ts1: u64,
+    pub ts2: u64,
+    pub tv1: u64,
+    pub tv2: u64,
+    pub shared: u64,
+}
+
+impl PairCategories {
+    pub(crate) fn n_eligible(&self) -> u64 {
+        self.ts1 + self.ts2 + self.tv1 + self.tv2 + self.shared
+    }
 }
 
 /// Raw p-distance restricted to loci that are single-copy in **both**
diff --git a/src/obikphylo/src/siblings/algorithms/snp_distance.rs b/src/obikphylo/src/siblings/algorithms/snp_distance.rs
new file mode 100644
index 00000000..7d6c0877
--- /dev/null
+++ b/src/obikphylo/src/siblings/algorithms/snp_distance.rs
@@ -0,0 +1,398 @@
+//! `snp-*` `--distance` values — closed-form (method-of-moments) pairwise
+//! corrections over the central-position SNP substitution spectrum, all
+//! derived from one shared [`PairwiseTally`] (see its own module docs: a
+//! single `O(n²)` structure, built once, from which every correction below
+//! is cheap post-processing, never a second index scan). See
+//! `DevDocMD/theory/evolutionary_distances.md`, "`--distance` unification"
+//! for the design discussion and the literature this implements
+//! ([`ape`](https://github.com/emmanuelparadis/ape)'s `dist.dna` C source,
+//! `src/dist_dna.c`, verified formula-by-formula against the upstream
+//! implementation rather than re-derived from memory).
+//!
+//! Deliberately **not** included here: LogDet/paralinear (needs the true
+//! *directional* per-pair base co-occurrence matrix — `PairwiseTally`
+//! only keeps the symmetrised substitution counts `BasePairTally` itself
+//! wants, which loses exactly the compositional-asymmetry information
+//! LogDet exists to detect), Tajima-Nei (needs each genome's *own* base
+//! composition, not the pair-pooled estimate `base_freq` below provides),
+//! F84 and the ML-fit HKY85 tree usage (no formula independently verified
+//! against a primary source at implementation time). Adding any of these
+//! later is a new function in this module plus, for LogDet/Tajima-Nei, a
+//! new field on `PairStats`/a per-genome accumulator — not an architecture
+//! change.
+
+use ndarray::Array2;
+use obikidxcache::index_cache::IndexCache;
+use obikindex::OKIResult;
+
+use super::pairwise::PairwiseTally;
+use super::sibling_family_size_histogram;
+use super::subsample::{EntropyBias, sample_index};
+
+/// One `snp-*` `--distance` value. `pub`: part of
+/// [`crate::siblings::extensions::SiblingExt::snp_distance`]'s public
+/// signature.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum SnpDistanceKind {
+    /// Uncorrected p-distance (`snp / (snp + shared)`).
+    Raw,
+    /// Jukes-Cantor (JC69): corrects for multiple substitutions per site,
+    /// equal rates, equal base frequencies.
+    Jc,
+    /// Kimura 2-parameter (K80): + transition/transversion rate bias.
+    K2p,
+    /// Kimura 3-parameter (K81/K3ST): splits transversions into the two
+    /// categories A↔C/G↔T and A↔T/C↔G, each its own rate.
+    K81,
+    /// Felsenstein 81 (F81): JC69 + unequal (pair-pooled empirical) base
+    /// frequencies, no ts/tv distinction.
+    F81,
+    /// Tamura 3-parameter (T92): K80 + GC-content bias.
+    T92,
+    /// Tamura-Nei (TN93): unequal base frequencies + separate purine
+    /// (A↔G) / pyrimidine (C↔T) transition rates + a transversion rate —
+    /// the richest closed-form correction implemented here.
+    Tn93,
+    /// Transversions-only p-distance — diagnostic/deep-divergence variant,
+    /// deliberately uncorrected (dropping transitions, which saturate
+    /// first, is itself the correction).
+    Tv,
+}
+
+impl SnpDistanceKind {
+    /// Whether `--gamma-shape` applies to this correction — every
+    /// rate-based model here except the two with no "-ln(x)" term to
+    /// reinterpret as a gamma mixture (`Raw`, nothing to correct; `Tv`,
+    /// deliberately uncorrected).
+    pub fn supports_gamma(self) -> bool {
+        !matches!(self, SnpDistanceKind::Raw | SnpDistanceKind::Tv)
+    }
+}
+
+/// `-ln(x)`, gamma-mixture corrected when `alpha` is given: replaces the
+/// single-rate `-ln(x)` term with the standard Jin-Nei (1990) gamma
+/// substitution `alpha * (x^(-1/alpha) - 1)` — the mechanical term-wise
+/// rewrite every gamma-corrected distance in the literature (JC+Γ, K80+Γ,
+/// F81+Γ, ...) applies to each log term of its base formula; verified here
+/// against `ape`'s own JC69/K80/F81 gamma branches (`dist_dna.c`), then
+/// generalised the same way to every other log term this module uses
+/// (K81/T92/TN93 each being multi-term sums of exactly this shape).
+fn corrected_log(x: f64, alpha: Option) -> f64 {
+    match alpha {
+        Some(alpha) => alpha * (x.powf(-1.0 / alpha) - 1.0),
+        None => -x.ln(),
+    }
+}
+
+pub(crate) fn raw(tally: &PairwiseTally, i: usize, j: usize) -> f64 {
+    let c = tally.categories(i, j);
+    let l = c.n_eligible();
+    if l == 0 {
+        return f64::NAN;
+    }
+    (c.ts1 + c.ts2 + c.tv1 + c.tv2) as f64 / l as f64
+}
+
+pub(crate) fn tv_only(tally: &PairwiseTally, i: usize, j: usize) -> f64 {
+    let c = tally.categories(i, j);
+    let l = c.n_eligible();
+    if l == 0 {
+        return f64::NAN;
+    }
+    (c.tv1 + c.tv2) as f64 / l as f64
+}
+
+pub(crate) fn jc(tally: &PairwiseTally, i: usize, j: usize, alpha: Option) -> f64 {
+    let p = raw(tally, i, j);
+    0.75 * corrected_log(1.0 - 4.0 * p / 3.0, alpha)
+}
+
+pub(crate) fn k2p(tally: &PairwiseTally, i: usize, j: usize, alpha: Option) -> f64 {
+    let c = tally.categories(i, j);
+    let l = c.n_eligible();
+    if l == 0 {
+        return f64::NAN;
+    }
+    let p = (c.ts1 + c.ts2) as f64 / l as f64;
+    let q = (c.tv1 + c.tv2) as f64 / l as f64;
+    let a1 = 1.0 - 2.0 * p - q;
+    let a2 = 1.0 - 2.0 * q;
+    0.5 * corrected_log(a1, alpha) + 0.25 * corrected_log(a2, alpha)
+}
+
+pub(crate) fn k81(tally: &PairwiseTally, i: usize, j: usize, alpha: Option) -> f64 {
+    let c = tally.categories(i, j);
+    let l = c.n_eligible();
+    if l == 0 {
+        return f64::NAN;
+    }
+    let p = (c.ts1 + c.ts2) as f64 / l as f64;
+    let q = c.tv1 as f64 / l as f64;
+    let r = c.tv2 as f64 / l as f64;
+    let a1 = 1.0 - 2.0 * p - 2.0 * q;
+    let a2 = 1.0 - 2.0 * p - 2.0 * r;
+    let a3 = 1.0 - 2.0 * q - 2.0 * r;
+    0.25 * (corrected_log(a1, alpha) + corrected_log(a2, alpha) + corrected_log(a3, alpha))
+}
+
+pub(crate) fn f81(tally: &PairwiseTally, i: usize, j: usize, alpha: Option) -> f64 {
+    let p = raw(tally, i, j);
+    let freq = tally.base_freq(i, j);
+    let e = 1.0 - freq.iter().map(|f| f * f).sum::();
+    match alpha {
+        Some(a) => e * a * ((1.0 - p / e).powf(-1.0 / a) - 1.0),
+        None => -e * (1.0 - p / e).ln(),
+    }
+}
+
+pub(crate) fn t92(tally: &PairwiseTally, i: usize, j: usize, alpha: Option) -> f64 {
+    let c = tally.categories(i, j);
+    let l = c.n_eligible();
+    if l == 0 {
+        return f64::NAN;
+    }
+    let p = (c.ts1 + c.ts2) as f64 / l as f64;
+    let q = (c.tv1 + c.tv2) as f64 / l as f64;
+    let freq = tally.base_freq(i, j);
+    let gc = freq[1] + freq[2]; // C + G
+    let wg = 2.0 * gc * (1.0 - gc);
+    let a1 = 1.0 - p / wg - q;
+    let a2 = 1.0 - 2.0 * q;
+    wg * corrected_log(a1, alpha) + 0.5 * (1.0 - wg) * corrected_log(a2, alpha)
+}
+
+pub(crate) fn tn93(tally: &PairwiseTally, i: usize, j: usize, alpha: Option) -> f64 {
+    let c = tally.categories(i, j);
+    let l = c.n_eligible();
+    if l == 0 {
+        return f64::NAN;
+    }
+    let freq = tally.base_freq(i, j);
+    let g_r = freq[0] + freq[2]; // A + G (purines)
+    let g_y = freq[1] + freq[3]; // C + T (pyrimidines)
+    let k1 = 2.0 * freq[0] * freq[2] / g_r;
+    let k2 = 2.0 * freq[1] * freq[3] / g_y;
+    let k3 = 2.0 * (g_r * g_y - freq[0] * freq[2] * g_y / g_r - freq[1] * freq[3] * g_r / g_y);
+    let p1 = c.ts1 as f64 / l as f64; // A<->G
+    let p2 = c.ts2 as f64 / l as f64; // C<->T
+    let q = (c.tv1 + c.tv2) as f64 / l as f64;
+    let w1 = 1.0 - p1 / k1 - q / (2.0 * g_r);
+    let w2 = 1.0 - p2 / k2 - q / (2.0 * g_y);
+    let w3 = 1.0 - q / (2.0 * g_r * g_y);
+    k1 * corrected_log(w1, alpha) + k2 * corrected_log(w2, alpha) + k3 * corrected_log(w3, alpha)
+}
+
+fn formula(kind: SnpDistanceKind) -> fn(&PairwiseTally, usize, usize, Option) -> f64 {
+    match kind {
+        SnpDistanceKind::Raw => |t, i, j, _| raw(t, i, j),
+        SnpDistanceKind::Tv => |t, i, j, _| tv_only(t, i, j),
+        SnpDistanceKind::Jc => jc,
+        SnpDistanceKind::K2p => k2p,
+        SnpDistanceKind::K81 => k81,
+        SnpDistanceKind::F81 => f81,
+        SnpDistanceKind::T92 => t92,
+        SnpDistanceKind::Tn93 => tn93,
+    }
+}
+
+/// Build (or reuse) a [`PairwiseTally`] and reduce it to one `n×n` distance
+/// matrix under `kind`. `n`: `Some(target)` samples proportionally per
+/// layer exactly like `--sankoff`/`--pseudo-alignment` (see
+/// `algorithms::subsample::sample_index`'s own docs); `None` means
+/// exhaustive — every non-monomorphic minorant of the whole index, achieved
+/// by handing `sample_index` a quota equal to the index-wide total (from
+/// [`sibling_family_size_histogram`], annex-bits-only, no cross-partition
+/// resolution), which makes every per-layer quota equal to that layer's own
+/// full eligible count, i.e. Bernoulli `p = 1` everywhere — no separate
+/// exhaustive-only driver needed. Entropy-biased sampling
+/// (`entropy_bias`) only makes sense for a genuine subsample, so it's
+/// ignored (forced to `None`) when `n` is `None`, regardless of what the
+/// caller passes.
+///
+/// `gamma_shape`: `--gamma-shape`, `None` disables the correction. Rejected
+/// with [`obikindex::OKIError::InvalidInput`] if given alongside a `kind`
+/// that doesn't support it ([`SnpDistanceKind::supports_gamma`]) — checked
+/// here rather than left to silently no-op.
+pub(crate) fn snp_distance(
+    cache: &IndexCache,
+    kind: SnpDistanceKind,
+    n: Option,
+    free_loss: bool,
+    no_ambiguity: bool,
+    excluded: &[bool],
+    entropy_bias: Option,
+    gamma_shape: Option,
+) -> OKIResult> {
+    if gamma_shape.is_some() && !kind.supports_gamma() {
+        return Err(obikindex::OKIError::InvalidInput(
+            "--gamma-shape has no effect on this --distance value".into(),
+        ));
+    }
+
+    let n_genomes = cache.meta().genomes().len();
+    let mut tally = PairwiseTally::new(n_genomes);
+
+    let (target, entropy_bias) = match n {
+        Some(target) => (target, entropy_bias),
+        None => {
+            let counts = sibling_family_size_histogram(cache)?;
+            let total_eligible = (counts[1] + counts[2] + counts[3]) as usize;
+            (total_eligible, None)
+        }
+    };
+
+    if target > 0 {
+        sample_index(
+            cache,
+            target,
+            free_loss,
+            no_ambiguity,
+            excluded,
+            entropy_bias,
+            |_partition, _layer, survivors| {
+                super::pairwise::reduce_pairwise(&survivors, &mut tally);
+            },
+        )?;
+    }
+
+    let f = formula(kind);
+    Ok(Array2::from_shape_fn((n_genomes, n_genomes), |(i, j)| {
+        if i == j { 0.0 } else { f(&tally, i, j, gamma_shape) }
+    }))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::siblings::FamilyMask;
+    use crate::siblings::algorithms::SurvivingFamily;
+    use super::super::pairwise::reduce_pairwise;
+
+    // Base bit positions, matching `PairwiseTally::categories`' convention
+    // (0=A, 1=C, 2=G, 3=T).
+    const A: u8 = 1 << 0;
+    const C: u8 = 1 << 1;
+    const G: u8 = 1 << 2;
+    const T: u8 = 1 << 3;
+
+    /// Two-genome synthetic tally: one variable family per `(genome0,
+    /// genome1)` pair in `pairs`, plus `n_shared` fully-agreeing families to
+    /// pad the eligible-locus denominator — cycled over all 4 bases (not a
+    /// single one) so the pooled base composition stays non-degenerate
+    /// (`F81`/`T92`/`TN93`'s denominators — `E`/`wg` — are exactly `0` for a
+    /// single-base composition, an edge case irrelevant to real data, not
+    /// something these tests need to exercise).
+    fn tally_2genomes(pairs: &[(u8, u8)], n_shared: usize) -> PairwiseTally {
+        let mut families = Vec::new();
+        for &(a, b) in pairs {
+            families.push(SurvivingFamily {
+                family_idx: families.len(),
+                mask: FamilyMask::EMPTY.with(a.trailing_zeros() as u8).with(b.trailing_zeros() as u8),
+                genome_mask: vec![a, b],
+            });
+        }
+        const BASES: [u8; 4] = [A, C, G, T];
+        for i in 0..n_shared {
+            let b = BASES[i % 4];
+            let other = BASES[(i + 1) % 4];
+            families.push(SurvivingFamily {
+                family_idx: families.len(),
+                mask: FamilyMask::EMPTY.with(b.trailing_zeros() as u8).with(other.trailing_zeros() as u8),
+                genome_mask: vec![b, b],
+            });
+        }
+        let mut tally = PairwiseTally::new(2);
+        reduce_pairwise(&families, &mut tally);
+        tally
+    }
+
+    #[test]
+    fn categories_classify_each_substitution_type() {
+        // A<->G (ts1), C<->T (ts2), A<->C (tv1), A<->T (tv2), + 1 shared.
+        let tally = tally_2genomes(&[(A, G), (C, T), (A, C), (A, T)], 1);
+        let c = tally.categories(0, 1);
+        assert_eq!((c.ts1, c.ts2, c.tv1, c.tv2, c.shared), (1, 1, 1, 1, 1));
+        assert_eq!(c.n_eligible(), 5);
+    }
+
+    #[test]
+    fn identical_genomes_give_zero_distance_under_every_correction() {
+        // No substitutions at all — every formula's "-ln(1)" term is 0.
+        let tally = tally_2genomes(&[], 10);
+        for kind in [
+            SnpDistanceKind::Raw,
+            SnpDistanceKind::Jc,
+            SnpDistanceKind::K2p,
+            SnpDistanceKind::K81,
+            SnpDistanceKind::F81,
+            SnpDistanceKind::T92,
+            SnpDistanceKind::Tn93,
+            SnpDistanceKind::Tv,
+        ] {
+            let d = formula(kind)(&tally, 0, 1, None);
+            assert!(d.abs() < 1e-9, "{kind:?} gave {d}, expected ~0");
+        }
+    }
+
+    #[test]
+    fn raw_matches_hand_computed_ratio() {
+        // 1 substitution (A<->G) out of 5 eligible loci (1 subst + 4 shared).
+        let tally = tally_2genomes(&[(A, G)], 4);
+        assert!((raw(&tally, 0, 1) - 0.2).abs() < 1e-12);
+    }
+
+    #[test]
+    fn jc_matches_hand_computed_correction() {
+        // p = 1/5 = 0.2 -> d = -0.75 * ln(1 - 4*0.2/3) = -0.75 * ln(1 - 4/15).
+        let tally = tally_2genomes(&[(A, G)], 4);
+        let expected = -0.75 * (1.0 - 4.0f64 * 0.2 / 3.0).ln();
+        assert!((jc(&tally, 0, 1, None) - expected).abs() < 1e-12);
+    }
+
+    #[test]
+    fn jc_uncorrected_undershoots_raw_p_distance() {
+        // JC always corrects *upward* from the raw p-distance (multiple-hit
+        // correction only adds distance, never removes it) for any p in its
+        // domain.
+        let tally = tally_2genomes(&[(A, G), (C, T)], 20);
+        let p = raw(&tally, 0, 1);
+        let d = jc(&tally, 0, 1, None);
+        assert!(d > p, "JC-corrected {d} should exceed raw p-distance {p}");
+    }
+
+    #[test]
+    fn gamma_correction_converges_to_uncorrected_as_alpha_grows() {
+        // The Jin-Nei substitution alpha*(x^(-1/alpha) - 1) -> -ln(x) as
+        // alpha -> infinity — a basic sanity check on `corrected_log`
+        // independent of any hand-checked literature value.
+        let tally = tally_2genomes(&[(A, G), (A, C)], 20);
+        let uncorrected = jc(&tally, 0, 1, None);
+        let large_alpha = jc(&tally, 0, 1, Some(1.0e6));
+        assert!(
+            (uncorrected - large_alpha).abs() < 1e-3,
+            "uncorrected {uncorrected} vs large-alpha gamma {large_alpha}"
+        );
+    }
+
+    #[test]
+    fn gamma_shape_rejected_for_unsupported_kind() {
+        assert!(!SnpDistanceKind::Raw.supports_gamma());
+        assert!(!SnpDistanceKind::Tv.supports_gamma());
+        assert!(SnpDistanceKind::Jc.supports_gamma());
+    }
+
+    #[test]
+    fn base_freq_sums_to_one_and_matches_pooled_counts() {
+        // 2 A/G substitutions (genome0=A, genome1=G each time) + 1 shared
+        // A/A family -> pooled bases across both genomes: A appears 4 times
+        // (1 per substitution's genome0 side, 2 from the shared locus's
+        // both genomes), G appears 2 times (1 per substitution's genome1
+        // side), total 6 (2 genomes * 3 loci).
+        let tally = tally_2genomes(&[(A, G), (A, G)], 1);
+        let freq = tally.base_freq(0, 1);
+        assert!((freq.iter().sum::() - 1.0).abs() < 1e-12);
+        assert!((freq[0] - 4.0 / 6.0).abs() < 1e-12); // A
+        assert!((freq[2] - 2.0 / 6.0).abs() < 1e-12); // G
+        assert!((freq[1] - 0.0).abs() < 1e-12); // C
+    }
+}
diff --git a/src/obikphylo/src/siblings/extensions/sibling_ext.rs b/src/obikphylo/src/siblings/extensions/sibling_ext.rs
index 6f4d9140..cc9e1177 100644
--- a/src/obikphylo/src/siblings/extensions/sibling_ext.rs
+++ b/src/obikphylo/src/siblings/extensions/sibling_ext.rs
@@ -9,15 +9,16 @@
 use std::io::{BufWriter, Write};
 use std::path::Path;
 
+use ndarray::Array2;
 use obikidxcache::index_cache::IndexCache;
 use obikindex::{OKIError, OKIResult};
 use obisys::progress_bar;
 
 use crate::siblings::algorithms::{
-    EntropyBias, SankoffBundle, Selection, SiblingAnnexStats, SnpAlignment,
+    EntropyBias, SankoffBundle, Selection, SiblingAnnexStats, SnpAlignment, SnpDistanceKind,
     build_layer_sibling_annex, family_entropy, family_entropy_4, is_fast_mode,
     sankoff_bundle, scan_layer_families, sibling_annex_stats, sibling_family_size_histogram,
-    snp_pseudo_alignment,
+    snp_distance, snp_pseudo_alignment,
 };
 use crate::siblings::extensions::SiblingBuilder;
 use crate::siblings::ENTROPY_ANNEX_FILE_NAME;
@@ -121,6 +122,38 @@ pub trait SiblingExt {
         entropy_bias: Option,
         ratio_ceiling: f64,
     ) -> OKIResult;
+
+    /// A `snp-*` `--distance` matrix — one of the closed-form corrections
+    /// in [`SnpDistanceKind`], over the central-position SNP substitution
+    /// spectrum (requires an already-built sibling annex, run
+    /// [`build_sibling_annex`](Self::build_sibling_annex) first).
+    ///
+    /// `n`: `Some(target)` samples proportionally per layer exactly like
+    /// [`sankoff_bundle`](Self::sankoff_bundle)/
+    /// [`snp_pseudo_alignment`](Self::snp_pseudo_alignment) (`--subsample`);
+    /// `None` means exhaustive — every non-monomorphic minorant of the
+    /// whole index, not an approximation (see
+    /// `algorithms::snp_distance::snp_distance`'s own docs for how this
+    /// reuses the same sampling machinery at `p = 1` rather than a second
+    /// driver). `free_loss`/`no_ambiguity`/`excluded`/`entropy_bias` — same
+    /// meaning as `snp_pseudo_alignment`'s own (the last forced to `None`
+    /// when `n` is `None`: entropy-biasing only makes sense for a genuine
+    /// subsample).
+    ///
+    /// `gamma_shape`: `--gamma-shape`, the Jin-Nei rate-heterogeneity
+    /// correction — `None` disables it, `Some(alpha)` is rejected with
+    /// [`OKIError::InvalidInput`] for a `kind` that doesn't support it
+    /// ([`SnpDistanceKind::supports_gamma`]).
+    fn snp_distance(
+        &self,
+        kind: SnpDistanceKind,
+        n: Option,
+        free_loss: bool,
+        no_ambiguity: bool,
+        excluded: &[bool],
+        entropy_bias: Option,
+        gamma_shape: Option,
+    ) -> OKIResult>;
 }
 
 impl SiblingExt for IndexCache {
@@ -265,4 +298,17 @@ impl SiblingExt for IndexCache {
     ) -> OKIResult {
         sankoff_bundle(self, n, free_loss, no_ambiguity, excluded, entropy_bias, ratio_ceiling)
     }
+
+    fn snp_distance(
+        &self,
+        kind: SnpDistanceKind,
+        n: Option,
+        free_loss: bool,
+        no_ambiguity: bool,
+        excluded: &[bool],
+        entropy_bias: Option,
+        gamma_shape: Option,
+    ) -> OKIResult> {
+        snp_distance(self, kind, n, free_loss, no_ambiguity, excluded, entropy_bias, gamma_shape)
+    }
 }
diff --git a/src/obikphylo/src/siblings/mod.rs b/src/obikphylo/src/siblings/mod.rs
index 35ec8b91..42c1d809 100644
--- a/src/obikphylo/src/siblings/mod.rs
+++ b/src/obikphylo/src/siblings/mod.rs
@@ -35,8 +35,8 @@ pub use siblingannex::{FamilyMask, SiblingAnnex, SiblingAnnexBuilder};
 
 pub use algorithms::{
     BasePairTally, CardinalityTally, EntropyBias, RawSnpDistanceOutput, SankoffBundle,
-    SiblingAnnexStats, SnpAlignment, cardinality_transition_probs, composition_transition_probs,
-    pairwise_cost_matrix,
+    SiblingAnnexStats, SnpAlignment, SnpDistanceKind, cardinality_transition_probs,
+    composition_transition_probs, pairwise_cost_matrix,
 };
 pub use extensions::SiblingExt;