1 ## multi2di.R (2008-04-09)
3 ## Collapse and Resolve Multichotomies
5 ## Copyright 2005-2008 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 nextnode <- length(phy$tip.label) + phy$Nnode + 1
17 new.edge <- edge2delete <- NULL
19 if (!is.null(phy$edge.length)) {
21 new.edge.length <- NULL
24 for (node in target) {
25 ind <- which(phy$edge[, 1] == node)
27 desc <- phy$edge[ind, 2]
29 ## if we shuffle the descendants, we need to eventually
30 ## reorder the corresponding branch lenghts (see below)
31 ## so we store the result of sample()
32 tmp <- sample(length(desc))
36 res <- matrix(0, 2*N - 2, 2)
37 res[, 1] <- N + rep(1:(N - 1), each = 2)
38 res[, 2] <- N + rep(2:N, each = 2)
39 res[seq(1, by = 2, length.out = N - 1), 2] <- 1:(N - 1)
43 ## keep the branch lengths coming from `node'
44 el <- numeric(dim(res)[1]) # initialized with 0's
46 if (random) phy$edge.length[ind][tmp] else phy$edge.length[ind]
48 ## now substitute the nodes in `res'
49 ## `node' stays at the "root" of these new
50 ## edges whereas their "tips" are `desc'
51 Nodes <- c(node, seq(from = nextnode, length.out = N - 2))
52 res[, 1] <- Nodes[res[, 1] - N]
54 res[tmp, 2] <- Nodes[res[tmp, 2] - N]
55 res[!tmp, 2] <- desc[res[!tmp, 2]]
56 new.edge <- rbind(new.edge, res)
57 edge2delete <- c(edge2delete, ind)
58 if (wbl) new.edge.length <- c(new.edge.length, el)
59 nextnode <- nextnode + N - 2
60 phy$Nnode <- phy$Nnode + N - 2
62 phy$edge <- rbind(phy$edge[-edge2delete, ], new.edge)
64 phy$edge.length <- c(phy$edge.length[-edge2delete], new.edge.length)
65 if (!is.null(attr(phy, "order"))) attr(phy, "order") <- NULL
66 if (!is.null(phy$node.label))
68 c(phy$node.label, rep("", phy$Nnode - length(phy$node.label)))
70 ##read.tree(text = write.tree(phy))
73 di2multi <- function(phy, tol = 1e-8)
75 if (is.null(phy$edge.length)) stop("the tree has no branch length")
76 ## We select only the internal branches which are
77 ## significantly small:
78 ind <- which(phy$edge.length < tol & phy$edge[, 2] > length(phy$tip.label))
81 ## recursive function to `propagate' node #'s in case
82 ## there is a series of consecutive edges to remove
83 foo <- function(ancestor, des2del) {
84 wh <- which(phy$edge[, 1] == des2del)
86 if (phy$edge[k, 2] %in% node2del) foo(ancestor, phy$edge[k, 2])
87 else phy$edge[k, 1] <<- ancestor
90 node2del <- phy$edge[ind, 2]
91 anc <- phy$edge[ind, 1]
93 if (anc[i] %in% node2del) next
94 foo(anc[i], node2del[i])
96 phy$edge <- phy$edge[-ind, ]
97 phy$edge.length <- phy$edge.length[-ind]
98 phy$Nnode <- phy$Nnode - n
99 ## Now we renumber the nodes that need to be:
100 sel <- phy$edge > min(node2del)
101 for (i in which(sel))
102 phy$edge[i] <- phy$edge[i] - sum(node2del < phy$edge[i])
103 if (!is.null(phy$node.label))
104 phy$node.label <- phy$node.label[-(node2del - length(phy$tip.label))]