fix: enforce user constraints at the start-tree, pool-capture and collapse boundaries

NEWS.md

@@ -1,5 +1,24 @@
 # To integrate into 2.0.0 notes
 
+- `constraint` now binds the trees `MaximizeParsimony()` returns, at three
+  boundaries where it did not.  A starting tree supplied through `tree` was
+  never checked against the constraint; because a constrained search rejects
+  every rearrangement away from a violating tree, the replicate froze on it and
+  reported a score no constraint-satisfying tree could reach, which then evicted
+  the compliant trees other replicates had found.  A violating start is now
+  rearranged until it complies before the search begins, **with a warning**.
+  Separately, a replicate's own tree entered the pool unchecked, and the final
+  collapse of unsupported branches could contract the very branch that displayed
+  an enforced grouping -- so under the default `collapse = TRUE` a returned tree
+  could break the constraint outright.  Both paths are now checked.
+
+  **Constrained results may therefore differ from previous versions**: scores
+  can rise to the true constrained optimum, and returned trees will display the
+  constrained groupings.  `MaximizeParsimony()` also warns if any replicate
+  ended on a tree that could not be made to satisfy the constraint, and now
+  raises an error rather than returning an unverified tree if no
+  constraint-satisfying tree was found at all.
+
 - `inapplicable = "xform"` scores are now reported at a canonical rooting, so a
   reported score is reproducible.  The x-transformation's step matrix is
   asymmetric -- a gain costs one more than the number of secondary characters it

R/AdditionTree.R

@@ -142,7 +142,9 @@ AdditionTree <- function(dataset, concavity = Inf, constraint, sequence) {
     addition_order = addition_order,
     concavity = as.double(concavity)
   )
-  result <- do.call(ts_wagner_tree, c(searchArgs, consArgs, profileArgs))
+  result <- do.call(ts_wagner_tree,
+                    c(searchArgs, .KernelConstraintArgs(consArgs),
+                      profileArgs))
 
   # Reconstruct phylo from edge matrix
   tree <- list(

R/MaximizeParsimony.R

@@ -123,6 +123,17 @@
   nConsStates <- ncol(consContrast)
   if (nConsStates < 2L) return(list())
 
+  # Constraints are enforced as bipartitions, so only the two extreme states of
+  # a character are read: taxa carrying an intermediate state are in neither
+  # group and go unconstrained.  Say so rather than let the caller infer, from
+  # `@param constraint`'s "compatible with each character", that a third state
+  # groups its taxa too.
+  if (nConsStates > 2L) {
+    warning("`constraint` characters with more than two states are enforced ",
+            "as the split between their first and last state only; taxa in ",
+            "any intermediate state are left unconstrained.", call. = FALSE)
+  }
+
   consMat <- matrix(unlist(constraint, use.names = FALSE),
                     nrow = length(constraint), byrow = TRUE)
   # For each constraint character, record the tips unambiguously in the "1"
@@ -190,6 +201,7 @@
 
   list(
     consSplitMatrix = consSplits,
+    consZero = consZero,
     consContrast = consContrast,
     consTipData = consTipData,
     consWeight = as.integer(consWeight),
@@ -198,6 +210,63 @@
   )
 }
 
+# Constraint fields the flat `ts_*` kernels declare as formals, in contrast to
+# the list-config entry points, which ignore anything they do not name.  A
+# `do.call()` onto a flat kernel has to be filtered through this, or a field
+# added for the list-config path becomes an unused-argument error there.
+.kernelConsFields <- c("consSplitMatrix", "consContrast", "consTipData",
+                       "consWeight", "consLevels", "consExpectedScore")
+
+.KernelConstraintArgs <- function(consArgs) {
+  consArgs[intersect(names(consArgs), .kernelConsFields)]
+}
+
+# Does `tree` display a split separating a constraint character's "1" group
+# from its "0" group?  This is the phyDat reading `constraint` is documented
+# in: tips ambiguous for the character sit on either side, so the test is
+# "some edge separates the two groups", not the stricter "the 1 group is
+# exactly a clade" that the search's locked-node machinery enforces
+# internally.  `consOne` / `consZero` are .PrepareConstraint()'s matrices, in
+# `tip_data` column order; `tree`'s tips must already be renumbered to match.
+#
+# The two groups are the character's extreme states, so this answers for
+# exactly what the engine enforces -- an intermediate state's taxa are in
+# neither group here and are unconstrained there too (.PrepareConstraint()
+# warns about that at input).
+.ConstraintViolated <- function(tree, consOne, consZero) {
+  edge <- Postorder(tree)[["edge"]]
+  parent <- edge[, 1L]
+  child <- edge[, 2L]
+  nTip <- ncol(consOne)
+  nRow <- nrow(consOne)
+  # One accumulation pass carries every group at once: rows 1..nRow are the
+  # "1" groups, the rest the "0" groups.  Nodes index the COLUMNS, so each
+  # accumulation touches one contiguous stretch of a column-major matrix.
+  counts <- matrix(0L, nrow = 2L * nRow, ncol = max(edge))
+  counts[, seq_len(nTip)] <- rbind(consOne, consZero)
+  for (i in seq_along(parent)) {
+    counts[, parent[i]] <- counts[, parent[i]] + counts[, child[i]]
+  }
+  # Postorder lists every node before its parent, so the first node holding a
+  # whole group is that group's MRCA; the groups are separated iff one MRCA
+  # holds none of the other group.
+  nodes <- c(child, parent[length(parent)])
+  for (r in seq_len(nRow)) {
+    one <- counts[r, ]
+    zero <- counts[nRow + r, ]
+    nOne <- sum(consOne[r, ])
+    nZero <- sum(consZero[r, ])
+    mrcaOne <- nodes[one[nodes] == nOne][1]
+    mrcaZero <- nodes[zero[nodes] == nZero][1]
+    displayed <- (!is.na(mrcaOne) && zero[mrcaOne] == 0L) ||
+      (!is.na(mrcaZero) && one[mrcaZero] == 0L)
+    if (!displayed) {
+      return(TRUE)
+    }
+  }
+  FALSE
+}
+
 # Ratchet depth for implied weights under `thorough`/`large`, applied after the
 # strategy preset (see MaximizeParsimony()). Kept out of `.StrategyPresets()` so
 # the preset table stays scorer-agnostic: this depth is calibrated for implied
@@ -615,6 +684,12 @@
 #'   says so in a warning.  Raise `targetHits` as well as `maxReplicates` to
 #'   use more of it.
 #'   If unspecified, all replicates start from random Wagner trees.
+#'   A start tree that does not satisfy `constraint` is rearranged until it
+#'   does before the search begins, with a warning: `constraint` is a
+#'   guarantee about the trees returned, whereas `tree` only says where to
+#'   begin, so when the two conflict the guarantee wins.  A taxon coded `?`
+#'   for a constraint character is unconstrained by it and may start on
+#'   either side of that split.
 #'   Edge lengths are not supported and will be deleted.
 #'   Rooted and unrooted trees are both accepted; an unrooted tree is rooted
 #'   arbitrarily (on its first tip) before the search begins, which may
@@ -696,6 +771,9 @@
 #' in any output tree.
 #' Constraint searches are supported natively: all tree rearrangements
 #' are filtered to respect the constraint topology.
+#' Each constraint character is enforced as a single split, so one with more
+#' than two states is read as the split between its first and last state
+#' alone: taxa in an intermediate state are left unconstrained, with a warning.
 #' @param effort Integer: how much search effort to spend, **relative to the
 #'   amount chosen automatically** for this dataset.  `0` (the default) accepts
 #'   the automatic choice; `1` asks for one notch more, `-1` one notch less.
@@ -1520,6 +1598,23 @@ MaximizeParsimony <- function(
     cli_alert_info("Constraint: {nrow(consArgs$consSplitMatrix)} split{?s}")
   }
 
+  # A start tree that breaks the constraint is not something the search can
+  # rearrange its way out of -- every constrained move from it is rejected, so
+  # it would freeze the replicate on a tree scoring better than any legal one.
+  # The engine repairs such a start before scoring it, but the conflict is the
+  # caller's to know about: either `tree` or `constraint` is not what they
+  # meant, and the tree they get back will not be the one they supplied.
+  if (userTree && length(consArgs) > 0L) {
+    violating <- vapply(startTrees, .ConstraintViolated, logical(1),
+                        consArgs[["consSplitMatrix"]], consArgs[["consZero"]])
+    if (any(violating)) {
+      warning(sum(violating), " of the ", length(startTrees),
+              " tree(s) supplied to `tree` do not satisfy `constraint`; ",
+              "they will be rearranged to comply before the search starts, ",
+              "or replaced if that fails.", call. = FALSE)
+    }
+  }
+
   # --- Profile parsimony: extract info_amounts ---
   profileArgs <- list()
   if (useProfile) {
@@ -1654,12 +1749,20 @@ MaximizeParsimony <- function(
     # matrix doesn't capture, so it stays visible even at zero length, while the
     # unsupported non-constraint branches still collapse.  consSplitMatrix rows
     # are the enforced bipartitions in tip_data order (see .PrepareConstraint).
+    # `consZero` names the tips the constraint places on the far side of the
+    # split; tips ambiguous for the character are in neither group.  Without it
+    # the kernel can only recognise a node whose tip set is the 1 group exactly,
+    # and a split realised by any wider node goes unprotected -- collapsing the
+    # enforced grouping out of the returned tree.
     consSplits <- if (!is.null(constraintConfig)) {
       constraintConfig[["consSplitMatrix"]]
     }
+    consZero <- if (!is.null(constraintConfig)) {
+      constraintConfig[["consZero"]]
+    }
     collapsed <- ts_collapse_pool(
       bestTrees, contrast, tip_data, weight, levels,
-      scoringConfig, hsjConfig, xformConfig, consSplits
+      scoringConfig, hsjConfig, xformConfig, consSplits, consZero
     )
     outTrees <- lapply(collapsed$trees, function(edgeMat) {
       tr <- list(
@@ -1680,6 +1783,18 @@ MaximizeParsimony <- function(
     })
   }
   if (length(outTrees) == 0L) {
+    # `treeTpl` is a starting tree, so under a constraint it is exactly what may
+    # not be handed back unchecked: an empty pool means no replicate produced a
+    # tree the constraint gate accepted -- or, benignly, that the time limit
+    # expired before the first one finished.  Check rather than assume, so a
+    # short budget still returns a tree when the fallback happens to comply.
+    if (!is.null(constraintConfig) &&
+        .ConstraintViolated(treeTpl, constraintConfig[["consSplitMatrix"]],
+                            constraintConfig[["consZero"]])) {
+      stop("The search returned no tree satisfying `constraint`. Check that ",
+           "the constraint is compatible with the data, and allow more search ",
+           "with `maxReplicates` or `maxSeconds`.")
+    }
     outTrees <- list(treeTpl)
   }
 

R/RcppExports.R

@@ -197,8 +197,8 @@ ts_driven_search <- function(contrast, tip_data, weight, levels, searchControl,
     .Call(`_TreeSearch_ts_driven_search`, contrast, tip_data, weight, levels, searchControl, runtimeConfig, scoringConfig, constraintConfig, hsjConfig, xformConfig)
 }
 
-ts_collapse_pool <- function(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig = NULL, xformConfig = NULL, consSplitMatrix = NULL) {
-    .Call(`_TreeSearch_ts_collapse_pool`, edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix)
+ts_collapse_pool <- function(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig = NULL, xformConfig = NULL, consSplitMatrix = NULL, consZero = NULL) {
+    .Call(`_TreeSearch_ts_collapse_pool`, edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix, consZero)
 }
 
 ts_resample_search <- function(contrast, tip_data, weight, levels, bootstrap = FALSE, jackProportion = 2.0 / 3.0, maxReplicates = 5L, targetHits = 2L, tbrMaxHits = 1L, ratchetCycles = 3L, ratchetPerturbProb = 0.04, driftCycles = 0L, min_steps = integer(), concavity = -1.0, consSplitMatrix = NULL, consContrast = NULL, consTipData = NULL, consWeight = NULL, consLevels = NULL, consExpectedScore = 0L, infoAmounts = NULL, xpiwe = FALSE, xpiwe_r = 0.5, xpiwe_max_f = 5.0, obs_count = integer()) {

R/Resample.R

@@ -369,7 +369,7 @@ Resample <- function(dataset, tree, method = "jack", proportion = 2 / 3,
     # Batch mode: run all replicates at once (optionally in parallel)
     batchArgs <- c(searchArgs,
                    list(nReplicates = nReplicates, nThreads = nThreads),
-                   consArgs, profileArgs)
+                   .KernelConstraintArgs(consArgs), profileArgs)
     result <- do.call(ts_parallel_resample, batchArgs)
 
     trees <- vector("list", nReplicates)
@@ -393,7 +393,9 @@ Resample <- function(dataset, tree, method = "jack", proportion = 2 / 3,
   }
 
   # Single-replicate path (original behavior)
-  result <- do.call(ts_resample_search, c(searchArgs, consArgs, profileArgs))
+  result <- do.call(ts_resample_search,
+                    c(searchArgs, .KernelConstraintArgs(consArgs),
+                      profileArgs))
 
   if (nrow(result$edge) == 0L) {
     tr <- if (!missing(tree) && inherits(tree, "phylo")) tree

R/SuccessiveApproximations.R

@@ -104,7 +104,9 @@ SuccessiveApproximations <- function (tree, dataset, outgroup = NULL, k = 3,
     xpiwe_max_f = as.double(xpiwe_max_f),
     obs_count = if (useXpiwe) obsCount else integer(0)
   )
-  result <- do.call(ts_successive_approx, c(searchArgs, consArgs, profileArgs))
+  result <- do.call(ts_successive_approx,
+                    c(searchArgs, .KernelConstraintArgs(consArgs),
+                      profileArgs))
 
   if (result$converged && verbosity > 0) {
     message("Successive approximations converged after ",

src/RcppExports.cpp

@@ -590,8 +590,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // ts_collapse_pool
-List ts_collapse_pool(List edges, NumericMatrix contrast, IntegerMatrix tip_data, IntegerVector weight, CharacterVector levels, List scoringConfig, Nullable<List> hsjConfig, Nullable<List> xformConfig, Nullable<IntegerMatrix> consSplitMatrix);
-RcppExport SEXP _TreeSearch_ts_collapse_pool(SEXP edgesSEXP, SEXP contrastSEXP, SEXP tip_dataSEXP, SEXP weightSEXP, SEXP levelsSEXP, SEXP scoringConfigSEXP, SEXP hsjConfigSEXP, SEXP xformConfigSEXP, SEXP consSplitMatrixSEXP) {
+List ts_collapse_pool(List edges, NumericMatrix contrast, IntegerMatrix tip_data, IntegerVector weight, CharacterVector levels, List scoringConfig, Nullable<List> hsjConfig, Nullable<List> xformConfig, Nullable<IntegerMatrix> consSplitMatrix, Nullable<IntegerMatrix> consZero);
+RcppExport SEXP _TreeSearch_ts_collapse_pool(SEXP edgesSEXP, SEXP contrastSEXP, SEXP tip_dataSEXP, SEXP weightSEXP, SEXP levelsSEXP, SEXP scoringConfigSEXP, SEXP hsjConfigSEXP, SEXP xformConfigSEXP, SEXP consSplitMatrixSEXP, SEXP consZeroSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -604,7 +604,8 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< Nullable<List> >::type hsjConfig(hsjConfigSEXP);
     Rcpp::traits::input_parameter< Nullable<List> >::type xformConfig(xformConfigSEXP);
     Rcpp::traits::input_parameter< Nullable<IntegerMatrix> >::type consSplitMatrix(consSplitMatrixSEXP);
-    rcpp_result_gen = Rcpp::wrap(ts_collapse_pool(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix));
+    Rcpp::traits::input_parameter< Nullable<IntegerMatrix> >::type consZero(consZeroSEXP);
+    rcpp_result_gen = Rcpp::wrap(ts_collapse_pool(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix, consZero));
     return rcpp_result_gen;
 END_RCPP
 }