1 ## multi2di.R (2010-01-23)
3 ## Collapse and Resolve Multichotomies
5 ## Copyright 2005-2010 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 multi2di <- function(phy, random = TRUE)
12 degree <- tabulate(phy$edge[, 1])
13 target <- which(degree > 2)
14 if (!length(target)) return(phy)
15 nb.edge <- dim(phy$edge)[1]
16 n <- length(phy$tip.label)
17 nextnode <- n + phy$Nnode + 1
18 new.edge <- edge2delete <- NULL
20 if (!is.null(phy$edge.length)) {
22 new.edge.length <- NULL
25 for (node in target) {
26 ind <- which(phy$edge[, 1] == node)
28 desc <- phy$edge[ind, 2]
30 ## if we shuffle the descendants, we need to eventually
31 ## reorder the corresponding branch lenghts (see below)
32 ## so we store the result of sample()
33 tmp <- sample(length(desc))
37 res <- matrix(0, 2*N - 2, 2)
38 res[, 1] <- N + rep(1:(N - 1), each = 2)
39 res[, 2] <- N + rep(2:N, each = 2)
40 res[seq(1, by = 2, length.out = N - 1), 2] <- 1:(N - 1)
44 ## keep the branch lengths coming from `node'
45 el <- numeric(dim(res)[1]) # initialized with 0's
47 if (random) phy$edge.length[ind][tmp] else phy$edge.length[ind]
49 ## now substitute the nodes in `res'
50 ## `node' stays at the "root" of these new
51 ## edges whereas their "tips" are `desc'
52 Nodes <- c(node, seq(from = nextnode, length.out = N - 2))
53 res[, 1] <- Nodes[res[, 1] - N]
55 res[tmp, 2] <- Nodes[res[tmp, 2] - N]
56 res[!tmp, 2] <- desc[res[!tmp, 2]]
57 new.edge <- rbind(new.edge, res)
58 edge2delete <- c(edge2delete, ind)
59 if (wbl) new.edge.length <- c(new.edge.length, el)
60 nextnode <- nextnode + N - 2
61 phy$Nnode <- phy$Nnode + N - 2
63 phy$edge <- rbind(phy$edge[-edge2delete, ], new.edge)
65 phy$edge.length <- c(phy$edge.length[-edge2delete], new.edge.length)
66 if (!is.null(attr(phy, "order"))) attr(phy, "order") <- NULL
67 if (!is.null(phy$node.label))
69 c(phy$node.label, rep("", phy$Nnode - length(phy$node.label)))
72 ## the node numbers are not in increasing order in edge[, 2]: this
73 ## will confuse drop.tip and other functions (root), so renumber them
74 newNb <- integer(phy$Nnode)
76 sndcol <- phy$edge[, 2] > n
78 ## reorder node labels before changing edge:
79 if (!is.null(phy$node.label)) {
80 o <- 1 + rank(phy$edge[sndcol, 2])
81 ## the root's label is not changed:
82 phy$node.label <- phy$node.label[c(1, o)]
85 ## executed from right to left, so newNb is modified before phy$edge:
86 phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2] - n] <-
88 phy$edge[, 1] <- newNb[phy$edge[, 1] - n]
92 di2multi <- function(phy, tol = 1e-8)
94 if (is.null(phy$edge.length)) stop("the tree has no branch length")
95 ## We select only the internal branches which are
96 ## significantly small:
97 ind <- which(phy$edge.length < tol & phy$edge[, 2] > length(phy$tip.label))
100 ## recursive function to `propagate' node #'s in case
101 ## there is a series of consecutive edges to remove
102 foo <- function(ancestor, des2del) {
103 wh <- which(phy$edge[, 1] == des2del)
105 if (phy$edge[k, 2] %in% node2del) foo(ancestor, phy$edge[k, 2])
106 else phy$edge[k, 1] <<- ancestor
109 node2del <- phy$edge[ind, 2]
110 anc <- phy$edge[ind, 1]
112 if (anc[i] %in% node2del) next
113 foo(anc[i], node2del[i])
115 phy$edge <- phy$edge[-ind, ]
116 phy$edge.length <- phy$edge.length[-ind]
117 phy$Nnode <- phy$Nnode - n
118 ## Now we renumber the nodes that need to be:
119 sel <- phy$edge > min(node2del)
120 for (i in which(sel))
121 phy$edge[i] <- phy$edge[i] - sum(node2del < phy$edge[i])
122 if (!is.null(phy$node.label))
123 phy$node.label <- phy$node.label[-(node2del - length(phy$tip.label))]