-## bind.tree.R (2010-02-12)
+## bind.tree.R (2010-03-15)
## Bind Trees
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-bind.tree <- function(x, y, where = "root", position = 0)
+`+.phylo` <- function(x, y)
{
- nb.tip <- length(x$tip.label)
- nb.node <- x$Nnode
- ROOT <- nb.tip + 1
- if (where == 0 || where == "root")
- where <- ROOT
- if (position < 0) position <- 0
- if (where > nb.tip + nb.node) stop("node number out of range for tree 'x'")
- nb.edge <- dim(x$edge)[1]
- yHasNoRootEdge <- is.null(y$root.edge)
- xHasNoRootEdge <- is.null(x$root.edge)
+ p <- if (is.null(x$root.edge)) 0 else x$root.edge
+ bind.tree(x, y, position = p)
+}
- ## check whether both trees have branch lengths:
- wbl <- TRUE
- noblx <- is.null(x$edge.length)
- nobly <- is.null(y$edge.length)
- if (noblx && nobly) wbl <- FALSE
- if (xor(noblx, nobly)) {
- if (nobly) x$edge.length <- NULL
- else y$edge.length <- NULL
- wbl <- FALSE
- warning("one tree has no branch lengths, they will be ignored")
- }
+bind.tree <- function(x, y, where = "root", position = 0, interactive = FALSE)
+{
+ nx <- length(x$tip.label)
+ mx <- x$Nnode
+ ROOTx <- nx + 1L
+ ny <- length(y$tip.label)
+ my <- y$Nnode
- ## To avoid problems with tips or nodes with indentical
- ## labels we substitute the one where `y' is grafted:
- if (where <= nb.tip) {
- Tip.Label.where <- x$tip.label[where]
- x$tip.label[where] <- "TheTipWhereToGraftY"
+ if (interactive) {
+ lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
+ if (lastPP$type != "phylogram" || lastPP$direction != "rightwards")
+ stop("you must plot tree 'x' as a 'rightward phylogram'")
+ cat("Click where you want to graft tree 'y'...\n")
+ xy <- locator(1)
+ d <- abs(xy$y - lastPP$yy)
+ d[lastPP$xx - xy$x < 0] <- Inf
+ where <- which.min(d)
+ position <- lastPP$xx[where] - xy$x
+ if (position < 0) position <- 0
+ cat("The following parameters are used:\n")
+ cat(" where =", where, " position =", position, "\n")
+ } else {
+ if (where == 0 || where == "root") where <- ROOTx
+ if (position < 0) position <- 0
+ if (where > nx + mx)
+ stop("argument 'where' out of range for tree 'x'")
}
- if (where > ROOT) {
- xHasNoNodeLabel <- TRUE
- if (is.null(x$node.label)) {
- x$node.label <- paste("NODE", 1:nb.node, sep = "")
- x$node.label[where - nb.tip] <- "TheNodeWhereToGraftY"
- } else {
- Node.Label.where <- x$node.label[where - nb.tip]
- x$node.label[where - nb.tip] <- "TheNodeWhereToGraftY"
- xHasNoNodeLabel <- FALSE
- }
+
+ ## check whether both trees have branch lengths:
+ switch(is.null(x$edge.length) + is.null(y$edge.length) + 1L,
+ wbl <- TRUE, {
+ x$edge.length <- y$edge.length <- NULL
+ wbl <- FALSE
+ warning("one tree has no branch lengths, they will be ignored")
+ },
+ wbl <- FALSE)
+
+ yHasNoRootEdge <- is.null(y$root.edge)
+ xHasNoRootEdge <- is.null(x$root.edge)
+
+ ## find the row of 'where' before renumbering
+ if (where == ROOTx) case <- 1 else {
+ i <- which(x$edge[, 2] == where)
+ case <- if (where <= nx) 2 else 3
}
+ ## case = 1 -> y is bound on the root of x
+ ## case = 2 -> y is bound on a tip of x
+ ## case = 3 -> y is bound on a node of x
- ## if we bind `y' under a node or tip of `y', we first
- ## adjust the edge lengths if needed
- if (position && wbl) {
- if (where == ROOT) {
- if (xHasNoRootEdge) stop("tree 'x' has no root edge")
- if (x$root.edge < position)
- stop("argument 'position' is larger than the root edge.")
- x$root.edge <- x$root.edge - position
+ ## check that 'position' is correct
+ if (position) {
+ if (!wbl)
+ stop("'position' is non-null but trees have no branch lengths")
+ if (case == 1) {
+ if (xHasNoRootEdge)
+ stop("tree 'x' has no root edge")
+ if (position > x$root.edge)
+ stop("'position' is larger than x's root edge")
} else {
- i <- which(x$edge[, 2] == where)
if (x$edge.length[i] < position)
- stop("argument 'position' is larger than the specified edge.")
- x$edge.length[i] <- x$edge.length[i] - position
+ stop("'position' is larger than the branch length")
}
- y$root.edge <- if (yHasNoRootEdge) position else y$root.edge + position
}
- if (is.null(y$root.edge) && where > nb.tip) y$root.edge <- 0
+ x <- reorder(x)
+ y <- reorder(y)
- X <- write.tree(x)
- Y <- write.tree(y)
- Y <- substr(Y, 1, nchar(Y) - 1)
+### because in all situations internal nodes need to be renumbered,
+### they are changed to negatives first, and nodes that eventually
+### will be added will be numbered sequentially
- if (where <= nb.tip) {
- if (position)
- X <- gsub("TheTipWhereToGraftY",
- paste("(", "TheTipWhereToGraftY", ",", Y, ")",
- sep = ""), X)
- else
- X <- gsub("TheTipWhereToGraftY", Y, X)
- }
- if (where == ROOT) {
- rmvx <- if (xHasNoRootEdge) "\\);$" else ";$"
- X <- gsub(rmvx, "", X)
- Y <- gsub("^\\(", "", Y)
- if (!xHasNoRootEdge) X <- paste("(", X, sep = "")
- X <- paste(X, ",", Y, ";", sep = "")
+ nodes <- x$edge > nx
+ x$edge[nodes] <- -(x$edge[nodes] - nx) # -1, ..., -mx
+ nodes <- y$edge > ny
+ y$edge[nodes] <- -(y$edge[nodes] - ny + mx) # -(mx+1), ..., -(mx+my)
+ ROOT <- -1L # may change later
+ next.node <- -(mx + my) - 1L
+
+ ## renumber now the tips in y:
+ new.nx <- if (where <= nx && !position) nx - 1L else nx
+ y$edge[!nodes] <- y$edge[!nodes] + new.nx
+
+ ## if 'y' as a root edge, use it:
+ if (!yHasNoRootEdge) {
+ y$edge <- rbind(c(0, y$edge[1]), y$edge)
+ ## ^ will be filled later
+ next.node <- next.node - 1L
+ if (wbl) y$edge.length <- c(y$root.edge, y$edge.length)
}
- if (where > ROOT) {
+
+ switch(case, { # case = 1
+ if (position) {
+ x$root.edge <- x$root.edge - position
+ x$edge <- rbind(c(next.node, x$edge[1]), x$edge)
+ ROOT <- next.node
+ if (wbl) x$edge.length <- c(position, x$edge.length)
+ }
+ if (yHasNoRootEdge) {
+ j <- which(y$edge[, 1] == y$edge[1])
+ y$edge[j, 1] <- ROOT
+ } else y$edge[1] <- ROOT
+ x$edge <- rbind(x$edge, y$edge)
+ if (wbl)
+ x$edge.length <- c(x$edge.length, y$edge.length)
+ }, { # case = 2
if (position) {
- ## find where is the node in `X':
- ## below 19 is: nchar("TheNodeWhereToGraftY") - 1
- for (i in 1:nchar(X)) {
- if ("TheNodeWhereToGraftY" == substr(X, i, i + 19))
- break
- i <- i + 1
+ x$edge[i, 2] <- next.node
+ x$edge <- rbind(x$edge[1:i, ], c(next.node, where), x$edge[-(1:i), ])
+ if (wbl) {
+ x$edge.length[i] <- x$edge.length[i] - position
+ x$edge.length <- c(x$edge.length[1:i], position, x$edge.length[-(1:i)])
}
- ## now go back to find the left matching parentheses
- n.paren <- 1
- i <- i - 2
- while (n.paren > 0) {
- if (substr(X, i, i) == ")") n.paren <- n.paren + 1
- if (substr(X, i, i) == "(") n.paren <- n.paren - 1
- i <- i - 1
+ i <- i + 1L
+ if (yHasNoRootEdge) {
+ j <- which(y$edge[, 1] == y$edge[1])
+ y$edge[j, 1] <- x$edge[i, 1]
+ } else y$edge[1] <- x$edge[i, 1]
+ } else {
+ if (yHasNoRootEdge) x$edge[i, 2] <- y$edge[1]
+ else {
+ ## the root edge of y is fused with the terminal edge of x
+ if (wbl) y$edge.length[1] <- y$edge.length[1] + x$edge.length[i]
+ y$edge[1] <- x$edge[i, 1]
+ ## delete i-th edge in x:
+ x$edge <- x$edge[-i, ]
+ i <- i - 1L
+ if (wbl) x$edge.length <- x$edge.length[-i]
+ }
+ x$tip.label <- x$tip.label[-where]
+ ## renumber the tips that need to:
+ ii <- which(x$edge[, 2] > where & x$edge[, 2] <= nx)
+ x$edge[ii, 2] <- x$edge[ii, 2] - 1L
+ }
+ x$edge <- rbind(x$edge[1:i, ], y$edge, x$edge[-(1:i), ])
+ if (wbl)
+ x$edge.length <- c(x$edge.length[1:i], y$edge.length, x$edge.length[-(1:i)])
+ }, { # case = 3
+ if (position) {
+ if (yHasNoRootEdge) {
+ j <- which(y$edge[, 1] == y$edge[1])
+ y$edge[j, 1] <- next.node
+ } else y$edge[1] <- next.node
+ x$edge <- rbind(x$edge[1:i, ], c(next.node, x$edge[i, 2]), x$edge[-(1:i), ])
+ x$edge[i, 2] <- next.node
+ if (wbl) {
+ x$edge.length[i] <- x$edge.length[i] - position
+ x$edge.length <- c(x$edge.length[1:i], position, x$edge.length[-(1:i)])
}
- ## insert the left parenthesis:
- ## here 21 is: nchar("TheNodeWhereToGraftY") + 1
- X <- paste(substr(X, 1, i - 1), "(",
- substr(X, i, 21), sep = "")
- ## and insert `y':
- X <- gsub("TheNodeWhereToGraftY",
- paste("TheNodeWhereToGraftY", ",", Y,
- sep = ""), X)
+ i <- i + 1L
} else {
- xx <- paste(")", "TheNodeWhereToGraftY", sep = "")
- X <- gsub(xx, paste(",", Y, xx, sep = ""), X)
+ if (yHasNoRootEdge) {
+ j <- which(y$edge[, 1] == y$edge[1])
+ y$edge[j, 1] <- x$edge[i, 2]
+ } else y$edge[1] <- x$edge[i, 2]
}
+ x$edge <- rbind(x$edge[1:i, ], y$edge, x$edge[-(1:i), ])
+ if (wbl)
+ x$edge.length <- c(x$edge.length[1:i], y$edge.length, x$edge.length[-(1:i)])
+ })
+
+ x$tip.label <- c(x$tip.label, y$tip.label)
+
+ if (is.null(x$node.label)) {
+ if (!is.null(y$node.label))
+ x$node.label <- c(rep(NA, mx), y$node.label)
+ } else {
+ x$node.label <-
+ if (is.null(y$node.label)) c(x$node.label, rep(NA, my))
+ else c(x$node.label, y$node.label)
}
- phy <- read.tree(text = X)
- ## restore the labels:
- if (where <= nb.tip)
- phy$tip.label[which(phy$tip.label == "TheTipWhereToGraftY")] <-
- Tip.Label.where
- if (where > ROOT) {
- if (xHasNoNodeLabel) phy$node.label <- NULL
- else
- phy$node.label[which(phy$node.label == "TheNodeWhereToGraftY")] <-
- Node.Label.where
- }
- phy
+
+ n <- length(x$tip.label)
+ x$Nnode <- dim(x$edge)[1] + 1L - n
+
+ ## update the node labels before renumbering (this adds NA for
+ ## the added nodes, and drops the label for those deleted)
+ if (!is.null(x$node.label))
+ x$node.label <- x$node.label[-unique(x$edge[, 1])]
+
+ ## renumber nodes:
+ newNb <- integer(x$Nnode)
+ newNb[-ROOT] <- n + 1L
+ sndcol <- x$edge[, 2] < 0
+ ## executed from right to left, so newNb is modified before x$edge:
+ x$edge[sndcol, 2] <- newNb[-x$edge[sndcol, 2]] <-
+ (n + 2):(n + x$Nnode)
+ x$edge[, 1] <- newNb[-x$edge[, 1]]
+
+ if (!is.null(x$node.label))
+ x$node.label <- x$node.label[newNb[newNb == 0] - n]
+
+ x
}
projects / ape.git / blobdiff