Score MorphyLength via native kernel for inapplicable data (Stage 1 of MorphyLib removal)

DESCRIPTION

@@ -10,14 +10,18 @@ Authors@R: c(
     ),
     person(
       given = "Martin", family = "Brazeau",
-      role = c("cph"),
-      comment = c(ORCID = "0000-0002-0650-1282", "Morphy library")
+      role = c("ctb"),
+      comment = c(ORCID = "0000-0002-0650-1282",
+                  "Brazeau-Guillerme-Smith inapplicable-data algorithm")
+    ),
+    person(
+      given = "Thomas", family = "Guillerme",
+      role = c("ctb"),
+      comment = c(ORCID = "0000-0003-4325-1275",
+                  "Brazeau-Guillerme-Smith inapplicable-data algorithm")
     )
   )
 License: GPL (>= 3)
-Copyright: Incorporates C/C++ code from
-   Morphy Phylogenetic Library by Martin Brazeau 
-   <https://github.com/mbrazeau/MorphyLib> (GPL3)
 Description: Reconstruct phylogenetic trees from discrete data using a
   high-performance C++ search engine with multi-replicate driven search.
   Supports equal weights, implied weights (Goloboff, 1993)
@@ -127,8 +131,8 @@ Collate:
   'data_manipulation.R'
   'fractional-weights.R'
   'length_range.R'
-  'mpl_morphy_objects.R'
-  'mpl_morphyex.R'
+  'PrepareData.R'
+  'morphy-deprecated.R'
   'pp_info_extra_step.r'
   'tree_length.R'
   'tree_rearrangement.R'
@@ -144,3 +148,4 @@ Encoding: UTF-8
 Language: en-GB
 VignetteBuilder: knitr
 Roxygen: list(markdown = TRUE)
+RoxygenNote: 7.3.3

NAMESPACE

@@ -21,11 +21,10 @@ S3method(TreeLength,phylo)
 S3method(print,CharacterHierarchy)
 S3method(print,ScoreSpectrum)
 S3method(print,SearchControl)
-S3method(summary,morphyPtr)
 export(.NonDuplicateRoot)
 export(.UniqueExceptHits)
 export(AdditionTree)
-export(C_MorphyLength)
+export(BootstrapTree)
 export(Carter1)
 export(CharacterHierarchy)
 export(CharacterLength)
@@ -38,14 +37,13 @@ export(Consistency)
 export(DoNothing)
 export(EasyTrees)
 export(EasyTreesy)
+export(EdgeListScore)
 export(EdgeListSearch)
 export(EmptyPhyDat)
 export(Evaluate)
 export(ExpectedLength)
 export(Fitch)
 export(FitchSteps)
-export(GapHandler)
-export(GetMorphyLength)
 export(IWScore)
 export(IWTreeSearch)
 export(JackLabels)
@@ -61,12 +59,9 @@ export(MaximizeParsimony)
 export(MaximizeParsimony2)
 export(MaximumLength)
 export(MinimumLength)
-export(Morphy)
 export(MorphyBootstrap)
-export(MorphyErrorCheck)
 export(MorphyLength)
 export(MorphyTreeLength)
-export(MorphyWeights)
 export(MostContradictedFreq)
 export(MultiRatchet)
 export(MutualClusteringConcordance)
@@ -78,6 +73,7 @@ export(PhyDat2Morphy)
 export(PhylogeneticConcordance)
 export(PlotCharacter)
 export(PolEscapa)
+export(PrepareData)
 export(PrepareDataIW)
 export(PrepareDataProfile)
 export(PresCont)
@@ -87,10 +83,12 @@ export(QCol)
 export(QuartetConcordance)
 export(QuartetResolution)
 export(RandomMorphyTree)
+export(RandomPostorderTree)
 export(RandomTreeScore)
 export(Ratchet)
 export(RatchetConsensus)
 export(RearrangeEdges)
+export(ReleaseData)
 export(Resample)
 export(RootedNNI)
 export(RootedNNISwap)
@@ -104,8 +102,8 @@ export(SPRSwap)
 export(SPRWarning)
 export(ScoreSpectrum)
 export(SearchControl)
-export(SetMorphyWeights)
 export(SharedPhylogeneticConcordance)
+export(SingleCharData)
 export(SingleCharMorphy)
 export(StepInformation)
 export(StopUnlessBifurcating)
@@ -118,6 +116,7 @@ export(TBRSwap)
 export(TBRWarning)
 export(TaxonInfluence)
 export(TreeLength)
+export(TreeScore)
 export(TreeSearch)
 export(UnloadMorphy)
 export(WhenFirstHit)
@@ -130,31 +129,9 @@ export(hierarchy_chars)
 export(hierarchy_controlling)
 export(hierarchy_from_names)
 export(hierarchy_to_blocks)
+export(is.ParsimonyData)
 export(is.morphyPtr)
 export(mc_fitch_scores)
-export(mpl_apply_tipdata)
-export(mpl_attach_rawdata)
-export(mpl_attach_symbols)
-export(mpl_delete_Morphy)
-export(mpl_first_down_recon)
-export(mpl_first_up_recon)
-export(mpl_get_charac_weight)
-export(mpl_get_gaphandl)
-export(mpl_get_num_charac)
-export(mpl_get_num_internal_nodes)
-export(mpl_get_numtaxa)
-export(mpl_get_symbols)
-export(mpl_init_Morphy)
-export(mpl_new_Morphy)
-export(mpl_second_down_recon)
-export(mpl_second_up_recon)
-export(mpl_set_charac_weight)
-export(mpl_set_gaphandl)
-export(mpl_set_num_internal_nodes)
-export(mpl_set_parsim_t)
-export(mpl_translate_error)
-export(mpl_update_lower_root)
-export(mpl_update_tip)
 export(non_hierarchy_weights)
 export(recode_hierarchy)
 export(validate_hierarchy)
@@ -167,19 +144,14 @@ importFrom(TreeDist,Entropy)
 importFrom(TreeDist,MutualClusteringInfo)
 importFrom(TreeDist,SharedPhylogeneticInfo)
 importFrom(TreeDist,entropy_int)
-importFrom(TreeTools,AddUnconstrained)
-importFrom(TreeTools,CharacterInformation)
 importFrom(TreeTools,CladeSizes)
 importFrom(TreeTools,CladisticInfo)
 importFrom(TreeTools,CompatibleSplits)
 importFrom(TreeTools,Consensus)
-importFrom(TreeTools,ConstrainedNJ)
 importFrom(TreeTools,DescendantEdges)
 importFrom(TreeTools,DescendantTips)
 importFrom(TreeTools,DoubleFactorial)
-importFrom(TreeTools,DropTip)
 importFrom(TreeTools,EdgeAncestry)
-importFrom(TreeTools,ImposeConstraint)
 importFrom(TreeTools,KeepTip)
 importFrom(TreeTools,LnUnrooted)
 importFrom(TreeTools,Log2DoubleFactorial)
@@ -198,7 +170,6 @@ importFrom(TreeTools,NexusTokens)
 importFrom(TreeTools,PaintTree)
 importFrom(TreeTools,PectinateTree)
 importFrom(TreeTools,PhyDatToMatrix)
-importFrom(TreeTools,PhyToString)
 importFrom(TreeTools,PolarizeSplits)
 importFrom(TreeTools,Postorder)
 importFrom(TreeTools,PostorderOrder)
@@ -229,13 +200,11 @@ importFrom(ape,read.nexus)
 importFrom(ape,root)
 importFrom(ape,write.nexus)
 importFrom(cli,cli_alert)
-importFrom(cli,cli_alert_danger)
 importFrom(cli,cli_alert_info)
 importFrom(cli,cli_alert_success)
 importFrom(cli,cli_alert_warning)
 importFrom(cli,cli_h1)
 importFrom(cli,cli_progress_bar)
-importFrom(cli,cli_progress_done)
 importFrom(cli,cli_progress_update)
 importFrom(colorspace,hex)
 importFrom(colorspace,max_chroma)
@@ -254,7 +223,6 @@ importFrom(stats,cophenetic)
 importFrom(stats,cutree)
 importFrom(stats,dnorm)
 importFrom(stats,median)
-importFrom(stats,runif)
 importFrom(stats,sd)
 importFrom(stats,setNames)
 importFrom(stats,weighted.mean)

NEWS.md

@@ -1,5 +1,39 @@
 # To integrate into 2.0.0 notes
 
+- **MorphyLib removed.** The vendored Morphy Phylogenetic Library (C/C++) has
+  been dropped; all parsimony scoring now runs through the native C++ kernel,
+  which implements the Brazeau, Guillerme & Smith (2019) inapplicable-state
+  algorithm correctly — including ambiguous-with-inapplicable tokens such as
+  `{1-}`, which MorphyLib mis-scored.
+
+- **`concavity` argument for `TreeSearch()`, `Ratchet()` and `Jackknife()`.**
+  Implied-weights and profile searches with the custom-search functions no
+  longer need a hand-written scorer: pass `concavity = k` (a finite constant)
+  for implied weights, `concavity = "profile"` for profile parsimony, or the
+  default `Inf` for equal weights.  (Adapted from the parallel T-200 work,
+  PR #216.)
+
+- **Custom-search scoring layer renamed** to drop the now-meaningless "Morphy"
+  branding.  `PhyDat2Morphy()` → `PrepareData()`; `UnloadMorphy()` →
+  `ReleaseData()`; `is.morphyPtr()` → `is.ParsimonyData()`; `SingleCharMorphy()`
+  → `SingleCharData()`; `MorphyLength()` → `EdgeListScore()`;
+  `MorphyTreeLength()` → `TreeScore()`; `MorphyBootstrap()` → `BootstrapTree()`;
+  `RandomMorphyTree()` → `RandomPostorderTree()`.  The old names remain as
+  deprecated aliases and will be removed in a future release.  `PrepareData()`
+  returns a lightweight, garbage-collected `ParsimonyData` object; only the
+  default `"inapplicable"` gap treatment is supported (recode data for the
+  missing/extra-state treatments).
+
+- `Jackknife()` and `BootstrapTree()` (formerly `MorphyBootstrap()`) now
+  resample characters natively, scoring the resampled weights through the
+  native kernel rather than by mutating a MorphyLib object — fixing a case
+  where resampled weights could be silently ignored.
+
+- The low-level MorphyLib bindings (`mpl_*()`), together with the
+  `MorphyWeights()`, `SetMorphyWeights()`, `GapHandler()`, `MorphyErrorCheck()`,
+  `GetMorphyLength()` and `C_MorphyLength()` helpers and the
+  `summary.morphyPtr()` method, have been removed.
+
 - `WideSample()` now dispatches to the appropriate Max-Min diversity (MMDP)
   solver from the \pkg{MaxMin} package, choosing the tier automatically
   from `length(trees)`: the exact node-packing optimum (Sayyady & Fathi, 2016)

R/Bootstrap.R

@@ -1,51 +1,34 @@
 #' @inheritParams TreeSearch
 #' @inheritParams EdgeListSearch
-#' @inheritParams MorphyTreeLength
+#' @param dataset A `ParsimonyData` object from [`PrepareData()`].
 #' @param maxIter Numeric specifying maximum number of iterations to perform in
 #' tree search.
 #' @param maxHits Numeric specifying maximum number of hits to accomplish in
 #' tree search.
 #' @param stopAtScore stop search as soon as this score is hit or beaten.
 #' @param \dots further parameters to send to `TreeScorer()`
 #'
-#' @return `MorphyBootstrap()` returns a tree that is optimal under a random
+#' @return `BootstrapTree()` returns a tree that is optimal under a random
 #' sampling of the original characters.
-#' 
+#'
 #' @rdname Ratchet
 #' @export
-MorphyBootstrap <- function (edgeList, morphyObj, EdgeSwapper = NNISwap, 
-                             maxIter, maxHits, verbosity = 1L, 
-                             stopAtPeak = FALSE, stopAtPlateau=0L, ...) {
-  startWeights <- MorphyWeights(morphyObj)["exact", ]
+BootstrapTree <- 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))
-  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 (mpl_apply_tipdata(morphyObj) -> error) {
-    stop("Error applying tip data: ", mpl_translate_error(error))
-  } # nocov end
-
-  res <- EdgeListSearch(edgeList[1:2], morphyObj, EdgeSwapper = EdgeSwapper,
+  # R copy-on-modify: the caller's `dataset` is unchanged.
+  dataset[["weight"]] <- as.integer(resampling)
+
+  res <- EdgeListSearch(edgeList[1:2], dataset, EdgeSwapper = EdgeSwapper,
                         maxIter = maxIter, maxHits = maxHits,
                         stopAtPeak = stopAtPeak, stopAtPlateau = stopAtPlateau,
                         verbosity = verbosity - 1L, ...)
-  errors <- vapply(eachChar, function (i) 
-         mpl_set_charac_weight(i, startWeights[i], morphyObj), integer(1))
-  if (any(errors)) { # nocov start
-    stop ("Error resampling morphy object: ",
-          mpl_translate_error(unique(errors[errors < 0L])))
-  }
-  if (mpl_apply_tipdata(morphyObj) -> error) {
-    stop("Error applying tip data: ", mpl_translate_error(error))
-  } # nocov end
-
+
   # Return:
   res[1:2]
 }

R/CustomSearch.R

@@ -12,7 +12,7 @@
 #' @keywords internal
 #' @export
 EdgeListSearch <- function (edgeList, dataset,
-                            TreeScorer = MorphyLength,
+                            TreeScorer = EdgeListScore,
                             EdgeSwapper = RootedTBRSwap,
                             maxIter = 100, maxHits = 20, 
                             bestScore = NULL, stopAtScore = NULL, 
@@ -128,13 +128,19 @@ EdgeListSearch <- function (edgeList, dataset,
 #' Will be overridden if a passed function has an attribute `stopAtPlateau` set
 #' by `attr(FunctionName, "stopAtPlateau") <- TRUE`.
 #'
-#' @param InitializeData Function that sets up data object to prepare for tree search. 
-#'        The function will be passed the `dataset` parameter.
+#' @param concavity Determines the optimality criterion when the default native
+#' scorer is used: `Inf` (the default) for equal weights; a finite, positive
+#' number for implied weighting with that concavity constant
+#' \insertCite{Goloboff1993}{TreeSearch}; or `"profile"` for profile parsimony.
+#' Ignored if a custom `InitializeData` is supplied (prepare the data yourself).
+#' @param InitializeData Function that sets up data object to prepare for tree
+#'        search.  The default, [`PrepareData()`], is called with `concavity`;
+#'        a custom function will be passed only the `dataset` parameter.
 #'        Its return value will be passed to `TreeScorer()` and `CleanUpData()`.
-#' @param CleanUpData Function to destroy data object on function exit.
+#' @param CleanUpData Function to release the data object on function exit.
 #'        The function will be passed the value returned by `InitializeData()`.
 #' @param TreeScorer function to score a given tree.
-#'        The function will be passed three parameters, corresponding to the 
+#'        The function will be passed three parameters, corresponding to the
 #'        `parent` and `child` entries of a tree's edge list, and a dataset.
 #'
 #' @param verbosity Numeric specifying level of detail to display in console: 
@@ -166,15 +172,18 @@ EdgeListSearch <- function (edgeList, dataset,
 #'   TreeSearch(njtree, Lobo.phy, maxIter = 20, EdgeSwapper = NNISwap)
 #'   TreeSearch(njtree, Lobo.phy, maxIter = 20, EdgeSwapper = RootedSPRSwap)
 #'   TreeSearch(njtree, Lobo.phy, maxIter = 20, EdgeSwapper = TBRSwap)
+#'
+#'   # Implied weighting (concavity constant k = 10):
+#'   TreeSearch(njtree, Lobo.phy, concavity = 10, maxIter = 20)
 #' }
 #' @template MRS
 #' @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 = PrepareData,
+                        CleanUpData    = ReleaseData,
+                        TreeScorer     = EdgeListScore,
                         EdgeSwapper    = RootedTBRSwap,
                         maxIter = 100L, maxHits = 20L,
                         stopAtPeak = FALSE, stopAtPlateau = 0L,
@@ -186,7 +195,13 @@ TreeSearch <- function (tree, dataset,
   edgeList <- tree[["edge"]]
   edgeList <- RenumberEdges(edgeList[, 1], edgeList[, 2])
 
-  initializedData <- InitializeData(dataset)
+  # `concavity` is honoured only by the default native `PrepareData`; a custom
+  # `InitializeData` receives the dataset unchanged.
+  if (identical(InitializeData, PrepareData)) {
+    initializedData <- PrepareData(dataset, concavity = concavity)
+  } else {
+    initializedData <- InitializeData(dataset)
+  }
   on.exit(initializedData <- CleanUpData(initializedData))
 
   bestScore <- attr(tree, "score")