//These default search parameters may be changed by assigning new values to the variables
BS_SPIKES_MIN = 3
BS_SI_MIN = 3


objref facTable
facTable = new Vector(171)
facTable.x(0) = 1
tmp = 1
for i=1,170 {
	tmp = tmp*i
	facTable.x(i) = tmp
}

func factorial() { //$1:n
	if ($1 > 170) {
//		print "factorial over limit ", $1
		return facTable.x(170)
	} else {
		return facTable.x($1)
	}
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
proc surpriseIndex() {local i,j,tmp,SI_MAX,lambda,explambda,poisscdf //$1:nSpikes, $2:timeInterval, $3:meanIsi, $&4:si, $&5:prb
////////////////
//function [si prb]=surpriseIndex(nSpikes, timeInterval, meanIsi)

SI_MAX = 100
prbMIN = exp(-SI_MAX)

//prb = 1 - poisscdf(nSpikes-1, timeInterval/ meanIsi);

lambda = $2/$3
if (lambda < 0) {
	print "surpriseIndex() lambda < 0"
	explambda = 0
} else {
	explambda = exp(-lambda)
}

//poisscdf = 0
//for i=0,$1-1 {
//	poisscdf = poisscdf + (explambda*lambda^i)/factorial(i)
//}

poisscdf = explambda
tmp = explambda
for i=1,$1-1 {
	tmp = tmp*(lambda/i)
	poisscdf = poisscdf + tmp
}


$&5 = 1 - poisscdf


if ($&5 <= prbMIN) {
	$&4 = SI_MAX
} else {
	$&4 = -log($&5)
}

}



////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
func findBurst() {local i,j,isiMean,bStart,bEnd,numSpikes,siMax,found,si,prb localobj x,ebIsi,b0bf //$o1:x, $o2:ebIsi, $3:isiMean, $4:bStart, $o5:b0bf
////////////////
//function [b_0, b_f]=findBurst(x, params, ebIsi, bStart)

x = $o1
ebIsi = $o2
isiMean = $3
bStart = $4
b0bf = $o5

b0bf.x(0) = bStart
b0bf.x(1) = bStart

bEnd = ebIsi.indwhere(">", bStart)

if (bEnd < 0) {
	bEnd = x.size()-1
} else {
	bEnd = ebIsi.x(bEnd)
}

//reject this burst segment if not enough spikes
numSpikes = bEnd - bStart + 1
if (numSpikes < BS_SPIKES_MIN) {
    	b0bf.x(1) = b0bf.x(0)
    	return -1
}

//find least probable (highest surprise) string of spikes in burst segment

siMax = 0
numSpikes = 0


for i=bStart,bEnd-BS_SPIKES_MIN+1 {
	for j=i+BS_SPIKES_MIN-1,bEnd {
		numIsi = j - i
		dur = x.x(j) - x.x(i)
		surpriseIndex(numIsi, dur, isiMean, &si, &prb)

		if (si > siMax) {
			found = 1
		} else {
			if ((si == siMax) && ((j-i+1) > numSpikes)) {
				found = 1
			} else {
				found = 0
			}
		}

		if (found > 0) {
			b0bf.x(0) = i
			b0bf.x(1) = j
			siMax = si
			numSpikes = j - i + 1
		}		
	}
}


if (siMax < BS_SI_MIN) {
    //No burst above surprise minimum in current burst segment
	b0bf.x(1) = b0bf.x(0)
	return -1
}


return 1

}


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
proc burstSearch() {local isiMean,flag,tmp,n localobj x,xIsi,res,bIsi,ebIsi,b0bf	//$o1:x, $o2:xIsi, $o3:result
////////////////
//function res=burstSearch(x, xIsi, isiMean)

x = $o1
xIsi = $o2
res = $o3
bIsi = new Vector()
ebIsi = new Vector()
b0bf = new Vector(2,0)

res.resize(1,1)
res.x[0][0] = 0


if (xIsi.size > 0) {
	isiMean = xIsi.mean()
} else {
	print "burstSearch xIsi.size() == 0"
	return
}


bIsi.indvwhere(xIsi, "()", 0, isiMean)
ebIsi.indvwhere(xIsi, "[)", isiMean, 1e9)

if (bIsi.size() > 0) {
	nextBurst = bIsi.x(0)
	flag = 1
} else {
	flag = 0
}

//keep finding bursts until starting points are exhausted 

while (flag > 0) {
    
	tmp = findBurst(x, ebIsi, isiMean, nextBurst, b0bf)

	if (tmp > 0) {
		n = res.nrow
		if (res.x[0][0] == 0) {
			res.resize(1,2)
		} else {
			n = n + 1
			res.resize(n, 2)
		}
		res.x[n-1][0] = b0bf.x(0)
		res.x[n-1][1] = b0bf.x(1)
	}

    	//next potential burst start

	tmp = ebIsi.indwhere(">=", b0bf.x(1))
	if (tmp >= 0) {
		tmp = ebIsi.x(tmp)
		tmp = bIsi.indwhere(">", tmp)
		if (tmp >= 0) {
			nextBurst = bIsi.x(tmp)
		} else {
			nextBurst = -1
		}
	} else {
		nextBurst = -1
	}

	flag = (nextBurst > 0)

}
   
return


}