Nothing
sim_t_comp<-function(phylo,pars,root.value,Nsegments=1000,model="MC,DDexp,DDlin"){
#return error if non-ultrametric tree
if(!is.ultrametric(phylo)){stop("phylo object must be ultrametric")}
if(length(pars)!=2){stop("pars must be a vector with a value for sig2 and either S for MC model, r for DDexp model, or b for DDlin model")}
if(model=="MC"){sterm=pars[2]}
if(model=="DDlin"){slope=pars[2]}
if(model=="DDexp"){r=pars[2]}
sig2=pars[1]
paste(rep(LETTERS,each=26),LETTERS,sep="")->TWOLETTERS
paste(rep(TWOLETTERS,each=26),LETTERS,sep="")->THREELETTERS
ROOT<-root.value
##define a few useful variables
##nodeDist, the distance between each node and the root
##nodeDiff, the distance between each node and the previous node (times for integration)
nodeDist<-vector(mode = "numeric", length = phylo$Nnode)
root <- length(phylo$tip.label) + 1
heights<-nodeHeights(phylo)
for (i in 1:dim(phylo$edge)[1]){
nodeDist[[phylo$edge[i, 1] - length(phylo$tip.label)]] <- heights[i]
}
nodeDist<-c(nodeDist,max(heights))
nodeDiff<-diff(nodeDist)
###label the branches for each segment of tree to be integrated and identify the node at which the branch terminates
if(sum(nodeDiff<0)>0){ ##this loop renumbers the nodes if trees nodes are not placed in sequential order
node.order<-match(rank(heights[,1],ties.method="min"),seq(1, by = 2, len = phylo$Nnode))
node.order<-node.order+length(phylo$tip.label)
old.edge<-phylo$edge
phylo$edge[,1]<-node.order
for(j in 1:length(phylo$edge[,2])){
if(phylo$edge[j,2]>length(phylo$tip.label)){
#match number order in old edge
#lookup value in new edge
#replace with value
phylo$edge[j,2]<-phylo$edge[,1][match(phylo$edge[j,2],old.edge[,1])]
}
}
nodeDist<-vector()
for (i in 1:dim(phylo$edge)[1]){
nodeDist[[phylo$edge[i, 1] - length(phylo$tip.label)]] <- heights[i]
}
nodeDist<-c(nodeDist,max(heights))
nodeDiff<-diff(nodeDist)
}
mat<-matrix(nrow=0, ncol=3)
counter_three_letters <- 0
for(i in 1:phylo$Nnode){
other<-phylo$edge[phylo$edge[,1]==i+length(phylo$tip.label), 2]
for(b in other){
int<-matrix(ncol=3)
int[1]<-i+length(phylo$tip.label)
if(b>length(phylo$tip.label)){
counter_three_letters <- counter_three_letters + 1
int[2]<-paste(".",THREELETTERS[counter_three_letters],sep="")
int[3]<-b
} else {
int[2]<-phylo$tip.label[[b]]
int[3]<-0 ##NOTE :I am considering tips to be "0" here and use this below
}
mat<-rbind(mat,int)
}
}
##now come up with a list of branches that exist for each time segment (+1 at each branching in this version, which means tree can't have any polytomies or speciation events at exactly the same time)
nat<-list()
for(i in 1:length(nodeDiff)){
if(i==1){
nat[[i]]<-list(mat[mat[,1]==(length(phylo$tip.label)+i),2])} else {
IN<-vector()
P<-mat[as.numeric(mat[,1])<=(length(phylo$tip.label)+i),c(2,3)]
IN<-c(IN, P[P[,2]=="0",1],P[as.numeric(P[,2])>(length(phylo$tip.label)+i),1])
nat[[i]]<-list(IN)
}
}
for(i in 2:length(nodeDiff)){ ##THIS LOOP checks for an error
if(length(unlist(nat[[i]]))!=(length(unlist(nat[[i-1]]))+1)){
print(paste("ERROR at node",i+length(phylo$tip.label)))
}
}
if(is.na(match(model,c("MC","DDexp","DDlin")))){stop("model not specified correctly, must be 'MC','DDexp', or'DDlin'")}
if(model=="DDexp"){
if((sig2/(sig2*exp(r*length(phylo$tip.label))))>=1e3){warning("r parameter leads to sig2 at tips that is 1000x smaller than root value, consider changing")}
simvalueDDexp<-function(sig2,start.value,seglen,r,branches){
x<-start.value+rnorm(1,0,sqrt(sig2*exp(r*branches)*seglen))
return(x)
}
}
if(model=="DDlin"){
if(sig2+(slope*length(phylo$tip.label))<=0){warning("b parameter leads to sig2=0 by the tips, consider changing")}
simvalueDDlin<-function(sig2,start.value,seglen,slope,branches){
x<-start.value+rnorm(1,0,sqrt(max((sig2+(slope*branches)),0)*seglen))
return(x)
}
}
if(model=="MC"){
simvalueMC<-function(sig2,sterm,mu,start.value,seglen){
x<-start.value+sterm*(mu-start.value)*seglen+rnorm(1,0,sqrt(sig2*seglen))
return(x)
}
}
N<-Nsegments #number of smaller segments to divide a tree into
seglength<-tail(nodeDist,n=1)/N
masterbranch<-list()
masterseg<-list()
for(i in 1:phylo$Nnode){ ##for each node interval
if(i==1){mu<-ROOT} #initialize mu value at ROOT (change if another root value is desired)
branchespresent<-length(unlist(nat[[i]]))
masterbranch[[i]]<-as.list(rep(ROOT,branchespresent))
names(masterbranch[[i]])<-unlist(nat[[i]])
if(i>1){
#update starting values to be the ending values from last time
for(l in 1:length(masterbranch[[i]])){
if(!is.na(match(names(masterbranch[[i]][l]),names(masterbranch[[i-1]])))){ #if the name of the element matches branch present in previous iteration
masterbranch[[i]][l]<-sapply(masterbranch[[i-1]][match(names(masterbranch[[i]][l]),names(masterbranch[[i-1]]))],tail,n=1)
} else {
##look up branch from which descended
prev.branch<-mat[mat[,3]==mat[mat[,2]==names(masterbranch[[i]][l]),1],2]
##look up starting value
masterbranch[[i]][l]<-sapply(masterbranch[[i-1]][match(prev.branch,names(masterbranch[[i-1]]))],tail,n=1)
}
}
}
masterseg[[i]]<-as.list(rep(0,branchespresent))
names(masterseg[[i]])<-unlist(nat[[i]])
for(m in 1:length(masterseg[[i]])){
if(i==1){
if(seglength<nodeDiff[i]){
masterseg[[i]][[m]]<-c(masterseg[[i]][[m]],seq(seglength,nodeDiff[i],seglength))} else{
masterseg[[i]][[m]]<-c(masterseg[[i]][[m]],nodeDiff[i])
}} else{
masterseg[[i]][[m]]<-seq(nodeDist[i],nodeDist[i+1],seglength)
}
if(tail(masterseg[[i]][[m]],n=1)!=nodeDiff[i]){
masterseg[[i]][[m]]<-c(masterseg[[i]][[m]],nodeDist[i+1])
}
}
for(k in 1:length(diff(masterseg[[i]][[1]]))){ ##for each of the segments between intervals
for(j in 1:branchespresent){## for each of the branches present
##extract starting value
start.value<-tail(masterbranch[[i]][[j]],n=1) #replace 1 with matching function
#run simvalue function with mu (where to get value of mu and update?)
segsize<-diff(masterseg[[i]][[1]])[k]
mu.value<-tail(mu,n=1)
if(model=="DDexp" ){
sv<-simvalueDDexp(sig2,start.value,segsize,r,branchespresent)
}
if(model=="DDlin" ){
sv<-simvalueDDlin(sig2,start.value,segsize,slope,branchespresent)
}
if(model=="MC"){
sv<-simvalueMC(sig2,sterm,mu.value,start.value,segsize)}
masterbranch[[i]][[j]]<-c(masterbranch[[i]][[j]],sv)
}
#update mu, add to VECTOR (ie mu=c(mu,UPDATE))
mu=c(mu,mean(sapply(masterbranch[[i]],tail,n=1)))# the last element in each element of branchlist/branches present ##this assumes that the first element is updated for each i
}
}
return(sapply(masterbranch[[phylo$Nnode]],tail,n=1))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.