Faster tree MRCA calculator in R (alt to APE's getMRCA)

get_mrca_fast.R

## written by:
# Joseph W. Brown (josephwb, @j0sephwb)
# Klaus Schliep (KlausVigo, @Klaus_Schlp)
get_mrca_fast <- 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
    tnd <- if (is.character(tip)) match(tip, phy$tip.label) else tip
    
    # keep track of traversed nodes
    done_v <- logical(max(phy$edge));
    
    ## 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
    mrcind <- integer(max(pars));
    mrcind[pars] <- 1:k;
    mrcand <- pars[1]; # initialized MRCA node

    ## traverse lineages for remaining tips, stop if hit common ancestor
    for (i in 2:length(tnd)) {
        cnd <- tnd[i];
        done <- done_v[cnd];
        while (!done) {
            done_v[cnd] <- TRUE;
            cpar <- pvec[cnd]; # get immediate parent
            done <- done_v[cpar]; # early exit if TRUE (already traversed this node, tho not in pars)
            if (cpar %in% pars) {
                if (cpar == rootnd) {
                    return(rootnd); # early exit (can't get shallower that the root)
                }
                if (mrcind[cpar] > mrcind[mrcand]) {
                    mrcand <- cpar;
                }
                done_v[cpar] <- TRUE;
                done <- TRUE;
            }
            cnd <- cpar; # keep going!
        }
    }
    return(mrcand);
}

This is now used in ape.