###This code is very confusing.  I tried to clean it up.  The source file is much clearer
###The gist of our algorithm is that we fit a poisson regression with a quadratically increaing 
###Mean from the change-point.  Final is our algorithm's estimate.    


###REAL DATA###
source("functions.r")

bw=4
num<-1
tau<-NA
toD<-c()
Cnext<-1
cum.est<-0
xxx=c()

thetaIM<-m[num,] ###This is the spike data
mann<-thetaIM
mids<-c(-500:199)  ###This is the index of the data (There are 700 data points)


nn=1
for(ii in 1:length(thetaIM))
{
	if(thetaIM[ii]>0)
	{
		for(jj in 1:thetaIM[ii])
		{
			xxx[nn]=ii
			nn=nn+1
		}
	}

}


xxxx=c()

place=1

for(ii in 1:length(xxx))
{
	if(xxx[ii]<=560)
	if(xxx[ii]>=500)
	{
		xxxx[place]=xxx[ii]
		place=place+1
	}
}





###Histogram of the data

pdf(file="RealData1.pdf")
par(mar=(c(5,4.7,4,2)+.1))
par(font=1)
par(cex.axis=1.25)
hist(xxxx,ylim = c(0,6), axes=FALSE, cex.lab=1.5, cex.main=1.75, xlim=c(500,560), main=NULL,xlab="Time (s)",ylab="Spikes per s",col="white",border="black",breaks=100)
axis(1,pos=0)
axis(2,pos=500)
dev.off()


ppp<-man(mann)
cum.est<-as.numeric(ppp[2]-500)

for(i in 1:700)
{
	thetaIM
	j<-as.numeric(thetaIM[i])	
	while(j>0)
	{
		toD[Cnext]<-i
		j<-j-1
		Cnext<-Cnext+1
	}
}


##########We improved the fit of the algorithm by truncating the data after peak spike
##########We detected this peak by kernel density estimation

ddd<-density(toD,bw=bw)  


y<-as.matrix(ddd$y)
x<-as.matrix(ddd$x)
z<-matrix(nrow=dim(y)[1],ncol=2)

for(i in 1:dim(y)[1])
{
	z[i,1]<-y[i,1]
	z[i,2]<-x[i,1]
}

for(i in 1:dim(y)[1])
{
	if(z[i,1]==max(z[,1]))
		max<-ceiling(z[i,2])

}

countsss<-find.counts2(max,thetaIM)

midsNEW<-mids[c(1:length(countsss))]


######################################Here is the important part

final<-as.numeric(optimise(over,c(-100,length(midsNEW)-500),mids=midsNEW,counts=countsss,maximum=TRUE))[1]  ###Our Algorithm


##################################################


bins=length(midsNEW)
col<-c(nrow<-(bins+1))
pch<-c(nrow<-length(col))

for(i in 1:length(countsss))
{
	col[i]<-4
	pch[i]<-1
	if(mids[i]<final)
		{
			col[i]<-6
			pch[i]<-1
		}
}

xpoints<-midsNEW
ypoints<-countsss

xpoints<-c(midsNEW,final)
ypoints<-c(countsss,-2)

col[bins+1]<-3
pch[bins+1]<-2

xpoints<-c(xpoints,tau)
ypoints<-c(ypoints,-2)

col[bins+2]<-1
pch[bins+2]<-2

xpoints<-c(xpoints,cum.est)
ypoints<-c(ypoints,-2)





col[bins+3]<-5
pch[bins+3]<-2

xpoints<-xpoints[c(450:(length(countsss)+3))]
ypoints<-ypoints[c(450:(length(countsss)+3))]
col<-col[c(450:(length(countsss)+3))]
pch<-pch[c(450:(length(countsss)+3))]

pdf(file="RealResults1.pdf")


par(font=1)
par(cex.axis=1.25)

plot(xpoints,ypoints, xlim = c(-40,40), ylim = c(-1,8), col=col,xlab="Time (ms)",ylab="Spikes per ms",main=NULL,pch=pch)



###Fit for the graph

z<-regress.poiss(matr(midsNEW,final),countsss)
lines(c(-50,final),c(z$coeff[1],z$coeff[1]))
curve(z$coef[1]+z$coef[3]*(x-final)^2,final,-500+bins,add=TRUE)
dev.off()

final=final+500

xx=c()

place=1

for(ii in 1:length(xxx))
{
	if(xxx[ii]<=536)
	if(xxx[ii]>=500)
	{
		xx[place]=xxx[ii]
		place=place+1
	}
}



###Graphing the results
pdf(file="RealResults1.pdf")
par(mar=(c(5,4.7,4,2)+.1))
par(cex.axis=1.25)
par(cex.lab=1.5)
hist(xx,breaks=35, axes=FALSE, main=NULL,xlab="Time (s)",ylab="Spikes per s",col="white",xlim=c(500,536),border="black")
axis(1,pos=0,at=c(500,505,510,515,520,525,530,535))
axis(2,pos=500)
lines(c(500,final),c(z$coeff[1],z$coeff[1]))
curve(z$coef[1]+z$coef[3]*(x-final)^2,final,536,add=TRUE)
points(final,-.1,pch=24,bg="black",col="black")
dev.off()