From: paradis <paradis@6e262413-ae40-0410-9e79-b911bd7a66b7>
Date: Wed, 9 Sep 2009 13:24:53 +0000 (+0000)
Subject: a collection of bug fixes
X-Git-Url: https://git.donarmstrong.com/?p=ape.git;a=commitdiff_plain;h=3ad385892d75db5c646c92f0f631ae9c5e3da4f6

a collection of bug fixes


git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@87 6e262413-ae40-0410-9e79-b911bd7a66b7
---

diff --git a/ChangeLog b/ChangeLog
index 4402627..275ad2f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -10,6 +10,17 @@ BUG FIXES
     o write.nexus() failed to write correctly trees with a "TipLabel"
       attribute.
 
+    o rcoal() failed to compute branch lengths with very large n.
+
+    o A small bug was fixed in compar.cheverud() (thanks to Michael
+      Phelan for the fix).
+
+    o seg.sites() failed when passing a vector.
+
+    o drop.tip() sometimes shuffled tip labels.
+
+    o root() shuffled node labels with 'resolve.root = TRUE'.
+
 
 
 		CHANGES IN APE VERSION 2.3-2
diff --git a/DESCRIPTION b/DESCRIPTION
index e09314b..85aee0b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ape
 Version: 2.3-3
-Date: 2009-07-27
+Date: 2009-09-09
 Title: Analyses of Phylogenetics and Evolution
 Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Gangolf Jobb, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer, Damien de Vienne
 Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
diff --git a/R/Cheverud.R b/R/Cheverud.R
index ab7e24b..6b70a12 100644
--- a/R/Cheverud.R
+++ b/R/Cheverud.R
@@ -12,7 +12,9 @@
 # Evolution 55(11): 2143-2160
 compar.cheverud <- function(y, W, tolerance=1e-6, gold.tol=1e-4)
 {
-  W <- W - diag(W) # ensure diagonal is zero
+  ## fix by Michael Phelan
+  diag(W) <- 0 # ensure diagonal is zero
+  ## end of fix
   y <- as.matrix(y)
   if(dim(y)[2] != 1) stop("Error: y must be a single column vector.")
   D <- solve(diag(apply(t(W),2,sum)))
diff --git a/R/DNA.R b/R/DNA.R
index d74d188..dd9c60b 100644
--- a/R/DNA.R
+++ b/R/DNA.R
@@ -1,4 +1,4 @@
-## DNA.R (2009-05-19)
+## DNA.R (2009-09-06)
 
 ##   Manipulations and Comparisons of DNA Sequences
 
@@ -274,9 +274,13 @@ GC.content <- function(x) sum(base.freq(x)[2:3])
 seg.sites <- function(x)
 {
     if (is.list(x)) x <- as.matrix(x)
-    n <- dim(x)
-    s <- n[2]
-    n <- n[1]
+    if (is.vector(x)) n <- 1
+    else { # 'x' is a matrix
+        n <- dim(x)
+        s <- n[2]
+        n <- n[1]
+    }
+    if (n == 1) return(integer(0))
     ans <- .C("SegSites", x, n, s, integer(s),
               DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
     which(as.logical(ans[[4]]))
diff --git a/R/compar.ou.R b/R/compar.ou.R
index 6512e63..ada36d7 100644
--- a/R/compar.ou.R
+++ b/R/compar.ou.R
@@ -40,7 +40,9 @@ compar.ou <- function(x, phy, node = NULL, alpha = NULL)
     }
     W <- cophenetic.phylo(phy)
     dev <- function(p) {
-        M <- rowSums(exp(-p[1] * Wstart) - exp(-p[1] * Wend) * p[-(1:2)])
+        ##M <- rowSums(exp(-p[1] * Wstart) - exp(-p[1] * Wend) * p[-(1:2)])
+        ##M <- -rowSums(exp(-p[1] * Wstart) - exp(-p[1] * Wend) * p[-(1:2)])
+        M <- rowSums((exp(-p[1] * Wend) - exp(-p[1] * Wstart)) * p[-(1:2)])
         V <- exp(-p[1]*W) * (1 - exp(-2*p[1]*(Tmax - W/2)))
         nb.tip*log(2*pi*p[2]) + log(det(V)) +
           (t(x - M) %*% chol2inv(V) %*% (x - M)) / p[2]
diff --git a/R/dist.topo.R b/R/dist.topo.R
index 1248c49..f82dd65 100644
--- a/R/dist.topo.R
+++ b/R/dist.topo.R
@@ -11,7 +11,7 @@
 dist.topo <- function(x, y, method = "PH85")
 {
     if (method == "BHV01" && (is.null(x$edge.length) || is.null(y$edge.length)))
-      stop("trees must have branch lengths for Billera et al.'s distance.")
+        stop("trees must have branch lengths for Billera et al.'s distance.")
     n <- length(x$tip.label)
     bp1 <- .Call("bipartition", x$edge, n, x$Nnode, PACKAGE = "ape")
     bp1 <- lapply(bp1, function(xx) sort(x$tip.label[xx]))
diff --git a/R/drop.tip.R b/R/drop.tip.R
index acdbe68..2d2a8f0 100644
--- a/R/drop.tip.R
+++ b/R/drop.tip.R
@@ -1,4 +1,4 @@
-## drop.tip.R (2009-07-06)
+## drop.tip.R (2009-09-09)
 
 ##   Remove Tips in a Phylogenetic Tree
 
@@ -167,8 +167,9 @@ drop.tip <-
 
     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 of their
+    ## in the 2nd col of edge, so no need to reorder $tip.label
+    phy$edge[TERMS, 2] <- rank(phy$edge[TERMS, 2])
 
     ## make new tip labels if necessary:
     if (subtree || !trim.internal) {
diff --git a/R/root.R b/R/root.R
index 3abdf65..aadd6c1 100644
--- a/R/root.R
+++ b/R/root.R
@@ -1,4 +1,4 @@
-## root.R (2009-07-06)
+## root.R (2009-09-09)
 
 ##   Root of Phylogenetic Trees
 
@@ -297,14 +297,19 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE)
     phy$edge[, 1] <- newNb[phy$edge[, 1]]
 
     if (!is.null(phy$node.label)) {
+        #browser()
         newNb <- newNb[-(1:n)]
         if (fuseRoot) {
             newNb <- newNb[-1]
             phy$node.label <- phy$node.label[-1]
         }
         phy$node.label <- phy$node.label[order(newNb)]
-        if (resolve.root)
-             phy$node.label <- c(phy$node.label[1], NA, phy$node.label[-1])
+        if (resolve.root) {
+            phy$node.label[is.na(phy$node.label)] <- phy$node.label[1]
+            phy$node.label[1] <- NA
+            ##phy$node.label <- c(phy$node.label[1], NA, phy$node.label[-1])
+            ##phy$node.label <- c("NA", phy$node.label)
+        }
     }
     phy
 }
diff --git a/R/rtree.R b/R/rtree.R
index f331ca3..6829b02 100644
--- a/R/rtree.R
+++ b/R/rtree.R
@@ -106,7 +106,7 @@ rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
     nbr <- 2*n - 2
     edge <- matrix(NA, nbr, 2)
     ## coalescence times by default:
-    x <- if (is.character(br)) 2*rexp(n - 1)/(n:2 * (n - 1):1)
+    x <- if (is.character(br)) 2*rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1))
     else br(n - 1, ...)
     if (n == 2) {
         edge[] <- c(3L, 3L, 1:2)
diff --git a/Thanks b/Thanks
index c9fea53..da83d4f 100644
--- a/Thanks
+++ b/Thanks
@@ -6,8 +6,9 @@ comments, or bug reports: thanks to all of you!
 
 Significant bug fixes were provided by Cécile Ané, James Bullard,
 Éric Durand, Olivier François, Rich FitzJohn, Bret Larget, Nick Matzke,
-Elizabeth Purdom, Dan Rabosky, Klaus Schliep, Tim Wallstrom, Li-San
-Wang, Yan Wong, and Peter Wragg. Contact me if I forgot someone.
+Michael Phelan, Elizabeth Purdom, Dan Rabosky, Klaus Schliep, Tim
+Wallstrom, Li-San Wang, Yan Wong, and Peter Wragg. Contact me if I
+forgot someone.
 
 Kurt Hornik, of the R Core Team, helped in several occasions to
 fix some problems and bugs.
diff --git a/src/nj.c b/src/nj.c
index 355c68f..bc52b2a 100644
--- a/src/nj.c
+++ b/src/nj.c
@@ -135,7 +135,7 @@ void nj(double *D, int *N, int *edge1, int *edge2, double *edge_length)
 			for (i = OTU1 - 1; i > 0; i--)
 				otu_label[i] = otu_label[i - 1];
 		if (OTU2 != n)
-			for (i = OTU2; i <= n; i++)
+			for (i = OTU2; i < n; i++)
 				otu_label[i - 1] = otu_label[i];
 		otu_label[0] = cur_nod;