New version of fastBM() for multiple simulations

In a previous post I described a fast Brownian simulation function, fastBM(), and I showed that it was indeed quite fast. I have since added some additional options to this function - for instance, simulating with a trend or simulating within hard boundaries.

One teeny-tiny problem though. Although this function is much faster than sim.char() in {geiger} for a single simulation, e.g.:

> tree<-rtree(1000)
> system.time(x<-sim.char(phy=tree,as.matrix(1)))
user system elapsed
80.87 0.00 80.99

compared to:

> system.time(x<-fastBM(tree))
user system elapsed
0.18 0.00 0.19

[using v0.2 of fastBM() on an Intel Core i5 @ 3.20GHz], it can actually be much slower for multiple simulations. For instance, sim.char() takes basically the same time to run 100 simulations as it does to run 1:

> system.time(X<-sim.char(phy=tree,as.matrix(1),nsim=100))
user system elapsed
80.76 0.02 80.87

whereas the only way to run multiple simulations on the same tree using fastBM() has been to loop the function many times, e.g.

> X<-matrix(NA,length(tree$tip.label),100,
   dimnames=list(tree$tip.label))
> system.time(X<-apply(X,2,function(x) x<-fastBM(tree)))
user system elapsed
19.74 0.00 19.74

Although fastBM() still does better in this example, it has gone from being 400 times faster - to a mere 4 times faster! That is because when we just loop calls to fastBM(), computation time increases in direct proportion to the number of simulations (in this case ~0.2s x 100 = 20s). By contrast, sim.char() only has one very hard calculation to do - the eigen-decomposition of the matrix given by vcv(tree). Once this calculation is complete, the function needs very little additional time to run more simulations.

It is not difficult to see where this is going. For any 400 simulations sim.char() and fastBM() take the same amount of time, and for any more than that sim.char() becomes faster.

Today, I decided to add the capacity to run multiple simulations - without looping over the tree many times - into fastBM(). The new version (tagged v0.4) is available on my R-phylogenetics page (direct link to code here). This was really pretty easy to do. The only tricky parts involved evaluating boundary violations across all the layers of the simulated data array, and then reflecting if node values were outside the bound. To do this, I first created an internal function reflect():

reflect<-function(yy,bounds){
    while(yy < bounds[1]||yy > bounds[2]){
        if(yy < bounds[1]) yy<-2*bounds[1]-yy
        if(yy > bounds[2]) yy<-2*bounds[2]-yy
    }
    return(yy)
}

and then I called this function using apply(), i.e.:

y[i,2,]<-apply(as.matrix(y[i,2,]),1,function(yy)
   reflect(yy,bounds))

I also created a logical variable, no.bounds, to avoid calling the internal function reflect() if the bounds are trivial (in other words, for bounds=c(-Inf,Inf), the default). I.e.:

if(!no.bounds) y[i,2,]<-apply(as.matrix(y[i,2,]),1,function(yy)
   reflect(yy,bounds))

With these problems solved, we can now evaluate the performance of our new version of fastBM():

> system.time(X<-fastBM(tree,nsim=100))
user system elapsed
0.28 0.00 0.28

The instinct to limit our calls to reflect() was indeed a good idea, because simulating with bounds is quite a bit slower than this, e.g.:

> system.time(X<-fastBM(tree,nsim=100,bounds=c(-2,4)))
user system elapsed
2.91 0.00 2.92

Neat.

Since I wrote this code this afternoon, please let me know if you are able to get it to run (either post a comment - or email me directly).