From ea74613cfdd9dde27614ded873adf3f61f1e285f Mon Sep 17 00:00:00 2001 From: R script Date: Thu, 30 Jul 2026 16:22:44 +0100 Subject: [PATCH 01/16] feat(search): deepen perturbation when targetHits is doubled Extends the existing `targetHits` escalation beyond ratchet depth. A caller who at least doubles `targetHits` has asked to keep searching well past ordinary convergence, so under `thorough`/`large` also deepen the per-replicate perturbation itself: drift 25 cycles, the auto/deep reweighting kick, a post-ratchet sectorial pass, and (internally) near-optimal pool retention for intraFuse. The ratio `targetHits / defaultHits` is factored out into `.TargetHitsEscalation()` and `.IwRatchetDepth()` now reads it too, so the two escalations share one signal rather than each deriving its own. Scoped to `thorough`/`large`, matching `.IwRatchetDepth`: `sprint`/`default` document themselves as shallow and have their own implied-weights operating point (`.iwStopPackage`), which this must not disturb. Unlike the ratchet deepening, this applies under any scoring regime. Calibrated on 25 training matrices x 5 seeds = 125 cells, both arms at the same raised `targetHits` so only these levers differ, no cell truncated at its wall cap. Strict paired wins by tier: small 0/0/35 tie, medium 0/0/35 tie, large 1 better/1 worse/33 tie, xlarge 5 better/0 worse/15 tie. All five xlarge wins are one matrix, project4284 at 4062 tips, improved on all five seeds by 1-9 steps; the other xlarge matrices (125/131/173 tips) tied. So the benefit is datasets too large to converge in an ordinary budget, not a size class. Cost: median wall x3.56 and time-to-best x2.49, but replicates-to-best x1.00 -- entirely work per replicate, not slower convergence. That is why these stay out of the presets and behind an explicit signal. `ratchetCycles` is deliberately excluded: depth belongs to `.IwRatchetDepth`, which scales it continuously against a 36-matrix calibration, and a flat value here would clobber it under implied weights. Equal weights loses nothing -- a 68-matrix comparison found no equal-weights reach gain from ratchet depth. `poolSuboptimal` is raised only as an internal aid (intraFuse needs suboptimal recipients). Because `collapse = FALSE` returns the pool verbatim, the return path now drops back to the best score when the escalation raised it, so an escalated search cannot silently return trees worse than `attr(, "score")`. A caller's own `poolSuboptimal` is untouched. Caveat recorded deliberately: the A/B measured seven levers ungated via `strategy = "auto"`; six ship, gated. The evidence for the shipped configuration is a tier decomposition of that run -- the xlarge win is under `large` and preserved, and the excluded tiers measured 0 better / 0 worse, so the gate should remove only cost -- not a direct measurement of it. The A/B also ran under equal weights only, so use under implied weights and profile parsimony is a deliberate but unmeasured extrapolation. Co-Authored-By: Claude Opus 4.8 --- NEWS.md | 19 +++ R/MaximizeParsimony.R | 173 ++++++++++++++++++++-- man/MaximizeParsimony.Rd | 20 ++- tests/testthat/test-reach-escalation.R | 191 +++++++++++++++++++++++++ 4 files changed, 392 insertions(+), 11 deletions(-) create mode 100644 tests/testthat/test-reach-escalation.R diff --git a/NEWS.md b/NEWS.md index 463268405..e9dd9d066 100644 --- a/NEWS.md +++ b/NEWS.md @@ -51,6 +51,25 @@ profile parsimony are unchanged, as is `thorough`/`large`, and setting any of these fields yourself overrides all of it. +- Doubling `targetHits` or more, under `strategy = "thorough"` or `"large"`, now + also deepens the per-replicate perturbation itself, extending the existing + `targetHits` escalation beyond ratchet depth: more drifting, a larger + reweighting kick, a second sectorial pass after the ratchet, and internal + retention of near-optimal trees to fuse against. Unlike the ratchet deepening + this applies under any scoring regime, though it was measured only under equal + weights. It targets datasets big or difficult enough that an ordinary search + stops short of the optimum: across 25 datasets spanning 20 to 4062 tips it + found shorter trees only on the 4062-tip matrix (on all five seeds tried, by + 1–9 steps), while from 20 to 173 tips it found trees of the same length and + simply took about 3.5× as long — a cost incurred as extra work per replicate, + not as slower convergence. Because most searches would pay for depth they do + not need, it is offered only on that explicit signal and only on those two + presets; `sprint` and `default` keep the implied-weights operating point + described above, `ratchetCycles` remains governed by the implied-weights + ratchet deepening, and any control field you set yourself is preserved. Note + that the documented large-`targetHits` idiom for collecting the full set of + most-parsimonious trees also engages this on those two presets. + - Fixed: a large `targetHits` combined with a large `perturbStopFactor` stopped the search after two replicates and silently returned a worse tree. The no-improvement rule computes `(targetHits / hits) * nTip * perturbStopFactor`, diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 4fed58204..2449c02b3 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -208,6 +208,20 @@ # the two in step if this changes. .iwRatchetMaxCycles <- 115L +# How far the caller has raised `targetHits` above its size-scaled default -- the +# single "search harder" signal, shared by every escalation below so that they +# read one quantity rather than each deriving its own. 1 means "not raised"; +# never less, so lowering `targetHits` (e.g. the documented `targetHits = 4` +# idiom) only stops the search sooner and never weakens the search itself. +.TargetHitsEscalation <- function(targetHits, defaultHits) { + if (length(defaultHits) == 1L && is.finite(defaultHits) && defaultHits > 0 && + length(targetHits) == 1L && is.finite(targetHits)) { + max(1, targetHits / defaultHits) + } else { + 1 + } +} + # Ratchet depth to impose for this call, or NULL to leave the preset's value. # `userSet` names the fields the caller set themselves (never overridden). # See the call site in MaximizeParsimony() for the calibration behind it. @@ -225,13 +239,7 @@ if ("ratchetCycles" %in% userSet) { return(NULL) } - escalation <- if (length(defaultHits) == 1L && is.finite(defaultHits) && - defaultHits > 0 && length(targetHits) == 1L && - is.finite(targetHits)) { - max(1, targetHits / defaultHits) - } else { - 1 - } + escalation <- .TargetHitsEscalation(targetHits, defaultHits) min(.iwRatchetMaxCycles, as.integer(round(.iwRatchetCycles * escalation))) } @@ -268,6 +276,88 @@ if (!length(out)) NULL else out } +# Escalation ratio at or above which the deeper-perturbation bundle below +# engages. The ratchet depth above scales continuously because cycle counts +# interpolate; this bundle is mostly switches, which cannot, so it needs one +# trigger point -- and 2 is the ratio that was actually measured (see the +# evidence in `.ReachEscalationDeltas`). +.reachEscalationMinRatio <- 2 + +# Deeper per-replicate perturbation for a caller who has at least doubled +# `targetHits` -- i.e. asked to keep searching well past ordinary convergence. +# +# EVIDENCE (general-pool A/B, 2026-07-28; 25 training matrices x 5 seeds = 125 +# cells; both arms ran at the SAME raised `targetHits`, so only these levers +# differ; no cell in either arm hit its wall cap, so the comparison is not a +# budget artefact). Strict paired final-score wins, by size tier: +# small (n=35): 0 better, 0 worse, 35 tie +# medium (n=35): 0 better, 0 worse, 35 tie +# large (n=35): 1 better, 1 worse, 33 tie (a wash) +# xlarge (n=20): 5 better, 0 worse, 15 tie +# Read that last row carefully: ALL FIVE wins are the SAME matrix, project4284 +# at 4062 tips, which improved on every one of its five seeds (by 1-9 steps). +# The other three xlarge matrices (125, 131 and 173 tips) all tied. So the +# demonstrated benefit is NOT "datasets over 120 tips" -- it is datasets far too +# large to converge within an ordinary budget, plus the hard-reach tail (on the +# 482-tip project5432 these move paired seeds 1946->1945 and 1946->1944, though +# a single call still floors one step above the best known tree). +# Cost: median total wall x3.56 and time-to-best x2.49, but replicates-to-best +# x1.00 -- the wall gap is entirely cost-PER-replicate (~15x candidates +# evaluated), not slower convergence. +# So over the tested 20-173 tip range these buy no reach and cost ~3.5x the wall +# (the sample jumps from 173 tips straight to 4062, so the range in between is +# untested), which is exactly why they are kept OUT of `.StrategyPresets()`: as +# blanket defaults drift especially is per-replicate overhead that, at a fixed +# budget, completes fewer replicates and so reaches the optimum LESS reliably. +# Gating them on a raised `targetHits` puts that cost only where the caller asked +# for it. Note the A/B ran under EQUAL weights, so applying these under implied +# weights and profile parsimony is an extrapolation -- deliberate (the levers are +# scorer-agnostic search machinery) but unmeasured there. +# +# NB `ratchetCycles` is deliberately NOT here. It was part of the tested bundle, +# but ratchet depth is owned by `.IwRatchetDepth` above, which scales it +# continuously against a 36-matrix calibration; a flat value here would clobber +# that under implied weights. Equal weights loses nothing by the omission -- a +# 68-matrix comparison found ratchet depth gives no equal-weights reach gain +# (0.970 vs 0.965), so it is not what the A/B above was measuring. +.ReachEscalationDeltas <- function() { + list( + ratchetPerturbMaxMoves = 0L, # 0 => auto/deep kick + driftCycles = 25L, + postRatchetSectorial = TRUE, + stallEscalateFactor = 1.5, + intraFuse = TRUE, + poolSuboptimal = 3 + ) +} + +# Apply the deeper-perturbation bundle, preserving every field the caller set. +# `userSet` names those fields, exactly as for `.IwRatchetDepth`. +# +# Scoped to `thorough`/`large`, matching `.IwRatchetDepth`: those are the presets +# the A/B's benefit came from (auto selects `large` for the 4062-tip matrix and +# `thorough` for project5432), and on smaller data the same A/B measured 0 better +# / 0 worse across 70 small- and medium-tier cells -- pure wall cost. Escalating +# `sprint`/`default` would also contradict their documented character ("Fast +# search: 3 ratchet cycles, no drift"), so they are left alone. +.ApplyReachEscalation <- function(control, strategy, escalation, + userSet = character(0)) { + if (!length(strategy) || !strategy %in% c("thorough", "large")) { + return(control) + } + if (!length(escalation) || !is.finite(escalation) || + escalation < .reachEscalationMinRatio) { + return(control) + } + deltas <- .ReachEscalationDeltas() + for (nm in names(deltas)) { + if (!(nm %in% userSet)) { + control[[nm]] <- deltas[[nm]] + } + } + control +} + # Strategy presets for adaptive search (Phase 6E). # Wrapped in a function to avoid load-order dependency on SearchControl(). .StrategyPresets <- function() { @@ -675,6 +765,24 @@ #' does not make the ratchet shallower than its default depth (fewer cycles #' were slower to the optimum on every matrix tested), and setting #' `ratchetCycles` yourself overrides this entirely. +#' +#' Doubling `targetHits` or more goes one step further and, under +#' `strategy = "thorough"` or `"large"`, also deepens the per-replicate +#' perturbation itself: more drifting, a larger reweighting kick, a second +#' sectorial pass after the ratchet, and (internally) retention of near-optimal +#' trees to fuse against. Unlike the ratchet deepening above this applies under +#' any scoring regime, though it was measured only under equal weights. It is +#' aimed at datasets big or difficult enough that an ordinary search stops short +#' of the optimum: in testing it found shorter trees only on a 4062-tip matrix +#' (on all five seeds tried), while from 20 to 173 tips it found trees of the +#' same length and simply took around 3.5 times as long -- so it is offered on +#' this explicit signal rather than enabled by default. Any of these values you +#' set yourself is left untouched, and the trees returned are still only the +#' best found. +#' Note that the large-`targetHits` idiom for collecting the full set of +#' most-parsimonious trees, above, therefore also engages this deeper search on +#' those two presets; set `driftCycles`, `intraFuse` and so on yourself if you +#' want the wider sampling without the extra per-replicate cost. #' @param maxSeconds Numeric: maximum wall-clock time in seconds for the #' search. When reached, the current replicate finishes and the search #' stops. `0` (default) means no time limit. @@ -951,6 +1059,8 @@ MaximizeParsimony <- function( } } } + # Set when the reach escalation below raises `poolSuboptimal` itself; see there. + escalatedPool <- FALSE # --- Apply strategy preset --- if (!is.null(strategy) && !identical(strategy, "none")) { @@ -1011,9 +1121,9 @@ MaximizeParsimony <- function( # Scoped to `thorough`/`large`, whose other knobs match the grid; `default` # and `sprint` co-tuned their ratchet with different sectorial settings and # are untouched. + userSet <- union(names(controlDots), attr(control, "explicit")) iwCycles <- .IwRatchetDepth( - strategy, concavity, targetHits, defaultHits, - userSet = union(names(controlDots), attr(control, "explicit")) + strategy, concavity, targetHits, defaultHits, userSet = userSet ) if (!is.null(iwCycles)) { control[["ratchetCycles"]] <- iwCycles @@ -1064,6 +1174,39 @@ MaximizeParsimony <- function( for (.f in names(iwStop)) { control[[.f]] <- iwStop[[.f]] } + + # The same `targetHits` signal, one step further: a caller who has at least + # DOUBLED it has asked to keep searching well past ordinary convergence, so + # also deepen the per-replicate perturbation itself (drift, the auto kick, + # a post-ratchet sectorial re-search, near-optimal pool retention). Unlike + # the ratchet depth above this applies under any scorer, but it is not + # scaled continuously -- it is mostly switches, and 2x is the ratio that was + # measured. See `.ReachEscalationDeltas` for the evidence and the reason + # `ratchetCycles` is left to `.IwRatchetDepth` alone. + # + # Three escalations now read `strategy`, and they are pairwise disjoint by + # design: `.IwRatchetDepth` (ratchet depth, thorough/large, implied weights) + # and `.IwStopPackage` (sprint/default, implied weights) never both fire, + # and this bundle is scoped to thorough/large so it cannot disturb the + # sprint/default operating point measured above. + escalation <- .TargetHitsEscalation(targetHits, defaultHits) + poolBefore <- control[["poolSuboptimal"]] + control <- .ApplyReachEscalation(control, strategy, escalation, + userSet = userSet) + # `poolSuboptimal` is raised here only as an internal aid: `intraFuse` needs + # suboptimal recipients to fuse against. It must not leak into the RESULT -- + # `collapse = FALSE` returns the pool verbatim, so without this the caller + # would silently receive trees up to 3 steps worse than `attr(, "score")` + # from a result documented as the best trees found. Recorded here (rather + # than compared later) so a caller's own `poolSuboptimal` is untouched: they + # asked for those trees and still get them. + escalatedPool <- !identical(control[["poolSuboptimal"]], poolBefore) + if (verbosity >= 1L && escalation >= .reachEscalationMinRatio && + strategy %in% c("thorough", "large")) { + cli::cli_alert_info( + "Deep search: {.field targetHits} raised {round(escalation, 1)}x" + ) + } } else if (!identical(strategy, "auto")) { warning("Unknown strategy '", strategy, "'; using default parameters.") } @@ -1497,7 +1640,17 @@ MaximizeParsimony <- function( }) nTopologies <- collapsed$n_topologies } else { - outTrees <- lapply(resultTrees, function(edgeMat) { + # The pool is returned verbatim here, suboptimal entries included -- which is + # what a caller who set `poolSuboptimal` themselves asked for. But when the + # reach escalation raised it on their behalf (an internal aid for `intraFuse`, + # see .ApplyReachEscalation), they did not: drop back to the best score so the + # result matches its documentation, "the best tree(s) found". + keepTrees <- if (escalatedPool) { + resultTrees[result$scores == result$best_score] + } else { + resultTrees + } + outTrees <- lapply(keepTrees, function(edgeMat) { tr <- treeTpl tr[["edge"]] <- edgeMat # C++ edge order may differ from template; renumber to valid preorder diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 62d7c21e4..127dfdfe3 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -224,7 +224,25 @@ much character reweighting a matrix needs, so a raised \code{targetHits} is take as the user's own signal that this one needs more. Lowering \code{targetHits} does not make the ratchet shallower than its default depth (fewer cycles were slower to the optimum on every matrix tested), and setting -\code{ratchetCycles} yourself overrides this entirely.} +\code{ratchetCycles} yourself overrides this entirely. + +Doubling \code{targetHits} or more goes one step further and, under +\code{strategy = "thorough"} or \code{"large"}, also deepens the per-replicate +perturbation itself: more drifting, a larger reweighting kick, a second +sectorial pass after the ratchet, and (internally) retention of near-optimal +trees to fuse against. Unlike the ratchet deepening above this applies under +any scoring regime, though it was measured only under equal weights. It is +aimed at datasets big or difficult enough that an ordinary search stops short +of the optimum: in testing it found shorter trees only on a 4062-tip matrix +(on all five seeds tried), while from 20 to 173 tips it found trees of the +same length and simply took around 3.5 times as long -- so it is offered on +this explicit signal rather than enabled by default. Any of these values you +set yourself is left untouched, and the trees returned are still only the +best found. +Note that the large-\code{targetHits} idiom for collecting the full set of +most-parsimonious trees, above, therefore also engages this deeper search on +those two presets; set \code{driftCycles}, \code{intraFuse} and so on yourself if you +want the wider sampling without the extra per-replicate cost.} \item{maxSeconds}{Numeric: maximum wall-clock time in seconds for the search. When reached, the current replicate finishes and the search diff --git a/tests/testthat/test-reach-escalation.R b/tests/testthat/test-reach-escalation.R new file mode 100644 index 000000000..176bd0b6a --- /dev/null +++ b/tests/testthat/test-reach-escalation.R @@ -0,0 +1,191 @@ +# Deep-search escalation: raising `targetHits` deepens per-replicate perturbation. +library("TreeTools", quietly = TRUE) +data("inapplicable.phyData", package = "TreeSearch") +ds <- inapplicable.phyData[["Vinther2008"]] # 23 tips + +test_that(".TargetHitsEscalation reports the raise ratio, floored at 1", { + TE <- TreeSearch:::.TargetHitsEscalation + expect_equal(TE(10L, 10L), 1) # at the default + expect_equal(TE(20L, 10L), 2) # doubled + expect_equal(TE(96L, 96L), 1) + expect_equal(TE(192L, 96L), 2) + expect_equal(TE(48L, 96L), 1) # LOWERED -> never below 1 + expect_equal(TE(4L, 10L), 1) # the documented "one tree, quickly" idiom + # Degenerate inputs fall back to "not raised" rather than erroring. + expect_equal(TE(NA_integer_, 10L), 1) + expect_equal(TE(10L, 0L), 1) + expect_equal(TE(Inf, 10L), 1) +}) + +test_that(".ReachEscalationDeltas is the exact expected lever set", { + expect_identical( + TreeSearch:::.ReachEscalationDeltas(), + list( + ratchetPerturbMaxMoves = 0L, + driftCycles = 25L, + postRatchetSectorial = TRUE, + stallEscalateFactor = 1.5, + intraFuse = TRUE, + poolSuboptimal = 3 + ) + ) +}) + +test_that("every escalation delta is a real SearchControl field", { + expect_true(all(names(TreeSearch:::.ReachEscalationDeltas()) %in% + names(SearchControl()))) +}) + +test_that("escalation does NOT touch ratchetCycles", { + # Ratchet depth belongs to .IwRatchetDepth, which scales it continuously + # against a 36-matrix calibration (48 cycles, up to 115). A flat value here + # would silently clobber that under implied weights, so the bundle must never + # carry ratchetCycles. + expect_false("ratchetCycles" %in% names(TreeSearch:::.ReachEscalationDeltas())) + ctrl <- TreeSearch:::.ApplyReachEscalation(SearchControl(ratchetCycles = 48L), + "thorough", escalation = 10) + expect_identical(ctrl[["ratchetCycles"]], 48L) +}) + +test_that("the two escalations compose without fighting over ratchetCycles", { + # Mirrors the call-site order: .IwRatchetDepth sets the depth, then the bundle + # is applied on top. The depth must survive, at its calibrated value. + ctrl <- TreeSearch:::.ApplyStrategyPreset( + SearchControl(), TreeSearch:::.StrategyPresets()[["thorough"]] + ) + iw <- TreeSearch:::.IwRatchetDepth("thorough", concavity = 10, + targetHits = 20L, defaultHits = 10L) + expect_equal(iw, 96L) # 48 * escalation(2), under the 115 cap + ctrl[["ratchetCycles"]] <- iw + ctrl <- TreeSearch:::.ApplyReachEscalation(ctrl, "thorough", escalation = 2) + expect_identical(ctrl[["ratchetCycles"]], 96L) # NOT clobbered by the bundle + expect_identical(ctrl[["driftCycles"]], 25L) # bundle still applied +}) + +test_that(".ApplyReachEscalation applies all deltas at or above the ratio", { + deltas <- TreeSearch:::.ReachEscalationDeltas() + for (strat in c("thorough", "large")) { + for (esc in c(2, 2.5, 20)) { + ctrl <- TreeSearch:::.ApplyReachEscalation(SearchControl(), strat, + escalation = esc) + for (nm in names(deltas)) { + expect_identical(ctrl[[nm]], deltas[[nm]], info = paste(strat, nm)) + } + } + } +}) + +test_that(".ApplyReachEscalation is inert below the ratio", { + stock <- SearchControl() + for (esc in c(1, 1.5, 1.99)) { + expect_identical(TreeSearch:::.ApplyReachEscalation(stock, "thorough", + escalation = esc), + stock) + } + # Degenerate escalation must not escalate. + for (esc in list(NA_real_, numeric(0), Inf)) { + expect_identical(TreeSearch:::.ApplyReachEscalation(stock, "thorough", + escalation = esc), + stock) + } +}) + +test_that(".ApplyReachEscalation is scoped to thorough/large", { + # `sprint` and `default` document themselves as fast/shallow ("3 ratchet + # cycles, no drift"), and the A/B measured 0 better / 0 worse across their + # 70 small- and medium-tier cells -- pure wall cost. They must not escalate. + stock <- SearchControl() + for (strat in c("sprint", "default", "none", NA_character_, character(0))) { + expect_identical( + TreeSearch:::.ApplyReachEscalation(stock, strat, escalation = 10), + stock, info = paste("strategy", strat) + ) + } +}) + +test_that(".ApplyReachEscalation preserves caller-set fields", { + ctrl <- TreeSearch:::.ApplyReachEscalation( + SearchControl(), "thorough", escalation = 4, + userSet = c("driftCycles", "intraFuse") + ) + expect_identical(ctrl[["driftCycles"]], SearchControl()[["driftCycles"]]) + expect_identical(ctrl[["intraFuse"]], SearchControl()[["intraFuse"]]) + # a field the caller did not set still escalates + expect_identical(ctrl[["postRatchetSectorial"]], TRUE) +}) + +test_that("the default targetHits never escalates", { + # Off-by-default guarantee: the size-scaled default is its own reference, so an + # un-tuned search sits at ratio 1 -- below the trigger -- for any tip count. + TE <- TreeSearch:::.TargetHitsEscalation + for (n in c(5L, 23L, 50L, 88L, 200L, 482L, 4062L)) { + defHits <- max(10L, as.integer(n / 5)) + expect_lt(TE(defHits, defHits), TreeSearch:::.reachEscalationMinRatio) + } +}) + +test_that("escalation measurably deepens the search end to end", { + # Not a smoke test: asserts a signal that DISAPPEARS if the feature is removed. + # The bundle multiplies per-replicate work (drift 2 -> 25, deep kick, an extra + # sectorial pass), so at matched replicates the escalated run must evaluate + # substantially more candidates. Vinther2008 (23 tips): default targetHits = 10, + # so 20 is exactly 2x. + runCand <- function(hits) { + set.seed(4242) + r <- MaximizeParsimony(ds, strategy = "thorough", maxReplicates = 2L, + targetHits = hits, maxSeconds = 0, verbosity = 0L) + list(cand = as.double(attr(r, "candidates_evaluated")), res = r) + } + ordinary <- runCand(10L) + escalated <- runCand(20L) + expect_s3_class(escalated$res, "multiPhylo") + expect_true(is.finite(attr(escalated$res, "score"))) + expect_true(is.finite(ordinary$cand) && ordinary$cand > 0) + # Same seed and same replicate count, so with the feature removed the two runs + # would be identical and the ratio exactly 1 (as the `sprint` test below shows). + # Measured here at ~1.49; 1.2 leaves room for stochastic drift while still + # failing outright if the escalation stops engaging. NB the ratio is much + # larger against a shallower preset -- `thorough` already drifts, so the + # increment over it is smaller than over `sprint`. + expect_gt(escalated$cand, 1.2 * ordinary$cand) +}) + +test_that("sprint is NOT escalated end to end", { + # The scope guard, observably. Under a preset outside the escalation's scope, + # raising targetHits cannot change per-replicate work at all: with the replicate + # count fixed the two runs are identical, so the candidate counts must match + # EXACTLY. This fails the moment the strategy gate is loosened. + runCand <- function(hits) { + set.seed(99L) + r <- MaximizeParsimony(ds, strategy = "sprint", maxReplicates = 2L, + targetHits = hits, maxSeconds = 0, verbosity = 0L) + as.double(attr(r, "candidates_evaluated")) + } + expect_identical(runCand(20L), runCand(10L)) +}) + +test_that("an escalated search still returns only best-score trees", { + # The bundle raises `poolSuboptimal` internally (intraFuse needs recipients). + # That must not leak into the result: with collapse = FALSE the pool is returned + # verbatim, so without the guard the caller would silently get trees up to 3 + # steps worse than attr(, "score") from a result documented as the best found. + set.seed(31L) + r <- MaximizeParsimony(ds, strategy = "thorough", maxReplicates = 3L, + targetHits = 20L, collapse = FALSE, verbosity = 0L) + best <- attr(r, "score") + expect_true(is.finite(best)) + scores <- vapply(r, function(t) TreeLength(t, ds, concavity = Inf), double(1)) + expect_true(all(scores == best)) +}) + +test_that("a caller's own poolSuboptimal is still honoured", { + # The guard above must not steal the documented behaviour from someone who + # asked for suboptimal trees themselves. + set.seed(31L) + r <- MaximizeParsimony(ds, strategy = "thorough", maxReplicates = 3L, + targetHits = 20L, poolSuboptimal = 3, + collapse = FALSE, verbosity = 0L) + expect_s3_class(r, "multiPhylo") + scores <- vapply(r, function(t) TreeLength(t, ds, concavity = Inf), double(1)) + expect_true(all(scores <= attr(r, "score") + 3)) +}) From 59e8c32ee0ac504257aa6f5defe45efe6cf99a62 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 31 Jul 2026 08:09:38 +0100 Subject: [PATCH 02/16] docs(benchmarks): record the targetHits-escalation ship gate and reach study Commits the A/B harness and its findings so the evidence behind the escalation is reproducible rather than living in scratch. Raw per-cell CSVs stay gitignored, as elsewhere in dev/benchmarks; the markdown is the durable record. The findings state plainly what the study does and does not support: the general-pool win is one 4062-tip matrix (all five of its seeds, 1-9 steps), not a property of the >=121-tip tier, and the 20-173 tip range measured 0 better / 0 worse at ~3.5x the wall -- which is the argument for gating rather than defaulting. Also records that the shipped six-lever gated configuration was never A/B'd as such (the run tested seven, ungated), and that the whole study was equal-weights only. The hard-tail section records the vehicle comparison for a future reach lever: an external block loop (5/16) and the in-engine TS_POOL_RESEED (4/16) are statistically indistinguishable, reaching the floor ~25-31% of the time either way, and a supplied starting tree seeds only replicate 0 -- so re-solving the incumbent is a restart-volume phenomenon that outerCycles cannot substitute for. Includes four methodological traps that cost time here: maxSeconds is not the search budget (enumTimeFraction reserves 10%, and identical elapsed times across seeds is the truncation signature); whole-suite runs need test_local() because several files call internals unqualified; a top-level skip_on_cran() reports "0 pass 0 fail" which is not a pass; and attributing a suite failure needs a pristine worktree at the parent commit, not a single-file swap. Co-Authored-By: Claude Opus 4.8 --- dev/benchmarks/reach_escalation_FINDINGS.md | 129 ++++++++++++ dev/benchmarks/reach_escalation_ab.R | 222 ++++++++++++++++++++ dev/benchmarks/reach_escalation_ab.sh | 21 ++ dev/benchmarks/reach_escalation_analyze.R | 137 ++++++++++++ 4 files changed, 509 insertions(+) create mode 100644 dev/benchmarks/reach_escalation_FINDINGS.md create mode 100644 dev/benchmarks/reach_escalation_ab.R create mode 100644 dev/benchmarks/reach_escalation_ab.sh create mode 100644 dev/benchmarks/reach_escalation_analyze.R diff --git a/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md new file mode 100644 index 000000000..ddbffb5b0 --- /dev/null +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -0,0 +1,129 @@ +# `targetHits` reach escalation — ship gate and reach study (2026-07-24 … 07-31) + +Durable record for the escalation shipped in `MaximizeParsimony()`: doubling `targetHits` +or more, under `strategy = "thorough"`/`"large"`, additionally deepens the per-replicate +perturbation. Harness: `reach_escalation_ab.R` + `reach_escalation_analyze.R` + +`reach_escalation_ab.sh` (this directory). Raw per-cell CSVs are gitignored; the numbers +below and the committed harness are the reproducible record. + +## The levers + +Six `SearchControl` fields, applied as a flat bundle when +`targetHits / max(10, nTip/5) >= 2`, skipping anything the caller set: + + ratchetPerturbMaxMoves = 0 # auto/deep reweighting kick + driftCycles = 25 + postRatchetSectorial = TRUE + stallEscalateFactor = 1.5 + intraFuse = TRUE + poolSuboptimal = 3 # internal only; filtered out of the returned trees + +`ratchetCycles` is deliberately NOT here — ratchet depth belongs to `.IwRatchetDepth()`, +which scales it continuously against its own 36-matrix calibration. A flat value here +would clobber that under implied weights, and a separate 68-matrix comparison found ratchet +depth gives no equal-weights reach gain (0.970 vs 0.965), so equal weights loses nothing. + +## Ship gate — general-pool A/B + +**Design.** `MBANK_FIXED_SAMPLE` (25 training matrices, 20–4062 tips) × 5 seeds = 125 +cells, one SLURM task per cell with both arms on the same node. Validation split +sequestered (asserted per matrix). Regime EW Fitch, gaps→missing. `strategy = "auto"`, +so the escalation layers on whatever preset auto picks — the faithful test. + +**The confound this avoids.** Comparing default-`targetHits` against raised-`targetHits` +would confound the levers with "raised `targetHits` searches longer anyway". So *both* +arms run at the SAME raised `targetHits` (exactly the trigger threshold) and differ ONLY in +the six levers. Run on a gate-free engine with the levers passed as dots, so the result +does not depend on the gate implementation being correct. + +**Target.** Union-best final score across arms within each cell (the established mbank +convention — there is no canonical best-known table for these matrices). Note this is +self-referential: if one arm alone attains a score, the other "misses" by construction, so +the reach fractions restate the paired counts rather than measuring absolute optimality. +The paired win counts below are the honest statistic. + +**Result.** No cell in either arm hit its wall cap, so nothing is a budget artefact. + +Strict paired final-score wins (deltas vs base): + +| tier | n | deltas better | base better | tie | +|------|---|---------------|-------------|-----| +| small (≤30 tips) | 35 | 0 | 0 | 35 | +| medium (31–60) | 35 | 0 | 0 | 35 | +| large (61–120) | 35 | 1 | 1 | 33 | +| xlarge (≥121) | 20 | 5 | 0 | 15 | + +**Read the xlarge row carefully: all five wins are the same matrix**, `project4284` at 4062 +tips, which improved on every one of its five seeds (by 1–9 steps). The other three xlarge +matrices (125, 131, 173 tips) all tied. So the demonstrated benefit is *datasets far too +large to converge within an ordinary budget*, **not** a property of "over 120 tips". + +Cost, median over cells: total wall ×3.56, time-to-best ×2.49, but **replicates-to-best +×1.00**. The wall gap is entirely work per replicate (~15× candidates evaluated against a +`sprint` baseline; ~1.5× against `thorough`), not slower convergence. Reporting wall alone +would read as a 2.5× regression; `rep2hit` is what shows it is not. (Same lesson as +`kick_anytime_FINDINGS.md`.) + +The one loss: `project2771` (large, seed 5821), base 911 → deltas 912. Against 6 better / +1 worse overall and a tied `large` tier, read as stochastic — but recorded, not swept away. + +**Pre-registered rule** (fixed before results): ship iff reach does not regress and no tier +regresses; wall cost is accepted by construction, since the escalation only fires when the +user has asked for it. → **SHIP.** + +## Why it is gated, and gated to `thorough`/`large` + +Below ~200 tips the levers buy no reach and cost ~3.5× the wall, so they must not be +defaults. Gating on `targetHits` puts the cost only where the user asked for it. Scoping +to `thorough`/`large` matches `.IwRatchetDepth`, keeps `sprint`/`default` at their own +measured implied-weights operating point (`.iwStopPackage`), and costs nothing measured — +the small and medium tiers were 0 better / 0 worse across 70 cells. + +## Hard-tail reach study (project5432, 482 tips, EW; the levers' origin) + +Best known 1943 (shared TS/TNT floor; a 1942 tree exists from a rare TNT run). All arms +`targetHits = 999`, `maxReplicates = 500`, matched seeds, current engine. + +| arm | configuration | reach 1943 | +|-----|---------------|------------| +| A | stock `thorough`, single call | 0/3 (floor 1944) | +| B | + the levers, single call | 0/3 (floor 1944) | +| C | + the levers, R block loop (`tree = best` carry-forward, ~48 re-entries) | **5/16** | +| D | + the levers, `TS_POOL_RESEED=0.5`, single call — **time-truncated** | 2/16 | +| E | as D, un-starved (`enumTimeFraction = 0`, 68 h) | **4/16** | + +**Mechanism.** A supplied `tree =` seeds only replicate 0 (`src/ts_driven.cpp`), and +nothing inside a single call ever re-optimises the global best through the deep pipeline +again — the restart strategies are all from-scratch. So a single call re-solves the +incumbent once; the block loop does so once per block. Reaching the floor is a +restart-volume phenomenon, and `outerCycles`/`maxOuterResets` are not a substitute (they +were identical in every arm). + +**Vehicle verdict.** Arm C 5/16 vs arm E 4/16 is statistically indistinguishable (Fisher +p ≈ 1.0), so the in-engine `TS_POOL_RESEED` matches an external block loop while keeping +the pool, conflict table and fuse donors that each re-entry discards, and paying the +MPT-enumeration reserve once instead of per block. **Reach is ~25–31% on any vehicle: the +recipe improves the odds, it does not guarantee the floor.** + +## Methodological traps worth not re-learning + +- **`maxSeconds` is not the search budget.** `main_deadline = maxSeconds × (1 − + enumTimeFraction)`, default 10% reserve. Identical elapsed times across seeds is the + time-truncation signature; arm B and arm D were truncated exactly there while appearing + replicate-capped. For pure reach runs set `enumTimeFraction = 0`. +- **Whole-suite test runs need `test_local()`/`devtools::test()`**, not `test_dir()` + + `library()`: several test files call internals unqualified and error otherwise. +- **A top-level `skip_on_cran()` makes a file report "0 pass 0 fail"** — that is *not* a + pass, it never ran. Use `NOT_CRAN=true`. +- **To attribute a suite failure, use a pristine detached worktree at the parent commit.** + Swapping a single source file is invalid once you have added a test file, because your own + tests then error on the missing functions and the totals stop being comparable. + +## Not measured + +- The shipped configuration was never A/B'd *as such*: the run above tested seven levers + ungated; six ship, gated to `thorough`/`large`. The inference is a tier decomposition of + that run (the xlarge win is under `large` and preserved; the excluded tiers measured 0 + better / 0 worse, so the gate should remove only cost), not a direct measurement. +- Equal weights only. Use under implied weights and profile parsimony is a deliberate but + unmeasured extrapolation; ratchet depth there is governed separately. diff --git a/dev/benchmarks/reach_escalation_ab.R b/dev/benchmarks/reach_escalation_ab.R new file mode 100644 index 000000000..a0268c90a --- /dev/null +++ b/dev/benchmarks/reach_escalation_ab.R @@ -0,0 +1,222 @@ +#!/usr/bin/env Rscript +# v1 REACH-ESCALATION GATE -- general-pool anytime A/B (the SHIP/NO-SHIP gate). +# +# WHAT IS BEING DECIDED: `MaximizeParsimony()` gains a gate that, when the user raises +# `targetHits` to >= 2 * max(10, nTip/5), applies 7 deeper per-replicate perturbation +# deltas (ratchetCycles 40, kick 0/auto-deep, driftCycles 25, postRatchetSectorial, +# stallEscalateFactor 1.5, intraFuse, poolSuboptimal 3). On the hard-reach tail +# (project5432) these measurably improve the score distribution. This A/B asks the +# OTHER question: when the gate fires on ORDINARY data, does it help or HURT? +# +# DESIGN (advisor-reviewed): +# * BUILDLESS + GATE-FREE ENGINE. Runs on plain cpp-search (curlib) and passes the +# deltas as top-level dots -- byte-equivalent to what the gate applies, so no gated +# build is needed and the result is not contingent on my patch being correct. +# * THE CONFOUND THIS AVOIDS: comparing default-targetHits vs raised-targetHits would +# confound the deltas with "raised targetHits searches longer anyway". So BOTH arms +# use the SAME raised targetHits (exactly the gate threshold) and differ ONLY in the +# 7 deltas. This isolates the escalation. +# base = strategy auto, targetHits = 2*max(10, nTip/5), stock preset +# deltas = identical + the 7 escalation deltas <- what the gate does when it fires +# * strategy = "auto" (NOT pinned): the gate layers on whatever preset auto picks, so +# auto is the faithful test -- and sprint/default (driftCycles 0-2) is exactly where +# driftCycles=25 is the largest relative change, i.e. where a regression would show. +# * METRIC = per-replicate improvement trace via progressCallback (the step function). +# Report REACH + tt_hit (wall) + **rep2hit (replicates)**. rep2hit is mandatory: the +# 2026-07-14 kick study's apparent 7-12% wall2hit "regression" was overturned by +# rep2hit = 1.000 -- same replicate to optimum, just costlier per replicate. Reporting +# wall alone manufactures a false regression verdict. +# * TARGET = union-best final score across arms within each (matrix, seed) cell -- the +# established mbank convention (there is NO canonical best-known table for mbank). +# * BUDGET: tier caps 3x the kick_anytime caps. driftCycles 2->25 is ~10x the drift +# work; reusing the old caps would starve the deltas arm and guarantee a spurious +# regression -- exactly the trap that time-truncated 5432 arm B (it stopped at +# maxSeconds*(1-enumTimeFraction), never spending its replicate budget). +# * VALIDATION SPLIT SEQUESTERED: asserts split == "training" (project_id %% 5 == 0 is +# validation and is a one-way door). +# * REGIME = EW Fitch, gaps->missing. +# +# PRE-REGISTERED DECISION RULE (fixed BEFORE any result was seen): +# The gate is OPT-IN and fires only when the user has said "time is no issue" by +# raising targetHits. So a wall cost is ACCEPTED by construction; what is NOT +# acceptable is the escalation finding WORSE trees. +# SHIP : reach(deltas) >= reach(base) overall AND no size tier shows a reach +# regression. Wall cost reported honestly alongside. +# NO-SHIP : reach(deltas) < reach(base) on the general pool (the deltas actively +# hurt when they fire) -> restrict the gate or drop v1. +# Evidence format = paired per-cell better/worse/tie counts + median paired ratios +# (the repo convention), NOT p-values. +# +# One SLURM array task = one (matrix, seed) cell, both arms on the same node (fair +# same-CPU wall comparison). TASK_ID (1-based) selects the manifest row. +# +# Env: TS_LIB, NEOTRANS_DIR, CAT_CSV, OUT_DIR, TASK_ID / SLURM_ARRAY_TASK_ID, +# N_SEEDS (default 5), TS_MAXREP (default 300), SMOKE (1 matrix, tiny budget). + +suppressMessages({ + ts_lib <- Sys.getenv("TS_LIB", "") + if (nzchar(ts_lib)) { + library(TreeSearch, lib.loc = normalizePath(ts_lib, winslash = "/", mustWork = TRUE)) + } else { + library(TreeSearch) + } + library(TreeTools) +}) + +neo_dir <- Sys.getenv("NEOTRANS_DIR", "") +cat_csv <- Sys.getenv("CAT_CSV", "") +out_dir <- Sys.getenv("OUT_DIR", ".") +if (!nzchar(neo_dir)) stop("NEOTRANS_DIR unset") +if (!nzchar(cat_csv)) stop("CAT_CSV unset") +dir.create(out_dir, showWarnings = FALSE, recursive = TRUE) + +# Fixed 25-matrix training sample (bench_datasets.R MBANK_FIXED_SAMPLE). +# "Do not modify: results are only comparable when the same sample is used." +MBANK_FIXED_SAMPLE <- c( + "project532", "project2346", "project2451", "project4501", + "project944", "project971_(1)", "project2762", + "project826", "project561", "project571", "project4146_(3)", + "project3688", "project4049", "project423", + "project4286", "project4359", "project4397", "project2084_(1)", + "project2771", "project2184", "project3938", + "syab07201", "project4133", "project804", "project4284" +) + +catalogue <- read.csv(cat_csv, stringsAsFactors = FALSE) +rownames(catalogue) <- catalogue$key + +to_fitch <- function(pd) { + m <- PhyDatToMatrix(pd, ambigNA = FALSE) + m[m == "-"] <- "?" + MatrixToPhyDat(m) +} + +load_matrix <- function(key) { + if (!key %in% catalogue$key) stop("key not in catalogue: ", key) + row <- catalogue[key, ] + # SEQUESTER: refuse validation-split matrices. + if (!identical(row$split, "training")) + stop(sprintf("key %s is split='%s' -- validation is SEQUESTERED", key, row$split)) + f <- file.path(neo_dir, row$filename) + if (!file.exists(f)) stop("matrix file not found: ", f) + pd <- suppressWarnings(TreeTools::ReadAsPhyDat(f)) + to_fitch(pd) +} + +`%||%` <- function(a, b) if (is.null(a)) b else a + +# The 7 escalation deltas -- must stay byte-identical to .ReachEscalationDeltas(). +ESCALATION_DELTAS <- list( + ratchetCycles = 40L, ratchetPerturbMaxMoves = 0L, driftCycles = 25L, + postRatchetSectorial = TRUE, stallEscalateFactor = 1.5, + intraFuse = TRUE, poolSuboptimal = 3 +) +# The gate's trigger, mirroring .ReachEscalationThreshold(nTip). +reach_threshold <- function(nTip) 2L * max(10L, as.integer(nTip / 5)) +ARMS <- c("base", "deltas") + +# ---- Anytime tracer: improvement events only (the step function) ---- +make_tracer <- function(t0) { + env <- new.env(parent = emptyenv()) + env$prev <- Inf + env$rows <- list() + cb <- function(info) { + if (!identical(info$phase, "replicate")) return(invisible()) + bs <- info$best_score + if (is.null(bs) || length(bs) != 1L || !is.finite(bs)) return(invisible()) + if (bs < env$prev - 1e-9) { + env$prev <- bs + env$rows[[length(env$rows) + 1L]] <- data.frame( + replicate = as.integer(info$replicate %||% NA_integer_), + elapsed_s = as.double(proc.time()["elapsed"] - t0), + engine_elapsed = if (!is.null(info$elapsed)) as.double(info$elapsed) else NA_real_, + best_score = as.double(bs), stringsAsFactors = FALSE) + } + invisible() + } + list(cb = cb, env = env) +} + +run_arm <- function(pd, nTip, arm, seed, maxrep, cap_s) { + set.seed(seed) + t0 <- proc.time()["elapsed"] + tr <- make_tracer(t0) + args <- list(pd, strategy = "auto", + maxReplicates = maxrep, maxSeconds = cap_s, + targetHits = reach_threshold(nTip), # SAME in both arms + nThreads = 1L, verbosity = 0L, + progressCallback = tr$cb) + if (identical(arm, "deltas")) args <- c(args, ESCALATION_DELTAS) + res <- suppressWarnings(do.call(MaximizeParsimony, args)) + wall <- as.double(proc.time()["elapsed"] - t0) + final_score <- as.double(attr(res, "score")) + reps <- attr(res, "replicates"); reps <- if (is.null(reps)) NA_integer_ else as.integer(reps) + cand <- attr(res, "candidates_evaluated"); cand <- if (is.null(cand)) NA_real_ else as.double(cand) + trace <- if (length(tr$env$rows)) do.call(rbind, tr$env$rows) else + data.frame(replicate = NA_integer_, elapsed_s = NA_real_, + engine_elapsed = NA_real_, best_score = final_score) + list(wall = wall, final_score = final_score, reps = reps, cand = cand, + trace = trace, n_events = length(tr$env$rows)) +} + +N_SEEDS <- as.integer(Sys.getenv("N_SEEDS", "5")) +BASE_SEED <- 5821L +SMOKE <- nzchar(Sys.getenv("SMOKE", "")) +maxrep <- as.integer(Sys.getenv("TS_MAXREP", if (SMOKE) "3" else "300")) + +keys <- MBANK_FIXED_SAMPLE +if (SMOKE) { keys <- keys[1]; N_SEEDS <- 1L } + +manifest <- expand.grid(key = keys, seed_idx = seq_len(N_SEEDS), + stringsAsFactors = FALSE) +manifest <- manifest[order(manifest$key, manifest$seed_idx), ] +rownames(manifest) <- NULL + +tid <- as.integer(Sys.getenv("TASK_ID", Sys.getenv("SLURM_ARRAY_TASK_ID", "1"))) +if (SMOKE) tid <- 1L +if (tid < 1 || tid > nrow(manifest)) + stop(sprintf("TASK_ID %d out of range 1..%d", tid, nrow(manifest))) + +key <- manifest$key[tid] +seed <- BASE_SEED + manifest$seed_idx[tid] - 1L +row <- catalogue[key, ] +nTip <- as.integer(row$ntax) +nChar <- as.integer(row$nchar) +tier <- cut(nTip, breaks = c(0, 30, 60, 120, Inf), + labels = c("small", "medium", "large", "xlarge")) +# 3x the kick_anytime caps: driftCycles 2->25 is ~10x the drift work, so the old caps +# would starve the deltas arm and manufacture a spurious regression. +cap_s <- if (SMOKE) 20 else switch(as.character(tier), + small = 180, medium = 360, large = 720, xlarge = 1440) + +cat(sprintf("=== reach_ab task %d/%d: %s (%dt, %dc, %s) seed=%d targetHits=%d cap=%gs maxrep=%d ===\n", + tid, nrow(manifest), key, nTip, nChar, tier, seed, + reach_threshold(nTip), cap_s, maxrep)) + +pd <- load_matrix(key) +stopifnot(length(pd) == nTip) + +all_rows <- list() +for (arm in ARMS) { + r <- run_arm(pd, nTip, arm, seed, maxrep, cap_s) + cat(sprintf(" %-6s final=%.0f reps=%s wall=%.1fs events=%d candM=%.2f\n", + arm, r$final_score, + ifelse(is.na(r$reps), "?", as.character(r$reps)), + r$wall, r$n_events, r$cand / 1e6)) + # SELF-CHECK: a completed search MUST emit >=1 improvement event (rep 1 from Inf). + if (r$n_events == 0L) + cat(sprintf(" WARN %s: progressCallback emitted 0 improvement events -- trace facility may be broken!\n", arm)) + tr <- r$trace + all_rows[[length(all_rows) + 1L]] <- data.frame( + dataset = key, nTip = nTip, nChar = nChar, tier = as.character(tier), + seed = seed, arm = arm, target_hits = reach_threshold(nTip), + event = "improve", replicate = tr$replicate, elapsed_s = tr$elapsed_s, + engine_elapsed = tr$engine_elapsed, best_score = tr$best_score, + final_score = r$final_score, reps_done = r$reps, wall_total_s = r$wall, + candidates = r$cand, cap_s = cap_s, stringsAsFactors = FALSE) +} + +D <- do.call(rbind, all_rows) +of <- file.path(out_dir, sprintf("cell_%03d_%s_s%d.csv", tid, gsub("[^A-Za-z0-9]", "", key), seed)) +write.csv(D, of, row.names = FALSE) +cat(sprintf("Wrote %s (%d rows)\n", of, nrow(D))) diff --git a/dev/benchmarks/reach_escalation_ab.sh b/dev/benchmarks/reach_escalation_ab.sh new file mode 100644 index 000000000..9ea07cb98 --- /dev/null +++ b/dev/benchmarks/reach_escalation_ab.sh @@ -0,0 +1,21 @@ +#!/bin/bash +#SBATCH --job-name=reach-ab +#SBATCH -p shared +#SBATCH -n 1 +#SBATCH --mem=8G +#SBATCH --time=2:00:00 +#SBATCH --array=1-125 +#SBATCH -o /nobackup/pjjg18/reach/logs/reachab_%A_%a.out +#SBATCH -e /nobackup/pjjg18/reach/logs/reachab_%A_%a.err +module load r/4.5.1 +export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 +# curlib = TreeSearch@cpp-search (49c7fcf, GATE-FREE) + TreeTools@head; deps from TreeSearch/lib. +# The deltas are passed as dots, so no gated build is needed. +export R_LIBS="/nobackup/pjjg18/curlib:/nobackup/pjjg18/TreeSearch/lib" +TS_LIB="/nobackup/pjjg18/curlib" \ +NEOTRANS_DIR="/nobackup/pjjg18/neotrans/inst/matrices" \ +CAT_CSV="/nobackup/pjjg18/reach/mbank_catalogue.csv" \ +OUT_DIR="/nobackup/pjjg18/reach/ab" \ +TASK_ID="${SLURM_ARRAY_TASK_ID:-1}" \ +N_SEEDS=5 \ + Rscript /nobackup/pjjg18/reach/reach_escalation_ab.R diff --git a/dev/benchmarks/reach_escalation_analyze.R b/dev/benchmarks/reach_escalation_analyze.R new file mode 100644 index 000000000..ec02bcd18 --- /dev/null +++ b/dev/benchmarks/reach_escalation_analyze.R @@ -0,0 +1,137 @@ +#!/usr/bin/env Rscript +# Analyse the reach-escalation A/B. Usage: Rscript reach_escalation_analyze.R +# +# Per (matrix, seed) cell the TARGET is the union-best final score across arms (the +# established mbank convention -- no canonical best-known table exists for mbank). +# Per arm: +# reached = final_score <= target +# tt_hit = first elapsed_s whose best_score <= target (NA if never reached) +# rep2hit = first replicate whose best_score <= target (NA if never reached) +# Anytime win = smaller tt_hit at equal-or-better reach; regression = larger tt_hit / +# worse reach. rep2hit is reported ALONGSIDE tt_hit because a wall-only read +# misattributes "costlier per replicate" as "slower to the optimum" (the 2026-07-14 +# kick study: wall2hit -7..-12% but rep2hit == 1.000). +# +# PRE-REGISTERED RULE (fixed before results): +# SHIP : reach(deltas) >= reach(base) overall AND no tier reach-regression. +# NO-SHIP : reach(deltas) < reach(base) on the general pool. +# Wall cost is accepted by construction (the gate is opt-in; the user raised +# targetHits meaning "time is no issue") but is reported honestly. + +args <- commandArgs(trailingOnly = TRUE) +dir <- if (length(args) >= 1L) args[[1]] else "." +files <- list.files(dir, pattern = "^cell_.*\\.csv$", full.names = TRUE) +if (!length(files)) stop("no cell_*.csv in ", dir) +D <- do.call(rbind, lapply(files, read.csv, stringsAsFactors = FALSE)) +cat(sprintf("Loaded %d rows from %d cell files\n", nrow(D), length(files))) + +cells <- unique(D[, c("dataset", "nTip", "tier", "seed")]) +out <- list() +for (i in seq_len(nrow(cells))) { + cl <- cells[i, ] + sub <- D[D$dataset == cl$dataset & D$seed == cl$seed, ] + target <- min(sub$final_score, na.rm = TRUE) # union-best across arms + for (a in unique(sub$arm)) { + sa <- sub[sub$arm == a, ] + sa <- sa[order(sa$elapsed_s), ] + hit <- which(sa$best_score <= target + 1e-9) + out[[length(out) + 1L]] <- data.frame( + dataset = cl$dataset, nTip = cl$nTip, tier = cl$tier, seed = cl$seed, arm = a, + target = target, final = sa$final_score[1], + reached = as.integer(sa$final_score[1] <= target + 1e-9), + tt_hit = if (length(hit)) sa$elapsed_s[hit[1]] else NA_real_, + rep2hit = if (length(hit)) sa$replicate[hit[1]] else NA_integer_, + wall_total = sa$wall_total_s[1], reps = sa$reps_done[1], + cap_s = sa$cap_s[1], + # TRUNCATION FLAG. The deltas arm costs ~15x the candidates per replicate, so at a + # fixed wall it completes far fewer reps. If an arm stopped AT the cap it did not + # converge -- its "reach" is a budget artefact, not a property of the config. This + # is exactly the error that made 5432 arm B look replicate-capped when it was + # time-truncated. A reach comparison is only honest on non-truncated cells. + truncated = as.integer(!is.na(sa$wall_total_s[1]) && + sa$wall_total_s[1] >= 0.95 * sa$cap_s[1]), + stringsAsFactors = FALSE) + } +} +P <- do.call(rbind, out) +write.csv(P, file.path(dir, "reach_ab_per_arm.csv"), row.names = FALSE) + +fmt <- function(x) if (all(is.na(x))) "NA" else sprintf("%.3g", median(x, na.rm = TRUE)) +cat("\n=== REACH by arm (fraction of cells attaining the cell's union-best) ===\n") +for (a in sort(unique(P$arm))) + cat(sprintf(" %-6s reach = %.3f (%d/%d)\n", a, + mean(P$arm == a & P$reached == 1L) / mean(P$arm == a), + sum(P$arm == a & P$reached == 1L), sum(P$arm == a))) + +cat("\n=== REACH by tier (the tier-regression check) ===\n") +for (tr in c("small", "medium", "large", "xlarge")) { + s <- P[P$tier == tr, ] + if (!nrow(s)) next + cat(sprintf(" %-7s", tr)) + for (a in sort(unique(P$arm))) + cat(sprintf(" %s=%.3f (%d/%d)", a, mean(s$reached[s$arm == a]), + sum(s$reached[s$arm == a]), sum(s$arm == a))) + cat("\n") +} + +cat("\n=== MEDIAN tt_hit (wall s) / rep2hit (replicates) / total wall ===\n") +for (a in sort(unique(P$arm))) + cat(sprintf(" %-6s tt_hit=%-8s rep2hit=%-6s wall_total=%-8s reps=%s\n", a, + fmt(P$tt_hit[P$arm == a]), fmt(P$rep2hit[P$arm == a]), + fmt(P$wall_total[P$arm == a]), fmt(P$reps[P$arm == a]))) + +# Paired per-cell comparison (repo convention: counts + median ratios, not p-values). +b <- P[P$arm == "base", ]; d <- P[P$arm == "deltas", ] +k <- intersect(paste(b$dataset, b$seed), paste(d$dataset, d$seed)) +b <- b[match(k, paste(b$dataset, b$seed)), ]; d <- d[match(k, paste(d$dataset, d$seed)), ] +cat(sprintf("\n=== PAIRED (n = %d cells) ===\n", length(k))) +cat(sprintf(" final score : %d better, %d worse, %d tie (deltas vs base)\n", + sum(d$final < b$final), sum(d$final > b$final), sum(d$final == b$final))) +cat(sprintf(" reach : base %d, deltas %d\n", sum(b$reached), sum(d$reached))) +ok <- !is.na(b$tt_hit) & !is.na(d$tt_hit) +if (any(ok)) { + cat(sprintf(" tt_hit ratio (deltas/base), median = %.3f [%d better, %d worse]\n", + median(d$tt_hit[ok] / b$tt_hit[ok]), + sum(d$tt_hit[ok] < b$tt_hit[ok]), sum(d$tt_hit[ok] > b$tt_hit[ok]))) + okr <- ok & !is.na(b$rep2hit) & !is.na(d$rep2hit) & b$rep2hit > 0 + if (any(okr)) + cat(sprintf(" rep2hit ratio (deltas/base), median = %.3f <- if ~1.0, any wall gap is COST-PER-REP, not slower reach\n", + median(d$rep2hit[okr] / b$rep2hit[okr]))) +} +cat(sprintf(" wall_total ratio (deltas/base), median = %.3f\n", + median(d$wall_total / b$wall_total, na.rm = TRUE))) + +worse <- d$final > b$final +if (any(worse)) { + cat("\n!! cells where deltas found a WORSE tree (the NO-SHIP signal):\n") + print(data.frame(dataset = d$dataset[worse], tier = d$tier[worse], seed = d$seed[worse], + base = b$final[worse], deltas = d$final[worse]), row.names = FALSE) +} else cat("\nNo cell where deltas found a worse tree.\n") + +cat(sprintf("\n=== TRUNCATION (stopped at the wall cap => did NOT converge) ===\n")) +cat(sprintf(" base %d/%d cells truncated\n deltas %d/%d cells truncated\n", + sum(b$truncated), nrow(b), sum(d$truncated), nrow(d))) +clean <- b$truncated == 0L & d$truncated == 0L +cat(sprintf(" cells where NEITHER arm truncated (the honest reach comparison): %d/%d\n", + sum(clean), length(clean))) +if (sum(d$truncated) > sum(b$truncated)) + cat(" NOTE: deltas truncated more often than base -- on those cells a reach gap is a\n", + " BUDGET artefact (deltas cost ~15x candidates/rep), not a config failure.\n") + +# The verdict is computed on NON-TRUNCATED cells only: a truncated arm never converged, +# so scoring its reach would repeat the arm-B error of reading a budget cut as a result. +reachB <- mean(b$reached[clean]); reachD <- mean(d$reached[clean]) +tierBad <- character(0) +for (tr in unique(P$tier)) { + ib <- clean & b$tier == tr; id <- clean & d$tier == tr + if (sum(ib) && mean(d$reached[id]) < mean(b$reached[ib])) tierBad <- c(tierBad, tr) +} +cat(sprintf("\n=== PRE-REGISTERED VERDICT (non-truncated cells, n = %d) ===\n", sum(clean))) +cat(sprintf(" reach base=%.3f deltas=%.3f; tier regressions: %s\n --> %s\n", + reachB, reachD, if (length(tierBad)) paste(tierBad, collapse = ",") else "none", + if (sum(clean) < 0.5 * length(clean)) + "INCONCLUSIVE -- too few non-truncated cells; re-run with larger caps" + else if (reachD >= reachB && !length(tierBad)) "SHIP v1" + else "NO-SHIP (restrict or drop)")) +cat("\n(All-cells reach, for reference only -- confounded by truncation: ") +cat(sprintf("base=%.3f deltas=%.3f)\n", mean(b$reached), mean(d$reached))) From 994ba410e468a1ed621d3ab59b0f629669289f0a Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 31 Jul 2026 08:23:41 +0100 Subject: [PATCH 03/16] fix(benchmarks): read A/B arm names from the data, not hard-coded The paired table selected arms by the literal names "base"/"deltas", so a variant harness with a differently-named arm produced an EMPTY comparison rather than an error -- a silent wrong answer. Names are now taken from the data, with "base" preferred as the baseline when present. Co-Authored-By: Claude Opus 4.8 --- dev/benchmarks/reach_escalation_analyze.R | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/dev/benchmarks/reach_escalation_analyze.R b/dev/benchmarks/reach_escalation_analyze.R index ec02bcd18..715d78520 100644 --- a/dev/benchmarks/reach_escalation_analyze.R +++ b/dev/benchmarks/reach_escalation_analyze.R @@ -81,13 +81,22 @@ for (a in sort(unique(P$arm))) fmt(P$wall_total[P$arm == a]), fmt(P$reps[P$arm == a]))) # Paired per-cell comparison (repo convention: counts + median ratios, not p-values). -b <- P[P$arm == "base", ]; d <- P[P$arm == "deltas", ] +# Arm names are READ FROM THE DATA rather than hard-coded, so a variant harness (e.g. a +# `deltas6` arm confirming the shipped six-lever form) analyses without an edit -- with +# fixed names the paired table silently comes back empty instead of erroring. +armNames <- sort(unique(P$arm)) +baseArm <- if ("base" %in% armNames) "base" else armNames[[1]] +testArm <- setdiff(armNames, baseArm)[[1]] +b <- P[P$arm == baseArm, ]; d <- P[P$arm == testArm, ] k <- intersect(paste(b$dataset, b$seed), paste(d$dataset, d$seed)) b <- b[match(k, paste(b$dataset, b$seed)), ]; d <- d[match(k, paste(d$dataset, d$seed)), ] -cat(sprintf("\n=== PAIRED (n = %d cells) ===\n", length(k))) -cat(sprintf(" final score : %d better, %d worse, %d tie (deltas vs base)\n", - sum(d$final < b$final), sum(d$final > b$final), sum(d$final == b$final))) -cat(sprintf(" reach : base %d, deltas %d\n", sum(b$reached), sum(d$reached))) +cat(sprintf("\n=== PAIRED (n = %d cells; baseline '%s' vs test '%s') ===\n", + length(k), baseArm, testArm)) +cat(sprintf(" final score : %d better, %d worse, %d tie (%s vs %s)\n", + sum(d$final < b$final), sum(d$final > b$final), sum(d$final == b$final), + testArm, baseArm)) +cat(sprintf(" reach : %s %d, %s %d\n", baseArm, sum(b$reached), + testArm, sum(d$reached))) ok <- !is.na(b$tt_hit) & !is.na(d$tt_hit) if (any(ok)) { cat(sprintf(" tt_hit ratio (deltas/base), median = %.3f [%d better, %d worse]\n", From e3f13c27a154830f5a997ab09eb8918d5e97e850 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 31 Jul 2026 09:27:11 +0100 Subject: [PATCH 04/16] bench(reach): confirm shipped 6-lever form; fix deadline detector Confirmation A/B of the form that actually ships (six levers, gated), job 18127149: 11 large/xlarge MBANK_FIXED_SAMPLE matrices x 5 seeds. Paired final score 8 better / 1 worse / 46 tie, no tier regression -- SHIP confirmed. Record what the numbers do NOT say: - The entire effect is 2 of 11 matrices (project4284 5/5; project2771 3/1/1). The "xlarge 20/20 vs 15/20" line is one matrix, as in the ship gate. Reach against a union-best-across-arms target is self-referential, so it inflates the loss count; paired counts are the statistic. - project4284 won having completed ZERO replicates -- one deep replicate, not more restarts. That is the opposite mechanism to the hard-tail study's restart-volume conclusion; both recorded, neither generalised. - project2771 is this battery's high-variance matrix (it also produced the ship gate's single loss); single cells there are noise. Fix a real bug in the committed analyzer: the truncation flag tested `wall >= 0.95 * cap_s`, but the engine stops at `maxSeconds * (1 - enumTimeFraction)` = 0.90 * cap_s, so it reported 1 of 110 deadline-bound ab6 cells when the true count was 59 -- the very trap documented three sections above it in the FINDINGS. The harness now records enum_time_fraction and the analyzer derives the deadline from it. Deadline-bound is not spoiled: when both arms stop at one deadline the cell is a valid equal-wall comparison (the cheaper-per-rep arm gets 2.24x the replicates and still has to win). Only asymmetry invalidates a cell. Discarding every deadline-bound cell, as the old code did, would have dropped 40 of 55 informative cells and printed INCONCLUSIVE over a clean result. The analyzer now classifies three ways and prints a per-matrix win/loss/tie table so a one-matrix effect cannot read as a tier property. Also: split FINDINGS "Not measured" into measured-on-score vs unit-tested-only (trigger threshold, strategy scoping, userSet skip, escalatedPool filter -- all inert in both benchmark runs, which passed levers as dots on a gate-free engine), and note that poolReseed (v2) cannot help on the regime carrying v1's biggest win, since it needs completed replicates to form a pool. Co-Authored-By: Claude Opus 4.8 --- dev/benchmarks/reach_escalation_FINDINGS.md | 93 ++++++++++++++- dev/benchmarks/reach_escalation_ab.R | 11 +- dev/benchmarks/reach_escalation_ab6.R | 107 ++++++++++++++++++ dev/benchmarks/reach_escalation_analyze.R | 118 +++++++++++++++----- dev/benchmarks/reach_recover4284.R | 83 ++++++++++++++ 5 files changed, 376 insertions(+), 36 deletions(-) create mode 100644 dev/benchmarks/reach_escalation_ab6.R create mode 100644 dev/benchmarks/reach_recover4284.R diff --git a/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md index ddbffb5b0..64fc3e28f 100644 --- a/dev/benchmarks/reach_escalation_FINDINGS.md +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -71,6 +71,48 @@ The one loss: `project2771` (large, seed 5821), base 911 → deltas 912. Agains regresses; wall cost is accepted by construction, since the escalation only fires when the user has asked for it. → **SHIP.** +## Confirmation A/B of the shipped form (job 18127149, 2026-07-31) + +The gap left by the ship gate — the shipped configuration is *six* levers, and was inferred +from a seven-lever run — is now closed on score. `reach_ab6.R`: the same design, restricted +to the 11 `large`+`xlarge` members of `MBANK_FIXED_SAMPLE` (on small/medium the gate is inert, +so those cells carry no information), 5 seeds, 55 cells, both arms at the same raised +`targetHits`, levers as dots on a gate-free engine. + +**Verdict: SHIP — confirmed.** Paired final score **8 better, 1 worse, 46 tie**; no tier +regression. Budget regime clean: 40 cells both-arms-at-deadline (equal-wall), 15 cells both +converged, **0 asymmetric**. + +**The whole effect is 2 of the 11 matrices.** Read this before quoting any tier number: + +| matrix | tips | win | loss | tie | +|--------|------|-----|------|-----| +| project4284 | 4062 | **5** | 0 | 0 | +| project2771 | 94 | 3 | 1 | 1 | +| the other nine (63–173 t) | | 0 | 0 | 45 | + +So the analyzer's `xlarge = 20/20 vs 15/20` is *one matrix*, exactly as in the ship gate. +Do not restate reach fractions as evidence: the union-best target is self-referential, so +base "misses" on 8 cells only because deltas6 got there. **The paired counts are the +statistic.** `project2771` is the *high-variance matrix of this battery* — it produced both +this run's single loss (911→912) and the ship gate's single loss, while deltas6 reached 911 on +4 of 5 seeds; treat any single 2771 cell as noise rather than signal. + +**project4284 won with ZERO completed replicates.** base completed 1–2 replicates in ~1330 s; +deltas6 completed **0** and still returned a tree 2–9 steps better on every seed. At 4062 +tips the escalation's benefit is *spending the whole budget deepening a single replicate*, not +more restarts. Note this sits crosswise to the hard-tail study below, which concluded reach +is a **restart-volume** phenomenon — different regimes (482 vs 4062 tips), not a +contradiction, but the two mechanisms are opposite and neither generalises to the other. + +**Cost.** In this run `wall_total` ratio is 1.000 — a tautology of the shared deadline, not a +measurement; the honest figure is that at equal wall the escalation completes **2.24× fewer +replicates**. On the 30 tied cells where both arms were deadline-bound (7 matrices, 86–173 +tips) that cost nothing in final score. The 16 tied cells that converged early (project4286 +at 5–9 s of a 720 s cap, project4359 stopping on `targetHits` at 28 replicates) carry **no** +cost information and must not be counted as evidence of cost-neutrality. The ×3.56 wall +figure belongs to the ship-gate run, which was not deadline-bound. + ## Why it is gated, and gated to `thorough`/`large` Below ~200 tips the levers buy no reach and cost ~3.5× the wall, so they must not be @@ -111,6 +153,24 @@ recipe improves the odds, it does not guarantee the floor.** enumTimeFraction)`, default 10% reserve. Identical elapsed times across seeds is the time-truncation signature; arm B and arm D were truncated exactly there while appearing replicate-capped. For pure reach runs set `enumTimeFraction = 0`. +- **…and writing that lesson down did not stop me re-committing it.** `reach_escalation_ + analyze.R` flagged truncation at `wall >= 0.95 × cap_s`, but cells stop at `0.90 × cap_s`, + so it reported **1 of 110** deadline-bound ab6 cells when the true count was 59 — three + sections below the paragraph above. Both the harness (records `enum_time_fraction`) and the + analyzer (derives the deadline from it, defaulting to 0.1) are fixed. A threshold that + encodes a magic number the engine owns will rot: derive it from recorded data. +- **Deadline-bound is not the same as spoiled.** When *both* arms stop at the same deadline + the cell is a valid equal-wall comparison — the stronger test, since the cheaper-per-replicate + arm gets ~2.2× more replicates and still has to win. What invalidates a cell is *asymmetry* + (one arm converged, one cut off). The analyzer now classifies cells three ways rather than + discarding every deadline-bound cell, which would have thrown away 40 of 55 informative + cells and printed "INCONCLUSIVE" over a clean result. +- **A tier win is not a tier property until you count matrices.** Twice now an `xlarge` + reach jump has been one matrix repeated across five seeds. The analyzer prints the + per-matrix win/loss/tie table and the "n of N matrices changed" line for this reason. +- **Union-best-across-arms targets make reach self-referential.** If one arm alone attains a + score the other misses *by construction*, so reach fractions inflate the loss count. + Report paired win/loss/tie; use reach only as a secondary description. - **Whole-suite test runs need `test_local()`/`devtools::test()`**, not `test_dir()` + `library()`: several test files call internals unqualified and error otherwise. - **A top-level `skip_on_cran()` makes a file report "0 pass 0 fail"** — that is *not* a @@ -121,9 +181,30 @@ recipe improves the odds, it does not guarantee the floor.** ## Not measured -- The shipped configuration was never A/B'd *as such*: the run above tested seven levers - ungated; six ship, gated to `thorough`/`large`. The inference is a tier decomposition of - that run (the xlarge win is under `large` and preserved; the excluded tiers measured 0 - better / 0 worse, so the gate should remove only cost), not a direct measurement. -- Equal weights only. Use under implied weights and profile parsimony is a deliberate but - unmeasured extrapolation; ratchet depth there is governed separately. +Measured as of job 18127149: the **six levers' effect on score**, under `auto`→`thorough`/ +`large`, equal weights, 61–4062 tips. + +Still **unit-tested only** — no benchmark evidence, because that run passed the levers as dots +on a gate-free engine and therefore never executed the gate: + +- the trigger threshold (`targetHits / defaultHits >= 2`); +- the `thorough`/`large` scoping; +- the `userSet` skip (caller-supplied fields must survive); +- the `escalatedPool` return filter. `escalatedPool` was FALSE throughout both benchmark + runs, so the filter that stops `poolSuboptimal = 3` leaking suboptimal trees into a + `collapse = FALSE` result has never run outside `test-reach-escalation.R`. + +Also not measured: + +- **Equal weights only.** Implied weights and profile parsimony are a deliberate but + unmeasured extrapolation; ratchet depth there is governed separately by `.IwRatchetDepth`. +- **No tree was retained for the project4284 result** the confirmation rests on — the harness + recorded `attr(res, "score")` only, so ~353 could not be independently re-scored at the + time of writing. A degenerate partial tree would score *worse*, so the win is very + unlikely to be an artefact, but "unlikely" is what the record says. Recovery run + (`recover4284.R`, job 18128376) re-runs the identical config and re-scores by label with + `TreeLength`; see the result note below. +- **`poolReseed` (v2) cannot help where v1 helped most.** It reseeds *replicates* from a + pool needing `size >= 2`; project4284's winning arm completed zero replicates. Any v2 + validation must therefore use matrices on which enough replicates complete for a pool to + form — 4062 tips is the wrong test bed for it. diff --git a/dev/benchmarks/reach_escalation_ab.R b/dev/benchmarks/reach_escalation_ab.R index a0268c90a..dcfb5d95a 100644 --- a/dev/benchmarks/reach_escalation_ab.R +++ b/dev/benchmarks/reach_escalation_ab.R @@ -137,12 +137,20 @@ make_tracer <- function(t0) { list(cb = cb, env = env) } +# The engine reserves a fraction of maxSeconds for MPT enumeration and stops the main search +# at maxSeconds * (1 - ENUM_TIME_FRACTION). Passed EXPLICITLY (at the engine default, so +# behaviour is unchanged) and recorded per row, so the analyzer can compute the real deadline +# instead of guessing at the nominal cap -- guessing is what let 59/110 deadline-bound ab6 +# cells be scored as "converged". +ENUM_TIME_FRACTION <- 0.1 + run_arm <- function(pd, nTip, arm, seed, maxrep, cap_s) { set.seed(seed) t0 <- proc.time()["elapsed"] tr <- make_tracer(t0) args <- list(pd, strategy = "auto", maxReplicates = maxrep, maxSeconds = cap_s, + enumTimeFraction = ENUM_TIME_FRACTION, targetHits = reach_threshold(nTip), # SAME in both arms nThreads = 1L, verbosity = 0L, progressCallback = tr$cb) @@ -213,7 +221,8 @@ for (arm in ARMS) { event = "improve", replicate = tr$replicate, elapsed_s = tr$elapsed_s, engine_elapsed = tr$engine_elapsed, best_score = tr$best_score, final_score = r$final_score, reps_done = r$reps, wall_total_s = r$wall, - candidates = r$cand, cap_s = cap_s, stringsAsFactors = FALSE) + candidates = r$cand, cap_s = cap_s, + enum_time_fraction = ENUM_TIME_FRACTION, stringsAsFactors = FALSE) } D <- do.call(rbind, all_rows) diff --git a/dev/benchmarks/reach_escalation_ab6.R b/dev/benchmarks/reach_escalation_ab6.R new file mode 100644 index 000000000..2be5472bc --- /dev/null +++ b/dev/benchmarks/reach_escalation_ab6.R @@ -0,0 +1,107 @@ +#!/usr/bin/env Rscript +# Confirmation A/B of the SHIPPED escalation form (lead 5 / "belt-and-braces"). +# +# The ship-gate A/B measured SEVEN levers, UNGATED, via strategy="auto" over all 25 +# MBANK_FIXED_SAMPLE matrices. What actually ships is SIX levers (ratchetCycles dropped to +# .IwRatchetDepth) GATED to strategy thorough/large. That configuration was never run as +# such; the ship argument is a tier decomposition of the earlier run. This closes that gap. +# +# Scope: only the matrices where auto resolves to thorough (61-120 tips) or large (>=121), +# i.e. the 11 large+xlarge members of MBANK_FIXED_SAMPLE -- on small/medium the gate now +# does nothing, so those cells would be identical in both arms and carry no information. +# 11 matrices x 5 seeds = 55 cells. Everything else matches the original harness: both arms +# at the SAME raised targetHits (the trigger threshold), differing ONLY in the six levers; +# gate-free engine with the levers as dots, so the result does not depend on the gate +# implementation; union-best target; training split asserted. +suppressMessages({ + ts_lib <- Sys.getenv("TS_LIB", "") + if (nzchar(ts_lib)) { + library(TreeSearch, lib.loc = normalizePath(ts_lib, winslash = "/", mustWork = TRUE)) + } else library(TreeSearch) + library(TreeTools) +}) +neo_dir <- Sys.getenv("NEOTRANS_DIR"); cat_csv <- Sys.getenv("CAT_CSV") +out_dir <- Sys.getenv("OUT_DIR", "."); dir.create(out_dir, showWarnings = FALSE, recursive = TRUE) + +# large (61-120) + xlarge (>=121) members of MBANK_FIXED_SAMPLE +KEYS <- c("project4286", "project4359", "project4397", "project2084_(1)", + "project2771", "project2184", "project3938", + "syab07201", "project4133", "project804", "project4284") + +# The SHIPPED six (ratchetCycles deliberately absent -- owned by .IwRatchetDepth). +DELTAS6 <- list(ratchetPerturbMaxMoves = 0L, driftCycles = 25L, + postRatchetSectorial = TRUE, stallEscalateFactor = 1.5, + intraFuse = TRUE, poolSuboptimal = 3) +reach_threshold <- function(nTip) 2L * max(10L, as.integer(nTip / 5)) + +catalogue <- read.csv(cat_csv, stringsAsFactors = FALSE) +rownames(catalogue) <- catalogue$key +to_fitch <- function(pd) { m <- PhyDatToMatrix(pd, ambigNA = FALSE); m[m == "-"] <- "?"; MatrixToPhyDat(m) } +load_matrix <- function(key) { + row <- catalogue[key, ] + if (!identical(row$split, "training")) + stop(sprintf("key %s is split='%s' -- validation is SEQUESTERED", key, row$split)) + to_fitch(suppressWarnings(TreeTools::ReadAsPhyDat(file.path(neo_dir, row$filename)))) +} +`%||%` <- function(a, b) if (is.null(a)) b else a +make_tracer <- function(t0) { + env <- new.env(parent = emptyenv()); env$prev <- Inf; env$rows <- list() + cb <- function(info) { + if (!identical(info$phase, "replicate")) return(invisible()) + bs <- info$best_score + if (is.null(bs) || length(bs) != 1L || !is.finite(bs)) return(invisible()) + if (bs < env$prev - 1e-9) { env$prev <- bs + env$rows[[length(env$rows) + 1L]] <- data.frame( + replicate = as.integer(info$replicate %||% NA_integer_), + elapsed_s = as.double(proc.time()["elapsed"] - t0), + best_score = as.double(bs), stringsAsFactors = FALSE) } + invisible() + } + list(cb = cb, env = env) +} +run_arm <- function(pd, nTip, arm, seed, maxrep, cap_s) { + set.seed(seed); t0 <- proc.time()["elapsed"]; tr <- make_tracer(t0) + args <- list(pd, strategy = "auto", maxReplicates = maxrep, maxSeconds = cap_s, + targetHits = reach_threshold(nTip), nThreads = 1L, verbosity = 0L, + progressCallback = tr$cb) + if (identical(arm, "deltas6")) args <- c(args, DELTAS6) + res <- suppressWarnings(do.call(MaximizeParsimony, args)) + wall <- as.double(proc.time()["elapsed"] - t0) + trace <- if (length(tr$env$rows)) do.call(rbind, tr$env$rows) else + data.frame(replicate = NA_integer_, elapsed_s = NA_real_, + best_score = as.double(attr(res, "score"))) + list(wall = wall, final_score = as.double(attr(res, "score")), + reps = attr(res, "replicates") %||% NA_integer_, + cand = attr(res, "candidates_evaluated") %||% NA_real_, + trace = trace, n_events = length(tr$env$rows)) +} +N_SEEDS <- as.integer(Sys.getenv("N_SEEDS", "5")); BASE_SEED <- 7731L +maxrep <- as.integer(Sys.getenv("TS_MAXREP", "300")) +manifest <- expand.grid(key = KEYS, seed_idx = seq_len(N_SEEDS), stringsAsFactors = FALSE) +manifest <- manifest[order(manifest$key, manifest$seed_idx), ]; rownames(manifest) <- NULL +tid <- as.integer(Sys.getenv("TASK_ID", Sys.getenv("SLURM_ARRAY_TASK_ID", "1"))) +if (tid < 1 || tid > nrow(manifest)) stop(sprintf("TASK_ID %d out of 1..%d", tid, nrow(manifest))) +key <- manifest$key[tid]; seed <- BASE_SEED + manifest$seed_idx[tid] - 1L +row <- catalogue[key, ]; nTip <- as.integer(row$ntax) +tier <- if (nTip >= 121L) "xlarge" else "large" +cap_s <- if (identical(tier, "xlarge")) 1440 else 720 +cat(sprintf("=== ab6 task %d/%d: %s (%dt, %s) seed=%d targetHits=%d cap=%gs ===\n", + tid, nrow(manifest), key, nTip, tier, seed, reach_threshold(nTip), cap_s)) +pd <- load_matrix(key); stopifnot(length(pd) == nTip) +rows <- list() +for (arm in c("base", "deltas6")) { + r <- run_arm(pd, nTip, arm, seed, maxrep, cap_s) + cat(sprintf(" %-8s final=%.0f reps=%s wall=%.1fs events=%d\n", + arm, r$final_score, r$reps, r$wall, r$n_events)) + if (r$n_events == 0L) cat(" WARN: 0 improvement events -- trace may be broken\n") + tr <- r$trace + rows[[length(rows) + 1L]] <- data.frame( + dataset = key, nTip = nTip, tier = tier, seed = seed, arm = arm, + target_hits = reach_threshold(nTip), event = "improve", replicate = tr$replicate, + elapsed_s = tr$elapsed_s, best_score = tr$best_score, final_score = r$final_score, + reps_done = r$reps, wall_total_s = r$wall, candidates = r$cand, cap_s = cap_s, + stringsAsFactors = FALSE) +} +D <- do.call(rbind, rows) +of <- file.path(out_dir, sprintf("cell_%03d_%s_s%d.csv", tid, gsub("[^A-Za-z0-9]", "", key), seed)) +write.csv(D, of, row.names = FALSE); cat(sprintf("Wrote %s (%d rows)\n", of, nrow(D))) diff --git a/dev/benchmarks/reach_escalation_analyze.R b/dev/benchmarks/reach_escalation_analyze.R index 715d78520..fba5b6aff 100644 --- a/dev/benchmarks/reach_escalation_analyze.R +++ b/dev/benchmarks/reach_escalation_analyze.R @@ -25,6 +25,18 @@ if (!length(files)) stop("no cell_*.csv in ", dir) D <- do.call(rbind, lapply(files, read.csv, stringsAsFactors = FALSE)) cat(sprintf("Loaded %d rows from %d cell files\n", nrow(D), length(files))) +# The engine's real stopping deadline, not the nominal cap (see the at_deadline comment). +ENUM_TIME_FRACTION_DEFAULT <- 0.1 +deadline_of <- function(sa) { + etf <- if ("enum_time_fraction" %in% names(sa) && !is.na(sa$enum_time_fraction[1])) { + sa$enum_time_fraction[1] + } else ENUM_TIME_FRACTION_DEFAULT + sa$cap_s[1] * (1 - etf) +} +if (!"enum_time_fraction" %in% names(D)) + cat(sprintf("NOTE: no enum_time_fraction column; assuming the %.2f default for deadlines\n", + ENUM_TIME_FRACTION_DEFAULT)) + cells <- unique(D[, c("dataset", "nTip", "tier", "seed")]) out <- list() for (i in seq_len(nrow(cells))) { @@ -43,13 +55,21 @@ for (i in seq_len(nrow(cells))) { rep2hit = if (length(hit)) sa$replicate[hit[1]] else NA_integer_, wall_total = sa$wall_total_s[1], reps = sa$reps_done[1], cap_s = sa$cap_s[1], - # TRUNCATION FLAG. The deltas arm costs ~15x the candidates per replicate, so at a - # fixed wall it completes far fewer reps. If an arm stopped AT the cap it did not - # converge -- its "reach" is a budget artefact, not a property of the config. This - # is exactly the error that made 5432 arm B look replicate-capped when it was - # time-truncated. A reach comparison is only honest on non-truncated cells. - truncated = as.integer(!is.na(sa$wall_total_s[1]) && - sa$wall_total_s[1] >= 0.95 * sa$cap_s[1]), + # DEADLINE FLAG. The deltas arm costs ~2.3x the wall per replicate, so at a fixed + # budget it completes far fewer reps. An arm that stopped at the budget did not + # converge -- reading its "reach" as a property of the config is exactly the error + # that made 5432 arm B look replicate-capped when it was time-truncated. + # + # The budget is NOT maxSeconds. The engine stops the main search at + # main_deadline = maxSeconds * (1 - enumTimeFraction) [src/ts_driven.cpp] + # with enumTimeFraction defaulting to 0.1, so a deadline-bound cell lands at ~0.90 * + # cap_s and NOT at cap_s. This flag originally tested `>= 0.95 * cap_s` and therefore + # scored 59 of 110 genuinely deadline-bound ab6 cells as "converged" -- it missed the + # very trap documented in reach_escalation_FINDINGS.md. Read enumTimeFraction from the + # data when the harness records it, else assume the 0.1 default. + deadline_s = deadline_of(sa), + at_deadline = as.integer(!is.na(sa$wall_total_s[1]) && + sa$wall_total_s[1] >= 0.98 * deadline_of(sa)), stringsAsFactors = FALSE) } } @@ -117,30 +137,70 @@ if (any(worse)) { base = b$final[worse], deltas = d$final[worse]), row.names = FALSE) } else cat("\nNo cell where deltas found a worse tree.\n") -cat(sprintf("\n=== TRUNCATION (stopped at the wall cap => did NOT converge) ===\n")) -cat(sprintf(" base %d/%d cells truncated\n deltas %d/%d cells truncated\n", - sum(b$truncated), nrow(b), sum(d$truncated), nrow(d))) -clean <- b$truncated == 0L & d$truncated == 0L -cat(sprintf(" cells where NEITHER arm truncated (the honest reach comparison): %d/%d\n", - sum(clean), length(clean))) -if (sum(d$truncated) > sum(b$truncated)) - cat(" NOTE: deltas truncated more often than base -- on those cells a reach gap is a\n", - " BUDGET artefact (deltas cost ~15x candidates/rep), not a config failure.\n") +# BUDGET REGIME. Being deadline-bound is not automatically a spoiled cell: when BOTH arms +# stop at the same deadline the cell is a valid EQUAL-WALL comparison, which is the stronger +# test (the cheaper-per-rep arm gets more replicates and still has to win). What invalidates +# a cell is ASYMMETRY -- one arm converged and the other was cut off -- because then the +# score gap may be purely budget. So classify cells three ways instead of dropping them. +cat(sprintf("\n=== BUDGET REGIME (deadline = cap_s x (1 - enumTimeFraction)) ===\n")) +cat(sprintf(" base %d/%d cells stopped at the deadline\n", sum(b$at_deadline), nrow(b))) +cat(sprintf(" %-6s %d/%d cells stopped at the deadline\n", testArm, + sum(d$at_deadline), nrow(d))) +bothDL <- b$at_deadline == 1L & d$at_deadline == 1L +neither <- b$at_deadline == 0L & d$at_deadline == 0L +asym <- !bothDL & !neither +cat(sprintf(" both at deadline (valid, EQUAL-WALL) : %d\n", sum(bothDL))) +cat(sprintf(" neither (both converged, valid) : %d\n", sum(neither))) +cat(sprintf(" exactly one (ASYMMETRIC, suspect) : %d\n", sum(asym))) +if (any(asym)) + print(data.frame(dataset = d$dataset[asym], seed = d$seed[asym], + baseAtDL = b$at_deadline[asym], testAtDL = d$at_deadline[asym], + base = b$final[asym], test = d$final[asym]), row.names = FALSE) +if (any(bothDL)) { + rr <- b$reps[bothDL] / d$reps[bothDL] + rr <- rr[is.finite(rr)] + if (length(rr)) + cat(sprintf(" work per replicate on equal-wall cells: %s does %.2fx the replicates\n", + baseArm, median(rr))) +} -# The verdict is computed on NON-TRUNCATED cells only: a truncated arm never converged, -# so scoring its reach would repeat the arm-B error of reading a budget cut as a result. -reachB <- mean(b$reached[clean]); reachD <- mean(d$reached[clean]) +# The verdict is computed on the VALID cells (both-at-deadline plus both-converged) and is +# driven by PAIRED SCORE COUNTS, not by reach. Reach here is measured against the union-best +# across arms, which is self-referential -- if one arm alone attains a score the other +# "misses" by construction -- so a reach fraction restates the paired counts with the losses +# inflated. Both are printed; the counts are the statistic. +valid <- bothDL | neither +nBetter <- sum(d$final[valid] < b$final[valid]) +nWorse <- sum(d$final[valid] > b$final[valid]) +reachB <- mean(b$reached[valid]); reachD <- mean(d$reached[valid]) tierBad <- character(0) for (tr in unique(P$tier)) { - ib <- clean & b$tier == tr; id <- clean & d$tier == tr - if (sum(ib) && mean(d$reached[id]) < mean(b$reached[ib])) tierBad <- c(tierBad, tr) + iv <- valid & b$tier == tr + if (sum(iv) && sum(d$final[iv] > b$final[iv]) > sum(d$final[iv] < b$final[iv])) + tierBad <- c(tierBad, tr) } -cat(sprintf("\n=== PRE-REGISTERED VERDICT (non-truncated cells, n = %d) ===\n", sum(clean))) -cat(sprintf(" reach base=%.3f deltas=%.3f; tier regressions: %s\n --> %s\n", - reachB, reachD, if (length(tierBad)) paste(tierBad, collapse = ",") else "none", - if (sum(clean) < 0.5 * length(clean)) - "INCONCLUSIVE -- too few non-truncated cells; re-run with larger caps" - else if (reachD >= reachB && !length(tierBad)) "SHIP v1" +cat(sprintf("\n=== PRE-REGISTERED VERDICT (valid cells, n = %d of %d) ===\n", + sum(valid), length(valid))) +cat(sprintf(" paired score: %d better, %d worse; tier regressions: %s\n", + nBetter, nWorse, if (length(tierBad)) paste(tierBad, collapse = ",") else "none")) +cat(sprintf(" reach (union-best, self-referential): base=%.3f %s=%.3f\n", + reachB, testArm, reachD)) +cat(sprintf(" --> %s\n", + if (!sum(valid)) "INCONCLUSIVE -- no valid cells" + else if (nBetter >= nWorse && !length(tierBad)) "SHIP" else "NO-SHIP (restrict or drop)")) -cat("\n(All-cells reach, for reference only -- confounded by truncation: ") -cat(sprintf("base=%.3f deltas=%.3f)\n", mean(b$reached), mean(d$reached))) + +# How concentrated is the effect? A tier-level win can be a single matrix repeated across +# seeds -- that happened here (project4284) and was briefly written up as a tier property. +chg <- valid & d$final != b$final +if (any(chg)) { + cat("\n=== WHERE THE EFFECT LIVES (per matrix; a 1-matrix effect is NOT a tier property) ===\n") + for (ds in unique(P$dataset[P$dataset %in% b$dataset[chg]])) { + i <- valid & b$dataset == ds + cat(sprintf(" %-18s %5dt win %d loss %d tie %d\n", ds, b$nTip[i][1], + sum(d$final[i] < b$final[i]), sum(d$final[i] > b$final[i]), + sum(d$final[i] == b$final[i]))) + } + cat(sprintf(" matrices with any change: %d of %d in the battery\n", + length(unique(b$dataset[chg])), length(unique(b$dataset)))) +} diff --git a/dev/benchmarks/reach_recover4284.R b/dev/benchmarks/reach_recover4284.R new file mode 100644 index 000000000..367498d83 --- /dev/null +++ b/dev/benchmarks/reach_recover4284.R @@ -0,0 +1,83 @@ +#!/usr/bin/env Rscript +# Recover the TREE for the project4284 result that the ab6 confirmation A/B rests on. +# +# WHY: reach_ab6.R recorded only attr(res, "score"). project4284 (4062 tips) is the matrix +# carrying the entire ab6 win (5/5 seeds, 353-357 vs base 359-364), and its deltas6 arm +# completed ZERO replicates -- it returned a tree found mid-replicate-0 when the deadline +# fired. A score attribute with no tree cannot be re-scored, and the standing rule here is +# to persist every new best tree the same turn. Every project4284 best_score on record in +# dev/benchmarks/ is 1040-1411 (all 30-120s budget-starved), so ~353 is the best value seen +# for this matrix anywhere -- exactly the thing that must not exist as a bare number. +# +# Config is IDENTICAL to reach_ab6.R's deltas6 arm (same seeds, same cap, same levers, same +# preprocessing) so the scores should reproduce; reproducing them also re-validates +# determinism. The only change is that the tree is written out and independently re-scored +# BY LABEL (never by raw edge+tip_data index -- RenumberTips permutes). +suppressMessages({ + ts_lib <- Sys.getenv("TS_LIB", "") + if (nzchar(ts_lib)) { + library(TreeSearch, lib.loc = normalizePath(ts_lib, winslash = "/", mustWork = TRUE)) + } else library(TreeSearch) + library(TreeTools) +}) +neo_dir <- Sys.getenv("NEOTRANS_DIR"); cat_csv <- Sys.getenv("CAT_CSV") +out_dir <- Sys.getenv("OUT_DIR", "."); dir.create(out_dir, showWarnings = FALSE, recursive = TRUE) + +KEY <- "project4284" +`%||%` <- function(a, b) if (is.null(a)) b else a +DELTAS6 <- list(ratchetPerturbMaxMoves = 0L, driftCycles = 25L, + postRatchetSectorial = TRUE, stallEscalateFactor = 1.5, + intraFuse = TRUE, poolSuboptimal = 3) + +catalogue <- read.csv(cat_csv, stringsAsFactors = FALSE) +rownames(catalogue) <- catalogue$key +row <- catalogue[KEY, ] +if (!identical(row$split, "training")) + stop(sprintf("key %s is split='%s' -- validation is SEQUESTERED", KEY, row$split)) +to_fitch <- function(pd) { + m <- PhyDatToMatrix(pd, ambigNA = FALSE); m[m == "-"] <- "?"; MatrixToPhyDat(m) +} +pd <- to_fitch(suppressWarnings(TreeTools::ReadAsPhyDat(file.path(neo_dir, row$filename)))) +nTip <- as.integer(row$ntax); stopifnot(length(pd) == nTip) + +tid <- as.integer(Sys.getenv("TASK_ID", Sys.getenv("SLURM_ARRAY_TASK_ID", "1"))) +seed <- 7731L + tid - 1L # ab6 used BASE_SEED 7731 + seed_idx - 1 +cap_s <- 1440 # ab6 xlarge cap, unchanged +targetHits <- 2L * max(10L, as.integer(nTip / 5)) +cat(sprintf("=== recover4284 task %d: seed=%d targetHits=%d cap=%gs (%d tips) ===\n", + tid, seed, targetHits, cap_s, nTip)) + +set.seed(seed); t0 <- proc.time()["elapsed"] +res <- suppressWarnings(do.call(MaximizeParsimony, c( + list(pd, strategy = "auto", maxReplicates = 300L, maxSeconds = cap_s, + targetHits = targetHits, nThreads = 1L, verbosity = 0L), DELTAS6))) +wall <- as.double(proc.time()["elapsed"] - t0) +reported <- as.double(attr(res, "score")) + +# INDEPENDENT re-score, by label, of every returned tree. TreeLength() re-derives the score +# from the tree + dataset, so agreement with attr(res,"score") is a real check that the +# returned object is a valid tree scoring what the harness claimed -- the point of the run. +trees <- if (inherits(res, "phylo")) structure(list(res), class = "multiPhylo") else res +lengths_ <- vapply(trees, function(tr) as.double(TreeLength(tr, pd, concavity = Inf)), + double(1)) +cat(sprintf(" reported=%.0f nTrees=%d re-scored: min=%.0f max=%.0f wall=%.1fs\n", + reported, length(trees), min(lengths_), max(lengths_), wall)) +agree <- isTRUE(all.equal(min(lengths_), reported)) +cat(sprintf(" RE-SCORE AGREES WITH REPORTED SCORE: %s\n", agree)) +if (!agree) + cat(" !! MISMATCH -- the reported score is NOT reproduced by TreeLength on the tree.\n") +binaryOK <- vapply(trees, function(tr) length(tr$edge[, 1]) == 2L * nTip - 3L, logical(1)) +cat(sprintf(" fully resolved (unrooted binary): %d/%d\n", sum(binaryOK), length(trees))) + +best <- trees[lengths_ <= min(lengths_) + 1e-9] +tf <- file.path(out_dir, sprintf("project4284_deltas6_s%d_score%.0f.tre", seed, min(lengths_))) +ape::write.tree(best, file = tf) +cat(sprintf(" WROTE %d tree(s) -> %s\n", length(best), tf)) +write.csv(data.frame(dataset = KEY, nTip = nTip, seed = seed, arm = "deltas6", + reported_score = reported, rescored_min = min(lengths_), + rescore_agrees = agree, n_trees = length(trees), + n_best = length(best), all_binary = all(binaryOK), + reps = attr(res, "replicates") %||% NA_integer_, + wall_s = wall, cap_s = cap_s, tree_file = basename(tf), + stringsAsFactors = FALSE), + file.path(out_dir, sprintf("recover_%02d_s%d.csv", tid, seed)), row.names = FALSE) From d269c6a6870cc9ad0b20403351865ac6d3c1dc2f Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 31 Jul 2026 09:30:28 +0100 Subject: [PATCH 05/16] docs(reach): correct two ship-gate claims the fixed detector falsified Re-ran the corrected analyzer over the original 125-cell ship-gate data. Two statements in the FINDINGS were wrong and are now fixed: - "No cell in either arm hit its wall cap, so nothing is a budget artefact" -- false, and a direct artefact of the 0.95*cap_s bug. Truth: 36 base / 40 deltas cells stopped at the deadline, and FOUR cells are asymmetric (project2184, deltas at the deadline while base converged). All four tied at 563, so no score comparison rests on an asymmetric budget -- but the clean-sweep claim does not stand and should not be repeated. - The ship-gate effect is also only 2 of 25 matrices: project4284 5/0/0 and project2771 1 win/1 loss/3 tie. So the `large` row's "1 better, 1 worse" is BOTH project2771 -- the matrix that is high-variance in both runs. That makes project4284 the only matrix with a clean win in either A/B, which is a narrower claim than the tier table alone suggests. The x3.56 wall figure is confirmed as belonging to this run (median 3.559). Co-Authored-By: Claude Opus 4.8 --- dev/benchmarks/reach_escalation_FINDINGS.md | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md index 64fc3e28f..a39070e75 100644 --- a/dev/benchmarks/reach_escalation_FINDINGS.md +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -42,7 +42,12 @@ self-referential: if one arm alone attains a score, the other "misses" by constr the reach fractions restate the paired counts rather than measuring absolute optimality. The paired win counts below are the honest statistic. -**Result.** No cell in either arm hit its wall cap, so nothing is a budget artefact. +**Result.** Budget regime (re-derived 2026-07-31 with the *fixed* deadline detector — the +original write-up said "no cell in either arm hit its wall cap", which was **wrong**, an +artefact of the `0.95 × cap_s` bug): 36 base and 40 deltas cells stopped at the deadline, 85 +cells converged in both arms, and **4 cells are asymmetric** (project2184 seeds 5821/5823/ +5824/5825 — deltas at the deadline, base not). All four tied at 563, so no score comparison +here rests on an asymmetric budget, but the clean-sweep claim does not stand. Strict paired final-score wins (deltas vs base): @@ -58,6 +63,12 @@ tips, which improved on every one of its five seeds (by 1–9 steps). The other matrices (125, 131, 173 tips) all tied. So the demonstrated benefit is *datasets far too large to converge within an ordinary budget*, **not** a property of "over 120 tips". +Per-matrix concentration for this run (added 2026-07-31): only **2 of the 25 matrices changed +at all** — `project4284` 5 win / 0 loss / 0 tie, and `project2771` **1 win / 1 loss / 3 tie**. +So the `large` row's "1 better, 1 worse" is *both* project2771, i.e. the one matrix that is +demonstrably high-variance in both runs. **project4284 is the only matrix with a clean win in +either A/B.** + Cost, median over cells: total wall ×3.56, time-to-best ×2.49, but **replicates-to-best ×1.00**. The wall gap is entirely work per replicate (~15× candidates evaluated against a `sprint` baseline; ~1.5× against `thorough`), not slower convergence. Reporting wall alone From 896f1dd64fa0fb14d219de532c53de51a53175d9 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 31 Jul 2026 09:32:19 +0100 Subject: [PATCH 06/16] docs(reach): project2771 is noise, not a second winning matrix Pooled across both A/Bs project2771 is 4 win / 2 loss / 4 tie over 10 seeds -- a coin flip -- and it supplied the single loss in each run. The ab6 table listed it beside project4284 in a way that reads as two contributing matrices, which is inconsistent with the correction made to the ship-gate section in d269c6a6. State the surviving claim precisely instead: project4284 is 10 win / 0 loss over 10 seeds across both runs, and no matrix regresses net. That is the whole of the positive evidence. It is enough to ship, and it is one matrix -- the largest in the battery. Also note in reach_escalation_ab.R's header that enumTimeFraction was added to its argument list post-hoc (same value as the engine default, so behaviour is unchanged) and that the 125-cell CSVs therefore predate the enum_time_fraction column. The analyzer's "assuming the 0.10 default" note on that data is correct behaviour, not a defect for a later session to "fix". Co-Authored-By: Claude Opus 4.8 --- dev/benchmarks/reach_escalation_FINDINGS.md | 13 ++++++++++--- dev/benchmarks/reach_escalation_ab.R | 7 +++++++ 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md index a39070e75..78d1bee56 100644 --- a/dev/benchmarks/reach_escalation_FINDINGS.md +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -105,9 +105,16 @@ converged, **0 asymmetric**. So the analyzer's `xlarge = 20/20 vs 15/20` is *one matrix*, exactly as in the ship gate. Do not restate reach fractions as evidence: the union-best target is self-referential, so base "misses" on 8 cells only because deltas6 got there. **The paired counts are the -statistic.** `project2771` is the *high-variance matrix of this battery* — it produced both -this run's single loss (911→912) and the ship gate's single loss, while deltas6 reached 911 on -4 of 5 seeds; treat any single 2771 cell as noise rather than signal. +statistic.** + +**`project2771` is noise, not a second winning matrix.** Pooled across both A/Bs it is +**4 win / 2 loss / 4 tie** over 10 seeds — a coin flip, and it supplied the single loss in each +run. Read the two-row table above as *one* winning matrix plus one high-variance matrix, and +treat any single 2771 cell as noise. + +**The claim that survives both runs:** `project4284` is **10 win / 0 loss over 10 seeds**, and +**no matrix regresses net** in either A/B. That is the whole of the positive evidence, and it +is enough — but it is one matrix, and it is the largest one in the battery. **project4284 won with ZERO completed replicates.** base completed 1–2 replicates in ~1330 s; deltas6 completed **0** and still returned a tree 2–9 steps better on every seed. At 4062 diff --git a/dev/benchmarks/reach_escalation_ab.R b/dev/benchmarks/reach_escalation_ab.R index dcfb5d95a..4074ccd5f 100644 --- a/dev/benchmarks/reach_escalation_ab.R +++ b/dev/benchmarks/reach_escalation_ab.R @@ -1,6 +1,13 @@ #!/usr/bin/env Rscript # v1 REACH-ESCALATION GATE -- general-pool anytime A/B (the SHIP/NO-SHIP gate). # +# POST-HOC EDIT, 2026-07-31: `enumTimeFraction = ENUM_TIME_FRACTION` was ADDED to the argument +# list and `enum_time_fraction` to the output columns, so the analyzer can compute the engine's +# real deadline instead of guessing. The value equals the engine default, so search behaviour +# is unchanged -- but this file is no longer textually identical to the script that produced +# the 125-cell result, and those CSVs therefore lack the column (the analyzer says so, and +# falls back to the 0.1 default; that note is correct, not a defect to "fix"). +# # WHAT IS BEING DECIDED: `MaximizeParsimony()` gains a gate that, when the user raises # `targetHits` to >= 2 * max(10, nTip/5), applies 7 deeper per-replicate perturbation # deltas (ratchetCycles 40, kick 0/auto-deep, driftCycles 25, postRatchetSectorial, From 58ae9f72ea21d7f1d6ca5ad957b643e6f344d1ba Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 31 Jul 2026 09:47:57 +0100 Subject: [PATCH 07/16] bench(reach): fix tree recovery (collapse=FALSE), measure poolSuboptimal leak Job 18128376 died after 22 minutes of search with "`tree` must be binary": MaximizeParsimony defaults to collapse = TRUE, which contracts zero-length branches, and TreeLength refuses a non-binary tree. Neither is at fault -- the check was mis-specified. Fixes, and one lesson worth more than the fix: - Ask for collapse = FALSE, so returned trees are binary and re-scorable. - WRITE THE TREES BEFORE VERIFYING THEM. v1 lost a whole search because it died in the verification step before persisting anything. The tree is the deliverable; the check is not. - Score each tree in a tryCatch and report per-tree, so one bad tree cannot take down the run. The corrected form also buys a measurement for free. A collapse = FALSE return is the pool verbatim, and under a gate-free harness escalatedPool is FALSE, so the filter that strips poolSuboptimal trees is inert -- comparing max against min of the re-scored set therefore tests directly whether poolSuboptimal = 3 leaks suboptimal trees to the caller. That was listed as unit-tested-only in "Not measured". Re-run is job 18128526. While checking that: confirmed by reading R/MaximizeParsimony.R that the collapse = TRUE branch already restricts to scores == best_score before collapsing, so the default path never leaked and the v1 return-path fix sits on the only branch that needed it. Recorded, since "did the fix cover both branches?" is the obvious question to ask of it. Co-Authored-By: Claude Opus 4.8 --- dev/benchmarks/reach_escalation_FINDINGS.md | 20 ++++- dev/benchmarks/reach_recover4284.R | 98 +++++++++++++++------ 2 files changed, 89 insertions(+), 29 deletions(-) diff --git a/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md index 78d1bee56..d48951fab 100644 --- a/dev/benchmarks/reach_escalation_FINDINGS.md +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -189,6 +189,16 @@ recipe improves the odds, it does not guarantee the floor.** - **Union-best-across-arms targets make reach self-referential.** If one arm alone attains a score the other misses *by construction*, so reach fractions inflate the loss count. Report paired win/loss/tie; use reach only as a secondary description. +- **You cannot `TreeLength()` a tree that `MaximizeParsimony()` returned by default.** + `collapse = TRUE` (the default) contracts zero-length branches, and `TreeLength` errors with + "`tree` must be binary". To re-score a returned tree, ask for `collapse = FALSE`. This + killed the first tree-recovery run (job 18128376) after 22 minutes of search, in the + *verification* step — so **persist the deliverable before verifying it**: write the trees + out first, then check them, or a failed check throws away the search too. +- **A `collapse = FALSE` return is the pool verbatim.** That makes it the natural place to + test whether `poolSuboptimal` leaks suboptimal trees to the caller: re-score every returned + tree and compare max against min. Under a gate-free harness `escalatedPool` is FALSE, so + the filter is inert and the leak (if real) is visible. - **Whole-suite test runs need `test_local()`/`devtools::test()`**, not `test_dir()` + `library()`: several test files call internals unqualified and error otherwise. - **A top-level `skip_on_cran()` makes a file report "0 pass 0 fail"** — that is *not* a @@ -220,8 +230,14 @@ Also not measured: recorded `attr(res, "score")` only, so ~353 could not be independently re-scored at the time of writing. A degenerate partial tree would score *worse*, so the win is very unlikely to be an artefact, but "unlikely" is what the record says. Recovery run - (`recover4284.R`, job 18128376) re-runs the identical config and re-scores by label with - `TreeLength`; see the result note below. + `reach_recover4284.R` re-runs the identical config, writes the trees and re-scores by label + with `TreeLength` (job 18128376 died in the check — see the `collapse` trap above; job + **18128526** is the corrected `collapse = FALSE` run). +- **The `escalatedPool` filter is only needed on the `collapse = FALSE` path** — verified by + reading, not benchmark: the `collapse = TRUE` branch already restricts to + `scores == best_score` before collapsing (`R/MaximizeParsimony.R`), which the `collapse` + roxygen also documents. So the default path never leaked; the fix sits on the one branch + that did. - **`poolReseed` (v2) cannot help where v1 helped most.** It reseeds *replicates* from a pool needing `size >= 2`; project4284's winning arm completed zero replicates. Any v2 validation must therefore use matrices on which enough replicates complete for a pool to diff --git a/dev/benchmarks/reach_recover4284.R b/dev/benchmarks/reach_recover4284.R index 367498d83..1fd0321fd 100644 --- a/dev/benchmarks/reach_recover4284.R +++ b/dev/benchmarks/reach_recover4284.R @@ -1,5 +1,6 @@ #!/usr/bin/env Rscript -# Recover the TREE for the project4284 result that the ab6 confirmation A/B rests on. +# Recover the TREE for the project4284 result that the ab6 confirmation A/B rests on, +# and measure the poolSuboptimal leak while we are here. # # WHY: reach_ab6.R recorded only attr(res, "score"). project4284 (4062 tips) is the matrix # carrying the entire ab6 win (5/5 seeds, 353-357 vs base 359-364), and its deltas6 arm @@ -9,10 +10,20 @@ # dev/benchmarks/ is 1040-1411 (all 30-120s budget-starved), so ~353 is the best value seen # for this matrix anywhere -- exactly the thing that must not exist as a bare number. # -# Config is IDENTICAL to reach_ab6.R's deltas6 arm (same seeds, same cap, same levers, same -# preprocessing) so the scores should reproduce; reproducing them also re-validates -# determinism. The only change is that the tree is written out and independently re-scored -# BY LABEL (never by raw edge+tip_data index -- RenumberTips permutes). +# v2 of this script. v1 (job 18128376) died in TreeLength with "`tree` must be binary": +# MaximizeParsimony defaults to collapse = TRUE, which contracts zero-length branches, and +# TreeLength refuses a non-binary tree. Not a bug in either -- a mis-specified check. +# +# So run with collapse = FALSE, which is both re-scorable AND a free measurement. With the +# levers passed as DOTS (gate-free engine) escalatedPool is FALSE, so the return-path filter +# that normally strips poolSuboptimal trees does NOT fire -- meaning the returned set should +# contain suboptimal pool trees if poolSuboptimal = 3 really leaks. That leak is currently +# listed as unit-tested-only in reach_escalation_FINDINGS.md "Not measured"; comparing +# max(rescored) against min(rescored) measures it directly, on real data. +# +# Config is otherwise IDENTICAL to reach_ab6.R's deltas6 arm (same seeds, cap, levers, +# preprocessing), so scores should reproduce; reproducing them re-validates determinism. +# collapse is post-processing on the returned pool and does not alter the search. suppressMessages({ ts_lib <- Sys.getenv("TS_LIB", "") if (nzchar(ts_lib)) { @@ -44,40 +55,73 @@ tid <- as.integer(Sys.getenv("TASK_ID", Sys.getenv("SLURM_ARRAY_TASK_ID", "1"))) seed <- 7731L + tid - 1L # ab6 used BASE_SEED 7731 + seed_idx - 1 cap_s <- 1440 # ab6 xlarge cap, unchanged targetHits <- 2L * max(10L, as.integer(nTip / 5)) -cat(sprintf("=== recover4284 task %d: seed=%d targetHits=%d cap=%gs (%d tips) ===\n", +cat(sprintf("=== recover4284 v2 task %d: seed=%d targetHits=%d cap=%gs (%d tips) ===\n", tid, seed, targetHits, cap_s, nTip)) set.seed(seed); t0 <- proc.time()["elapsed"] res <- suppressWarnings(do.call(MaximizeParsimony, c( list(pd, strategy = "auto", maxReplicates = 300L, maxSeconds = cap_s, - targetHits = targetHits, nThreads = 1L, verbosity = 0L), DELTAS6))) + targetHits = targetHits, nThreads = 1L, verbosity = 0L, + collapse = FALSE), # <- binary trees, and exposes the pool verbatim + DELTAS6))) wall <- as.double(proc.time()["elapsed"] - t0) reported <- as.double(attr(res, "score")) - -# INDEPENDENT re-score, by label, of every returned tree. TreeLength() re-derives the score -# from the tree + dataset, so agreement with attr(res,"score") is a real check that the -# returned object is a valid tree scoring what the harness claimed -- the point of the run. trees <- if (inherits(res, "phylo")) structure(list(res), class = "multiPhylo") else res -lengths_ <- vapply(trees, function(tr) as.double(TreeLength(tr, pd, concavity = Inf)), - double(1)) -cat(sprintf(" reported=%.0f nTrees=%d re-scored: min=%.0f max=%.0f wall=%.1fs\n", - reported, length(trees), min(lengths_), max(lengths_), wall)) -agree <- isTRUE(all.equal(min(lengths_), reported)) +cat(sprintf(" reported=%.0f nTrees=%d reps=%s wall=%.1fs\n", reported, length(trees), + attr(res, "replicates") %||% "?", wall)) + +# WRITE THE TREES FIRST. v1 of this script lost a 22-minute search because it died in the +# verification step before persisting anything; the tree is the deliverable, the check is not. +allf <- file.path(out_dir, sprintf("project4284_deltas6_s%d_ALL.tre", seed)) +ape::write.tree(trees, file = allf) +cat(sprintf(" wrote all %d returned tree(s) -> %s\n", length(trees), basename(allf))) + +# Now verify: resolution first (TreeLength demands binary), then re-score BY LABEL. +nEdge <- vapply(trees, function(tr) nrow(tr$edge), integer(1)) +binaryOK <- nEdge == 2L * nTip - 3L +cat(sprintf(" fully resolved (unrooted binary, %d edges): %d/%d edge counts seen: %s\n", + 2L * nTip - 3L, sum(binaryOK), length(trees), + paste(sort(unique(nEdge)), collapse = ","))) +lengths_ <- rep(NA_real_, length(trees)) +for (i in seq_along(trees)) { + lengths_[i] <- tryCatch(as.double(TreeLength(trees[[i]], pd, concavity = Inf)), + error = function(e) { cat(sprintf(" TreeLength failed on tree %d: %s\n", + i, conditionMessage(e))); NA_real_ }) +} +ok <- !is.na(lengths_) +agree <- any(ok) && isTRUE(all.equal(min(lengths_[ok]), reported)) +cat(sprintf(" re-scored %d/%d: min=%s max=%s\n", sum(ok), length(trees), + if (any(ok)) sprintf("%.0f", min(lengths_[ok])) else "NA", + if (any(ok)) sprintf("%.0f", max(lengths_[ok])) else "NA")) cat(sprintf(" RE-SCORE AGREES WITH REPORTED SCORE: %s\n", agree)) -if (!agree) - cat(" !! MISMATCH -- the reported score is NOT reproduced by TreeLength on the tree.\n") -binaryOK <- vapply(trees, function(tr) length(tr$edge[, 1]) == 2L * nTip - 3L, logical(1)) -cat(sprintf(" fully resolved (unrooted binary): %d/%d\n", sum(binaryOK), length(trees))) +if (any(ok) && !agree) + cat(" !! MISMATCH -- reported score NOT reproduced by TreeLength on any returned tree.\n") -best <- trees[lengths_ <= min(lengths_) + 1e-9] -tf <- file.path(out_dir, sprintf("project4284_deltas6_s%d_score%.0f.tre", seed, min(lengths_))) -ape::write.tree(best, file = tf) -cat(sprintf(" WROTE %d tree(s) -> %s\n", length(best), tf)) +# The poolSuboptimal-leak measurement: any returned tree scoring worse than the best is a +# suboptimal pool tree that reached the caller, which is what the escalatedPool filter exists +# to prevent (it is inert here by design -- levers as dots, so escalatedPool == FALSE). +nAbove <- if (any(ok)) sum(lengths_[ok] > min(lengths_[ok]) + 1e-9) else NA_integer_ +cat(sprintf(" SUBOPTIMAL trees in the returned set: %s of %d -> leak %s\n", + nAbove, sum(ok), + if (is.na(nAbove)) "UNKNOWN" else if (nAbove > 0) "REAL (filter is load-bearing)" + else "not observed on this cell")) + +if (any(ok)) { + best <- trees[ok][lengths_[ok] <= min(lengths_[ok]) + 1e-9] + bf <- file.path(out_dir, sprintf("project4284_deltas6_s%d_BEST_score%.0f.tre", + seed, min(lengths_[ok]))) + ape::write.tree(best, file = bf) + cat(sprintf(" wrote %d best tree(s) -> %s\n", length(best), basename(bf))) +} write.csv(data.frame(dataset = KEY, nTip = nTip, seed = seed, arm = "deltas6", - reported_score = reported, rescored_min = min(lengths_), + reported_score = reported, + rescored_min = if (any(ok)) min(lengths_[ok]) else NA_real_, + rescored_max = if (any(ok)) max(lengths_[ok]) else NA_real_, rescore_agrees = agree, n_trees = length(trees), - n_best = length(best), all_binary = all(binaryOK), + n_rescored = sum(ok), n_suboptimal = nAbove, + n_binary = sum(binaryOK), all_binary = all(binaryOK), reps = attr(res, "replicates") %||% NA_integer_, - wall_s = wall, cap_s = cap_s, tree_file = basename(tf), + wall_s = wall, cap_s = cap_s, collapse_arg = FALSE, stringsAsFactors = FALSE), file.path(out_dir, sprintf("recover_%02d_s%d.csv", tid, seed)), row.names = FALSE) +cat(" CSV written.\n") From 13c6728cf252d3db865c0fcec0a7e65cd8002e25 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 31 Jul 2026 10:17:35 +0100 Subject: [PATCH 08/16] bench(reach): recover project4284 trees; margin is not a stable effect size Job 18128526 re-ran the identical deltas6 config on all 5 seeds with collapse = FALSE, wrote the trees, and re-scored every returned tree by label. CONFIRMED: on all five seeds the re-scored minimum exactly equals the reported attr(res,"score"), over 100 returned trees per seed (500 trees, all scored, all fully resolved). The ab6 headline is not an artefact of a degenerate partial tree. Trees are now in Hamilton floors/. Best tree HELD = 354 (seed 7732); best score ever observed = 353 (ab6 seed 7731), whose tree is lost. QUALIFIED: the per-seed scores did NOT reproduce -- 353/356/357 came back as 355/360/363, mean 355.0 -> 357.6. Not a determinism bug: these cells complete zero replicates and stop on the wall clock, so node speed moves the answer. Elsewhere in the battery 46 cells tied exactly, so the instability is confined to the matrix that never converges. Consequences, stated in the FINDINGS: - The 10 win / 0 loss count STANDS. Both A/Bs ran base and deltas6 sequentially in one R session on one node, so the pairing controls node speed; the re-run had no base arm and cannot add or remove wins. - But the 2-9 step margin must NOT be quoted as an effect size. Cross-run spread is up to 6 steps, comparable to the effect, and the re-run's seed 7734 would have lost to four of five original base cells. The paired within-task design is load-bearing. - The poolSuboptimal leak was not observed, and this was the wrong cell to test it on: all 100 trees shared one score, but with zero completed replicates the retention path plausibly never engaged. Still unmeasured, not disproved. Also fix the resolution check: trees come back ROOTED, so fully resolved is 2n-2 edges, not the unrooted 2n-3. The old test reported "0/100 binary" for trees that TreeLength -- which refuses non-binary input -- had just scored without complaint. Observed 8122 = 2*4062-2. Co-Authored-By: Claude Opus 4.8 --- dev/benchmarks/reach_escalation_FINDINGS.md | 57 ++++++++++++++++++--- dev/benchmarks/reach_recover4284.R | 11 ++-- 2 files changed, 57 insertions(+), 11 deletions(-) diff --git a/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md index d48951fab..02107d1db 100644 --- a/dev/benchmarks/reach_escalation_FINDINGS.md +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -114,7 +114,49 @@ treat any single 2771 cell as noise. **The claim that survives both runs:** `project4284` is **10 win / 0 loss over 10 seeds**, and **no matrix regresses net** in either A/B. That is the whole of the positive evidence, and it -is enough — but it is one matrix, and it is the largest one in the battery. +is enough — but it is one matrix, and it is the largest one in the battery. Read the +reproducibility note below before quoting the *size* of that win. + +### project4284 tree recovery, and what it says about the margin (job 18128526) + +Re-ran the identical deltas6 config on all 5 seeds with `collapse = FALSE`, wrote the trees and +re-scored every returned tree by label with `TreeLength`. + +**Good news — the result is real.** On all five seeds the re-scored minimum **exactly equals** +the reported `attr(res, "score")`, across 100 returned trees per seed (500 trees, all scored, +all fully resolved: 8122 edges = 2n−2, the *rooted* binary count). So the ab6 headline is not +an artefact of a degenerate partial tree, and the trees now exist in Hamilton `floors/`. +Best tree **held** = **354** (seed 7732). Best score ever *observed* = 353 (ab6 seed 7731) — +**that tree is lost**, so 354 is the best recoverable tree for this matrix. + +**Bad news — the per-seed scores did not reproduce.** + +| seed | 7731 | 7732 | 7733 | 7734 | 7735 | mean | +|------|------|------|------|------|------|------| +| ab6 deltas6 | 353 | 354 | 356 | 357 | 355 | 355.0 | +| re-run deltas6 | 355 | 354 | **360** | **363** | 356 | 357.6 | +| ab6 base | 361 | 360 | 361 | 359 | 364 | 361.0 | + +Only seed 7732 matched. **This is not a determinism bug**: these cells complete *zero* +replicates and stop on the wall clock, so the search halts wherever the clock happened to be +and node speed and load move the answer. (Contrast the 46 exact ties elsewhere in the battery +— on matrices that actually converge, the score is stable. The instability is confined to the +matrix that never converges.) + +What this does and does not change: + +- **The 10/0 count stands.** Both A/Bs ran base and deltas6 sequentially *in one R session on + one node*, so the pairing controls node speed. The re-run had no base arm and so cannot add + or remove wins. +- **Do not quote the 2–9 step margin as an effect size.** The cross-run spread on this matrix + is up to 6 steps — comparable to the effect — and the re-run's seed 7734 (363) would have + *lost* to four of the five original base cells. The paired within-task design is + load-bearing: the effect is only cleanly visible because of it. +- **The `poolSuboptimal` leak was not observed here, and this was the wrong cell to test it + on.** All 100 returned trees shared one score (the pool sits at `poolMaxSize` = 100), but + with zero completed replicates the suboptimal-retention path plausibly never engaged — the + same structural reason `poolReseed` cannot work on this matrix. Treat the leak as still + unmeasured, not as disproved. **project4284 won with ZERO completed replicates.** base completed 1–2 replicates in ~1330 s; deltas6 completed **0** and still returned a tree 2–9 steps better on every seed. At 4062 @@ -226,13 +268,12 @@ Also not measured: - **Equal weights only.** Implied weights and profile parsimony are a deliberate but unmeasured extrapolation; ratchet depth there is governed separately by `.IwRatchetDepth`. -- **No tree was retained for the project4284 result** the confirmation rests on — the harness - recorded `attr(res, "score")` only, so ~353 could not be independently re-scored at the - time of writing. A degenerate partial tree would score *worse*, so the win is very - unlikely to be an artefact, but "unlikely" is what the record says. Recovery run - `reach_recover4284.R` re-runs the identical config, writes the trees and re-scores by label - with `TreeLength` (job 18128376 died in the check — see the `collapse` trap above; job - **18128526** is the corrected `collapse = FALSE` run). +- ~~No tree was retained for the project4284 result.~~ **RESOLVED** by job 18128526 — see the + recovery section above. Trees are in Hamilton `floors/`; the reported score is exactly + reproduced by `TreeLength` on the returned tree. What the recovery *added* to the + "not measured" list is that **the per-seed score on this matrix is not reproducible across + nodes** (deadline-truncated, zero completed replicates), so the margin is not a stable + effect size. - **The `escalatedPool` filter is only needed on the `collapse = FALSE` path** — verified by reading, not benchmark: the `collapse = TRUE` branch already restricts to `scores == best_score` before collapsing (`R/MaximizeParsimony.R`), which the `collapse` diff --git a/dev/benchmarks/reach_recover4284.R b/dev/benchmarks/reach_recover4284.R index 1fd0321fd..37cf5848b 100644 --- a/dev/benchmarks/reach_recover4284.R +++ b/dev/benchmarks/reach_recover4284.R @@ -77,10 +77,15 @@ ape::write.tree(trees, file = allf) cat(sprintf(" wrote all %d returned tree(s) -> %s\n", length(trees), basename(allf))) # Now verify: resolution first (TreeLength demands binary), then re-score BY LABEL. +# MaximizeParsimony returns trees ROOTED (the engine roots at tip 0), so fully resolved means +# 2n-2 edges, NOT the unrooted 2n-3. Testing 2n-3 reported "0/100 binary" on trees that +# TreeLength -- which itself refuses non-binary input -- had just scored without complaint; +# the observed count was 8122 = 2*4062-2. Accept either, and print the counts so a future +# mismatch is diagnosable instead of mysterious. nEdge <- vapply(trees, function(tr) nrow(tr$edge), integer(1)) -binaryOK <- nEdge == 2L * nTip - 3L -cat(sprintf(" fully resolved (unrooted binary, %d edges): %d/%d edge counts seen: %s\n", - 2L * nTip - 3L, sum(binaryOK), length(trees), +binaryOK <- nEdge %in% c(2L * nTip - 2L, 2L * nTip - 3L) +cat(sprintf(" fully resolved (%d rooted / %d unrooted edges): %d/%d edge counts seen: %s\n", + 2L * nTip - 2L, 2L * nTip - 3L, sum(binaryOK), length(trees), paste(sort(unique(nEdge)), collapse = ","))) lengths_ <- rep(NA_real_, length(trees)) for (i in seq_along(trees)) { From 3913e0fd32678aaef01dcbb2aa0bfca7a514a121 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Fri, 31 Jul 2026 20:42:29 +0100 Subject: [PATCH 09/16] docs(reach): upstream `effort` breaks the gate's premise; correct the 4284 record Two findings from a verification sweep, both checked against origin/cpp-search and the archived March library rather than inferred. 1. SHIPPING BLOCKER. 419168d4 removed `strategy` in favour of `effort`, a relative offset on an internal ladder. Traced on the tip: .AutoRung(nTip>=120, nChar>=100) = 4; effort=+1 -> rung 5; hitMultiplier = 2^(5-4) = 2; line ~1124 multiplies targetHits by it ONLY when !userSetHits. So on any >=120-tip dataset, effort=+1 makes targetHits/defaultHits exactly 2.0, and .reachEscalationMinRatio = 2 tests >= -- the bundle fires on a mechanical ladder artefact, which is the precise opposite of the "user's own signal" premise it was justified by. Rung 5 has already doubled maxReplicates, so layering a ~3.5x-per-replicate bundle makes one notch of effort cost ~7x the work against docs promising about twice -- and the bundle's only positive evidence is one 4062-tip zero-replicate matrix while the 125/131/173-tip cells all tied. Recommends gating on `userSetHits`, which the tip already computes and whose doctrine its own roxygen states. 2. Correct the project4284 record. 354 is confirmed by five independent scorers and is the best RECOVERABLE tree; 353 has no artefact on disk and must not be quoted as attained. Two further caveats now recorded: - Gap treatment must be normalised across harnesses: "-"->"?" (Fitch) vs "-" as a sixth level (BGS) is worth +61 on the SAME tree (354 vs 415), over 7 of 27 characters. - The t252 CSVs are not a comparable baseline, and NOT because of budget as I previously wrote. Bare AdditionTree, same data and seed, no search: 1590 under the archived March library vs 409 under the current one -- worse and 5x faster, the signature of the since-fixed union-of-finals insertion-cost bug. The current engine returns the same score at maxSeconds 1/5/25/45 and the entire 30s->1440s span is x1.068, so time was never the driver. Co-Authored-By: Claude Opus 4.8 --- dev/benchmarks/reach_escalation_FINDINGS.md | 60 ++++++++++++++++++++- 1 file changed, 58 insertions(+), 2 deletions(-) diff --git a/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md index 02107d1db..cfc672e19 100644 --- a/dev/benchmarks/reach_escalation_FINDINGS.md +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -126,8 +126,25 @@ re-scored every returned tree by label with `TreeLength`. the reported `attr(res, "score")`, across 100 returned trees per seed (500 trees, all scored, all fully resolved: 8122 edges = 2n−2, the *rooted* binary count). So the ab6 headline is not an artefact of a degenerate partial tree, and the trees now exist in Hamilton `floors/`. -Best tree **held** = **354** (seed 7732). Best score ever *observed* = 353 (ab6 seed 7731) — -**that tree is lost**, so 354 is the best recoverable tree for this matrix. +Best tree **held** = **354** (seed 7732), independently confirmed by *five* scorers (current +`TreeLength`, the archived March `TreeLength`, phangorn Fitch, phangorn Sankoff, and a +hand-rolled Fitch). Best score ever *observed* = 353 (ab6 seed 7731) — **that tree is lost**, so +354 is the best recoverable tree and **353 should not be quoted as an attained score**: no +artefact on disk realises it. + +**Normalise the gap treatment before any cross-harness comparison.** These harnesses map +`"-" → "?"` (plain Fitch) via `to_fitch()`; the older `t252` sweeps keep `"-"` as a sixth level +(BGS three-pass). On the *same tree* that is worth **+61** — 415 under BGS versus 354 under +Fitch — localised to 7 of the 27 characters. Two numbers for this matrix are therefore not +comparable unless the gap handling matches. + +**Do not compare against the `t252` CSVs at all.** Their project4284 values (1040–1411) came +from an engine whose bare Wagner addition was several-fold worse: same data and seed, no search, +`AdditionTree` scores **1590** under the archived March library +(`/nobackup/pjjg18/TreeSearch/lib-t252`) against **409** under the current one — worse *and* +five times faster, the signature of the since-fixed union-of-finals insertion-cost bug. Budget +is not the explanation: the current engine returns the same score at `maxSeconds` 1, 5, 25 and +45, and the whole 30 s → 1440 s span is only ×1.068. **Bad news — the per-seed scores did not reproduce.** @@ -173,6 +190,45 @@ at 5–9 s of a 720 s cap, project4359 stopping on `targetHits` at 28 replicates cost information and must not be counted as evidence of cost-neutrality. The ×3.56 wall figure belongs to the ship-gate run, which was not deadline-bound. +## 🚨 The upstream API change breaks the gate's premise (2026-07-31, `419168d4`) + +`origin/cpp-search` **removed `strategy`** and replaced it with `effort = 0L`, a *relative* +offset on an internal ladder (`.effortLadder = sprint, default, thorough, large`, then rungs +that double the replicate cap). This is not a rename — it invalidates the gate's trigger. + +Traced on the tip: + + .AutoRung(nTip >= 120, nChar >= 100) -> 4 # `large` + effort = +1 -> rung 5 + spec$hitMultiplier at rung 5 = 2^(5-4) = 2 + line ~1124: if (!userSetHits && hitMultiplier > 1) + targetHits <- targetHits * hitMultiplier + +So on **any dataset of ≥120 tips, `effort = +1` makes `targetHits / defaultHits` exactly 2.0**, +and `.reachEscalationMinRatio = 2` tests `>=` — **the bundle fires**. Note the multiplication is +applied only when `!userSetHits`, i.e. precisely when the user did *not* raise the hit target. +The gate's entire justification is "the user's own signal that this dataset needs more"; under +the new API it would fire on a **mechanical ladder artefact** instead. + +The cost of getting this wrong is concrete. Rung 5 has already doubled `maxReplicates` +(500 → 1000); layering a bundle measured at ~3.5× wall per replicate on top makes one notch of +`effort` cost roughly **7× the work**, against documentation promising about twice. And the +bundle's only positive evidence is a single 4062-tip matrix that completes zero replicates — +the 125-, 131- and 173-tip cells all tied. **At 120–200 tips this is pure cost.** + +Options, with the recommendation: + +- **(a) leave it** — rejected: pure cost across the tier where it would newly fire. +- **(b) raise `.reachEscalationMinRatio` above 2** — only postpones the misfire; `effort = +2` + gives a multiplier of 4 and trips any threshold below that. +- **(c) gate on `userSetHits` (recommended)** — the tip already computes exactly that flag, and + its own roxygen states the doctrine ("a raised `targetHits` is taken as the user's own + signal"). This restores the premise precisely and is immune to future ladder changes. + +Open question for the maintainer: at high `effort` (say ≥ +3) the user arguably *does* want the +bundle. If so, the trigger becomes `userSetHits && ratio >= 2` **or** an explicit `effort` +threshold — a deliberate second door rather than an accident of the multiplier. + ## Why it is gated, and gated to `thorough`/`large` Below ~200 tips the levers buy no reach and cost ~3.5× the wall, so they must not be From d5039ea0d7c6f504c25ff392b3347c68dc5a4cbf Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Mon, 3 Aug 2026 15:55:18 +0100 Subject: [PATCH 10/16] docs(reach): drop the "1-9 steps" magnitude claim from NEWS Unsupported. Re-running the identical config on the 4062-tip matrix returned 355/354/360/363/356 against the original 353/354/356/357/355: the cross-run spread is up to 6 steps, comparable to the quoted 1-9 step effect, because those cells complete zero replicates and stop on the wall clock. The DIRECTION is sound (paired within one session on one node, 10 win / 0 loss over both A/Bs), so the sentence keeps "found shorter trees ... on all five seeds tried" and simply stops quoting a step count. The x3.5 wall figure is retained: it comes from the ship-gate run, which was not deadline-bound. NOT fixed here, deliberately: this entry and the corresponding roxygen still say `strategy = "thorough"` / `"large"`, and upstream has REMOVED `strategy` in favour of `effort`. That rewrite is deferred because the correct wording depends on an unmade decision about the gate's trigger (userSetHits vs an effort door) -- rewriting now would mean rewriting again. Both files auto-merge without a conflict marker, so this staleness is recorded in the campaign hub as part of the blocked rebase rather than left to be discovered. Co-Authored-By: Claude Opus 4.8 --- NEWS.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index e9dd9d066..f8f55b02a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -59,8 +59,8 @@ this applies under any scoring regime, though it was measured only under equal weights. It targets datasets big or difficult enough that an ordinary search stops short of the optimum: across 25 datasets spanning 20 to 4062 tips it - found shorter trees only on the 4062-tip matrix (on all five seeds tried, by - 1–9 steps), while from 20 to 173 tips it found trees of the same length and + found shorter trees only on the 4062-tip matrix (on all five seeds tried), + while from 20 to 173 tips it found trees of the same length and simply took about 3.5× as long — a cost incurred as extra work per replicate, not as slower convergence. Because most searches would pay for depth they do not need, it is offered only on that explicit signal and only on those two From 34cfc07aba0e276261d37973608e480d262783ae Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Mon, 3 Aug 2026 16:26:15 +0100 Subject: [PATCH 11/16] fix(reach): gate the deep-search bundle on a CALLER-SET targetHits The escalation triggered on `targetHits / defaultHits >= 2` alone, which upstream's `effort` ladder now satisfies by itself: .RungSpec()'s `hitMultiplier` doubles `targetHits` at rung 5 precisely WHEN THE USER DID NOT SET IT, so the ratio lands on exactly 2.0 and trips the gate. On any dataset over 120 tips .AutoRung() gives rung 4, so `effort = 1` alone fired a bundle that costs ~3.5x the per-replicate wall -- and rung 5 has already doubled `maxReplicates`, making one notch ~7x the work against docs that promise about 2x. .ApplyReachEscalation() now requires `userSetHits`. That axis is also measured flat by two independent lines: this feature's own A/B ties on every 125-, 131- and 173-tip cell, and the NA certify/effort panel (60 cells) found notch +2 matching notch +1 on every matrix while spending 798 replicates against 500. `.IwRatchetDepth()` still follows the ladder-raised value -- upstream states that is the point of the multiplier -- so only this bundle changes. Also corrects two claims in the evidence comment that my own fix to the deadline detector falsified, and which survived in R/ after being fixed in the FINDINGS: - "no cell in either arm hit its wall cap" is false. The engine stops at `maxSeconds * (1 - enumTimeFraction)`, i.e. 0.9x: 85 cells converged in both arms, 36/40 stopped at that deadline, and 4 are asymmetric. Both arms at one deadline is still a valid equal-wall comparison and all four asymmetric cells tied at 563, so no conclusion moves -- but the clean-sweep claim does not stand. - the "1-9 steps" margin is not an effect size. Re-running those cells moved scores several steps in both directions, because a deadline-truncated run completing zero replicates stops wherever the clock lands. The win/loss count replicates; the gap does not. Records the confirmation A/B (55 cells, shipped form, 0 asymmetric: 8 better / 1 worse / 46 tie) and states plainly that the effect is ONE matrix -- project4284 10/0 over both runs, project2771 4/2/4 (noise), nine matrices 0/0/45 -- and that project4284 won having completed zero replicates, so what pays is depth within a replicate, not restart volume. Docs: purges the removed `strategy = ` argument from the NEWS entry and the roxygen (it auto-merges with no conflict marker, so it needed catching by hand), including one pre-existing upstream entry that named it. Tests: new provenance test holding every VALUE equal -- rung 5 reaching targetHits = 20 by ladder vs a caller naming 20 -- so only the source of the number differs; the ladder run must match the rung below it and the caller's must deepen. Four end-to-end tests migrated off `strategy = `. Verified: 16/16 test_that blocks pass (0 failed, 0 error), plus test-iw-ratchet-depth.R 18/18, via temp-lib R CMD INSTALL. Co-Authored-By: Claude Opus 5 --- NEWS.md | 47 +++++++------ R/MaximizeParsimony.R | 96 ++++++++++++++++++++------ man/MaximizeParsimony.Rd | 36 +++++----- tests/testthat/test-reach-escalation.R | 86 +++++++++++++++++++---- 4 files changed, 194 insertions(+), 71 deletions(-) diff --git a/NEWS.md b/NEWS.md index ac643d539..32280a80c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -164,8 +164,8 @@ `nThreads > 1` it is evaluated when the coordinating thread polls, so it fires later and less predictably. -- Implied-weights searches under `strategy = "sprint"` or `"default"` now run a - deeper ratchet paid for by that flat patience: `sprint` takes +- Implied-weights searches at `effort` rung 1 (`sprint`) or 2 (`default`) now run + a deeper ratchet paid for by that flat patience: `sprint` takes `ratchetCycles = 12`, `ratchetPerturbProb = 0.25` and `stopPatience = 20`; `default` takes `ratchetCycles = 20` and `stopPatience = 15`. The two knobs ship together because each fails on its own — the deeper ratchet improves the @@ -182,24 +182,31 @@ profile parsimony are unchanged, as is `thorough`/`large`, and setting any of these fields yourself overrides all of it. -- Doubling `targetHits` or more, under `strategy = "thorough"` or `"large"`, now - also deepens the per-replicate perturbation itself, extending the existing - `targetHits` escalation beyond ratchet depth: more drifting, a larger - reweighting kick, a second sectorial pass after the ratchet, and internal - retention of near-optimal trees to fuse against. Unlike the ratchet deepening - this applies under any scoring regime, though it was measured only under equal - weights. It targets datasets big or difficult enough that an ordinary search - stops short of the optimum: across 25 datasets spanning 20 to 4062 tips it - found shorter trees only on the 4062-tip matrix (on all five seeds tried), - while from 20 to 173 tips it found trees of the same length and - simply took about 3.5× as long — a cost incurred as extra work per replicate, - not as slower convergence. Because most searches would pay for depth they do - not need, it is offered only on that explicit signal and only on those two - presets; `sprint` and `default` keep the implied-weights operating point - described above, `ratchetCycles` remains governed by the implied-weights - ratchet deepening, and any control field you set yourself is preserved. Note - that the documented large-`targetHits` idiom for collecting the full set of - most-parsimonious trees also engages this on those two presets. +- Setting `targetHits` yourself to at least twice its default, at `effort` rung 3 + (`thorough`) or above, now also deepens the per-replicate perturbation itself, + extending the existing `targetHits` escalation beyond ratchet depth: more + drifting, a larger reweighting kick, a second sectorial pass after the ratchet, + and internal retention of near-optimal trees to fuse against. Unlike the + ratchet deepening this applies under any scoring regime, though it was measured + only under equal weights. It targets datasets big or difficult enough that an + ordinary search stops short of the optimum: across 25 datasets spanning 20 to + 4062 tips it found shorter trees only on the 4062-tip matrix — there on all ten + seeds tried across two runs — while from 20 to 173 tips it found trees of the + same length and simply took about 3.5× as long, a cost incurred as extra work + per replicate rather than as slower convergence. Because most searches would + pay for depth they do not need, it is offered only on that explicit signal and + only at those rungs; `sprint` and `default` keep the implied-weights operating + point described above, `ratchetCycles` remains governed by the implied-weights + ratchet deepening, and any control field you set yourself is preserved. + + It follows a `targetHits` that *you* set, and only that. Raising `effort` also + raises `targetHits` from rung 5, but that is a change of budget rather than a + statement about the dataset, and deepening the perturbation on top of it was + measured to buy nothing: every 125-, 131- and 173-tip cell tied, and a separate + 60-cell panel found `effort = 2` matching `effort = 1` on every matrix while + spending 798 replicates against 500. Note that the documented + large-`targetHits` idiom for collecting the full set of most-parsimonious trees + does engage this, since you set the number. - Fixed: a large `targetHits` combined with a large `perturbStopFactor` stopped the search after two replicates and silently returned a worse tree. The diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index a38060aed..577b34b07 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -288,14 +288,24 @@ # # EVIDENCE (general-pool A/B, 2026-07-28; 25 training matrices x 5 seeds = 125 # cells; both arms ran at the SAME raised `targetHits`, so only these levers -# differ; no cell in either arm hit its wall cap, so the comparison is not a -# budget artefact). Strict paired final-score wins, by size tier: +# differ). Budget regime, re-derived once the deadline detector was fixed (the +# engine stops at `maxSeconds * (1 - enumTimeFraction)`, i.e. 0.9x by default, +# NOT at `maxSeconds`): 85 cells converged in both arms, 36 base and 40 deltas +# cells stopped at that deadline, and 4 cells are ASYMMETRIC (project2184, +# deltas at the deadline and base not). Both arms stopping at one shared +# deadline is a valid equal-wall comparison; only the asymmetric cells are +# suspect, and all four of those tied at 563, so no claim below rests on one. +# Strict paired final-score wins, by size tier: # small (n=35): 0 better, 0 worse, 35 tie # medium (n=35): 0 better, 0 worse, 35 tie # large (n=35): 1 better, 1 worse, 33 tie (a wash) # xlarge (n=20): 5 better, 0 worse, 15 tie # Read that last row carefully: ALL FIVE wins are the SAME matrix, project4284 -# at 4062 tips, which improved on every one of its five seeds (by 1-9 steps). +# at 4062 tips, which improved on every one of its five seeds. The MARGIN is +# not a quotable effect size: re-running those cells moved scores by several +# steps in both directions, because a deadline-truncated run that completes zero +# replicates stops wherever the clock lands. The win/loss COUNT replicates; the +# gap does not. # The other three xlarge matrices (125, 131 and 173 tips) all tied. So the # demonstrated benefit is NOT "datasets over 120 tips" -- it is datasets far too # large to converge within an ordinary budget, plus the hard-reach tail (on the @@ -304,6 +314,16 @@ # Cost: median total wall x3.56 and time-to-best x2.49, but replicates-to-best # x1.00 -- the wall gap is entirely cost-PER-replicate (~15x candidates # evaluated), not slower convergence. +# +# CONFIRMED in exactly the form shipped here (2026-07-30; 11 large/xlarge +# matrices x 5 seeds = 55 cells, the six levers below and nothing else, 0 +# asymmetric cells): paired 8 better / 1 worse / 46 tie, no tier regression. +# The effect is still ONE matrix, though. Pooling both runs: project4284 is +# 10 win / 0 loss over 10 seeds, project2771 is 4 win / 2 loss / 4 tie (noise), +# and the other nine matrices are 0 win / 0 loss / 45 tie. Note too that +# project4284 won having completed ZERO replicates, so what pays there is depth +# WITHIN one replicate -- the opposite lever to the restart volume that a raised +# `maxReplicates` buys, and a reason not to expect these to compose. # So over the tested 20-173 tip range these buy no reach and cost ~3.5x the wall # (the sample jumps from 173 tips straight to 4062, so the range in between is # untested), which is exactly why they are kept OUT of `.StrategyPresets()`: as @@ -334,14 +354,33 @@ # Apply the deeper-perturbation bundle, preserving every field the caller set. # `userSet` names those fields, exactly as for `.IwRatchetDepth`. # +# `userSetHits` must be TRUE: this fires only where the CALLER named +# `targetHits`, never where the effort ladder raised it. The ratio alone will +# not do, because .RungSpec()'s `hitMultiplier` doubles `targetHits` at rung 5 +# precisely WHEN THE USER DID NOT SET IT -- landing the ratio on exactly 2 and +# tripping the `>=` below on a mechanical ladder artefact, on any dataset over +# 120 tips, where `effort = 1` reaches rung 5 unaided. That is not what the A/B +# measured; it measured a caller asking for a deeper search than the defaults. +# +# Riding the ladder would also mean riding a signal that is measured FLAT. Two +# independent lines say the returns die above rung 4: the A/B above ties on every +# 125-, 131- and 173-tip cell, and the NA certify/effort panel (60 cells) found +# notch +2 matching notch +1 on every matrix while burning 798 replicates against +# 500. Rung 5 has already doubled `maxReplicates`; layering ~3.5x the +# per-replicate cost on top would make one notch roughly 7x the work for no +# measured gain, against documentation that promises about 2x. +# # Scoped to `thorough`/`large`, matching `.IwRatchetDepth`: those are the presets # the A/B's benefit came from (auto selects `large` for the 4062-tip matrix and # `thorough` for project5432), and on smaller data the same A/B measured 0 better # / 0 worse across 70 small- and medium-tier cells -- pure wall cost. Escalating # `sprint`/`default` would also contradict their documented character ("Fast # search: 3 ratchet cycles, no drift"), so they are left alone. -.ApplyReachEscalation <- function(control, strategy, escalation, +.ApplyReachEscalation <- function(control, strategy, escalation, userSetHits, userSet = character(0)) { + if (!isTRUE(userSetHits)) { + return(control) + } if (!length(strategy) || !strategy %in% c("thorough", "large")) { return(control) } @@ -911,23 +950,27 @@ #' were slower to the optimum on every matrix tested), and setting #' `ratchetCycles` yourself overrides this entirely. #' -#' Doubling `targetHits` or more goes one step further and, under -#' `strategy = "thorough"` or `"large"`, also deepens the per-replicate -#' perturbation itself: more drifting, a larger reweighting kick, a second -#' sectorial pass after the ratchet, and (internally) retention of near-optimal -#' trees to fuse against. Unlike the ratchet deepening above this applies under -#' any scoring regime, though it was measured only under equal weights. It is -#' aimed at datasets big or difficult enough that an ordinary search stops short -#' of the optimum: in testing it found shorter trees only on a 4062-tip matrix -#' (on all five seeds tried), while from 20 to 173 tips it found trees of the -#' same length and simply took around 3.5 times as long -- so it is offered on -#' this explicit signal rather than enabled by default. Any of these values you -#' set yourself is left untouched, and the trees returned are still only the -#' best found. +#' Setting `targetHits` yourself to at least twice its default goes one step +#' further and, at `effort` rung 3 (`thorough`) or above, also deepens the +#' per-replicate perturbation itself: more drifting, a larger reweighting kick, +#' a second sectorial pass after the ratchet, and (internally) retention of +#' near-optimal trees to fuse against. Unlike the ratchet deepening above, this +#' applies under any scoring regime, though it was measured only under equal +#' weights. It is aimed at datasets big or difficult enough that an ordinary +#' search stops short of the optimum: in testing it found shorter trees on one +#' 4062-tip matrix, on all ten seeds tried across two runs, while from 20 to 173 +#' tips it found trees of the same length and simply took around 3.5 times as +#' long -- so it is offered on this explicit signal rather than enabled by +#' default. It follows a `targetHits` that *you* set, and only that: raising +#' `effort` also raises `targetHits` from rung 5, but that is a change of budget +#' rather than a statement about the dataset, and deepening the perturbation on +#' top of it was measured to buy nothing. Any of these values you set yourself +#' is left untouched, and the trees returned are still only the best found. +#' #' Note that the large-`targetHits` idiom for collecting the full set of -#' most-parsimonious trees, above, therefore also engages this deeper search on -#' those two presets; set `driftCycles`, `intraFuse` and so on yourself if you -#' want the wider sampling without the extra per-replicate cost. +#' most-parsimonious trees, above, therefore also engages this deeper search at +#' those rungs; set `driftCycles`, `intraFuse` and so on yourself if you want +#' the wider sampling without the extra per-replicate cost. #' @param maxSeconds Numeric: maximum wall-clock time in seconds for the #' search. When reached, the current replicate finishes and the search #' stops. `0` (default) means no time limit. @@ -1347,8 +1390,8 @@ MaximizeParsimony <- function( control[[.f]] <- iwStop[[.f]] } - # The same `targetHits` signal, one step further: a caller who has at least - # DOUBLED it has asked to keep searching well past ordinary convergence, so + # The same `targetHits` signal, one step further: a caller who has THEMSELVES + # at least doubled it has asked to keep searching past ordinary convergence, so # also deepen the per-replicate perturbation itself (drift, the auto kick, # a post-ratchet sectorial re-search, near-optimal pool retention). Unlike # the ratchet depth above this applies under any scorer, but it is not @@ -1361,9 +1404,15 @@ MaximizeParsimony <- function( # implied weights) and `.IwStopPackage` (sprint/default, implied weights) # never both fire, and this bundle is scoped to thorough/large so it # cannot disturb the sprint/default operating point measured above. + # `userSetHits` is passed on deliberately: the rung's `hitMultiplier` above + # raised `targetHits` only when the caller did NOT set it, so the ratio on + # its own cannot distinguish "search harder, please" from a rung change. + # .IwRatchetDepth() above *should* follow the ladder-raised value -- that is + # stated to be the point of the multiplier -- and this bundle should not. escalation <- .TargetHitsEscalation(targetHits, defaultHits) poolBefore <- control[["poolSuboptimal"]] control <- .ApplyReachEscalation(control, strategy, escalation, + userSetHits = userSetHits, userSet = userSet) # `poolSuboptimal` is raised here only as an internal aid: `intraFuse` needs # suboptimal recipients to fuse against. It must not leak into the RESULT -- @@ -1373,7 +1422,8 @@ MaximizeParsimony <- function( # than compared later) so a caller's own `poolSuboptimal` is untouched: they # asked for those trees and still get them. escalatedPool <- !identical(control[["poolSuboptimal"]], poolBefore) - if (verbosity >= 1L && escalation >= .reachEscalationMinRatio && + if (verbosity >= 1L && userSetHits && + escalation >= .reachEscalationMinRatio && strategy %in% c("thorough", "large")) { cli::cli_alert_info( "Deep search: {.field targetHits} raised {round(escalation, 1)}x" diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index d953c8a86..316726a4d 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -267,23 +267,27 @@ does not make the ratchet shallower than its default depth (fewer cycles were slower to the optimum on every matrix tested), and setting \code{ratchetCycles} yourself overrides this entirely. -Doubling \code{targetHits} or more goes one step further and, under -\code{strategy = "thorough"} or \code{"large"}, also deepens the per-replicate -perturbation itself: more drifting, a larger reweighting kick, a second -sectorial pass after the ratchet, and (internally) retention of near-optimal -trees to fuse against. Unlike the ratchet deepening above this applies under -any scoring regime, though it was measured only under equal weights. It is -aimed at datasets big or difficult enough that an ordinary search stops short -of the optimum: in testing it found shorter trees only on a 4062-tip matrix -(on all five seeds tried), while from 20 to 173 tips it found trees of the -same length and simply took around 3.5 times as long -- so it is offered on -this explicit signal rather than enabled by default. Any of these values you -set yourself is left untouched, and the trees returned are still only the -best found. +Setting \code{targetHits} yourself to at least twice its default goes one step +further and, at \code{effort} rung 3 (\code{thorough}) or above, also deepens the +per-replicate perturbation itself: more drifting, a larger reweighting kick, +a second sectorial pass after the ratchet, and (internally) retention of +near-optimal trees to fuse against. Unlike the ratchet deepening above, this +applies under any scoring regime, though it was measured only under equal +weights. It is aimed at datasets big or difficult enough that an ordinary +search stops short of the optimum: in testing it found shorter trees on one +4062-tip matrix, on all ten seeds tried across two runs, while from 20 to 173 +tips it found trees of the same length and simply took around 3.5 times as +long -- so it is offered on this explicit signal rather than enabled by +default. It follows a \code{targetHits} that \emph{you} set, and only that: raising +\code{effort} also raises \code{targetHits} from rung 5, but that is a change of budget +rather than a statement about the dataset, and deepening the perturbation on +top of it was measured to buy nothing. Any of these values you set yourself +is left untouched, and the trees returned are still only the best found. + Note that the large-\code{targetHits} idiom for collecting the full set of -most-parsimonious trees, above, therefore also engages this deeper search on -those two presets; set \code{driftCycles}, \code{intraFuse} and so on yourself if you -want the wider sampling without the extra per-replicate cost.} +most-parsimonious trees, above, therefore also engages this deeper search at +those rungs; set \code{driftCycles}, \code{intraFuse} and so on yourself if you want +the wider sampling without the extra per-replicate cost.} \item{maxSeconds}{Numeric: maximum wall-clock time in seconds for the search. When reached, the current replicate finishes and the search diff --git a/tests/testthat/test-reach-escalation.R b/tests/testthat/test-reach-escalation.R index 176bd0b6a..1f6026272 100644 --- a/tests/testthat/test-reach-escalation.R +++ b/tests/testthat/test-reach-escalation.R @@ -1,4 +1,6 @@ -# Deep-search escalation: raising `targetHits` deepens per-replicate perturbation. +# Deep-search escalation: a CALLER-SET `targetHits` at 2x its default deepens +# per-replicate perturbation. Caller-set is half the contract: the effort ladder +# raises `targetHits` itself from rung 5, and that must NOT engage the bundle. library("TreeTools", quietly = TRUE) data("inapplicable.phyData", package = "TreeSearch") ds <- inapplicable.phyData[["Vinther2008"]] # 23 tips @@ -43,7 +45,8 @@ test_that("escalation does NOT touch ratchetCycles", { # carry ratchetCycles. expect_false("ratchetCycles" %in% names(TreeSearch:::.ReachEscalationDeltas())) ctrl <- TreeSearch:::.ApplyReachEscalation(SearchControl(ratchetCycles = 48L), - "thorough", escalation = 10) + "thorough", escalation = 10, + userSetHits = TRUE) expect_identical(ctrl[["ratchetCycles"]], 48L) }) @@ -57,7 +60,8 @@ test_that("the two escalations compose without fighting over ratchetCycles", { targetHits = 20L, defaultHits = 10L) expect_equal(iw, 96L) # 48 * escalation(2), under the 115 cap ctrl[["ratchetCycles"]] <- iw - ctrl <- TreeSearch:::.ApplyReachEscalation(ctrl, "thorough", escalation = 2) + ctrl <- TreeSearch:::.ApplyReachEscalation(ctrl, "thorough", escalation = 2, + userSetHits = TRUE) expect_identical(ctrl[["ratchetCycles"]], 96L) # NOT clobbered by the bundle expect_identical(ctrl[["driftCycles"]], 25L) # bundle still applied }) @@ -67,7 +71,8 @@ test_that(".ApplyReachEscalation applies all deltas at or above the ratio", { for (strat in c("thorough", "large")) { for (esc in c(2, 2.5, 20)) { ctrl <- TreeSearch:::.ApplyReachEscalation(SearchControl(), strat, - escalation = esc) + escalation = esc, + userSetHits = TRUE) for (nm in names(deltas)) { expect_identical(ctrl[[nm]], deltas[[nm]], info = paste(strat, nm)) } @@ -79,17 +84,41 @@ test_that(".ApplyReachEscalation is inert below the ratio", { stock <- SearchControl() for (esc in c(1, 1.5, 1.99)) { expect_identical(TreeSearch:::.ApplyReachEscalation(stock, "thorough", - escalation = esc), + escalation = esc, + userSetHits = TRUE), stock) } # Degenerate escalation must not escalate. for (esc in list(NA_real_, numeric(0), Inf)) { expect_identical(TreeSearch:::.ApplyReachEscalation(stock, "thorough", - escalation = esc), + escalation = esc, + userSetHits = TRUE), stock) } }) +test_that(".ApplyReachEscalation requires the CALLER to have set targetHits", { + # The ratio alone cannot distinguish "search harder" from a rung change: + # .RungSpec()'s hitMultiplier doubles `targetHits` at rung 5 exactly when the + # user did NOT set it, landing the ratio on 2.0 and tripping the `>=` gate. + # Two independent lines measure that axis flat (see .ApplyReachEscalation), so + # the bundle must stay shut unless the caller named the number themselves. + stock <- SearchControl() + for (notSet in list(FALSE, NA, NULL, logical(0))) { + expect_identical( + TreeSearch:::.ApplyReachEscalation(stock, "thorough", escalation = 10, + userSetHits = notSet), + stock, info = paste("userSetHits", format(notSet)) + ) + } + # ... and open when they did, at the same ratio. + expect_identical( + TreeSearch:::.ApplyReachEscalation(stock, "thorough", escalation = 10, + userSetHits = TRUE)[["driftCycles"]], + 25L + ) +}) + test_that(".ApplyReachEscalation is scoped to thorough/large", { # `sprint` and `default` document themselves as fast/shallow ("3 ratchet # cycles, no drift"), and the A/B measured 0 better / 0 worse across their @@ -97,7 +126,8 @@ test_that(".ApplyReachEscalation is scoped to thorough/large", { stock <- SearchControl() for (strat in c("sprint", "default", "none", NA_character_, character(0))) { expect_identical( - TreeSearch:::.ApplyReachEscalation(stock, strat, escalation = 10), + TreeSearch:::.ApplyReachEscalation(stock, strat, escalation = 10, + userSetHits = TRUE), stock, info = paste("strategy", strat) ) } @@ -105,7 +135,7 @@ test_that(".ApplyReachEscalation is scoped to thorough/large", { test_that(".ApplyReachEscalation preserves caller-set fields", { ctrl <- TreeSearch:::.ApplyReachEscalation( - SearchControl(), "thorough", escalation = 4, + SearchControl(), "thorough", escalation = 4, userSetHits = TRUE, userSet = c("driftCycles", "intraFuse") ) expect_identical(ctrl[["driftCycles"]], SearchControl()[["driftCycles"]]) @@ -130,9 +160,11 @@ test_that("escalation measurably deepens the search end to end", { # sectorial pass), so at matched replicates the escalated run must evaluate # substantially more candidates. Vinther2008 (23 tips): default targetHits = 10, # so 20 is exactly 2x. + # Vinther2008 is 23 tips, so .AutoRung() gives rung 1 (`sprint`); `effort = 2` + # is rung 3 (`thorough`), which is in scope and below the rung-5 hitMultiplier. runCand <- function(hits) { set.seed(4242) - r <- MaximizeParsimony(ds, strategy = "thorough", maxReplicates = 2L, + r <- MaximizeParsimony(ds, effort = 2L, maxReplicates = 2L, targetHits = hits, maxSeconds = 0, verbosity = 0L) list(cand = as.double(attr(r, "candidates_evaluated")), res = r) } @@ -157,7 +189,7 @@ test_that("sprint is NOT escalated end to end", { # EXACTLY. This fails the moment the strategy gate is loosened. runCand <- function(hits) { set.seed(99L) - r <- MaximizeParsimony(ds, strategy = "sprint", maxReplicates = 2L, + r <- MaximizeParsimony(ds, effort = 0L, maxReplicates = 2L, targetHits = hits, maxSeconds = 0, verbosity = 0L) as.double(attr(r, "candidates_evaluated")) } @@ -170,7 +202,7 @@ test_that("an escalated search still returns only best-score trees", { # verbatim, so without the guard the caller would silently get trees up to 3 # steps worse than attr(, "score") from a result documented as the best found. set.seed(31L) - r <- MaximizeParsimony(ds, strategy = "thorough", maxReplicates = 3L, + r <- MaximizeParsimony(ds, effort = 2L, maxReplicates = 3L, targetHits = 20L, collapse = FALSE, verbosity = 0L) best <- attr(r, "score") expect_true(is.finite(best)) @@ -178,11 +210,41 @@ test_that("an escalated search still returns only best-score trees", { expect_true(all(scores == best)) }) +test_that("the effort ladder's own targetHits rise does NOT deepen the search", { + # The provenance test, with every VALUE held equal. At rung 5 the ladder + # doubles the 23-tip default of 10 to 20 by itself; the second run names 20, + # so the ladder skips its multiplier and leaves it at 20. Both runs therefore + # search with targetHits = 20, the same preset (`large`), the same replicate + # cap and the same seed -- the ONLY difference is who set the number. Only the + # caller's version may deepen the perturbation. With the gate reading the + # ratio alone (as it first did) both runs escalate and the counts match, which + # would fire this bundle on every dataset over 120 tips at `effort = 1`. + runCand <- function(...) { + set.seed(808L) + r <- MaximizeParsimony(ds, effort = 4L, maxReplicates = 2L, + maxSeconds = 0, verbosity = 0L, ...) + as.double(attr(r, "candidates_evaluated")) + } + ladder <- runCand() + asked <- runCand(targetHits = 20L) + expect_true(is.finite(ladder) && ladder > 0) + expect_gt(asked, 1.2 * ladder) + # And the ladder run is indistinguishable from the rung below it, whose ratio + # is 1: rung 5 with a user-set `maxReplicates` differs only in the hit target. + set.seed(808L) + rung4 <- as.double(attr( + MaximizeParsimony(ds, effort = 3L, maxReplicates = 2L, maxSeconds = 0, + verbosity = 0L), + "candidates_evaluated" + )) + expect_identical(ladder, rung4) +}) + test_that("a caller's own poolSuboptimal is still honoured", { # The guard above must not steal the documented behaviour from someone who # asked for suboptimal trees themselves. set.seed(31L) - r <- MaximizeParsimony(ds, strategy = "thorough", maxReplicates = 3L, + r <- MaximizeParsimony(ds, effort = 2L, maxReplicates = 3L, targetHits = 20L, poolSuboptimal = 3, collapse = FALSE, verbosity = 0L) expect_s3_class(r, "multiPhylo") From 088f6c985dd03844bb633f2964a4536116db834e Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Mon, 3 Aug 2026 16:27:27 +0100 Subject: [PATCH 12/16] docs(reach): the effort-ladder blocker is resolved; no effort door Records the decision and, more importantly, WHY the open question closed NO rather than by preference: notch +2 is measured inert on a second, independent panel (60 cells, NA certify/effort), which corroborates this study's own 125-/131-/173-tip ties from a different regime. Co-Authored-By: Claude Opus 5 --- dev/benchmarks/reach_escalation_FINDINGS.md | 24 ++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md index cfc672e19..0f04b76f5 100644 --- a/dev/benchmarks/reach_escalation_FINDINGS.md +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -190,7 +190,28 @@ at 5–9 s of a 720 s cap, project4359 stopping on `targetHits` at 28 replicates cost information and must not be counted as evidence of cost-neutrality. The ×3.56 wall figure belongs to the ship-gate run, which was not deadline-bound. -## 🚨 The upstream API change breaks the gate's premise (2026-07-31, `419168d4`) +## ✅ RESOLVED — the upstream API change broke the gate's premise (`419168d4`) + +**Fixed 2026-08-03 in `34cfc07a`: option (c), `userSetHits`, with no `effort` door.** +`.ApplyReachEscalation()` now takes `userSetHits` and returns `control` untouched unless it is +`TRUE`. The open question below — whether high `effort` should open a second door — is +**answered NO**, and by evidence rather than taste: notch +2 is measured INERT (NA +certify/effort panel, `1958f211`, 60 cells: +2 matched +1 on every matrix while spending 798 +replicates against 500; *"further tuning above rung 4 is measured flat"*), which corroborates +this study's own 125-/131-/173-tip ties from a different regime. Two independent lines now say +the returns die after the first notch, so an `effort` trigger would buy nothing at ~3.5× wall. + +Guarded by a provenance test that holds every *value* equal: rung 5 reaching `targetHits = 20` +via the ladder, against a caller naming 20. Same preset, same seed, same replicate cap, same +hit target — only the *source* of the number differs. The ladder run must match the rung below +it; the caller's must deepen. Reading the ratio alone makes the two identical — the bug exactly. + +`.IwRatchetDepth()` deliberately still follows the ladder-raised value: upstream's own comment +states that coupling is the point of the multiplier. Only this bundle changed. + +The analysis that found it, kept because the mechanism is the reusable part: + +### 🚨 The mechanism (2026-07-31) `origin/cpp-search` **removed `strategy`** and replaced it with `effort = 0L`, a *relative* offset on an internal ladder (`.effortLadder = sprint, default, thorough, large`, then rungs @@ -228,6 +249,7 @@ Options, with the recommendation: Open question for the maintainer: at high `effort` (say ≥ +3) the user arguably *does* want the bundle. If so, the trigger becomes `userSetHits && ratio >= 2` **or** an explicit `effort` threshold — a deliberate second door rather than an accident of the multiplier. +**Answered NO, 2026-08-03 — see the resolution at the head of this section.** ## Why it is gated, and gated to `thorough`/`large` From 00ea838253629be236c64355e8125b8c374ad977 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Tue, 4 Aug 2026 07:10:11 +0100 Subject: [PATCH 13/16] bench(t253): probe to decide annotate-vs-retract on the gap analysis `t253_conv_gap_mbank.csv` and `t253_gap_characterization.md` publish an n=23 Spearman analysis derived from t252 cells, and t252 ran on an engine whose Wagner addition was several-fold worse (project4284: bare AdditionTree, same data and seed, no search -- 1590 March vs 409 current). If that gap is broad, every t252 start tree is systematically bad and the analysis must be RETRACTED; if it is project4284's alone, one row needs ANNOTATING. Nothing on record distinguishes those, so this measures bare `AdditionTree` on all 25 MBANK_FIXED_SAMPLE matrices under both engines. Three confounds are designed out, each of which would fake a result: - PREPROCESSING: the libraries ship different TreeTools versions, so `prep` builds every dataset ONCE and saves RDS; both engines read the same object. - THE SCORER: `score` re-scores every tree under ONE library, so numbers differ only by tree quality -- not by scorer version. - GAP TREATMENT: both arms map "-" -> "?", since the two conventions differ by +61 on one tree and would dwarf the effect being measured. Reporting is per MATRIX (medians over 3 seeds), not per (matrix, seed): seeds within a matrix are not independent evidence about an engine. The ratio is deliberately OLD/NEW so that "> 1" means "March built a longer tree", which is what every caption asserts. Alphabetical engine ordering put `cur` in the numerator and inverted every count against its own caption -- caught in the smoke test, and worth the explicit comment it now carries. Smoke-tested end to end locally on 3 built-in matrices, including a same-engine control that returns ratio 1.000 on all three (which also confirms the seeding is deterministic and the scorer is applied consistently). NOT YET RUN on Hamilton: the VPN dropped mid-session. Co-Authored-By: Claude Opus 5 --- dev/benchmarks/t253_wagner_era_probe.R | 239 ++++++++++++++++++++++++ dev/benchmarks/t253_wagner_era_probe.sh | 55 ++++++ 2 files changed, 294 insertions(+) create mode 100644 dev/benchmarks/t253_wagner_era_probe.R create mode 100644 dev/benchmarks/t253_wagner_era_probe.sh diff --git a/dev/benchmarks/t253_wagner_era_probe.R b/dev/benchmarks/t253_wagner_era_probe.R new file mode 100644 index 000000000..bb7978b77 --- /dev/null +++ b/dev/benchmarks/t253_wagner_era_probe.R @@ -0,0 +1,239 @@ +#!/usr/bin/env Rscript +# t253 DECISION PROBE: is the March->current Wagner improvement GENERAL, or is it +# project4284's alone? +# +# WHY THIS EXISTS. `dev/benchmarks/t252_mbank_*.csv` (2026-03-27) were produced by +# an engine whose Wagner addition was several-fold worse than the current one: on +# project4284, bare `AdditionTree` with the same data and seed and NO SEARCH scores +# 1590 under the March library against 409 under the current one (worse AND 5x +# faster -- the signature of the insertion-cost bug). `t253_conv_gap_mbank.csv` and +# `t253_gap_characterization.md` publish an n=23 Spearman analysis derived from those +# t252 cells, so the question is not academic: +# +# * if the improvement is BROAD, every t252 start tree is systematically bad and +# the t253 analysis rests on them -> RETRACT it. +# * if it is essentially project4284 alone, one row is contaminated -> ANNOTATE. +# +# Nothing else decides this. A search-level comparison would confound the addition +# tree with everything the search does afterwards, which is why this measures BARE +# `AdditionTree` and nothing else. +# +# THREE CONFOUNDS THIS DESIGN REMOVES, each of which would fake a result: +# +# 1. PREPROCESSING. The two libraries ship different TreeTools versions, so running +# `ReadAsPhyDat` + `PhyDatToMatrix` under each would let a preprocessing +# difference masquerade as an engine difference. STEP `prep` therefore builds +# every dataset ONCE, under one library, and saves it as RDS; both engines read +# the identical object. +# 2. THE SCORER. Asking each engine to score its own tree compares scorers as well +# as builders. STEP `score` re-scores EVERY saved tree under ONE library, so the +# numbers differ only by tree quality. Cf. na-validation-alignment-gotcha. +# 3. GAP TREATMENT. t252 kept "-" as a sixth level (BGS); this maps "-" -> "?" +# (plain Fitch) for BOTH arms, matching `reach_escalation_ab.R`. That is worth +# +61 on one tree, so mixing the two conventions across arms would dwarf the +# effect being measured. Both arms use the SAME convention, which is what makes +# this an engine comparison. +# +# Env: STEP (prep|addition|score), OUT_DIR, ENGINE (label, for `addition`), +# NEOTRANS_DIR + CAT_CSV (for `prep`), N_SEEDS (default 3). +# The library selection is the CALLER's job, via R_LIBS -- see the .sh. + +step <- Sys.getenv("STEP", "") +outDir <- Sys.getenv("OUT_DIR", "") +engine <- Sys.getenv("ENGINE", "") +nSeeds <- as.integer(Sys.getenv("N_SEEDS", "3")) +if (!nzchar(outDir)) stop("OUT_DIR unset") +dir.create(outDir, showWarnings = FALSE, recursive = TRUE) +pdDir <- file.path(outDir, "pd") +treeDir <- file.path(outDir, "trees") + +# The fixed 25-matrix training sample, verbatim from bench_datasets.R / +# reach_escalation_ab.R: "results are only comparable when the same sample is used". +MBANK_FIXED_SAMPLE <- c( + "project532", "project2346", "project2451", "project4501", + "project944", "project971_(1)", "project2762", + "project826", "project561", "project571", "project4146_(3)", + "project3688", "project4049", "project423", + "project4286", "project4359", "project4397", "project2084_(1)", + "project2771", "project2184", "project3938", + "syab07201", "project4133", "project804", "project4284" +) + +BASE_SEED <- 1301L +seedsFor <- function(i) BASE_SEED + seq_len(nSeeds) - 1L + (i - 1L) * 100L + +safeKey <- function(k) gsub("[^A-Za-z0-9]", "_", k) + +# ---------------------------------------------------------------- STEP: prep ---- +if (identical(step, "prep")) { + suppressPackageStartupMessages(library("TreeTools")) + neoDir <- Sys.getenv("NEOTRANS_DIR", "") + catCsv <- Sys.getenv("CAT_CSV", "") + if (!nzchar(neoDir)) stop("NEOTRANS_DIR unset") + if (!nzchar(catCsv)) stop("CAT_CSV unset") + catalogue <- read.csv(catCsv, stringsAsFactors = FALSE) + rownames(catalogue) <- catalogue$key + dir.create(pdDir, showWarnings = FALSE, recursive = TRUE) + + toFitch <- function(pd) { + m <- PhyDatToMatrix(pd, ambigNA = FALSE) + m[m == "-"] <- "?" + MatrixToPhyDat(m) + } + + rows <- list() + for (key in MBANK_FIXED_SAMPLE) { + if (!key %in% catalogue$key) { + message("SKIP (not in catalogue): ", key) + next + } + row <- catalogue[key, ] + # SEQUESTER the validation split: a one-way door. + if (!identical(row$split, "training")) { + message("SKIP (split=", row$split, ", SEQUESTERED): ", key) + next + } + f <- file.path(neoDir, row$filename) + if (!file.exists(f)) { + message("SKIP (file missing): ", key) + next + } + pd <- toFitch(suppressWarnings(ReadAsPhyDat(f))) + saveRDS(pd, file.path(pdDir, paste0(safeKey(key), ".rds"))) + rows[[length(rows) + 1L]] <- data.frame( + key = key, nTip = length(pd), nChar = attr(pd, "nr"), + stringsAsFactors = FALSE + ) + cat(sprintf("prep %-18s nTip=%5d nChar=%4d\n", key, length(pd), attr(pd, "nr"))) + } + write.csv(do.call(rbind, rows), file.path(outDir, "manifest.csv"), + row.names = FALSE) + cat("prep done:", length(rows), "datasets\n") +} + +# ------------------------------------------------------------ STEP: addition ---- +# Bare `AdditionTree`, no search. Trees are saved as RDS rather than Newick: a +# round-trip through Newick could renumber or reorder, and STEP `score` must see the +# tree the engine actually produced. +if (identical(step, "addition")) { + if (!nzchar(engine)) stop("ENGINE unset") + suppressPackageStartupMessages(library("TreeSearch")) + dir.create(treeDir, showWarnings = FALSE, recursive = TRUE) + tsVer <- as.character(utils::packageVersion("TreeSearch")) + tsLib <- dirname(system.file(package = "TreeSearch")) + cat("engine=", engine, " TreeSearch ", tsVer, " from ", tsLib, "\n", sep = "") + + manifest <- read.csv(file.path(outDir, "manifest.csv"), stringsAsFactors = FALSE) + rows <- list() + for (i in seq_len(nrow(manifest))) { + key <- manifest$key[i] + pd <- readRDS(file.path(pdDir, paste0(safeKey(key), ".rds"))) + for (sd in seedsFor(i)) { + set.seed(sd) + t0 <- Sys.time() + tr <- tryCatch(AdditionTree(pd), error = function(e) e) + wall <- as.double(difftime(Sys.time(), t0, units = "secs")) + if (inherits(tr, "error")) { + cat(sprintf("FAIL %-18s seed=%d: %s\n", key, sd, conditionMessage(tr))) + rows[[length(rows) + 1L]] <- data.frame( + key = key, engine = engine, tsVersion = tsVer, seed = sd, + wallS = wall, ok = FALSE, stringsAsFactors = FALSE + ) + next + } + # Write the deliverable BEFORE anything that could fail on it: a verification + # error must not take the computed tree with it (the 2026-07-31 lesson). + saveRDS(tr, file.path(treeDir, sprintf("%s__%s__s%d.rds", + safeKey(key), engine, sd))) + rows[[length(rows) + 1L]] <- data.frame( + key = key, engine = engine, tsVersion = tsVer, seed = sd, + wallS = wall, ok = TRUE, stringsAsFactors = FALSE + ) + cat(sprintf("add %-18s %-6s seed=%d %.2fs\n", key, engine, sd, wall)) + } + } + write.csv(do.call(rbind, rows), + file.path(outDir, paste0("addition_", engine, ".csv")), + row.names = FALSE) +} + +# --------------------------------------------------------------- STEP: score ---- +# ONE scorer for every tree, whichever engine built it. +if (identical(step, "score")) { + suppressPackageStartupMessages(library("TreeSearch")) + cat("scorer: TreeSearch ", as.character(utils::packageVersion("TreeSearch")), + "\n", sep = "") + manifest <- read.csv(file.path(outDir, "manifest.csv"), stringsAsFactors = FALSE) + files <- list.files(treeDir, pattern = "\\.rds$", full.names = TRUE) + rows <- list() + for (f in files) { + parts <- strsplit(sub("\\.rds$", "", basename(f)), "__", fixed = TRUE)[[1]] + if (length(parts) != 3L) { + cat("SKIP unparseable filename:", basename(f), "\n") + next + } + keySafe <- parts[[1]] + idx <- match(keySafe, safeKey(manifest$key)) + if (is.na(idx)) { + cat("SKIP no manifest row:", basename(f), "\n") + next + } + key <- manifest$key[idx] + pd <- readRDS(file.path(pdDir, paste0(keySafe, ".rds"))) + tr <- readRDS(f) + sc <- tryCatch(TreeLength(tr, pd, concavity = Inf), + error = function(e) { + cat("SCORE FAIL", basename(f), ":", conditionMessage(e), "\n") + NA_real_ + }) + rows[[length(rows) + 1L]] <- data.frame( + key = key, engine = parts[[2]], seed = as.integer(sub("^s", "", parts[[3]])), + nTip = manifest$nTip[idx], nChar = manifest$nChar[idx], + score = as.double(sc), stringsAsFactors = FALSE + ) + } + scores <- do.call(rbind, rows) + write.csv(scores, file.path(outDir, "scores.csv"), row.names = FALSE) + + # ------- the decision table: per MATRIX, not per (matrix, seed) ------- + # Pairing on cells would be pseudo-replicated -- seeds within a matrix are not + # independent evidence about the ENGINE. + # ORIENTATION IS LOAD-BEARING. The ratio must be OLD / NEW, so that "> 1" means + # "the March engine built a worse (longer) tree" -- which is what every sentence + # below asserts. Alphabetical order would put `cur` in the numerator and silently + # invert every count against its own caption. + engines <- unique(scores$engine) + engines <- c(intersect(c("t252"), engines), sort(setdiff(engines, "t252"))) + if (length(engines) == 2L) { + medOf <- function(k, e) { + v <- scores$score[scores$key == k & scores$engine == e] + if (!length(v) || all(is.na(v))) NA_real_ else median(v, na.rm = TRUE) + } + keys <- unique(scores$key) + tab <- data.frame( + key = keys, + nTip = manifest$nTip[match(keys, manifest$key)], + a = vapply(keys, medOf, double(1), e = engines[[1]]), + b = vapply(keys, medOf, double(1), e = engines[[2]]), + stringsAsFactors = FALSE + ) + names(tab)[3:4] <- engines + tab$ratio <- tab[[3]] / tab[[4]] + tab <- tab[order(-tab$ratio), ] + write.csv(tab, file.path(outDir, "decision_table.csv"), row.names = FALSE) + cat("\n==== per-MATRIX medians (", engines[[1]], " vs ", engines[[2]], + ") ====\n", sep = "") + print(tab, row.names = FALSE) + ok <- !is.na(tab$ratio) + cat("\nmatrices where", engines[[1]], "is WORSE (ratio > 1):", + sum(tab$ratio[ok] > 1), "of", sum(ok), "\n") + cat("matrices within 1%:", sum(abs(tab$ratio[ok] - 1) < 0.01), "\n") + cat("matrices >10% worse:", sum(tab$ratio[ok] > 1.10), "\n") + cat("median ratio:", median(tab$ratio[ok]), "\n") + cat("\nREAD THIS AS: broad ratios > 1 => the t253 n=23 analysis rests on", + "systematically bad start trees (RETRACT). A ratio > 1 on project4284", + "alone => ANNOTATE that row.\n") + } else { + cat("\nOnly", length(engines), "engine(s) present; run both arms.\n") + } +} diff --git a/dev/benchmarks/t253_wagner_era_probe.sh b/dev/benchmarks/t253_wagner_era_probe.sh new file mode 100644 index 000000000..1f5b8a6fd --- /dev/null +++ b/dev/benchmarks/t253_wagner_era_probe.sh @@ -0,0 +1,55 @@ +#!/bin/bash +#SBATCH --job-name=t253probe +#SBATCH -p shared +#SBATCH -n 1 +#SBATCH --mem=24G +#SBATCH --time=8:00:00 +#SBATCH -o /nobackup/pjjg18/reach/logs/t253probe_%j.out +#SBATCH -e /nobackup/pjjg18/reach/logs/t253probe_%j.err +# +# Decides annotate-vs-retract on the t253 gap analysis. See the .R for the design. +# +# NOT an array, deliberately: the three steps are SEQUENTIAL (prep feeds both +# addition arms, and both arms feed score), and the whole job is ~150 bare +# `AdditionTree` calls -- an array would need a barrier for no gain. +# +# The archived March engine is the point of this job. It needs r/4.5.1 (4.4.1 +# refuses it: "built under R version 4.5.1") and its own dep path for Rcpp/TreeTools. +set -u +module load r/4.5.1 +export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 + +OUT_DIR=/nobackup/pjjg18/reach/t253probe +PROBE=/nobackup/pjjg18/reach/t253_wagner_era_probe.R + +CURLIB="/nobackup/pjjg18/curlib:/nobackup/pjjg18/TreeSearch/lib" +# lib-t252 = TreeSearch 2.0.0 packaged 2026-03-27 10:12, i.e. 32 min before +# t252_mbank_30s_20260327_1044.csv was written. lib-baseline supplies TreeTools 2.2.0 +# (and Rcpp, which curlib also needs from TreeSearch/lib). +T252LIB="/nobackup/pjjg18/TreeSearch/lib-t252:/nobackup/pjjg18/ts-bench/lib-baseline" + +mkdir -p "$OUT_DIR" /nobackup/pjjg18/reach/logs + +echo "=== STEP 1/4: prep (one preprocessing, shared by both arms) ===" +R_LIBS="$CURLIB" \ +STEP=prep OUT_DIR="$OUT_DIR" \ +NEOTRANS_DIR=/nobackup/pjjg18/neotrans/inst/matrices \ +CAT_CSV=/nobackup/pjjg18/reach/mbank_catalogue.csv \ + Rscript "$PROBE" || { echo "PREP FAILED"; exit 1; } + +echo "=== STEP 2/4: AdditionTree under the MARCH engine (lib-t252) ===" +R_LIBS="$T252LIB" \ +STEP=addition ENGINE=t252 OUT_DIR="$OUT_DIR" N_SEEDS=3 \ + Rscript "$PROBE" + +echo "=== STEP 3/4: AdditionTree under the CURRENT engine (curlib) ===" +R_LIBS="$CURLIB" \ +STEP=addition ENGINE=cur OUT_DIR="$OUT_DIR" N_SEEDS=3 \ + Rscript "$PROBE" + +echo "=== STEP 4/4: score EVERY tree with ONE scorer (curlib) ===" +R_LIBS="$CURLIB" \ +STEP=score OUT_DIR="$OUT_DIR" \ + Rscript "$PROBE" + +echo "=== done; results in $OUT_DIR ===" From f5f6dae0db7c0a473dc9eac1ac6216df864498aa Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Tue, 4 Aug 2026 07:13:26 +0100 Subject: [PATCH 14/16] bench(t253): refuse to report unless the two arms used different libraries Both engines report `Version: 2.0.0`, so the recorded version string cannot tell them apart. If `R_LIBS` failed to take, both arms would run the SAME library, every ratio would be 1.000, and that reads exactly like the "annotate, not retract" answer -- a null result obtainable from a broken arm. The addition step now records the resolved library path, and the score step refuses to produce a decision table unless the two arms provably came from different paths. Verified by running it against the deliberate same-library control, where it stops with PROVENANCE FAILURE instead of reporting ratio 1.000 on all three matrices. Submitted as Hamilton job 18183132. Co-Authored-By: Claude Opus 5 --- dev/benchmarks/t253_wagner_era_probe.R | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/dev/benchmarks/t253_wagner_era_probe.R b/dev/benchmarks/t253_wagner_era_probe.R index bb7978b77..c550ab30d 100644 --- a/dev/benchmarks/t253_wagner_era_probe.R +++ b/dev/benchmarks/t253_wagner_era_probe.R @@ -136,7 +136,7 @@ if (identical(step, "addition")) { if (inherits(tr, "error")) { cat(sprintf("FAIL %-18s seed=%d: %s\n", key, sd, conditionMessage(tr))) rows[[length(rows) + 1L]] <- data.frame( - key = key, engine = engine, tsVersion = tsVer, seed = sd, + key = key, engine = engine, tsVersion = tsVer, tsLib = tsLib, seed = sd, wallS = wall, ok = FALSE, stringsAsFactors = FALSE ) next @@ -146,7 +146,7 @@ if (identical(step, "addition")) { saveRDS(tr, file.path(treeDir, sprintf("%s__%s__s%d.rds", safeKey(key), engine, sd))) rows[[length(rows) + 1L]] <- data.frame( - key = key, engine = engine, tsVersion = tsVer, seed = sd, + key = key, engine = engine, tsVersion = tsVer, tsLib = tsLib, seed = sd, wallS = wall, ok = TRUE, stringsAsFactors = FALSE ) cat(sprintf("add %-18s %-6s seed=%d %.2fs\n", key, engine, sd, wall)) @@ -164,6 +164,27 @@ if (identical(step, "score")) { cat("scorer: TreeSearch ", as.character(utils::packageVersion("TreeSearch")), "\n", sep = "") manifest <- read.csv(file.path(outDir, "manifest.csv"), stringsAsFactors = FALSE) + + # ---- PROVENANCE ASSERTION: both engines report Version 2.0.0, so the version + # string CANNOT distinguish them. If R_LIBS failed to take, both arms would run + # the SAME library and every ratio would be 1.000 -- which reads exactly like the + # "annotate, not retract" answer. A null result must not be obtainable from a + # broken arm, so refuse to report unless the two arms provably used different + # library paths. + addFiles <- list.files(outDir, pattern = "^addition_.*\\.csv$", full.names = TRUE) + if (length(addFiles) >= 2L) { + prov <- do.call(rbind, lapply(addFiles, read.csv, stringsAsFactors = FALSE)) + byEngine <- unique(prov[, c("engine", "tsVersion", "tsLib")]) + cat("\n---- provenance ----\n") + print(byEngine, row.names = FALSE) + if (anyDuplicated(byEngine$tsLib)) { + stop("PROVENANCE FAILURE: two engine labels resolved to the SAME library ", + "path, so the arms are not independent. R_LIBS did not take. Any ratio ", + "of 1.000 here would be an artefact, not a finding.") + } + cat("provenance OK: arms used distinct libraries\n") + } + files <- list.files(treeDir, pattern = "\\.rds$", full.names = TRUE) rows <- list() for (f in files) { From b4e9c37f7b7eb857881c81b780e58ee4ea2d37c4 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Tue, 4 Aug 2026 07:18:36 +0100 Subject: [PATCH 15/16] =?UTF-8?q?bench(t253):=20RETRACT=20the=20MorphoBank?= =?UTF-8?q?=20gap=20analysis=20=E2=80=94=2025=20of=2025=20matrices=20affec?= =?UTF-8?q?ted?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Hamilton job 18183151 answers the annotate-vs-retract question, and it answers retract. Bare `AdditionTree`, no search, identical preprocessing, 3 seeds, every tree re-scored by one scorer, provenance asserted: matrices where the March engine built a LONGER tree 25 of 25 matrices >10% longer 23 of 25 median ratio (March / current) 1.365 range 1.050 -> 3.232 So project4284 (3.23x) is the EXTREME of a universal effect, not a special case to be excused in a footnote. A per-matrix "did it converge in 30 s" proxy built on start trees 5-223% too long cannot support the rho values t253 quotes, so the MorphoBank half of `t253_gap_characterization.md` is retracted and `t253_conv_gap_mbank.csv` gets a sidecar (a header comment would break every existing reader). Scoped deliberately: the t265 half — 8 named datasets, TNT vs TreeSearch at 120 s — does not touch the t252 engine and stands. The CONCLUSION that ntax predicts difficulty is corroborated elsewhere and is not called false; only this document's evidence for it is withdrawn. Two probe defects found by running it, both of which would have shipped a partial answer that looked complete: - the `__` filename separator collided with `safeKey()` output for parenthesised keys, silently dropping 3 of 25 matrices from the first run. The separator is now a character safeKey cannot emit, and the score step reports coverage and names any missing matrix. - both libraries report `Version: 2.0.0`, so a failed `R_LIBS` would have run both arms on one engine, produced ratio 1.000 everywhere, and read exactly like "annotate, not retract". The score step now refuses to report unless the arms provably resolved to different library paths. Result CSVs are force-added past `.gitignore`'s blanket `dev/benchmarks/*.csv`, matching the 73 result CSVs already tracked there — this one is the retraction's evidence. Co-Authored-By: Claude Opus 5 --- .../t253_conv_gap_mbank.RETRACTED.md | 29 ++++ dev/benchmarks/t253_gap_characterization.md | 32 ++++ dev/benchmarks/t253_wagner_era_decision.csv | 26 +++ dev/benchmarks/t253_wagner_era_probe.R | 28 +++- dev/benchmarks/t253_wagner_era_scores.csv | 151 ++++++++++++++++++ 5 files changed, 263 insertions(+), 3 deletions(-) create mode 100644 dev/benchmarks/t253_conv_gap_mbank.RETRACTED.md create mode 100644 dev/benchmarks/t253_wagner_era_decision.csv create mode 100644 dev/benchmarks/t253_wagner_era_scores.csv diff --git a/dev/benchmarks/t253_conv_gap_mbank.RETRACTED.md b/dev/benchmarks/t253_conv_gap_mbank.RETRACTED.md new file mode 100644 index 000000000..1fb763f91 --- /dev/null +++ b/dev/benchmarks/t253_conv_gap_mbank.RETRACTED.md @@ -0,0 +1,29 @@ +# ⛔ `t253_conv_gap_mbank.csv` IS RETRACTED (2026-08-03) + +A sidecar rather than a header, because a comment line inside the CSV would break every +reader that currently parses it. + +**Do not use `t253_conv_gap_mbank.csv` as evidence.** Its rows are derived from the +`t252_mbank_*` CSVs (2026-03-27), produced by an engine whose Wagner addition was worse on +**every matrix in the sample** — so the "convergence gap" it tabulates is substantially the +addition bug rather than a property of the datasets. + +Measured, Hamilton job 18183151 (`t253_wagner_era_probe.R`, results in +`t253_wagner_era_decision.csv`): bare `AdditionTree`, no search, identical preprocessing for +both engines, 3 seeds each, every tree re-scored by a single scorer, and a provenance +assertion confirming the two arms really used different libraries. + +| | | +|---|---| +| matrices where the March engine built a longer tree | **25 of 25** | +| matrices >10% longer | **23 of 25** | +| median ratio (March ÷ current) | **1.365** | +| range | 1.050 → 3.232 | + +The probe was built to decide *annotate one row* vs *retract the analysis*. 25 of 25 answers +retract: `project4284` (3.23×) is the extreme of a universal effect, not an outlier to be +excused. + +Full reasoning, and what survives, in the retraction box at the top of +`t253_gap_characterization.md`. Note that its `t265` half — 8 named datasets, TNT vs +TreeSearch at 120 s — does **not** depend on the t252 engine and stands. diff --git a/dev/benchmarks/t253_gap_characterization.md b/dev/benchmarks/t253_gap_characterization.md index c1ed84fbb..2d4259c27 100644 --- a/dev/benchmarks/t253_gap_characterization.md +++ b/dev/benchmarks/t253_gap_characterization.md @@ -1,5 +1,37 @@ # T-253: Gap Characterization by Dataset Features +> # ⛔ RETRACTED 2026-08-03 — DO NOT CITE THE MorphoBank HALF OF THIS DOCUMENT +> +> The `t252_mbank_*` half of this analysis is built on an engine whose Wagner addition was +> **worse on every matrix in the sample**, so its "convergence proxy" measures the addition +> bug, not dataset difficulty. +> +> **Measured** (Hamilton job 18183151, `t253_wagner_era_probe.R`; bare `AdditionTree`, no +> search, identical preprocessing, 3 seeds, every tree re-scored by ONE scorer; +> `t253_wagner_era_decision.csv`): +> +> | | | +> |---|---| +> | matrices where the March engine built a **longer** tree | **25 of 25** | +> | matrices >10% longer | **23 of 25** | +> | median ratio (March ÷ current) | **1.365** | +> | range | 1.050 (project561) → **3.232** (project4284) | +> +> This was run to decide *annotate one row* vs *retract*, and it answers **retract**: +> `project4284` is the extreme of a **universal** effect, not a special case. A per-matrix +> "did it converge in 30 s" proxy derived from start trees that are 5–223% too long cannot +> support the ρ values quoted below. +> +> **What does NOT fall with it:** the `t265` half (8 named datasets, TNT vs TreeSearch at +> 120 s) is independent of the t252 engine and is untouched. The *conclusion* that ntax +> predicts difficulty is also corroborated elsewhere and is not being called false — only +> **this document's evidence for it** is withdrawn. `t253_conv_gap_mbank.csv` carries the +> same defect; `headtohead_phase0.csv` was checked and contains no project4284 row, so no +> canonical target was contaminated. +> +> Cf. `reach_escalation_FINDINGS.md`; the archived-library technique that made this +> measurable at all is the reusable part. + **Date:** 2026-03-27 **Agent:** F **Data sources:** diff --git a/dev/benchmarks/t253_wagner_era_decision.csv b/dev/benchmarks/t253_wagner_era_decision.csv new file mode 100644 index 000000000..70fe3f000 --- /dev/null +++ b/dev/benchmarks/t253_wagner_era_decision.csv @@ -0,0 +1,26 @@ +"key","nTip","t252","cur","ratio" +"project4284",4062,1267,392,3.23214285714286 +"project804",173,2710,1315,2.06083650190114 +"syab07201",125,27551,15168,1.81638976793249 +"project3938",119,6111,3492,1.75 +"project4133",131,3855,2425,1.58969072164948 +"project826",33,705,452,1.55973451327434 +"project3688",60,1356,885,1.53220338983051 +"project4286",63,414,280,1.47857142857143 +"project4146_(3)",59,371,260,1.42692307692308 +"project4049",60,7505,5362,1.39966430436404 +"project4359",71,265,190,1.39473684210526 +"project4397",75,2333,1708,1.36592505854801 +"project2346",23,411,301,1.36544850498339 +"project423",60,672,497,1.35211267605634 +"project2084_(1)",86,33315,25150,1.32465208747515 +"project2771",94,1270,964,1.31742738589212 +"project2451",24,973,741,1.31309041835358 +"project2184",114,755,595,1.26890756302521 +"project944",25,168,135,1.24444444444444 +"project2762",29,306,261,1.17241379310345 +"project571",42,751,659,1.13960546282246 +"project971_(1)",26,180,159,1.13207547169811 +"project532",21,1282,1144,1.12062937062937 +"project4501",24,129,120,1.075 +"project561",34,1215,1157,1.05012964563526 diff --git a/dev/benchmarks/t253_wagner_era_probe.R b/dev/benchmarks/t253_wagner_era_probe.R index c550ab30d..e759bb34c 100644 --- a/dev/benchmarks/t253_wagner_era_probe.R +++ b/dev/benchmarks/t253_wagner_era_probe.R @@ -64,6 +64,13 @@ seedsFor <- function(i) BASE_SEED + seq_len(nSeeds) - 1L + (i - 1L) * 100L safeKey <- function(k) gsub("[^A-Za-z0-9]", "_", k) +# Filename field separator. It must be a character `safeKey()` can NEVER emit, or +# keys like "project2084_(1)" -- which safeKey renders "project2084___1_" -- split +# into the wrong number of fields and get silently dropped from the results. That +# cost 3 of 25 matrices on the first run: a partial answer that still LOOKED +# complete. "_" is unusable for exactly that reason; "@" is not in [A-Za-z0-9]. +FS <- "@@" + # ---------------------------------------------------------------- STEP: prep ---- if (identical(step, "prep")) { suppressPackageStartupMessages(library("TreeTools")) @@ -143,8 +150,8 @@ if (identical(step, "addition")) { } # Write the deliverable BEFORE anything that could fail on it: a verification # error must not take the computed tree with it (the 2026-07-31 lesson). - saveRDS(tr, file.path(treeDir, sprintf("%s__%s__s%d.rds", - safeKey(key), engine, sd))) + saveRDS(tr, file.path(treeDir, paste0( + paste(safeKey(key), engine, paste0("s", sd), sep = FS), ".rds"))) rows[[length(rows) + 1L]] <- data.frame( key = key, engine = engine, tsVersion = tsVer, tsLib = tsLib, seed = sd, wallS = wall, ok = TRUE, stringsAsFactors = FALSE @@ -187,16 +194,19 @@ if (identical(step, "score")) { files <- list.files(treeDir, pattern = "\\.rds$", full.names = TRUE) rows <- list() + nSkipped <- 0L for (f in files) { - parts <- strsplit(sub("\\.rds$", "", basename(f)), "__", fixed = TRUE)[[1]] + parts <- strsplit(sub("\\.rds$", "", basename(f)), FS, fixed = TRUE)[[1]] if (length(parts) != 3L) { cat("SKIP unparseable filename:", basename(f), "\n") + nSkipped <- nSkipped + 1L next } keySafe <- parts[[1]] idx <- match(keySafe, safeKey(manifest$key)) if (is.na(idx)) { cat("SKIP no manifest row:", basename(f), "\n") + nSkipped <- nSkipped + 1L next } key <- manifest$key[idx] @@ -216,6 +226,18 @@ if (identical(step, "score")) { scores <- do.call(rbind, rows) write.csv(scores, file.path(outDir, "scores.csv"), row.names = FALSE) + # ---- COMPLETENESS: say out loud what is missing. A table covering 22 of 25 + # matrices reads exactly like a table covering all of them, and the first run of + # this probe did precisely that (3 parenthesised keys lost to the field separator). + covered <- unique(scores$key) + missing <- setdiff(manifest$key, covered) + cat("\ncoverage:", length(covered), "of", nrow(manifest), "matrices;", + nSkipped, "tree files skipped\n") + if (length(missing)) { + cat("!! MISSING MATRICES (the result below is PARTIAL):", + paste(missing, collapse = ", "), "\n") + } + # ------- the decision table: per MATRIX, not per (matrix, seed) ------- # Pairing on cells would be pseudo-replicated -- seeds within a matrix are not # independent evidence about the ENGINE. diff --git a/dev/benchmarks/t253_wagner_era_scores.csv b/dev/benchmarks/t253_wagner_era_scores.csv new file mode 100644 index 000000000..83a42bb07 --- /dev/null +++ b/dev/benchmarks/t253_wagner_era_scores.csv @@ -0,0 +1,151 @@ +"key","engine","seed","nTip","nChar","score" +"project2084_(1)","cur",3001,86,3570,25193 +"project2084_(1)","cur",3002,86,3570,25150 +"project2084_(1)","cur",3003,86,3570,25082 +"project2084_(1)","t252",3001,86,3570,33265 +"project2084_(1)","t252",3002,86,3570,33315 +"project2084_(1)","t252",3003,86,3570,34335 +"project2184","cur",3201,114,168,595 +"project2184","cur",3202,114,168,613 +"project2184","cur",3203,114,168,586 +"project2184","t252",3201,114,168,755 +"project2184","t252",3202,114,168,772 +"project2184","t252",3203,114,168,705 +"project2346","cur",1401,23,141,299 +"project2346","cur",1402,23,141,301 +"project2346","cur",1403,23,141,301 +"project2346","t252",1401,23,141,421 +"project2346","t252",1402,23,141,411 +"project2346","t252",1403,23,141,393 +"project2451","cur",1501,24,367,737 +"project2451","cur",1502,24,367,741 +"project2451","cur",1503,24,367,744 +"project2451","t252",1501,24,367,1027 +"project2451","t252",1502,24,367,973 +"project2451","t252",1503,24,367,947 +"project2762","cur",1901,29,171,258 +"project2762","cur",1902,29,171,261 +"project2762","cur",1903,29,171,262 +"project2762","t252",1901,29,171,319 +"project2762","t252",1902,29,171,305 +"project2762","t252",1903,29,171,306 +"project2771","cur",3101,94,123,967 +"project2771","cur",3102,94,123,964 +"project2771","cur",3103,94,123,957 +"project2771","t252",3101,94,123,1279 +"project2771","t252",3102,94,123,1270 +"project2771","t252",3103,94,123,1183 +"project3688","cur",2401,60,245,870 +"project3688","cur",2402,60,245,894 +"project3688","cur",2403,60,245,885 +"project3688","t252",2401,60,245,1263 +"project3688","t252",2402,60,245,1356 +"project3688","t252",2403,60,245,1428 +"project3938","cur",3301,119,677,3466 +"project3938","cur",3302,119,677,3492 +"project3938","cur",3303,119,677,3493 +"project3938","t252",3301,119,677,6052 +"project3938","t252",3302,119,677,6212 +"project3938","t252",3303,119,677,6111 +"project4049","cur",2501,60,719,5301 +"project4049","cur",2502,60,719,5362 +"project4049","cur",2503,60,719,5398 +"project4049","t252",2501,60,719,7739 +"project4049","t252",2502,60,719,7505 +"project4049","t252",2503,60,719,7100 +"project4133","cur",3501,131,349,2432 +"project4133","cur",3502,131,349,2382 +"project4133","cur",3503,131,349,2425 +"project4133","t252",3501,131,349,3797 +"project4133","t252",3502,131,349,3971 +"project4133","t252",3503,131,349,3855 +"project4146_(3)","cur",2301,59,130,260 +"project4146_(3)","cur",2302,59,130,266 +"project4146_(3)","cur",2303,59,130,258 +"project4146_(3)","t252",2301,59,130,371 +"project4146_(3)","t252",2302,59,130,359 +"project4146_(3)","t252",2303,59,130,420 +"project423","cur",2601,60,212,497 +"project423","cur",2602,60,212,496 +"project423","cur",2603,60,212,510 +"project423","t252",2601,60,212,672 +"project423","t252",2602,60,212,682 +"project423","t252",2603,60,212,600 +"project4284","cur",3701,4062,27,388 +"project4284","cur",3702,4062,27,393 +"project4284","cur",3703,4062,27,392 +"project4284","t252",3701,4062,27,1267 +"project4284","t252",3702,4062,27,1264 +"project4284","t252",3703,4062,27,1309 +"project4286","cur",2701,63,135,280 +"project4286","cur",2702,63,135,278 +"project4286","cur",2703,63,135,280 +"project4286","t252",2701,63,135,418 +"project4286","t252",2702,63,135,376 +"project4286","t252",2703,63,135,414 +"project4359","cur",2801,71,114,193 +"project4359","cur",2802,71,114,190 +"project4359","cur",2803,71,114,187 +"project4359","t252",2801,71,114,251 +"project4359","t252",2802,71,114,265 +"project4359","t252",2803,71,114,266 +"project4397","cur",2901,75,222,1715 +"project4397","cur",2902,75,222,1679 +"project4397","cur",2903,75,222,1708 +"project4397","t252",2901,75,222,2306 +"project4397","t252",2902,75,222,2373 +"project4397","t252",2903,75,222,2333 +"project4501","cur",1601,24,41,121 +"project4501","cur",1602,24,41,120 +"project4501","cur",1603,24,41,120 +"project4501","t252",1601,24,41,132 +"project4501","t252",1602,24,41,129 +"project4501","t252",1603,24,41,126 +"project532","cur",1301,21,420,1144 +"project532","cur",1302,21,420,1140 +"project532","cur",1303,21,420,1145 +"project532","t252",1301,21,420,1420 +"project532","t252",1302,21,420,1282 +"project532","t252",1303,21,420,1213 +"project561","cur",2101,34,329,1157 +"project561","cur",2102,34,329,1159 +"project561","cur",2103,34,329,1151 +"project561","t252",2101,34,329,1215 +"project561","t252",2102,34,329,1242 +"project561","t252",2103,34,329,1194 +"project571","cur",2201,42,125,659 +"project571","cur",2202,42,125,652 +"project571","cur",2203,42,125,667 +"project571","t252",2201,42,125,814 +"project571","t252",2202,42,125,723 +"project571","t252",2203,42,125,751 +"project804","cur",3601,173,565,1313 +"project804","cur",3602,173,565,1316 +"project804","cur",3603,173,565,1315 +"project804","t252",3601,173,565,2710 +"project804","t252",3602,173,565,2795 +"project804","t252",3603,173,565,2676 +"project826","cur",2001,33,213,452 +"project826","cur",2002,33,213,492 +"project826","cur",2003,33,213,451 +"project826","t252",2001,33,213,799 +"project826","t252",2002,33,213,705 +"project826","t252",2003,33,213,672 +"project944","cur",1701,25,72,135 +"project944","cur",1702,25,72,137 +"project944","cur",1703,25,72,130 +"project944","t252",1701,25,72,166 +"project944","t252",1702,25,72,168 +"project944","t252",1703,25,72,202 +"project971_(1)","cur",1801,26,73,161 +"project971_(1)","cur",1802,26,73,159 +"project971_(1)","cur",1803,26,73,158 +"project971_(1)","t252",1801,26,73,180 +"project971_(1)","t252",1802,26,73,190 +"project971_(1)","t252",1803,26,73,166 +"syab07201","cur",3401,125,2813,15168 +"syab07201","cur",3402,125,2813,15333 +"syab07201","cur",3403,125,2813,15041 +"syab07201","t252",3401,125,2813,27499 +"syab07201","t252",3402,125,2813,27551 +"syab07201","t252",3403,125,2813,27660 From 7e105ee5522729be15aab84fa2b3905bcaf19d0c Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Tue, 4 Aug 2026 07:20:19 +0100 Subject: [PATCH 16/16] docs(reach): the March-engine gap is every matrix, not project4284's Closes the loop between this study's incidental finding and the t253 retraction: the 1590-vs-409 observation that started as an explanation for one anomalous row turns out to hold on 25 of 25 matrices (median 1.365). Co-Authored-By: Claude Opus 5 --- dev/benchmarks/reach_escalation_FINDINGS.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md index 0f04b76f5..95e0d7a5c 100644 --- a/dev/benchmarks/reach_escalation_FINDINGS.md +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -146,6 +146,16 @@ five times faster, the signature of the since-fixed union-of-finals insertion-co is not the explanation: the current engine returns the same score at `maxSeconds` 1, 5, 25 and 45, and the whole 30 s → 1440 s span is only ×1.068. +**And it is not project4284's peculiarity either — it is EVERY matrix** (job 18183151, +`t253_wagner_era_probe.R`, `t253_wagner_era_decision.csv`). Extending the same bare-`AdditionTree` +comparison to all 25 MBANK_FIXED_SAMPLE matrices, with one shared preprocessing, 3 seeds each and +every tree re-scored by a single scorer: the March engine built a **longer tree on 25 of 25**, +**23 of 25 by more than 10%**, median ratio **1.365**, range 1.050 → 3.232. project4284 (3.23×) +is the extreme of a universal effect. This is what retired the last hope that one row needed a +footnote: **`t253_gap_characterization.md`'s MorphoBank half is now retracted**, since a "did it +converge in 30 s" proxy cannot be built on start trees 5–223% too long. Its `t265` half is +independent of that engine and stands. + **Bad news — the per-seed scores did not reproduce.** | seed | 7731 | 7732 | 7733 | 7734 | 7735 | mean |