obikmer/docmd/theory/evolutionary_distances.md

# Central-position SNP distance (discussion)

Not implemented. Design discussion for a substitution-rate estimator that
observes SNPs directly from paired-genome k-mer comparison, as an alternative
to Mash's Poisson-Jaccard inference (see [obicompactvec](../implementation/obicompactvec.md)
for the implemented Jaccard/Mash distances).

## Motivation

**Primary intent: restrict the comparison to what is actually comparable.**
Mash's Jaccard is computed over the **union** of both genomes' k-mer content:
anything not identically shared is folded into a single undifferentiated
mass, whether the cause is a point substitution, a genuinely absent
homologous region (lineage-specific content, gene-family expansion, HGT,
genome-size asymmetry), or a diverged paralogous copy. The model then
back-infers a single mutation rate from that mass, silently attributing
non-homology to mutation. The central-SNP approach instead conditions every
comparison on local, positive evidence of homology: a locus only enters the
statistic if its `2m` flanking bases (`m = (k-1)/2`) are found intact in
*both* genomes — genuinely absent or non-homologous content is excluded from
the comparison entirely (neither numerator nor denominator), rather than
silently counted as divergence. This is a conditioning on comparability, not
just a richer summary statistic — see "Statistic and correspondence with
`shared`" below for how it plays out against genome-size asymmetry and
diverged gene families, and "Heterozygosity, ploidy, and consensus-assembly
inputs" for the corresponding paralogy/heterozygosity filter.

**Secondary benefit: access to the substitution's nature.** Because the
central base of an odd-k window is directly observable once the flanks are
confirmed conserved, this also yields more than a rate — the
transition/transversion split — enabling classical corrected distances
(Jukes-Cantor, Kimura 2-parameter, LogDet) that a single Jaccard scalar
cannot support.

## Statistic and correspondence with `shared`

A genomic position `p` is covered by `k` overlapping k-mer windows. Requiring
the substitution to sit at the window's **center** makes exactly one window
per SNP eligible — a 1:1 correspondence between SNP and center-neighbor k-mer
pair, avoiding the ~k-fold overcount of an any-position neighbor search.

A locus with a fully conserved `k`-window (flanks **and** center) is an
exact-shared k-mer at that locus; a locus with conserved flanks but a
substituted center is a "central SNP". Both count each locus exactly once, in
matching units:

```
p_hat[i,j] = SNP[i,j] / (SNP[i,j] + shared[i,j])
```

`p_hat` is `P(center substituted | 2m flanks conserved)`. `shared[i,j]` here
is **not** the general-purpose `shared_kmers` matrix used by Jaccard/Mash
(`--shared-kmers`, `BitPartials::partial_jaccard` /
`CountPartials::partial_threshold_jaccard`) — that matrix counts raw k-mer
identity with no per-genome copy-number constraint, whereas `p_hat`'s
denominator applies the eligibility rule defined below (raw or
paralogy-filtered). Both `SNP` and `shared` are accumulated by the same
sweep, from the same per-locus candidate set (source k-mer + 3 variants),
under the same eligibility rule — see "Locus eligibility" below, and
"Heterozygosity, ploidy, and consensus-assembly inputs" for why the
copy-number constraint matters and what it costs.

**Canonical invariance**: for odd k, the central position maps to itself under
reverse-complement (`m -> k-1-m = m`, base complemented). A transition maps to
a transition, a transversion to a transversion — the transition/transversion
split is well-defined in canonical space.

### Definitions: family, and the canonical form of a family

**Family.** The family of a k-mer `x` is the set of (up to) 4 k-mers sharing
`x`'s `2m` flanking bases, differing only at the central base `m`. Membership
is a property of the flank pattern, not of `x` itself: any of the 4 possible
central substitutions belongs to the same family.

**`central_canonical_neighbors()`** (`obikseq`, `CanonicalKmerOf::central_canonical_neighbors`)
generates all 4 members from any one of them (observed or not), each
independently canonicalised (`.canonical()`, i.e. `min(kmer, revcomp(kmer))`).
This independent canonicalisation is necessary because a central substitution
can flip which orientation is lexicographically smaller — two members of the
same family can end up canonicalised in *different* orientations. Despite
that, the **set** of 4 resulting canonical k-mers is invariant: calling
`central_canonical_neighbors()` on any member of a family — present in the
index or not — yields the same 4 values. This is relied upon throughout the
rest of this document.

**Canonical form of a family.** Because orientation can differ member to
member, "which of the 4 is the reference" cannot be defined relative to
*whichever member happened to be visited first*, nor relative to the
minorant (see below) — both are data-dependent (they depend on what is
actually observed), so using either as the reference would make the
reference itself vary depending on what happens to be present in a given
index. Instead: **the canonical form of a family is, by definition, the
member whose own central base — read in its own already-canonical
orientation — is `A`.** This is well-defined for every family, computed
purely from the flank pattern, whether or not that specific member (or any
member at all) is actually observed anywhere in the index. Concretely: call
`central_canonical_neighbors()` on any member (observed or not) to get the
family's 4 canonical forms; the one among them whose own centre nucleotide is
`A` is the family's canonical form. The other 3 (`C`, `G`, `T`) are labelled
relative to *that* fixed reference, not relative to the calling member's own
orientation.

**Consequence for the minorant.** With this fixed A-referenced labelling,
`minorant` (the smallest raw encoding among the family's *observed* members,
introduced further below) becomes directly computable rather than needing to
be tracked as extra state: regenerate the family's 4 canonical forms from
any member's own k-mer (cheap, no lookup), compare the raw encodings of
whichever are marked present, and take the smallest. No separate stored bit
is required — see Step 2b below, where this replaces the earlier
minorant-bit design.

## Locus eligibility: raw definition vs. paralogy filter

For each k-mer `x` observed in genome A (source, one MPHF slot; the 3
central-position variants generated as in the sweep below): check whether
A's locus (flanks fixed) is resolvable in genome B under one of the 4
central forms.

**Raw / no model.** The locus counts in the denominator iff at least one of
the 4 forms is present in B; it counts in the numerator iff the form found in
B differs from A's own. No constraint on A's or B's own copy number at this
locus. Open question, not resolved: what if **more than one** of the 4 forms
is present in B simultaneously (ambiguous target — count once arbitrarily,
count all, or drop)? The stringent filter below sidesteps the question by
construction rather than answering it.

**Stringent / paralogy-aware.** The locus counts only if exactly one of the
4 forms is present in A **and** exactly one is present in B (`count == 1` at
that slot too, when a count index is available, to also exclude same-allele
duplicates that presence alone cannot see). This drops the raw definition's
ambiguous-B case automatically, at the cost of also dropping heterozygous
sites indiscriminately alongside true duplications (see "Heterozygosity,
ploidy, and consensus-assembly inputs" below).

**Rejected: parsimony-based multiset pairing for multiplicity > 1.** Rather
than dropping ambiguous loci, pair identical alleles between A and B first
(0-mutation explanation preferred), then take `min(unmatched_A, unmatched_B)`
as inferred SNP pairs. Rejected on two grounds: (1) circularity — selecting
pairs by minimal apparent divergence, then measuring divergence on those same
pairs, deflates the estimate by construction, not a neutral heuristic; (2)
the discriminating signal is a single base among 4 possible values, and the
flanks are *already* guaranteed identical for every candidate by
construction (that is how the locus was selected) — no information remains
in a k-mer window to tell which copy in B truly corresponds to which copy in
A once multiplicity > 1 on either side. Any pairing rule invents a
correspondence the data cannot support. Multiplicity > 1 is treated as
non-identifiable, not as a puzzle to solve with a heuristic.

## Multi-genome framing: family as pseudo-alignment column

**Idea.** Instead of resolving locus eligibility and correspondence one
genome pair at a time, treat a family as a column of a pseudo multiple
alignment across *all* genomes simultaneously: for each family, each genome
has either a net single-copy state (`A`/`C`/`G`/`T`, when the genome carries
exactly one of the 4 forms) or "missing" (`?`, multi-copy or absent). Flank
conservation (the `2m` bases fixed by construction) supplies positional
homology for free — the same role a real MSA would play, without alignment
software, gap penalties, or progressive-alignment approximations. Stacking
one such column per family, genomes as rows, produces a genuine SNP
pseudo-alignment matrix, not just a bag of pairwise distances.

**Precedent.** This is the same principle behind reference-free
k-mer-based phylogenomics tools — SKA (Split K-mer Analysis, Harris 2018) and
kSNP: split the k-mer around a variable center, use flank identity to call
homologous columns across arbitrarily many genomes with no reference and no
MSA step, then feed the resulting pseudo-alignment to standard phylogenetic
tools. Landing on the same design independently is a good sign, not a
coincidence.

**Resolves the pairwise-correspondence problem, properly.** The "Rejected:
parsimony-based multiset pairing" case above failed because, with only two
genomes' cardinalities to look at, there is no external constraint to justify
picking one correspondence between leftover alleles over another — `min(a,b)`
is a lower bound dressed up as a point estimate (see the follow-up discussion
on Felsenstein-style parsimony inconsistency: minimum-event explanations are
systematically biased low whenever homoplasy/multiplicity is real, not
noise-cancelling). With `N` genomes and many families jointly, the same
question can be answered the way real phylogenetics answers it: ancestral
state reconstruction / ML mapping over a tree estimated from the whole
column set. The tree supplies the missing constraint that two isolated
columns cannot — this is the principled way out, not a heuristic replacement
for one.

**Relation to what's already implemented.** `KmerIndex::raw_snp_distance`
already computes, internally, per family, exactly this row — `single_form:
Vec<Option<u8>>`, one entry per genome, `None` where ambiguous/absent —
before immediately collapsing it into pairwise `snp[i,j]`/`shared[i,j]`
tallies. The pivot this section proposes is small at the implementation
level: stop collapsing early, and surface the per-family row as a first-class
artifact (a `families x genomes` matrix). Pairwise raw p-distance becomes one
projection of that matrix (what's computed today), not the primary object;
downstream, the matrix itself could feed real phylogenetic tools (parsimony/
ML, e.g. RAxML/IQ-TREE-style) instead of only NJ/UPGMA on a homemade
pairwise-distance matrix.

**Caveat: column completeness shrinks with `N`.** The probability that a
family's flanks stay intact simultaneously across all `N` genomes decays with
`N` (same ascertainment-bias mechanism as Bias 1 above, compounded over more
genomes) — fully-resolved columns (no `?` anywhere) become rare as more
genomes are added. Same missing-data situation any real multi-species
alignment faces, and phylogenetic tools already handle it well; the practical
implication is that columns should be allowed partial coverage (>=2 resolved
genomes, not unanimous) rather than requiring every genome to be net
single-copy at that locus.

## Context, detectability, and a 3-way ordinal distance per pair

Empirical follow-up to the pseudo-alignment idea above: `obikmer distance
--snp` was run on a real 20-genome benchmark index and the resulting FASTA
fed to `raxml-ng`. Two problems surfaced, both traced back to conflating
distinct notions under one symbol.

**"Context", precisely.** Sharing a central base between two genomes is not
just sharing a nucleotide — it is sharing a **context**: the `2m` flanking
bases, identical, which is a homology claim about that flanked window
(guaranteed non-coincidental by k-specificity, Bias 4 above), *not* a claim
about orthology or paralogy of the copy each genome carries. `A` opposite `C`
= same context, divergent centre. `A` opposite nothing = **this context is
not observed in one of the two genomes** — informative, not neutral.

**Why the IUPAC/DNA encoding used for the first `--snp` test was wrong.**
Feeding IUPAC-coded ambiguity into a standard DNA model (`raxml-ng --model
GTR+G`) is a semantic mismatch: Felsenstein-pruning ML treats an ambiguous
tip as "exactly one true state, unknown which" (a uniform partial-likelihood
vector over compatible bases), not "these states are simultaneously
present". The two encodings look identical (same IUPAC letters) but the
software reads them backwards from what was intended — this invalidates the
literal branch lengths from that first experiment (topology-level groupings
by genus were still informative, see the worked example further down).

**Why `-` (absence) must not be scored as similarity, but also must not be
scored as a shared character between two absences.** Two genomes both
lacking a context are not observed to resemble each other at that locus —
neither is observed to differ from the other either. It is a symmetric
non-observation, uninformative for that pair, and should contribute nothing
(not a small positive nor a small negative signal) to their distance. A
genome carrying a state (`A`) against one carrying none is a different case
entirely: informative, and should not be scored as neutral "missing data"
the way a generic DNA/ML pipeline would.

**Detectability vs existence — a deliberate simplification, accepted.**
"Context not observed" conflates two different biological events: (1) true
loss of the locus, (2) the locus still exists but a mutation/indel *outside*
the centre, anywhere in the `2m` flanks, broke k-mer recognition. The design
adopts a rigorist stance on purpose: any flank-breaking mutation counts as
"this context no longer exists", full stop — because both causes (1) and (2)
independently require *at least* one more mutational event than a lone
central substitution would. This licenses treating "context absent in one of
the two genomes" as a **lower bound** on distance strictly greater than a
plain central SNP, without needing to know which of the two causes applies.
Coarser than a true event count, and accepted as such (fine substitution-type
modelling, e.g. transition/transversion weighting, is a secondary
refinement, not required for this to be useful).

**Resulting ordinal distance between two genomes at one family/context:**

| Comparison | Distance | Meaning |
|---|---|---|
| same centre (`A`/`A`) | `0` | identical |
| different centre, both single-copy (`A`/`C`) | `1` | plain central SNP |
| one genome has a state, the other has none | `>1` (lower bound) | context undetectable in one genome — at least one extra mutational event, of unknown type |
| neither genome has any state (`∅`/`∅`) | excluded | symmetric non-observation, not comparable, contributes nothing |

This is a direct extension of `KmerIndex::raw_snp_distance` (`obikindex/src/siblings.rs`),
which today only implements the `0`/`1` rows and silently drops everything
else (including the informative `>1` row) rather than scoring it.

**Open, not yet resolved:**
- Calibrating `>1` to a real number for tools expecting continuous distances
  (NJ/UPGMA, ML branch lengths), rather than an arbitrary placeholder.
  Natural route: estimate `p_hat` from the resolved (`0`/`1`) sites first,
  then use the already-derived ascertainment formula (`P(usable window
  showing a central SNP) = p * (1-p)^(2m)`, Bias 1 above) to derive a
  model-consistent value for the `>1` bucket instead of guessing a constant.
- Where multi-copy/ambiguous states (the IUPAC case: a genome carrying more
  than one form) fit into this ordinal scheme — plausibly also `>1` by the
  same "at least one extra event" argument (a second form appearing is a
  gain, itself an event), but not yet worked out.

**Practical alternative validated for the pseudo-alignment output itself**
(orthogonal to the ordinal-distance question above, useful regardless of how
`>1` ends up calibrated): re-encode each family as 4 independent binary
presence/absence characters (`A`,`C`,`G`,`T` columns) instead of one IUPAC
column, feed to a `BIN`-type model instead of `DNA`. `∅` becomes an explicit
`0000` state (identity with another `0000`, not missing data) rather than a
gap — removes the semantic mismatch above by construction. Known cost,
accepted for now: a plain central substitution (`A` -> `C`) becomes 2 binary
flips (`1000` -> `0100`), overweighting substitutions relative to true
gain/loss events, and the 4 sub-characters of one family are not
statistically independent the way a generic `BIN` model assumes. A proper
fix (single 16-state alphabet, i.e. the powerset of `{A,C,G,T}`, with a
substitution-rate structure that respects the subset lattice rather than a
fully general 16x16 GTR-analogue) is very likely not expressible in
`raxml-ng`'s `MULTI` datatype as-is (Mk or fully-general rates only) and a
fully general 16-state rate matrix is almost certainly unidentifiable here
(states of cardinality >=3 are ~2% of sites in the benchmark run). Treated as
a longer-term research question, not a near-term implementation target.

## Sankoff parsimony as the resolution of the 16-state model problem

The "longer-term research question" just above (a 16-state alphabet — the
powerset of `{A,C,G,T}` — with a substitution structure that respects the
subset lattice) turns out to have a near-term answer, once the *unification*
question below is worked through.

**Distance methods (NJ/UPGMA/ME) vs. character methods (parsimony/ML): not a
deep philosophical divide, but a real practical distinction for this
project.** Historically "phenetic" (characters -> distances -> tree) and
"cladistic" (characters -> tree directly) approaches were presented as
opposed schools; the modern view is a mathematical continuity, not a
dichotomy — Minimum Evolution (ME: find the tree minimizing total branch
length from a distance matrix) and Maximum Parsimony (MP: find the tree
minimizing total character-state changes) are both instances of "minimise a
global explanatory cost", and coincide under simple encodings (see Farris
1983, "The Logical Basis of Phylogenetic Analysis"; the MP/ME connection is
developed in the Minimum Evolution / Balanced Minimum Evolution literature,
e.g. Nei and colleagues — citations not independently re-verified here, flag
before quoting further). NJ's own agglomeration step already uses the whole
distance matrix jointly (the Q-matrix), not just the pair being merged — an
earlier claim in this discussion that distance methods are "blind" to
cross-taxon structure at every stage was too strong.

What *does* remain a real, structural distinction for this project: in a
character method, a given character's cost is **re-evaluated per candidate
topology** during tree search (the same family can cost 1 change under one
topology, 2 under another). In a pairwise-distance pipeline
(`raw_snp_distance` as it exists today), each family's contribution to
`d(i,j)` is computed **once**, independent of any candidate topology, before
NJ/UPGMA ever runs — so a question like "does this shared `∅` look like a
synapomorphy under topology T" can never be posed in that pipeline, for any
T. That question is only answerable by a method that tests candidate
topologies and re-scores characters under each — i.e. a character method.

**Sankoff parsimony directly resolves the `∅`/gain-loss/substitution
question, without the identifiability problem of a fitted 16-state model.**
Sankoff's algorithm generalises Fitch parsimony to an arbitrary
user-supplied cost matrix between states (`obikseq`/`obikindex` would treat
each family as a `2^{4}`-state character, state = subset of `{A,C,G,T}`
observed in that genome, `∅` included as a real state, not a gap). The
previous 16-state idea failed specifically because *fitting* a full 16x16
rate matrix by ML is unidentifiable at this data volume; Sankoff sidesteps
that because the cost matrix is **fixed a priori from domain knowledge**, not
estimated — e.g. `c({A},{C}) = 1` (a substitution), `c({A},{A,C}) = 1` (a
gain), `c({A,C},{A}) = 1` (a loss), `c({A,C},{G,T}) = 2` (two changes) — no
estimation, no overparameterisation. This reframes gain/loss and central
substitution as two cost categories with independently chosen weights,
exactly the "two families of parameters" (`mu_substitution`, `mu_gain/loss`)
floated earlier in this discussion, now with an actual algorithmic home.

**Caveat, not blocking for this project's scope.** Sankoff is still
parsimony: in principle exposed to Felsenstein's statistical-inconsistency
result under long-branch attraction (already invoked earlier against
`D_F = min(a,b)`) — parsimony and ML only provably coincide in the
short-branch regime. This is not a practical concern here because it is
exactly this estimator's declared target (closely related genomes, short
branches) — the regime where parsimony's known failure mode does not apply —
but worth stating explicitly as a scope guard rather than leaving it
implicit.

**Cheapest next experiment: don't write a Sankoff tree-search engine, use
one that exists.** The hard part of a from-scratch implementation is not the
Sankoff DP itself (a straightforward dynamic program over a *fixed* tree) but
the topology search (SPR/NNI with incremental re-scoring) that comes for
free with `raxml-ng` on the ML side. **TNT** (Tree analysis using New
Technology, free, standard in morphological cladistics) already implements
Sankoff parsimony with a custom cost matrix plus topology search — the
family-state matrix (already close to what `--snp` produces, minus the
IUPAC/DNA-model mismatch) could be fed there directly, no new code required,
before considering a bespoke engine.

**The concrete comparison this unlocks:** run both pipelines on the same
family data —
`k-mer families -> D_ij -> NJ/ME/UPGMA` (phenetic, what exists today) vs.
`k-mer families -> characters -> argmin_T Sankoff-cost(T)` (cladistic, via
TNT) — and compare the resulting topologies. Agreement would validate that
the pairwise-distance projection preserves the phylogenetic signal;
disagreement would pinpoint exactly what the projection to a single number
per pair loses. Not yet run.

### A concrete Sankoff cost matrix for the 16-state alphabet

**Why `c_gl` and `c_ctx` are the same constant, not two.** An earlier version
of this design used two independent parameters: `c_gl` for an ordinary
gain/loss step between two nonempty states (e.g. `{A,C}->{A}`), and `c_ctx`
only for a state collapsing all the way to `∅`. That distinction doesn't
survive contact with what the index actually observes. Presence/absence
tracking never sees "the flanks" separately from "the centre" — it sees
whole, distinct, homologous 31-mers (two 31-mers sharing >=30 bases are
homologous by construction, though orthology vs. paralogy is undecidable
from that alone — settled earlier, see "Locus eligibility"). Losing the `A`
form of a family while `C` remains (`{A,C}->{C}`) and losing the last
remaining form (`{A}->∅`) are the *same kind of event*: a specific, complete
31-mer that used to be observed no longer is, because at least one of its
31 positions mutated. There is no separate "gain/loss of a still-recognised
allele" mechanism distinct from "loss of context" — both are exactly the
event `c_ctx(p)` (below) already computes the expected cost of. So: one
constant, `c_ctx`, used everywhere a member is gained or lost — `X -> ∅`
included, and (see next point) not even as a special case there.

**Set-edit-distance formula.** For two states `X, Y ⊆ {A,C,G,T}`, split into
`seulement_X = X \ Y` (size `a`) and `seulement_Y = Y \ X` (size `b`).
Elements present in both cost nothing. Pair up to `min(a,b)` of the
remaining elements as **substitutions** (cheaper than treating them as an
unrelated loss + gain whenever `c_sub < 2*c_ctx`, which any sane parameter
choice satisfies); whatever is left over after pairing is a pure
**gain/loss**:

```
cost(X,Y) = min(a,b)*c_sub + |a-b|*c_ctx
```

Worked examples (`c_sub = c_ctx = 1`): `{A}->{C}` = 1 (one substitution);
`{A}->{A,C}` = 1 (one gain, no substitution pair available since nothing is
only-in-Y that matches an only-in-X element after the shared `A` is
excluded); `{A,C}->{G,T}` = 2 (two substitution pairs, `A/C` vs `G/T`, both
same-size sets share nothing); `{A,C,G}->{A,C,T}` = 1 (`G`/`T` is the only
mismatched pair, `A,C` shared).

**`∅` is *not* a flat-cost special case — a corrected position, reversing
an earlier draft of this design.** An earlier version charged
`cost(X, ∅) = c_ctx` flat, independent of `|X|`, on the grounds that every
member of a family shares identical flanking sequence, so one mutation
breaking that context should unrecognise all of them at once rather than
`|X|` separate times. That argument doesn't survive comparison with how the
rest of the matrix already works: `{A,C,T} -> {A}` (losing two members,
one remaining) is charged `2*c_ctx` via the general formula above, with no
equivalent "maybe it was one shared event" discount — nothing distinguishes
that case from `{A,C} -> ∅` (losing two members, none remaining) other than
which state happens to be the target. Singling out `∅` for special
treatment was arbitrary, not principled: the uncertainty about "one event
or several" is identical in both cases, and the model already resolves it
uniformly elsewhere by simply counting elements. So `∅` is now an ordinary
node in the state graph below like any other, with no patch applied
afterwards:

```
cost(X, ∅) = cost(∅, X) = |X| * c_ctx   (via |X| single-element graph edges)
cost(∅, ∅) = 0
```

**Computing `c_ctx`.** Not guessed — built from the value already derived
for the `>1` bucket in "Context, detectability, and a 3-way ordinal
distance per pair" above:

```
c_ctx(p) = mean_sub_cost * [(2m*p) / (1 - (1-p)^(2m)) + p]
```

`p` is `p_hat`, the calibrated per-site mutation probability (see
`calibrate_p_hat` / the "Experiment" section below); `m = (k-1)/2`, the
flank length on each side of the centre (`k=31` → `m=15`). Read the
bracketed term, well, term by term: `(1-p)^(2m)` is the probability that
*none* of the `2m` flanking positions mutated, so `1-(1-p)^(2m)` is the
probability that *at least one* did — i.e. the probability that this
specific 31-mer stops being observable at all, which is exactly the event
`c_ctx` prices. `2m*p` is the unconditional expected mutation count over
those `2m` positions. Their ratio is the conditional expectation
`E[mutations | at least one occurred]` — provable directly: for any
nonnegative integer random variable `X`, `X · 1{X>=1} = X` always (both
sides are `0` when `X=0`, both are `X` otherwise), so
`E[X | X>=1] = E[X·1{X>=1}]/P(X>=1) = E[X]/P(X>=1)`, which is
`(2m*p)/(1-(1-p)^(2m))` here. The trailing `+p` adds the (much smaller,
first-order, not itself a conditional expectation) marginal contribution
of the centre position's own mutation probability.

**The `mean_sub_cost` multiplier — a correction, not part of the original
derivation.** The bracketed term is a *count* of expected mutations, not a
cost — an earlier version used it directly as `c_ctx`, implicitly pricing
every one of those mutations at a flat `1` regardless of type. That stopped
being defensible once substitution costs were calibrated per base-pair
category (`substitution_costs_from_tally`): transitions are markedly more
frequent than transversions in this project's own data (roughly 5-6x, e.g.
130k vs. 20-28k observed instances per category in the eubacteria run
below) and correspondingly cheaper, so "one mutation" isn't worth a flat
unit. `mean_sub_cost` (`mean_substitution_cost`) is the empirical average
substitution cost, weighted by each category's observed frequency — it
converts the bracketed term's expected mutation *count* into an actual
expected *cost*. Concretely, in the eubacteria run: the bracketed term
alone is `~1.17` (close to `1`, i.e. context loss is *usually* attributable
to a single mutation, since `p_hat` is small); `mean_sub_cost ~= 1.48`
(pulled up from the cheapest transition cost of `1.0` by the substantial
minority of transversions in the mix); final `c_ctx ~= 1.73` — sensibly
between the cheapest transition (`~1.0`) and the transversion costs
(`~2.5`-`2.9`), rather than coinciding almost exactly with the cheapest
transition purely by construction accident, which is what the unmultiplied
version did and which is what exposed this gap in the first place.

**Better construction method: shortest path in a small state graph, not the
closed-form formula directly.** Build a graph on all 16 states, `∅`
included, with two edge types — substitution edges between same-cardinality
sets differing by one element (weight `c_sub`, or `c_ts`/`c_tv` if split
further below), and gain/loss edges between sets whose cardinality differs
by one (weight `c_ctx`, per the unification above; `∅` connects to each
singleton state this way, being a strict subset of it) — then define
`cost(X,Y)` as shortest-path distance in that graph, uniformly for every
`X,Y` including `∅`, precomputed once (16 nodes, trivial) into a dense
16x16 matrix before feeding it to Sankoff/TNT. Verified equivalent to the
closed-form formula above in the
uniform-cost case (checked by hand on `{A}->{C,G,T}`: both give
`c_sub + 2*c_ctx`). The graph construction is not just a reformulation for
its own sake: it is the version that generalises correctly once substitution
costs stop being uniform (next point) — the closed-form's `min(a,b)`
counting silently assumes *any* pairing costs the same, which breaks the
moment `c_sub` depends on which two bases are involved.

**Substitution refinement — implemented as the full 6-category symmetric
matrix, not just Ts/Tv.** `c_sub` was originally going to split into just
`c_ts` (A<->G or C<->T) and `c_tv` (the other four pairs) — already
well-defined in canonical space (see "Canonical invariance" above: a
transition maps to a transition, a transversion to a transversion,
regardless of orientation). A symmetric cost matrix allows finer
resolution than that 2-category split, though: each of the 6 distinct
undirected base pairs (AC, AG, AT, CG, CT, GT) can be costed
independently, with no grouping at all — `SankoffWeights::sub_cost` is a
full `[[f64;4];4]` table (`SankoffWeights::ts_tv` remains as a convenience
constructor for the coarser 2-category case, if ever wanted). This turns
"pick `min(a,b)` substitution pairs" into a genuine (tiny, <=4 elements per
side, trivially enumerable) minimum-cost bipartite matching problem instead
of a plain count — the state-graph shortest-path construction handles this
automatically, no separate logic needed. Calibration is not a guess
either: `KmerIndex::base_pair_tally` collects the pooled `(centre_i,
centre_j)` distribution over resolved (SNP) sites — restricted to the same
ratio-ceiling-included genome pairs as `p_hat`'s own calibration, for the
same saturation-exclusion reason — and `substitution_costs_from_tally`
turns that into `cost(a,b) = -ln(rate(a,b))`, normalised so the most
frequent category costs `1.0` (the standard generalised-parsimony
step-weighting heuristic, generalised from 2 categories to 6). A second
full pass over the annex is required for this (see `base_pair_tally`'s own
docs for why it can't share `raw_snp_distance`'s single pass) — both now
share one traversal helper, `scan_family_pairs`, rather than duplicating
the per-family reconciliation logic.

**Caveat carried over from the `D_F = min(a,b)` rejection earlier:** this
cost matrix is only valid as **Sankoff step-cost input**, re-evaluated for
every branch of every candidate topology during tree search. Reusing
`cost(leaf_A, leaf_B)` directly as a standalone pairwise distance (bypassing
the tree) would reintroduce the exact circularity already rejected — the
`min(a,b)` pairing here is a locally-defined edit distance between two
states, not a claim about the true evolutionary history between two
specific genomes.

**Feasibility confirmed:** TNT's `costs` command accepts custom step
matrices for multistate characters, so this whole construction (16x16
matrix derived from the state graph, `c_ts`/`c_tv`/`c_ctx` as the three
tunable parameters) is directly usable there — no new tooling required
before testing it.

### Experiment: TNT run on real data (2026-08-11)

**Status:** validated at genus/family/order scale within Bacteria; not
informative across domains with this character type. Exploratory only —
run entirely outside the repo (`/tmp/tnt_run`, TNT installed locally under
`TNT/`), no new Rust code. Kept here as the record of what was learned.

**Pipeline.** `obikmer distance --snp` emits one IUPAC-coded pseudo-alignment
row per genome (`snp_pseudo_alignment`, one column per family with
`family_size() >= 2`). A small external Python script decodes IUPAC back to
the 16-state bitmask, applies the set-edit-distance formula above, and
emits a complete TNT script (`xread` matrix + `smatrix` step-matrix + `hold`/
`mult` search). No new Rust code was needed for this pass.

**Calibration ("quick option").** Rather than modifying
`write_raw_snp_distance_csv` to emit raw counts, `c_ctx` was approximated as
the unweighted mean of the pairwise `snp/(snp+shared)` ratios already
present in `rawsnp.csv`, restricted to the relevant taxon subset. Used
values: `mu = gamma = 10`, `c_ctx = 18` (this run predates the `c_gl`/`c_ctx`
unification above — `gamma` here is what `c_gl` was called before it turned
out to just be `c_ctx`; ratio `c_ctx/mu = 1.8`, consistent
with the observed pairwise ratios for genomes within Enterobacteriaceae/
Eubacteria, which cluster around 1.5-2 and barely move as the taxon set
widens — the context signal is stable, not sensitive to which subset is
chosen). The more principled route (raw counts, weighted `p_hat`, Ts/Tv
split) remains a follow-up, not yet done.

**Run 1 — Enterobacteriaceae (11 taxa: 4 *E. coli*, 4 *Salmonella enterica*,
3 *Klebsiella pneumoniae*).** All three genera recovered as monophyletic.
Initial read of the exported (unrooted, TNT/Nexus `[&U]`) tree as showing a
genus-arrangement disagreement with known systematics (*Escherichieae*:
*Escherichia*+*Salmonella* sister vs. more distant *Klebsielleae*) was
**wrong** — diagnosed via `force = (taxa);` monophyly-constraint test
(identical score constrained vs. free ⟹ no real disagreement). Root cause of
the misreading: with exactly 3 clades and no outgroup, an unrooted tree has
only **one possible topology** (a single trifurcation) — there is no
internal arrangement to get right or wrong. This run cannot test
inter-genus relationships at all; it can only test intra-genus monophyly
(which held).

**Run 2 — Eubacteria (18 taxa: run 1 + *Acidobacterium capsulatum*,
*Opitutus terrae*, *Bacillus subtilis*, *Shouchella clausii*, *Wolbachia*
endosymbiont, *Proteus mirabilis*, *Yersinia ruckeri*).** The
Enterobacteriaceae substructure from run 1 is reproduced identically, now
correctly rooted by real outgroups, resolving the tribal arrangement left
undetermined in run 1: *Escherichia*+*Salmonella* sister, *Klebsiella* more
distant — matching known systematics. Outgroup placement: `(Bacillus,
Shouchella)` sister pair (Firmicutes/*Bacillales*) splits from all
Proteobacteria — a correct phylum-level split; `Wolbachia`
(Alphaproteobacteria) splits from the Gammaproteobacteria block
(`Proteus`, `Yersinia`, Enterobacteriaceae) — a correct class-level split.
Lower confidence: the fine nested order `(Proteus, (Yersinia,
Enterobacteriaceae))` and the relative position of `Acidobacterium` vs.
`Opitutus` — plausible, not independently verified against current
Enterobacterales family-level literature.

**Run 3 — full domain set (20 taxa: run 2 + *Candidozyma auris* [yeast] and
*Saccharolobus islandicus* [archaeon]).** The bacterial clade from run 2 is
reproduced **unchanged and intact** — a real robustness signal, the method
does not fragment the ingroup when unrelated deep taxa are added. But the
result carries **no information on Bacteria/Archaea/Eukarya relationships**:
with exactly one eukaryote, one archaeon and one bacterial clade, the
unrooted tree is again forced into the single 3-clade trifurcation from run
1's caveat — there is no second representative of either outgroup domain to
resolve internal arrangement, so nothing about their relative position can
be read from the topology (a "ladder" ordering in the exported tree is
serialization, not signal). Independently, `rawsnp.csv` shows *why* this
character type cannot reach further: pairwise ratios involving
*Candidozyma*/*Saccharolobus* are almost all `NA` (no central-position
family shared at all) or saturated at `1.0` (every shared family differs) —
central-position families require literal 31 bp context conservation, which
simply does not survive domain-level divergence. Cross-domain placement
would need conserved-marker characters (rRNA, ribosomal proteins), not this
estimator.

### Native `--sankoff --tnt`/`--phyg` export (2026-08-12), superseding the external scripts above

The ad hoc Python glue from the previous section is superseded: `obikmer
distance --sankoff` now calibrates the matrix natively (`p_hat`, 6-category
substitution costs, `c_ctx` weighted by `mean_sub_cost` — see the worked
example above) and `--tnt`/`--phyg` each write a ready-to-run script from
it, no external script needed. `∅` is an ordinary 16th state throughout
(never `-`), specifically to avoid gap-semantics confusion in downstream
tools — see "A concrete Sankoff cost matrix" above for why.

**TNT (`--tnt`).** `write_sankoff_tnt` (`obikmer/src/cmd/distance.rs`)
recodes to TNT's own `0-9A-F` xread alphabet (its default reader rejects
the wider IUPAC set otherwise), scales and rounds costs to integers
(`smatrix`/`cost` reject decimals), then re-runs integer Floyd-Warshall on
the rounded matrix (`scaled_metric_matrix`) — independently rounding each
cell of an already-metric real-valued matrix can break the triangle
inequality (e.g. two real costs of `1.734` round to `173` each, summing to
`346`, while their real sum `3.468` rounds to `347`), which TNT otherwise
silently "fixes" itself with an unreproducible correction. Verified against
the real 20-genome benchmark index: zero triangle-inequality violations
after the fix, TNT loads the file without its "triangle inequality
violated... Fixed" warning.

Two syntax facts worth recording because they're wrong in intuitive
guesses and contradicted actual TNT behavior when tested: TNT's plain
command stream has **no comment syntax** of its own — `/* */` and `[ ]`
only work inside the (separately-enabled, off by default) macro scripting
language, and error with "No command!" otherwise. The working substitute
is `quote TEXT ;` (prints the text, doesn't affect parsing) — but the text
itself can't contain a literal `;` (TNT's universal terminator); the
manual's own escape (`.,`) exists but the script here just avoids
semicolons in the text instead.

The default search command embedded in the script is `mult` (traditional:
random addition sequences + TBR), not `xmult` (New Technology search:
ratchet/drift/tree-fusion). `xmult` with TNT's default `mxram` (16 MB,
must be set *before* `xread` if changed) ran out of RAM on the real
908k-character dataset ("`xmult - out of ram`"); `mult` does not, matching
what had already been validated by hand outside this session.

**PhyG (`--phyg`).** `write_sankoff_phyg` writes a `tcm:` custom-alphabet
matrix (same scale+round+metric-closure treatment as TNT) and reuses
`--sankoff`'s own `_sankoff.fasta` as-is via `prefasta:` — PhyG's `tcm:`
alphabet is read from the matrix file's own first line, so (unlike TNT) no
recoding is needed. PhyG auto-adds its own indel/gap state as an
`(n+1)`-th row/column of the tcm; inert here since the alignment encodes
absence as `0`, never `-`.

`report("file", newick, overwrite)` — exactly as shown in PhyG's own
manual — triggers `Unrecognized/missing report option ... defaulting to
'graphs'` on the locally installed binary (1.3, commit `3c1a1fa`); the
working form adds `graphs` explicitly: `report("file", graphs, newick,
overwrite)`. Manual/binary mismatches like this (also true of `criterion:`
— the binary accepts `parsimony`/`ml`/`pmdl`, the manual instead documents
`mapa`/`ncm`/`parsimony`/`pmdl`/`si`) mean command syntax against this PhyG
build should be verified empirically, not trusted from the PDF alone.
`instances:N` (not a separate CPU flag) is what actually parallelises the
search across cores — PhyG uses all physical cores by default but only
across as many instances as are running, so raise it to the physical core
count to use them all (the CLI-level `+RTS -NX -RTS` flag also exists but
controls something else: capping/limiting cores, not requesting more).

Both scripts share one `--sankoff-cost-scale` (default `100`), not two
separate flags — they scale the same calibrated matrix for the same
reason (integer-only cost commands) and no PhyG-specific
accumulator-width constraint was ever found to justify a different
default from TNT's (TNT: hinted 32-bit accumulators in its own manual;
PhyG: no such hint found — Haskell's native `Int` is typically 64-bit).

**Open problem: PhyG reports all branch lengths as `0.0`.** The graph-level
parsimony cost is correct (`3.3286×10⁸` on the real dataset, consistent
with TNT's `328574911` on the same calibrated matrix), but every individual
edge in the exported Newick shows `:0.0`, with the total cost only ever
shown as a whole-tree annotation (`[3.32860377e8]`). Not fixed, not fully
diagnosed — PhyG's manual describes per-edge branch length as computed by
ancestral-state (HTU) backtracking, well documented for sequence/standard
character types, but nothing found (the term "Sankoff" doesn't even appear
in the manual) confirming this backtracking is wired up for a custom
`tcm:` matrix character. Switching `criterion:` to a likelihood-family
option (`ml`, or the manual's `mapa`/`ncm`/`si`) was considered as a
possible fix but is very unlikely to be one: those criteria are
information-theoretic reparametrisations of the *same* step-counting
machinery as parsimony (`ncm` in particular is known in the literature to
be numerically equivalent to weighted parsimony), not classical
continuous-time-Markov ML with a real rate matrix — so they wouldn't
change how branch length is attributed per edge either.

**Export format note (not a bug in the generator).** Neither script's
`.tre` output opens in PearTree (FigTree's successor) via File > Open —
association/Launch-Services quirks were ruled out (the file was opened
directly through the app, not by double-click). Likely cause, not yet
confirmed: TNT's export is a minimal NEXUS `begin trees;` block with no
preceding `Taxa` block and bare numeric (untranslated) leaf labels; PhyG's
is multiple raw Newick trees concatenated with no NEXUS wrapper at all
plus a trailing `[cost]` bracket tag after the root label. Both differ
from a "normal" single-tree, fully-declared NEXUS file; this is PhyG/TNT's
own export format, not something `write_sankoff_tnt`/`write_sankoff_phyg`
could fix without post-processing the *other* program's output after the
fact.

### Next direction: genuine ML branch lengths, not parsimony (open, 2026-08-12)

Decided: parsimony (the whole `--sankoff`/`--tnt`/`--phyg` pipeline above)
is a stopgap, not the destination. The goal is maximum likelihood with
real, calibrated branch lengths (expected substitutions/site), which
parsimony step-counts were never going to give directly (see the open
"branch lengths are `0.0`" problem above — even if fixed, TNT/PhyG-style
parsimony branch length is a step count, not a continuous ML estimate).

**Model choices, settled:**
- **The exchangeability `R` is symmetric; the rate matrix `Q` is not.**
  (Superseded an earlier, wrong framing here that treated "symmetric
  model" as one thing — see the resolution below on `R` vs `Q` vs `π` for
  the full reasoning.) `R(a,b) = R(b,a)` because `BasePairTally` never
  captured direction — a fact about the data, not a modelling choice.
  `Q(a,b) = R(a,b)·π_b` is asymmetric whenever the real state frequencies
  `π` are (which they are, empirically) — biology drives this via `π`,
  not via `R`.
- **`∅` stays an ordinary 16th state**, as already established for
  TNT/PhyG — same reasoning applies to any ML tool: encode as a real
  alphabet symbol, never as `-`/gap, or the RAxML-era failure (empty set
  silently treated as missing data) repeats.

**Stationary frequencies for the 16 states — resolved (2026-08-12).**
A CTMC needs a rate matrix `Q`, generally asymmetric. `Q(i,j) = R(i,j) ·
π_j`, where `R` (exchangeability) is symmetric and `π` (stationary
frequencies) need not be — this factoring is what makes `Q` reversible
(satisfies detailed balance, `π_i·Q(i,j) = π_j·Q(j,i)`) for *any* `π`, not
just uniform, as long as `R` is symmetric. Two separate, both-easy
quantities, not one hard inverse problem:
- **`R` is already calibrated**: `sub_cost` (`-ln(observed rate)` per pair,
  from `BasePairTally`) *is* `R` up to a log transform — recover it as
  `R(a,b) = exp(-sub_cost(a,b))`. Symmetric by construction, because the
  tally itself never distinguished direction (unordered-pair counts only)
  — not a modelling choice, a fact about what the data can say.
- **`π` is a direct count**: empirical marginal frequency of each of the
  16 states across the whole alignment (same kind of scan already used to
  confirm `N` occurs 1383 times in the real 20-genome benchmark). With
  ~908k sites × 20 genomes, the counts are large enough that this is
  precise on its own — no need to spend ML degrees of freedom
  re-estimating it via IQ-TREE's `+FO`.

Checked and ruled out along the way: IQ-TREE's `+F` (empirical, "compute
from the alignment") does **not** work as a shortcut for this — for a
custom-file morphology model, `readParameters` always requires the file's
own frequency line unconditionally; omitting it and passing `+F` instead
just fails (`ERROR: State frequencies could not be read`). `π` has to be
computed by `obikmer` and written into the file, not left to IQ-TREE.

Net effect: no free-rate ML estimation needed for this piece at all (the
mistaken assumption that motivated most of this discussion — that
building an asymmetric `Q` a priori would require solving a linear system
from `Q` itself — doesn't apply, because `R`, the only piece that's
genuinely hard to get directionally, is symmetric and already in hand).

**Candidate tool: IQ-TREE**, because it supports user-defined multistate
models (unlike RAxML's `MULTI` data type, which is limited to the
equal-rate Mk model and can't take a custom rate matrix at all — a genuine
tool limitation, not a gap-symbol encoding problem this time). IQ-TREE 3
(3.0.1) is now installed locally (Homebrew, `iqtree3`).

**IQ-TREE custom-model format — verified empirically against the local
binary (2026-08-12).** The web docs' `-mdef` NEXUS `begin models;
frequency NAME = ...; model NAME = ...; end;` mechanism (initially assumed
to apply directly, see history below) turned out to be for **named
components used inside `MIX{...}`/`FMIX{...}` mixture models only** — it
does **not** apply to a single, non-mixture custom morphology matrix, and
using it that way fails (`ERROR: File not found <name>` — traced in
IQ-TREE 3's own source, `model/modelmorphology.cpp`: any `-m` string that
isn't `MK`/`ORDERED`/`GTR`/`GTRX` is passed straight to
`ModelMarkov::readParameters()`, which opens it **as a literal file path**,
never consulting the `-mdef` models block at all for this data type).

**The confirmed working recipe** (built a tiny 5-taxon/3-state toy dataset
and rate file, ran it end to end with `iqtree3`, got a real ML tree with
non-zero branch lengths and an optimized log-likelihood — ground truth,
not documentation):
- No `-mdef` needed. Write one plain file (any name) containing, as
  whitespace/newline-separated numbers, in order: the **lower-triangular
  rate matrix** (`N(N-1)/2` values, PAML row-major order — for 16 states,
  120 values, the same count and layout already produced for TNT's
  `smatrix`), immediately followed by the **N state frequencies** on the
  same stream (no header, no blank line required — confirmed by reading
  `ModelMorphology::readRates`/`ModelMarkov::readStateFreq` directly, which
  just pull tokens off the stream in sequence).
- Invoke with `-m <path-to-that-file>+ASC` (`+ASC` for the no-constant-site
  correction, as before). An explicit `+F{f1,...,fN}` on the command line
  overrides the file's own frequency line if given (confirmed in
  `ModelMorphology::init`) — useful once real calibrated stationary
  frequencies exist, a placeholder equal-frequency line works meanwhile
  (the still-open gap noted above).
- `--seqtype MORPH` (alphabet `0`-`9`,`A`-`V`, ≤31 states) — reuse the same
  `0-9A-F` recoding already built for TNT (`TNT_STATE_SYMBOL`).

**Risk, confirmed, precisely characterised, and resolved by design
(2026-08-12).**
`--seqtype MORPH{16}` **does not force the state count** for real ML
analysis — tested directly (`--seqtype MORPH{4}` on the 3-symbol toy
alignment gave the byte-for-byte identical 3-state result as no `{4}` at
all) and confirmed in source: the value it sets
(`params.alisim_num_states_morph`, `utils/tools.cpp`) is consumed only by
the `--alisim` simulator; the main analysis path always calls
`getDataBlockMorphStates`/an equivalent scan (`alignment/alignment.cpp`),
for both FASTA/PHYLIP and NEXUS input (a NEXUS `symbols=` declaration
doesn't change this either — checked, same code path). No CLI flag or
NEXUS declaration overrides it.

The precise rule (from `getDataBlockMorphStates`, `alignment.cpp:1058`):
`N` = **one plus the highest state ordinal actually observed anywhere in
the alignment**, ordinal being the symbol's position in IQ-TREE's own
fixed table `"0123456789ABCDEFGHIJKLMNOPQRSTUV"` — not a count of distinct
symbols seen. So the risk is narrower than "any missing symbol breaks it":
concretely, it's whether the symbol mapped to state index 15 (`F` in the
`0-9A-F` recoding already used for TNT, i.e. our `N` = "all four bases
ambiguous") occurs **at least once anywhere** in the real alignment — if
it does, `N` correctly comes out to 16 regardless of which lower-index
symbols (including `0`/`∅`) are rare or absent; if it doesn't, `N` silently
undercounts and misaligns every value in a 16-entry rate/frequency file,
with no error to catch it. **Checked against the real biological alignment**
(`/tmp/msg_test/eub_sankoff.fasta`, the 20-genome benchmark index): the
symbol `N` (IUPAC "all four bases ambiguous," state index 15) occurs 1383
times across 13 of the 20 sequences — present, so this specific real
dataset is not at risk. Still worth a real, general presence check inside
`obikmer` before this is wired in, rather than assuming every future
dataset will have `N` too (nothing in IQ-TREE would catch it if not).

**Resolution: subset + compact-renumber, not rely on all 16 appearing.**
Since IQ-TREE always infers `N` from the alignment's own content and
nothing overrides that, the fix is to make the file `obikmer` writes match
that inference *by construction*, for every run, rather than hope the
16th (or any particular) state happens to occur:
1. Scan the real alignment for which of the 16 canonical states actually
   occur anywhere (not per-column — anywhere in the whole alignment).
2. Renumber the occurring states to a **compact, consecutive** `0..k-1`
   range, preserving their relative order (the original bitmask/`STATE_
   SYMBOL` ordering) — not just filtering, since a *gap* in the ordinal
   sequence (e.g. keeping states `{0,1,2,4}` numbered as-is instead of
   `{0,1,2,3}`) reproduces the exact same "highest observed ordinal"
   miscount this was meant to fix.
3. Recode the alignment itself with this new compact `k`-symbol alphabet
   (same recoding mechanism already used for TNT's `0-9A-F`, just over a
   possibly-smaller symbol set).
4. Extract the matching `k×k` submatrix (rows/columns for the kept
   states only) from the full calibrated 16×16 cost matrix, in the same
   lower-triangular order the rate-matrix file needs — and, later, the
   matching `k`-length subset of stationary frequencies once those are
   calibrated (still the open gap noted earlier in this section).

Consequence, and why nothing is lost: a state that never occurs in a given
alignment can, by definition, never contribute a transition to score in
that same alignment — dropping it from that run's matrix costs nothing.
The subset (and therefore `k`) can differ from one dataset/run to the
next; this has to be done freshly per alignment, not computed once and
reused.

*(Superseded reasoning, kept for the record: the `-mdef` NEXUS route
below was the original plan, based on IQ-TREE's own web documentation for
protein mixture models, before the local install allowed testing it —
`GTRX` combined with a `-mdef`-referenced custom model, `+Fname` frequency
reference. Both pieces exist and parse without error individually, but
`GTRX`/`GTR` are IQ-TREE's own fixed built-in equal-structure multistate
model, not a hook for an arbitrary custom matrix; a custom matrix is a
file path in `-m` directly, no `-mdef` or `GTRX` involved.)*

Source for the empirical findings above: `model/modelmorphology.cpp` and
`model/modelmarkov.cpp` in the local `iqtree/iqtree3` source (cloned to
inspect the exact parsing logic after documentation didn't resolve the
`+Fname` reference error) — more reliable here than the PDF/web manual,
which (like TNT/PhyG) doesn't always match this specific binary. Original
(partially superseded) sources: [Substitution
Models](https://iqtree.github.io/doc/Substitution-Models), [Complex
Models](https://iqtree.github.io/doc/Complex-Models).

**Relation to the existing calibration.** `sub_cost[a][b] = -ln(rate)`
(see "A concrete Sankoff cost matrix" above) is already a log-rate — a
genuine CTMC rate matrix `Q` could plausibly be recovered as
`rate(a,b) = exp(-cost(a,b))`, renormalised so each row sums to zero, once
the stationary-frequency gap above is closed. Superseded by the two
sections below, which implement and then substantially revise this.

### `R`/`π` implemented; `--exclude-genome` added; rogue-taxon test negative (2026-08-12)

`write_iqtree` (`obikmer/src/cmd/distance/iqtree.rs`) implements exactly
the `R = exp(-cost)` / empirical-`π` design above: writes
`<prefix>_iqtree.model` (lower-triangular `R`, PAML order, then `π`) and
`<prefix>_iqtree.fasta` (alignment recoded to the compact `0..k-1`
alphabet of states actually present), states/frequencies restricted to
whichever of the 16 canonical states actually occur, compactly
renumbered — see the `MORPH{N}` risk above. Verified end to end on the
20-genome benchmark: real, non-zero, varied branch lengths (`Total tree
length: 6.773`), converged log-likelihood, `State frequencies:` in
IQ-TREE's own output matching the computed `π` exactly.

**Rogue-taxon hypothesis, tested and refuted.** The backbone (*Yersinia*,
*Proteus*, *Opitutus*, *Shouchella*, *Wolbachia*...) resolves as a
near-linear comb with ~8 near-zero branch lengths — visible on the real
tree, and independently on a Mash+NJ tree built from an entirely
different signal (whole-genome k-mer distance, no relation to the Sankoff
pipeline), which shows the same comb shape. Hypothesis: `Saccharolobus`/
`Candidozyma` (near-zero real signal — see below) destabilise the
heuristic tree search enough to also blur resolution elsewhere ("rogue
taxa", a documented phenomenon). Tested directly: reran `iqtree3` after
removing both taxa (`--exclude-genome`, `+ASC` recomputed on the
resulting variable-sites-only alignment since removing taxa turns some
columns invariant) — **still exactly 8 near-zero backbone branches**,
identical comb shape. Refuted for this dataset: the backbone's weak
resolution is a property of the character system's signal at that
divergence depth (matches "Run 3" above), not rogue-taxon interference.

**Diagnosed why `Saccharolobus`/`Candidozyma` place so poorly**, using a
new diagnostic (`--raw-snp-counts`, `<prefix>_rawsnp_counts.csv`: one row
per genome pair, `n_snp,n_shared,n_eligible,ratio` — the counts
`--raw-snp-distance`'s ratio-only matrix discards, needed because
`ratio=0.0` from 2 eligible loci and from 2000 look identical in the
ratio alone). `Saccharolobus` has 21,020 non-`∅` sites in the real
alignment (comparable to other taxa) — but **21,017 of them (100.0%) are
private**: no other genome has a non-`∅` state at the same site. Only 3
sites are shared with any other genome at all (1, 2, and 14
co-occurring genomes respectively). Real informativeness for placement
tracks shared sites, not raw non-`∅` count — with no other archaeon (or
even archaea-adjacent bacterium) in the dataset, there is nothing to
anchor `Saccharolobus`'s position to, regardless of how much of its own
data exists. `Candidozyma` shows the same pattern, more extreme (almost
all `NA`, the 3 non-`NA` pairs all exactly `0.0`, never `1.0` — itself a
tell: with `n_eligible=1`, the ratio can only be exactly `0` or `1`, so
3-for-3 landing on `0.0` is more than sampling noise alone would predict;
possibly ascertainment bias — see below — or possibly a few genuinely
ultra-conserved loci; not resolved).

**`--exclude-genome LABEL`** (repeatable, `obikmer distance`) added for
exactly this kind of test: zeroes the excluded genome's row/column in
`RawSnpDistanceOutput` after `raw_snp_distance` runs (a pair with zero
counts is already skipped by `calibrate_p_hat`/`base_pair_tally` — no
`obikindex` traversal change needed) and drops its row from
`SnpAlignment` before any output is written. Deliberately *not* index
surgery (a new, smaller on-disk index) — genome sets to exclude are
expected to change between quick tests, so an in-memory filter is the
right tool, not a new index-rewriting subsystem. Scoped to the sibling-annex
family of computations (`--raw-snp-distance`/`--raw-snp-counts`/`--snp`/
`--sankoff` and everything it implies) — does not affect the plain
`--metric` distance matrix/NJ/UPGMA path (a different, unrelated
computation on `idx.distance()`, not touched).

**Caveat surfaced while reusing a stale `π`**: rerunning IQ-TREE on a
genome-reduced alignment while keeping the *original* (20-genome)
model file is inconsistent — `π`'s composition shifts once low-cardinality-0
columns that were only variable because of the removed taxa drop out.
Measured directly: `π(∅)` `0.9048` (20 genomes) → `0.8869` (18 genomes,
variable sites only) — real (~378k affected cells, from `Saccharolobus`'s
~21k private sites × 18 remaining genomes) but modest in *proportion*
(~2 points) because the alignment was already monomorphic-filtered
before `Saccharolobus` was ever added, so removing it only drops the
subset of columns that were variable *because of* it specifically, not
every column it appears in.

### `R` via `exp(-cost)` is wrong for a CTMC; cardinality/composition decomposition (open, 2026-08-12)

**The flaw in `R = exp(-cost)`, precisely.** `cost` (`build_cost_matrix`'s
output) is a *shortest-path closure* over an elementary-edit graph
(Floyd-Warshall) — correct and required for Sankoff parsimony, where
`cost(a,b)` must be a metric. But a CTMC's own matrix exponential
(`exp(Qt)`, computed internally by IQ-TREE) *already* sums over
paths of every length through the elementary rates — that's the whole
mechanism by which a CTMC generates indirect transitions. Feeding it a
pre-summed, multi-hop shortest-path cost and exponentiating that *again*
as if each entry were a direct edge double-applies the "compose multiple
steps" logic once in log-space (Floyd-Warshall, additive) and once more
inside IQ-TREE's own exponential — systematically over-penalising
non-adjacent state pairs (e.g. `∅→{A,C}` priced as two chained edges,
`2×c_ctx`, when it should be one direct lookup).

**Resolution, in two parts — both estimated directly from the real
alignment, not smoothed through a small parametric formula:**

1. **Cardinality model**: a 5-state (`0,1,2,3,4`) first-order Markov
   chain, estimated from the empirical cardinality co-occurrence table
   (pooled across all *included* genome pairs — same
   saturated/`NA`-pair exclusion discipline as `calibrate_p_hat`, 14
   saturated + 18 `NA` pairs excluded of 190 in the benchmark; the
   result is materially different from the unfiltered version and more
   internally consistent, not just "cleaner"). Diagonal included (the
   probability of a family *staying* at the same cardinality is part of
   the model, not assumed away). Measured, cost `= -ln(observed/expected
   under independence)`, on the 20-genome benchmark (`π` here from the
   marginal cardinality distribution: `90.481%, 5.261%, 4.183%, 0.067%,
   0.008%` for `c=0..4` respectively — note `c=2` is *not* rare, almost
   as common as `c=1`):

   ```
   c=0/c=0: 0.166   c=0/c=1: 0.160   c=0/c=2: 0.157   c=0/c=3: 0.059   c=0/c=4: 0.006
   c=1/c=1: -0.803 (enriched)        c=1/c=2: 0.381   c=1/c=3: 1.027   c=1/c=4: 1.794
   c=2/c=2: 4.032 (sharply suppressed)                c=2/c=3: 2.620   c=2/c=4: 2.416
   c=3/c=3: -0.254  c=3/c=4: -0.411 (too few observations to trust)
   ```
   Two robust findings, both stable under the saturation filter: (a) all
   `∅`-involving costs are low and close to each other (`0.006`–`0.166`)
   regardless of how many members are gained/lost at once — sharply at
   odds with the current model's implicit `×2`/`×3` multi-hop scaling;
   (b) `c=2/c=2` (two genomes both showing an ambiguous 2-member state at
   the same site) is dramatically under-represented (~1.5–1.8% of the
   independence expectation) — a real, robust anomaly, not explained.

2. **Composition model**: unchanged — `sub_cost` (the existing 6-category
   Ts/Tv-biased calibration), estimated **only from unambiguous sites**
   (cardinality-1 ↔ cardinality-1 pairs), because composition bias is a
   substitution phenomenon and only means something when cardinality is
   conserved. Checked whether composition bias also appears in pure
   gain/loss events (no substitution involved, so no bias expected a
   priori): single-base "gain" (`∅→`single base) is close to uniform
   (`24.35/25.92/26.22/23.51%` for A/C/G/T) — consistent with "no
   mutational mechanism, so no bias" as expected. Two-base "gain"
   (`∅→`2-member state) is *not* uniform even after correcting for the
   real (non-25/25/25/25) marginal base frequencies: `{A,G}`/`{C,T}`
   (the transition-linked pairs) mildly enriched (`obs/exp` `1.12`/
   `1.16`), `{C,G}` sharply suppressed (`obs/exp 0.563`, expected to be
   the *most* common pair under independence since C and G are
   individually the two most frequent bases, observed the least) — an
   unexplained anomaly, deliberately **not** built into the model (no
   mechanistic story for why gain/loss would carry a `{C,G}`-specific
   bias), left as an open puzzle rather than fit.

**Composing the two into a full pairwise cost — the actual replacement
for both the elementary-edge graph and its Floyd-Warshall closure.** For
any two states `A`, `B` (not just the "clean" same-cardinality or
pure-subset cases the current graph handles directly): let `shared = A∩B`
(free), `lost = A\B`, `gained = B\A`. Pair off `k = min(|lost|,|gained|)`
elements between `lost` and `gained` as substitution events, choosing the
pairing that minimises total `sub_cost` (a trivial assignment problem —
at most 4 elements per side). The `|lost|-k` (or `|gained|-k`) leftover,
unpaired elements are a *pure* cardinality change, `|A|→|B|`, costed by
**one direct lookup** in the cardinality model above — no chaining, no
Floyd-Warshall. Example: `{A,C}→{G}` (cardinality 2→1, no shared base):
pair 1 substitution (cheaper of `A→G`, `C→G`), 1 base left over unpaired
→ cost `= sub_cost(chosen pair) + cardinality_cost(2→1)`. This is a
direct, closed-form cost for *every* pair of the 16 states, replacing
`build_cost_matrix`'s graph-plus-shortest-path construction outright —
and specifically fixes the CTMC double-counting problem, since every
entry is now a single decomposed lookup, never a sum of chained edges.

**Reframed as a likelihood (product of probabilities), not a cost (sum of
`-ln`s) — same content, but forces the diagonals (no-change cases) to be
kept rather than implicitly dropped.** Both sub-models are proper
transition *probability* matrices, diagonal included — `P_cardinality`
includes "stay at the same cardinality", `P_composition` includes "stay
the same base" (e.g. `P(G→G)`, not just `sub_cost`'s off-diagonal
entries). Composing:

```
P(A→B) = P_cardinality(|A|→|B|)
        × ∏_{x ∈ A∩B} P_composition(x→x)         (shared bases: "stayed")
        × ∏_{(x,y) paired}  P_composition(x→y)    (parsimony-paired substitutions)
```

with the unpaired leftover `lost`/`gained` elements (if `|lost|≠|gained|`)
contributing *nothing further* beyond the `P_cardinality` term already
counted — consistent with the finding above that pure gain/loss carries
no separate composition bias worth modelling. Worked example,
`{A,C}→{A,G}` (`shared={A}`, one paired substitution `C→G`, nothing left
over): `P = P_cardinality(2→2) · P_composition(A→A) · P_composition(C→G)`.

The row-wise product of these two independently-calibrated models isn't
guaranteed to already sum to exactly `1` across all `B` for a fixed `A`
(the two aren't perfectly independent in reality) — so each row of the
resulting 16×16 matrix is renormalised (divided by its own sum) after
composition, **not** the matrix as a whole (which would produce a joint
distribution over `(A,B)` pairs, the wrong object — a transition matrix
needs each row, "given I start in `A`", to be a valid distribution over
where I end up).

This now covers every pair of the 16 states with no unhandled case
identified. Gives 120 parameters, but derived from two small,
well-estimated pieces (a 5×5 cardinality model, a 4×4 composition model)
rather than fit or smoothed independently per pair.

### Implemented (2026-08-12): `pairwise_cost_matrix` replaces `build_cost_matrix` entirely

New module `obikindex/src/cardcomp.rs`, replacing `sankoff::build_cost_matrix`
and the `c_ctx`/`SankoffWeights`/`PHatEstimate`/`calibrate_p_hat`/
`c_ctx_from_p_hat`/`mean_substitution_cost`/`substitution_costs_from_tally`
machinery it depended on outright — not kept in parallel as a fallback
(all now unreferenced outside their own tests; `sankoff.rs` itself is a
pending removal, not yet done).

**New primitives, `obikindex/src/siblings.rs`:**
- `BasePairTally` gained a `same: [u64; 4]` field (diagonal — "both
  genomes at the same single base", pooled from
  [`base_pair_tally`](obikindex::KmerIndex::base_pair_tally)'s existing
  traversal, extended to also tally the `bi == bj` case it previously
  discarded).
- `CardinalityTally { counts: [[u64; 5]; 5] }` and
  `KmerIndex::cardinality_tally`, a new traversal (same shape as
  `snp_pseudo_alignment`'s — needs full per-genome presence masks, not
  `scan_family_pairs`'s single-resolved-form view, since cardinality 2-4
  is exactly the signal being tallied, not noise to drop). Same
  saturated/no-data pair exclusion as `base_pair_tally`. Restricted to
  variable families (`family_size() >= 2`), matching
  `snp_pseudo_alignment`'s own scope. Verified against the 20-genome
  benchmark: counts match the earlier hand-rolled Python analysis exactly
  (e.g. `c=0/c=0: 117,158,166`, `c=0/c=1: 13,707,223` — the same numbers
  this whole investigation started from).

**`cardcomp.rs`:**
- `cardinality_transition_probs`/`composition_transition_probs`: row-
  normalise the two tallies into proper transition probability matrices,
  diagonal included ("stay the same" is a real, calibrated outcome).
- `pairwise_cost_matrix`: for every pair of the 16 states, `shared = A∩B`
  contributes `∏ P_composition(x→x)`; `lost = A\B`, `gained = B\A` are
  parsimony-paired (`best_pairing_cost`, brute-force over the ≤4!
  injections — small enough that hand-rolling beats a dependency) into
  substitution events on `P_composition`, minimising total `-ln`; the
  cardinality-difference leftover is priced once via
  `P_cardinality(|A|→|B|)`, never chained. Row-normalised, `-ln`'d, then
  **symmetrised**: `cost_sym(A,B) = (cost(A,B)+cost(B,A))/2` — equivalent
  to taking the *geometric* mean of the two raw probabilities
  (`-ln(√(P(A,B)·P(B,A))) = (-ln P(A,B) - ln P(B,A))/2`), not their
  arithmetic mean. Required, not just convenient for IQ-TREE's
  lower-triangular file format: Sankoff parsimony's score is independent
  of where an *unrooted* tree (what TNT/PhyG actually search over) gets
  rooted only if the cost matrix is symmetric — the discrete-parsimony
  analogue of CTMC reversibility, established by direct reasoning, not
  assumed. Bonus of the same decision: 120 free parameters instead of the
  240 a fully asymmetric matrix would need.

**Verified on the 20-genome benchmark**: resulting matrix symmetric
(checked numerically, zero asymmetric cells), zero diagonal, no NaN/Inf.
`--tnt` output still loads into TNT with no triangle-inequality warning
(`scaled_metric_matrix`'s rounding-metric-closure step still needed and
still applied — nothing in the new construction guarantees the *rounded
integer* matrix stays a metric, even though the real-valued one is exact
by construction here, unlike the old Floyd-Warshall-closed matrix which
needed it for a different reason). IQ-TREE loads the new model file and
reports the same `π` as before (only `R` changed).

### Two consistency bugs found and fixed post-implementation (2026-08-13)

**`--exclude-genome` didn't drop columns that become monomorphic once the
excluded genome(s) are gone.** `snp_pseudo_alignment`'s "variable family"
test (`family_size() >= 2`) is a property of the annex computed over
*every* genome in the index — unaffected by the CLI-level exclusion, which
only dropped the excluded genome's *row*. A family variable only because
of the excluded genome stayed in the alignment as a now-constant column —
silently wrong data for TNT/PhyG, and a hard failure for IQ-TREE's `+ASC`
(verified: excluding 2 taxa on the benchmark left 116,351 such columns —
matches the manual `+ASC` failures hit earlier in this same investigation,
before `--exclude-genome` existed). Fixed in `drop_excluded`
(`obikmer/src/cmd/distance/mod.rs`): after dropping excluded rows,
re-scan each column among the *surviving* sequences and drop any that are
now constant. Verified: 908,723 → 792,372 sites after excluding 2 taxa,
zero monomorphic columns remain, `π` recomputed from the corrected
alignment matches an independent recount exactly. The compact-alphabet
renumbering (`iqtree::compact_alphabet`) needed no equivalent fix — it
already recomputes which of the 16 states occur fresh on every call, from
whatever alignment it's actually handed, so a symbol disappearing (e.g.
excluding every genome that carries `N`) is already handled correctly;
verified directly (excluded 13 genomes to force `N` out: "15 of 16 states"
reported, correctly-shaped model file).

**`cardinality_tally`'s `family_size() >= 2` filter looked inconsistent
with `base_pair_tally` — removing it was tried, and was wrong; reverted.**
`cardinality_tally` (modelled after `snp_pseudo_alignment`) had the
filter; `base_pair_tally` didn't (it visits every family via
`scan_family_pairs` unconditionally, folding fully-invariant loci into its
own `same` diagonal). Read as `cardinality_tally` under-counting its
diagonal relative to `base_pair_tally`, and — independently — as another
angle on the `--exclude-genome` drift (`family_size()` being global-only
meant a family kept here post-exclusion could differ from what the
now-correctly-filtered alignment kept). First fix tried: drop
`cardinality_tally`'s filter entirely, matching `base_pair_tally`'s
whole-annex scope.

**That fix was empirically wrong, confirmed by a real IQ-TREE run, not
just a hunch.** Log-likelihood dropped from the earlier correct run's
`-8,364,671`/`-8,371,082` to `-9,170,228` (worse fit, not better), with
repeated `NNI search needs unusual large number of steps (20) to
converge!` warnings — and the completed run's **total tree length came
out at 67.644**, roughly 30× the earlier correct runs' ~2.0, i.e. branches
blowing up/saturating. Root cause, only clear in hindsight: `+ASC`
("ascertainment bias correction") exists specifically because the
likelihood only ever sees *variable* sites — the alignment fed to
IQ-TREE, by construction, contains not one invariant column. Calibrating
`R` from a population overwhelmingly dominated by genome-wide invariant
background (family_size()<2 loci outnumber the ~908k variable ones by
orders of magnitude) describes a completely different population than the
one `+ASC` and the alignment actually model — the "consistency" argument
for matching `base_pair_tally`'s scope was real, but pointed the wrong
way: `base_pair_tally`'s own unrestricted `same` diagonal turned out to
have the *identical* latent bug (only unmasked once its diagonal existed
at all, which happened earlier the same day when `same` was added), not a
correct baseline to match `cardinality_tally` to.

**Final fix**: restored `cardinality_tally`'s `family_size() >= 2` filter,
and gave `base_pair_tally`'s `same` diagonal the equivalent restriction —
`scan_family_pairs` (shared with `raw_snp_distance`, which legitimately
*does* want fully-invariant families counted as `shared`) now passes an
extra `variable: bool` (the family's own `family_size() >= 2`) to its
`on_pair` callback; `base_pair_tally` only increments `same` when
`variable` is true, `raw_snp_distance`'s callback ignores the new
argument. Both tallies now describe the same variable-families-only
population the `+ASC`-corrected alignment does. Verified: calibration
counts back to their original values exactly (`c=0/c=0`: 117,158,166,
matching the pre-regression run bit for bit), and a full IQ-TREE rerun
converged normally — log-likelihood `-8,389,106.273` (same order as the
two earlier correct runs), **total tree length 2.044** (was 67.644), no
NNI convergence warnings.

## Heterozygosity, ploidy, and consensus-assembly inputs

A within-genome multiplicity signal (more than one of the 4 central forms
present at a locus) is produced identically by two distinct causes:
paralogous duplication and diploid/polyploid heterozygosity. K-mer data alone
cannot distinguish them. The one real discriminator is sequencing depth
(heterozygous site: total depth of the present forms ~= the genome's
single-copy average; duplication: ~2x or more) — but that signal only exists
if genome "counts" are raw-read depth (FASTQ input), not occurrence counts in
an assembled FASTA, where per-locus depth is not preserved.

**Magnitude is taxon- and mating-system-dependent, not universal.**
Heterozygosity density: mammals ~1 site / 1-1.5 kb (~0.1%); highly
outcrossing plants (maize, poplar) reported an order of magnitude higher
(~1%); self-fertilising plants (*Arabidopsis thaliana*) near zero — but with
a documented failure mode where segmental duplication masquerades as
"pseudo-heterozygosity"; fungi split between haploid vegetative stages
(non-issue) and dikaryotic Basidiomycetes, where two long-diverged haploid
nuclei coexist without fusing. The estimator's target use case (closely
related genomes, k=31) is exactly where the stringent filter above costs the
least for low-heterozygosity taxa and the most for outcrossing/dikaryotic
ones — no universal threshold; this is a scope caveat to document, not a
problem to solve generically.

**Why assembled-consensus inputs don't make measured distances wrong.**
Phylogenetic inputs are near-universally assemblies, not raw reads, and
assemblers collapse heterozygous sites to one consensus allele per
position — effectively an arbitrary, largely uncorrelated-between-assemblies
choice at each het site. This does not inject unbounded noise: standard
population genetics gives `d_xy = d_a + (pi_A + pi_B)/2` — the expected
pairwise difference between a random allele of population A and a random
allele of population B equals the net (fixed) divergence `d_a` plus the
average of the two populations' own within-population diversity `pi`.
Consensus flattening realises exactly this random-allele draw, so the
measured genome-to-genome distance is a `d_xy`-like quantity, not `d_a` —
inflated by heterozygosity by a well-characterised additive term, not
distorted unpredictably. The term is negligible when `pi << d_xy` (the common
case for cross-species comparisons), and becomes material precisely in the
two cases already flagged above: very closely related genomes (this
estimator's explicit target) and highly heterozygous outcrossing organisms,
where `pi` and `d_xy` are the same order of magnitude.

Caveat: this assumes the flattening is uncorrelated with the phylogenetic
signal — plausible for de novo assembly, not guaranteed for reference-guided
assembly biased toward one allele (e.g. the reference's) at each het site,
which would turn the noise term into a systematic bias toward the reference
lineage. Not evaluated here.

**Forward-looking implication, not part of the current design.** The
multiplicity > 1 signal discarded by the stringent filter is a crude
per-genome proxy for `pi` (under low background paralogy). If a `pi_hat` per
genome were tallied alongside `SnpTally`, a `d_a` correction
(`p_hat - mean(pi_hat_i, pi_hat_j)/2`, roughly) could recover an estimate
closer to net divergence instead of `d_xy` — a possible extension, not
scoped here.

## Sufficient statistic: 4x4 base-pair tally

Tabulating the joint distribution of `(center_i, center_j)` over conserved-flank
loci, per genome pair, is sufficient for every downstream correction:

| Estimator | Input | Formula |
|---|---|---|
| Raw p-distance | total off-diagonal / total | `p = SNP / (SNP + shared)` |
| Jukes-Cantor | p | `d = -3/4 * ln(1 - 4p/3)` |
| Kimura 2-parameter | transition rate P, transversion rate Q | `d = 1/2 ln(1/(1-2P-Q)) + 1/4 ln(1/(1-2Q))` |
| LogDet/paralinear | full 4x4 + base-composition margins | `d ~= -1/4 ln det(F)`, robust to non-stationary base composition |

JC/K2P need only the total and the transition/transversion split (the
diagonal collapses to a single "shared" total). LogDet needs the full 4x4,
already populated at no extra cost (see Step 1/2 below).

Memory for the 4x4 tally: `n^2 * 16` counters. Trivial for the project's
genome-scale use case (tens to hundreds of genomes); ~13 GB at n=10^4 — outside
scope but worth flagging if n grows.

## Biases (properties of the estimator, not defects)

1. **Conserved-flank ascertainment bias.** Only SNPs with intact `2m`-base
   flanks are visible; window-intact probability decays as `(1-p)^{2m}`. For
   k=31 (2m=30): 0.74 at p=1%, 0.21 at p=5%, 0.04 at p=10%. This estimator
   targets **closely related genomes**. Under rate heterogeneity across sites
   (universal in practice), conserved flanks correlate with slow centers, so
   `p_hat` underestimates the genome-wide average rate — it specifically
   estimates the substitution rate of **conserved regions**.
   Two distinct factors are at play here, not one: `P(centre of a given
   window is a SNP) = p` exactly, **independent of k** — a direct restatement
   of the raw per-site rate via the bijective window<->centre-position
   correspondence (Statistic section above), not a k-dependent quantity.
   `(1-p)^{2m}` is the *separate*, genuinely k-dependent ascertainment factor
   (are the flanks also intact). The two multiply:
   `P(usable window showing a central SNP) = p * (1-p)^{2m}` — e.g. at
   p=1/31 (~3.2%), k=31: `p * (1-p)^30 ~= 0.0323 * 0.374 ~= 1.2%`, i.e. about
   1 window in 83, not 1 in 31 (which is only the centre-mutated fraction,
   before requiring intact flanks).
2. **Bias toward isolated SNPs.** Two SNPs within k of each other disqualify
   each other's flanks. Hypervariable regions are invisible by construction.
3. **Indels are invisible.** A frameshift destroys k-mer matches in a block;
   this channel captures substitutions only. Indel divergence shows up as lost
   shared k-mers (lower Jaccard/Mash), not as SNP signal.
4. **k-dependent specificity.** "A k-mer match implies common ancestry" is
   quantitative. For a 3 Gbp genome, expected random flank-30 collisions
   (k=31): `(3e9)^2 / 4^30 ~= 8` — negligible. At k=21: `(3e9)^2 / 4^20 ~= 2e6`
   — no longer negligible. k=31 is safe; k<=21 is marginal to unreliable for
   large genomes. The large k that guarantees homology is the same k that
   shrinks the detectable-divergence window — an inherent tension.

## Implementation: avoid materializing a de Bruijn graph

A central-SNP pair is topologically a simple bubble in the colored de Bruijn
graph (source/sink k-mer shared, two length-k branches differing only at the
midpoint). Classical bubble-calling (Cortex/discoSNP-style) finds these, but
requires the graph — nodes plus adjacency for ~10^9 colored k-mers — resident
in memory. **Rejected**: prohibitive RAM for this project's scale.

A naive per-pair generalisation of variant lookup across n genomes (query each
non-shared k-mer's 3 central variants against every counterpart genome's
index) costs `O(n^2 . N . 3)` random lookups, with the same k-mer's 3 variants
regenerated and requeried once per counterpart genome — pure redundant work.
**Rejected** as the basis for an n-genome design.

## Implementation: sequential per-partition sweep (no scratch, no graph)

`KmerIndex::distance()` already opens every partition's `presence_store`/
`count_store` simultaneously, memory-mapped, into one `LayeredStore`
(`distance.rs:73-77`). "Querying another partition" is therefore not a new
I/O pattern to design — it is the same O(1) MPHF+evidence lookup the `query`
command already performs at scale. This lets the SNP tally be computed with
**no scratch files and no auxiliary graph**, by sweeping partitions once each
as a source:

1. For source partition `p`, enumerate its **distinct** k-mers (one per MPHF
   slot; each already carries its full multi-genome presence/count vector —
   no need to explode per (k-mer, genome) occurrence).
2. For each, generate the 3 central-substitution variants and **canonicalise
   each independently** (`min(kmer, revcomp)`, exactly as any normal query) —
   this avoids the orientation edge case a masked-flank grouping would have
   (a substitution that flips canonical orientation is handled correctly
   because each variant is canonicalised on its own, not inferred from a
   fixed-orientation flank key).
3. Compute each variant's target partition `q` via its minimizer; batch/sort
   the partition's outgoing variant queries by `q` for locality.
4. Look up each variant in `q`'s already-mmap'd MPHF+evidence; on a hit,
   combine the source's presence vector (base `a`) with the variant's
   presence vector (base `b`): for every `i` carrying `a` and `j` carrying
   `b`, `tally[i,j][a,b] += 1`.

**Deduplication needs no persisted state.** Sweeping partitions in a fixed
order `p = 0, 1, ..., P-1` and only acting on a variant when its target
partition `q >= p` guarantees each unordered SNP pair is counted exactly
once: a pair with `q < p` was already resolved earlier, when `q` was itself
the source partition and `p` (being `>= q`) was a valid forward target. No
cross-partition flag array is needed — the sweep order *is* the
deduplication rule. Within the same partition (`q == p`), a lightweight
transient tie-break suffices: either a `#slots(p)`-bit scratch flag reset per
partition, or simply comparing the two k-mers' raw `u64` encodings and only
counting when `kmer_source < kmer_variant` — no storage at all.

This is a distinct computation stage, not a `partial_*` in the existing
additive-by-partition sense: step 3-4 read across partition boundaries by
construction, unlike the row-local `partial_jaccard`/`partial_threshold_jaccard`
primitives. But it requires no new index files, permanent or scratch:
`unitigs.bin`, `mphf.bin`, `evidence.bin`, and the presence/count columns are
read as-is, and the only extra memory is the current partition's small
outgoing-query batch (`#kmers(p) * 3`, released once `p` is done) plus the
persistent `tally` accumulator (`n^2 * 16` counters, see above).

**Outer loop (over source partitions `p`) must stay sequential.** Two
independent reasons, not just one: (a) memory — the bounded-footprint claim
above only holds with one partition's outgoing-query batch in flight; running
`T` source partitions concurrently multiplies that batch by `T`, exactly the
blowup the design avoids; (b) correctness — the `q >= p` deduplication rule
requires partitions to be claimed as sources in a fixed order; running `p1 <
p2` concurrently gives no guarantee `p1` has finished claiming its `q >= p1`
targets before `p2` starts claiming its own, breaking the "counted exactly
once" property.

**Inner loop (target-partition lookups for a fixed `p`) parallelises safely.**
Each lookup is an O(1) read against an already-mmap'd structure, independent
of the others, with no growing allocation — no memory blowup, no ordering
dependency between different `q`. The only shared mutable state is `tally`;
give each worker thread a **thread-local partial tally** (fixed `n^2 * 16`
size, independent of partition size) and merge into the global `tally` once
`p`'s inner loop completes — the same reduce-then-merge pattern Rayon already
uses elsewhere in this codebase to open partitions in parallel. Extra memory:
`#threads * n^2 * 16`, negligible (~512 MB at n~1000, 32 threads) and
unrelated to partition size.

**Cost**: `3 * N_distinct` MPHF lookups total across the whole index (each
partition swept once as source) — the same order of magnitude and the same
operation as running `query` over the index's entire k-mer content against
itself, three times. This is the tool's already-optimized regime, not a new
I/O profile to validate.

## Cheaper: subsampling

Since the target is a ratio, restricting the source-partition sweep to a
bottom-`s` hash sketch (only enumerate k-mers with `hash < threshold` as
sources) divides the lookup count by the sampling factor without biasing
`p_hat`. Mash-like tradeoff: rate estimated from a sample, not the full
k-mer set.

## Recommendation

Sequential per-partition sweep (Route D): reuse the already-mmap'd
per-partition MPHF/evidence/presence structures for O(1) variant lookups,
dedup via the fixed sweep-order rule (`q >= p`, plus an in-partition
tie-break), no scratch files, no graph materialisation. Both the SNP
(off-diagonal) and shared (diagonal) counts are accumulated by this same
sweep, under the locus-eligibility rule chosen (raw or paralogy-filtered) —
not reused from the general-purpose `shared_kmers` matrix, whose raw-identity
definition does not apply the same copy-number constraint. Distances (p, JC,
K2P, LogDet) as finalisations of the resulting 4x4 tally, mirroring the
`partial_* -> *_dist_matrix` pattern used for Jaccard/Mash/Bray-Curtis/etc.

## Detailed implementation plan

Grounded in the current codebase. File/type references are anchors, not
prescriptions; adjust to reality when implementing.

### Step 0 — new low-level primitives (`obikseq`)

Two helpers do not yet exist and are prerequisites:

1. **Central neighbours.** `CanonicalKmerOf<L>` already exposes
   `left_canonical_neighbors()` / `right_canonical_neighbors()`
   (`obikseq/src/kmer.rs`), each returning the 4 canonicalised neighbours at
   an end position. Add `central_canonical_neighbors()` returning the 4
   variants at position `m = (k-1)/2` (each independently canonicalised via
   `.canonical()`). The 3 that differ from the source are the query variants;
   skip the identity. Building on `nucleotide(i)` / the raw 2-bit layout keeps
   it O(1).
2. **Lone-k-mer minimiser.** Routing a *synthetic* variant to its partition
   needs its minimiser, but `RollingStat` (`obiskbuilder/src/rolling_stat.rs`)
   only computes minimisers incrementally along a sequence. Add a standalone
   `minimizer(kmer) -> Minimizer` that scans the `k-m+1` m-mer windows
   (`PackedSeq::mmer`, `obikseq/src/packed_seq.rs`), canonicalises each, and
   takes the min by `seq_hash()` — the same selection `RollingStat` performs,
   evaluated once. Partition index is then
   `(minimizer.seq_hash() & (n_partitions - 1)) as usize`, exactly as
   `QueryBatch::from_records` (`obikmer/src/cmd/query.rs:142`); `n_partitions`
   is a power of two so the mask is valid.

### Step 1 — the tally accumulator (`obikindex`)

A `SnpTally` holding, per genome pair, the 4x4 joint count of central bases:
`n * n * 4 * 4` `u64` (or a packed lower-triangular form since it is
symmetric). Provide `merge(&mut self, other: &SnpTally)` for the thread-local
reduce, and accessors yielding, per pair `(i,j)`: total off-diagonal (SNP),
diagonal (shared, i.e. `p_hat`'s denominator minus SNP), transition count
`P`, transversion count `Q`. The diagonal is always populated — it is not an
optional LogDet-only extra, since `p_hat`'s denominator is no longer sourced
from the external `shared_kmers` matrix (see "Locus eligibility" and
"Statistic and correspondence with `shared`" above): the source k-mer's own
presence/count vector, already in hand when it is enumerated, supplies the
diagonal entry directly, at no extra lookup cost.

### Step 2 — the sweep (`obikindex`, new `snp.rs`)

Mirror `distance.rs`: open the presence or count store per partition. But
instead of a per-partition `partial_*`, run the sequential source sweep:

```text
for p in 0..n_partitions:                       # OUTER — sequential
    open source partition p's layers (QueryLayer-style, obikpartitionner)
    enumerate distinct canonical k-mers of p (one per MPHF slot) with their
        presence/count vectors                   # column-major, as query stage 2
    par_iter over these source k-mers:           # INNER — rayon, thread-local tally
        apply eligibility rule to the source's own vector (raw: none;         # diagonal
            stringent: exactly one of the 4 forms present in each genome)     # gate
        for i in genomes eligible with source base a:
            for j in genomes eligible with source base a:
                thread_tally[i,j][a,a] += 1        # diagonal — no extra lookup
        for each of the 3 central variants:
            q = partition_of(variant)
            if q < p: continue                    # dedup: forward targets only
            if q == p and variant <= source.raw(): continue   # in-partition tie-break
            slot = layers[q].find_slot(variant)   # MphfLayer::find, mmap'd
            if hit:
                vb = variant presence/count vector
                apply eligibility rule to vb (as above)
                for i in eligible genomes with source base a:
                    for j in eligible genomes with variant base b:
                        thread_tally[i,j][a,b] += 1
    merge thread-local tallies into global SnpTally
```

The inner lookup is precisely `QueryLayer::find_slot` +
`col_value(g, slot)` (`obikpartitionner/src/query_layer.rs`) — reuse or factor
out that path rather than reimplementing MPHF access. Enumerating "all distinct
k-mers of a partition with their vectors" is the `dump`/`query` stage-2
column-major scan already implemented in `dump_layer.rs` /
`query_partition_with`; factor a reusable iterator if none fits.

`presence_threshold` applies exactly as elsewhere: a genome "carries base b"
iff its count at that slot is `>= presence_threshold` (trivially `>= 1` for
presence indexes).

### Open problem (unresolved, session end — not yet fully convinced)

The `q >= p` / tie-break dedup rule in Step 2's pseudocode above is **flawed**
for the stringent (paralogy-filtered) eligibility rule: it only ever brings
two family members into view at once (the source and one looked-up variant),
never all four simultaneously, and which subset gets compared depends on
partition sweep order. "Exactly one of the 4 forms present in genome A" is a
whole-family property and cannot be decided correctly from a sequence of
pairwise, order-dependent glimpses — the pseudocode above needs revision, not
just the eligibility gate bolted onto it as written.

Direction discussed, **not yet settled**:

1. **Every distinct source k-mer looks up all 3 variants unconditionally**
   (drop the `q < p` skip entirely) so that every observed family member
   independently gathers all 4 vectors (its own + whichever of the 3
   variants exist) at once — a whole-family, order-independent view, computed
   redundantly once per observed member. Same total lookup order of
   magnitude as already budgeted (`3 * N_distinct`), just organised
   differently (no lookup actually skipped, versus the original rule which
   skipped roughly half).
2. **Tie-break after gathering, not before**: only the member whose own
   canonical encoding is the smallest *among the members actually observed*
   (now known, since all were just looked up) writes to `SnpTally`; the
   others silently discard their redundant computation. Deterministic,
   order-independent — as a side effect this also removes the "outer loop
   must stay sequential" constraint from the cost/parallelism discussion
   above, since no step depends on partition processing order any more.
3. **Proposed optimisation**: precompute, once at index build time, a
   compact global (not per-genome) annex per MPHF slot — the count of
   *other* family members observed anywhere in the dataset (0-3). Slots with
   count 0 (majority under low divergence and few genomes, but see the
   scaling caveat below) need no cross-lookup at all: eligibility reduces to
   a local `count == 1` check at that single slot, and only slots with count
   >= 1 enter the 3-lookup sweep machinery above. Revised (see Step 2b
   below): minorant status *is* stored alongside the count after all, on 3
   bits rather than 2 — it comes for free from the same lookups needed to
   count siblings, and storing it lets the sweep discard non-minorant slots
   without re-fetching anything.

**Minorant/sibling-count relationship, worked out precisely.** "Minorant" is
a one-way implication from sibling count, not an equivalence: `0 siblings
=> minorant` (trivially — with no other observed member, the k-mer is by
definition the smallest of the observed set, itself alone), and its
contrapositive `not minorant => >= 1 sibling`. The converse does not hold:
being the minorant says nothing about sibling count — a minorant can have 0,
1, 2 or 3 siblings, all with larger encodings than itself. Consequence: this
confirms, as a logical necessity rather than a heuristic, that a 0-sibling
slot can always write its diagonal contribution with zero ambiguity and no
lookup (it is unconditionally its own minorant) — but it gives no shortcut
for the >= 1-sibling case, where minorant status still requires the actual
comparison of gathered encodings; sibling count alone never determines it.

**When to compute the annex, and cache invalidation.** Sibling count is a
property of the whole set of columns (genomes/groups) currently in the
index, not of any single genome — it cannot be computed correctly at
mono-genome build time (a family may gain siblings, or its minorant may
change, once more genomes are merged in later). Computing it eagerly at
every `merge` would also waste work on intermediate merged states nobody
ever queries. Instead: compute it lazily, on first `distance` call against a
given index, and persist the result alongside that index for subsequent
calls — the same lazy-derived-cache pattern `PersistentBitMatrix` already
uses for `Columnar` -> `Packed`. This requires no explicit invalidation for
`merge` or `filter` (`obikindex/src/merge.rs`, `obikmer/src/cmd/filter.rs`):
both only ever write to a fresh `--output` directory, never mutate an input
index in place, so a re-merged/re-filtered index is simply a new state with
no annex yet. `select --in-place` (`select_layer.rs:139-235`) is the
exception: it aggregates genome columns into groups (Any/All/None/Sum/Min/
Max) by mutating the existing index's files without changing its location.
It does not remove k-mer rows, but it can still change eligibility and
sibling counts derived from those rows (e.g. a `Sum` over several
single-copy genomes can read as multi-copy at the group level). Because it
mutates in place, **`select --in-place` must explicitly invalidate (delete
or mark stale) any cached sibling-count annex for that index** — the one
operation in the current pipeline where this doesn't happen for free.

Not yet convinced this is the right shape, and Step 2's pseudocode above has
not been rewritten to match — flagged for the next pass rather than resolved
here.

### Step 2b — sibling-count / minorant annex (consolidated plan)

Scope: only the precursor annex — not the SNP tally itself, whose Step 2
sweep remains unresolved above. This piece is simpler than the sweep,
because it writes to an independent per-slot value, not a shared
cross-k-mer accumulator, so it needs no dedup/ownership logic at all at this
stage.

**Revised annex encoding — 4-bit presence mask, not 3-bit (minorant +
count).** Superseded after settling the "canonical form of a family"
definition above. The 3-bit design (1 minorant bit + 2-bit sibling count,
§ below, kept for the historical record) had two problems: it discards
*which* variants are present (only how many), so any future consumer
(the SNP sweep, or a stats pass — see below) that needs to know which bases
exist still has to regenerate and blindly re-query all 3 candidates; and
the minorant bit's meaning was tied to whichever member was visited, not to
a fixed reference. Storing instead a **4-bit mask** — one bit per base
(A/C/G/T), set iff that member of the family (labelled relative to the
family's fixed canonical form, i.e. the member with `A` at the centre — see
above) is observed anywhere in the index — fixes both:
- **Sibling count is derived, not stored**: `siblings = popcount(mask) - 1`.
- **Minorant is derived, not stored**: regenerate the family's 4 canonical
  forms from the slot's own k-mer (cheap, no lookup — see above), compare
  the raw encodings of whichever bits are set in the mask, take the
  smallest.
- **A future consumer knows exactly which variants to (re-)query** —
  `popcount(mask) - 1` lookups instead of always 3, and it knows *which*
  3 (or fewer) to issue, not just how many hits to expect.
- The all-zero value (no base present at all) is still logically
  unreachable as a real result — the slot's *own* base is always present in
  its own family — so it remains available as a free "not yet computed"
  sentinel, exactly as before.

1. **Primitive.** Reuse `central_canonical_neighbors()` from Step 0
   unchanged — the 3 canonicalised central-substitution variants of a k-mer
   (plus the identity, i.e. all 4 members of the family — see "Definitions"
   above).
2. **New annex type** (`obicompactvec`, alongside `bitmatrix.rs`): a 4-bit-
   per-slot packed array (the presence mask above), one per partition — same
   on-disk shape family as `PersistentBitMatrix`'s `Packed` variant, but
   simpler (no per-genome columns, a single derived read-only value per
   slot).
   <details><summary>Superseded 3-bit design (historical)</summary>
   3 bits, storing minorant status alongside sibling count directly, since
   it came for free from the same lookups (point 3 below) — 5 real states
   (not-minorant; minorant with 0/1/2/3 siblings) fit in 3 bits (8 states,
   3 unused). This let the future SNP sweep discard a non-minorant slot
   instantly, with no lookup at all. The otherwise-unreachable combination
   "not-minorant + 0 siblings" (0 siblings always implies minorant) doubled
   as the "not yet computed" sentinel. Replaced by the 4-bit mask above,
   which subsumes this benefit (minorant still derivable, now for free at
   read time rather than stored) while also fixing the "which variant"
   blindness.
   </details>
3. **Computation pass** (`obikindex`, new `siblings.rs`): **one
   `obipipeline` run per layer, iterated sequentially over the index's
   layers** — settled after two false starts, worth recording both.
   - *False start 1*: "fully parallel over every partition/slot at once,
     no ordering at all". Correctness is fine with this (sibling count and
     minorant are order-independent, unlike the old `q >= p` dedup they
     replace), but it reintroduces, at a larger scale, exactly the
     memory-blowup the original Step 2 sweep's sequential-outer-loop
     constraint existed to prevent: scattering every source partition at
     once multiplies the in-flight outgoing-query volume by the number of
     partitions.
   - *False start 2*: push the layer loop itself into the pipeline (source
     = the index's layers, a first `Flat` stage expands each layer into
     its k-mers). `obipipeline`'s scheduler already bounds memory on its
     own — it dispatches every item through a **shared** worker pool at
     each stage boundary (`scheduler.rs:217-372`, `dispatch()` into a
     common `worker_tx` queue, any free worker picks up any pending item;
     not "one worker owns a chunk end to end"), with a biased `Select`
     that prioritises draining items already advanced in the chain over
     admitting new source items (`scheduler.rs:271-282`: stage results
     outrank the source, "vider le pipeline en priorité" / "dernier
     recours" for new data) — so bounded channel `capacity` plus this
     drain-first bias already caps in-flight work without any external
     sequential discipline. Correct, but it means k-mers from several
     layers can be completing concurrently, so the sink would need to
     track several open per-layer annex-file writers at once — real,
     avoidable complexity.
   - **Settled design**: keep the layer loop external and sequential —
     not for memory (the pipeline's own `capacity`/priority mechanism
     already provides that, for free, regardless), but so each pipeline
     run's sink targets exactly one layer's annex file, no concurrent
     multi-writer bookkeeping. Per layer: source = that layer's distinct
     k-mers; a `Flat` (1->N) stage generates the 3 central variants of a
     k-mer, each tagged with its origin (local slot); a transform stage
     routes each variant to its target partition (unchanged per-k-mer
     minimiser); a transform stage performs the lookup (existence-only —
     `find_slot` hit/miss, cheaper than the SNP sweep's full column
     fetch); a final stage/sink folds each answer into its origin's
     running state (below) and, once a layer's k-mers are all resolved,
     flushes the completed array to that layer's annex file. Many small,
     single-purpose stages on purpose, to let the scheduler interleave
     them finely across many in-flight items — this deliberately does
     **not** mirror how `obipipeline` is used elsewhere today: `query.rs`'s
     `process_chunk` lumps parse+route+query+serialise into one closure
     (`query.rs:325,743-758`), and `scatter.rs` only pipelines file-
     reading/superkmer construction, routing partitions afterwards in a
     plain sequential loop (`KmerPartition::write_batch`,
     `partition.rs:140`) — both under-use the fine-grained scheduling the
     mechanism offers, so they are not precedents to copy, only existing
     (and arguably improvable, out of scope here) usages. Cross-partition
     lookups (querying another layer's MPHF for a variant) remain
     necessary as before — only the *output* side is kept single-layer.
   - **Reconciliation**: processed at the granularity of one *answer batch
     per destination partition*, not one source k-mer at a time — this is a
     proper shuffle, not a per-k-mer wait. Each source partition `p` holds a
     small array of running states `(minorant = true, siblings = 0)`, one
     per local slot, initialised at scatter time and **persisting across
     however many destination-partition batches answer it** (up to 3, one
     per variant, not necessarily all from the same `q`). Every scattered
     query carries an origin tag (source partition + local slot) so its
     answer can be routed back. When target partition `q` returns its batch
     (all answers for every query that named `q`, regardless of which source
     k-mer or which source partition they came from), that batch is walked
     once, locally, and each answer updates — via its origin tag — the
     matching entry in *its* source partition's array: a miss changes
     nothing; a hit does `siblings += 1`, and if the found sibling's own
     encoding is smaller than the source's, `minorant = false`. A given
     source k-mer's state is final only once every destination batch
     concerning it has been folded in; its partition's array is flushed to
     the persistent annex once complete. Commutative per entry, so the order
     in which destination batches arrive and get folded in doesn't matter.

   **Open optimisation, not adopted yet — real tradeoff, not a free win.**
   Since looking up sibling `y` from `x`'s visit already yields everything
   needed to fill `y`'s own annex entry too, one visit per *family* could in
   principle replace one visit per *observed family member* — cutting this
   pass's cost roughly by the average family size instead of paying
   `3 * N_distinct` regardless. But it means threads processing different
   source k-mers can end up writing the *same* sibling's slot concurrently —
   the fully independent, ownership-free parallelism of the plan above is
   deliberately traded away for this gain. It stays safe only because the
   computed value for a given slot is deterministic regardless of who
   computes it, so redundant concurrent writes converge to the same
   value — correct as long as each write is atomic, no locking needed — but
   it is a real design complexity increase over "every member redoes its
   own 3 lookups independently," not a strict improvement to adopt by
   default.
4. **Trigger and caching** (`obikindex::KmerIndex`/`distance.rs`): compute
   lazily on first `distance` call for an SNP-family metric against a given
   index; check for an existing annex file first (mirrors
   `PersistentBitMatrix::open()`'s auto-detect-and-fall-back,
   `bitmatrix.rs:264-287`); if absent, run step 3 and persist; if present,
   mmap and reuse.
5. **Invalidation.** `merge` and `filter` always write to a fresh `--output`
   directory (`obikindex/src/merge.rs`, `obikmer/src/cmd/filter.rs`) so a
   re-merged/re-filtered index simply has no annex yet — nothing to
   invalidate. `select --in-place` (`select_layer.rs:139-235`) mutates
   columns of an existing index without changing its location, which can
   change sibling counts without removing rows — it must explicitly delete
   any cached annex for that index as part of its in-place rewrite.
6. **Testing**: hand-built tiny indexes with known sibling counts (0-3);
   order-independence (recompute twice on a static index, identical
   result, given the fully-parallel no-ownership design); invalidation
   (annex absent/correctly recomputed after `select --in-place`); once
   Step 2's sweep is fixed, a regression check that sibling_count == 0
   slots are never looked up cross-partition during the sweep.

Cost: `3 * N_distinct` existence-only lookups, computed once per index
state and amortised over every subsequent `distance` call that reuses the
cached annex — cheaper per-lookup than the sweep itself (hit/miss only, no
column fetch).

### Step 3 — finalisation (`obikindex`)

From the global `SnpTally` alone (diagonal and off-diagonal both populated by
the sweep, see Step 1/2 — no dependency on the external `shared_kmers`
matrix), derive n x n distance matrices, each a pure function of the
accumulated counts (same shape as `jaccard_to_mash`):

- `p_hat[i,j] = SNP / (SNP + shared)`
- Jukes-Cantor, Kimura-2P (from `P`, `Q`), optionally LogDet (needs the
  diagonal + base-composition margins).

Guard the singularities (`p >= 3/4` for JC, `1-2P-Q <= 0` or `1-2Q <= 0` for
K2P) by clamping to a max distance, as `jaccard_to_mash` clamps `J <= 0`.

### Step 4 — surfacing (`obikindex` + `obikmer` CLI)

These metrics do not fit `DistanceMetric`'s current `LayeredStore`-partial
dispatch (they need the cross-partition sweep and produce a different
intermediate). Two options, to decide:

- **(a)** New `DistanceMetric` variants (`Pdistance`, `JukesCantor`,
  `Kimura2P`, `LogDet`) whose `KmerIndex::distance` arm calls the sweep
  (`snp.rs`) instead of the partial path, still returning `DistanceOutput`.
  Keeps one CLI surface (`--metric jukes-cantor`), at the cost of a branch in
  `distance()` that ignores the `LayeredStore` it built.
- **(b)** A dedicated pathway (`KmerIndex::snp_distance`) and a distinct CLI
  entry, if mixing a cross-partition sweep into the partition-local `distance`
  command is judged architecturally muddy.

Recommendation: (a) for user ergonomics (all pairwise distances under
`distance`, all feeding NJ/UPGMA/`--shared-kmers` unchanged), but compute the
sweep lazily only when an SNP-family metric is requested, so the existing
metrics keep their partition-local fast path untouched.

### Step 5 — subsampling flag

Add `--snp-sample <fraction>` (or a bottom-`s` hash threshold): restrict the
source-k-mer enumeration in Step 2 to `seq_hash(kmer) < threshold`. Divides
lookups proportionally; `p_hat` is unbiased. Off by default (exact).

### Testing

- **Primitive unit tests**: `central_canonical_neighbors` on hand-checked
  k-mers incl. palindrome-boundary cases; lone-k-mer `minimizer` against
  `RollingStat`'s incremental result on the same k-mer.
- **End-to-end tiny index**: two 1-genome indexes differing by a handful of
  known isolated SNPs (transitions and transversions placed by hand), assert
  exact `SNP`, `P`, `Q` counts and the resulting JC/K2P values.
- **Dedup invariant**: assert the tally is identical regardless of genome/
  partition order and that no pair is double-counted (compare against a
  brute-force all-pairs reference on a small index).
- **Subsampling**: `p_hat` within sampling error of the exact run.

### Suggested phasing

1. Step 0 primitives + their unit tests (self-contained, no distance wiring).
   This also unblocks the long-declared-but-unimplemented `query --mismatch`
   (`obikmer/src/cmd/query.rs:676`, currently a warning), which needs the same
   neighbour + routing machinery.
2. `SnpTally` + finalisation math with a brute-force (non-swept) reference
   backend, validated on a tiny index.
3. The real per-partition sweep (Step 2) behind the same finalisation; assert
   it matches the brute-force backend.
4. CLI surfacing (Step 4a) and NJ/UPGMA integration (already generic over the
   matrix).
5. Subsampling (Step 5).

## References

The Mash mutation-rate model this discussion contrasts with:
[@Mash-distances-doc; @Fan2015-mash-formula].