From 639ef27e02f423fc53d1c08d2a84f06dd1907d2e Mon Sep 17 00:00:00 2001
From: "Martin R. Smith" <1695515+ms609@users.noreply.github.com>
Date: Tue, 10 Dec 2024 15:01:55 +0000
Subject: [PATCH] Improve SPR root behaviour (#165)

---
 DESCRIPTION               |  2 +-
 NEWS.md                   |  3 ++-
 R/SPR.R                   | 30 +++++++++++++++++++-----------
 tests/testthat/test-SPR.R | 18 ++++++++++++++++++
 4 files changed, 40 insertions(+), 13 deletions(-)
 create mode 100644 tests/testthat/test-SPR.R

diff --git a/DESCRIPTION b/DESCRIPTION
index 3c817547c..ec63cccb2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: TreeSearch
 Title: Phylogenetic Analysis with Discrete Character Data
-Version: 1.5.1.9000
+Version: 1.5.1.9001
 Authors@R: c(
     person(
       "Martin R.", 'Smith',
diff --git a/NEWS.md b/NEWS.md
index ac6f3f200..62ca4d8eb 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,7 @@
-# TreeSearch 1.5.1.9000 (2024-10)
+# TreeSearch 1.5.1.9001 (2024-12)
 
 - Fix display of state labels in app
+- Fix SPR behaviour when move is close to root
 
 
 # TreeSearch 1.5.1 (2024-05-23)
diff --git a/R/SPR.R b/R/SPR.R
index c409a8e10..c9b481c85 100644
--- a/R/SPR.R
+++ b/R/SPR.R
@@ -43,11 +43,11 @@ SPRWarning <- function (parent, child, error) {
   rightSide <- DescendantEdges(edge = 1, parent, child, nEdge = nEdge)
   nEdgeRight <- sum(rightSide)
   if (nEdgeRight == 1) {
-    notDuplicateRoot[2] <- FALSE
+    notDuplicateRoot[[2]] <- FALSE
   } else if (nEdgeRight == 3) {
-    notDuplicateRoot[4] <- FALSE
+    notDuplicateRoot[[4]] <- FALSE
   } else {
-    notDuplicateRoot[1] <- FALSE
+    notDuplicateRoot[[1]] <- FALSE
   }
   notDuplicateRoot
 }
@@ -97,12 +97,13 @@ SPR <- function(tree, edgeToBreak = NULL, mergeEdge = NULL) {
   if (!is.null(edgeToBreak) && edgeToBreak == -1) {
     child <- edge[, 2]
     nEdge <- length(parent)
-    stop("Negative edgeToBreak not yet supported; on TODO list for next release")
+    stop("Negative edgeToBreak not yet supported; please request on GitHub")
     notDuplicateRoot <- .NonDuplicateRoot(parent, child, nEdge)
     # Return:
     unique(unlist(lapply(which(notDuplicateRoot), AllSPR,
-      parent=parent, child=child, nEdge=nEdge, notDuplicateRoot=notDuplicateRoot),
-      recursive=FALSE)) # TODO the fact that we need to use `unique` indicates that 
+      parent = parent, child = child, nEdge = nEdge, 
+      notDuplicateRoot = notDuplicateRoot),
+      recursive = FALSE)) # TODO the fact that we need to use `unique` indicates that 
                          #      we're being inefficient here.
   } else {
     newEdge <- SPRSwap(parent, edge[, 2], edgeToBreak = edgeToBreak,
@@ -180,7 +181,7 @@ SPRSwap <- function (parent, child, nEdge = length(parent), nNode = nEdge / 2L,
   
   if (is.null(edgeToBreak)) {
     # Pick an edge at random
-    edgeToBreak <- SampleOne(which(notDuplicateRoot), len=nEdge - 1L)
+    edgeToBreak <- SampleOne(which(notDuplicateRoot), len = nEdge - 1L)
   } else if (edgeToBreak > nEdge) {
     return(SPRWarning(parent, child, "edgeToBreak > nEdge"))
   } else if (edgeToBreak < 1) {
@@ -199,10 +200,17 @@ SPRSwap <- function (parent, child, nEdge = length(parent), nNode = nEdge / 2L,
   brokenEdgeSister <- parent == brokenEdge.parentNode & !brokenEdge
   brokenEdgeDaughters <- parent == brokenEdge.childNode
   nearBrokenEdge <- brokenEdge | brokenEdgeSister | brokenEdgeParent | brokenEdgeDaughters
-  if (breakingRootEdge <- !any(brokenEdgeParent)) { 
+  breakingRootEdge <- !any(brokenEdgeParent)
+  if (breakingRootEdge) {
+    if (edgeToBreak != 1 && all(edgesCutAdrift[-1])) {
+      return(SPRWarning(parent, child, "No rearrangement possible with this root position."))
+    }
+    
     # Edge to break is the Root Node.
+    # These daughters are going to have the root as a parent.
     brokenRootDaughters <- parent == child[brokenEdgeSister]
-    nearBrokenEdge <- nearBrokenEdge | brokenRootDaughters
+    # Why did I do this?  Breaks SPR(BalancedTree(4), 1)
+    # nearBrokenEdge <- nearBrokenEdge | brokenRootDaughters
   }
   
   if (!is.null(mergeEdge)) { # Quick sanity checks
@@ -233,9 +241,9 @@ SPRSwap <- function (parent, child, nEdge = length(parent), nNode = nEdge / 2L,
   if (breakingRootEdge) {
     parent[brokenRootDaughters] <- brokenEdge.parentNode
     spareNode <- child[brokenEdgeSister]
-    child [brokenEdgeSister] <- child[mergeEdge]
+    child[brokenEdgeSister] <- child[[mergeEdge]]
     parent[brokenEdge | brokenEdgeSister] <- spareNode
-    child[mergeEdge] <- spareNode
+    child[[mergeEdge]] <- spareNode
   } else {
     parent[brokenEdgeSister] <- parent[brokenEdgeParent]
     parent[brokenEdgeParent] <- parent[[mergeEdge]]
diff --git a/tests/testthat/test-SPR.R b/tests/testthat/test-SPR.R
new file mode 100644
index 000000000..35eff2823
--- /dev/null
+++ b/tests/testthat/test-SPR.R
@@ -0,0 +1,18 @@
+library("TreeTools")
+
+test_that("SPR handles root rearrangement", {
+  # Soft tests.
+  expect_true(inherits(SPR(BalancedTree(4), edgeToBreak = 1), "phylo"))
+  expect_true(inherits(SPR(BalancedTree(5), edgeToBreak = 1), "phylo"))
+  
+  expect_equal(sum(sapply(c(
+    # hard-coded 1 %in% 2 due to 
+    # https://github.com/r-lib/testthat/issues/1661
+    SPR(PectinateTree(4), edgeToBreak = 1, mergeEdge = 3),
+    SPR(PectinateTree(4), edgeToBreak = 1, mergeEdge = 4),
+    SPR(PectinateTree(4), edgeToBreak = 1, mergeEdge = 5),
+    SPR(PectinateTree(4), edgeToBreak = 1, mergeEdge = 6)),
+    all.equal, SPR(PectinateTree(4), edgeToBreak = 1))), 1)
+  expect_warning(SPR(PectinateTree(4), edgeToBreak = 2),
+                 "No rearrangement possible with this root position")
+})