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.
Enum + dispatch (intmatrix.rs): PersistentCompactIntMatrix::Sparse
+ (PersistentSparseCompactIntMatrix), detected in open/detect_storage
+ via a singleton_values.pciv marker (mirrors PersistentBitMatrix's own
+ sparse_meta.json check), reported via storage_kind(). col/
+ col_view/col_persist panic/Unsupported on Sparse, same convention
+ as the bit side. sub_matrix/fill_sub_matrix and nonzero_iter
+ unified the same way PersistentBitMatrix's already are (drain
+ nonzero_iter, one traversal per format — see "Implemented
+ (2026-08-20)" above); nonzero_iter had to become Box<dyn Iterator<...>>
+ for the same reason (Columnar/Packed/Sparse are different concrete
+ types). No change needed in obikindex at all — KmerLayer::Count
+ already only ever holds TypedLayer<PersistentCompactIntMatrix>, so the
+ enum absorbing Sparse fixes the unreadable-layer bug for free, same as
+ PersistentBitMatrix::Sparse already did on the presence side.
CountPartials, non-naive (sparse_intmatrix.rs): unlike
+ PersistentSparseBitMatrix's dict-driven col_weights_and_pair_counts,
+ values here aren't deduplicated (two rows can share the same non-zero
+ column set via the same dict_id while carrying different counts), so
+ the "weight by how many rows share a dict entry" shortcut doesn't carry
+ over. What does: a single row-major pass (row_major_pairwise, decodes
+ each row once via for_each_cell_in_row, nests over that row's own
+ co-present columns) — O(Σ k̄²) over populated rows instead of the naive
+ O(n_cols² × n) column-pair rescan, same complexity class as the bit
+ side minus the dict multiplicity discount. Kernels used: min(a,b)
+ (bray, relfreq-bray — both vanish when either side is absent, so no
+ correction needed), a·b and √(a·b) (euclidean/relfreq-euclidean and
+ hellinger — these do need a correction, reconstructed from per-column
+ marginals via Σ(a-b)² = Σa²+Σb²-2Σab, since (a-0)² = a² ≠ 0 unlike
+ the min-based formulas). threshold_jaccard(1) shortcuts straight to
+ support's own BitPartials::partial_jaccard (threshold 1 is exactly
+ presence); threshold_jaccard(0) is closed-form (every u32 is ≥ 0).
Two pre-existing bugs found and fixed while wiring the threshold==1
+ shortcut (bitmatrix/sparse.rs, BitPartials for
+ PersistentSparseBitMatrix, present since the 2026-08-15 implementation
+ above, never caught because no test compared Sparse's raw partial_*
+ output against dense on real data — only the diagonal-blind
+ jaccard_dist_matrix/hamming_dist_matrix finalisations were tested):
partial_jaccard's diagonal was (0, 2×col_weights[i]) instead of a
+ genuine self-comparison (col_weights[i], col_weights[i]) —
+ col_weights_and_pair_counts's inter never pairs a column with
+ itself by construction.partial_hamming's off-diagonal formula itself was wrong: total -
+ union (count of rows where neither column is present) instead of
+ the actual Hamming distance col_weights[i] + col_weights[j] -
+ 2×inter[i,j] (symmetric-difference size). Only coincides with the
+ correct value when col_weights[i] + col_weights[j] == total, so
+ small/synthetic test data could easily have hidden it.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.
tests::intmatrix::sparse_roundtrip_matches_columnar/
+ sparse_roundtrip_from_packed (the open-dispatch fix, both build
+ paths); tests::intmatrix::sparse_count_partials_match_dense (all six
+ CountPartials formulas, thresholds 0/1/2/3, against Columnar on
+ asymmetric-presence data — this is what caught the diagonal gap in the
+ int side's own new code before it shipped, the same way it exposed the
+ two pre-existing bit-side bugs above); obikindex's
+ count_layer_transparently_reads_sparse_after_pack — the actual
+ end-to-end regression test for the original "layer unreadable after
+ pack --sparse" bug, built → packed sparse → reopened, compared against
+ the pre-pack dense read. cargo test -p obicompactvec -p obikindex:
+ green, no regressions (180 + 12 tests).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:
Packed |
-single matrix.pbmx mmap file |
+single 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-major |
+Sparse |
+bit: 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)
+
+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).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.
+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 f3981d08..03a78e93 100644
Binary files a/DevDoc/sitemap.xml.gz and b/DevDoc/sitemap.xml.gz differ
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.
+
+
+
+value
+corrects for
+inputs 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 bias
+ts/tv split
+
+
+snp-k81 (Kimura 3-parameter, K3ST)+ splits transversions into 2 categories
+ts/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 behavior
+empirical base freqs
+
+
+snp-t92 (Tamura 3-parameter)K2P + GC-content bias
+ts/tv split, GC content
+
+
+snp-f84 (Felsenstein 84)full empirical base freqs + single ts/tv rate
+empirical base freqs, ts/tv split
+
+
+snp-hky85 (Hasegawa-Kishino-Yano, pairwise estimator)same inputs as F84, different formula
+empirical base freqs, ts/tv split
+
+
+snp-tn93 (Tamura-Nei)full empirical base freqs + separate purine/pyrimidine transition rates + transversion rate
+empirical base freqs, purine-ts/pyrimidine-ts/tv split
+
+
+snp-logdet (LogDet / paralinear)no shared-model or stationarity assumption at all — general divergence-matrix determinant
+full 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;