X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Froot.R;h=5655873a75398cdc16d8e84de39a53904ebc2e45;hb=f5c4abe6ac31486e821d82788bf66b5db2be51d1;hp=5691bb518f01e9b1105b7b886c5ca48952b02e61;hpb=3ece2ec76da287a8a86339827cc44e193fe16cdd;p=ape.git

diff --git a/R/root.R b/R/root.R
index 5691bb5..5655873 100644
--- a/R/root.R
+++ b/R/root.R
@@ -1,4 +1,4 @@
-## root.R (2011-04-29)
+## root.R (2011-08-05)
 
 ##   Root of Phylogenetic Trees
 
@@ -21,43 +21,73 @@ unroot <- function(phy)
 {
     if (!inherits(phy, "phylo"))
         stop('object "phy" is not of class "phylo"')
-    if (dim(phy$edge)[1] < 3)
+    N <- dim(phy$edge)[1]
+    if (N < 3)
         stop("cannot unroot a tree with less than three edges.")
+
     ## delete FIRST the root.edge (in case this is sufficient to
     ## unroot the tree, i.e. there is a multichotomy at the root)
     if (!is.null(phy$root.edge)) phy$root.edge <- NULL
     if (!is.rooted(phy)) return(phy)
-    ## We remove one of the edges coming from the root, and
-    ## eventually adding the branch length to the other one
-    ## also coming from the root.
-    ## In all cases, the node deleted is the 2nd one (numbered
-    ## nb.tip+2 in `edge'), so we simply need to renumber the
-    ## nodes by adding 1, except the root (this remains the
-    ## origin of the tree).
-    nb.tip <- length(phy$tip.label)
-    ROOT <- nb.tip + 1L
-    EDGEROOT <- which(phy$edge[, 1] == ROOT)
-    ## j: the target where to stick the edge
-    ## i: the edge to delete
-    if (phy$edge[EDGEROOT[1], 2] == ROOT + 1L) {
-        j <- EDGEROOT[2]
-        i <- EDGEROOT[1]
+
+    n <- length(phy$tip.label)
+    ROOT <- n + 1L
+
+### EDGEROOT[1]: the edge to delete
+### EDGEROOT[2]: the target where to stick the edge to delete
+
+### If the tree is in pruningwise order, then the last two edges
+### are those connected to the root; the node situated in
+### phy$edge[N - 2L, 1L] will be the new root...
+
+    ophy <- attr(phy, "order")
+    if (!is.null(ophy) && ophy == "pruningwise") {
+        NEWROOT <- phy$edge[N - 2L, 1L]
+        EDGEROOT <- c(N, N - 1L)
     } else {
-        j <- EDGEROOT[1]
-        i <- EDGEROOT[2]
+
+### ... otherwise, we remove one of the edges coming from
+### the root, and eventually adding the branch length to
+### the other one also coming from the root.
+### In all cases, the node deleted is the 2nd one (numbered
+### nb.tip+2 in 'edge'), so we simply need to renumber the
+### nodes by adding 1, except the root (this remains the
+### origin of the tree).
+
+        EDGEROOT <- which(phy$edge[, 1L] == ROOT)
+        NEWROOT <- ROOT + 1L
     }
-    ## This should work whether the tree is in pruningwise or
-    ## cladewise order.
-    phy$edge <- phy$edge[-i, ]
-    nodes <- phy$edge > ROOT # renumber all nodes except the root
-    phy$edge[nodes] <- phy$edge[nodes] - 1L
+
+    ## make sure EDGEROOT is ordered as described above:
+    if (phy$edge[EDGEROOT[1L], 2L] != NEWROOT)
+        EDGEROOT <- EDGEROOT[2:1]
+
+    phy$edge <- phy$edge[-EDGEROOT[1L], ]
+
+    s <- phy$edge == NEWROOT # renumber the new root
+    phy$edge[s] <- ROOT
+
+    s <- phy$edge > NEWROOT # renumber all nodes greater than the new root
+    phy$edge[s] <- phy$edge[s] - 1L
+
     if (!is.null(phy$edge.length)) {
-        phy$edge.length[j] <- phy$edge.length[j] + phy$edge.length[i]
-        phy$edge.length <- phy$edge.length[-i]
+        phy$edge.length[EDGEROOT[2L]] <-
+            phy$edge.length[EDGEROOT[2L]] + phy$edge.length[EDGEROOT[1L]]
+        phy$edge.length <- phy$edge.length[-EDGEROOT[1L]]
     }
+
     phy$Nnode <- phy$Nnode - 1L
-    if (!is.null(phy$node.label))
-      phy$node.label <- phy$node.label[-2]
+
+    if (!is.null(phy$node.label)) {
+        if (NEWROOT == n + 2L)
+            phy$node.label <- phy$node.label[-1]
+        else {
+            lbs <- phy$node.label
+            tmp <- lbs[NEWROOT - n]
+            lbs <- lbs[-c(1, NEWROOT)]
+            phy$node.label <- c(tmp, lbs)
+        }
+    }
     phy
 }
 
@@ -91,34 +121,21 @@ root <- function(phy, outgroup, node = NULL,
         ## First check that the outgroup is monophyletic--
         ## unless there's only one tip specified of course
         if (length(outgroup) > 1) {
-            seq.nod <- .Call("seq_root2tip", phy$edge, n,
-                             phy$Nnode, PACKAGE = "ape")
-            sn <- seq.nod[outgroup]
-            ## We go from the root to the tips: the sequence of nodes
-            ## is identical until the MRCA:
-            newroot <- ROOT
-            i <- 2 # we start at the 2nd position since the root
-                   # of the tree is a common ancestor to all tips
-            repeat {
-                x <- unique(unlist(lapply(sn, "[", i)))
-                if (length(x) != 1) break
-                newroot <- x
-                i <- i + 1
-            }
-            ## Check that all descendants of this node
-            ## are included in the outgroup.
-            ## (below is slightly faster than calling "bipartition")
-            desc <- which(unlist(lapply(seq.nod,
-                                        function(x) any(x %in% newroot))))
-            msg <- "the specified outgroup is not monophyletic"
+            pp <- prop.part(phy)
             ingroup <- (1:n)[-outgroup]
-            ## 'outgroup' and 'desc' are already sorted:
-            if (newroot != ROOT) {
-                if (!identical(outgroup, desc) && !identical(ingroup, desc))
-                    stop(msg)
-            } else { # otherwise check monophyly of the ingroup
-                if (!is.monophyletic(phy, ingroup)) stop(msg)
+            newroot <- 0L
+            for (i in 2:phy$Nnode) {
+                if (identical(pp[[i]], ingroup)) {
+                    newroot <- i + n
+                    break
+                }
+                if (identical(pp[[i]], outgroup)) {
+                    newroot <- phy$edge[which(phy$edge[, 2] == i + n), 1]
+                    break
+                }
             }
+            if (!newroot)
+                stop("the specified outgroup is not monophyletic")
         } else newroot <- phy$edge[which(phy$edge[, 2] == outgroup), 1]
     }
     N <- Nedge(phy)
@@ -134,8 +151,8 @@ root <- function(phy, outgroup, node = NULL,
                 phy$edge <- rbind(phy$edge[a, ], c(ROOT, newnod),
                                   phy$edge[b, ])
                 if (!is.null(phy$edge.length))
-                phy$edge.length <-
-                    c(phy$edge.length[a], 0, phy$edge.length[b])
+                    phy$edge.length <-
+                        c(phy$edge.length[a], 0, phy$edge.length[b])
                 phy$Nnode <- phy$Nnode + 1L
                 ## node renumbering (see comments below)
                 newNb <- integer(n + oldNnode)
@@ -297,8 +314,7 @@ root <- function(phy, outgroup, node = NULL,
     newNb[newroot] <- n + 1L
     sndcol <- phy$edge[, 2] > n
     ## executed from right to left, so newNb is modified before phy$edge:
-    phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
-        (n + 2):(n + phy$Nnode)
+    phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <- n + 2:phy$Nnode
     phy$edge[, 1] <- newNb[phy$edge[, 1]]
 
     if (!is.null(phy$node.label)) {