- if (class(phy) != "phylo")
- stop('object "phy" is not of class "phylo"')
- if (is.character(outgroup))
- outgroup <- which(phy$tip.label %in% outgroup)
- nb.tip <- length(phy$tip.label)
- if (length(outgroup) == nb.tip) return(phy)
-
- ## First check that the outgroup is monophyletic--
- ## unless there's only one tip specified of course
- ROOT <- nb.tip + 1
- if (length(outgroup) > 1) {
- msg <- "the specified outgroup is not monophyletic"
- ## If all tips in `tip' are not contiguous, then
- ## no need to go further:
- if (!all(diff(outgroup) == 1)) stop(msg)
- seq.nod <- .Call("seq_root2tip", phy$edge, nb.tip,
- 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
+ if (!inherits(phy, "phylo"))
+ stop('object "phy" is not of class "phylo"')
+ phy <- reorder(phy)
+ n <- length(phy$tip.label)
+ 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")
+ outgroup <- NULL
+ newroot <- node
+ } else {
+ if (is.numeric(outgroup)) {
+ if (any(outgroup > n))
+ stop("incorrect taxa#: should not be greater than the number of taxa")
+ outgroup <- sort(outgroup) # used below
+ }
+ if (is.character(outgroup))
+ outgroup <- which(phy$tip.label %in% outgroup)
+ if (length(outgroup) == n) return(phy)
+
+ ## 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"
+ 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]
+ }
+ 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
+
+ start <- which(phy$edge[, 1] == ROOT)
+ end <- c(start[-1] - 1, N)
+ o <- integer(N)
+ INV <- logical(N)
+
+ w <- which(phy$edge[, 2] == newroot)
+ nod <- phy$edge[w, 1]
+ i <- w
+ NEXT <- 1L
+
+ ## The loop below starts from the new root and goes up in the edge
+ ## matrix reversing the edges that need to be as well as well
+ ## inverting their order. The other edges must not be changed, so
+ ## their indices are stored in `stack'.
+ ## We then bind the other edges in a straightforward way.
+
+ if (nod != ROOT) {
+ ## it is important that the 3 next lines
+ ## are inside this "if" statement
+ o[NEXT] <- w
+ NEXT <- NEXT + 1L
+ INV[w] <- TRUE
+ i <- w - 1L
+ stack <- 0L