T-207/T-210: Multi-cycle PCSA perturbation phase

R/MaximizeParsimony.R

@@ -172,7 +172,7 @@
     ratchetAdaptive = TRUE,
     nniPerturbCycles = 0L,
     driftCycles = 0L,
-    annealPhases = 5L, annealTStart = 20, annealTEnd = 0,
+    annealCycles = 3L, annealPhases = 5L, annealTStart = 20, annealTEnd = 0,
     xssRounds = 3L, xssPartitions = 6L,
     rssRounds = 2L, cssRounds = 1L, cssPartitions = 6L,
     sectorMinSize = 8L, sectorMaxSize = 100L,

R/SearchControl.R

@@ -103,14 +103,18 @@
 #'   an outer cycle, the counter resets up to this many times, allowing
 #'   productive re-exploration.  Set to \eqn{-1} for unlimited resets.
 #'   Strategy presets (`"default"`, `"thorough"`) set 2–3.
-#' @param annealPhases Integer; number of simulated annealing temperature
-#'   steps (default 0 = disabled).  When > 0, runs a linear cooling
-#'   schedule from `annealTStart` to `annealTEnd` using stochastic TBR
-#'   with Boltzmann acceptance.  Runs between drift and final TBR polish.
-#' @param annealTStart Numeric; initial Boltzmann temperature for annealing
-#'   (default 20).  Higher temperatures accept more suboptimal moves.
+#' @param annealCycles Integer; number of simulated annealing perturbation
+#'   cycles (PCSA) per replicate.  Each cycle perturbs the current best tree
+#'   via scheduled SA cooling, then reconverges with TBR.  If the result
+#'   improves on the best, it becomes the new starting point.  Effective at
+#'   escaping deep basins under equal-weights parsimony at \eqn{\ge}100 tips.
+#'   0 (default) disables SA perturbation.
+#' @param annealPhases Integer; number of temperature steps in the linear
+#'   cooling schedule per SA cycle (default 5).
+#' @param annealTStart Numeric; initial Boltzmann temperature for SA cooling
+#'   schedule (default 20).  Higher temperatures accept more suboptimal moves.
 #' @param annealTEnd Numeric; final Boltzmann temperature (default 0 =
-#'   strict hill-climbing at end).
+#'   strict hill-climbing at end of each cycle).
 #' @param annealMovesPerPhase Integer; stochastic TBR moves per temperature
 #'   step (default 0 = number of tips).
 #' @param enumTimeFraction Numeric between 0 and 0.5; fraction of `maxSeconds`
@@ -190,8 +194,9 @@ SearchControl <- function(
     consensusStableReps = 0L,
     adaptiveLevel = FALSE,
     consensusConstrain = FALSE,
-    # Simulated annealing (linear cooling schedule)
-    annealPhases = 0L,
+    # Simulated annealing perturbation (PCSA, T-207)
+    annealCycles = 0L,
+    annealPhases = 5L,
     annealTStart = 20,
     annealTEnd = 0,
     annealMovesPerPhase = 0L,
@@ -238,6 +243,7 @@ SearchControl <- function(
       consensusStableReps = as.integer(consensusStableReps),
       adaptiveLevel = as.logical(adaptiveLevel),
       consensusConstrain = as.logical(consensusConstrain),
+      annealCycles = as.integer(annealCycles),
       annealPhases = as.integer(annealPhases),
       annealTStart = as.double(annealTStart),
       annealTEnd = as.double(annealTEnd),
@@ -259,8 +265,8 @@ print.SearchControl <- function(x, ...) {
                    "ratchetTaper"),
     "NNI Perturbation" = c("nniPerturbCycles", "nniPerturbFraction"),
     "Drift" = c("driftCycles", "driftAfdLimit", "driftRfdLimit"),
-    "Annealing" = c("annealPhases", "annealTStart", "annealTEnd",
-                     "annealMovesPerPhase"),
+    "Annealing" = c("annealCycles", "annealPhases", "annealTStart",
+                     "annealTEnd", "annealMovesPerPhase"),
     "Sectorial" = c("xssRounds", "xssPartitions", "rssRounds",
                      "cssRounds", "cssPartitions",
                      "sectorMinSize", "sectorMaxSize"),

R/ts-driven-compat.R

@@ -131,7 +131,11 @@ ts_driven_search <- function(
 
   # Anneal config: fold into SearchControl if provided
   if (!is.null(annealConfig)) {
-    sc$annealPhases <- as.integer(annealConfig$phases %||% 0L)
+    phases <- as.integer(annealConfig$phases %||% 5L)
+    # Backward compat: if phases > 0 but cycles not specified, default to 1
+    sc$annealCycles <- as.integer(annealConfig$cycles %||%
+                                    if (phases > 0L) 1L else 0L)
+    sc$annealPhases <- phases
     sc$annealTStart <- as.double(annealConfig$tStart %||% 20)
     sc$annealTEnd <- as.double(annealConfig$tEnd %||% 0)
     sc$annealMovesPerPhase <- as.integer(annealConfig$movesPerPhase %||% 0L)

src/ts_driven.cpp

@@ -366,19 +366,45 @@ ReplicateResult run_single_replicate(
       return result;
     }
 
-    // 5b. Simulated annealing (linear cooling schedule)
-    if (params.anneal_phases > 0) {
+    // 5b. SA perturbation (multi-cycle PCSA with best-tree restart).
+    // Each cycle: perturb best tree via scheduled SA cooling -> TBR
+    // reconverge -> keep if improved (T-207).
+    if (params.anneal_cycles > 0) {
       AnnealParams ap;
       ap.t_start = params.anneal_t_start;
       ap.t_end = params.anneal_t_end;
       ap.n_phases = params.anneal_phases;
       ap.moves_per_phase = params.anneal_moves_per_phase;
-      anneal_search(result.tree, ds, ap, cd, check_timeout);
 
+      double best_sa_score = score_tree(result.tree, ds);
+      TreeState best_sa_tree = result.tree;
+
+      for (int cyc = 0; cyc < params.anneal_cycles; ++cyc) {
+        if (cyc > 0) result.tree = best_sa_tree;
+
+        anneal_search(result.tree, ds, ap, cd, check_timeout);
+
+        // TBR reconverge after SA perturbation
+        {
+          TBRParams tp;
+          tp.tabu_size = params.tabu_size;
+          tbr_search(result.tree, ds, tp, cd, nullptr, nullptr, check_timeout);
+        }
+
+        double cyc_score = score_tree(result.tree, ds);
+        if (cyc_score < best_sa_score - 1e-10) {
+          best_sa_score = cyc_score;
+          best_sa_tree = result.tree;
+        }
+
+        if (ts::check_interrupt() || check_timeout()) break;
+      }
+
+      result.tree = best_sa_tree;
       result.timings.anneal_ms += ph_lap();
       if (verbosity >= 2) {
         Rprintf("  %s score: %.5g [%.0f ms total]\n",
-                outer_label("Anneal").c_str(),
+                outer_label("SA").c_str(),
                 score_tree(result.tree, ds), result.timings.anneal_ms);
       }
     }

src/ts_driven.h

@@ -58,13 +58,16 @@ struct DrivenParams {
   int drift_afd_limit = 3;
   double drift_rfd_limit = 0.1;
 
-  // Simulated annealing: single-chain linear cooling schedule using
-  // stochastic TBR.  Runs between drift and final TBR polish.
-  // n_phases = 0 disables annealing.
-  int anneal_phases = 0;           // 0 = disabled
-  double anneal_t_start = 20.0;    // initial Boltzmann temperature
-  double anneal_t_end = 0.0;       // final temperature
-  int anneal_moves_per_phase = 0;  // 0 = n_tip
+  // Simulated annealing perturbation (PCSA: post-convergence SA).
+  // Multi-cycle SA with best-tree restart, inserted after drift phase.
+  // Each cycle: SA cooling schedule -> TBR reconverge -> keep if improved.
+  // Effective at escaping deep basins under EW at >=100 tips.
+  // 0 = disabled (default). Typical: 3-5 cycles for large trees.
+  int anneal_cycles = 0;
+  int anneal_phases = 5;             // temperature steps per cycle
+  double anneal_t_start = 20.0;      // initial Boltzmann temperature
+  double anneal_t_end = 0.0;         // final temperature
+  int anneal_moves_per_phase = 0;    // 0 = n_tip
 
   // Sectorial search
   int xss_rounds = 3;

src/ts_rcpp.cpp

@@ -1240,7 +1240,8 @@ static void unpack_search_control(List ctrl, ts::DrivenParams& params) {
   params.adaptive_start        = as<bool>(ctrl["adaptiveStart"]);
   params.enum_time_fraction    = as<double>(ctrl["enumTimeFraction"]);
 
-  // Simulated annealing (folded into searchControl)
+  // Simulated annealing perturbation (PCSA)
+  params.anneal_cycles          = as<int>(ctrl["annealCycles"]);
   params.anneal_phases          = as<int>(ctrl["annealPhases"]);
   params.anneal_t_start         = as<double>(ctrl["annealTStart"]);
   params.anneal_t_end           = as<double>(ctrl["annealTEnd"]);

src/ts_temper.cpp

@@ -385,9 +385,16 @@ AnnealResult anneal_search(
   double score = temper_full_rescore(tree, ds);
   result.best_score = score;
 
+  // Track best tree across phases (T-210). Boltzmann acceptance can
+  // displace a good tree found in an earlier phase; without tracking,
+  // PCSA reconverges from the final (possibly worse) SA state.
+  TreeState best_tree = tree;
+  double best_tree_score = score;
+
   for (int phase = 0; phase < np; ++phase) {
     // Linear temperature schedule: t_start -> t_end
-    double frac = (np == 1) ? 1.0
+    // Single phase: run at t_start (hot perturbation, not t_end cold)
+    double frac = (np == 1) ? 0.0
                   : static_cast<double>(phase) / (np - 1);
     double temp = params.t_start + (params.t_end - params.t_start) * frac;
 
@@ -404,10 +411,19 @@ AnnealResult anneal_search(
       result.best_score = tr.best_score;
     }
 
+    // Save tree at phase boundary if it improved
+    double phase_score = temper_full_rescore(tree, ds);
+    if (phase_score < best_tree_score - 1e-10) {
+      best_tree_score = phase_score;
+      best_tree = tree;
+    }
+
     if (ts::check_interrupt()) break;
     if (check_timeout && check_timeout()) break;
   }
 
+  // Restore best tree found across all phases
+  tree = best_tree;
   result.final_score = temper_full_rescore(tree, ds);
   return result;
 }

tests/testthat/test-ts-anneal.R

@@ -4,8 +4,10 @@ skip_on_cran()
 # --- Helpers ---
 # Build annealConfig list from SearchControl fields
 make_anneal_config <- function(ctrl) {
-  if (isTRUE(as.integer(ctrl$annealPhases) > 0L)) {
+  cycles <- as.integer(ctrl$annealCycles %||% 0L)
+  if (cycles > 0L) {
     list(
+      cycles = cycles,
       phases = as.integer(ctrl$annealPhases),
       tStart = as.double(ctrl$annealTStart),
       tEnd = as.double(ctrl$annealTEnd),
@@ -20,7 +22,7 @@ anneal_search <- function(ds, ctrl_overrides = list(), maxSeconds = 5,
                           maxReplicates = 3L, verbosity = 0L,
                           concavity = -1.0) {
   ctrl <- SearchControl(
-    annealPhases = 5L, annealTStart = 20, annealTEnd = 0,
+    annealCycles = 1L, annealPhases = 5L, annealTStart = 20, annealTEnd = 0,
     ratchetCycles = 2L, driftCycles = 0L,
     xssRounds = 0L, rssRounds = 0L, cssRounds = 0L,
     fuseInterval = 0L
@@ -67,7 +69,7 @@ test_that("anneal_ms > 0 when annealing is enabled", {
 
 test_that("anneal_ms = 0 when annealing is disabled", {
   result <- anneal_search(ds,
-    ctrl_overrides = list(annealPhases = 0L),
+    ctrl_overrides = list(annealCycles = 0L),
     maxSeconds = 3, maxReplicates = 2L
   )
   expect_equal(result$timings[["anneal_ms"]], 0)
@@ -81,8 +83,10 @@ test_that("Annealing with IW scoring works", {
 })
 
 test_that("SearchControl accepts annealing parameters", {
-  ctrl <- SearchControl(annealPhases = 5L, annealTStart = 15,
-                        annealTEnd = 1, annealMovesPerPhase = 50L)
+  ctrl <- SearchControl(annealCycles = 3L, annealPhases = 5L,
+                        annealTStart = 15, annealTEnd = 1,
+                        annealMovesPerPhase = 50L)
+  expect_equal(ctrl$annealCycles, 3L)
   expect_equal(ctrl$annealPhases, 5L)
   expect_equal(ctrl$annealTStart, 15)
   expect_equal(ctrl$annealTEnd, 1)
@@ -91,13 +95,15 @@ test_that("SearchControl accepts annealing parameters", {
 
 test_that("SearchControl defaults disable annealing", {
   ctrl <- SearchControl()
-  expect_equal(ctrl$annealPhases, 0L)
+  expect_equal(ctrl$annealCycles, 0L)
+  expect_equal(ctrl$annealPhases, 5L)
 })
 
 test_that("Large preset enables annealing and disables drift", {
   presets <- TreeSearch:::.AutoStrategy(200L, 200L)
   expect_equal(presets, "large")
   large_ctrl <- TreeSearch:::.StrategyPresets()[["large"]]
+  expect_equal(large_ctrl$annealCycles, 3L)
   expect_gt(large_ctrl$annealPhases, 0L)
   expect_equal(large_ctrl$driftCycles, 0L)
 })
@@ -113,12 +119,37 @@ test_that("Annealing with T=0 acts as strict hill-climbing", {
 })
 
 test_that("MaximizeParsimony respects annealing in SearchControl", {
-  ctrl <- SearchControl(annealPhases = 3L, annealTStart = 10,
-                        annealTEnd = 0, ratchetCycles = 1L,
+  ctrl <- SearchControl(annealCycles = 1L, annealPhases = 3L,
+                        annealTStart = 10, annealTEnd = 0,
+                        ratchetCycles = 1L,
                         driftCycles = 0L, xssRounds = 0L,
                         rssRounds = 0L, cssRounds = 0L)
   result <- MaximizeParsimony(dataset, concavity = Inf,
                               maxSeconds = 3, maxReplicates = 2L,
                               control = ctrl)
   expect_s3_class(result[[1]], "phylo")
 })
+
+test_that("Multi-cycle PCSA runs and reports sa_ms", {
+  result <- anneal_search(ds,
+    ctrl_overrides = list(annealCycles = 3L),
+    maxSeconds = 5, maxReplicates = 2L
+  )
+  expect_gt(result$pool_size, 0)
+  expect_gt(result$timings[["anneal_ms"]], 0)
+  validate_result(result, length(dataset))
+})
+
+test_that("Multi-cycle PCSA score <= single-cycle", {
+  set.seed(7418)
+  single <- anneal_search(ds,
+    ctrl_overrides = list(annealCycles = 1L),
+    maxSeconds = 5, maxReplicates = 3L
+  )
+  multi <- anneal_search(ds,
+    ctrl_overrides = list(annealCycles = 3L),
+    maxSeconds = 5, maxReplicates = 3L
+  )
+  # Multi-cycle should find scores at least as good (may tie)
+  expect_lte(multi$best_score, single$best_score + 1)
+})