X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Fdrop.tip.R;h=113c448eaeaf8805f51460e8ce3a511175e6b017;hb=c488b74490ee3d9d200de0e471881f002a18fe4f;hp=6e497cb9c0245dd66e39bbd79d291bd484cab0af;hpb=039cecc7ab5174412f1be3164f991f1ddc637f83;p=ape.git

diff --git a/R/drop.tip.R b/R/drop.tip.R
index 6e497cb..113c448 100644
--- a/R/drop.tip.R
+++ b/R/drop.tip.R
@@ -1,35 +1,37 @@
-## drop.tip.R (2009-06-23)
+## drop.tip.R (2012-11-29)
 
 ##   Remove Tips in a Phylogenetic Tree
 
-## Copyright 2003-2009 Emmanuel Paradis
+## Copyright 2003-2012 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
 
-extract.clade <- function(phy, node, root.edge = 0)
+extract.clade <- function(phy, node, root.edge = 0, interactive = FALSE)
 {
     Ntip <- length(phy$tip.label)
     ROOT <- Ntip + 1
     Nedge <- dim(phy$edge)[1]
     wbl <- !is.null(phy$edge.length)
-    if (length(node) > 1) {
-        node <- node[1]
-        warning("only the first value of 'node' has been considered")
-    }
-    if (is.character(node)) {
-        if (is.null(phy$node.label))
-            stop("the tree has no node labels")
-        node <- which(phy$node.label %in% node) + Ntip
+    if (interactive) node <- identify(phy)$nodes else {
+        if (length(node) > 1) {
+            node <- node[1]
+            warning("only the first value of 'node' has been considered")
+        }
+        if (is.character(node)) {
+            if (is.null(phy$node.label))
+                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)
     start <- root.node + 1 # start of the clade looked for
     anc <- phy$edge[root.node, 1] # the ancestor of 'node'
-    next.anc <- which(phy$edge[-(1:start), 1] == anc) # find the next occurence of 'anc'
+    next.anc <- which(phy$edge[-(1:start), 1] <= anc) # find the next occurence of 'anc' or an 'older' node
 
     keep <- if (length(next.anc)) start + 0:(next.anc[1] - 1) else start:Nedge
 
@@ -52,7 +54,11 @@ extract.clade <- function(phy, node, root.edge = 0)
     if (wbl) phy$edge.length <- phy$edge.length[keep]
     TIPS <- phy$edge[, 2] <= Ntip
     tip <- phy$edge[TIPS, 2]
-    phy$tip.label <- phy$tip.label[tip]
+    ## Fix by Ludovic Mallet and Mahendra Mariadassou (2011-11-21):
+    name <- vector("character", length(tip))
+    name[order(tip)] <- phy$tip.label[tip]
+    phy$tip.label <- name
+    ## End of fix
     ## keep the ordering so no need to reorder tip.label:
     phy$edge[TIPS, 2] <- order(tip)
     if (!is.null(phy$node.label))
@@ -73,15 +79,31 @@ extract.clade <- function(phy, node, root.edge = 0)
 
 drop.tip <-
     function(phy, tip, trim.internal = TRUE, subtree = FALSE,
-             root.edge = 0, rooted = is.rooted(phy))
+             root.edge = 0, rooted = is.rooted(phy), interactive = FALSE)
 {
     if (!inherits(phy, "phylo"))
         stop('object "phy" is not of class "phylo"')
+    if (!length(tip)) return(phy)
 
     Ntip <- length(phy$tip.label)
     ## find the tips to drop:
-    if (is.character(tip))
-        tip <- which(phy$tip.label %in% tip)
+    if (interactive) {
+        cat("Left-click close to the tips you want to drop; right-click when finished...\n")
+        xy <- locator()
+        nToDrop <- length(xy$x)
+        tip <- integer(nToDrop)
+        lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
+        for (i in 1:nToDrop) {
+            d <- sqrt((xy$x[i] - lastPP$xx)^2 + (xy$y[i] - lastPP$yy)^2)
+            tip[i] <- which.min(d)
+        }
+    } else {
+        if (is.character(tip))
+            tip <- which(phy$tip.label %in% tip)
+    }
+    if (!length(tip)) return(phy)
+    if (any(tip > Ntip))
+        warning("some tip numbers were higher than the number of tips")
 
     if (!rooted && subtree) {
         phy <- root(phy, (1:Ntip)[-tip][1])
@@ -105,15 +127,6 @@ drop.tip <-
     edge2 <- phy$edge[, 2] #
     keep <- !logical(Nedge)
 
-    ## find the tips to drop:
-    if (is.character(tip))
-        tip <- which(phy$tip.label %in% tip)
-
-    if (!rooted && subtree) {
-        phy <- root(phy, (1:Ntip)[-tip][1])
-        root.edge <- 0
-    }
-
     ## delete the terminal edges given by `tip':
     keep[match(tip, edge2)] <- FALSE
 
@@ -165,39 +178,42 @@ drop.tip <-
     ## get the old No. of the nodes and tips that become tips:
     oldNo.ofNewTips <- phy$edge[TERMS, 2]
 
+    ## in case some tips are dropped but kept because of 'subtree = TRUE':
+    if (subtree) {
+        i <- which(tip %in% oldNo.ofNewTips)
+        if (length(i)) {
+            phy$tip.label[tip[i]] <- "[1_tip]"
+            tip <- tip[-i]
+        }
+    }
+
     n <- length(oldNo.ofNewTips) # the new number of tips in the tree
 
-    ## assumes that the ordering of tips is unchanged:
-    phy$edge[TERMS, 2] <- 1:n
+    ## the tips may not be sorted in increasing order in the
+    ## 2nd col of edge, so no need to reorder $tip.label
+    phy$edge[TERMS, 2] <- rank(phy$edge[TERMS, 2])
+    phy$tip.label <- phy$tip.label[-tip]
 
     ## make new tip labels if necessary:
     if (subtree || !trim.internal) {
-        ## get the logical indices of the tips that are kept within 'oldNo.ofNewTips':
-        tips.kept <- oldNo.ofNewTips <= Ntip & !(oldNo.ofNewTips %in% tip)
-        new.tip.label <- character(n)
-        new.tip.label[tips.kept] <- phy$tip.label[-tip]
         ## get the numbers of the nodes that become tips:
-        node2tip <- oldNo.ofNewTips[!tips.kept]
-        new.tip.label[!tips.kept] <- if (subtree) {
+        node2tip <- oldNo.ofNewTips[oldNo.ofNewTips > Ntip]
+        new.tip.label <- if (subtree) {
             paste("[", N[node2tip], "_tips]", sep = "")
         } else {
             if (is.null(phy$node.label)) rep("NA", length(node2tip))
             else phy$node.label[node2tip - Ntip]
         }
-        if (!is.null(phy$node.label))
-            phy$node.label <- phy$node.label[-(node2tip - Ntip)]
-        phy$tip.label <- new.tip.label
-    } else phy$tip.label <- phy$tip.label[-tip]
-
-    ## update node.label if needed:
-    if (!is.null(phy$node.label))
-        phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
+#        if (!is.null(phy$node.label))
+#            phy$node.label <- phy$node.label[-(node2tip - Ntip)]
+        phy$tip.label <- c(phy$tip.label, new.tip.label)
+    }
 
     phy$Nnode <- dim(phy$edge)[1] - n + 1L # update phy$Nnode
 
     ## The block below renumbers the nodes so that they conform
-    ## to the "phylo" format -- same as in root()
-    newNb <- integer(n + phy$Nnode)
+    ## to the "phylo" format, same as in root()
+    newNb <- integer(Ntip + Nnode)
     newNb[NEWROOT] <- n + 1L
     sndcol <- phy$edge[, 2] > n
     ## executed from right to left, so newNb is modified before phy$edge:
@@ -205,5 +221,8 @@ drop.tip <-
         (n + 2):(n + phy$Nnode)
     phy$edge[, 1] <- newNb[phy$edge[, 1]]
     storage.mode(phy$edge) <- "integer"
+    if (!is.null(phy$node.label)) # update node.label if needed
+        phy$node.label <- phy$node.label[which(newNb > 0) - Ntip]
     collapse.singles(phy)
 }
+