diff --git a/NEWS.md b/NEWS.md index 3492cf2d4..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,6 +182,32 @@ profile parsimony are unchanged, as is `thorough`/`large`, and setting any of these fields yourself overrides all of it. +- 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 no-improvement rule computes `(targetHits / hits) * nTip * perturbStopFactor`, diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 3ccce1f44..577b34b07 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,127 @@ 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). 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. 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 +# 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. +# +# 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 +# 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`. +# +# `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, userSetHits, + userSet = character(0)) { + if (!isTRUE(userSetHits)) { + return(control) + } + 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() { @@ -820,6 +949,28 @@ #' 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. +#' +#' 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 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. @@ -1117,6 +1268,8 @@ MaximizeParsimony <- function( } rungName <- .rung } + # Set when the reach escalation below raises `poolSuboptimal` itself; see there. + escalatedPool <- FALSE # --- Apply the rung --- if (!identical(rungName, "none")) { @@ -1183,9 +1336,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 @@ -1236,6 +1389,46 @@ MaximizeParsimony <- function( for (.f in names(iwStop)) { control[[.f]] <- iwStop[[.f]] } + + # 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 + # 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 the rung's preset, 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. + # `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 -- + # `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 && userSetHits && + escalation >= .reachEscalationMinRatio && + strategy %in% c("thorough", "large")) { + cli::cli_alert_info( + "Deep search: {.field targetHits} raised {round(escalation, 1)}x" + ) + } } } @@ -1672,7 +1865,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/dev/benchmarks/reach_escalation_FINDINGS.md b/dev/benchmarks/reach_escalation_FINDINGS.md new file mode 100644 index 000000000..95e0d7a5c --- /dev/null +++ b/dev/benchmarks/reach_escalation_FINDINGS.md @@ -0,0 +1,373 @@ +# `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.** 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): + +| 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". + +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 +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.** + +## 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 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. 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), 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. + +**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 | +|------|------|------|------|------|------|------| +| 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 +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. + +## ✅ 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 +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. +**Answered NO, 2026-08-03 — see the resolution at the head of this section.** + +## 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`. +- **…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. +- **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 + 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 + +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.~~ **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` + 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 + 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 new file mode 100644 index 000000000..4074ccd5f --- /dev/null +++ b/dev/benchmarks/reach_escalation_ab.R @@ -0,0 +1,238 @@ +#!/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, +# 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) +} + +# 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) + 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, + enum_time_fraction = ENUM_TIME_FRACTION, 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_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 new file mode 100644 index 000000000..fba5b6aff --- /dev/null +++ b/dev/benchmarks/reach_escalation_analyze.R @@ -0,0 +1,206 @@ +#!/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))) + +# 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))) { + 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], + # 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) + } +} +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). +# 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; 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", + 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") + +# 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 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)) { + 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 (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)")) + +# 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..37cf5848b --- /dev/null +++ b/dev/benchmarks/reach_recover4284.R @@ -0,0 +1,132 @@ +#!/usr/bin/env Rscript +# 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 +# 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. +# +# 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)) { + 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 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, + collapse = FALSE), # <- binary trees, and exposes the pool verbatim + DELTAS6))) +wall <- as.double(proc.time()["elapsed"] - t0) +reported <- as.double(attr(res, "score")) +trees <- if (inherits(res, "phylo")) structure(list(res), class = "multiPhylo") else res +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. +# 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 %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)) { + 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 (any(ok) && !agree) + cat(" !! MISMATCH -- reported score NOT reproduced by TreeLength on any returned tree.\n") + +# 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 = 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_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, collapse_arg = FALSE, + stringsAsFactors = FALSE), + file.path(out_dir, sprintf("recover_%02d_s%d.csv", tid, seed)), row.names = FALSE) +cat(" CSV written.\n") 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 new file mode 100644 index 000000000..e759bb34c --- /dev/null +++ b/dev/benchmarks/t253_wagner_era_probe.R @@ -0,0 +1,282 @@ +#!/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) + +# 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")) + 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, tsLib = tsLib, 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, 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 + ) + 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) + + # ---- 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() + nSkipped <- 0L + for (f in files) { + 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] + 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) + + # ---- 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. + # 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 ===" 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 diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 917751868..316726a4d 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -265,7 +265,29 @@ 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. + +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 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 new file mode 100644 index 000000000..1f6026272 --- /dev/null +++ b/tests/testthat/test-reach-escalation.R @@ -0,0 +1,253 @@ +# 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 + +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, + userSetHits = TRUE) + 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, + userSetHits = TRUE) + 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, + userSetHits = TRUE) + 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, + 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, + 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 + # 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, + userSetHits = TRUE), + stock, info = paste("strategy", strat) + ) + } +}) + +test_that(".ApplyReachEscalation preserves caller-set fields", { + ctrl <- TreeSearch:::.ApplyReachEscalation( + SearchControl(), "thorough", escalation = 4, userSetHits = TRUE, + 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. + # 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, effort = 2L, 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, effort = 0L, 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, effort = 2L, 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("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, effort = 2L, 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)) +})