#-------------------------------------------------------
#let ssimul1 the data set associated to the neuronal activity of the neuron 1
#let ssimul2 the data set associated to the neuronal activity of the neuron 2
#Both data sets are matrix with the same number of rows and the row i is associated to the trial i.
#The first column of ssimul1 (resp. ssimul2) contains the number of spikes of each trial.
#From column 2, we have the time of spikes. Thus, on the row i, from column 2 to column (1+ssimul1[i]) (resp. (1+ssimul2[i]))
#we have the times of spikes for the trial i of neuron 1 (resp. neuron 2) sorted in increasing order. The end of the row contains just 0. 

#the function creaXY1Y2e allows to compute:
#   - Xb: the average observed coincidence count for a delay delta on a window [a,b] associated to the datasets ssimul1 and ssimul2
#   - Y1b: (b-a)*the estimation of the firing rate of neuron 1
#   - Y2b: (b-a)*the estimation of the firing rate of neuron 2
 
#the output of this function is a vector of length 3 whose elements are Xb, Y1b and Y2b.
#------------------------------------------------------


creaXY1Y2e<-function(ssimul1,ssimul2,a,b,delta)
{
nessai=nrow(ssimul1)	
a1=0
a2=0
x=0
#creer N1j et N2j
	for (j in 1:nessai)
	{
	N1j=ssimul1[j,]	
	N2j=ssimul2[j,]
	n1=N1j[1]								
	n2=N2j[1]
	N1j=N1j[2:(n1+1)]							
	N2j=N2j[2:(n2+1)]
#creer N1j[a,b] et N2j[a,b]	
	linf=a
	lsup=b
	S1j=(linf<=N1j)*(N1j<=lsup)						
	S2j=(linf<=N2j)*(N2j<=lsup)						
	n1j=sum(S1j)						
	n2j=sum(S2j)
	a1=a1+n1j
	a2=a2+n2j
		if (n1!=0)
		{
			if (n2!=0)
			{
			SS1j=which(S1j==1)
			SS2j=which(S2j==1)
			nn1=length(SS1j)
			nn2=length(SS2j)
				if (nn1*nn2>=1)
				{	
				NN1j=N1j[SS1j]	
				NN2j=N2j[SS2j]
				NN1jm=NN1j-delta
				NN1jp=NN1j+delta
				nn2=length(NN2j)
					for (i in 1:nn2)
					{
					v=(NN1jm<=NN2j[i])*(NN2j[i]<=NN1jp)
					w=sum(v)
					x=x+w
					}
				}
			}
		}
	}
Xb=x/nessai
Y1b=a1/nessai
Y2b=a2/nessai
creaXY1Y2e<-c(Xb,Y1b,Y2b)
}


#-----------------------------------------
#The function testnous computes the p-value associated to the GAUE symmetric test and the numerator of this test.
#The arguments of this function are:
#  - C: a matrix with 3 rows. It contains the output of the function creaXY1Y2e 
#  - nessai: the number of trials of the datasets used to compute the average coincidence count with delay delta
#  - a and b: the lower limit and the upper limit of the time window on with C is computed
#  - delta: the delay

#The output of this function is a list with two elements, the first one is the p-value associated to the GAUE symmetric test
#and the second one is th enumerator of the statistic test.
#------------------------------------------

testnous<-function(C,nessai,a,b,delta)					
{
c=2*delta*(b-a)-delta^2
lambda1chapo=C[2,]/(b-a)
lambda2chapo=C[3,]/(b-a)
T1=sqrt(nessai)*(C[1,]-c*lambda1chapo*lambda2chapo)
T2=lambda1chapo*lambda2chapo*(c+(2/3*delta^3-delta^4/(b-a))*(lambda1chapo+lambda2chapo))
T=T1/sqrt(T2)
pv=pnorm(abs(T),mean = 0, sd = 1)					
alpha=2*(1-pv)									
testnous<-list(alpha,T1)							
}



#-----------------------------------------
#This function BHb allows to identify the small time interval detected by the MTGAUE procedure.
#The argument of this function are:
#  - A: a matrix with two rows. Each column is associated to a small time window [a,b] with a and the first row and b on the second.
#  - q: the maximal value for the false positive
#  - delta : the delay
#  - sssimul1 (resp. sssimul2): the matrix associated to the neuronal activity of neuron 1 (resp. neuron 2). The format of sssimul1 (resp. sssimul2) is the same as the
# one of ssimul 1 (resp. ssimul2) defined in creaXY1Y2e.

#The output ids a list with 4 elements:
#  - the number of detected small time window
#  - a matrix which identifies th edetected small windows. The first two rows od this matrix is A. The third row is a sequence of 0, 1 and -1, 0 meaning that
#	the small time window is not detected, -1 meaning that the small time window is detected and that the average coincidence count is smaller than the expected one
#	1 meaning that the small time window is detected and that the average coincidence count is bigger than the expected one (in case of independency) 
#  - the vector of p-values associated to each GAUE symmetric test
#  - the numerator of each GAUE symmetric statistic test
#------------------------------------------

BHb<-function(A,q,delta,sssimul1,sssimul2)				
						
{
nessai=nrow(sssimul1)
m=ncol(A)										 
pval=c()
cre=c()										
valeurtest=c()									
Sortie=matrix(data=0,ncol=m,nrow=3)
Sortie[1:2,]=A

#the GAUE symmetric test is performed on each small time window [a,b] defined in A

for (i in 1:m)									
  {
  a=A[1,i]									
  b=A[2,i]
  crea=creaXY1Y2e(sssimul1,sssimul2,a,b,delta)			
  crea=matrix(crea,ncol=1)
  test=testnous(crea,nessai,a,b,delta)
  alpha=test[[1]]		                        		
  valeurtest=c(valeurtest,test[[2]])
  pval=c(pval,alpha)							
  cre=c(cre,list(crea))	
}

lNa=which(pval=='NaN')							

llNa=length(lNa)

if (llNa==0)
  {
  pvalordd=sort(pval,index.return=TRUE)
  pvalord=pvalordd$x
  }
else
  {
  qqq=max(max(pval[-lNa]),q*1.001)
  pval[lNa]=rep(qqq,llNa)
  pvalordd=sort(pval,index.return=TRUE)
  pvalord=pvalordd$x
  }

j=m
										
while (pvalord[j]>j*q/m )						
  {
  j=j-1									
  if (j==0){break}
  }

m=j
if (m==0)
  {
  BHb=list(nbreplage=0,plage=Sortie,pvaleur=pval,test=valeurtest)
  }
else
  {
  pos=pvalordd$ix[1:m]
  Sortie[3,pos]=rep(1,length(pos))
  Sortie[3,pos]=Sortie[3,pos]-2*(valeurtest[pos]<0)
  BHb=list(m=m,plage=Sortie,pvaleur=pval,test=valeurtest) 	
  }														
}