Merge pull request 'ci: disable numa feature, bump obikmer, and document Sankoff costs' (#65) from push-oruynkvporsn into main

Reviewed-on: #65

.gitea/workflows/ci.yml

@@ -25,17 +25,19 @@ jobs:
             ~/.cargo/registry
             ~/.cargo/git
             src/target
-          # v2: bump this salt to force a clean cache when a stuck/killed job
-          # 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-
 
+      # Both `obikmer` and `obikindex` default to the `numa` feature
+      # (hwloc-based topology detection + CPU pinning), which is only useful
+      # on bare-metal multi-socket indexing hosts. Under this runner's
+      # container/cgroup setup it deadlocks at startup — confirmed live
+      # (2026-08-11): the same test binary hangs indefinitely with `numa` on
+      # and passes instantly, repeatedly, with it off, on the same
+      # container. Disable it for CI; it has nothing to do with test
+      # correctness.
       - name: Build
-        run: cargo build --release
+        run: cargo build --release --no-default-features
 
       - name: Test
-        run: cargo test --release
+        run: cargo test --release --no-default-features

docmd/theory/evolutionary_distances.md

@@ -377,6 +377,102 @@ 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
+
+**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_gl`, 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_gl
+```
+
+Worked examples (`c_sub = c_gl = 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 a flat-cost special case, not a instance of the general formula.**
+Applying the formula naively to e.g. `{A,C,G} -> ∅` would charge `3*c_gl`
+(three independent losses). Reject that: total context disappearance
+(flank-breaking, or true structural loss) is plausibly **one** event, not
+`|X|` of them, so:
+
+```
+cost(X, ∅) = cost(∅, X) = c_ctx   (constant, independent of |X|)
+cost(∅, ∅) = 0
+```
+
+`c_ctx` should not be guessed — it is exactly the value already derived for
+the `>1` bucket in "Context, detectability, and a 3-way ordinal distance per
+pair" above (`(2m*p_hat)/(1-(1-p_hat)^(2m)) + p_hat`), reused rather than
+invented.
+
+**Better construction method: shortest path in a small state graph, not the
+closed-form formula directly.** Build a graph on the 16 states 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_gl`) — then define `cost(X,Y)` as shortest-path distance in that graph,
+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_gl`). 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.
+
+**Transition/transversion refinement.** Split `c_sub` into `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). 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:
+the "Sufficient statistic: 4x4 base-pair tally" section below already plans
+to collect the joint `(centre_i, centre_j)` distribution over resolved (`0`
+or `1`) sites — that tally directly gives the empirical Ts/Tv ratio, from
+which `c_ts`/`c_tv` follow (e.g. `cost ∝ -log(observed rate)`, the standard
+generalised-parsimony step-weighting heuristic), reusing a statistic already
+planned rather than adding a new one.
+
+**`gamma`/`mu` (i.e. `c_gl`/`c_sub`) sensitivity sweep before calibration.**
+No strong prior on whether gain/loss events should cost more or less than a
+point substitution — duplication/deletion rates are not a priori equal to
+point-mutation rates, but the direction isn't obvious, and setting it too
+high risks the parsimony search *eliminating* exactly the
+heterozygosity/paralogy signal the design is meant to tolerate (see
+"Heterozygosity, ploidy..." below). Cheap first step: run the topology at a
+handful of ratios (`0.5, 1, 2, 5`) and check whether it's stable — a robust
+topology across that range is far more trustworthy than one built on a
+single, unvalidated guess. Empirical calibration of `c_gl` from the
+family-size distribution already available (`sibling_annex_stats`) is a
+natural follow-up once the sensitivity sweep shows the topology is worth
+refining further.
+
+**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_gl`/`c_ctx` as
+tunable parameters) is directly usable there — no new tooling required
+before testing it.
+
 ## 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.42"
+version = "1.1.43"
 dependencies = [
  "clap",
  "csv",

src/obikmer/Cargo.toml

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