ccoptim<-function(datloc,outloc=getwd(),graph=FALSE) {
	require(ucminf)
	
	tdat<-t(read.table(datloc,header=FALSE,sep='\t'))
	dat<-array(dim=c((length(t(tdat[,1]))-1),length(t(tdat[1,]))))
	colnames(dat)<-tdat[1,]
	tdat<-tdat[-1,]
	for(i in 1:length(tdat[,1])) dat[i,]<-as.numeric(tdat[i,])
	
	# Covariance function
	covar<-function(idat,est,boolean,boolean2) {# Inputs are individual gene data and parameter estimates
		pft1<-function(idat) {					# Partial derivative for theta 1
			return(1)
			}
		pft2<-function(idat) {					# Partial derivative for theta 2
			return(idat)
			}
		pft3<-function(idat,est) {				# Patrial derivative for theta 3
			return(1/sqrt(2*pi)*integrate(
			function(z) {
				cos((2*pi*idat)/(est[5]*exp(est[6]*z))+est[4])*exp(-.5*z^2)
				},
			-10,10,stop.on.error=FALSE,subdivisions=100000)$val)
			}
		pft4<-function(idat,est) {				# Partial derivative for theta 4
			return(-1/sqrt(2*pi)*integrate(
			function(z) {
				sin((2*pi*idat)/(est[5]*exp(est[6]*z))+est[4])*exp(-.5*z^2)
				},
			-10,10,stop.on.error=FALSE,subdivisions=100000)$val)
			}
		pft5<-function(idat,est) {				# Partial derivative for theta 5
			return(1/sqrt(2*pi)*(integrate(
			function(z) {
				(sin((2*pi*idat)/(est[5]*exp(est[6]*z))+est[4])*exp(-.5*z^2)*2*pi*idat)/(est[5]^2*exp(est[6]*z))
				},
			-10,10,stop.on.error=FALSE,subdivisions=100000)$val))
			}
		pft6<-function(idat,est) {				# Partial derivative for theta 6
			return(1/sqrt(2*pi)*(integrate(
			function(z) {
				(sin((2*pi*idat)/(est[5]*exp(est[6]*z))+est[4])*exp(-.5*z^2)*2*pi*idat*z)/(est[5]*exp(est[6]*z))
				},
			-10,10,stop.on.error=FALSE,subdivisions=100000)$val))
			}
		
		A<-array(dim=c(length(idat[,1]),6))		# Matrix A of partial derivatives evaluated at timepoints t for derived estimates
		A[,1]<-apply(as.array(idat[,1]),1,pft1)
		A[,2]<-apply(as.array(idat[,1]),1,pft2)
		A[,3]<-apply(as.array(idat[,1]),1,pft3,est)
		A[,4]<-apply(as.array(idat[,1]),1,pft4,est)
		A[,5]<-apply(as.array(idat[,1]),1,pft5,est)
		A[,6]<-apply(as.array(idat[,1]),1,pft6,est)
		
		# Sigma^2 calculation
		sigmasq<-sum((idat[,2]-apply(as.array(idat[,1]),1,
		function(t) {
			est[1]+est[2]*t+est[3]/sqrt(2*pi)*(integrate(
			function(z)	cos(2*pi*t/(est[5]*exp(est[6]*z))+est[4])*exp(-5*z^2),-10,10,subdivisions=10000)$val)
			}
		))^2)/(length(idat[,1])-6)
		ifelse(boolean,
			ifelse(boolean2,return(sigmasq*solve(t(A)%*%A)),return(1/sigmasq*(t(A)%*%A))),
			ifelse(boolean2,return(sigmasq*solve(t(A[,-2])%*%A[,-2])),return(1/sigmasq*(t(A[,-2])%*%A[,-2]))))
		}
		
	# Least absolute error random period model function
	f<-function(est,idat) {
		sum<-0
		for(i in 1:length(idat[,1])) {
			eval<-abs(idat[i,2]-(est[1]+est[2]*idat[i,1]+est[3]/sqrt(2*pi)*(integrate(
			function(z) {
				ifelse((est[5]*exp(est[6]*z))==0,return(0),return(cos(2*pi*idat[i,1]/(est[5]*exp(est[6]*z))+est[4])*exp(-.5*z^2))) 
				}
			,-10,10,subdivisions=10000,stop.on.error=FALSE)$val)))
			sum<-sum+eval
			i=i+1
			}
		return(sum)	
		}
	
	curdir=getwd()				# Get current directory
	if(outloc!=curdir) setwd(outloc)
		
	# Output function
	out<-function(idat,est,A,B,C,D,f,outloc) {
		cat(
		dimnames(idat)[[2]][2],"\t",	# Gene name
		est[1],"\t",			# Six estimated parameters
		est[2],"\t",
		est[3],"\t",
		est[4],"\t",
		est[5],"\t",
		est[6],"\t",
		f(est,idat),"\n",		# Function evaluated at estimated parameters
		file='output.txt',		# output file name
		append=TRUE				# appends current output.txt file to show results procedurally
		)
		cat(
		dimnames(idat)[[2]][2],"\t",	# Gene name
		A[1,1],"\t",			# 15 relevant covariance values with the first 6 being the diagonals
		A[2,2],"\t",
		A[3,3],"\t",
		A[4,4],"\t",
		A[5,5],"\t",
		A[6,6],"\t",
		A[1,2],"\t",
		A[1,3],"\t",
		A[1,4],"\t",
		A[1,5],"\t",
		A[1,6],"\t",
		A[2,3],"\t",
		A[2,4],"\t",
		A[2,5],"\t",
		A[2,6],"\t",
		A[3,4],"\t",
		A[3,5],"\t",
		A[3,6],"\t",
		A[4,5],"\t",
		A[4,6],"\t",
		A[5,6],"\t",
		kappa(A),"\n",
		file='covarm.txt',
		append=TRUE
		)
		cat(
		dimnames(idat)[[2]][2],"\t",	# Gene name
		B[1,1],"\t",			# 15 relevant covariance values with the first 6 being the diagonals
		B[2,2],"\t",
		B[3,3],"\t",
		B[4,4],"\t",
		B[5,5],"\t",
		B[1,2],"\t",
		B[1,3],"\t",
		B[1,4],"\t",
		B[1,5],"\t",
		B[2,3],"\t",
		B[2,4],"\t",
		B[2,5],"\t",
		B[3,4],"\t",
		B[3,5],"\t",
		B[4,5],"\t",
		kappa(B),"\n",
		file='covarmwoslp.txt',
		append=TRUE
		)
		cat(
		dimnames(idat)[[2]][2],"\t",	# Gene name
		C[1,1],"\t",			# 15 relevant covariance values with the first 6 being the diagonals
		C[2,2],"\t",
		C[3,3],"\t",
		C[4,4],"\t",
		C[5,5],"\t",
		C[6,6],"\t",
		C[1,2],"\t",
		C[1,3],"\t",
		C[1,4],"\t",
		C[1,5],"\t",
		C[1,6],"\t",
		C[2,3],"\t",
		C[2,4],"\t",
		C[2,5],"\t",
		C[2,6],"\t",
		C[3,4],"\t",
		C[3,5],"\t",
		C[3,6],"\t",
		C[4,5],"\t",
		C[4,6],"\t",
		C[5,6],"\t",
		kappa(C),"\n",
		file='infom.txt',
		append=TRUE
		)
		cat(
		dimnames(idat)[[2]][2],"\t",	# Gene name
		D[1,1],"\t",			# 15 relevant covariance values with the first 6 being the diagonals
		D[2,2],"\t",
		D[3,3],"\t",
		D[4,4],"\t",
		D[5,5],"\t",
		D[1,2],"\t",
		D[1,3],"\t",
		D[1,4],"\t",
		D[1,5],"\t",
		D[2,3],"\t",
		D[2,4],"\t",
		D[2,5],"\t",
		D[3,4],"\t",
		D[3,5],"\t",
		D[4,5],"\t",
		kappa(D),"\n",
		file='infomwoslp.txt',
		append=TRUE
		)
		}
		
	plotf<-function(est,idat) {
		plotter<-function(t) {
			est[1]+est[2]*t+est[3]/sqrt(2*pi)*(integrate(function(z) cos(2*pi*t/(est[5]*exp(est[6]*z))+est[4])*exp(-.5*z^2),-10,10,subdivisions=10000)$val)
			}	
		p<-idat[1,1]:idat[length(idat[,1]),1]
		ep<-lapply(p,plotter)
		plot(idat,main=dimnames(idat)[[2]][2],xlab="Timepoints",ylab="Gene expression values",xlim=c(idat[1,1],idat[length(idat[,2]),1]),ylim=c(min(c(idat[,2],as.numeric(ep))),max(c(idat[,2]),as.numeric(ep))))
		points(p,ep,"l")	
		}

	pop<-expand.grid(			# Initialize population values
		c(-2,-1,0,1,2),
		c(-0.01,0,0.01),
		c(1,2,4),
		c(0,pi/6,pi/4,pi/3,pi/2,2*pi/3,3*pi/4,5*pi/6,pi,7*pi/6,5*pi/4,4*pi/3,3*pi/2,5*pi/3,7*pi/4,11*pi/6),
		c(60,70,80,90,100,110,120,130),
		c(0.001,0.01,0.1,0.2,0.3))
		
	if(graph) jpeg(filename="Rplot%05d.jpeg")	

		# Run for all genes
		for(n in 2:length(dat[1,])) {
			idat<-dat[,c(1,n)]				# Create individual gene data table
			idat<-idat[!is.na(idat[,2]),]	# Remove NA values from individual table
			
			evpop<-apply(pop,1,f,idat)		# Evaluate population values
			bpop<-pop[which.min(evpop),]	# Store population value with lowest cost as best value
			
			# Optimize best value
			est<-ucminf(bpop,f,gr=NULL,idat,control=list(maxeval=25000,xtol=1e-25))$par
			est1=est
			est1[4]=2*pi - est[4]
			# Calculate covariance matrix
			A<-covar(idat,est,TRUE,TRUE)
			B<-covar(idat,est,FALSE,TRUE)
			C<-covar(idat,est,TRUE,FALSE)
			D<-covar(idat,est,FALSE,FALSE)
			
			# Output results
			out(idat,est1,A,B,C,D,f,outloc)

			# Plot results
			if(graph) plotf(est,idat)
			}
		setwd(curdir)
		if(graph) graphics.off()
	}