From 0f717ab8186382152b581c9dd67bf85779112e7f Mon Sep 17 00:00:00 2001 From: paradis Date: Wed, 10 Jun 2009 09:49:21 +0000 Subject: [PATCH] improved ace() git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@76 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 4 ++++ DESCRIPTION | 2 +- R/ace.R | 27 ++++++++++++++------------- Thanks | 2 +- 4 files changed, 20 insertions(+), 15 deletions(-) diff --git a/ChangeLog b/ChangeLog index 8ad799a..5705d9a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -22,6 +22,10 @@ BUG FIXES o A small bug was fixed in CDAM.global(). + o ace() failed with large data sets. Thanks to Rich FitzJohn for + the fix. With other improvements, this function is now about 6 + times faster. + OTHER CHANGES diff --git a/DESCRIPTION b/DESCRIPTION index fac51bd..c12694c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape Version: 2.3-1 -Date: 2009-06-08 +Date: 2009-06-10 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 diff --git a/R/ace.R b/R/ace.R index b45a88e..91994be 100644 --- a/R/ace.R +++ b/R/ace.R @@ -1,4 +1,4 @@ -## ace.R (2009-05-10) +## ace.R (2009-06-10) ## Ancestral Character Estimation @@ -149,10 +149,13 @@ as the number of categories in `x'") rate[rate == 0] <- np + 1 # to avoid 0's since we will use this an numeric indexing liks <- matrix(0, nb.tip + nb.node, nl) - for (i in 1:nb.tip) liks[i, x[i]] <- 1 + TIPS <- 1:nb.tip + liks[cbind(TIPS, x)] <- 1 phy <- reorder(phy, "pruningwise") Q <- matrix(0, nl, nl) + ## from Rich FitzJohn: + comp <- numeric(nb.tip + nb.node) # Storage... dev <- function(p, output.liks = FALSE) { Q[] <- c(p, 0)[rate] diag(Q) <- -rowSums(Q) @@ -161,14 +164,14 @@ as the number of categories in `x'") anc <- phy$edge[i, 1] des1 <- phy$edge[i, 2] des2 <- phy$edge[j, 2] - tmp <- eigen(Q * phy$edge.length[i], symmetric = FALSE) - P1 <- tmp$vectors %*% diag(exp(tmp$values)) %*% solve(tmp$vectors) - tmp <- eigen(Q * phy$edge.length[j], symmetric = FALSE) - P2 <- tmp$vectors %*% diag(exp(tmp$values)) %*% solve(tmp$vectors) - liks[anc, ] <- P1 %*% liks[des1, ] * P2 %*% liks[des2, ] + v.l <- matexpo(Q * phy$edge.length[i]) %*% liks[des1, ] + v.r <- matexpo(Q * phy$edge.length[j]) %*% liks[des2, ] + v <- v.l * v.r + comp[anc] <- sum(v) + liks[anc, ] <- v/comp[anc] } - if (output.liks) return(liks[-(1:nb.tip), ]) - - 2 * log(sum(liks[nb.tip + 1, ])) + if (output.liks) return(liks[-TIPS, ]) + -2 * sum(log(comp[-TIPS])) } out <- nlminb(rep(ip, length.out = np), function(p) dev(p), lower = rep(0, np), upper = rep(Inf, np)) @@ -184,10 +187,8 @@ as the number of categories in `x'") else obj$se <- sqrt(diag(solve(h))) obj$index.matrix <- index.matrix if (CI) { - lik.anc <- dev(obj$rates, TRUE) - lik.anc <- lik.anc / rowSums(lik.anc) - colnames(lik.anc) <- lvls - obj$lik.anc <- lik.anc + obj$lik.anc <- dev(obj$rates, TRUE) + colnames(obj$lik.anc) <- lvls } } obj$call <- match.call() diff --git a/Thanks b/Thanks index 57e3cb2..157d4cc 100644 --- a/Thanks +++ b/Thanks @@ -5,7 +5,7 @@ Many users gave important feed-back with their encouragements, 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, Bret Larget, Nick Matzke, +Éric Durand, Olivier François, Rich FitzJohn Bret Larget, Nick Matzke, Elizabeth Purdom, Dan Rabosky, Klaus Schliep, Li-San Wang, Yan Wong, and Peter Wragg. Contact me if I forgot someone. -- 2.39.2