R/multi2di.R

   1 ## multi2di.R (2010-01-23)
   2
   3 ##   Collapse and Resolve Multichotomies
   4
   5 ## Copyright 2005-2010 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 multi2di <- function(phy, random = TRUE)
  11 {
  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     n <- length(phy$tip.label)
  17     nextnode <- n + phy$Nnode + 1
  18     new.edge <- edge2delete <- NULL
  19     wbl <- FALSE
  20     if (!is.null(phy$edge.length)) {
  21         wbl <- TRUE
  22         new.edge.length <- NULL
  23     }
  24
  25     for (node in target) {
  26         ind <- which(phy$edge[, 1] == node)
  27         N <- length(ind)
  28         desc <- phy$edge[ind, 2]
  29         if (random) {
  30           ## if we shuffle the descendants, we need to eventually
  31           ## reorder the corresponding branch lenghts (see below)
  32           ## so we store the result of sample()
  33             tmp <- sample(length(desc))
  34             desc <- desc[tmp]
  35             res <- rtree(N)$edge
  36         } else {
  37             res <- matrix(0, 2*N - 2, 2)
  38             res[, 1] <- N + rep(1:(N - 1), each = 2)
  39             res[, 2] <- N + rep(2:N, each = 2)
  40             res[seq(1, by = 2, length.out = N - 1), 2] <- 1:(N - 1)
  41             res[length(res)] <- N
  42         }
  43         if (wbl) {
  44             ## keep the branch lengths coming from `node'
  45             el <- numeric(dim(res)[1]) # initialized with 0's
  46             el[res[, 2] <= N] <-
  47               if (random) phy$edge.length[ind][tmp] else phy$edge.length[ind]
  48         }
  49         ## now substitute the nodes in `res'
  50         ## `node' stays at the "root" of these new
  51         ## edges whereas their "tips" are `desc'
  52         Nodes <- c(node, seq(from = nextnode, length.out = N - 2))
  53         res[, 1] <- Nodes[res[, 1] - N]
  54         tmp <- res[, 2] > N
  55         res[tmp, 2] <- Nodes[res[tmp, 2] - N]
  56         res[!tmp, 2] <- desc[res[!tmp, 2]]
  57         new.edge <- rbind(new.edge, res)
  58         edge2delete <- c(edge2delete, ind)
  59         if (wbl) new.edge.length <- c(new.edge.length, el)
  60         nextnode <- nextnode + N - 2
  61         phy$Nnode <- phy$Nnode + N - 2
  62     }
  63     phy$edge <- rbind(phy$edge[-edge2delete, ], new.edge)
  64     if (wbl)
  65         phy$edge.length <- c(phy$edge.length[-edge2delete], new.edge.length)
  66     if (!is.null(attr(phy, "order"))) attr(phy, "order") <- NULL
  67     if (!is.null(phy$node.label))
  68         phy$node.label <-
  69             c(phy$node.label, rep("", phy$Nnode - length(phy$node.label)))
  70     phy <- reorder(phy)
  71
  72     ## the node numbers are not in increasing order in edge[, 2]: this
  73     ## will confuse drop.tip and other functions (root), so renumber them
  74     newNb <- integer(phy$Nnode)
  75     newNb[1] <- n + 1L
  76     sndcol <- phy$edge[, 2] > n
  77
  78     ## reorder node labels before changing edge:
  79     if (!is.null(phy$node.label)) {
  80         o <- 1 + rank(phy$edge[sndcol, 2])
  81         ## the root's label is not changed:
  82         phy$node.label <- phy$node.label[c(1, o)]
  83     }
  84
  85     ## executed from right to left, so newNb is modified before phy$edge:
  86     phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2] - n] <-
  87         n + 2:phy$Nnode
  88     phy$edge[, 1] <- newNb[phy$edge[, 1] - n]
  89     phy
  90 }
  91
  92 di2multi <- function(phy, tol = 1e-8)
  93 {
  94     if (is.null(phy$edge.length)) stop("the tree has no branch length")
  95     ## We select only the internal branches which are
  96     ## significantly small:
  97     ind <- which(phy$edge.length < tol & phy$edge[, 2] > length(phy$tip.label))
  98     n <- length(ind)
  99     if (!n) return(phy)
 100     ## recursive function to `propagate' node #'s in case
 101     ## there is a series of consecutive edges to remove
 102     foo <- function(ancestor, des2del) {
 103         wh <- which(phy$edge[, 1] == des2del)
 104         for (k in wh) {
 105             if (phy$edge[k, 2] %in% node2del) foo(ancestor, phy$edge[k, 2])
 106             else phy$edge[k, 1] <<- ancestor
 107         }
 108     }
 109     node2del <- phy$edge[ind, 2]
 110     anc <- phy$edge[ind, 1]
 111     for (i in 1:n) {
 112         if (anc[i] %in% node2del) next
 113         foo(anc[i], node2del[i])
 114     }
 115     phy$edge <- phy$edge[-ind, ]
 116     phy$edge.length <- phy$edge.length[-ind]
 117     phy$Nnode <- phy$Nnode - n
 118     ## Now we renumber the nodes that need to be:
 119     sel <- phy$edge > min(node2del)
 120     for (i in which(sel))
 121       phy$edge[i] <- phy$edge[i] - sum(node2del < phy$edge[i])
 122     if (!is.null(phy$node.label))
 123         phy$node.label <- phy$node.label[-(node2del - length(phy$tip.label))]
 124     phy
 125 }