X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Froot.R;h=ac5ec27b5b45d27470b8ca322adc787e86dc6193;hb=d7105238f98acefc5fc0ad6278a515e682453b8b;hp=e7c916b55e61679f5088bf892c36918d08b4e888;hpb=1d0651b1374592d87400614a03b34b4e0cc63aae;p=ape.git

diff --git a/R/root.R b/R/root.R
index e7c916b..ac5ec27 100644
--- a/R/root.R
+++ b/R/root.R
@@ -1,29 +1,28 @@
-## root.R (2008-02-12)
+## root.R (2010-02-11)
 
 ##   Root of Phylogenetic Trees
 
-## Copyright 2004-2008 Emmanuel Paradis
+## Copyright 2004-2010 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
 
 is.rooted <- function(phy)
 {
-    if (!"phylo" %in% class(phy))
-      stop('object "phy" is not of class "phylo"')
-    if (!is.null(phy$root.edge)) return(TRUE)
+    if (!inherits(phy, "phylo"))
+        stop('object "phy" is not of class "phylo"')
+    if (!is.null(phy$root.edge)) TRUE
     else
-      if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2)
-        return(FALSE)
-      else return(TRUE)
+        if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2)
+            FALSE else TRUE
 }
 
 unroot <- function(phy)
 {
-    if (class(phy) != "phylo")
-      stop('object "phy" is not of class "phylo"')
+    if (!inherits(phy, "phylo"))
+        stop('object "phy" is not of class "phylo"')
     if (dim(phy$edge)[1] < 3)
-      stop("cannot unroot a tree with two edges.")
+        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
@@ -36,11 +35,11 @@ unroot <- function(phy)
     ## nodes by adding 1, except the root (this remains the
     ## origin of the tree).
     nb.tip <- length(phy$tip.label)
-    ROOT <- nb.tip + 1
+    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 + 1) {
+    if (phy$edge[EDGEROOT[1], 2] == ROOT + 1L) {
         j <- EDGEROOT[2]
         i <- EDGEROOT[1]
     } else {
@@ -51,25 +50,29 @@ unroot <- function(phy)
     ## cladewise order.
     phy$edge <- phy$edge[-i, ]
     nodes <- phy$edge > ROOT # renumber all nodes except the root
-    phy$edge[nodes] <- phy$edge[nodes] - 1
+    phy$edge[nodes] <- phy$edge[nodes] - 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$Nnode <- phy$Nnode - 1
+    phy$Nnode <- phy$Nnode - 1L
     if (!is.null(phy$node.label))
       phy$node.label <- phy$node.label[-2]
     phy
 }
 
-root <- function(phy, outgroup, node = NULL, resolve.root = FALSE)
+root <- function(phy, outgroup, node = NULL,
+                 resolve.root = FALSE, interactive = FALSE)
 {
-    if (class(phy) != "phylo")
-      stop('object "phy" is not of class "phylo"')
-    ord <- attr(phy, "order")
-    if (!is.null(ord) && ord == "pruningwise") phy <- reorder(phy)
+    if (!inherits(phy, "phylo"))
+        stop('object "phy" is not of class "phylo"')
+    phy <- reorder(phy)
     n <- length(phy$tip.label)
-    ROOT <- n + 1
+    ROOT <- n + 1L
+    if (interactive) {
+        node <- identify(phy)$nodes
+        cat("You have set resolve.root =", resolve.root, "\n")
+    }
     if (!is.null(node)) {
         if (node <= n)
             stop("incorrect node#: should be greater than the number of taxa")
@@ -88,7 +91,6 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE)
         ## First check that the outgroup is monophyletic--
         ## unless there's only one tip specified of course
         if (length(outgroup) > 1) {
-            msg <- "the specified outgroup is not monophyletic"
             seq.nod <- .Call("seq_root2tip", phy$edge, n,
                              phy$Nnode, PACKAGE = "ape")
             sn <- seq.nod[outgroup]
@@ -108,19 +110,50 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE)
             ## (below is slightly faster than calling "bipartition")
             desc <- which(unlist(lapply(seq.nod,
                                         function(x) any(x %in% newroot))))
-            if (length(outgroup) != length(desc)) stop(msg)
-            ## both vectors below are already sorted:
-            if (!all(outgroup == desc)) stop(msg)
+            msg <- "the specified outgroup is not monophyletic"
+            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)
+            }
         } else newroot <- phy$edge[which(phy$edge[, 2] == outgroup), 1]
     }
-    if (newroot == ROOT) return(phy)
+    N <- Nedge(phy)
+    oldNnode <- phy$Nnode
+    if (newroot == ROOT) {
+        if (resolve.root) {
+            snw <- which(phy$edge[, 1] == newroot)
+            if (length(snw) > 2) {
+                a <- snw[1]:(snw[2] - 1)
+                b <- snw[2]:N
+                newnod <- oldNnode + n + 1
+                phy$edge[snw[-1], 1] <- newnod
+                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$Nnode <- phy$Nnode + 1L
+                ## node renumbering (see comments below)
+                newNb <- integer(n + oldNnode)
+                newNb[newroot] <- n + 1L
+                sndcol <- phy$edge[, 2] > n
+                phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
+                    (n + 2):(n + phy$Nnode)
+                phy$edge[, 1] <- newNb[phy$edge[, 1]]
+            }
+        }
+        return(phy)
+    }
 
     phy$root.edge <- NULL # just in case...
     Nclade <- tabulate(phy$edge[, 1])[ROOT] # degree of the root node
     ## if only 2 edges connect to the root, we have to fuse them:
     fuseRoot <- Nclade == 2
 
-    N <- Nedge(phy)
     start <- which(phy$edge[, 1] == ROOT)
     end <- c(start[-1] - 1, N)
     o <- integer(N)
@@ -218,7 +251,6 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE)
     ne <- length(s)
     o[NEXT:(NEXT + ne - 1L)] <- s
 
-    oldNnode <- phy$Nnode
     if (fuseRoot) {
         phy$Nnode <- oldNnode - 1
         N <- N - 1L
@@ -270,14 +302,19 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE)
     phy$edge[, 1] <- newNb[phy$edge[, 1]]
 
     if (!is.null(phy$node.label)) {
+        #browser()
         newNb <- newNb[-(1:n)]
         if (fuseRoot) {
             newNb <- newNb[-1]
             phy$node.label <- phy$node.label[-1]
         }
         phy$node.label <- phy$node.label[order(newNb)]
-        if (resolve.root)
-             phy$node.label <- c(phy$node.label[1], NA, phy$node.label[-1])
+        if (resolve.root) {
+            phy$node.label[is.na(phy$node.label)] <- phy$node.label[1]
+            phy$node.label[1] <- NA
+            ##phy$node.label <- c(phy$node.label[1], NA, phy$node.label[-1])
+            ##phy$node.label <- c("NA", phy$node.label)
+        }
     }
     phy
 }