Thursday, May 2, 2013

Bug fix in fastAnc

A couple of days ago a phytools user reported a bug in geomorph (which reverse depends on phytools) that seemed to be due to a problem with fastAnc, which is used internally.

Well, I finally got around to looking into it. It turns out that there was a bug in fastAnc, but I'd overlooked it for good reason - it only occurs under the somewhat idiosyncratic circumstances of when is.binary.tree(tree)=FALSE and the user-specified option vars=FALSE. This is because fastAnc works by re-rooting the tree at all internal nodes of a binary tree and computing the PIC ancestral state & variance. If the input tree is not binary, then it uses multi2di, but then has to back-translate to the original tree which it does using phytools matchNodes.

The problem arose because:

if(!is.binary.tree(tree)){
  ancNames<-matchNodes(tree,btree)
  anc<-anc[as.character(ancNames[,2])]
  names(anc)<-ancNames[,1]
  if(vars) v[as.character(ancNames[,2])]
  names(v)<-ancNames[,1]
}
should have been:
if(!is.binary.tree(tree)){
  ancNames<-matchNodes(tree,btree)
  anc<-anc[as.character(ancNames[,2])]
  names(anc)<-ancNames[,1]
  if(vars||CI){
    v[as.character(ancNames[,2])]
    names(v)<-ancNames[,1]
  }
}
in the code that is executed to back-translate nodes betweeen trees. Basically, if is.binary.tree(tree)=FALSE and vars=FALSE (but under no other circumstances), fastAnc will try to assign names to a vector of zero length. Oops.

The fixed source code for this very simple function is here, but I also posted a new phytools build (phytools 0.2-55).

3 comments:

  1. Thank you for the update! I had this really cool car in high school where as soon as I fixed one thing something else broke, which may be happening here. I get a different error now.

    p1 <- physignal(phy=MSB2brl, Pat$coords, iter=100)
    Error: internal error -3 in R_decompress1

    ReplyDelete
    Replies
    1. I like the metaphor!

      However, it looks like this is a problem with geomorph, not phytools.

      Basically, it looks like physignal assumes the tree is binary without checking.

      Your options are to make your tree fully bifurcating using multi2di (not recommended because this will be make ancestral state estimation unnecessarily slow); or modify the source of physignal replacing every instance of nrow(x) - 1 with phy$Nnode.

      Let us know if that works - and perhaps report to the maintainer of geomorph if it does.

      All the best, Liam

      Delete
  2. That did it! Thanks and I'll pass along the note to Dean Adams

    ReplyDelete