feat(spr,drift): exact directional EW scorer — spr_search default + drift/sector opt-in + findings

dev/benchmarks/hamilton_drift_exactness_array.sh

@@ -0,0 +1,44 @@
+#!/bin/bash
+# Authoritative matched-wall A/B gate for the exact-directional EW drift scorer.
+# One array task per dataset; each runs the committed harness
+# dev/profiling/drift-exactness-gate-bench.R with production seeds/budgets,
+# union (default) vs TS_DRIFT_EXACT, from identical per-seed local-opt starts.
+# Submit chained: sbatch --dependency=afterok:<build-jobid> this.sh
+#SBATCH --job-name=ts-driftexact-gate
+#SBATCH -p shared
+#SBATCH -n 1
+#SBATCH --mem=4G
+#SBATCH --time=4:00:00
+#SBATCH --array=0-2
+#SBATCH --output=/nobackup/%u/TreeSearch/logs/driftexact_gate_%A_%a.out
+#SBATCH --error=/nobackup/%u/TreeSearch/logs/driftexact_gate_%A_%a.err
+
+set -e
+module load r/4.5.1
+module load gcc/14.2
+export OMP_NUM_THREADS=1
+export OPENBLAS_NUM_THREADS=1
+# Dependencies + TreeTools come from the shared lib; our TreeSearch is loaded
+# explicitly from LIB via the harness's lib.loc=GATE_LIB (so this shared path
+# resolves deps without shadowing the branch build).
+export R_LIBS_USER=/nobackup/$USER/TreeSearch/lib
+
+LIB=/nobackup/$USER/TreeSearch/lib-driftexact
+REPO=/nobackup/$USER/TreeSearch-driftexact
+OUTDIR=/nobackup/$USER/TreeSearch/drift-exactness/out
+mkdir -p "$OUTDIR"
+cd "$REPO"                       # data-raw/*.nex + the harness live here
+
+DATASETS=(Zhu2013 Zanol2014 Dikow2009)
+DS=${DATASETS[$SLURM_ARRAY_TASK_ID]}
+
+# 30 seeds; nCycles sweep extended into the saturated regime to resolve the
+# plateau-vs-climb crossover the 8-seed local pilot showed.
+GATE_LIB="$LIB" \
+GATE_SEEDS=30 \
+GATE_NCYC=8,16,32,64,128,256 \
+GATE_DATA="$DS" \
+GATE_OUT="$OUTDIR/${DS}.csv" \
+  Rscript dev/profiling/drift-exactness-gate-bench.R
+
+echo "DONE $DS -> $OUTDIR/${DS}.csv"

dev/benchmarks/hamilton_drift_exactness_build.sh

@@ -0,0 +1,41 @@
+#!/bin/bash
+# Build TreeSearch from the drift-exactness branch into a DEDICATED shared lib
+# (does NOT clobber the cpp-search shared lib other sessions use). Chain the
+# gate array on afterok of this job.
+#SBATCH --job-name=ts-driftexact-build
+#SBATCH -p shared
+#SBATCH -n 4
+#SBATCH --mem=8G
+#SBATCH --time=0:45:00
+#SBATCH --output=/nobackup/%u/TreeSearch/logs/driftexact_build_%j.out
+#SBATCH --error=/nobackup/%u/TreeSearch/logs/driftexact_build_%j.err
+
+set -e
+module load r/4.5.1
+module load gcc/14.2
+
+BRANCH=claude/happy-jennings-e7a165
+LIB=/nobackup/$USER/TreeSearch/lib-driftexact
+REPO=/nobackup/$USER/TreeSearch-driftexact
+# Dependencies (Rcpp, TreeTools, ape, ...) live in the shared TreeSearch/lib;
+# resolve them from there while installing OUR TreeSearch into the dedicated LIB
+# (so we never clobber the shared lib other sessions use).
+export R_LIBS_USER=/nobackup/$USER/TreeSearch/lib
+mkdir -p "$LIB" /nobackup/$USER/TreeSearch/logs
+
+if [ ! -d "$REPO/.git" ]; then
+  git clone https://github.com/ms609/TreeSearch.git "$REPO"
+fi
+cd "$REPO"
+git fetch origin "$BRANCH"
+git checkout "$BRANCH"
+git reset --hard "origin/$BRANCH"
+echo "Git HEAD: $(git log --oneline -1)"
+
+rm -f src/*.o src/*.so
+R CMD build --no-build-vignettes --no-manual --no-resave-data .
+R CMD INSTALL --library="$LIB" TreeSearch_*.tar.gz
+rc=$?
+rm -f TreeSearch_*.tar.gz
+echo "INSTALL exit: $rc; version: $(R_LIBS_USER=$LIB Rscript -e 'cat(as.character(packageVersion("TreeSearch")))' 2>/dev/null)"
+exit $rc

dev/benchmarks/hamilton_drift_exactness_recipe.sh

@@ -0,0 +1,28 @@
+#!/bin/bash
+# Recipe-level whole-search test (route B'): drift scorer choice in wall-limited
+# thorough on TRAINING-split MorphoBank (EW-recoded). One catalogue key per task.
+# Reuses lib-driftexact (no rebuild). TRAINING split only (sequestration).
+#SBATCH --job-name=ts-driftexact-recipe
+#SBATCH -p shared
+#SBATCH -n 1
+#SBATCH --mem=8G
+#SBATCH --time=12:00:00
+#SBATCH --array=0-7
+#SBATCH --output=/nobackup/%u/TreeSearch/logs/driftexact_recipe_%A_%a.out
+#SBATCH --error=/nobackup/%u/TreeSearch/logs/driftexact_recipe_%A_%a.err
+set -e
+module load r/4.5.1; module load gcc/14.2
+export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1
+export R_LIBS_USER=/nobackup/$USER/TreeSearch/lib
+LIB=/nobackup/$USER/TreeSearch/lib-driftexact
+REPO=/nobackup/$USER/TreeSearch-driftexact
+OUTDIR=/nobackup/$USER/TreeSearch/drift-exactness/out
+mkdir -p "$OUTDIR"; cd "$REPO"
+# Deterministic mid-size (65-135t) training-split keys, spread across sizes.
+KEYS=(project3617 project589 project4624 project2124 project5201 project3807 project3602 project4614)
+KEY=${KEYS[$SLURM_ARRAY_TASK_ID]}
+GATE_LIB="$LIB" GATE_KEY="$KEY" GATE_CAT=/nobackup/$USER/floor/mbank_catalogue.csv \
+  GATE_SEEDS=6 GATE_REPS=2,4,8 \
+  GATE_OUT="$OUTDIR/recipe_${KEY}.csv" \
+  Rscript dev/profiling/drift-exactness-recipe-bench.R
+echo "DONE recipe $KEY"

dev/benchmarks/hamilton_drift_exactness_replicate.sh

@@ -0,0 +1,26 @@
+#!/bin/bash
+# Route-3 gap-closer (arm A): +replicate marginal-value. Same datasets as route3
+# for direct frontier comparison. Reuses lib-driftexact (no rebuild).
+#SBATCH --job-name=ts-driftexact-repl
+#SBATCH -p shared
+#SBATCH -n 1
+#SBATCH --mem=6G
+#SBATCH --time=8:00:00
+#SBATCH --array=0-3
+#SBATCH --output=/nobackup/%u/TreeSearch/logs/driftexact_repl_%A_%a.out
+#SBATCH --error=/nobackup/%u/TreeSearch/logs/driftexact_repl_%A_%a.err
+set -e
+module load r/4.5.1; module load gcc/14.2
+export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1
+export R_LIBS_USER=/nobackup/$USER/TreeSearch/lib
+LIB=/nobackup/$USER/TreeSearch/lib-driftexact
+REPO=/nobackup/$USER/TreeSearch-driftexact
+OUTDIR=/nobackup/$USER/TreeSearch/drift-exactness/out
+mkdir -p "$OUTDIR"; cd "$REPO"
+DATASETS=(Zhu2013 Zanol2014 Dikow2009 mbank_X30754)
+DS=${DATASETS[$SLURM_ARRAY_TASK_ID]}
+SEEDS=20; if [ "$DS" = "mbank_X30754" ]; then SEEDS=10; fi
+GATE_LIB="$LIB" GATE_SEEDS=$SEEDS GATE_KMAX=16 GATE_DATA="$DS" \
+  GATE_OUT="$OUTDIR/repl_${DS}.csv" \
+  Rscript dev/profiling/drift-exactness-replicate-bench.R
+echo "DONE repl $DS"

dev/benchmarks/hamilton_drift_exactness_route3.sh

@@ -0,0 +1,32 @@
+#!/bin/bash
+# Route-3 marginal-value experiment. One array task per dataset; runs the
+# committed harness dev/profiling/drift-exactness-route3-bench.R (arms:
+# ratchet / drift_exact / drift_union) from a shared per-seed local optimum,
+# cycles knob, matched-wall frontier. mbank_X30754 (180t) is the large/hard
+# anchor where a TBR local-opt has real room; the 74-88t sets are easy controls.
+# Reuses lib-driftexact (no rebuild).
+#SBATCH --job-name=ts-driftexact-route3
+#SBATCH -p shared
+#SBATCH -n 1
+#SBATCH --mem=6G
+#SBATCH --time=8:00:00
+#SBATCH --array=0-3
+#SBATCH --output=/nobackup/%u/TreeSearch/logs/driftexact_route3_%A_%a.out
+#SBATCH --error=/nobackup/%u/TreeSearch/logs/driftexact_route3_%A_%a.err
+set -e
+module load r/4.5.1
+module load gcc/14.2
+export OMP_NUM_THREADS=1
+export OPENBLAS_NUM_THREADS=1
+export R_LIBS_USER=/nobackup/$USER/TreeSearch/lib
+LIB=/nobackup/$USER/TreeSearch/lib-driftexact
+REPO=/nobackup/$USER/TreeSearch-driftexact
+OUTDIR=/nobackup/$USER/TreeSearch/drift-exactness/out
+mkdir -p "$OUTDIR"; cd "$REPO"
+DATASETS=(Zhu2013 Zanol2014 Dikow2009 mbank_X30754)
+DS=${DATASETS[$SLURM_ARRAY_TASK_ID]}
+SEEDS=20; if [ "$DS" = "mbank_X30754" ]; then SEEDS=10; fi
+GATE_LIB="$LIB" GATE_SEEDS=$SEEDS GATE_CYC=1,2,4,8,16 GATE_DATA="$DS" \
+  GATE_OUT="$OUTDIR/route3_${DS}.csv" \
+  Rscript dev/profiling/drift-exactness-route3-bench.R
+echo "DONE route3 $DS"

dev/benchmarks/hamilton_drift_exactness_sector.sh

@@ -0,0 +1,33 @@
+#!/bin/bash
+# Sector-godrift isolation gate (route B). One array task per dataset; runs the
+# committed harness dev/profiling/drift-exactness-sector-bench.R (3 arms:
+# none/union/exact) at production seeds/replicates. Reuses lib-driftexact (code
+# unchanged since f7f17e16 -- docs-only commits since), so NO rebuild: just
+# `git pull` the repo for the new harness, then submit this.
+#SBATCH --job-name=ts-driftexact-sector
+#SBATCH -p shared
+#SBATCH -n 1
+#SBATCH --mem=6G
+#SBATCH --time=8:00:00
+#SBATCH --array=0-3
+#SBATCH --output=/nobackup/%u/TreeSearch/logs/driftexact_sector_%A_%a.out
+#SBATCH --error=/nobackup/%u/TreeSearch/logs/driftexact_sector_%A_%a.err
+set -e
+module load r/4.5.1
+module load gcc/14.2
+export OMP_NUM_THREADS=1
+export OPENBLAS_NUM_THREADS=1
+export R_LIBS_USER=/nobackup/$USER/TreeSearch/lib
+LIB=/nobackup/$USER/TreeSearch/lib-driftexact
+REPO=/nobackup/$USER/TreeSearch-driftexact
+OUTDIR=/nobackup/$USER/TreeSearch/drift-exactness/out
+mkdir -p "$OUTDIR"
+cd "$REPO"
+DATASETS=(Zhu2013 Zanol2014 Dikow2009 mbank_X30754)
+DS=${DATASETS[$SLURM_ARRAY_TASK_ID]}
+# Fewer seeds for the 180-tip mbank (each thorough search is far slower).
+SEEDS=15; if [ "$DS" = "mbank_X30754" ]; then SEEDS=8; fi
+GATE_LIB="$LIB" GATE_SEEDS=$SEEDS GATE_REPS=1,2,4 GATE_DATA="$DS" \
+  GATE_OUT="$OUTDIR/sector_${DS}.csv" \
+  Rscript dev/profiling/drift-exactness-sector-bench.R
+echo "DONE sector $DS"

dev/profiling/drift-exactness-census.R

@@ -0,0 +1,64 @@
+#!/usr/bin/env Rscript
+# Drift EW decision-flip census: quantify how far the deployed union-of-finals
+# drift scorer diverges from the exact directional path, on the three gate
+# datasets recoded to pure Fitch (EW).  Prints the DRIFT-SCANCHK line for the
+# union (default) and exact paths.  Local, bounded -- severity probe only.
+suppressMessages({
+  library(TreeSearch, lib.loc = ".agent-hj")
+  library(TreeTools)
+})
+
+fitchPhy <- function(p) {
+  m <- PhyDatToMatrix(p, ambigNA = FALSE)
+  m[m == "-"] <- "?"                       # inapplicable -> missing => standard Fitch
+  MatrixToPhyDat(m)
+}
+
+makeDs <- function(dataset) {
+  at <- attributes(dataset)
+  list(contrast = at$contrast,
+       tip_data = matrix(unlist(dataset, use.names = FALSE),
+                         nrow = length(dataset), byrow = TRUE),
+       weight = at$weight,
+       levels = at$levels)
+}
+
+datasets <- c("Zhu2013", "Zanol2014", "Dikow2009")
+
+for (nm in datasets) {
+  phy <- fitchPhy(ReadAsPhyDat(file.path("data-raw", paste0(nm, ".nex"))))
+  ds <- makeDs(phy)
+  nTip <- length(phy)
+  hasGap <- "-" %in% ds$levels
+  cat(sprintf("\n=== %s : %d tips, %d chars(weighted), levels={%s} hasGap=%s ===\n",
+              nm, nTip, sum(ds$weight), paste(ds$levels, collapse = ""), hasGap))
+
+  # Production-representative start: random tree -> TBR to a local optimum, then
+  # census drift from there (drift runs AFTER a hill-climb in real searches).
+  set.seed(1L)
+  start <- RandomTree(nTip, root = TRUE)
+  # Align edge tip index i with tip_data row i (= names(phy)[i]).
+  start$tip.label <- names(phy)
+  tbr <- TreeSearch:::ts_tbr_search(start$edge, ds$contrast, ds$tip_data,
+                                    ds$weight, ds$levels, maxHits = 5L)
+  startTree <- list(edge = tbr$edge)
+  cat(sprintf("  local-opt score after TBR = %.0f\n", tbr$score))
+
+  runDrift <- function(exact) {
+    if (exact) Sys.setenv(TS_DRIFT_EXACT = "1") else Sys.unsetenv("TS_DRIFT_EXACT")
+    Sys.setenv(TS_DRIFT_SCANCHK = "1")
+    set.seed(42L)
+    r <- TreeSearch:::ts_drift_search(startTree$edge, ds$contrast, ds$tip_data,
+                                      ds$weight, ds$levels,
+                                      nCycles = 12L, afdLimit = 3L,
+                                      rfdLimit = 0.1, maxHits = 1L)
+    Sys.unsetenv("TS_DRIFT_SCANCHK"); Sys.unsetenv("TS_DRIFT_EXACT")
+    cat(sprintf("  [%s] final drift score = %.0f (drift_moves=%d tbr_moves=%d)\n",
+                if (exact) "EXACT" else "UNION", r$score,
+                r$total_drift_moves, r$total_tbr_moves))
+    invisible(r)
+  }
+  runDrift(FALSE)  # union (deployed default) -- severity
+  runDrift(TRUE)   # exact -- wiring validation (applied_mismatch must be 0)
+}
+cat("\nDONE\n")

dev/profiling/drift-exactness-frontier-combined.R

@@ -0,0 +1,27 @@
+# Combined frontier: original nc{8..256} + extension nc{512,1024}, matched wall.
+dir <- "dev/profiling"
+read1 <- function(f) read.csv(file.path(dir, f))
+scoreAtWall <- function(d, W) {
+  sd <- unique(d$seed)
+  v <- sapply(sd, function(s){x<-d[d$seed==s & d$wall_s<=W,]; if(nrow(x)==0) NA else min(x$score)})
+  mean(v, na.rm=TRUE)
+}
+for (nm in c("Zhu2013","Zanol2014")) {
+  df <- rbind(read1(sprintf("drift-exactness-gate-hamilton-%s.csv", nm)),
+              read1(sprintf("drift-exactness-gate-hamilton-%s-ext.csv", nm)))
+  cat(sprintf("\n== %s : mean score / mean wall by nCycles x scorer ==\n", nm))
+  ag <- aggregate(cbind(score,wall_s)~nCycles+scorer, df, mean)
+  ag <- ag[order(ag$scorer, ag$nCycles),]; ag$score<-round(ag$score,2); ag$wall_s<-round(ag$wall_s,2)
+  print(ag, row.names=FALSE)
+  cat(sprintf("\n-- %s matched-WALL frontier (mean over 30 seeds) --\n", nm))
+  cat(sprintf("  %6s %9s %9s   %s\n","wall","union","exact","winner(margin)"))
+  for (W in c(0.8,1.6,3.2,6.4,9,12,16,20)) {
+    u<-scoreAtWall(df[df$scorer=="union",],W); e<-scoreAtWall(df[df$scorer=="exact",],W)
+    win <- if(any(is.na(c(u,e)))) "-" else if(abs(u-e)<1e-9) "tie" else if(e<u) sprintf("EXACT (%.2f)",u-e) else sprintf("union (%.2f)",e-u)
+    cat(sprintf("  %6.1f %9.2f %9.2f   %s\n",W,u,e,win))
+  }
+  du<-df[df$scorer=="union",]; de<-df[df$scorer=="exact",]
+  cat(sprintf("  best-of-30 any budget: union=%.0f exact=%.0f ; mean@maxNc(%d): union=%.2f exact=%.2f\n",
+      min(du$score),min(de$score),max(df$nCycles),
+      mean(du$score[du$nCycles==max(du$nCycles)]),mean(de$score[de$nCycles==max(de$nCycles)])))
+}

dev/profiling/drift-exactness-frontier.R

@@ -0,0 +1,48 @@
+#!/usr/bin/env Rscript
+# Authoritative matched-wall frontier analysis for the Hamilton gate (30 seeds).
+# Score achievable by wall <= W: per seed+scorer, min score among runs with
+# wall_s <= W (score is monotone non-increasing in nCycles => in wall); mean and
+# best over the 30 seeds. Reports union vs exact and the crossover per dataset.
+dir <- "dev/profiling"
+files <- list.files(dir, pattern = "^drift-exactness-gate-hamilton-.*\\.csv$",
+                    full.names = TRUE)
+df <- do.call(rbind, lapply(files, read.csv))
+cat(sprintf("rows=%d datasets=%s seeds=%d nCycles=%s\n\n",
+            nrow(df), paste(unique(df$dataset), collapse=","),
+            length(unique(df$seed)), paste(sort(unique(df$nCycles)), collapse=",")))
+
+# Raw mean score + mean wall by dataset x nCycles x scorer (matched-cycles view)
+cat("== mean score / mean wall(s) by dataset x nCycles x scorer ==\n")
+ag <- aggregate(cbind(score, wall_s) ~ dataset + nCycles + scorer, df, mean)
+ag <- ag[order(ag$dataset, ag$nCycles, ag$scorer), ]
+ag$score <- round(ag$score, 2); ag$wall_s <- round(ag$wall_s, 3)
+print(ag, row.names = FALSE)
+
+# Matched-WALL frontier
+scoreAtWall <- function(d, W) {
+  seeds <- unique(d$seed)
+  v <- sapply(seeds, function(s){
+    x <- d[d$seed==s & d$wall_s<=W, ]; if (nrow(x)==0) NA else min(x$score)
+  })
+  c(mean=mean(v, na.rm=TRUE), best=min(v, na.rm=TRUE), cov=mean(!is.na(v)))
+}
+Ws <- c(0.1,0.2,0.4,0.8,1.6,3.2,6.4,12.8)
+for (nm in unique(df$dataset)) {
+  cat(sprintf("\n== %s : matched-WALL frontier (mean score over 30 seeds; * = coverage<1) ==\n", nm))
+  cat(sprintf("  %6s %10s %10s   %s\n","wall","union","exact","winner(margin)"))
+  for (W in Ws) {
+    u <- scoreAtWall(df[df$dataset==nm & df$scorer=="union",], W)
+    e <- scoreAtWall(df[df$dataset==nm & df$scorer=="exact",], W)
+    mk <- function(z) if (z["cov"]<1) "*" else " "
+    win <- if (any(is.na(c(u["mean"],e["mean"])))) "-" else
+           if (abs(u["mean"]-e["mean"])<1e-9) "tie" else
+           if (e["mean"]<u["mean"]) sprintf("EXACT (%.2f)", u["mean"]-e["mean"]) else
+                                    sprintf("union (%.2f)", e["mean"]-u["mean"])
+    cat(sprintf("  %6.2f %9.2f%s %9.2f%s   %s\n", W, u["mean"], mk(u), e["mean"], mk(e), win))
+  }
+  # saturated reach (best score at max budget, over seeds)
+  du <- df[df$dataset==nm & df$scorer=="union",]; de <- df[df$dataset==nm & df$scorer=="exact",]
+  cat(sprintf("  saturated best-of-30 (any budget): union=%.0f exact=%.0f | mean@maxNc: union=%.2f exact=%.2f\n",
+      min(du$score), min(de$score),
+      mean(du$score[du$nCycles==max(du$nCycles)]), mean(de$score[de$nCycles==max(de$nCycles)])))
+}