getGroups w/ edges

getGroups.R

getGroups <- function(tree, get_edges=TRUE) {
  .fillInRes <- function(i) {
    res <- vector('list', length=2)
    names(res) <- c('g1', 'g2')
    res[['g1']] <- which(edge.tips[i, ] == 1)
    res[['g2']] <- which(edge.tips[i, ] == 0)
    res
  }
  .addEdges <- function(i) {
    res[[i]][['edglngth']] <<- edglngths[[i]]
    NULL
  }
  edge.tips <- .Call("bifurcations", PACKAGE="phylofactor",
                     as.integer(tree$edge[, 1]),
                     as.integer(tree$edge[, 2]),
                     as.integer(length(tree$tip.label)))
  res <- lapply(1:nrow(edge.tips), .fillInRes)
  names(res) <- 1:nrow(tree$edge) # remove this line after testing
  if(is.rooted(tree)) {
    root <- length(tree$tip.label) + 1
    rootedges <- which(tree$edge[ ,1] == root)
    edglngths <- tree$edge.length
    rootedglngth <- sum(edglngths[rootedges])
    edglngths[rootedges[1]] <- rootedglngth
    pull <- 1:nrow(tree$edge) %in% rootedges[-1]
    edglngths <- edglngths[!pull]
    res <- res[!pull]
  } else {
    edglngths <- tree$edge.length
  }
  if(get_edges) {
    lapply(1:length(res), .addEdges)
  }
  res
}
# n.b. having labels in 'res' will make the code marginally slower

# TESTING
plotWEdgeInf <- function() {
  # plot w/ edge number and length
  plot(tree, no.margin=TRUE)
  edgelabels(text=round(tree$edge.length, 2), adj=-.5)
  edgelabels()
}
polytomise <- function(tree, n) {
  # make polytomies for random n internal nodes
  intnds <- length(tree$tip.label):nrow(tree$edge)
  edges <- which(tree$edge[ ,2] %in% intnds)
  tree$edge.length[sample(edges, 3)] <- 0
  tree <- di2multi(tree)
  tree
}
library(phylofactor)
# test bifurcating and rooted tree
tree <- rtree(10)
tree$tip.label <- 1:10
plotWEdgeInf()
res <- getGroups(tree)
res[['1']]  # res contains edge numbers as names
# test polytmous and unrooted tree
tree <- unroot(tree)
plotWEdgeInf()
res <- getGroups(tree)
res[['1']]
# test tree with multiple polytomies
tree <- polytomise(rtree(10), 3)
tree$tip.label <- 1:10
plotWEdgeInf()
res <- getGroups(tree)

You might want to remove l18 in the future