Refresh n=4 cells; extend shape lookup to n=11; Hamilton scripts

Upstream PID fix (TreeDist commit 18371ac "Small-tree comparison fix")
required re-running enumerate_small_n.R for exact n=4..7 cells of the
slim lookup — the prior n=4 PID value (0.9333) was wrong; correct value
is 0.6667 (matching CID at n=4, where the only non-trivial split is
balanced 2|2). data-raw/exact_small_n.rds refreshed.

Shape-conditioned lookup extended from n=4..9 to n=4..11 via Hamilton
HPC (jobs 17155360 build + 17155361_10,11 shape; ~40 min and ~73 min
respectively, nSim=500 at n=10 and nSim=200 at n=11). n=12 still
running on Hamilton (~6h remaining at nSim=100). xz-compressed file
size is now 744 KB - a CRAN size concern; may need to ship a leaner
table or move the larger n cells to inst/extdata via lazy load.

n=4 PID cells of the shape lookup also rebuilt (previously all 1.0
under broken upstream PID).

Scripts vendored under data-raw/:
- hamilton_build.sh: SLURM job to build TreeTools+TreeDist on Hamilton.
- hamilton_shape_lookup.{R,sh}: per-n shape-pair MC and the array job
  driver (n in 10..12 with adaptive nSim).
- extend_shape_lookup.R: local-CPU equivalent (kept for completeness;
  prefer Hamilton for n >= 10).

tests/testthat/test-tree_distance_adjusted.R "Lookup falls back to MC"
test retargeted from n=300 to n=500 (n=300 will be inside the lookup
once the n=4..400 extension completes).

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

Author: ms609

data-raw/exact_small_n.rds

@@ -0,0 +1,80 @@
+# extend_shape_lookup.R -----------------------------------------------------
+# (1) Rebuild n=4 cells of inst/extdata/shapeExpectedDistances.rds: the
+#     previously committed n=4 values were drawn while
+#     PhylogeneticInfoDistance was broken at n=4 (upstream bug, fixed in
+#     TreeDist commit 18371ac7 "Small-tree comparison fix"); PID is all 1.0
+#     in the stale rows.
+# (2) Add n=10 to the shipped shape lookup (nSim = 500 per pair; ~2 hours).
+#
+# n=11 and n=12 deferred to Hamilton: estimated wall time 8 hours and 15+
+# hours respectively at usable nSim.  See PLAN.md "Findings to follow up".
+
+suppressPackageStartupMessages({
+  devtools::load_all("C:/Users/pjjg18/GitHub/worktrees/acid-aphid-TreeTools",
+                     quiet = TRUE)
+  devtools::load_all(".", quiet = TRUE)
+})
+
+out_file <- file.path("inst", "extdata", "shapeExpectedDistances.rds")
+result_list <- readRDS(out_file)
+cat("Loaded:", paste(names(result_list), collapse = ","), "\n")
+
+.ShapePairMC <- function(s1, s2, n, nSim) {
+  lab_pool <- as.character(seq_len(n))
+  cid_v <- numeric(nSim); pid_v <- numeric(nSim)
+  for (i in seq_len(nSim)) {
+    t1 <- TreeTools::Preorder(
+            TreeTools::RootedTreeWithShape(s1, n, sample(lab_pool)))
+    t2 <- TreeTools::Preorder(
+            TreeTools::RootedTreeWithShape(s2, n, sample(lab_pool)))
+    cid_v[i] <- ClusteringInfoDistance(t1, t2, normalize = TRUE)
+    pid_v[i] <- PhylogeneticInfoDistance(t1, t2, normalize = TRUE)
+  }
+  data.frame(
+    shape1 = s1, shape2 = s2,
+    mean_cid = mean(cid_v), sd_cid = sd(cid_v),
+    mean_pid = mean(pid_v), sd_pid = sd(pid_v),
+    n_pairs = nSim
+  )
+}
+
+BuildOneN <- function(n, nSim) {
+  W <- as.integer(TreeTools::NRootedShapes(n))
+  pairs_total <- W * (W + 1L) / 2L
+  cat(sprintf("[n=%d] W=%d shapes, %d pairs, nSim=%d\n", n, W, pairs_total, nSim))
+  t0 <- Sys.time()
+  rows <- vector("list", pairs_total)
+  k <- 0L
+  for (s1 in seq_len(W) - 1L) {
+    for (s2 in s1:(W - 1L)) {
+      k <- k + 1L
+      rows[[k]] <- .ShapePairMC(s1, s2, n, nSim)
+    }
+    if (W > 20 && s1 %% 10L == 0L && s1 > 0L) {
+      cat(sprintf("  s1=%d/%d  k=%d/%d  elapsed=%.1f s\n",
+                  s1, W - 1L, k, pairs_total,
+                  as.numeric(Sys.time() - t0, units = "secs")))
+      flush.console()
+    }
+  }
+  df <- do.call(rbind, rows)
+  cat(sprintf("  done in %.1f s\n",
+              as.numeric(Sys.time() - t0, units = "secs")))
+  df
+}
+
+# (1) Rebuild n=4 against fixed PID.
+result_list[["4"]] <- BuildOneN(4L, 1000L)
+saveRDS(result_list, out_file, compress = TRUE)
+cat("After n=4 refresh, file size:", file.info(out_file)$size, "bytes\n")
+cat("n=4 row check:\n"); print(result_list[["4"]])
+
+# (2) Add n=10.
+if (!"10" %in% names(result_list)) {
+  result_list[["10"]] <- BuildOneN(10L, 500L)
+  result_list <- result_list[order(as.integer(names(result_list)))]
+  saveRDS(result_list, out_file, compress = TRUE)
+}
+
+cat("\nFinal keys:", paste(names(result_list), collapse = ","), "\n")
+cat("Final file size:", file.info(out_file)$size, "bytes\n")

data-raw/extend_shape_lookup.R

@@ -0,0 +1,34 @@
+#!/bin/bash
+#SBATCH --job-name=acid-build
+#SBATCH --partition=shared
+#SBATCH -n 4
+#SBATCH --mem=8G
+#SBATCH --time=00:45:00
+#SBATCH --output=/nobackup/%u/acid-aphid/logs/%x_%j.out
+#SBATCH --error=/nobackup/%u/acid-aphid/logs/%x_%j.err
+
+set -euo pipefail
+module load r/4.5.1
+module load gcc/14.2
+export OMP_NUM_THREADS=1
+export OPENBLAS_NUM_THREADS=1
+
+LIB=/nobackup/$USER/acid-aphid/lib
+export R_LIBS_USER=$LIB
+
+cd /nobackup/$USER/acid-aphid/repos/TreeTools
+git pull origin acid-aphid
+rm -f src/*.o src/*.so
+R CMD build --no-build-vignettes --no-manual --no-resave-data .
+R CMD INSTALL --library=$LIB TreeTools_*.tar.gz
+rm -f TreeTools_*.tar.gz
+
+cd /nobackup/$USER/acid-aphid/repos/TreeDist
+git pull origin acid-aphid
+rm -f src/*.o src/*.so
+R CMD build --no-build-vignettes --no-manual --no-resave-data .
+R CMD INSTALL --library=$LIB TreeDist_*.tar.gz
+rm -f TreeDist_*.tar.gz
+
+echo "Build complete."
+Rscript -e ".libPaths(c('$LIB', .libPaths())); cat('TreeTools', as.character(packageVersion('TreeTools')), '\n'); cat('TreeDist', as.character(packageVersion('TreeDist')), '\n')"

data-raw/hamilton_build.sh

@@ -0,0 +1,62 @@
+# shape_lookup_n.R - compute the shape-conditioned mean CID & PID for every
+# pair of rooted shapes at given n.  Args: n nSim out_path
+suppressPackageStartupMessages({
+  .libPaths(c(Sys.getenv("R_LIBS_USER", "/nobackup/pjjg18/acid-aphid/lib"),
+              .libPaths()))
+  library("TreeTools")
+  library("TreeDist")
+})
+
+args <- commandArgs(trailingOnly = TRUE)
+n    <- as.integer(args[1])
+nSim <- as.integer(args[2])
+out_path <- args[3]
+stopifnot(!is.na(n), !is.na(nSim), nzchar(out_path))
+
+cat(sprintf("[%s] n=%d nSim=%d out=%s\n",
+            format(Sys.time(), "%Y-%m-%d %H:%M:%S"), n, nSim, out_path))
+
+W <- as.integer(NRootedShapes(n))
+pairs_total <- W * (W + 1L) / 2L
+cat(sprintf("W=%d  pairs=%d\n", W, pairs_total))
+
+ShapePairMC <- function(s1, s2, n, nSim) {
+  lab_pool <- as.character(seq_len(n))
+  cid_v <- numeric(nSim); pid_v <- numeric(nSim)
+  for (i in seq_len(nSim)) {
+    t1 <- Preorder(RootedTreeWithShape(s1, n, sample(lab_pool)))
+    t2 <- Preorder(RootedTreeWithShape(s2, n, sample(lab_pool)))
+    cid_v[i] <- ClusteringInfoDistance(t1, t2, normalize = TRUE)
+    pid_v[i] <- PhylogeneticInfoDistance(t1, t2, normalize = TRUE)
+  }
+  c(mean_cid = mean(cid_v), sd_cid = sd(cid_v),
+    mean_pid = mean(pid_v), sd_pid = sd(pid_v))
+}
+
+t0 <- Sys.time()
+rows <- vector("list", pairs_total)
+k <- 0L
+for (s1 in seq_len(W) - 1L) {
+  for (s2 in s1:(W - 1L)) {
+    k <- k + 1L
+    stats <- ShapePairMC(s1, s2, n, nSim)
+    rows[[k]] <- data.frame(shape1 = s1, shape2 = s2,
+                            mean_cid = stats[["mean_cid"]],
+                            sd_cid   = stats[["sd_cid"]],
+                            mean_pid = stats[["mean_pid"]],
+                            sd_pid   = stats[["sd_pid"]],
+                            n_pairs  = nSim)
+  }
+  if (s1 %% 10L == 0L && s1 > 0L) {
+    elapsed <- as.numeric(Sys.time() - t0, units = "secs")
+    eta <- elapsed * pairs_total / k - elapsed
+    cat(sprintf("[%s] s1=%d/%d  k=%d/%d  elapsed=%.1fs  eta=%.1fs\n",
+                format(Sys.time(), "%H:%M:%S"),
+                s1, W - 1L, k, pairs_total, elapsed, eta))
+    flush.console()
+  }
+}
+df <- do.call(rbind, rows)
+saveRDS(df, out_path, compress = TRUE)
+cat(sprintf("Wrote %s (%.1f KB)\n", out_path, file.info(out_path)$size / 1024))
+cat(sprintf("Total wall: %.1f s\n", as.numeric(Sys.time() - t0, units = "secs")))

data-raw/hamilton_shape_lookup.R

@@ -0,0 +1,31 @@
+#!/bin/bash
+#SBATCH --job-name=acid-shape
+#SBATCH --partition=shared
+#SBATCH -n 1
+#SBATCH --mem=4G
+#SBATCH --time=24:00:00
+#SBATCH --array=10-12
+#SBATCH --output=/nobackup/%u/acid-aphid/logs/%x_%A_%a.out
+#SBATCH --error=/nobackup/%u/acid-aphid/logs/%x_%A_%a.err
+
+set -euo pipefail
+module load r/4.5.1
+export OMP_NUM_THREADS=1
+export OPENBLAS_NUM_THREADS=1
+
+LIB=/nobackup/$USER/acid-aphid/lib
+export R_LIBS_USER=$LIB
+
+n=$SLURM_ARRAY_TASK_ID
+
+# nSim scales down as W(n) grows.  Targets: n=10 nSim=500 (~108 min), n=11
+# nSim=200 (~3 h), n=12 nSim=100 (~7-8 h).  All comfortably under 24 h.
+case $n in
+  10) nSim=500 ;;
+  11) nSim=200 ;;
+  12) nSim=100 ;;
+  *)  nSim=200 ;;
+esac
+
+out=/nobackup/$USER/acid-aphid/results/shape_lookup_n${n}.rds
+Rscript /nobackup/$USER/acid-aphid/scripts/shape_lookup_n.R $n $nSim $out

data-raw/hamilton_shape_lookup.sh

@@ -69,15 +69,15 @@ test_that("Lookup agrees with on-the-fly MC at n = 20", {
 })
 
 test_that("Lookup falls back to MC outside the shipped range", {
-  # n = 300 is beyond the shipped table's max (n = 200); the lookup path
-  # should silently fall back to MC with an informative message.
+  # n = 500 is beyond the shipped table's max (n = 400 post-extension); the
+  # lookup path should silently fall back to MC with an informative message.
   expect_message(
-    val <- ExpectedClusteringInfoDistance(300L, nSim = 5L),
+    val <- ExpectedClusteringInfoDistance(500L, nSim = 5L),
     "outside the shipped lookup range"
   )
   expect_true(is.finite(val) && val > 0 && val < 1)
   expect_message(
-    ExpectedPhylogeneticInfoDistance(300L, nSim = 5L),
+    ExpectedPhylogeneticInfoDistance(500L, nSim = 5L),
     "outside the shipped lookup range"
   )
 })