rm(list=ls())
library(TreeSim)

## Loading required package: ape

## Loading required package: geiger

# plot tree using path space
PlotTreeHistory<-function(phy=phy,pathdata=pathdata,main=main,colors=colors){
  ntips<-length(phy$tip.label)
  edge.number<-dim(phy$edge)[1]
  start<-rep(0,length(edge.number))
  end<-rep(0,length(edge.number))
  
  mtx<-(cbind(phy$edge,ceiling(phy$edge.length)+1,start,end))
  mtx<-mtx[dim(mtx)[1]:1,]
  colnames(mtx)<-c("anc","des","step","start","end")
  mtx<-as.data.frame(mtx)
  mtx$start[1:2]<-1
  mtx$end[1:2]<-mtx$step[1:2]
  for(edgeIndex in 3:dim(mtx)[1]){
    mtx$start[edgeIndex] <-  mtx$end[mtx$anc[edgeIndex]==mtx$des]
    mtx$end[edgeIndex] <-    mtx$end[mtx$anc[edgeIndex]==mtx$des]+mtx$step[edgeIndex]-1
  }
  xlim<-c(1, max(mtx$end) )
  unlistpath<-unlist(pathdata)
  ylim<-c(min(unlistpath),max(unlistpath))
  plot(NA,xlim=xlim,ylim=ylim,ylab="Trait value", xlab="Times Steps",main=main)
  points(1, pathdata[[1]][1] ,lwd=5,col= 1)
  for(edgeIndex in 1:edge.number){
    #edgeIndex<-7
    path<-unlist(pathdata[edgeIndex])
    l1<-length(path)
    l2<-length(mtx$start[edgeIndex]:mtx$end[edgeIndex])
    if(l1<l2){   
      addupath<-l2-l1
      path<-c(path,  rep(path[l1],l2-l1))
    }
    points(mtx$start[edgeIndex]:mtx$end[edgeIndex], path,type="l",lwd=2,col=colors[edgeIndex])
    points(mtx$end[edgeIndex]:mtx$end[edgeIndex], path[length(path)],cex=2,lwd=1,col=colors[edgeIndex])
  }
}#end of plot history

#Brownian bridge
bb<-function(sigma,a=0,b=0,t0=0,T=1,N=100){
  if(T<=t0){stop("wrong times")}
  dt <- (T-t0)/N
  t<- seq(t0,T,length=N+1)
  X<-c(0,sigma*cumsum(rnorm(N)*sqrt(dt)))
  bbpath<-a+X-(t-t0)/(T-t0)*(X[N+1]-b+a)
  return(bbpath)
}

bbtree<-function(sigma,phy=phy,tips=tips){
  anclist<-ace(x=tips,phy=phy)
  fullnodedata<-c(tips,anclist$ace)
  anc<-phy$edge[,1]
  des<-phy$edge[,2]
  ntips<-length(phy$tip.label)
  edge.number<-dim(phy$edge)[1]
  edge.length<-phy$edge.length
  pathdata<-list()
  for(edgeIndex in edge.number:1){
    brlen<-edge.length[edgeIndex]
    N<-ceiling(brlen)
    b<-fullnodedata[des[edgeIndex]]
    a<-fullnodedata[anc[edgeIndex]]
    assign(paste("path",edgeIndex,sep=""),bb(sigma,a=a,b=b,t0=0,T=brlen,N=N))
    pathdata<-c(pathdata, list(get(paste("path",edgeIndex,sep=""))) )
  }
  return(pathdata)
}

library(ape)
set.seed(443)
#tree<-"((z:1,(y:0.42,x:0.42):0.58):0);"
#phy<-read.newick(text=tree)
phy<-sim.bdtree(b=1, d=0, stop=c("taxa", "time"), n=3, t=4, seed=0, extinct=TRUE)
phy<-reorder(phy,"postorder")
colors<-c("red","orange","blue","black")
phy$edge.length

## [1] 0.5899223 0.5899223 0.7257700 0.1358477

min.length<-20
while(min(phy$edge.length)<min.length){
  phy$edge.length<-1.001*phy$edge.length
}
phy$edge.length

## [1]  86.89841  86.89841 106.90943  20.01102

phy$tip.label<-c("   x","   y","   z")

C<-vcv(phy)
X<-c(23.4,26.7,24.5) 
sigma<-0.1

par(mfrow=c(1,2))
plot(phy,edge.width=5,cex=2,edge.color=rev(colors))
title("A 3 Taxa Time-calibrated Tree ")
axisPhylo(1, root.time = NULL, backward = F)

bbpathdata<-bbtree(sigma,phy=phy,tips=X)
PlotTreeHistory(phy=phy,pathdata=bbpathdata,main="A Gaussian Process Trait Evolution",colors=colors)

Star Tree

set.seed(56)
phy <- compute.brlen(stree(3))
phy$tip.label<-c("   x","   y","   z")
phy$edge.length<-rep(100,3)
colors<-c("red","blue","black")
par(mfrow=c(1,3))
plot(phy,edge.width=5,cex=2,edge.color=rev(colors))
title("A 3 Taxa Star Tree ")
axisPhylo(1, root.time = NULL, backward = F)
X<-c(23.4,26.7,24.5) 
sigma<-0.01
T<-1e3
N<-1000
bb1<-bb(sigma=sigma,a=mean(X),b=X[1],t0=0,T=T,N=N)
bb2<-bb(sigma=sigma,a=mean(X),b=X[2],t0=0,T=T,N=N)
bb3<-bb(sigma=sigma,a=mean(X),b=X[3],t0=0,T=T,N=N)

plot(bb1,type="l",col="black",lwd=2,ylim=range(c(bb1,bb2,bb3)),xlab="Time",ylab="Trait Value", main="Trajectories",cex.lab=1.5)
points(1:1001,bb2,type="l",col="red",lwd=2)
points(1:1001,bb3,type="l",col="blue",lwd=2)

id.path.data.x<-list(bb1,bb2,bb3)
bb1.y<-bb(sigma=sigma,a=mean(X),b=X[1],t0=0,T=T,N=N)
bb2.y<-bb(sigma=sigma,a=mean(X),b=X[2],t0=0,T=T,N=N)
bb3.y<-bb(sigma=sigma,a=mean(X),b=X[3],t0=0,T=T,N=N)
id.path.data.y<-list(bb1.y,bb2.y,bb3.y)

# ID 2D tree plot 
xlims<-range(unlist(id.path.data.x))
ylims<-range(unlist(id.path.data.y))

plot(id.path.data.y[[1]]~id.path.data.x[[1]],col=colors[1],xlim=xlims,ylim=ylims,type="l",xlab="Trait 1",ylab="Trait 2",lwd=2,main="",cex.lab=1.5)
text(24,26.5, "ID 2D",cex=2)

for(Index in 2:3){
  points(id.path.data.y[[Index]]~id.path.data.x[[Index]],col=colors[Index],type="l",lwd=2)
}

A 3 taxa tree and the Gaussian process replicate

Star Tree