fix: loosen constraint matching to the documented free-taxa contract

.AGENTS/memory/architecture.md

@@ -94,8 +94,34 @@ Profile mode sets `ds.concavity = 1.0` (finite sentinel) so existing
 
 ## Constraint enforcement
 
+- A constraint split names **two disjoint groups** plus a FREE remainder. A tree
+  satisfies it iff some edge separates the groups; free tips may fall either
+  side. This is what `?MaximizeParsimony`'s `constraint` documents and, since
+  agent-issues/TreeSearch#54, what every entry point enforces. `ts::node_displays_split()`
+  (`ts_constraint.h`) is THE shared predicate — `map_constraint_nodes()`,
+  `wagner_tree_displays_constraint()` and `ts_collapse_pool()` all call it.
+  Reintroducing an exact-clade test at any one of them freezes replicates.
 - `build_constraint()` reads R split matrix with **column-major** indexing:
-  `split_matrix[s + n_splits * t]`.
+  `split_matrix[s + n_splits * t]`. Values: `1` = together-group, `0` =
+  apart-group, **anything else (`NA_INTEGER`) = free**. A hand-built 0/1 matrix
+  therefore means "no free tips" and reduces to the exact-clade behaviour, which
+  is what `build_constraint_from_bitsets()` (consensus constraints) relies on.
+- `ConstraintData` carries `split_zeros` (the apart-group) alongside
+  `split_tips`, and both ends of the displaying-node chain:
+  `constraint_node` (tightest, used for "must land outside") and
+  `constraint_node_hi` (highest, "must land inside"). Any writer of one must
+  write the other — `ts_wagner.cpp` pins hi to the tight anchor.
+- `.PrepareConstraint()` drops (and warns about) a character whose `1` **or**
+  `0` group holds fewer than two taxa: vacuous under the documented contract,
+  since every tree separates such a group from the rest.
+- A user constraint binds at three boundaries besides the per-move filter
+  (agent-issues/TreeSearch#59): the start tree is repaired by `impose_constraint()`
+  before it is scored, each replicate's finished tree is gated by
+  `capture_satisfies_constraint()` on its way into the pool, and the enforced
+  split is kept out of the final collapse. `impose_constraint()` is heuristic and
+  can fail, so every caller re-verifies. The collapse protects a realising node
+  only when no other realising node already survives — protecting
+  unconditionally would resolve a branch the constraint never asked for.
 - Wagner uses LCA-based constraint mapping (`wagner_map_constraint_nodes`)
   since splits aren't fully present during incremental construction.
 - Wagner has a posthoc retry loop (up to 100 random addition orders) as a

NEWS.md

@@ -18,6 +18,29 @@
   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.
+- Every part of the search now reads `constraint` the way it is documented: a
+  tree is compatible with a constraint character when some edge separates the
+  taxa coded `1` from those coded `0`, with `?`-coded and unmentioned taxa free
+  to fall on either side.  The locked-node filter that screens individual
+  rearrangements, the constrained Wagner build and the collapse pass previously
+  required the `1` group to be a clade *exactly*, free taxa excluded.  That is
+  strictly stronger, so the search never accepted a rearrangement that broke the
+  documented constraint; but a start tree that satisfied the documented
+  constraint without making either group an exact clade matched no node, every
+  rearrangement was rejected, and the replicate returned its start unimproved.
+  Constrained searches with `?`-coded taxa therefore reach better scores.
+  The exact match also blunted the collapse protection described above: with
+  free taxa it matched no branch, so the separating edge could still be
+  contracted away -- the one route by which a *returned* tree could break the
+  constraint.
+- A constraint character whose `1` or `0` group holds fewer than two taxa now
+  warns and is ignored, rather than being enforced as a clade.  Every tree
+  separates such a group from the rest, so the character constrains nothing
+  under the documented reading.  The test is symmetric in the two groups, which
+  the old one was not: `c(a = 1, b = 1, c = 0)` and `c(a = 0, b = 0, c = 1)`
+  state the same constraint and are now treated the same way.  Code the taxa
+  that must fall outside a group as `0`, rather than leaving them `?`, to keep
+  it enforced.
 - `TreeLength()`, `CharacterLength()`, `TreeScore()` and `EdgeListScore()` -- and
   so `Consistency()`, `ExpectedLength()`, `ConcordantInformation()`,
   `LengthAdded()` and `SuccessiveApproximations()`, which score trees through

R/MaximizeParsimony.R

@@ -121,7 +121,18 @@
 
   consContrast <- attr(constraint, "contrast")
   nConsStates <- ncol(consContrast)
-  if (nConsStates < 2L) return(list())
+  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)
+    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
@@ -156,10 +167,30 @@
     }
   }
 
-  keep <- apply(consSplits, 1, function(row) {
-    s <- sum(row)
-    s >= 1 && s < length(constraint) - 1
-  })
+  # 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.
+  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)
+  }
+  keep <- !inert
   consSplits <- consSplits[keep, , drop = FALSE]
   consZero   <- consZero[keep, , drop = FALSE]
   if (nrow(consSplits) == 0L) return(list())
@@ -199,6 +230,15 @@
   consTipData <- matrix(unlist(constraint, use.names = FALSE),
                         nrow = length(constraint), byrow = TRUE)
 
+  # Fold the two groups into the single membership matrix the C++ engine reads:
+  # 1 = "together", 0 = "apart", NA = free to fall on either side.  A tip that
+  # is in neither group must not be coded 0, or the engine would enforce the
+  # stricter "the 1 group is an exact clade" reading and refuse to move a start
+  # tree that already satisfies the documented one (agent-issues/TreeSearch#54).
+  # build_constraint() (src/ts_constraint.cpp) treats any value that is neither
+  # 1 nor 0 as free, so a plain 0/1 matrix still means "no free tips".
+  consSplits[consSplits == 0L & consZero == 0L] <- NA_integer_
+
   list(
     consSplitMatrix = consSplits,
     consZero = consZero,
@@ -234,6 +274,11 @@
 # neither group here and are unconstrained there too (.PrepareConstraint()
 # warns about that at input).
 .ConstraintViolated <- function(tree, consOne, consZero) {
+  # `consOne` is the membership matrix the C++ kernels read, so it uses their
+  # coding: 1 = in the group, anything else -- including the NA that marks a
+  # free tip -- out of it.  Reduce it to 0/1 here rather than let an NA
+  # propagate through the accumulation below and turn every comparison NA.
+  consOne <- (!is.na(consOne) & consOne == 1L) * 1L
   edge <- Postorder(tree)[["edge"]]
   parent <- edge[, 1L]
   child <- edge[, 2L]
@@ -769,6 +814,14 @@
 #' returned trees will be perfectly compatible with each character in
 #' `constraint`; or a tree of class `phylo`, all of whose nodes will occur
 #' in any output tree.
+#' A returned tree is compatible with a constraint character when some edge
+#' separates the taxa coded `1` from those coded `0`.  Taxa coded `?`, and taxa
+#' that `constraint` does not mention, are unconstrained: they may fall on
+#' either side of that edge, and are not required to join either group.
+#' A character whose `1` or `0` group contains fewer than two taxa therefore
+#' constrains nothing -- every tree separates such a group from the rest -- and
+#' is ignored with a warning.  To group taxa, code the taxa they must be
+#' separated from as `0` rather than leaving them `?`.
 #' 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
@@ -1748,21 +1801,17 @@ MaximizeParsimony <- function(
     # enforced clade"): a constraint is external evidence for a grouping the
     # 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.
+    # are the enforced bipartitions in tip_data order, carrying both groups
+    # (1 = together, 0 = apart, NA = free; see .PrepareConstraint).  The kernel
+    # needs both: a tree with free tips generally realises the split at a node
+    # whose tip set is wider than the 1 group, which no exact match reaches, so
+    # the enforced grouping would collapse 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, consZero
+      scoringConfig, hsjConfig, xformConfig, consSplits
     )
     outTrees <- lapply(collapsed$trees, function(edgeMat) {
       tr <- list(

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, consZero = NULL) {
-    .Call(`_TreeSearch_ts_collapse_pool`, edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix, consZero)
+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_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()) {

inst/WORDLIST

@@ -266,6 +266,7 @@ rearranger
 reconverged
 reconverges
 regraft
+regrafted
 regrafting
 regrafts
 reoptimisation

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, 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) {
+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) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -604,8 +604,7 @@ 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::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));
+    rcpp_result_gen = Rcpp::wrap(ts_collapse_pool(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix));
     return rcpp_result_gen;
 END_RCPP
 }

src/TreeSearch-init.c

@@ -60,7 +60,7 @@ extern SEXP _TreeSearch_ts_ev_cache_key_probe(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP
 extern SEXP _TreeSearch_ts_ls_fit(SEXP, SEXP, SEXP, SEXP);
 extern SEXP _TreeSearch_ts_ls_search(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP _TreeSearch_ts_collapsed_flags_debug(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
-extern SEXP _TreeSearch_ts_collapse_pool(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
+extern SEXP _TreeSearch_ts_collapse_pool(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 
 static const R_CallMethodDef callMethods[] = {
   {"_TreeSearch_nni",               (DL_FUNC) &_TreeSearch_nni, 3},
@@ -116,7 +116,7 @@ static const R_CallMethodDef callMethods[] = {
   {"_TreeSearch_ts_ls_fit", (DL_FUNC) &_TreeSearch_ts_ls_fit, 4},
   {"_TreeSearch_ts_ls_search", (DL_FUNC) &_TreeSearch_ts_ls_search, 6},
   {"_TreeSearch_ts_collapsed_flags_debug", (DL_FUNC) &_TreeSearch_ts_collapsed_flags_debug, 6},
-  {"_TreeSearch_ts_collapse_pool", (DL_FUNC) &_TreeSearch_ts_collapse_pool, 10},
+  {"_TreeSearch_ts_collapse_pool", (DL_FUNC) &_TreeSearch_ts_collapse_pool, 9},
   {NULL, NULL, 0}
 };