improve fastBM

R/fastBM.R

@@ -27,13 +27,14 @@ fastBM<-function(tree,a=0,mu=0,sig2=1,bounds=c(-Inf,Inf),internal=FALSE,nsim=1,.
 		cat("Warning: OU with bounds not permitted. Bounds will be ignored.\n")
 	## if BM
 	if(is.null(alpha)) x<-simBM(tree,a,mu,sig2,bounds,internal,nsim)
-	else x<-if(nsim==1) simOU(tree,alpha,sig2,theta,a,internal) else replicate(nsim,simOU(tree,alpha,sig2,theta,a,internal))
+	else x <- simOU(tree, alpha, sig2, theta, a, internal, nsim) 
 	x
 }
 
 ## internal function does BM simulation
 ## written by Liam J. Revell 2011, 2013
 simBM<-function(tree,a,mu,sig2,bounds,internal,nsim){
+    tree <- reorder(tree)
 	if(bounds[2]<bounds[1]){
 		warning("bounds[2] must be > bounds[1]. Simulating without bounds.")
 		bounds<-c(-Inf,Inf)
@@ -56,25 +57,38 @@ simBM<-function(tree,a,mu,sig2,bounds,internal,nsim){
 		}
 		return(yy)
 	}
-	# how many species?
-	n<-length(tree$tip)
-	# first simulate changes along each branch
-	x<-matrix(data=rnorm(n=length(tree$edge.length)*nsim,mean=rep(mu*tree$edge.length,nsim),sd=rep(sqrt(sig2*tree$edge.length),nsim)),length(tree$edge.length),nsim)
-	# now add them up
-	y<-array(0,dim=c(nrow(tree$edge),ncol(tree$edge),nsim))
-	for(i in 1:nrow(x)){
-		if(tree$edge[i,1]==(n+1))
-			y[i,1,]<-a
-		else
-			y[i,1,]<-y[match(tree$edge[i,1],tree$edge[,2]),2,]
+# the following code is replaced and makes use of the ordering of the tree	
+#	# how many species?
+#	n<-length(tree$tip)
+#	# first simulate changes along each branch
+#	x<-matrix(data=rnorm(n=length(tree$edge.length)*nsim,mean=rep(mu*tree$edge.length,nsim),sd=rep(sqrt(sig2*tree$edge.length),nsim)),length(tree$edge.length),nsim)
+#	# now add them up
+#	y<-array(0,dim=c(nrow(tree$edge),ncol(tree$edge),nsim))
+#	for(i in 1:nrow(x)){
+#		if(tree$edge[i,1]==(n+1))
+#			y[i,1,]<-a
+#		else
+#			y[i,1,]<-y[match(tree$edge[i,1],tree$edge[,2]),2,]
+#
+#		y[i,2,]<-y[i,1,]+x[i,]
+#		if(!no.bounds) y[i,2,]<-apply(as.matrix(y[i,2,]),1,function(yy) reflect(yy,bounds))
+#	}
+#	rm(x); x<-matrix(data=rbind(y[1,1,],as.matrix(y[,2,])),length(tree$edge.length)+1,nsim)
+#	rownames(x)<-c(n+1,tree$edge[,2])
+#	x<-as.matrix(x[as.character(1:(n+tree$Nnode)),])
+#	rownames(x)[1:n]<-tree$tip.label
 
-		y[i,2,]<-y[i,1,]+x[i,]
-		if(!no.bounds) y[i,2,]<-apply(as.matrix(y[i,2,]),1,function(yy) reflect(yy,bounds))
+## additions by Klaus Schliep Dec 2015	
+	nTips <- Ntip(tree)
+	x<-matrix(0,max(tree$edge),nsim)
+	x[tree$edge[1,1], ] <- a
+	for(i in 1:nrow(tree$edge)){
+	    x[tree$edge[i,2],] <- x[tree$edge[i,1],] + rnorm(nsim, mean= mu * tree$edge.length[i], sd=sqrt(sig2*tree$edge.length[i]))
+	    if(!no.bounds) x[tree$edge[i,2],]<-apply(as.matrix(x[tree$edge[i,2],]),1,function(yy) reflect(yy,bounds))
 	}
-	rm(x); x<-matrix(data=rbind(y[1,1,],as.matrix(y[,2,])),length(tree$edge.length)+1,nsim)
-	rownames(x)<-c(n+1,tree$edge[,2])
-	x<-as.matrix(x[as.character(1:(n+tree$Nnode)),])
-	rownames(x)[1:n]<-tree$tip.label
+	rownames(x) <- c(tree$tip.label, (nTips+1L):max(tree$edge))
+# 	
+
 	# return simulated data
 	if(internal==TRUE)
 		return(x[1:nrow(x),]) # include internal nodes
@@ -84,22 +98,37 @@ simBM<-function(tree,a,mu,sig2,bounds,internal,nsim){
 
 ## internal function does BM simulation
 ## written by Liam J. Revell 2013
-simOU<-function(tree,alpha,sig2,theta,a0,internal){
+## additions by Klaus Schliep Dec 2015
+
+# simOU<-function(tree,alpha,sig2,theta,a0,internal){
+simOU <- function(tree, alpha, sig2, theta, a0, internal, nsim){    
 	tree<-reorder(tree,"cladewise")
-	X<-matrix(0,nrow(tree$edge),ncol(tree$edge))
-	root<-length(tree$tip.label)+1
-	X[which(tree$edge[,1]==root),1]<-a0
-	for(i in 1:nrow(X)){
-		t<-tree$edge.length[i]
-		s2<-sig2*(1-exp(-2*alpha*t))/(2*alpha)
-		X[i,2]<-exp(-alpha*t)*X[i,1]+(1-exp(-alpha*t))*theta+rnorm(n=1,sd=sqrt(s2))
-		ii<-which(tree$edge[,1]==tree$edge[i,2])
-		if(length(ii)>0) X[ii,1]<-X[i,2]
+#	X<-matrix(0,nrow(tree$edge),ncol(tree$edge))
+#	root<-length(tree$tip.label)+1
+#	X[which(tree$edge[,1]==root),1]<-a0
+#	for(i in 1:nrow(X)){
+#		t<-tree$edge.length[i]
+#		s2<-sig2*(1-exp(-2*alpha*t))/(2*alpha)
+#		X[i,2]<-exp(-alpha*t)*X[i,1]+(1-exp(-alpha*t))*theta+rnorm(n=1,sd=sqrt(s2))
+#		ii<-which(tree$edge[,1]==tree$edge[i,2])
+#		if(length(ii)>0) X[ii,1]<-X[i,2]
+#	}
+#	x<-sapply(1:max(tree$edge),function(x,y,tree) y[which(tree$edge==x)[1]],y=X,tree=tree)
+#	x<-setNames(x,c(tree$tip.label,1:tree$Nnode+length(tree$tip.label)))
+
+## new 		
+	nTips <- Ntip(tree)
+	x<-matrix(0,max(tree$edge),nsim)
+	x[tree$edge[1,1], ]<-a0
+	for(i in 1:nrow(tree$edge)){
+	    t<-tree$edge.length[i]
+	    s2<-sig2*(1-exp(-2*alpha*t))/(2*alpha)
+	    x[tree$edge[i,2],] <- exp(-alpha*t) * x[tree$edge[i,1],] + (1-exp(-alpha*t))*theta + rnorm(n=nsim, sd=sqrt(s2))
 	}
-	x<-sapply(1:max(tree$edge),function(x,y,tree) y[which(tree$edge==x)[1]],y=X,tree=tree)
-	x<-setNames(x,c(tree$tip.label,1:tree$Nnode+length(tree$tip.label)))
+	rownames(x) <- c(tree$tip.label, (nTips+1L):max(tree$edge))
+
 	if(internal==TRUE)
-		return(x) # include internal nodes
+		return(x[1:nrow(x),]) # include internal nodes
 	else
-		return(x[1:length(tree$tip.label)]) # tip nodes only
+		return(x[1:length(tree$tip.label),]) # tip nodes only
 }