numerical_r.R

library(qfa)

numerical_r=function(obsdat,mkPlots=FALSE,span=0.5,nBrute=1000,cDiffDelta=0.0001){
	# Generate numerical (model-free) estimate for maximum intrinsic growth rate
	tims=obsdat$Expt.Time
	gdat=obsdat$Growth
	tmax=max(tims)

	# Smooth data, find slope as function of time and maximum slope
	gdat=log(gdat)
	tims=tims[!is.na(gdat)]
	gdat=gdat[!is.na(gdat)]
	la=NA
	try(la<-loapproxfun(tims,gdat,span=span),silent=TRUE)
	if(!is.function(la)) return(list(objective=0,maximum=NA))
	centralDiff=function(f,delta) return(function(x) (f(x+delta/2.0)-f(x-delta/2.0))/delta)
	slp=centralDiff(la,cDiffDelta)

	# Brute force optimization
	stimes=seq(min(tims),max(tims),length.out=nBrute)
	vals=la(stimes)
	slps=slp(stimes)
	opt=which.max(slps)
	maxval=list(objective=slps[opt],maximum=stimes[opt])
	maxslope=maxval$objective
	
	if(mkPlots){
		mainlab=paste("Max. intrinsic growth rate =",formatC(maxslope,3),"(estimated on log scale)")
		# Plot synthetic data, Loess approximation and estimated slope
		op=par(mar=c(5,4,4,5)+.1,mfrow=c(1,2))
		plot(NULL,xlab="Time (d)",ylab="",xlim=c(-0.1*tmax,tmax),ylim=range(gdat),main=mainlab,cex.main=0.75)
		abline(v=maxval$maximum,lwd=3,col="green",lty=2)
		slope=maxval$objective
		intercept=la(maxval$maximum)-slope*maxval$maximum
		abline(a=intercept,b=slope,lwd=3,col="green")
		points(tims,gdat)
		mtext("Log Cell density (AU)",side=2,line=3,col="red")
		curve(la,from=-0.1*tmax,to=tmax,add=TRUE,col="red",lwd=2,xlim=c(-0.1*tmax,tmax))
		par(new=TRUE)
		
		curve(slp(x),from=-0.1*tmax,to=tmax,col="blue",xlab="",ylab="",xaxt="n",yaxt="n",lwd=2,xlim=c(-0.1*tmax,tmax))
		axis(4)
		mtext("Slope of Log Cell Density",side=4,line=3,col="blue")

		mainlab=paste("Max. intrinsic growth rate =",formatC(maxslope,3),"(estimated on log scale)")
		# Plot synthetic data, Loess approximation and estimated slope
		op=par(mar=c(5,4,4,5)+.1)
		plot(NULL,xlab="Time (d)",ylab="",xlim=c(-0.1*tmax,tmax),ylim=range(obsdat$Growth),main=mainlab,cex.main=0.75)
		abline(v=maxval$maximum,lwd=3,col="green",lty=2)
		slope=exp(la(maxval$maximum))*maxval$objective
		intercept=exp(la(maxval$maximum))-slope*maxval$maximum
		abline(a=intercept,b=slope,lwd=3,col="green",untf=FALSE)
		points(obsdat$Expt.Time,obsdat$Growth)
		mtext("Cell density (AU)",side=2,line=3,col="red")
		curve(exp(la(x)),from=-0.1*tmax,to=tmax,add=TRUE,col="red",lwd=2,xlim=c(-0.1*tmax,tmax))
		par(new=TRUE)
		
		curve(x*slp(x),from=-0.1*tmax,to=tmax,col="blue",xlab="",ylab="",xaxt="n",yaxt="n",lwd=2,xlim=c(-0.1*tmax,tmax))
		axis(4)
		mtext("Slope of Cell Density",side=4,line=3,col="blue")

		par(op)
	}
	maxval$objective=maxslope
	return(maxval)
}

# Generate some synthetic data with simple error model
syndat=function(K,r,g,v,sigma,N=10,tmax=5.0,regularTime=FALSE,tsigma=0.01){
	tims=seq(0,tmax,length.out=N)
	if(!regularTime){tims=tims+rnorm(length(tims),0,tsigma)}
	gdat=Glogist(K,r,g,v,tims)+rnorm(length(tims),0,sigma)
	return(data.frame(Expt.Time=tims,Growth=gdat))
}

# Plot comparing linear and log-scale estimate
pdf("SyntheticGrowthCurves.pdf",width=10,height=5)
	K=0.3;r=16;g=9e-5;v=1.0;sigma=0.01;N=10;tmax=1.5;tsigma=0.01;regularTime=FALSE
	replicate(5,numerical_r(syndat(K,r,g,v,sigma,N=N,tmax=tmax,tsigma=tsigma,regularTime=FALSE),mkPlots=TRUE)$objective)
	K=0.3;r=16;g=9e-5;v=1.0;sigma=0.01;N=30;tmax=1.5;tsigma=0.01;regularTime=FALSE
	replicate(5,numerical_r(syndat(K,r,g,v,sigma,N=N,tmax=tmax,tsigma=tsigma,regularTime=FALSE),mkPlots=TRUE)$objective)
	K=0.3;r=8;g=9e-5;v=1.0;sigma=0.01;N=10;tmax=1.5;tsigma=0.01;regularTime=FALSE
	replicate(5,numerical_r(syndat(K,r,g,v,sigma,N=N,tmax=tmax,tsigma=tsigma,regularTime=FALSE),mkPlots=TRUE)$objective)
	K=0.3;r=8;g=9e-5;v=1.0;sigma=0.01;N=30;tmax=1.5;tsigma=0.01;regularTime=FALSE
	replicate(5,numerical_r(syndat(K,r,g,v,sigma,N=N,tmax=tmax,tsigma=tsigma,regularTime=FALSE),mkPlots=TRUE)$objective)
dev.off()

infer_r=function(K=0.3,r=16,g=9e-5,v=1.0,sigma=0.01,N=10,tmax=1.5,tsigma=0.01,regularTime=FALSE,Nsamples=1000){
	# Infer intrinsic growth rate using numeric_r and using logistic model
	numeric_r=replicate(Nsamples,numerical_r(syndat(K,r,g,v,sigma,N=N,tmax=tmax,regularTime=regularTime),mkPlots=FALSE)$objective)
	#qfa_r=replicate(Nsamples,makeFitness(data.frame(t(makefits(syndat(K,r,g,v,sigma,N=N,tmax=tmax,regularTime=regularTime),g))))$r)
	qfa_r=replicate(Nsamples,makefits(syndat(K,r,g,v,sigma,N=N,tmax=tmax,regularTime=regularTime),g,fixG=FALSE,glog=FALSE)[["r"]])
	summary(numeric_r)
	summary(qfa_r)

	mlab=paste("K=",K," r=",r," g=",g," v=",v,"\n sigma=",sigma," N=",N," tsigma=",tsigma," Nsamp=",Nsamples,sep="")
	plot(density(qfa_r),lwd=2,type="n",xlim=range(c(density(qfa_r)$x,density(numeric_r)$x)),xlab="Intrinsic growth rate (1/d)",ylab="Density",main=mlab,cex.main=0.75)
	abline(v=r,col="green",lwd=2)
	points(density(numeric_r),lwd=2,type="l")
	points(density(qfa_r),col="blue",type="l",lwd=2)
	legend("topright",col=c("black","blue","green"),legend=c("Numeric","QFA","True"),lwd=2)
}

pdf("Inferring_r.pdf")
op=par(mfrow=c(2,2))
	infer_r(K=0.3,r=16,g=9e-5,v=1.0,sigma=0.01,N=10,tmax=1.5,tsigma=0.01,regularTime=FALSE,Nsamples=100)
	infer_r(K=0.3,r=16,g=9e-5,v=1.0,sigma=0.01,N=10,tmax=1.5,tsigma=0.01,regularTime=FALSE,Nsamples=100)
	infer_r(K=0.3,r=8,g=9e-5,v=1.0,sigma=0.01,N=10,tmax=1.5,tsigma=0.01,regularTime=FALSE,Nsamples=100)
	infer_r(K=0.3,r=8,g=9e-5,v=1.0,sigma=0.01,N=10,tmax=1.5,tsigma=0.01,regularTime=FALSE,Nsamples=100)
par(op)
dev.off()