diff --git a/DevDocMD/theory/evolutionary_distances.md b/DevDocMD/theory/evolutionary_distances.md index 7b1a3506..d8d9bf33 100644 --- a/DevDocMD/theory/evolutionary_distances.md +++ b/DevDocMD/theory/evolutionary_distances.md @@ -2249,9 +2249,85 @@ F84/TN93 and is included below. `snp-` prefix on every CLI value. **`+Γ` rate-heterogeneity modifier, applicable to `snp-jc`, `snp-k2p`, `snp-k81`, `snp-t92`, `snp-f84`, `snp-hky85`, `snp-tn93`** (not `snp-raw`, nothing to correct; not `snp-logdet`, no standard gamma formulation) — same -formula as the base correction, weighted by a shape parameter `α` supplied -by the user (`--gamma-shape `), not estimated by ML. A modifier on -existing values, not a separate enum arm per distance. +formula as the base correction, weighted by a shape parameter `α` either +supplied by the user (`--gamma-shape `) or estimated from the data +(`--gamma-shape auto`/`estimate`, method-of-moments — not ML; see +"Automatic α estimation" below). A modifier on existing values, not a +separate enum arm per distance. + +### Automatic α estimation (`--gamma-shape auto`) + +**Correction (verified against the primary source, 2026-09-11):** Jin & +Nei (1990) itself (*"Limitations of the Evolutionary Parsimony Method of +Phylogenetic Analysis"*, Mol. Biol. Evol. 7(2):82–102 — the paper this +whole `+Γ` correction is cited from, confirmed algebraically to match this +codebase's `corrected_log`/`k2p` exactly against their eq. A4, general, and +A8, the `a = 1` case) contains **no** data-driven α-estimation procedure. +Their own recommendation (p. 98) is a fixed default: *"we suggest that the +gamma distance with a = 1 [eq. A8] be used. However, one may choose a +different gamma distance, estimating a from data. Wilson et al. (1989) +recently used a distance with a = 1/2 for restriction-site data of +mitochondrial DNA in hominoids."* — i.e. Jin & Nei explicitly punt +data-driven estimation to a *different* paper (Wilson et al. 1989), not +read/verified here. The estimator below is therefore **not** "Jin & Nei's +method" under any framing — that attribution (present in an earlier +revision of this section) was wrong, not just under-cited. + +**Implemented** (`PartitionDispersion`, `obikphylo/src/siblings/algorithms/pairwise.rs`) +as an independent method-of-moments estimator, unrelated to any specific +published procedure: pools substitution counts by **partition** rather +than by genome pair, during the same `reduce_pairwise` pass that builds +`PairwiseTally` (no second scan). + +For partition `i`: `n_i` = substitutions pooled over every genome pair, +`L_i` = eligible loci pooled over every genome pair, `R_i = n_i / L_i`. +Modeling among-site rate heterogeneity the same way as the `+Γ` correction +itself (a `Gamma(α, α)`-distributed, mean-1, multiplicative rate shared by +every locus in a partition — the classical Poisson-Gamma/negative-binomial +mixture, the general identity behind gamma-rate-heterogeneity corrections, +also behind Uzzell & Corbin's (1971) original observation that substitution +counts across sites/regions are over-dispersed relative to Poisson): + +\[ +\mathbb{E}[R_i] = \mu \qquad \mathrm{Var}[R_i] = \frac{\mu}{L_i} + \frac{\mu^2}{\alpha} +\] + +Weighting each partition's squared deviation by its own `L_i` cancels the +Poisson term before attributing what's left to `α`: + +\[ +\hat\mu = \frac{\sum_i n_i}{\sum_i L_i} \qquad +V = \frac{\sum_i L_i (R_i-\hat\mu)^2}{\sum_i L_i} \qquad +\hat\alpha = \frac{\hat\mu^2}{V - \hat\mu/\bar L} +\] + +where `\bar L` is the mean partition size. Returns "no estimate" (falls +back to the uncorrected formula, warns) when fewer than 2 partitions have +data, `\hat\mu \le 0`, or `V` doesn't exceed the Poisson floor +`\hat\mu/\bar L` — no detectable over-dispersion, `α` would be unbounded. + +**Caveat, stated explicitly rather than left implicit**: unlike every +closed-form correction in `snp_distance.rs` (each verified line-by-line +against `ape`'s `dist_dna.c`, and now also against Jin & Nei 1990 directly +for the base `+Γ` formula), this estimator is derived from first +principles (the general Poisson-Gamma/negative-binomial identity) with no +primary-source procedure behind it at all — not Jin & Nei's (confirmed +above), and Wilson et al. (1989), the paper they point to instead, hasn't +been read/verified either. Mathematically self-consistent (re-derived and +checked, not guessed), but a from-scratch method, not a literature +implementation. If `--gamma-shape` needs a value with a literature +pedigree rather than an estimated one, Jin & Nei's own stated default, +`α = 1` (`--gamma-shape 1`), is the better-supported choice today. + +Deliberately **not** gated by `--sankoff-ratio-ceiling` the way +`base_pair_tally` is (same precedent as `cardinality_tally` — see its own +doc comment): that filter excludes individual saturated *pairs* from a +composition estimate computed once at the very end, from the complete +`PairwiseTally`; the partition axis needed here only exists transiently, +one partition at a time, while `PairwiseTally` is still being built — long +before any pair's final SNP ratio (and thus its ratio_ceiling eligibility) +is known. `--exclude-genome` isn't applied either, matching +`reduce_pairwise`'s own raw per-pair fold. **Implemented now: `snp-raw`, `snp-jc`, `snp-k2p`, `snp-k81`, `snp-f81`, `snp-t92`, `snp-tn93`, `snp-tv`, all with `+Γ` except `raw`/`tv`** — see @@ -2391,6 +2467,397 @@ mechanical substitution but weren't independently checked against an `ape`-equivalent reference for those three specifically — flagged here, not silently assumed correct. +### α-estimation literature survey and design discussion (2026-09-11) + +Follow-up discussion after implementing `PartitionDispersion::estimate_alpha` +(previous section), triggered by the user pointing out that the literature +conflates two distinct problems: *estimating α from an alignment* versus +*using a given α in a Jin-Nei-style corrected distance*. Jin & Nei (1990) +itself is squarely in the second camp — α is a user-supplied parameter +there, not something their paper estimates (confirmed by reading the paper +directly, see previous section). This section catalogs the actual +α-estimation literature and records which approaches were considered and +why each was or wasn't adopted. + +**Bibliography** (methodological papers, not exhaustive users of `+Γ`): + +| Method | Reference | Needs a tree? | Estimates "true" α or a task-specific one? | Cost | +|---|---|---|---|---| +| ML, continuous Γ (foundational) | **Gu, Fu & Li (1995)**, MBE 12:546–557 — read in full | yes, topology fixed a priori | true α (+ θ jointly, invariant+Γ) | very high, ≤5–6 taxa in practice | +| ML, discrete Γ | Yang (1993/94), MBE 10:1396; JME 39:306 | yes | true α | high | +| Parsimony-based ML on inferred changes | Yang & Kumar (1996), MBE 13:650 | yes, equal branch lengths assumed | true α, biased when branches unequal | low | +| Corrected substitutions/site → NB fit | **Gu & Zhang (1997)**, MBE 14:1106 — read in full | yes, + ancestral reconstruction | true α, close to ML accuracy | low (given a tree) | +| Topology-optimizing grid search | **Guindon & Gascuel (2002)**, MBE 19:534 — read in full | yes, but built *from* the candidate distances themselves (BIONJ) | **not** true α — deliberately biased *upward* for topological accuracy | low–moderate (grid × tree build) | +| Free (non-parametric) rate distribution | Susko et al. (2003), Syst. Biol. 52:594 | yes | tests the Γ assumption itself | high | +| Parsimony counts/site → NB fit | **Wakeley (1993)**, JME 37:613 — read in full | yes, NJ tree + Fitch parsimony ancestral reconstruction | true α, biased *upward* by parsimony (conservative) | low (given a tree) | +| Poisson-Gamma/negative-binomial, foundational | Uzzell & Corbin (1971), Science 172:1089 — **not read**, only seen via citations | unknown — likely also count-per-site/tree-based, unverified | true α (indirect, via NB overdispersion) | unknown | +| Bayesian posterior over α | BPP/MrBayes/BEAST2/RevBayes | yes | posterior, not point estimate | very high | + +**Correction (2026-09-11, after reading Wakeley in full):** an earlier +revision of this table filed Wakeley (1993) under "no tree needed" — +wrong. Wakeley's method: build a tree (neighbor-joining, Saitou & Nei +1987, on the real mtDNA data; coalescent-simulated for the simulation +study), reconstruct ancestral states at internal nodes by Fitch (1971) +parsimony, count the *minimum* number of changes per site this implies, +then fit either a two-rate Poisson mixture or a gamma-distributed-rates +model (⟹ negative binomial, parameters by the methods of Cohen 1965 and +Johnson & Kotz 1969) to the resulting per-site count distribution — the +same general "tree + ancestral reconstruction + per-site counts + NB fit" +shape as Gu & Zhang (1997), just using raw parsimony counts directly +rather than Gu & Zhang's multiple-hit-corrected `k̂`. Confirmed by +simulation (Wakeley's own Fig. 1) that parsimony underestimates both the +mean *and* (more severely) the variance of per-site change counts, which +biases `α` **upward** (toward apparent uniformity) — consistent with Gu & +Zhang's own citation of Wakeley for exactly this bias, and with the +general "parsimony overestimates α" pattern noted throughout this +literature. Real hypervariable-region-1 mtDNA data gave `α̂ ≈ 0.44–0.60` +per data set, `0.47` combined — despite the conservative bias, still +solidly in the "strong heterogeneity" range, which is the paper's own +point: the method is biased but usefully conservative, not useless. + +**Where `PartitionDispersion::estimate_alpha` actually sits**: nowhere in +this table's tree-dependent rows. It shares only the *abstract +mathematical identity* (Poisson-Gamma mixture ⟹ negative binomial +over-dispersion) with Wakeley/Uzzell-Corbin/Gu&Zhang, applied via a +genuinely different, tree-free and ancestral-reconstruction-free route: +pairwise genome-to-genome substitution counts pooled by *partition* +instead of per-site counts pooled across an alignment's columns after +ancestral state reconstruction. As far as this survey has established +(Uzzell & Corbin not read, so not fully ruled out), **no tree and no +ancestral reconstruction of any kind** appears to be a genuinely +distinguishing property of `PartitionDispersion`, not something it +inherits from prior art — consistent with it being an independent +derivation rather than a literature implementation (see previous +section's caveat, now on firmer footing). + +**Why a tree-free estimate is even valid — the identifiability question.** +Raised directly by the user after reading Gu, Fu & Li (1995), which states +explicitly: *"one cannot estimate the rate heterogeneity when the number +of sequences is <3."* This is a real non-identifiability result, not a +practical inconvenience: with only 2 sequences, a site's history is +summarized by one binary observation (differ / don't differ), collapsed +across all sites into two aggregate numbers (transition and transversion +proportions) — divergence time and among-site rate variance are +confounded in that pair of numbers, with no way to separate them. Every +method surveyed above (Gu,Fu&Li, Yang, Wakeley, Gu&Zhang) needs ≥3 +sequences *sharing one evolutionary history* (a tree) specifically to get +multiple independent looks at the *same* site's rate across different +lineages, which is what identifies its variance. + +`PartitionDispersion` doesn't violate this, because it isn't estimating +the same quantity. The literature above estimates rate variation **across +sites, within one shared tree**. `PartitionDispersion` estimates rate +variation **across partitions, pooled over every genome pair in the +index** — each genome pair stands in for one independent draw of a +partition's relative rate, the same identifying role multiple lineages +play in the tree-based methods, just substituting "genome pair" for +"lineage" and "partition" for "site." With `P` genome pairs (potentially +in the hundreds here) all contributing to every partition's pooled count, +the ≥3-samples identifiability requirement is met by the pair count, not +by tree depth — closer in spirit to Jin & Nei's own original motivating +scenario (comparing several genomic *regions* across one fixed panel of +taxa) than to single-tree site-rate estimation. + +This substitution carries its own assumption, which must be named rather +than left implicit: that a partition's *relative* rate (fast/slow) is +reasonably stable **across genome pairs**, not just across sites within +one pair — the direct analogue of the "no lineage-specific rate variation, +only site-specific" assumption every tree-based `+Γ` method already makes, +just moved from lineage→pair and site→partition. If this breaks — e.g. an +index mixing very closely related and very divergent genome pairs, where a +partition's saturation behavior differs qualitatively between the two — +the estimate could be misled in a way a tree-based method would at least +have the topology to detect and a tree-free, pool-everything method +cannot. Not yet tested against real data with strongly heterogeneous +pairwise divergence; worth keeping in mind as the main open validity +question for `--gamma-shape auto`, not the citation question (now settled) +this whole discussion started from. + +**Gu, Fu & Li (1995), read in full — the foundational ML method, not a +candidate to implement.** Models rate variation as invariant+Γ (`θ` = +proportion of invariant sites, `α` = gamma shape among the variable sites), +and derives the exact joint likelihood over nucleotide *configurations* +across all `n` sequences on a **fixed, a-priori-known tree topology** — +for 3 sequences this is closed-form (Jukes-Cantor), for `n` sequences it's +a sum over `4^n` configurations weighted by coefficients tied to the +specific topology's branching structure (their eq. 32–34), with branch +lengths, `α`, and `θ` jointly optimized by Newton-Raphson (Hessian over all +free parameters). The site-rate integral has the same closed form +`E = θ + (1-θ)(1+D/α)^{-α}` that underlies `+Γ` distance corrections +generally (the same mechanism as `corrected_log`, not a coincidence — Jin & +Nei's own formula is the `θ=0` special case applied to a pairwise `D` +rather than a whole-tree one). Explicitly stated by the authors as +infeasible beyond about 5–6 taxa (their own simulations stop at 5); not a +candidate for obikmer's typical genome counts. Notable contribution worth +keeping in mind regardless: `ρ = (1+θα)/(1+α)` is shown to be a more robust +summary of rate heterogeneity than `α` alone whenever `θ` isn't negligible +— `α` and `θ` are confounded (nearly-invariant sites can come from either a +small `α`'s left tail or from a nonzero `θ`), so `α` alone can be very +unstable while `ρ` stays well-behaved. `PartitionDispersion` has no `θ` +term at all (no invariant-site component), so this confound doesn't arise +for it the same way, but it's worth remembering if an invariant+Γ variant +is ever considered. + +**Gu & Zhang (1997), read in full — rejected for now.** Their procedure: +(1) a tree with least-squares branch lengths must already exist; (2) +ancestral states are reconstructed at every internal node (they use Zhang & +Nei 1997's likelihood method, but note plain parsimony would work, just +less accurately); (3) per site, branches are split into "changed"/"unchanged" +given the ancestral reconstruction, and the *expected* substitution count +`k̂` (corrected for multiple hits) is obtained by solving their eq. (7) +(or eq. 12 for the generalized model) numerically — a root-find per site, +not a closed form; (4) the site-level `k̂` values (real-valued, not +integers) are fit to a negative binomial (Uzzell & Corbin 1971) by ML to +get α — again no closed form, numerical optimization. This is a real +per-site pipeline requiring a tree *before* distances/α can be computed — +the reverse of obikmer's current dependency order (`snp_distance` computes +distances, which `--nj`/`--upgma` only turn into a tree *afterward*). Their +own answer to this chicken-and-egg problem is an iterative +distance→tree→α→distance loop, which they themselves flag as expensive. +Adopting this would mean a new ancestral-reconstruction subsystem (even a +parsimony/Fitch-only version) plus a per-family nonlinear solve plus an NB +ML fit — a much bigger addition than `PartitionDispersion`, not something +to build without a concrete need beyond "closer to Gu&Zhang's own +simulated accuracy than our moment estimator." + +**Guindon & Gascuel (2002), read in full — plausible future addition, not implemented.** +Their **Q criterion**: for a candidate α, build a tree (they use BIONJ) from +the α-corrected distance matrix; for every internal branch, group taxa into +the 4 subtrees it separates (A, B, C, D), compute mean inter-subtree +distances \(\bar\delta_{AB}\), \(\bar\delta_{CD}\), etc., and let `S ≤ M ≤ L` +be the three pairings' sums (`{AB,CD}`, `{AC,BD}`, `{AD,BC}`); the branch's +reliability score is `Q_branch = L − M` (zero when the four points are +perfectly tree-additive). The whole-tree `Q` is the mean over internal +branches (negative/zero branches excluded). `α* = argmin_α Q(α)` over a +grid (their run: ~60–100 points, 0.1 to 5000, finer spacing where +sensitivity is highest). Complexity: `O(n²l)` once for the base +frequencies, then `O(n²r)` for the `r` candidate distance matrices, then +`O(n³r)` worst case for building `r` trees and evaluating `Q` on each +(same order as the tree-building step itself, so "free" relative to it). + +Central finding, **the opposite direction from an initial misreading in +this discussion**: `α_opt` (Guindon-Gascuel's topology-optimal value) is +**always ≥ the true α**, not the reverse — underestimating rate +heterogeneity (i.e., picking a larger α than reality) reduces the variance +of distance estimates and thereby improves NJ/BIONJ topological accuracy, +especially when the molecular clock roughly holds (where `α_opt → ∞`, i.e. +no correction at all is topologically best). This is an *empirical +regularity* observed across their tested conditions (η ∈ {0.5, 2.0}, true +α ∈ {0.1, 0.7}, 20-taxon trees) — not a theorem — so treating it as a hard +guarantee would be another unverified-claim mistake of the kind this whole +discussion has been correcting. + +**Proposed integration (not yet implemented)**: since `α_true ≤ α_opt` +empirically, `PartitionDispersion::estimate_alpha()`'s output could seed +the *lower bound* of Guindon-Gascuel's grid (search `[α̂, ∞)` instead of +`[0.1, 5000]`), narrowing the grid substantially and reusing the existing +estimator rather than needing an independent "true α" source — this also +sidesteps Gu & Zhang's tree-first dependency problem entirely, since +Guindon-Gascuel's own grid already builds a tree per candidate. Still +requires: (a) a **BIONJ** implementation (see below), (b) the Q-criterion +subtree-grouping logic, (c) a grid-search driver. Not started. + +**BIONJ availability in the Rust ecosystem — checked 2026-09-11, none found.** +Searched crates.io directly (API query for "bionj": 0 results) and read the +docs.rs pages of the most plausible candidates: `phylotree` (builds/reads/ +manipulates trees but has no distance-matrix reconstruction — no NJ, no +BIONJ, no UPGMA despite web-search summaries claiming otherwise, corrected +after checking the actual docs), `speedytree` (NJ only — Canonical/RapidNJ/ +Hybrid variants, explicitly no BIONJ), `nj` (plain NJ only), `phylo` (no +distance-based reconstruction at all). `obikphylo::neighbor_joining` +(`tree.rs`) is plain NJ, not BIONJ — BIONJ (Gascuel 1997) differs only in +how branch lengths are weighted during agglomeration (accounts for +distance-estimate variance, not just the sum), a modest delta over an +existing NJ implementation rather than a rewrite, but would need to be +written in-house if the Guindon-Gascuel integration above is ever pursued +(their own paper notes other tree-building methods gave similar results in +their simulations, so plain NJ may be an acceptable substitute if BIONJ +proper is judged not worth the effort). + +### Sharing one sample across algorithms (`--sankoff`/`--tnt`/`--phyg`/`--iqtree` + `snp-*` `--distance`) — implemented 2026-09-11 + +**Correctness bug, not just a performance one.** `--sankoff`/`--tnt`/`--phyg`/ +`--iqtree` (via `SiblingExt::sankoff_bundle`) and a `snp-*` `--distance` +(via `SiblingExt::snp_distance`) always consume the *exact same* selection +parameters when requested together in one `obikmer phylo` invocation +(`n`/`free_loss`/`no_ambiguity`/excluded-set/`entropy_bias` — the CLI has +no way to give them different values in one run, see `cmd/phylo/args.rs`). +Before this fix, each independently called `sample_index`, and +`sample_layer`'s `rand::rng()` (`subsample.rs:235`) is a **thread-local +generator that advances across calls, not reseeded each time** — so the +second call silently drew a *different* random sample of sites than the +first, even with identical parameters. Verified empirically (throwaway +8-genome index): running `snp-jc` alone vs. combined with `--sankoff` +produced visibly different distance matrices from the same index/params. +This defeated the actual point of combining these flags — running several +algorithms (Sankoff/TNT/PhyG/IQ-TREE calibration, a `snp-*` distance) on +one *identical* site selection for direct comparison. + +**Fix**: `SankoffBundle` (`sankoff.rs`) now retains its `PairwiseTally`/ +`PartitionDispersion` internally (previously consumed into `raw`/ +`base_pair_tally`/`cardinality_tally` and dropped) and exposes +`SankoffBundle::snp_distance(kind, gamma_shape)`, which computes the +matrix from that *same* already-sampled tally — no second `sample_index` +call. The shared post-sampling logic (the `--gamma-shape` support check, +`alpha` resolution, final matrix build — previously all inline in the +`snp_distance()` free function) was factored into +`pub(crate) fn distance_matrix(tally, dispersion, kind, gamma_shape)` in +`snp_distance.rs`, called by both the standalone `snp_distance()` (after +its own fresh sampling) and `SankoffBundle::snp_distance` (reusing the +bundle's). `cmd/phylo/mod.rs` keeps the `Option` alive past +the Sankoff-family `if` block and, when a `snp-*` `--distance` is also +requested, calls `bundle.snp_distance(...)` instead of +`cache.snp_distance(...)` whenever a bundle was built — since both branches +are driven by the same `args.*` fields, the parameters trivially always +match when both fire; no runtime "do the params match" check needed. +Verified end-to-end: the reused path logs `(reusing the Sankoff bundle's +sample)` and the `snp_distance` stage timer reads `0ms` (`formula(kind)` +is `O(n²)` post-processing, no I/O), vs. ~150-230ms for a fresh sample on +the same tiny test index. + +#### Explicit `--session` (implemented 2026-09-12) + +The in-process fix above only covered reuse *within one command*. +`--session DIR` (`cmd/phylo/args.rs`) names a directory (outside the +index) that persists the `snp-*` `--distance` sample — the +`PairwiseTally`/`PartitionDispersion` pair `sample_index` would otherwise +rebuild from scratch — across *separate* `obikmer phylo` invocations. + +**Crate placement, as planned**: a new **`obiksession`** crate +(`src/obiksession/`), domain-agnostic — no `PairwiseTally`/"site" concept +anywhere in it. `Session::open(dir, params: &[u8], force: bool)` handles +directory lifecycle and the params-conflict check against an opaque byte +blob; `Session::store(name, bytes)`/`Session::restore(name) -> Option` +handle a checksummed artifact cache. `obikphylo` +(`siblings/algorithms/session_cache.rs`, shared by `snp_distance.rs` and +`sankoff.rs` — see below) is the consumer: it serializes its own +`PairwiseTally`/`PartitionDispersion`/`SnpAlignment` to bytes (via `rkyv`) +and hands them to `obiksession`, which never sees their type. + +**Serialization**: `rkyv` 0.8.18, added as planned. One deviation from the +original zero-copy pitch, scoped down deliberately: artifacts are restored +via `rkyv::from_bytes` (full owned deserialize) rather than `rkyv::access` +(zero-copy over the `mmap`). True zero-copy would require every +`PairwiseTally` query method (`categories`, `base_freq`, `pair`, ...) to +work generically over `Archived` as well as the owned type +— a separate, larger change not needed to get the actual win (skipping +`sample_index`'s expensive re-scan of the sibling annex; deserializing an +already-in-memory-sized count array is comparatively cheap CPU, not I/O). +`#[derive(Archive, Serialize, Deserialize)]` was added to `PairStats`, +`PairwiseTally`, and `PartitionDispersion` (`pairwise.rs`) for this. + +**Startup behavior, as designed**: `Session::open` creates `DIR` if +missing (fresh, params recorded); if `DIR` already holds different saved +params, `cmd/phylo` reports a hard error and exits rather than silently +preferring one side, unless `--session-force` is given (which discards the +directory's cached artifacts and starts over under the new params) — +implements exactly the "explicit escape hatch, not silent override" +decision from the design discussion. `obiksession` itself doesn't decode +either side's bytes to produce a field-by-field diff (true to "opaque +blob," left to whoever needs it); `cmd/phylo`'s conflict message is +currently a single generic sentence listing which flags could be the +cause, not a computed diff — an honest v1 simplification, not a limitation +of `obiksession` itself. + +**Per-artifact caching, as designed**: only `PairwiseTally`/ +`PartitionDispersion` are cached (`snp_distance.rs`'s `restore_tally`/ +`store_tally`) — `ratio_ceiling`/`gamma_shape` remain pure post-processing +over whatever tally is in hand, cached or fresh, exactly as planned. + +**Locking, reconsidered and kept as the whole-session `DirLock`.** First +written up as a corner cut ("per-entry locking was sketched, not +implemented"). On actually working through the concrete race, whole-session +locking turned out to be the *correct* choice, not a simplification: +`--session-force` resets a session by deleting its whole `artifacts/` +directory. If the directory lock were only held briefly (during +`Session::open`, then released so a long-running `Session` doesn't block +other processes from touching *different* artifacts), a second process +calling `open(..., force: true)` could acquire the lock, see mismatched +params, and wipe `artifacts/` **while a first process's already-open +`Session` is mid-checkpoint** — genuine corruption risk, not just wasted +work. Holding the lock for the `Session`'s entire lifetime makes that +sequence impossible: a `force` reset must wait for every currently-open +`Session` on that directory to be dropped first. The concurrency this +gives up (two *different* processes reading/writing *different* artifacts +of the same session at once) isn't exercised by anything today — only one +code path (`snp_distance`/`sankoff_bundle`, sequentially per invocation) +uses sessions at all. + +**Per-layer chunked/resumable dumps — implemented.** `subsample::sample_index` +gained a `skip_layers: usize` parameter (skip the first N (partition, +layer) pairs, in its own deterministic iteration order, entirely — no +I/O, no `on_layer` call) and its `on_layer` callback gained a `raw_index` +argument (that pair's 0-based position in the same order), so a caller can +tell exactly how far a run has gotten. `snp_distance`/`sankoff_bundle` +checkpoint their tally/dispersion (and, for `sankoff_bundle`, the +in-progress `SnpAlignment`) to the session every +`session_cache::CHECKPOINT_INTERVAL_LAYERS` (8) layers or +`CHECKPOINT_INTERVAL` (30s) of wall time, whichever comes first — bounding +how much completed scanning a crash can lose, without paying a full +tally-rewrite on every one of what can be thousands of layers. A `progress` +artifact records how many layers are reflected in the checkpoint, with a +`PROGRESS_COMPLETE` (`u64::MAX`) sentinel once every layer has been +processed (distinct from any real layer count, so "finished" is never +confused with "stopped after N layers where N happens to equal the total"). + +Verified with a real kill: built a 10-genome/64-partition index, ran +`--distance snp-jc --session DIR` under `timeout 0.4`, confirmed a +checkpoint at `progress=8` (exactly `CHECKPOINT_INTERVAL_LAYERS`) was on +disk, then re-ran the same command and confirmed it resumed from layer 8 +and completed successfully — 0-effort verification that the mechanism +itself is real, not just plausible-sounding. + +**Important correction from that same test — resume is *not* bit-reproducible, +and shouldn't be.** The output matrix from the resumed run differed from +an uninterrupted control run against the same index/params. Root cause: +`sample_layer`'s random draws come from `rand::rng()`, a per-process, +OS-seeded generator with no continuity across a process boundary — the +layers processed before a kill and the layers processed after a resume +come from two independent random streams. An earlier revision of this doc +(and of the code comments) claimed resume "reproduces the exact same +sample an uninterrupted run would have" — that's wrong, and was corrected +after the user pointed out the actual reasoning error: true randomness has +no notion of continuity in the first place, so there was never a +correctness requirement for a resumed run to match an uninterrupted one +bit-for-bit. The sample produced after a resume is exactly as legitimate a +random draw as one from an uninterrupted run — just a *different* one, +same as running the command twice without `--session` already gives two +different samples today. That variance is explicitly wanted (see below), +not a defect to fix. Comments in `subsample.rs`/`session_cache.rs` were +corrected to state this plainly instead of the false reproducibility claim. + +**Explicitly declined: deterministic/seeded sampling.** Raised as a +possible fix for the above (seed `sample_layer`'s RNG from a per-session +value plus `(partition, layer)`, making a resume reproduce bit-identical +results, and, if applied everywhere, making *any* two runs with identical +parameters produce identical samples). The user explicitly declined this +for now: repeated runs are relied on to measure sampling variance, and a +`--seed` flag would need its own explicit design if wanted later — this is +not an accidental gap, it's a stated preference to revisit only if asked. + +**`SankoffBundle` integration — implemented.** `sankoff_bundle` accepts +the same `session: Option<&Session>` and uses the *same* artifact names +(`session_cache::{TALLY_ARTIFACT, DISPERSION_ARTIFACT, PROGRESS_ARTIFACT}`) +as `snp_distance`, plus its own `ALIGNMENT_ARTIFACT` (`sankoff.rs`) folded +in lockstep with the tally at every checkpoint (same `on_layer` call +produces both, so they can never drift out of sync). Consequence verified +directly: a tally checkpointed by `--sankoff --session DIR` is restored by +a *later, separate* `--distance snp-k2p --session DIR` invocation with no +`--sankoff` at all, and vice versa — cross-consumer reuse, not just +same-command reuse, exactly the synergy this design aimed for. + +**Per-artifact caching, as designed**: only the tally/dispersion (plus, +for `sankoff_bundle`, the alignment) are cached — `ratio_ceiling`/ +`gamma_shape` remain pure post-processing over whatever tally is in hand, +cached or fresh, exactly as planned. + +**Still not done — genuine follow-up, not a justified tradeoff:** +**exhaustive-sample caching for `--subsample`-varying exploration.** The +design discussion's suggestion — cache the exhaustive (`n = None`) sample +once and subsample *from* it on each run instead of one session per exact +`n` — wasn't implemented; a session today is still scoped to one exact +parameter tuple, `--subsample` included. + ### Output format: PHYLIP-relaxed by default for the distance matrix **Implemented.** The primary distance-matrix output diff --git a/UserDocMD/usage/phylo.md b/UserDocMD/usage/phylo.md index f933ddd4..db40e84e 100644 --- a/UserDocMD/usage/phylo.md +++ b/UserDocMD/usage/phylo.md @@ -17,7 +17,7 @@ obikmer phylo INDEX [OPTIONS] | Option | Default | Description | |---|---|---| | `--distance` | `jaccard` | See the two tables below for the full list of accepted values | -| `--gamma-shape ALPHA` | none | Rate-heterogeneity correction, for `snp-*` values that support it (see below). No effect on the other values; rejected if given together with a value that doesn't support it | +| `--gamma-shape ALPHA\|auto` | none | Rate-heterogeneity correction, for `snp-*` values that support it (see below). Either a fixed $\alpha$ or `auto`/`estimate` to fit it from the data (see "Automatic $\alpha$ estimation" below). No effect on the other values; rejected if given together with a value that doesn't support it | | `--presence-threshold` | `1` | Minimum count for a kmer to be considered present, for `jaccard`/`mash` on a count index | | `--csv` | off | Write the matrix as plain CSV instead of the default relaxed-PHYLIP format | | `--shared-kmers` | off | Also write the shared-kmer count matrix. Only valid with a whole-index metric, not a `snp-*` value | @@ -121,6 +121,28 @@ $$d = Q$$ `--gamma-shape ALPHA` applies to every value above except `snp-raw` and `snp-tv`: each $-\ln(x)$ term in the formulas above is replaced by $\alpha\left(x^{-1/\alpha}-1\right)$ (the same weight, same $x$). +### Automatic $\alpha$ estimation (`--gamma-shape auto`) + +`--gamma-shape auto` (or the equivalent `--gamma-shape estimate`) fits $\alpha$ from the index itself instead of requiring a user-supplied value, using a method-of-moments estimator computed once, from the same sampling pass that builds the pairwise substitution tally — no extra scan of the index. + +The estimator pools substitution counts by **partition** rather than by genome pair: for partition $i$, let $n_i$ be the total number of substitutions observed across every genome pair, and $L_i$ the total number of eligible loci across every genome pair, in that partition. Define the partition's observed substitution rate: + +$$R_i = \frac{n_i}{L_i}$$ + +Under a single shared substitution rate with no among-site heterogeneity, each $R_i$ would vary only by Poisson sampling noise. Rate heterogeneity is modeled, as elsewhere in this correction, by a $\mathrm{Gamma}(\alpha,\alpha)$-distributed multiplicative rate (mean 1) shared by every locus in a partition — the classical Poisson–Gamma (negative-binomial) mixture. Under that model: + +$$\mathbb{E}[R_i] = \mu \qquad \mathrm{Var}[R_i] = \frac{\mu}{L_i} + \frac{\mu^2}{\alpha}$$ + +where $\mu$ is the pooled substitution rate across every partition. Weighting each partition's squared deviation by its own $L_i$ removes the first (Poisson) term before attributing what's left to genuine rate heterogeneity: + +$$\hat\mu = \frac{\sum_i n_i}{\sum_i L_i} \qquad V = \frac{\sum_i L_i\,(R_i-\hat\mu)^2}{\sum_i L_i} \qquad \bar L = \frac{\sum_i L_i}{\text{number of partitions}}$$ + +$$\hat\alpha = \frac{\hat\mu^2}{V - \hat\mu/\bar L}$$ + +If the measured variance $V$ doesn't exceed the Poisson floor $\hat\mu/\bar L$ (no detectable over-dispersion across partitions — the data are consistent with a single shared rate), $\alpha$ is left undefined: the correction is silently disabled for that run rather than applying a fabricated value, and a warning is logged. When an estimate is produced, it's logged at the `info` level before the distance matrix is computed. + +Note: this is a method-of-moments estimator derived from the standard Poisson–Gamma relationship between substitution counts and gamma-distributed rate variation, applied per-partition — it is not part of Jin & Nei's (1990) original publication, which only defines the `+Γ` distance formula itself and, absent an estimate, recommends the fixed default $\alpha = 1$ (`--gamma-shape 1`) rather than proposing a way to estimate it from data. `alpha < 1` indicates strong among-site rate heterogeneity (many near-invariant loci, a few fast ones); `alpha` growing large makes the correction converge to the uncorrected formula. + ### Output Without `-o`, the matrix goes to stdout in relaxed-PHYLIP format (`n` on the first line, then one `labelvalue...` row per genome). With `--csv`, the format is instead a header row `genome,,,...` followed by one `