Joseph, Klaus,
This is great. I slightly edited the code and renamed the argument
'tips' to 'tip' (as in the original definition). I run on a bunch of
random trees and the results are exactly identical with the current
version of getMRCA. I'm changing the code in ape now.
Cheers,
Emmanuel
getMRCA <- function(phy, tip)
{
if (!inherits(phy, "phylo"))
stop('object "phy" is not of class "phylo"')
if (length(tip) < 2) return(NULL)
Ntip <- length(phy$tip.label)
rootnd <- Ntip + 1L
pars <- integer(phy$Nnode) # worst case assignment, usually far too
long
mrcaind <- 1L # index in pars of the mrca node. will return highest
one traversed by other lineages
tnd <- if (is.character(tip)) match(tip, phy$tip.label) else tip
## build a lookup table to get parents faster
pvec <- integer(Ntip + phy$Nnode)
pvec[phy$edge[, 2]] <- phy$edge[, 1]
## get entire lineage for first tip
nd <- tnd[1]
for (k in 1:phy$Nnode) {
nd <- pvec[nd]
pars[k] <- nd
if (nd == rootnd) break
}
pars <- pars[1:k] # delete the rest
## traverse lineages for remaining tips, stop if hit common ancestor
for (i in 2:length(tnd)) {
done <- FALSE
cnd <- tnd[i]
while (!done) {
cpar <- pvec[cnd] # get immediate parent
if (cpar %in% pars) {
## early exit! if the mrca of any set of taxa is the
root then we are done
if (cpar == rootnd) return(rootnd)
idx <- which(pars == cpar)
if (idx > mrcaind) mrcaind <- idx
done <- TRUE
} else cnd <- cpar # keep going!
}
}
pars[mrcaind]
}
Le 10/06/2017 à 00:28, Klaus Schliep a écrit :
Hi Joseph
I guess I figured out where some room for improvement was.
You grow the vector pars <- c(pars, nd)
I preallocated pars, but you could even just grow pars using pars[k] <- nd
without pre-allocation. This operation is now (with R >= 3.4) far less
costly as mentioned in the NEWS:
"Assigning to an element of a vector beyond the current length now
over-allocates by a small fraction. The new vector is marked internally as
growable, and the true length of the new vector is stored in the truelength
field. This makes building up a vector result by assigning to the next
element beyond the current length more efficient, though pre-allocating is
still preferred. The implementation is subject to change and not intended
to be used in packages at this time."
One of the many little improvements recently making R quite a bit faster
and easier.
So the worst case behavior was pretty bad
tree <- stree(100000, "right")
system.time(get_mrca_fast(tree, c(1,2)))
user system elapsed
13.424 0.012 13.440
mrca
[1] 199999
with the changes it is much better:
system.time(mrca <- get_mrca_fast(tree, c(1,2)))
user system elapsed
0.012 0.000 0.012
mrca
[1] 199999
And here is the code:
get_mrca_fast <- function (phy, tips) {
rootnd <- length(phy$tip.label) + 1;
pars <- integer(phy$Nnode); # worst case assignment, usually far too
long
mrcaind <- 1; # index in pars of the mrca node. will return
highest one traversed by other lineages
tnd <- NULL; # node ids of tips
if ("character" %in% class(tips)) {
tnd <- match(tips, phy$tip.label);
} else {
tnd <- tips;
}
# build a lookup table to get parents faster
pvec <- integer(max(phy$edge));
pvec[phy$edge[,2]] <- phy$edge[,1];
# get entire lineage for first tip
nd <- tnd[1];
for(k in 1:phy$Nnode){
nd <- pvec[nd]
pars[k] <- nd
if (nd == rootnd) break
}
pars <- pars[1:k] #delete the rest
# traverse lineages for remaining tips, stop if hit common ancestor
for (i in 2:length(tnd)) {
done <- FALSE;
cnd <- tnd[i];
while (!done) {
cpar <- pvec[cnd]; # get immediate parent
if (cpar %in% pars) {
if (cpar == rootnd) {
# early exit! if the mrca of any set of taxa is the
root then we are done
return(rootnd);
}
idx <- which(pars == cpar);
if (idx > mrcaind) {
mrcaind <- idx;
}
done <- TRUE;
} else {
# keep going!
cnd <- cpar;
}
}
}
return(pars[mrcaind]);
}
Have a nice weekend,
Klaus
On Fri, Jun 9, 2017 at 3:59 PM, Joseph W. Brown <josep...@umich.edu> wrote:
Liam,
I put the (updated) code up as a gist <https://gist.github.com/josephwb/
ad61fd29ed4fb06e712e23d67422c813>. Feel free to use/improve/whatever.
Emmanuel sees room for improvement but frankly taking 0.3 seconds to find
the MRCA of 5000 taxa on a million-tip tree is good enough for my purposes.
JWB
________________________________________
Joseph W. Brown
Post-doctoral Researcher, Smith Laboratory
University of Michigan
Department of Ecology & Evolutionary Biology
Room 2071, Kraus Natural Sciences Building
Ann Arbor MI 48109-1079
josep...@umich.edu
On 9 Jun, 2017, at 15:20, Liam J. Revell <liam.rev...@umb.edu> wrote:
On the other hand, phytools does have a function - the somewhat
imprecisely named fastMRCA - which can find the MRCA of just a pair of
species much faster than getMRCA (however still slower than or only about
as fast as Joseph & Klaus's solutions).
Liam J. Revell, Associate Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://blog.phytools.org
On 6/9/2017 5:22 AM, Liam J. Revell wrote:
Juan. findMRCA was written before getMRCA existed, but the latter was
faster so now it just calls getMRCA internally. All the best, Liam
Liam J. Revell, Associate Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://blog.phytools.org
On 6/9/2017 1:57 AM, Juan Antonio Balbuena wrote:
Package phytools includes a function, findMRCA, that is supposed to
work
very efficiently with large trees. you may wish to compare it with your
function.
Cheers
Juan
--
Dr. Juan A. Balbuena
Cavanilles Institute of Biodiversity and Evolutionary Biology
University of Valencia
http://www.uv.es/~balbuena <http://www.uv.es/%7Ebalbuena>
P.O. Box 22085
http://www.uv.es/cophylpaco
<http://www.uv.es/cavanilles/zoomarin/index.htm>
46071 Valencia, Spain
e-mail: j.a.balbu...@uv.es <mailto:j.a.balbu...@uv.es> tel. +34 963
543 658 fax +34 963 543 733
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*NOTE!*For shipments by EXPRESS COURIER use the following street
address:
C/ Catedrático José Beltrán 2, 46980 Paterna (Valencia), Spain.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
<https://www.avast.com/sig-email?utm_medium=email&utm_
source=link&utm_campaign=sig-email&utm_content=emailclient>
Libre de virus. www.avast.com
<https://www.avast.com/sig-email?utm_medium=email&utm_
source=link&utm_campaign=sig-email&utm_content=emailclient>
<#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at
http://www.mail-archive.com/r-sig-phylo@r-project.org/
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at
http://www.mail-archive.com/r-sig-phylo@r-project.org/
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-
sig-ph...@r-project.org/
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
