3 ## Find Most Recent Common Ancestors Between Pairs
5 ## Copyright 2005-2006 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 mrca <- function(phy, full = FALSE)
12 if (class(phy) != "phylo") stop('object "phy" is not of class "phylo"')
13 ## if (!is.rooted(phy)) stop("the tree must be rooted.")
15 nb.tip <- length(phy$tip.label)
17 BP <- .Call("bipartition", phy$edge, nb.tip,
18 nb.node, PACKAGE = "ape")
21 ## In the following matrix, numeric indexing will be used:
25 ## We start at the root:
27 while (length(next.node)) {
29 for (anc in next.node) {
30 ## Find the branches which `anc' is the ancestor...:
31 id <- which(phy$edge[, 1] == anc)
32 ## ... and get their descendants:
33 desc <- phy$edge[id, 2]
34 ## `anc' is itself the MRCA of its direct descendants:
35 M[anc, desc] <- M[desc, anc] <- anc
36 ## Find all 2-by-2 combinations of `desc': `anc'
38 for (i in 1:length(desc))
39 M[cbind(desc[i], desc[-i])] <- anc
40 ## If one element of `desc' is a node, then the tips it
41 ## leads to and the other elements of `desc' have also
43 for (i in 1:length(desc)) {
44 if (desc[i] < ROOT) next
46 tips <- BP[[desc[i] - nb.tip]]
47 ## Same thing for the nodes...
48 node.desc <- numeric(0)
49 for (k in 1:nb.node) {
50 if (k == desc[i] - nb.tip) next
51 ## If the clade of the current node is a
52 ## subset of desc[i], then it is one of its
54 if (all(BP[[k]] %in% tips))
55 node.desc <- c(node.desc, k)
57 ## all nodes and tips which are descendants of
59 ALLDESC <- c(tips, node.desc + nb.tip)
60 M[ALLDESC, desc[-i]] <- M[desc[-i], ALLDESC] <- anc
61 for (j in 1:length(desc)) {
62 if (j == i || desc[j] < ROOT) next
63 tips2 <- BP[[desc[j] - nb.tip]]
64 node.desc <- numeric(0)
65 for (k in 1:nb.node) {
66 if (k == desc[j] - nb.tip) next
67 if (all(BP[[k]] %in% tips2))
68 node.desc <- c(node.desc, k)
70 ALLDESC2 <- c(tips2, node.desc + nb.tip)
71 M[ALLDESC, ALLDESC2] <- M[ALLDESC2, ALLDESC] <- anc
73 ## `anc' is also the MRCA of itself and its descendants:
74 M[ALLDESC, anc] <- M[anc, ALLDESC] <- anc
76 ## When it is done, `desc' i stored to become
77 ## the new `next.node', if they are nodes:
78 tmp <- c(tmp, desc[desc > nb.tip])
82 M[cbind(1:N, 1:N)] <- 1:N
84 dimnames(M)[1:2] <- list(as.character(1:N))
86 M <- M[1:nb.tip, 1:nb.tip]
87 dimnames(M)[1:2] <- list(phy$tip.label)