R/multi2di.R

   1 ## multi2di.R (2008-03-17)
   2
   3 ##   Collapse and Resolve Multichotomies
   4
   5 ## Copyright 2005-2008 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     nextnode <- length(phy$tip.label) + phy$Nnode + 1
  17     new.edge <- edge2delete <- NULL
  18     wbl <- FALSE
  19     if (!is.null(phy$edge.length)) {
  20         wbl <- TRUE
  21         new.edge.length <- NULL
  22     }
  23
  24     for (node in target) {
  25         ind <- which(phy$edge[, 1] == node)
  26         N <- length(ind)
  27         desc <- phy$edge[ind, 2]
  28         if (random) {
  29           ## if we shuffle the descendants, we need to eventually
  30           ## reorder the corresponding branch lenghts (see below)
  31           ## so we store the result of sample()
  32             tmp <- sample(length(desc))
  33             desc <- desc[tmp]
  34             res <- rtree(N)$edge
  35         } else {
  36             res <- matrix(0, 2*N - 2, 2)
  37             res[, 1] <- N + rep(1:(N - 1), each = 2)
  38             res[, 2] <- N + rep(2:N, each = 2)
  39             res[seq(1, by = 2, length.out = N - 1), 2] <- 1:(N - 1)
  40             res[length(res)] <- N
  41         }
  42         if (wbl) {
  43             ## keep the branch lengths coming from `node'
  44             el <- numeric(dim(res)[1]) # initialized with 0's
  45             el[res[, 2] <= N] <-
  46               if (random) phy$edge.length[ind][tmp] else phy$edge.length[ind]
  47         }
  48         ## now substitute the nodes in `res'
  49         ## `node' stays at the "root" of these new
  50         ## edges whereas their "tips" are `desc'
  51         Nodes <- c(node, seq(from = nextnode, length.out = N - 2))
  52         res[, 1] <- Nodes[res[, 1] - N]
  53         tmp <- res[, 2] > N
  54         res[tmp, 2] <- Nodes[res[tmp, 2] - N]
  55         res[!tmp, 2] <- desc[res[!tmp, 2]]
  56         new.edge <- rbind(new.edge, res)
  57         edge2delete <- c(edge2delete, ind)
  58         if (wbl) new.edge.length <- c(new.edge.length, el)
  59         nextnode <- nextnode + N - 2
  60         phy$Nnode <- phy$Nnode + N - 2
  61     }
  62     phy$edge <- rbind(phy$edge[-edge2delete, ], new.edge)
  63     if (wbl)
  64       phy$edge.length <- c(phy$edge.length[-edge2delete], new.edge.length)
  65     if (!is.null(attr(phy, "order"))) attr(phy, "order") <- NULL
  66     print(phy$node.label)
  67     if (!is.null(phy$node.label))
  68         phy$node.label <-
  69             c(phy$node.label, rep("", phy$Nnode - length(phy$node.label)))
  70     print(phy$node.label)
  71     reorder(phy)
  72     ##read.tree(text = write.tree(phy))
  73 }
  74
  75 di2multi <- function(phy, tol = 1e-8)
  76 {
  77     if (is.null(phy$edge.length)) stop("the tree has no branch length")
  78     ## We select only the internal branches which are
  79     ## significantly small:
  80     ind <- which(phy$edge.length < tol & phy$edge[, 2] > length(phy$tip.label))
  81     n <- length(ind)
  82     if (!n) return(phy)
  83     ## recursive function to `propagate' node #'s in case
  84     ## there is a series of consecutive edges to remove
  85     foo <- function(ancestor, des2del) {
  86         wh <- which(phy$edge[, 1] == des2del)
  87         for (k in wh) {
  88             if (phy$edge[k, 2] %in% node2del) foo(ancestor, phy$edge[k, 2])
  89             else phy$edge[k, 1] <<- ancestor
  90         }
  91     }
  92     node2del <- phy$edge[ind, 2]
  93     anc <- phy$edge[ind, 1]
  94     for (i in 1:n) {
  95         if (anc[i] %in% node2del) next
  96         foo(anc[i], node2del[i])
  97     }
  98     phy$edge <- phy$edge[-ind, ]
  99     phy$edge.length <- phy$edge.length[-ind]
 100     phy$Nnode <- phy$Nnode - n
 101     ## Now we renumber the nodes that need to be:
 102     sel <- phy$edge > min(node2del)
 103     for (i in which(sel))
 104       phy$edge[i] <- phy$edge[i] - sum(node2del < phy$edge[i])
 105     if (!is.null(phy$node.label))
 106         phy$node.label <- phy$node.label[-(node2del - length(phy$tip.label))]
 107     phy
 108 }