.gitea/workflows/ci.yml

@@ -25,8 +25,14 @@ jobs:
             ~/.cargo/registry
             ~/.cargo/git
             src/target
-          key: ${{ runner.os }}-cargo-${{ hashFiles('src/Cargo.lock') }}
+          # v2: bump this salt to force a clean cache when a stuck/killed job
-          restore-keys: ${{ runner.os }}-cargo-
+          # may have saved a corrupted incremental-compilation `src/target`
+          # (observed 2026-08-11: a stale test binary deadlocked in NUMA
+          # worker startup; a from-scratch rebuild of the same source fixed
+          # it instantly, pointing at cache corruption rather than a source
+          # bug).
+          key: ${{ runner.os }}-cargo-v2-${{ hashFiles('src/Cargo.lock') }}
+          restore-keys: ${{ runner.os }}-cargo-v2-
 
       - name: Build
         run: cargo build --release

docmd/theory/evolutionary_distances.md

@@ -203,6 +203,180 @@ 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.
+
 ## Heterozygosity, ploidy, and consensus-assembly inputs
 
 A within-genome multiplicity signal (more than one of the 4 central forms

src/Cargo.lock

@@ -1715,7 +1715,7 @@ dependencies = [
 
 [[package]]
 name = "obikmer"
-version = "1.1.41"
+version = "1.1.42"
 dependencies = [
  "clap",
  "csv",

src/obikmer/Cargo.toml

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