From 9ba3a61fec6c41318c6cdf64b8eba2e4c331e571 Mon Sep 17 00:00:00 2001 From: klash Date: Thu, 10 Dec 2015 19:43:10 -0500 Subject: [PATCH 1/2] improved fastBM --- R/fastBM.R | 89 ++++++++++++++++++++++++++++++++++++------------------ 1 file changed, 59 insertions(+), 30 deletions(-) diff --git a/R/fastBM.R b/R/fastBM.R index 3f18d8c4..e7b3eadd 100644 --- a/R/fastBM.R +++ b/R/fastBM.R @@ -25,13 +25,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]. Simulating without bounds.") bounds<-c(-Inf,Inf) @@ -54,25 +55,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 @@ -82,20 +96,35 @@ 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 else From eea678e1737c9291b0b9468ab2716e3a7e987fae Mon Sep 17 00:00:00 2001 From: klash Date: Thu, 10 Dec 2015 19:58:31 -0500 Subject: [PATCH 2/2] names were missing --- R/fastBM.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/fastBM.R b/R/fastBM.R index e7b3eadd..e32e5723 100644 --- a/R/fastBM.R +++ b/R/fastBM.R @@ -126,7 +126,7 @@ simOU <- function(tree, alpha, sig2, theta, a0, internal, nsim){ 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 }