Sync cpp-search from ms609 (reverse convergence, post-split)

NEWS.md

@@ -1,5 +1,25 @@
 # To integrate into 2.0.0 notes
 
+- New `SearchControl()` parameter `enumMaxTrees`: a retention ceiling applied to
+  the post-search MPT-enumeration phase alone.  `0` (the default)
+  keeps `poolMaxSize` throughout, so behaviour is unchanged unless you set it.
+  From `effort` rung 5 the ladder now doubles it each notch alongside the
+  replicate budget and the hit target, so asking for more effort also asks for a
+  more complete tree set — previously a search could be given eight times the
+  budget and still return only the default 100 trees.
+
+  `poolMaxSize` is deliberately **not** scaled, and the split is the point.
+  During the replicate loop the pool cap is not a ceiling on what is returned but
+  the size of the working set the search reads: fusing draws its donors from the
+  whole pool, conflict-guided sectorial search reads the pool's split
+  frequencies once per replicate, and `consensusConstrain` reads its consensus
+  splits.  Raising it therefore changes which trees the search *visits*, so the
+  anytime-dominance argument that licenses raising `maxReplicates` — a higher cap
+  only appends later replicates and can never delay an earlier improvement — does
+  not transfer to it.  Once the loop is over the pool is pure output, and there
+  the same argument does hold, which is why the ceiling is raised at that point
+  instead.  Raising `poolMaxSize` yourself still works and still governs both
+  phases; `enumMaxTrees` is the side-effect-free way to keep more trees.
 - Profile parsimony computes exactly for more multi-state characters, where it
   previously approximated nearly all of them.  The exact Maddison & Slatkin
   solver caches into fixed-capacity memo tables and bails out when one fills --

R/Concordance.R

@@ -177,20 +177,10 @@ ClusteringConcordance <- function(
   dataset <- dataset[keep]
 
   # Prepare data
-  # `tree` may carry tips absent from `dataset` (already dropped from `keep`).
-  # Restrict `splits` to the shared taxon set via `Subsplit()` rather than
-  # pruning `tree` itself: retaining extra tips computes bipartitions over a
-  # different taxon set than `dataset` describes (and column-indexing the
-  # unpruned splits matrix by `keep` cannot recover the correct, smaller set
-  # of splits), but `KeepTip()` would renumber nodes and break any caller
-  # (e.g. `ConcordanceTable()`, `PaintCharacters()`) that matches these split
-  # names against `tree$edge`.
+  # `tree` may carry tips absent from `dataset`; Subsplit restricts `splits` to the
+  # shared taxa without renumbering nodes.
   splits <- as.logical(Subsplit(as.Splits(tree), keep))
-  # `Subsplit()` drops row names entirely when exactly one split survives
-  # restriction (its own version of the drop-to-a-vector bug this file works
-  # around elsewhere) -- recover each surviving split's original node number
-  # by matching it (or its complement) against `tree`'s own splits, likewise
-  # restricted to `keep`'s columns.
+  # Recover each surviving split's original node number
   if (is.null(rownames(splits))) {
     fullRestricted <- as.logical(as.Splits(tree))[, TipLabels(tree) %in% keep,
                                                   drop = FALSE]

R/Consistency.R

@@ -156,19 +156,9 @@ ExpectedLength <- function(dataset, tree, nRelabel = 1000, compress = FALSE) {
     as.integer(intToBits(x)[1:nLevels])
   }, integer(nLevels)))
 
-  # Key on the unlabelled rooted shape, which is what the sampled distribution
-  # is a function of: leaf states are permuted uniformly, and relabelling
-  # composes with a uniform permutation to leave it uniform, so any two trees
-  # of the same shape are sampling the same distribution.  Keying on the
-  # labelled topology instead would be sound but strictly weaker -- identical
-  # topologies are a subset of identical shapes, so it would miss every reuse
-  # this catches and none of its own.  Rooting is part of the shape, as these
-  # characters may contain inapplicable tokens, whose lengths are not
-  # rooting-invariant.
+  # Key on the unlabelled rooted shape
   treeKey <- .ShapeKey(tree)
-  # Cache per shape, and within that per character, rather than pasting both
-  # into one key: that keeps the shape key out of every character's entry, and
-  # leaves no ambiguity about where the shape key ends and the counts begin.
+  # Cache per shape, and within that per character
   treeCache <- .CharLengthCache[[treeKey]]
   if (is.null(treeCache)) {
     treeCache <- new.env(hash = TRUE, parent = emptyenv())
@@ -218,11 +208,7 @@ ExpectedLength <- function(dataset, tree, nRelabel = 1000, compress = FALSE) {
 
 # Canonical identifier of a rooted tree's unlabelled shape, after
 # Aho, Hopcroft & Ullman: a leaf encodes as `01`, and an internal node wraps
-# its children's codes, sorted into a fixed order, in `0`...`1`.  Sorting is
-# what makes the code canonical, so it is already invariant to edge order and
-# to node rotation, and two rooted shapes are isomorphic exactly if their codes
-# agree.  Unlike `TreeTools::RootedTreeShape()`, which enumerates shapes into
-# an integer and so stops at 55 leaves, this is bounded only by string length.
+# its children's codes, sorted into a fixed order, in `0`...`1`.
 # @param tree A rooted, binary tree of class `phylo`.
 # @return A string identifying the shape of `tree`.
 #' @importFrom TreeTools NTip Postorder

R/MaximizeParsimony.R

@@ -122,15 +122,7 @@
   consContrast <- attr(constraint, "contrast")
   nConsStates <- ncol(consContrast)
   if (nConsStates < 2L) {
-    # One state means no taxon is coded `0`, so this is the extreme case of the
-    # inert character warned about below -- and the loudest one, because it is
-    # what `MatrixToPhyDat(c(a = 1, b = 1, c = 1))` produces: a user asking for
-    # a clade and getting no constraint at all.  Warn here rather than returning
-    # silently; the group-size test below never sees these characters.
-    warning("Constraint constrains nothing, and is ignored: no taxon is coded ",
-            "`0`, so every tree separates the `1` taxa from the (empty) `0` ",
-            "group.  Code the taxa that must fall outside the group as `0`.",
-            call. = FALSE)
+    warning("Ignoring empty constraint", call. = FALSE)
     return(list())
   }
 
@@ -150,8 +142,7 @@
   # For each constraint character, record the tips unambiguously in the "1"
   # group (derived state present, ancestral absent) and, separately, those in
   # the "0" group (ancestral present, derived absent).  Tips ambiguous for the
-  # character ("?", or unconstrained taxa) belong to neither group and are free
-  # to plot on either side of the split.
+  # character ("?", or unconstrained taxa) may plot anywhere.
   consSplits <- matrix(0L, nrow = ncol(consMat), ncol = length(constraint))
   consZero   <- matrix(0L, nrow = ncol(consMat), ncol = length(constraint))
   for (ch in seq_len(ncol(consMat))) {
@@ -167,43 +158,20 @@
     }
   }
 
-  # Every tree separates a group of fewer than two taxa from anything: the edge
-  # above a lone tip already does it, and an empty group needs no edge at all.
-  # Such a character constrains nothing under the documented contract, so
-  # enforcing its group as a clade would restrict the search for a guarantee it
-  # already has -- which is the over-strict reading agent-issues/TreeSearch#54
-  # is about.  The test is symmetric in the two groups because they are
-  # interchangeable: which one a user calls "1" is arbitrary, and
-  # build_constraint() swaps them freely to canonicalise.
-  #
-  # Warn rather than drop silently: a character coding only "1" and "?" almost
-  # certainly means "group these taxa", which is not what it says.
+  # Ignore trivial constraints
   nOne <- rowSums(consSplits)
   nZero <- rowSums(consZero)
   inert <- nOne < 2 | nZero < 2
   if (any(inert)) {
-    warning("Constraint character", if (sum(inert) > 1) "s" else "", " ",
-            paste(which(inert), collapse = ", "),
-            if (sum(inert) > 1) " constrain" else " constrains",
-            " nothing, and", if (sum(inert) > 1) " are" else " is",
-            " ignored: every tree separates a group of fewer than two taxa ",
-            "from the rest.  Taxa coded `?` join neither group; code those ",
-            "that must fall outside the group as `0`.", call. = FALSE)
+    warning("Ignoring trivial constraint character", if (sum(inert) > 1) "s" else "", " ",
+            paste(which(inert), collapse = ", "), call. = FALSE)
   }
   keep <- !inert
   consSplits <- consSplits[keep, , drop = FALSE]
   consZero   <- consZero[keep, , drop = FALSE]
   if (nrow(consSplits) == 0L) return(list())
 
   # Every returned tree must display all constraint splits simultaneously.
-  # Two splits are jointly displayable iff they are compatible in the
-  # four-gamete sense: treating each as a bipartition of the tips it
-  # constrains (its "1" group vs its "0" group, ambiguous tips excluded), the
-  # pair is compatible iff at least one of the four group intersections is
-  # empty.  A laminar (nested-or-disjoint) test alone is too strict: it rejects
-  # the case where the two "0" groups are disjoint -- i.e. the splits' "1"
-  # sides jointly cover the constrained tips -- which is perfectly displayable,
-  # e.g. ab | cef and abcd | ef coexist on ((a,b),(d,(c,(e,f)))).
   nSplits <- nrow(consSplits)
   if (nSplits > 1L) {
     for (i in seq_len(nSplits - 1L)) {
@@ -216,7 +184,7 @@
                       !any(aZero & bOne) || !any(aZero & bZero)
         if (!compatible) {
           stop("Constraint is impossible to satisfy: splits ", i, " and ", j,
-               " are incompatible (all four taxon groupings co-occur)")
+               " are incompatible")
         }
       }
     }
@@ -602,7 +570,32 @@
     # through .IwRatchetDepth()'s targetHits/defaultHits escalation (capped at
     # .iwRatchetMaxCycles), which is a genuine reach lever the equal-weights
     # measurement above cannot see.
-    hitMultiplier = if (rung <= 4L) 1L else as.integer(2^(rung - 4L))
+    hitMultiplier = if (rung <= 4L) 1L else as.integer(2^(rung - 4L)),
+    # `enumMaxTrees` multiplier: the size of the returned MPT set, relative to
+    # `poolMaxSize`.  Doubling in step with the other two knobs, from rung 5, so
+    # that one notch keeps meaning "roughly twice the work" on this axis too.
+    #
+    # This scales the ENUMERATION ceiling only, never `poolMaxSize` itself, and
+    # the distinction is the whole point.  During the replicate loop the pool cap
+    # is the size of the working set the search reads -- fuse donors are the
+    # entire pool (uncapped, and taken under the pool mutex on the parallel
+    # path), conflict-guided sector selection reads the pool's split frequencies
+    # once per replicate, and `consensusConstrain` reads its consensus splits --
+    # so scaling it would change which trees the search VISITS.  The
+    # anytime-dominance argument that licenses raising `maxReplicates` above
+    # therefore does NOT transfer to `poolMaxSize`: a bigger pool can delay
+    # every later improvement rather than merely appending to the result.
+    # After the loop, the pool is pure output and a bigger ceiling can only
+    # append equal-score topologies, so the same argument DOES hold there.
+    #
+    # It is bounded in practice without needing a cap: enumeration shares the
+    # `maxSeconds * enumTimeFraction` reserve, and its loop exits as soon as the
+    # pool fills, so an over-generous ceiling costs enumeration time, never a
+    # worse tree.  As with `hitMultiplier` the doubling shape is an OPERATING
+    # POINT rather than a measurement -- what is measured is that the July 2026
+    # 182-tip runs returned exactly `poolMaxSize` trees in all four analyses,
+    # i.e. the ceiling bound the answer rather than the MPT count doing so.
+    enumMultiplier = if (rung <= 4L) 1L else as.integer(2^(rung - 4L))
   )
 }
 
@@ -677,9 +670,13 @@
 #' most-parsimonious tree (\acronym{MPT}) is recovered.
 #' The size of the returned set is bounded by, in order:
 #' \enumerate{
-#'   \item **`poolMaxSize`** (default `100`) — a hard ceiling on the number of
-#'     trees retained.  Raise it (via [`SearchControl()`]) to keep more MPTs;
-#'     with the default you will never see more than 100.
+#'   \item **`enumMaxTrees`**, falling back to **`poolMaxSize`** (default `100`)
+#'     when `enumMaxTrees` is `0` — a hard ceiling on the number of trees
+#'     retained; with the default you will never see more than 100.  Prefer
+#'     raising `enumMaxTrees` (via [`SearchControl()`]): it applies only once the
+#'     search is over, so it cannot alter which trees are visited, whereas
+#'     `poolMaxSize` also sizes the working set that fusing and sectorial search
+#'     read.  From `effort` rung 5 the ladder raises `enumMaxTrees` for you.
 #'   \item **MPT-enumeration time.** After the main search, a TBR plateau walk
 #'     enumerates equal-score neighbours of each pool tree, within a time
 #'     reserve of `maxSeconds * enumTimeFraction`.  If this phase times out it
@@ -851,9 +848,13 @@
 #'       TBR-disconnected islands that random restarts alone miss.}
 #'     \item{4, `large`}{`thorough`'s provisioning with `maxReplicates` raised
 #'       to 500, to suit the higher per-replicate cost of big trees.}
-#'     \item{5 and up}{`thorough`'s provisioning, with both the replicate budget
-#'       and the hit target doubling each notch (1000, 2000, 4000 ...
-#'       replicates), so that one notch always means roughly twice the work.
+#'     \item{5 and up}{`thorough`'s provisioning, with the replicate budget, the
+#'       hit target and the \acronym{MPT}-enumeration ceiling (`enumMaxTrees`)
+#'       all doubling each notch (1000, 2000, 4000 ... replicates), so that one
+#'       notch always means roughly twice the work.  `poolMaxSize` is
+#'       deliberately *not* scaled: it sizes the working set that fusing and
+#'       sectorial search read during the run, so raising it would change which
+#'       trees are visited rather than only how many are returned.
 #'       There is no policy ceiling: extra replicates cannot cost reach, only
 #'       wall, which is what you asked to spend.  The ladder stops only at rung
 #'       26, where the replicate budget outgrows R's integer type.}
@@ -1283,6 +1284,18 @@ MaximizeParsimony <- function(
         targetHits <- as.integer(targetHits * spec[["hitMultiplier"]])
       }
 
+      # Rung-scaled MPT-enumeration ceiling (rung 5 and up).  Keyed off the
+      # POST-merge `poolMaxSize`, so a user who raised the pool gets a
+      # proportionally larger returned set rather than having their value
+      # ignored.  Skipped when the user named `enumMaxTrees` themselves.
+      # `poolMaxSize` is deliberately not touched -- see .RungSpec().
+      if (!("enumMaxTrees" %in% union(names(controlDots),
+                                      attr(control, "explicit"))) &&
+          spec[["enumMultiplier"]] > 1L) {
+        control[["enumMaxTrees"]] <-
+          as.integer(control[["poolMaxSize"]] * spec[["enumMultiplier"]])
+      }
+
       # Implied-weights ratchet depth. Under implied weights the optimum often
       # sits in a small basin at fine score resolution, separated from an
       # easy-to-find near-optimum by a fraction of a step; character reweighting

R/ParsSim.R

@@ -132,11 +132,7 @@ ParsSim <- function(tree,
   # --- Determine state counts per character ----------------------------------
   n_states_vec <- rep(seq_along(nChar) + 1L, times = nChar)
   if (any(n_states_vec > 31L)) {
-    stop("ParsSim() supports at most 31 states per character (state codes ",
-         "0:30): the internal Fitch bit-set representation packs states ",
-         "into a 32-bit integer via bitwShiftL(), which silently overflows ",
-         "to NA beyond that. Requested up to ", max(n_states_vec),
-         " states via `nChar`.")
+    stop("ParsSim() supports at most 31 states per character.")
   }
 
   # --- Validate and expand rootState ------------------------------------------

R/SearchControl.R

@@ -123,8 +123,22 @@
 #'   within each replicate, after TBR polish.  This approximates TNT's
 #'   within-replicate fusing pattern. Default: `FALSE`.
 #' @param poolMaxSize Integer; maximum trees retained in the pool.
+#'   This governs the pool throughout the search, so it is not only a ceiling on
+#'   the trees returned: fuse draws its donors from the pool, and conflict-guided
+#'   sector selection and `consensusConstrain` both read it.  Raising it
+#'   therefore changes the trajectory as well as the output; to keep more
+#'   most-parsimonious trees without that side effect, use `enumMaxTrees`.
 #' @param poolSuboptimal Numeric; retain trees that are this many steps
 #'   worse than the best tree.  0 (default) keeps only optimal trees.
+#' @param enumMaxTrees Integer; retention ceiling applied to the
+#'   \acronym{MPT}-enumeration phase alone, which runs after the last replicate.
+#'   `0` (default) keeps `poolMaxSize` throughout, reproducing the behaviour
+#'   before this argument existed.  Because enumeration happens once the search
+#'   is over, a larger ceiling here only appends further equal-score topologies:
+#'   it cannot change which trees the search visits.  This is the knob
+#'   [`MaximizeParsimony()`]'s `effort` scales; `poolMaxSize` is deliberately
+#'   left alone.  Values below `poolMaxSize` are ignored (the ceiling is only
+#'   ever raised).
 #' @param consensusStableReps Integer; stop when the strict consensus of
 #'   best-score pool trees has been unchanged for this many consecutive
 #'   replicates.
@@ -377,7 +391,8 @@ SearchControl <- function(
     # sampling from {Wagner-random, Wagner-Goloboff, Wagner-entropy,
     # random-tree, pool-ratchet, pool-NNI-perturb}. Overrides wagnerBias.
     adaptiveStart = FALSE,
-    enumTimeFraction = 0.1
+    enumTimeFraction = 0.1,
+    enumMaxTrees = 0L
 ) {
   # Record which fields the caller set explicitly (by name or position;
   # `match.call()` normalises positional args to their names).  This lets
@@ -398,6 +413,14 @@ SearchControl <- function(
       stop("`", .p, "` must be a single positive integer")
     }
   }
+  # `enumMaxTrees` takes 0 ("follow poolMaxSize") but never a negative: the
+  # kernel only ever RAISES the ceiling, so a negative would be silently inert
+  # rather than reported, hiding a sign typo.
+  .emt <- as.integer(enumMaxTrees)
+  if (length(.emt) != 1L || is.na(.emt) || .emt < 0L) {
+    stop("`enumMaxTrees` must be a single non-negative integer ",
+         "(0 follows `poolMaxSize`)")
+  }
   # `stopPatience` is a replicate count, so a negative value is meaningless; the
   # kernel treats anything <= 0 as "off", which would silently ignore a typo
   # such as -20 rather than honouring the obvious intent.
@@ -498,7 +521,8 @@ SearchControl <- function(
       annealTEnd = as.double(annealTEnd),
       annealMovesPerPhase = as.integer(annealMovesPerPhase),
       adaptiveStart = as.logical(adaptiveStart),
-      enumTimeFraction = as.double(enumTimeFraction)
+      enumTimeFraction = as.double(enumTimeFraction),
+      enumMaxTrees = .emt
     ),
     class = "SearchControl",
     explicit = .explicit
@@ -530,7 +554,7 @@ print.SearchControl <- function(x, ...) {
                      "sectorCombStarts", "sectorFuseRounds",
                      "postRatchetSectorial"),
     "Fuse/Pool" = c("fuseInterval", "fuseAcceptEqual", "intraFuse",
-                     "poolMaxSize", "poolSuboptimal"),
+                     "poolMaxSize", "poolSuboptimal", "enumMaxTrees"),
     "Stopping" = c("consensusStableReps", "perturbStopFactor", "stopPatience",
                     "adaptiveLevel",
                     "consensusConstrain", "adaptiveStart",

R/ts-driven-compat.R

@@ -42,6 +42,7 @@ ts_driven_search <- function(
     fuseAcceptEqual = FALSE,
     poolMaxSize = 100L,
     poolSuboptimal = 0.0,
+    enumMaxTrees = 0L,
     maxSeconds = 0.0,
     verbosity = 0L,
     min_steps = integer(0),
@@ -140,7 +141,8 @@ ts_driven_search <- function(
     pruneReinsertDrop = as.double(pruneReinsertDrop),
     pruneReinsertSelection = as.integer(pruneReinsertSelection),
     adaptiveStart = as.logical(adaptiveStart),
-    enumTimeFraction = as.double(enumTimeFraction)
+    enumTimeFraction = as.double(enumTimeFraction),
+    enumMaxTrees = as.integer(enumMaxTrees)
   )
 
   # Anneal config: fold into SearchControl if provided