// Simulate bar presentation, Papoutsi, 2017
//Incorporation of tag picking from Petousakis

//----------------Declare objects
//Spike trains
objref Stim_basal[total_syn_basal],Stim_apical[total_syn_apic], vc_stim_apical[total_syn_apic], vc_stim_basal[total_syn_basal]
//Stimulus synapses
objref stim_apical_ampa[total_syn_apic], stim_apical_nmda[total_syn_apic]
objref stim_basal_ampa[total_syn_basal], stim_basal_nmda[total_syn_basal]
//NetCons
objref nc_stim_apical_ampa[total_syn_apic], nc_stim_apical_nmda[total_syn_apic]
objref nc_stim_basal_ampa[total_syn_basal], nc_stim_basal_nmda[total_syn_basal]
// Random generators Random location of synapses
objref pstim, cut_off,s_time, basal_or
//Vectors 
objref X, sample, probs, mat_bas_st, mat_ap_st, tags_basal_syns[5], tag_dend_vec,tags_apical_syns, target

//Synapses accorind to length
mat_bas_st=new Vector()
mat_ap_st=new Vector()
tot_ap=0
tot_bas=0
for i=0, apical.count()-1 {tot_ap=tot_ap+apical.o(i).sec.L}
for i=0, apical.count()-1 {mat_ap_st.append(round((apical.o(i).sec.L/tot_ap)*total_syn_apic))}
for i=0, basal.count()-1 {tot_bas=tot_bas+basal.o(i).sec.L}
//Re-arrange to be primary dendrite-specific
for pd=0,4 {
	for i=0, list_primary[pd].count()-1 {
		mat_bas_st.append(round((list_primary[pd].o(i).sec.L/tot_bas)*total_syn_basal))
	}	
}
//Make activation function of each synapse 
func activation_vec(){
	pa=6
	a_Hz=(1/(pa*sqrt(2*PI)))*(exp((-0.5*(($1)^2))/(pa^2))+exp((-0.5*(($1-180)^2))/(pa^2))+exp((-0.5*(($1-360)^2))/(pa^2))+exp((-0.5*(($1+360)^2))/(pa^2))+exp((-0.5*(($1+180)^2))/(pa^2)))*15
	return a_Hz
}
//-----------------------------------------------------Set variables/functions for synapse tags----------------------------------------------------------------//
stmax=72  	//Maximum number of different synaptic orientation preferences (tags). 
func gauss5pdf(){  //This function describes the probability distribution function (pdf) for our gaussian distribution. We are using a five-term equation, wrapped around the unit circle.
		frac=(1/($1*sqrt(2*PI)))
		e1=(exp(-0.5*(($2-$3-360)^2)/($1^2)))  //2k*pi = -360, k=-1
		e2=(exp(-0.5*(($2-$3-180)^2)/($1^2)))  //2k*pi = -180, k=-0.5
		e3=(exp(-0.5*(($2-$3)^2)/($1^2)))      //2k*pi = 0, k=0
		e4=(exp(-0.5*(($2-$3+180)^2)/($1^2)))  //2k*pi = 180, k=0.5
		e5=(exp(-0.5*(($2-$3+360)^2)/($1^2)))  //2k*pi = 360, k=1
		prob=frac*(e1+e2+e3+e4+e5)  //This is a function of the synaptic tag ($2), and gives a probability value for every different tag.
		return prob
}
//Make reference stimulus vector
X=new Vector()
for i=0,stmax-1 {X.append(i*5)}

//-----------------------------------Stimulus procedure
proc stimulus() {
cut_off=new Random(num_seed2)	
cut_off.uniform(0,1)
s_time=new Random(num_seed)	
//Range of orientation preferences in basal dendrites 
basal_or=new Random(num_seed2-500000)	
tag_dend_vec=new Vector()
step_b=(disp*2)/4
if (disp==0) {
	target=new Vector(5,tag_basal)
} else {
	target=new Vector(5)
	target.indgen((tag_basal/10)-disp,(tag_basal/10)+disp, step_b)
}
for pd=0,4 { 
	basal_or.discunif(0, target.size()-1)	
	id_d=basal_or.repick()
	tag_basal_dend=target.x(id_d)*10
	target.remove(id_d)
	if (tag_basal_dend<=-10) {tag_basal_dend=180+tag_basal_dend}
	tag_dend_vec.append(tag_basal_dend)
}

syn=0
c_d=-1
//For list of primary dendrites and their child-all will have the same mean orientation
for pd=0,4 {

	//Traget orientation of this dendrite (and it's child)
	tag_basal_dend=tag_dend_vec.x(pd)
	// Fill in the probs and sums variables-linear arrays containing probability data.
	objref probs
	probs=new Vector()
	sums=0
	for k=0, stmax-1{  							//Iterates through all potential differences.
		st=k*5	
		probs.append(gauss5pdf(basal_width,st,tag_basal_dend))  	//Calls gauss5pdf() to calculate all individual probabilities, then stores them.
		sums=sums+probs.x(k) 						//Sums all individual probabilities and stores the sums.
	}	
	probs.div(sums)								//Normalize probabilities in 0-1 range
	tags_basal_syns[pd]=new Vector()							
	//For the basal dendrites of this branch (primary+child)
	for bd=0, list_primary[pd].count()-1 {
		matchcount=0
		c_d=c_d+1
		while (matchcount<= mat_bas_st.x(c_d)-1) { 	//For all the synapses in this dendrite
			matched=0  				//While "matched" is 0, the algorithm keeps trying to find a match for a given synapse.
			k=0  					//Variable to be used when iterating through all synaptic tags for a mc.
			mc=cut_off.repick() 			//Pick new random comparison value in the interval 0-1 (inclusive).
			while (matched==0){ 
				if(k==0){  
					rsum=probs.x(k)  //rsum is the sum of probabilities up to and including the current one.
				}
				if(k!=0){  //For all subsequent iterations.
					rsum=rsum+probs.x(k)  //As above, rsum is the sum of probabilities up to and including the current one.
				}
				if (mc<=rsum){  	//The value of rsum is the end point for the normalized probability of the current synaptic tag - beyond it lies the probability for the next one. 
					matched=1  	//As a match has been found, this breaks the inner loop, and returns to the outer loop.
					tag_syn=k*5					
					tags_basal_syns[pd].append(tag_syn)
					//Update this dendrite's and dendritic trees synapse counts
					matchcount=matchcount+1
					syn=syn+1
					//Create stimulation train of the synapse
					sample=new Vector()
					sample=X.c
					sample.sub(abs(tag_syn-stimulus_presentation))
					sample.abs()
					id=sample.min_ind()
					act_Hz=activation_vec(X.x(id))
					p_Hz=act_Hz*s_Hz
					s_time.poisson(p_Hz/1000)
					Stim_basal[syn-1]=new Vector()
					for tt=0, stimulus_duration-1 {
						if(s_time.repick()) {
							Stim_basal[syn-1].append(tt)
						}
					}

				}
				if(mc>rsum){  	//If the comparison value does not fall within the (nsum1, nsum2) range, we have a negative match.
					k=k+1  //Proceed to the next synaptic tag, as the one that was previously selected was not a match.
				}
			}
		}
	}
}	
matchcount=0
tags_apical_syns=new Vector()
objref probs
probs=new Vector()
sums=0
for k=0, stmax-1{  						//As before						
	st=k*5 	
	probs.append(gauss5pdf(apical_width,st,tag_apical))  	
	sums=sums+probs.x(k) 						
}	
probs.div(sums)	
while (matchcount<=total_syn_apic-1) { 		//For all the synapses in this tree
	matched=0  				
	k=0  					
	mc=cut_off.repick() 		
	while(matched==0){ 
		if(k==0){  
			rsum=probs.x(k)  
		}
		if(k!=0){ 
			rsum=rsum+probs.x(k)  
		}
		if(mc<=rsum){ 
			matched=1  		
			tag_syn=k*5					
			tags_apical_syns.append(tag_syn)
			//Update synapses
			matchcount=matchcount+1
			//Create stimulation of the synapse
			sample=new Vector()
			sample=X.c
			sample.sub(abs(tag_syn-stimulus_presentation))
			sample.abs()
			id=sample.min_ind()
			act_Hz=activation_vec(X.x(id))
			p_Hz=act_Hz*s_Hz
			s_time.poisson(p_Hz/1000)
			Stim_apical[matchcount-1]=new Vector()
			for tt=0, stimulus_duration-1 {
				if(s_time.repick()) {
					Stim_apical[matchcount-1].append(tt)
				}
			}
		}
		if(mc>rsum){  	
			k=k+1 
		}
	}
}	
pstim=new Random(num_seed2+400000)
pstim.uniform(0,1)	
//Connect basal synapses
syn=0
c_d=-1
for pd=0,4 {
	for num=0, list_primary[pd].count()-1 {
		c_d=c_d+1
		for i=0, mat_bas_st.x(c_d)-1 {
			vc_stim_basal[syn] = new VecStim(0.5)
			vc_stim_basal[syn].delay = 500
			vc_stim_basal[syn].play(Stim_basal[syn])
			PID=pstim.repick()
			list_primary[pd].o(num).sec stim_basal_ampa[syn]=new GLU(PID)
			list_primary[pd].o(num).sec stim_basal_nmda[syn]=new nmda(PID)

			if (!cut_basal) {
				nc_stim_basal_ampa[syn] = new NetCon(vc_stim_basal[syn], stim_basal_ampa[syn], -20, 0, ampa_g)
				nc_stim_basal_nmda[syn] = new NetCon(vc_stim_basal[syn], stim_basal_nmda[syn], -20, 0, nmda_g)
			} else if (pd!=dend_id1 && pd!=dend_id2) {
				nc_stim_basal_ampa[syn] = new NetCon(vc_stim_basal[syn], stim_basal_ampa[syn], -20, 0, ampa_g)
				nc_stim_basal_nmda[syn] = new NetCon(vc_stim_basal[syn], stim_basal_nmda[syn], -20, 0, nmda_g)
			}
			syn=syn+1
		}
	}
}
syn=-1
for num=0, apical.count()-1 {
	for i=0, mat_ap_st.x(num)-1 {
		syn=syn+1

		vc_stim_apical[syn] = new VecStim(0.5)
		vc_stim_apical[syn].delay = 500
		vc_stim_apical[syn].play(Stim_apical[syn])
	
		PID=pstim.repick()

		apical.o(num).sec stim_apical_ampa[syn]=new GLU(PID)
		apical.o(num).sec stim_apical_nmda[syn]=new nmda(PID)
		if (!ablated) {
			nc_stim_apical_ampa[syn] = new NetCon(vc_stim_apical[syn], stim_apical_ampa[syn], -20, 0, ampa_g)
			nc_stim_apical_nmda[syn] = new NetCon(vc_stim_apical[syn], stim_apical_nmda[syn], -20, 0, nmda_g)
		}
	
	}
}
}