R/which.edge.R

   1 ## which.edge.R (2009-05-10)
   2
   3 ##   Identifies Edges of a Tree
   4
   5 ## Copyright 2004-2009 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 getMRCA <- function(phy, tip)
  11 ### Find the MRCA of the tips given as `tip'
  12 ### (see `root.R' for comments on the code)
  13 {
  14     Ntip <- length(phy$tip.label)
  15     seq.nod <- .Call("seq_root2tip", phy$edge, Ntip,
  16                      phy$Nnode, PACKAGE = "ape")
  17     sn <- seq.nod[tip]
  18     MRCA <- Ntip + 1
  19     i <- 2
  20     repeat {
  21         x <- unique(unlist(lapply(sn, "[", i)))
  22         if (length(x) != 1) break
  23         MRCA <- x
  24         i <- i + 1
  25     }
  26     MRCA
  27 }
  28
  29 which.edge <- function(phy, group)
  30 {
  31     if (!inherits(phy, "phylo"))
  32       stop('object "phy" is not of class "phylo"')
  33     if (is.character(group))
  34       group <- which(phy$tip.label %in% group)
  35     if (length(group) == 1)
  36       return(match(group, phy$edge[, 2]))
  37     nb.tip <- length(phy$tip.label)
  38     MRCA <- getMRCA(phy, group)
  39     if (MRCA == nb.tip + 1) {
  40         from <- 1
  41         to <- dim(phy$edge)[1]
  42     } else {
  43         from <- which(phy$edge[, 2] == MRCA) + 1
  44         to <- max(which(phy$edge[, 2] %in% group))
  45     }
  46     wh <- from:to
  47     tmp <- phy$edge[wh, 2]
  48     ## check that there are no extra tips:
  49     ## (we do this by selecting the tips in `group' and the nodes
  50     ##  i.e., the internal edges)
  51     test <- tmp %in% group | tmp > nb.tip
  52     if (any(!test)) {
  53         wh <- wh[test] # drop the extra tips
  54         ## see if there are no extra internal edges:
  55         tmp <- phy$edge[wh, ]
  56         test <- !(tmp[, 2] %in% tmp[, 1]) & tmp[, 2] > nb.tip
  57         while (any(test)){
  58             wh <- wh[!test]
  59             tmp <- phy$edge[wh, ]
  60             test <- !(tmp[, 2] %in% tmp[, 1]) & tmp[, 2] > nb.tip
  61         }
  62     }
  63     wh
  64 }