# function animates branching random diffusion
# written by Liam Revell 2011

branching.diffusion<-function(sig2=1,b=0.0023,time.stop=1000,ylim=NULL,smooth=TRUE,pause=0.02,record=NULL,path=NULL){
	if(!is.null(record)) if(!require(animation)) stop("must first install 'animation' package")
	N<-1; Y<-matrix(0,1,N)
	if(is.null(ylim)) ylim<-c(-2*sqrt(sig2*time.stop),2*sqrt(sig2*time.stop))
	par(bg="white")
	plot(0,0,xlim=c(0,time.stop),ylim=ylim,xlab="time",ylab="phenotype",main="branching diffusion")	
	if(!is.null(record)){
		if(is.null(path)) path="C:/Program Files/ffmpeg/bin/ffmpeg.exe"
		ani.options(interval=0.01,ffmpeg=path,outdir=getwd())
		ani.record(reset=TRUE)
	}
	for(i in 2:time.stop){
		time<-1:i
		Y<-rbind(Y,Y[i-1,])
		for(j in 1:N){
			if(b>runif(n=1)){
				Y<-cbind(Y,Y[,j])
				N<-N+1
				Y[i,N]<-Y[i,j]+rnorm(n=1,sd=sqrt(sig2))
			}
			Y[i,j]<-Y[i,j]+rnorm(n=1,sd=sqrt(sig2))
		}
		if(smooth==FALSE) for(j in 1:N) lines(c(time[i-1],time[i]),c(Y[i-1,j],Y[i,j]))
		if(smooth==TRUE){
			plot(0,0,xlim=c(0,time.stop),ylim=ylim,xlab="time",ylab="phenotype",main="branching diffusion")
			for(j in 1:N) lines(time,Y[,j])
		}
		if(!is.null(record)) ani.record()
		Sys.sleep(pause)
	}
	if(!is.null(record)) saveVideo(ani.replay(),video.name=record,other.opts="-b 300k",clean=TRUE)
}