1 ## cophenetic.phylo.R (2007-01-23)
3 ## Pairwise Distances from a Phylogenetic Tree
5 ## Copyright 2006-2007 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 dist.nodes <- function(x)
12 if (is.null(x$edge.length))
13 stop("your tree has no branch lengths")
15 if (!is.binary.tree(x) || !is.rooted(x))
16 x <- multi2di(x, random = FALSE)
17 n <- length(x$tip.label)
20 x <- reorder(x, order = "pruningwise")
22 res <- matrix(NA, N, N)
23 res[cbind(1:N, 1:N)] <- 0 # implicit mode conversion
25 ## I like the simplicity of this one:
26 res[x$edge] <- res[x$edge[, 2:1]] <- x$edge.length
28 ## compute the distances ...
29 for (i in seq(from = 1, by = 2, length.out = n.node)) {
35 ## If `des1' is a node, we look for the nodes and tips for
36 ## which the distance up to `des1' has already been
37 ## computed, including `des1' itself. For all these, we can
38 ## compute the distance up to `anc' and all node(s) and
40 if (des1 > n) des1 <- which(!is.na(res[des1, ]))
43 if (des2 > n) des2 <- which(!is.na(res[des2, ]))
45 ## The following expression is vectorized only on `des2' and
46 ## not on `des1' because they may be of different lengths.
48 res[y, des2] <- res[des2, y] <- res[anc, y] + res[anc, des2]
49 ## compute the distances between the tip(s) and node(s)
50 ## in `des2' and the ancestor of `anc'; id. for `des2'
51 ## (only if it is not the root)
53 ind <- which(x$edge[, 2] == anc)
54 nod <- x$edge[ind, 1] # the ancestor of `anc'
55 l <- x$edge.length[ind]
56 res[des2, nod] <- res[nod, des2] <- res[anc, des2] + l
57 res[des1, nod] <- res[nod, des1] <- res[anc, des1] + l
60 dimnames(res)[1:2] <- list(1:N)
64 cophenetic.phylo <- function(x)
66 n <- length(x$tip.label)
67 ans <- dist.nodes(x)[1:n, 1:n]
68 dimnames(ans)[1:2] <- list(x$tip.label)