fix in birthdeath()

[ape.git] / R / read.tree.R

diff --git a/R/read.tree.R b/R/read.tree.R
index cbcfd3f9c4a6e418cdf208fcd4075f15874ace5b..cba51d10fd7a5641568b12480006bc6c7635078e 100644 (file)
--- a/R/read.tree.R
+++ b/R/read.tree.R
@@ -1,8 +1,8 @@
-## read.tree.R (2008-02-18)
+## read.tree.R (2010-09-27)
 
 ##   Read Tree Files in Parenthetic Format
 
-## Copyright 2002-2008 Emmanuel Paradis
+## Copyright 2002-2010 Emmanuel Paradis, Daniel Lawson and Klaus Schliep
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
@@ -12,11 +12,13 @@ tree.build <- function(tp)
     add.internal <- function() {
         edge[j, 1] <<- current.node
         edge[j, 2] <<- current.node <<- node <<- node + 1L
+        index[node] <<- j # set index
         j <<- j + 1L
     }
     add.terminal <- function() {
         edge[j, 1] <<- current.node
         edge[j, 2] <<- tip
+        index[tip] <<- j # set index
         X <- unlist(strsplit(tpc[k], ":"))
         tip.label[tip] <<- X[1]
         edge.length[j] <<- as.numeric(X[2])
@@ -25,7 +27,7 @@ tree.build <- function(tp)
         j <<- j + 1L
     }
     go.down <- function() {
-        l <- which(edge[, 2] == current.node)
+        l <- index[current.node]
         X <- unlist(strsplit(tpc[k], ":"))
         node.label[current.node - nb.tip] <<- X[1]
         edge.length[l] <<- as.numeric(X[2])
@@ -38,7 +40,7 @@ tree.build <- function(tp)
         obj$edge.length <- as.numeric(tp[3])
         obj$Nnode <- 1L
         obj$tip.label <- tp[2]
-        if (length(tp) == 4) obj$node.label <- tp[4]
+        if (tp[4] != "") obj$node.label <- tp[4]
         class(obj) <- "phylo"
         return(obj)
     }
@@ -60,7 +62,9 @@ tree.build <- function(tp)
     edge.length <- numeric(nb.edge)
     edge <- matrix(0L, nb.edge, 2)
     current.node <- node <- as.integer(nb.tip + 1) # node number
-    edge[nb.edge, 2] <- node #
+    edge[nb.edge, 2] <- node
+    index <- numeric(nb.edge + 1) # hash index to avoid which
+    index[node] <- nb.edge
 
     ## j: index of the line number of edge
     ## k: index of the line number of tpc
@@ -94,9 +98,19 @@ tree.build <- function(tp)
     obj
 }
 
-read.tree <- function(file = "", text = NULL, tree.names = NULL,
-                      skip = 0, comment.char = "#", ...)
+read.tree <- function(file = "", text = NULL, tree.names = NULL, skip = 0,
+    comment.char = "#", keep.multi = FALSE, ...)
 {
+    unname <- function(treetext) {
+        nc <- nchar(treetext)
+       tstart <- 1
+       while (substr(treetext, tstart, tstart) != "(" && tstart <= nc)
+            tstart <- tstart + 1
+       if (tstart > 1)
+            return(c(substr(treetext, 1, tstart - 1),
+                     substr(treetext, tstart, nc)))
+       return(c("", treetext))
+    }
     if (!is.null(text)) {
         if (!is.character(text))
           stop("argument `text' must be of mode character")
@@ -120,6 +134,13 @@ read.tree <- function(file = "", text = NULL, tree.names = NULL,
     STRING <- character(Ntree)
     for (i in 1:Ntree)
         STRING[i] <- paste(tree[x[i]:y[i]], sep = "", collapse = "")
+
+    tmp <- unlist(lapply(STRING, unname))
+    tmpnames <- tmp[c(TRUE, FALSE)]
+    STRING <- tmp[c(FALSE, TRUE)]
+    if (is.null(tree.names) && any(nzchar(tmpnames)))
+        tree.names <- tmpnames
+
     colon <- grep(":", STRING)
     if (!length(colon)) {
         obj <- lapply(STRING, clado.build)
@@ -137,9 +158,9 @@ read.tree <- function(file = "", text = NULL, tree.names = NULL,
         ## is a bifurcation at the root
         ROOT <- length(obj[[i]]$tip.label) + 1
         if(sum(obj[[i]]$edge[, 1] == ROOT) == 1 && dim(obj[[i]]$edge)[1] > 1)
-            stop(paste("There is apparently two root edges in your file: cannot read tree file.\n  Reading Newick file aborted at tree no.", i, sep = ""))
+            stop(paste("The tree has apparently singleton node(s): cannot read tree file.\n  Reading Newick file aborted at tree no.", i))
     }
-    if (Ntree == 1) obj <- obj[[1]] else {
+    if (Ntree == 1 && !keep.multi) obj <- obj[[1]] else {
         if (!is.null(tree.names)) names(obj) <- tree.names
         class(obj) <- "multiPhylo"
     }