T-200: Decouple R-loop search from MorphyLib; add concavity parameter

DESCRIPTION

@@ -91,6 +91,7 @@ Collate:
   'SearchControl.R'
   'MaximizeParsimony.R'
   'Morphy.R'
+  'NativeSearch.R'
   'NNI.R'
   'ParsSim.R'
   'PlotCharacter.R'

NAMESPACE

@@ -28,6 +28,7 @@ export(C_MorphyLength)
 export(Carter1)
 export(CharacterHierarchy)
 export(CharacterLength)
+export(CleanNativeData)
 export(ClusterStrings)
 export(ClusteringConcordance)
 export(ConcordanceTable)
@@ -69,13 +70,16 @@ export(MultiRatchet)
 export(MutualClusteringConcordance)
 export(NNI)
 export(NNISwap)
+export(NativeBootstrap)
+export(NativeLength)
 export(ParsSim)
 export(PhyDat2Morphy)
 export(PhylogeneticConcordance)
 export(PlotCharacter)
 export(PolEscapa)
 export(PrepareDataIW)
 export(PrepareDataProfile)
+export(PrepareNativeData)
 export(PresCont)
 export(QACol)
 export(QALegend)

R/Bootstrap.R

@@ -16,6 +16,7 @@
 MorphyBootstrap <- function (edgeList, morphyObj, EdgeSwapper = NNISwap, 
                              maxIter, maxHits, verbosity = 1L, 
                              stopAtPeak = FALSE, stopAtPlateau=0L, ...) {
+  .Deprecated("NativeBootstrap")
   startWeights <- MorphyWeights(morphyObj)["exact", ]
   eachChar <- seq_along(startWeights)
   deindexedChars <- rep.int(eachChar, startWeights)

R/CustomSearch.R

@@ -12,7 +12,7 @@
 #' @keywords internal
 #' @export
 EdgeListSearch <- function (edgeList, dataset,
-                            TreeScorer = MorphyLength,
+                            TreeScorer = NativeLength,
                             EdgeSwapper = RootedTBRSwap,
                             maxIter = 100, maxHits = 20, 
                             bestScore = NULL, stopAtScore = NULL, 
@@ -171,10 +171,10 @@ EdgeListSearch <- function (edgeList, dataset,
 #' @family custom search functions
 #' @importFrom TreeTools RenumberTips
 #' @export
-TreeSearch <- function (tree, dataset,
-                        InitializeData = PhyDat2Morphy,
-                        CleanUpData    = UnloadMorphy,
-                        TreeScorer     = MorphyLength,
+TreeSearch <- function (tree, dataset, concavity = Inf,
+                        InitializeData = PrepareNativeData,
+                        CleanUpData    = CleanNativeData,
+                        TreeScorer     = NativeLength,
                         EdgeSwapper    = RootedTBRSwap,
                         maxIter = 100L, maxHits = 20L,
                         stopAtPeak = FALSE, stopAtPlateau = 0L,
@@ -186,7 +186,11 @@ TreeSearch <- function (tree, dataset,
   edgeList <- tree[["edge"]]
   edgeList <- RenumberEdges(edgeList[, 1], edgeList[, 2])
 
-  initializedData <- InitializeData(dataset)
+  if (identical(InitializeData, PrepareNativeData)) {
+    initializedData <- PrepareNativeData(dataset, concavity = concavity)
+  } else {
+    initializedData <- InitializeData(dataset)
+  }
   on.exit(initializedData <- CleanUpData(initializedData))
 
   bestScore <- attr(tree, "score")

R/Jackknife.R

@@ -19,10 +19,11 @@
 #' @family split support functions
 #' @family custom search functions
 #' @export
-Jackknife <- function(tree, dataset, resampleFreq = 2 / 3,
-                      InitializeData = PhyDat2Morphy,
-                      CleanUpData    = UnloadMorphy,
-                      TreeScorer     = MorphyLength,
+Jackknife <- function(tree, dataset, concavity = Inf,
+                      resampleFreq = 2 / 3,
+                      InitializeData = PrepareNativeData,
+                      CleanUpData    = CleanNativeData,
+                      TreeScorer     = NativeLength,
                       EdgeSwapper    = TBRSwap,
                       jackIter = 5000L, searchIter = 4000L, searchHits = 42L,
                       verbosity = 1L, ...) {
@@ -34,19 +35,29 @@ Jackknife <- function(tree, dataset, resampleFreq = 2 / 3,
   edgeList <- tree[["edge"]]
   edgeList <- RenumberEdges(edgeList[, 1], edgeList[, 2])
 
-  morphyObj <- InitializeData(dataset)
-  on.exit(morphyObj <- CleanUpData(morphyObj))
+  if (identical(InitializeData, PrepareNativeData)) {
+    initializedData <- PrepareNativeData(dataset, concavity = concavity)
+  } else {
+    initializedData <- InitializeData(dataset)
+  }
+  on.exit(initializedData <- CleanUpData(initializedData))
 
-  startWeights <- MorphyWeights(morphyObj)["exact", ]
+  useMorphy <- inherits(initializedData, "morphyPtr")
+
+  if (useMorphy) {
+    startWeights <- MorphyWeights(initializedData)["exact", ]
+  } else {
+    startWeights <- initializedData[["original_weight"]]
+  }
   eachChar <- seq_along(startWeights)
   deindexedChars <- rep.int(eachChar, startWeights)
   charsToKeep <- ceiling(resampleFreq * length(deindexedChars))
 
   if (charsToKeep < 1L) {
-    stop("resampleFreq of ", resampleFreq, " is too low; can't keep 0 of ",
+    stop("resampleFreq of ", resampleFreq, " is too low; cannot keep 0 of ",
          length(deindexedChars), " characters.")
   } else if (charsToKeep >= length(deindexedChars)) {
-    stop("resampleFreq of ", resampleFreq, " is too high; can't keep all ",
+    stop("resampleFreq of ", resampleFreq, " is too high; cannot keep all ",
          length(deindexedChars), " characters.")
   }
 
@@ -61,16 +72,25 @@ Jackknife <- function(tree, dataset, resampleFreq = 2 / 3,
     } #nocov end
     resampling <- tabulate(sample(deindexedChars, charsToKeep, replace = FALSE),
                            nbins = length(startWeights))
-    errors <- vapply(eachChar, function (i) 
-      mpl_set_charac_weight(i, resampling[i], morphyObj), integer(1))
-    if (any(errors)) { #nocov start
-      stop ("Error resampling morphy object: ", 
-            mpl_translate_error(unique(errors[errors < 0L])))
+    if (useMorphy) {
+      errors <- vapply(eachChar, function (i)
+        mpl_set_charac_weight(i, resampling[i], initializedData), integer(1))
+      if (any(errors)) { #nocov start
+        stop("Error resampling morphy object: ",
+             mpl_translate_error(unique(errors[errors < 0L])))
+      }
+      if (mpl_apply_tipdata(initializedData) -> error) {
+        stop("Error applying tip data: ", mpl_translate_error(error))
+      } #nocov end
+      searchData <- initializedData
+    } else {
+      # Native: local copy with resampled weights
+      searchData <- initializedData
+      searchData[["weight"]] <- as.integer(resampling)
     }
-    if (mpl_apply_tipdata(morphyObj) -> error) {
-      stop("Error applying tip data: ", mpl_translate_error(error))
-    } #nocov end
-    res <- EdgeListSearch(edgeList[1:2], morphyObj, EdgeSwapper = EdgeSwapper,
+    res <- EdgeListSearch(edgeList[1:2], searchData,
+                          TreeScorer = TreeScorer,
+                          EdgeSwapper = EdgeSwapper,
                           maxIter = searchIter, maxHits = searchHits,
                           verbosity = verbosity - 1L, ...)
     res[1:2]

R/Morphy.R

@@ -253,7 +253,7 @@
 #'   c(0, 1, 1, 1, 0, 0,
 #'     1, 1, 1, 0, 0, 0), ncol = 2,
 #'   dimnames = list(letters[1:6], NULL)))
-#' Morphy(characters, constraint = constraint, verbosity = 0)
+#' suppressWarnings(Morphy(characters, constraint = constraint, verbosity = 0))
 #' 
 #' @template MRS
 #' 
@@ -578,8 +578,8 @@ Morphy <- function(dataset, tree,
 
   # Initialize constraints
   if (constrained) {
-    morphyConstr <- PhyDat2Morphy(constraint)
-    on.exit(morphyConstr <- UnloadMorphy(morphyConstr), add = TRUE)
+    morphyConstr <- suppressWarnings(PhyDat2Morphy(constraint))
+    on.exit(suppressWarnings(morphyConstr <- UnloadMorphy(morphyConstr)), add = TRUE)
     constraintWeight <- attr(constraint, "weight")
     if (any(constraintWeight > 1)) {
       cli_alert_warning("Some constraints are exact duplicates.")
@@ -651,8 +651,8 @@ Morphy <- function(dataset, tree,
   if ((!iw && !profile) || # Required for equal weights search
       (isTRUE(ratchEW) && ratchIter > 0) # For EW ratchet searches
   ) {
-    morphyObj <- PhyDat2Morphy(dataset)
-    on.exit(morphyObj <- UnloadMorphy(morphyObj), add = TRUE)
+    morphyObj <- suppressWarnings(PhyDat2Morphy(dataset))
+    on.exit(suppressWarnings(morphyObj <- UnloadMorphy(morphyObj)), add = TRUE)
   }
 
   if (iw || profile) {
@@ -662,8 +662,9 @@ Morphy <- function(dataset, tree,
     startWeights <- at[["weight"]]
     minLength <- MinimumLength(dataset, compress = TRUE)
     morphyObjects <- lapply(characters, SingleCharMorphy)
-    on.exit(morphyObjects <- vapply(morphyObjects, UnloadMorphy, integer(1)),
-            add = TRUE)
+    on.exit(suppressWarnings(
+      morphyObjects <- vapply(morphyObjects, UnloadMorphy, integer(1))),
+      add = TRUE)
 
     nLevel <- length(at[["level"]])
     nChar <- at[["nr"]]

R/NativeSearch.R

@@ -0,0 +1,134 @@
+#' Native C++ scoring for custom search functions
+#'
+#' These functions provide a native C++ scoring interface for use with
+#' [`TreeSearch()`], [`Ratchet()`], and [`Jackknife()`], replacing the
+#' MorphyLib-based defaults (`PhyDat2Morphy`, `UnloadMorphy`,
+#' `MorphyLength`, `MorphyBootstrap`).
+#'
+#' @name NativeSearch
+#' @family custom search functions
+NULL
+
+#' @describeIn NativeSearch Prepare a phyDat dataset for native C++ scoring.
+#' Replaces [`PhyDat2Morphy()`] as the `InitializeData` function.
+#'
+#' @param dataset A phyDat object.
+#' @param concavity Concavity constant for implied weighting, `Inf` for
+#'   equal weights, or `"profile"` for profile parsimony.
+#' @return `PrepareNativeData()` returns a list containing pre-extracted data
+#'   matrices suitable for repeated scoring with [`NativeLength()`].
+#' @export
+PrepareNativeData <- function(dataset, concavity = Inf) {
+  if (!inherits(dataset, "phyDat")) {
+    stop("`dataset` must be a phyDat object.")
+  }
+
+  useProfile <- identical(concavity, "profile")
+  if (useProfile) {
+    dataset <- PrepareDataProfile(dataset)
+    infoAmounts <- attr(dataset, "info.amounts")
+  } else {
+    infoAmounts <- NULL
+  }
+
+  iw <- !useProfile && is.finite(concavity)
+  if (iw) {
+    if (concavity <= 0) {
+      stop("`concavity` must be positive (or Inf for equal weights, ",
+           "or \"profile\" for profile parsimony).")
+    }
+    if (!("min.length" %in% names(attributes(dataset)))) {
+      dataset <- PrepareDataIW(dataset)
+    }
+    minSteps <- as.integer(attr(dataset, "min.length"))
+  } else {
+    minSteps <- integer(0)
+  }
+
+  at <- attributes(dataset)
+  weight <- at[["weight"]]
+
+  list(
+    contrast = at[["contrast"]],
+    tip_data = matrix(unlist(dataset, use.names = FALSE),
+                      nrow = length(dataset), byrow = TRUE),
+    weight = as.integer(weight),
+    levels = at[["levels"]],
+    min_steps = minSteps,
+    concavity = if (iw) as.double(concavity) else -1.0,
+    info_amounts = infoAmounts,
+    original_weight = as.integer(weight),
+    index = at[["index"]]
+  )
+}
+
+#' @describeIn NativeSearch No-op cleanup (native data has no external
+#'   resources). Replaces [`UnloadMorphy()`] as the `CleanUpData` function.
+#' @param dataset Prepared dataset returned by `PrepareNativeData()`.
+#' @return `CleanNativeData()` returns `dataset` invisibly.
+#' @export
+CleanNativeData <- function(dataset) {
+  invisible(dataset)
+}
+
+#' @describeIn NativeSearch Score a tree using the native C++ Fitch engine.
+#' Replaces [`MorphyLength()`] as the `TreeScorer` function.
+#'
+#' @param parent Integer vector of parent node indices for each edge.
+#' @param child Integer vector of child node indices for each edge.
+#' @param dataset Prepared dataset from [`PrepareNativeData()`].
+#' @param \dots Ignored (present for interface compatibility).
+#' @return `NativeLength()` returns a double: the parsimony score of the tree.
+#' @export
+NativeLength <- function(parent, child, dataset, ...) {
+  ts_fitch_score(
+    cbind(parent, child),
+    dataset[["contrast"]],
+    dataset[["tip_data"]],
+    dataset[["weight"]],
+    dataset[["levels"]],
+    dataset[["min_steps"]],
+    dataset[["concavity"]],
+    dataset[["info_amounts"]]
+  )
+}
+
+#' @describeIn NativeSearch Bootstrap resampling for use with [`Ratchet()`].
+#' Replaces [`MorphyBootstrap()`] as the `Bootstrapper` function.
+#'
+#' @param edgeList A list containing parent and child integer vectors
+#'   (and optionally a score and hit count).
+#' @param dataset Prepared dataset from [`PrepareNativeData()`].
+#' @param EdgeSwapper A function to rearrange edges, such as
+#'   [`NNISwap()`] or [`TBRSwap()`].
+#' @param maxIter Maximum rearrangement iterations.
+#' @param maxHits Maximum hits of best score before stopping.
+#' @param verbosity Verbosity level.
+#' @param stopAtPeak,stopAtPlateau Passed to [`EdgeListSearch()`].
+#' @param \dots Further parameters passed to the `TreeScorer`.
+#' @return `NativeBootstrap()` returns an edge list (list of parent and
+#'   child vectors).
+#' @export
+NativeBootstrap <- function(edgeList, dataset, EdgeSwapper = NNISwap,
+                            maxIter, maxHits, verbosity = 1L,
+                            stopAtPeak = FALSE, stopAtPlateau = 0L, ...) {
+  startWeights <- dataset[["original_weight"]]
+  eachChar <- seq_along(startWeights)
+  deindexedChars <- rep.int(eachChar, startWeights)
+  resampling <- tabulate(sample(deindexedChars, replace = TRUE),
+                         length(startWeights))
+
+  # R's copy-on-modify: local copy, caller's dataset unchanged
+  dataset[["weight"]] <- as.integer(resampling)
+
+  res <- EdgeListSearch(edgeList[1:2], dataset,
+                        TreeScorer = NativeLength,
+                        EdgeSwapper = EdgeSwapper,
+                        maxIter = maxIter, maxHits = maxHits,
+                        stopAtPeak = stopAtPeak,
+                        stopAtPlateau = stopAtPlateau,
+                        verbosity = verbosity - 1L, ...)
+
+  # Return:
+  res[1:2]
+}

R/Ratchet.R

@@ -69,11 +69,11 @@
 #' @family custom search functions
 #' @importFrom TreeTools RenumberEdges RenumberTips
 #' @export
-Ratchet <- function(tree, dataset, 
-                    InitializeData = PhyDat2Morphy,
-                    CleanUpData    = UnloadMorphy,
-                    TreeScorer     = MorphyLength,
-                    Bootstrapper   = MorphyBootstrap,
+Ratchet <- function(tree, dataset, concavity = Inf,
+                    InitializeData = PrepareNativeData,
+                    CleanUpData    = CleanNativeData,
+                    TreeScorer     = NativeLength,
+                    Bootstrapper   = NativeBootstrap,
                     swappers = list(TBRSwap, SPRSwap, NNISwap),
                     BootstrapSwapper = if (is.list(swappers))
                      swappers[[length(swappers)]] else swappers,
@@ -93,7 +93,11 @@ Ratchet <- function(tree, dataset,
   edgeList <- tree[["edge"]]
   edgeList <- RenumberEdges(edgeList[, 1], edgeList[, 2])
 
-  initializedData <- InitializeData(dataset)
+  if (identical(InitializeData, PrepareNativeData)) {
+    initializedData <- PrepareNativeData(dataset, concavity = concavity)
+  } else {
+    initializedData <- InitializeData(dataset)
+  }
   on.exit(initializedData <- CleanUpData(initializedData))
 
   bestScore <- TreeScorer(edgeList[[1]], edgeList[[2]], initializedData, ...)