Faster version of getStates, but why is it faster...?

I was tinkering with the function describe.simmap to try and speed it up, when I discovered that the very simple function getStates, which does nothing more than get the states on a mapped tree in memory for the nodes or tips, had run-time using lapply to iterate over the trees in an object of class "multiPhylo", that seems to rise non-linearly with the number of trees.

> ## simulate tree & data
> tree<-pbtree(n=200,scale=1)
> Q<-matrix(c(-1,1,1,-1),2,2)
> colnames(Q)<-rownames(Q)<-letters[1:2]
> trees<-sim.history(tree,Q,nsim=200)
> # now lets get the states for all trees at all nodes
> # using sapply
> system.time(XX1<-getStates(trees[[1]]))
  user  system elapsed
  0.02    0.00    0.02
> system.time(XX10<-sapply(trees[1:10],getStates))
  user  system elapsed
  0.05    0.00    0.04
> system.time(XX50<-sapply(trees[1:50],getStates))
  user  system elapsed
  0.28    0.00    0.29
> system.time(XX100<-sapply(trees[1:100],getStates))
  user  system elapsed
  1.16    0.00    1.16
> system.time(XX200<-sapply(trees,getStates))
  user  system elapsed
    6.6    0.0    6.6

Hmmm. What's going on here?

What I discovered (somehow - I'm not sure why I tried this) is that if I first split my list of 200 trees into, say, 20 lists of 10 trees; and then I ran sapply(...,getStates) on each of these lists; then recombined the results using cbind, this is much faster. So, for instance:

> g<-function(trees){
 ff<-as.factor(ceiling(1:length(trees)/10))
 aa<-lapply(split(trees,ff),function(x)   sapply(x,getStates))
 y<-if(length(trees)>10) aa[[1]] else aa
 for(i in 2:length(aa)) y<-cbind(y,aa[[i]])
 y
}
> # now run it
> system.time(YY200<-g(trees))
  user  system elapsed
  0.97    0.00    0.97
> # check all equal
> dim(YY200)
[1] 199 200
> dim(XX200)
[1] 199 200
> all(XX200==YY200)
[1] TRUE

Seriously - what's going on here? (I have some ideas - but I'm not sure. Feedback welcome.)

Here's a new version of getStates with this hack implemented internally:

# function to get node states from simmap style trees
# written by Liam J. Revell 2013
getStates<-function(tree,type=c("nodes","tips")){
  type<-type[1]
  if(class(tree)=="multiPhylo"){
    ff<-as.factor(ceiling(1:length(tree)/10))
    aa<-lapply(split(tree,ff),function(x)
     sapply(x,getStates))
    y<-if(length(tree)>10) aa[[1]] else aa
    for(i in 2:length(aa)) y<-cbind(y,aa[[i]])
  } else if(class(tree)=="phylo"){
    if(type=="nodes"){
      y<-setNames(sapply(tree$maps,function(x)
       names(x)[1]),tree$edge[,1])
      y<-y[as.character(length(tree$tip)+1:tree$Nnode)]
    } else if(type=="tips"){
      y<-setNames(sapply(tree$maps,function(x)
       names(x)[length(x)]),tree$edge[,2])
      y<-setNames(y[as.character(1:length(tree$tip))],
       tree$tip)
    }
  } else stop("tree should be an object of class 'phylo' or 'multiPhylo'")
  return(y)
}