new is.monophyletic() + various changes/fixes

[ape.git] / R / drop.tip.R

diff --git a/R/drop.tip.R b/R/drop.tip.R
index d819c757ce1c845d5165a555b9063668f5460962..a4a3d030207f2baa34d5a63ef535826fd11bce27 100644 (file)
--- a/R/drop.tip.R
+++ b/R/drop.tip.R
@@ -1,4 +1,4 @@
-## drop.tip.R (2009-01-07)
+## drop.tip.R (2009-05-10)
 
 ##   Remove Tips in a Phylogenetic Tree
 
@@ -22,7 +22,8 @@ extract.clade <- function(phy, node, root.edge = 0)
             stop("the tree has no node labels")
         node <- which(phy$node.label %in% node) + Ntip
     }
-    if (node <= Ntip) stop("node number must be greater than the number of tips")
+    if (node <= Ntip)
+        stop("node number must be greater than the number of tips")
     if (node == ROOT) return(phy)
     phy <- reorder(phy) # insure it is in cladewise order
     root.node <- which(phy$edge[, 2] == node)
@@ -73,8 +74,9 @@ extract.clade <- function(phy, node, root.edge = 0)
 drop.tip <-
     function(phy, tip, trim.internal = TRUE, subtree = FALSE, root.edge = 0)
 {
-    if (class(phy) != "phylo")
+    if (!inherits(phy, "phylo"))
         stop('object "phy" is not of class "phylo"')
+    phy <- reorder(phy)
     Ntip <- length(phy$tip.label)
     NEWROOT <- ROOT <- Ntip + 1
     Nnode <- phy$Nnode
@@ -94,24 +96,23 @@ drop.tip <-
     ## find the tips to drop:
     if (is.character(tip))
         tip <- which(phy$tip.label %in% tip)
-    trms <- edge2 <= Ntip
     ## delete the terminal edges given by `tip':
     keep[match(tip, edge2)] <- FALSE
 
     if (trim.internal) {
-        ## delete the internal edges that do not have descendants
-        ## anymore (ie, they are in the 2nd column of `edge' but
+        internals <- edge2 <= Ntip
+        ## delete the internal edges that do not have anymore
+        ## descendants (ie, they are in the 2nd col of `edge' but
         ## not in the 1st one)
         repeat {
-            sel <- !(edge2 %in% edge1[keep]) & !trms & keep
+            sel <- !(edge2 %in% edge1[keep]) & internals & keep
             if (!sum(sel)) break
             keep[sel] <- FALSE
         }
         if (subtree) {
             ## keep the subtending edge(s):
             subt <- edge1 %in% edge1[keep] & edge1 %in% edge1[!keep]
-            ## <FIXME> 'if (... ' needed below?
-            if (any(subt)) keep[which(subt)] <- TRUE
+            keep[subt] <- TRUE
         }
         if (root.edge && wbl) {
             degree <- tabulate(edge1[keep])
@@ -143,8 +144,11 @@ drop.tip <-
     ## 2) renumber the remaining tips now
     TIPS <- phy$edge[, 2] <= Ntip
     ## keep the ordering so no need to reorder tip.label:
-    phy$edge[TIPS, 2] <- order(phy$edge[TIPS, 2])
-    Ntip <- length(phy$tip.label) # update Ntip
+    phy$edge[TIPS, 2] <- rank(phy$edge[TIPS, 2])
+    ## 3) update node.label if needed
+    if (!is.null(phy$node.label))
+        phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
+    Ntip <- length(phy$tip.label) # 4) update Ntip
 
     ## make new tip labels if necessary
     if (subtree || !trim.internal) {
@@ -176,6 +180,6 @@ drop.tip <-
     phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
         (Ntip + 2):(Ntip + phy$Nnode)
     phy$edge[, 1] <- newNb[phy$edge[, 1]]
-
+    storage.mode(phy$edge) <- "integer"
     collapse.singles(phy)
 }