// $Id: control.g,v 16.1 2000/03/02 14:43:39 ajay Exp $
// control.g genesis - 2.1
// functions for sequence storage and retrieval and performance measure.

include bin_pyr_spikes.g
include bin_inta_spikes.g
include bin_intb_spikes.g

//-------------------------------------------------------------------------------------------------

// add_sawtooth_depol creates table which generates a ramp waveform from ramp_curr_min (0 Amps)
// to ramp_curr_max (= curr_val Amps) followed by steady state (plateau) current injection.
// The table entries are values of current (in Amps) which are 
// injected into the soma compartment of every non_ILN (background) pyr cell.

function add_sawtooth_depol(curr_val, plateau_period, ramp_period)

	float	curr_val
	float	theta_period = {aff_input_interval} // == 100 ms

	float	xmin = 0
	float	xmax = theta_period / dt // one entry for every clock step in theta period
	float	xdivs = xmax

	float	baseline_period = 15e-3 // zero current for 18 ms since aff. input is at 15 ms.

	float	ramp_curr_min = 0 // min value in table
	float	ramp_curr_max = curr_val // max value of current to be injected (in Amps)
	float	ramp_period // = theta_period / 25  // current injected for this duration
	int	ramp_idx_min = baseline_period / dt
	int	ramp_idx_max = (baseline_period + ramp_period) / dt // max index used in for loop
	float	curr_slope = (ramp_curr_max - ramp_curr_min) / (ramp_idx_max - ramp_idx_min)

	// plateau current follows ramp current injection 
	int	plateau_idx_min // starting time of plateau current = end of ramp phase
	int	plateau_idx_max // ending time of plateau current
	float	plateau_period // 0.010, duration of plateau period in sec
	float	plateau_curr // plateau current amplitude = peak of ramp current amplitude

	int	i, nmes

	if (!{exists /saw_tab})
		echo creating saw_tab to inject depol current into non_ILN cells...
		create table /saw_tab
	end
	call /saw_tab TABCREATE {xdivs} {xmin} {xmax} // sets values to 0

	// provide output of table to all non_ILN pyr cells. 
	// First, delete any previously created outgoing messages

	nmes = {getmsg /saw_tab -out -count}
	for (i=0; i < nmes; i=i+1)
	  deletemsg /saw_tab 0 -outgoing
	end

	for (i= {num_ILN}; i < {NPYR}; i=i+1)
		addmsg /saw_tab /pyr_layer/pyr[{i}]/soma INJECT output
	end

	setfield /saw_tab step_mode 1 stepsize 1 // TAB_LOOP

	// have linearly increasing values from ramp_curr_min to ramp_curr_max for time period of ramp_period,
	// starting after the baseline period. 

	if (ramp_idx_max > xmax)
		ramp_idx_max = xmax // ensure current duration does not exceed theta period
		ramp_period = theta_period - baseline_period // reset value of ramp_period
	end

	for (i = 0; i < ramp_idx_min; i = i+1)
		setfield /saw_tab table->table[{i}] {0} // baseline period
	end

	for (i = ramp_idx_min; i < ramp_idx_max; i = i+1)
		setfield /saw_tab table->table[{i}] {(i-ramp_idx_min) * curr_slope}
	end

	plateau_curr = (i-1) * curr_slope / 5 // 1/5 of ramp_curr_max
	plateau_idx_min = ramp_idx_max

	if ((baseline_period +  ramp_period + plateau_period) > theta_period)
		plateau_period = theta_period - ramp_period - baseline_period
	end
	plateau_idx_max = (baseline_period + ramp_period + plateau_period) / dt
	for (i = plateau_idx_min; i < plateau_idx_max; i = i+1)
		setfield /saw_tab table->table[{i}] {plateau_curr} 
	end

	// add hyperpol current at end of theta cycle to rest rmp
	if (curr_val > 0)
		for (i = xmax-30; i <= xmax-15; i=i+1)
//			setfield /saw_tab table->table[{i}] -0.3e-9
		end
	end

	// call reset to make changes effective
//	reset
end

 add_sawtooth_depol {bg_depol_current} {bg_depol_plateau_period} {bg_depol_ramp_period}
// 0.18e-9 Amps, 0.010 s, 4 ms 

//-------------------------------------------------------------------------------------------------


// Reset initVm of non-ILN pyr cells from c_start to c_end-1 to initvm_val (depolarize them)

function reset_Em_pyr(c_start, c_end, initvm_val)
  int   c_start, c_end
  float initvm_val
  int   i
  str   compname
  str   channame
  str   elt_name
  float Itot = 0
  float newEm = 0
  float tempvar = 0

  echo resetting initVm of pyr[{c_start}] - [{c_end}] to {initvm_val} and adjusting Em accordingly...
 
  for (i= {c_start}; i < {c_end}; i=i+1) 
    pushe /pyr_layer/pyr[{i}]

    // First,  set initVm to initvm_val for all compartments so that
    //  upon reset Vm starts at initVm. (EREST_ACT defined in pyramidal.p)
    foreach compname ({el ##[OBJECT=symcompartment]})
	setfield {compname} initVm  {initvm_val}   // {EREST_ACT - 0.01}
	call {compname} RESET
    end

    foreach elt_name ({el /pyr_layer/pyr[{i}]/##})
	call {elt_name} RESET
    end  

    // Now calculate and set Em for every compartment under {CellRootDir}
    foreach compname ({el ##[OBJECT=symcompartment]})
	Itot = 0
	foreach channame ({el {compname}/##[OBJECT=tabchannel]})
	  Itot = Itot + {getfield {channame} Ik}
	end
	foreach channame ({el {compname}/##[OBJECT=vdep_channel]})
	  Itot = Itot + {getfield {channame} Ik}
	end
	tempvar = Itot * {getfield {compname} Rm}
	newEm = {getfield {compname} Vm} - tempvar
	setfield {compname} Em {newEm}
	call {compname} RESET
    end

    pope
  end // for
end

// reset_Em_pyr {num_ILN} {NPYR} {bg_steady_depol} // depolarizes non-ILN cells during learning phase

//-------------------------------------------------------------------------------------------------

// change_ILN_Ka is used to increase  Gbar of K_A in soma of ILN pyr cells (defined to be 50 S/m^2 
// = 1.66112e-7 S in pyramidal.p).

function change_ILN_Ka(Ka_gbar)
  float	Ka_gbar // S/m^2
  float	Ka_gbar_total // S. Got by multiplying Ka_gbar with surface area of soma
  float	pi = 3.1416
  float	len // length of soma of pyr cells
  float	dia // dia of soma of pyr cells
  int	i

  len = {getfield /pyr_layer/pyr[0]/soma len}
  dia = {getfield /pyr_layer/pyr[0]/soma dia}
  Ka_gbar_total = Ka_gbar * {pi} * {len} * {dia}

  for(i = 0; i < {num_ILN}; i = i + 1)
    setfield /pyr_layer/pyr[{i}]/soma/K_A Gbar {Ka_gbar_total}
  end
end

change_ILN_Ka {ILN_K_A}

//-------------------------------------------------------------------------------------------------

// Fix weights of AMPA hebbsynchans of int_a and int_b cells so that
// they do not participate in learning.

function fix_int_weights
	int i, isyn
	echo fixing weights of hebbsynchans on interneurons ...	
	for(i = 0; i < {NINTA}; i = i + 1)
	      setfield /int_a_layer/int_a[{i}]/apical/AMPA_channel \
	               change_weights 0
	end
	for(i = 0; i < {NINTB}; i = i + 1)
	      setfield /int_b_layer/int_b[{i}]/apical/AMPA_channel \
	               change_weights 0
	end
end

fix_int_weights

//-------------------------------------------------------------------------------------------------

// Create table c_ILN which stores for each non-ILN pyr cell (pyr[{num_ILN} - {NPYR-1}]), a '0'
// if that cell is not connected to an ILN, or '1' otherwise.
// 'Not connected' implies does not receive input from, nor sends output to an ILN.

function make_c_ILN
  int   c, isyn, numsyn
  float num_non_ILN = {NPYR - num_ILN} // number of non-ILN cells in network
  float xmin = 0
  float xmax = {NPYR - 1}
  float xdivs = {xmax}
  str   src // src element of synapse
  str   dest // dest. elt. of synapse
  str   cell_type
  int   id_open_bracket 
  int   id_close_bracket
  int   id_first_underscore
  str   cell_id_str
  int   cell_id // index in c_ILN table

  if (num_non_ILN > 0)

    // create table called c_ILN to store 1 or 0 for each non-ILN that is either connected (1)
    // or not connected (0) to an ILN.
    if ({exists /c_ILN}) 
	echo /c_ILN already exists. Returning from make_c_ILN without any changes
	return
    end

    echo creating table /c_ILN to indicate if non-ILN cell is connected to an ILN ...
    create table /c_ILN
    call /c_ILN TABCREATE {xdivs} {xmin} {xmax} // resets elements to 0	

    for(c = 0; c < {num_ILN}; c = c + 1) // loop over ILNs
	// first loop over all its synapses to identify all non-ILNs that it receives input from.
	// Set value in /c_nc_ILN to 1 for those non-ILNs.
	numsyn = {getfield /pyr_layer/pyr[{c}]/apical/AMPA_channel nsynapses}
	for (isyn = 0; isyn < numsyn; isyn = isyn+1)
	  src = {getsynsrc /pyr_layer/pyr[{c}]/apical/AMPA_channel {isyn}}

	  // get index of cell number from name string
	  id_open_bracket = {findchar {src} "["}
	  id_close_bracket = {findchar {src} "]"}
	  if ((id_open_bracket != -1) && (id_close_bracket != -1))
		cell_id_str = {substring {src} {id_open_bracket + 1} {id_close_bracket - 1}}
		cell_id = cell_id_str
	  else // pyr[0]
		cell_id = 0
	  end
	  if (cell_id >= num_ILN) // non-ILN cell
		setfield /c_ILN table->table[{cell_id}] 1
	  end
	end

	// Now loop over its spike messages (outgoing synapses) to identify all non-ILNs that
	// it sends output to and set corresponding values to 1 in c_ILN.
	numsyn = {getsyncount /pyr_layer/pyr[{c}]/soma/spike}
	for (isyn = 0; isyn < numsyn; isyn = isyn+1)
	  dest = {getsyndest /pyr_layer/pyr[{c}]/soma/spike {isyn}}

	  // Note: outgoing spike messages also go to NMDA channels on pyr cells as well as to
	  // AMPA channels on int_a and int_b cells. We only want messages going to pyr cells.
	  // Differentiating between AMPA and NMDA channels on pyr cells is not necessary, since
	  // if it goes to AMPA channel, it also goes to NMDA channel.
	  id_first_underscore = {findchar {dest} "_"}
	  cell_type = {substring {dest} 1 {id_first_underscore - 1}}
	  if ({cell_type} == "pyr")
		// get index of cell number from name string
		id_open_bracket = {findchar {dest} "["}
		id_close_bracket = {findchar {dest} "]"}
		if ((id_open_bracket != -1) && (id_close_bracket != -1))
		  cell_id_str = {substring {src} {id_open_bracket + 1} {id_close_bracket - 1}}
		  cell_id = cell_id_str
		else // pyr[0]
		  cell_id = 0
		end
		if (cell_id >= num_ILN) // non-ILN cell
			setfield /c_ILN table->table[{cell_id}] 1
		end
	  end
	end

    end // for

  end // if num_non_ILN > 0
end

// make_c_ILN


//--------------------------------------------------------------------------------------------


// Fix weights of AMPA hebbsynchans of non-ILN pyr cells that are also 
// "not connected to Input Layer Pyr Cells" (ncILNs) so that they do not participate in learning.
// 'Not connected' implies does not receive input from, nor sends output to an ILN.

function fix_ncILN_weights
  float num_non_ILN = {NPYR - num_ILN} // number of non-ILN cells in network
  int   cell_id // index in c_ILN table

  if (num_non_ILN > 0)

    if (!{exists /c_ILN})
	make_c_ILN
    end

    // fix weights of all non-ILNs that are not connected to ILNs
    echo fixing weights of hebbsynchans on non-ILN pyr cells not connected to ILNs ...
    for (cell_id = {num_ILN}; cell_id < {NPYR}; cell_id = cell_id+1)
	if ({getfield /c_ILN table->table[{cell_id}]} == 0) // not connected to ILN
	  setfield /pyr_layer/pyr[{cell_id}]/apical/AMPA_channel change_weights 0
	  echo not connected non-ILN: pyr[{cell_id}]
	end
    end
  end // if num_non_ILN > 0
end

// fix_ncILN_weights

//--------------------------------------------------------------------------------------------

// initialize non-changeable weights to some fixed value.
// non-changeable weights are weights of synapses that do not undergo
// hebbian learning. These are NMDA, GABAa and GABAb for pyr cells,
// and AMPA, GABAa and GABAb for int_a and int_b cells.

function init_nonchg_weights
  int c, imes, i, idx

  echo initializing non-changeable weights on pyr, int_a and int_b cells...

  // pyr cells (apical - NMDA, GABAb; soma - GABAa)

  for(c = 0; c < {NPYR}; c = c + 1)
    imes = {getfield /pyr_layer/pyr[{c}]/apical/NMDA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/apical/NMDA_channel \
		synapse[{i}].weight {pyr_nmda_wt}
    end

    imes = {getfield /pyr_layer/pyr[{c}]/apical/GABAb_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/apical/GABAb_channel \
		synapse[{i}].weight {pyr_gabab_wt}
    end
    
    imes = {getfield /pyr_layer/pyr[{c}]/soma/GABAa_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/soma/GABAa_channel \
		synapse[{i}].weight {pyr_gabaa_wt}
    end

    // reset weight of GABAa synapse from feedforward interneuron back to {pyr_ffi_wt}
    idx = {getsynindex /ffi_layer/int_a[0]/soma/spike /pyr_layer/pyr[{c}]/soma/GABAa_channel}
    setfield /pyr_layer/pyr[{c}]/soma/GABAa_channel synapse[{idx}].weight {pyr_ffi_wt}

  end


  // int_a cells (apical - AMPA, GABAb; soma - GABAa)

  for(c = 0; c < {NINTA}; c = c + 1)
    imes={getfield /int_a_layer/int_a[{c}]/apical/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_a_layer/int_a[{c}]/apical/AMPA_channel \
		synapse[{i}].weight {inta_ampa_wt}
    end    

    imes={getfield /int_a_layer/int_a[{c}]/apical/GABAb_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_a_layer/int_a[{c}]/apical/GABAb_channel \
		synapse[{i}].weight {inta_gabab_wt}
    end

    imes = {getfield /int_a_layer/int_a[{c}]/soma/GABAa_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_a_layer/int_a[{c}]/soma/GABAa_channel \
		synapse[{i}].weight {inta_gabaa_wt}
    end
  end

  // int_b cells (apical - AMPA, GABAb; soma - GABAa)

  for(c = 0; c < {NINTB}; c = c + 1)
    imes={getfield /int_b_layer/int_b[{c}]/apical/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_b_layer/int_b[{c}]/apical/AMPA_channel \
		synapse[{i}].weight {intb_ampa_wt}
    end    

    imes={getfield /int_b_layer/int_b[{c}]/apical/GABAb_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_b_layer/int_b[{c}]/apical/GABAb_channel \
		synapse[{i}].weight {intb_gabab_wt}
    end
    imes = {getfield /int_b_layer/int_b[{c}]/soma/GABAa_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_b_layer/int_b[{c}]/soma/GABAa_channel \
		synapse[{i}].weight {intb_gabaa_wt}
    end
  end

end

// init_nonchg_weights

//--------------------------------------------------------------------------------------------

// connect_ILN_in_pat connects all the ILNs within each pattern (if not already connected),
// and sets their weights to high values. This ensures that all cells of a pattern fire
// nearly synchronously, as soon as even one of them gets activated.

function connect_ILN_in_pat
  int	patnum, prenum, postnum
  int	preidx, postidx
  int	nsyn, index
  float	ILN_wt = 25.0 // weight of synapse between ILNs of same pattern

  echo connecting all ILNs within each pattern and setting their weights to {ILN_wt}...

  for (patnum = 0; patnum < num_aff_pat; patnum=patnum+1)
    for (preidx = 0; preidx < 5; preidx=preidx+1) // 5 cells per pattern
	for (postidx = 0; postidx < 5; postidx=postidx+1)
	  prenum = preidx + (patnum*5)
	  postnum = postidx + (patnum*5)
	  nsyn = {getsyncount /pyr_layer/pyr[{prenum}]/soma/spike \
			/pyr_layer/pyr[{postnum}]/apical/AMPA_channel}
	  if (nsyn > 0) // connection exists
	  	index = {getsynindex /pyr_layer/pyr[{prenum}]/soma/spike \
			/pyr_layer/pyr[{postnum}]/apical/AMPA_channel}
		setfield /pyr_layer/pyr[{postnum}]/apical/AMPA_channel \
			synapse[{index}].weight {ILN_wt}
	  else
		if (prenum != postnum) // do not form self connections
			addmsg /pyr_layer/pyr[{prenum}]/soma/spike \
				/pyr_layer/pyr[{postnum}]/apical/AMPA_channel SPIKE
			setfield /pyr_layer/pyr[{postnum}]/apical/AMPA_channel \
				synapse[0].weight {ILN_wt}
		end
	  end
	end
    end
  end
end

// connect_ILN_in_pat

//--------------------------------------------------------------------------------------------

// initialize GABAa weights onto pyr cells (pyr_gabaa_wt).

function reset_pyr_gaba_weights
  float	gabaa_weight
  float	gabab_weight
  float	old_wt, new_wt
  int 	c, imes, i

  echo multiplying GABAa and GABAbweights on pyr cells by {recall_gaba_wt_mult}...

  for(c = 0; c < {NPYR}; c = c + 1)
    // GABA_A weights
    imes = {getfield /pyr_layer/pyr[{c}]/soma/GABAa_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	old_wt = {getfield /pyr_layer/pyr[{c}]/soma/GABAa_channel synapse[{i}].weight}
	new_wt = old_wt * recall_gaba_wt_mult
	setfield /pyr_layer/pyr[{c}]/soma/GABAa_channel \
		synapse[{i}].weight {new_wt}
    end


    // GABA_B weights
    imes = {getfield /pyr_layer/pyr[{c}]/apical/GABAb_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	old_wt = {getfield /pyr_layer/pyr[{c}]/apical/GABAb_channel synapse[{i}].weight}
	new_wt = old_wt * recall_gaba_wt_mult
	setfield /pyr_layer/pyr[{c}]/apical/GABAb_channel \
		synapse[{i}].weight {new_wt}
    end
  end

end

//---------------------------------------------------------------------------------------------------------

// init_syn_delays initializes the delays of all synapses to fixed values (set in constants.g).
// The synaptic delay includes the axon conduction time.

function init_syn_delays
  int c, imes, i, idx

  echo initializing synaptic delays...


  // pyr cells (apical - AMPA, NMDA, GABAb; soma - GABAa)

  for(c = 0; c < {NPYR}; c = c + 1)

    imes = {getfield /pyr_layer/pyr[{c}]/apical/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/apical/AMPA_channel \
		synapse[{i}].delay {p_to_p_syndelay}
    end

    imes = {getfield /pyr_layer/pyr[{c}]/apical/NMDA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/apical/NMDA_channel \
		synapse[{i}].delay {p_to_p_syndelay}
    end

    imes = {getfield /pyr_layer/pyr[{c}]/apical/GABAb_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/apical/GABAb_channel \
		synapse[{i}].delay {i_to_p_syndelay}
    end
    
    imes = {getfield /pyr_layer/pyr[{c}]/soma/GABAa_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/soma/GABAa_channel \
		synapse[{i}].delay {i_to_p_syndelay}
    end

  end


  // int_a cells (apical - AMPA, GABAb; soma - GABAa)

  for(c = 0; c < {NINTA}; c = c + 1)
    imes={getfield /int_a_layer/int_a[{c}]/apical/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_a_layer/int_a[{c}]/apical/AMPA_channel \
		synapse[{i}].delay {p_to_i_syndelay}
    end    

    imes={getfield /int_a_layer/int_a[{c}]/apical/GABAb_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_a_layer/int_a[{c}]/apical/GABAb_channel \
		synapse[{i}].delay {i_to_i_syndelay}
    end

    imes = {getfield /int_a_layer/int_a[{c}]/soma/GABAa_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_a_layer/int_a[{c}]/soma/GABAa_channel \
		synapse[{i}].delay {i_to_i_syndelay}
    end
  end

  // int_b cells (apical - AMPA, GABAb; soma - GABAa)

  for(c = 0; c < {NINTB}; c = c + 1)
    imes={getfield /int_b_layer/int_b[{c}]/apical/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_b_layer/int_b[{c}]/apical/AMPA_channel \
		synapse[{i}].delay {p_to_i_syndelay}
    end    

    imes={getfield /int_b_layer/int_b[{c}]/apical/GABAb_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_b_layer/int_b[{c}]/apical/GABAb_channel \
		synapse[{i}].delay {i_to_i_syndelay}
    end
    imes = {getfield /int_b_layer/int_b[{c}]/soma/GABAa_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /int_b_layer/int_b[{c}]/soma/GABAa_channel \
		synapse[{i}].delay {i_to_i_syndelay}
    end
  end

end

init_syn_delays

//---------------------------------------------------------------------------------------------------------

// reset_p2p_syndelay(p2pd) changes the delay between p2p connections just prior to recall
// phase, so that overlapping (jitter) of recalled patterns may be avoided.

function reset_p2p_syndelay(p2pd)
  float p2pd // synaptic delay
  int	c, imes, i

  // pyr cells (apical - AMPA, NMDA, GABAb; soma - GABAa)

  for(c = 0; c < {NPYR}; c = c + 1)

    imes = {getfield /pyr_layer/pyr[{c}]/apical/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/apical/AMPA_channel \
		synapse[{i}].delay {p2pd}
    end

    imes = {getfield /pyr_layer/pyr[{c}]/apical/NMDA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/apical/NMDA_channel \
		synapse[{i}].delay {p2pd}
    end
  end
end

//---------------------------------------------------------------------------------------------------------

// initialize changeable weights (AMPA of pyr cells) [random values ??]

function init_hebb_weights
  int c, imes, i

  echo initializing Hebbsynchans on pyr cells ...
  for(c = 0; c < {NPYR}; c = c + 1)
    imes = {getfield /pyr_layer/pyr[{c}]/apical/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/apical/AMPA_channel \
		synapse[{i}].weight {pyr_ampa_wt}
    end

    imes = {getfield /pyr_layer/pyr[{c}]/apical_2/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	setfield /pyr_layer/pyr[{c}]/apical_2/AMPA_channel \
		synapse[{i}].weight {pyr_EC_wt}
    end

  end
end

// init_hebb_weights

//-------------------------------------------------------------------------------------------

// multiply p->p weights OF POTENTIATED SYNAPSES onto ILN cells ONLY by recall_exc_wt_mult.
// This has to be done because of high K_A on ILN cells, preventing much learning.
// Alternative is ti decrease post_tau of ILN AMPA hebbsynchans to 1ms (from 3 ms)
// This is done prior to recall phase. 

function reset_pyr_ampa_wt
  int	c, imes, i
  float	old_wt, new_wt

  echo multiplying pyr_ampa_wt of potentiated synapses on ILN cells ONLY by recall_exc_wt_mult ...
  for(c = 0; c < {num_ILN}; c = c + 1)
    imes = {getfield /pyr_layer/pyr[{c}]/apical/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	old_wt = {getfield /pyr_layer/pyr[{c}]/apical/AMPA_channel synapse[{i}].weight}
	if (old_wt > {pyr_ampa_wt}) // synapse has gotten potentiated
		new_wt = old_wt * recall_exc_wt_mult
//		new_wt = pyr2ILN_maxweight
		setfield /pyr_layer/pyr[{c}]/apical/AMPA_channel \
			synapse[{i}].weight {new_wt}
	end
    end
  end
end

//-------------------------------------------------------------------------------------------

// set EC_to_pyr weights of POTENTIATED EC SYNAPSES ONLY to EC2pyr_recall_maxwt_{ILN/bg}
// This is done prior to recall phase

function  reset_EC2pyr_wt
  int	c, imes, i
  float	old_wt, new_wt

  echo resetting weights of potentiated EC to ILN synapses to  {EC2pyr_recall_maxwt_ILN}

  // change weights of synapses from EC to ILN pyr cells
  for(c = 0; c < {num_ILN}; c = c + 1)
    imes = {getfield /pyr_layer/pyr[{c}]/apical_2/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	old_wt = {getfield /pyr_layer/pyr[{c}]/apical_2/AMPA_channel synapse[{i}].weight}
	if (old_wt > {pyr_EC_wt}) // synapse has gotten potentiated
		new_wt = EC2pyr_recall_maxwt_ILN
		setfield /pyr_layer/pyr[{c}]/apical_2/AMPA_channel \
			synapse[{i}].weight {new_wt}
	end
    end
  end

  // change weights of synapses from EC to bg pyr cells
  for(c = num_ILN; c < NPYR; c = c + 1)
    imes = {getfield /pyr_layer/pyr[{c}]/apical_2/AMPA_channel nsynapses} - 1
    for(i = 0; i <= imes; i = i + 1)
	old_wt = {getfield /pyr_layer/pyr[{c}]/apical_2/AMPA_channel synapse[{i}].weight}
	if (old_wt > {pyr_EC_wt}) // synapse has gotten potentiated
		new_wt = EC2pyr_recall_maxwt_bg
		setfield /pyr_layer/pyr[{c}]/apical_2/AMPA_channel \
			synapse[{i}].weight {new_wt}
	end
    end
  end

end

//-------------------------------------------------------------------------------------------

// set_EC2pyr_weights set weights of synapses from EC[ecnum] to ILN cells of seq [1 or both 1&2],
// to the value in EC2pyr_maxweight
// Note, synapse[0] is synapse from EC[0] and synapse[1] from EC[1]

function set_EC2pyr_weights

  int	i

  echo Setting weights of synapses from EC cells to ILN cells...

  if (num_seq==1)
    for (i=5; i<num_ILN; i=i+1) // skip input to recall cue pattern
	setfield /pyr_layer/pyr[{i}]/apical_2/AMPA_channel synapse[0].weight {EC2pyr_maxweight}
    end 
  elif (num_seq==2)
    for (i=5; i<25; i=i+1)
	setfield /pyr_layer/pyr[{i}]/apical_2/AMPA_channel synapse[0].weight {EC2pyr_maxweight}
    end 
    for (i=25; i<45; i=i+1) // 2nd seq has cells pyr[0-4, 25-44]
	setfield /pyr_layer/pyr[{i}]/apical_2/AMPA_channel synapse[1].weight {EC2pyr_maxweight}
    end 
  end

end

//-------------------------------------------------------------------------------------------

// functions to allow or disallow weight changes of pyr AMPA_channels

function fix_hebb_weights
	int i, isyn
	echo Hebbsynchans fixed ...	
	for(i = 0; i < {NPYR}; i = i + 1)
	      // pyr_to_pyr connections
	      setfield /pyr_layer/pyr[{i}]/apical/AMPA_channel change_weights 0

	      // EC_to_pyr connections
	      setfield /pyr_layer/pyr[{i}]/apical_2/AMPA_channel change_weights 0
	end
end

function free_hebb_weights
	int i, isyn
	echo Hebbsynchans free to change ...
	for(i = 0; i < {NPYR}; i = i + 1)
	      // pyr_to_pyr connections
	      setfield /pyr_layer/pyr[{i}]/apical/AMPA_channel change_weights 1

	      // EC_to_pyr connections
	      setfield /pyr_layer/pyr[{i}]/apical_2/AMPA_channel change_weights 1
	end
end


// function below fixes (val == 0) or frees (val == 1) AMPA weights of only ILN cells
// and those non-ILN pyr cells that are connected to ILNs.

function change_ILN_hebb_weights(val)
  int   val
  int   i
  float num_non_ILN = {NPYR - num_ILN} // number of non-ILN cells in network
  int   cell_id // index in c_ILN table

  if ({val} == 0)
    echo Hebbsynchans of ILNs and 'connected non-ILNs' fixed ...	
  elif ({val} == 1)
    echo Hebbsynchans of ILNs and 'connected non-ILNs' free to change ...	
  else
    echo argument of function has to be 0 (fix) or 1 (free). Aborting function...
    return	
  end

  // fix or free  weights of ILNs
  for(i = 0; i < {num_ILN}; i = i + 1)
	setfield /pyr_layer/pyr[{i}]/apical/AMPA_channel change_weights {val}
  end

  // fix or free weights of 'connected non-ILNs
  if (num_non_ILN > 0)
    if (!{exists /c_ILN})
	make_c_ILN
    end

    for (cell_id = {num_ILN}; cell_id < {NPYR}; cell_id = cell_id+1)
	if ({getfield /c_ILN table->table[{cell_id}]} == 1) // connected to an ILN
		setfield /pyr_layer/pyr[{cell_id}]/apical/AMPA_channel change_weights {val}
	end
    end
  end // if num_non_ILN > 0

end


//----------------------------------------------------------------------------------------------------------

// ff_off turns off activation of feedforward inhibition. This is done just prior to start of recall phase

function ff_off
  int c, idx, i

  intb_ampa_wt = 0

  echo Turning off feedforward inhibition by int_a cell

  for(c = 0; c < {NPYR}; c = c + 1)
    idx = {getsynindex /ffi_layer/int_a[0]/soma/spike /pyr_layer/pyr[{c}]/soma/GABAa_channel}
    setfield /pyr_layer/pyr[{c}]/soma/GABAa_channel synapse[{idx}].weight 0
  end    

end

//-----------------------------------------------------------------------------------------


// Performance of network on sequence recogniion task is determined by
// calculating its 'Recall Accuracy' or RCL_ACC. To get this measure, neuronal
// activity is binned in windows of duration (.005 * num_bins) ms, (5 ms is bin_width
// used to store spike data). num_bins is varied from 1 to MAX_NUM_BINS.
// The start of a window lies between 1 bin_width after the start of the previous pattern's
// window (called new_pat_win_start) and (new_pat_win_start + 100ms) and is equal to the
// earliest spike time of a cell in the current pattern within this interval. If no appropriate
// ILN fires within this interval (one pattern skipped), the sequence recall is considered to be
// over, even though the following patterns may fire in sequence.
   
// RCL_ACC is calculated by looking at activities of input layer neurons (ILNs), 
// which are the pyr cells receiving afferent input, pyr[0 - (5*num_pattern - 1)], and which
// code the patterns in the sequnce. If even one ILN fires in a window, a prob_pattern
// value is determined for that pattern (window). This is given by:
// prob_pattern = (#ILN_supposed_to_fire that fire / 5) - (other_ILNs that fire / (num_ILN-5))
// (For e.g., pyr[0-4] are supposed to fire in the first pattern, and pyr[5- (num_ILN-1)]
// are the 'other_ILNs' that are not supposed to fire for that pattern.).
// The sum of prob_pattern over all the patterns (except first recall cue), divided by the
// sequence length (num_pattern-1), gives the RCL_ACC value.
// The recall_start_time argument should indicate time after recall_cue


function gpm3(recall_start_time, num_pattern)
  float recall_start_time // RCL_ACC measured after learning phase completed
  int   num_pattern    // number of patterns in sequence
  int   pattern_idx  // index that loops over patterns in the sequence
  int   cell_idx  // index that loops over appropriate cells in a pattern
  int   pat_cell_start  // start value of cell_idx in a pattern
  int   pat_cell_end  // end value of cell_idx in a pattern
  float num_fired_in_pattern  // num. of ILNs firing that are part of current pattern
  float num_fired_out_pattern // num. of ILNs firing that are not in current pattern
  int   i
  int   stab_idx = 0 // index into stab_i, got from stab_indices
  int   num_ILN_spikes_over = 0 // number of ILNs that do not spike anymore
  int   no_pattern_skipped = 1 // bool var to indicate whether a pattern is skipped during recall
  float cell_spike_time // value got from pyr_stab_i
  float prob_sequence = 0  // sum of prob_pattern over all patterns in sequence
  float prob_pattern = 0 // (# appropriate cells firing) / 5
  float max_prob_pattern // max over various win_lengths
  float xdivs = {num_ILN-1} // used for table /curr_stab_idx.
  float xmin = 0
  float xmax = {xdivs}
  float new_pat_win_start // start time of previous window
  float curr_win_start // start time of current window
  float	win_start_interval = 0.002 // 2 ms interval between start of curr. pat and start of next pat.
  float curr_win_end   // end time of current window
  float max_spike_time, min_spike_time
  int   num_pat_cell_fired
  int   nb_idx, num_bins
  float bin_int = 0.0003 // 0.3ms increments in interval of window
  int   MAX_NUM_BINS = 50 // maximum window duration is bin_int * MAX_NUM_BINS
  float MAX_PAT_SEPARATION = 0.1 // next pattern should occur within 100ms of new_pat_win_start 
  echo Measuring performance on sequence recall from time {recall_start_time} ....

  // table '/curr_stab_idx' stores for each ILN, the current position in the stab table
  // If that ILN has no more spikes (-1 reached in stab table), then /curr_stab_idx also
  // stores -1 for that cell.
  if (!{exists /curr_stab_idx}) 
	create table /curr_stab_idx
  end
  call /curr_stab_idx TABCREATE {xdivs} {xmin} {xmax} // resets elements to 0
  num_ILN_spikes_over = 0 // set to 1 when no more ILNs left to fire
  no_pattern_skipped = 1 // set to 0 if a pattern is not recalled in sequence. Terminates scoring.
  prob_sequence = 0

  new_pat_win_start = recall_start_time

  for (pattern_idx=1; pattern_idx < num_pattern; pattern_idx=pattern_idx+1) // start after cue pattern
    if ((num_ILN_spikes_over < num_ILN) && (no_pattern_skipped == 1)) 
	pat_cell_start = pattern_idx * 5 // 5 cells per pattern
	pat_cell_end = pat_cell_start + 4 

	// for cells in pattern get earliest spike time after new_pat_win_start.
	// Also get min and max of these spike times
	num_pat_cell_fired = 0
	min_spike_time = 1000 // initialized to a large value
	max_spike_time = -1  // initialized to below possible min

	for (cell_idx = pat_cell_start; cell_idx <= pat_cell_end; cell_idx=cell_idx+1)
	  stab_idx = {getfield /curr_stab_idx table->table[{cell_idx}]} 
	  if ( stab_idx != -1) // if cell still to fire
	      cell_spike_time = {getfield /stables/pyr_stab_{cell_idx} table->table[{stab_idx}]}
	      while ((stab_idx < (nbins-1)) && (cell_spike_time != -1) &&  \
			(cell_spike_time < new_pat_win_start)) // 
		stab_idx = stab_idx + 1
		cell_spike_time = {getfield /stables/pyr_stab_{cell_idx} table->table[{stab_idx}]}
	      end

	      // If -1 or end of table, then no more spikes >= new_pat_win_start
	      if ( (cell_spike_time == -1) || ( stab_idx == (nbins-1) ) )
		setfield /curr_stab_idx table->table[{cell_idx}] -1
		num_ILN_spikes_over = num_ILN_spikes_over + 1
//		echo no more spikes for pyr {cell_idx} after {new_pat_win_start}
	      else
//		echo cell {cell_idx}, stime = {cell_spike_time}
		setfield /curr_stab_idx table->table[{cell_idx}] {stab_idx} // start from same point
		num_pat_cell_fired = num_pat_cell_fired + 1
		if (cell_spike_time < min_spike_time)
			min_spike_time = cell_spike_time
		end
		if (cell_spike_time > max_spike_time)
			max_spike_time = cell_spike_time
		end
	      end
	  end // if stab_idx != -1
	end  // for cell_idx
//	echo pat {pattern_idx} min_spike_time {min_spike_time}
//	echo pat {pattern_idx} max_spike_time {max_spike_time}

	if ((num_pat_cell_fired == 0) || (min_spike_time > (new_pat_win_start + MAX_PAT_SEPARATION)))
	  // ***************************************************************
	  // pattern has been skipped. Terminate prob_pattern computations
	  // ***************************************************************
	  no_pattern_skipped = 0
	else // go on with computation of prob_pattern for this and following patterns
	  curr_win_start = min_spike_time - bin_int // include cells firing one bin_int prior
							// to min_spike_time
	  new_pat_win_start = curr_win_start + win_start_interval // 2 ms after curr_win_start
	  num_bins = {round {(max_spike_time - curr_win_start) / bin_int}} + 1
	  if (num_bins > MAX_NUM_BINS)
		num_bins = MAX_NUM_BINS
	  elif (num_bins == 0)
		num_bins = 1
	  end
//	  echo pat {pattern_idx} num_bins {num_bins}

	  // find prob_pattern for each win_length (nb_idx from 1 to  num_bins)
	  // and select max prob_pattern
	  max_prob_pattern = -1
	  for (nb_idx = 1; nb_idx <= num_bins; nb_idx = nb_idx+1)
	    curr_win_end = curr_win_start + (nb_idx * bin_int) 
//	    echo pat {pattern_idx} curr_win_start {curr_win_start}
//	    echo pat {pattern_idx} curr_win_end {curr_win_end}
	    num_fired_in_pattern = 0
	    num_fired_out_pattern = 0 

	    for (cell_idx = 0; cell_idx < num_ILN; cell_idx=cell_idx+1)
	      // get first spike time >= curr_win_start from /stables/pyr_stab_{cell_idx}
	      // Do this only if this cell has not already stopped firing
	      stab_idx = {getfield /curr_stab_idx table->table[{cell_idx}]} 
	      if ( stab_idx != -1)
	        cell_spike_time = {getfield /stables/pyr_stab_{cell_idx} table->table[{stab_idx}]}
	        while ((stab_idx < (nbins-1)) && (cell_spike_time != -1) &&  \
			(cell_spike_time < curr_win_start)) // -1  => end of spike data
		  stab_idx = stab_idx + 1
		  cell_spike_time = {getfield /stables/pyr_stab_{cell_idx} table->table[{stab_idx}]}
	        end

	        // If -1 or end of table, then last entry < curr_win_start. Else, check if < curr_win_end
	        if ((cell_spike_time == -1) || (stab_idx == (nbins-1)))
		  setfield /curr_stab_idx table->table[{cell_idx}] -1
		  num_ILN_spikes_over = num_ILN_spikes_over + 1
//		  echo last spike reached: cell {cell_idx}, win_start {curr_win_start}
	        else
//		  echo cell {cell_idx}, stime = {cell_spike_time}
		  setfield /curr_stab_idx table->table[{cell_idx}] {stab_idx} // start from same
		  // spike time for next window since it may occur after win_start for next window
		  if (cell_spike_time < curr_win_end) // note curr_win_end is open end of interval
		    if ((cell_idx >= pat_cell_start) && (cell_idx <= pat_cell_end))
			num_fired_in_pattern = num_fired_in_pattern + 1
		    else
			num_fired_out_pattern = num_fired_out_pattern + 1
//			echo out_pattern cell idx = {cell_idx}
		    end
		  end
	        end
	      end // if stab_idx != -1
	    end // for cell_idx

//	    echo num_fired_in_pattern = {num_fired_in_pattern}
//	    echo num_fired_out_pattern = {num_fired_out_pattern}
	    prob_pattern = (num_fired_in_pattern / 5) - (num_fired_out_pattern /(num_ILN-5))
	    if (prob_pattern > max_prob_pattern)
		max_prob_pattern = prob_pattern
	    end
	  end // for nb_idx
	  echo prob_pattern {pattern_idx} == {max_prob_pattern}
	  prob_sequence = prob_sequence + max_prob_pattern
//	  echo prob_sequence = {prob_sequence}
	end // else portion of if num_pat_cell_fired == 0
    end // if (num_ILN_spikes_over < num_ILN)..
  end // for pattern_idx

  RCL_ACC = prob_sequence / (num_pattern-1) // first (cue) pattern not included in measure 

  echo *********************************
  echo RECALL ACCURACY = {RCL_ACC}
  echo *********************************

end

//--------------------------------------------------------------------------------------------------

// compute_bg_perc saves measures of how even the distribution of background firing
// in between two successive patterns is, over the entire sequence during the learning phase. 
// The percentage of bg cells (called bg_perc), that fire in between two successive input patterns
// (this time period called bg_epoch), 
// is computed for all patterns in the sequence and the following measures are derived from it:
// 1. Overall mean(bg_perc) and SD(bg_perc).
// 2. Mean and SD of absolute difference of neighbouring bg_percs.

// Note: uses variable 'nbins' which is set (in make_<pyr/inta/intb>_stabs) after lrs has been called.
// ********Thus, compute_bg_perc should be called only AFTER lrs HAS BEEN CALLED!***********************

function compute_bg_perc

  int   pattern_idx  // index that loops over patterns in the sequence
  int   cnum  // pyr[cnum]
  int	cnum_idx // index into current_stab_idx table corresponding to cnum
  int   stab_idx = 0 // index into pyr_stab_i, got from stab_indices
  int   i
  int	xdivs = {NPYR - num_ILN - 1} // used for table /current_stab_idx.
  int	xmin = 0
  int	xmax = {xdivs}
  int	bg_perc_xdivs = {num_aff_pat - 2}
  int	bg_perc_xmin = 0
  int	bg_perc_xmax = {bg_perc_xdivs}

  float cell_spike_time // value got from pyr_stab_i
  float bg_epoch_start // start time of interval to measure bg_perc
  float bg_epoch_end   // end time of above interval
  float	num_bg_fired
  float	perc_bg_fired
  float	perc_val
  float	diff_perc_val

  // table '/current_stab_idx' stores for each background cell (bg), the current position in the stab table
  // If that bg has no more spikes (-1 reached in stab table), then /current_stab_idx also
  // stores -1 for that cell.
  if (!{exists /current_stab_idx}) 
	create table /current_stab_idx
  end
  call /current_stab_idx TABCREATE {xdivs} {xmin} {xmax} // resets elements to 0

  // table /bg_perc stores percentage of bg cells that fire in between successive aff. patterns
  // bg_perc[0] stores percent firing between patterns 1 (first pattern)  and 2, etc.
  if (!{exists /bg_perc}) 
	create table /bg_perc
  end
  call /bg_perc TABCREATE {bg_perc_xdivs} {bg_perc_xmin} {bg_perc_xmax} // resets elements to 0
  for (i=0; i<={num_aff_pat - 2}; i=i+1) // initialize elements to -1
	setfield /bg_perc table->table[{i}] {-1}
  end

  for (pattern_idx=1; pattern_idx < num_aff_pat; pattern_idx=pattern_idx+1) // first pat. = pat 1

    // derive bg_epoch_start and bg_epoch_end times
    bg_epoch_start	= (pattern_idx-1) * {aff_input_interval} + 0.015 // 15ms latency before aff input
    bg_epoch_end	= bg_epoch_start + aff_input_interval
 
    num_bg_fired = 0

    for (cnum = num_ILN; cnum < NPYR; cnum=cnum+1) // loop over all bc's
	cnum_idx = cnum - num_ILN // since current_stab_idx has indices that run from 0 to {NPYR-num_ILN-1}
	stab_idx = {getfield /current_stab_idx table->table[{cnum_idx}]} 
	if ( stab_idx != -1) // if cell still to fire
		cell_spike_time = {getfield /stables/pyr_stab_{cnum} table->table[{stab_idx}]}
		while ((stab_idx < (nbins-1)) && (cell_spike_time != -1) &&  \
		  (cell_spike_time < bg_epoch_start)) // 
			stab_idx = stab_idx + 1
			cell_spike_time = {getfield /stables/pyr_stab_{cnum} table->table[{stab_idx}]}
		end

		// If -1 or end of table, then no more spikes >= bg_epoch_start
		if ( (cell_spike_time == -1) || ( stab_idx == (nbins-1) ) )
			setfield /current_stab_idx table->table[{cnum_idx}] -1
//			echo no more spikes for pyr {cnum} after {bg_epoch_start}
		else
			if (cell_spike_time < bg_epoch_end) // increment num_fired and stab_idx
			  num_bg_fired = num_bg_fired + 1
			  setfield /current_stab_idx table->table[{cnum_idx}] {stab_idx+1}
			else // read the same spike time for next bin
			  setfield /current_stab_idx table->table[{cnum_idx}] {stab_idx}
			end
		end
	end // if stab_idx != -1
    end  // for cnum

  // calculate percentage and set value in /bg_perc
//  echo {pattern_idx}: numfired = {num_bg_fired}
  perc_bg_fired = num_bg_fired / {NPYR} * 100
  echo {pattern_idx}: percfired = {perc_bg_fired}
  setfield /bg_perc table->table[{pattern_idx-1}] {perc_bg_fired}
  end // for pattern_idx

  // calculate MEAN and SD of perc_bg_fired for all background epochs in the sequence
  mean_perc_bg = 0
  for (i=0; i <= {bg_perc_xmax}; i=i+1)
	mean_perc_bg = mean_perc_bg + {getfield /bg_perc table->table[{i}]}
  end
  mean_perc_bg = mean_perc_bg / (bg_perc_xmax + 1)
  echo mean_perc_bg = {mean_perc_bg}

  // SD = sqrt{ sum[ ( x-mean(x) )^2 ] / (n-1) }
  SD_perc_bg = 0
  for (i=0; i <= {bg_perc_xmax}; i=i+1)
	perc_val = {getfield /bg_perc table->table[{i}]}
	SD_perc_bg = SD_perc_bg + (perc_val - mean_perc_bg)**2
  end
  SD_perc_bg = (SD_perc_bg / bg_perc_xmax)**0.5
  echo SD_perc_bg = {SD_perc_bg}

  // calculate mean and SD of difference between adjacent perc_bg_fired
  mean_diff_perc_bg = 0
  for (i=0; i <= {bg_perc_xmax - 1}; i=i+1)
	diff_perc_val = {abs {{getfield /bg_perc table->table[{i+1}]} - {getfield /bg_perc table->table[{i}]}}}
	mean_diff_perc_bg = mean_diff_perc_bg + diff_perc_val
  end
  mean_diff_perc_bg = mean_diff_perc_bg / (bg_perc_xmax)
  echo mean_diff_perc_bg = {mean_diff_perc_bg}

  // SD = sqrt{ sum[ ( x-mean(x) )^2 ] / (n-1) }
  SD_diff_perc_bg = 0
  for (i=0; i <= {bg_perc_xmax - 1}; i=i+1)
	diff_perc_val = {abs {{getfield /bg_perc table->table[{i+1}]} - {getfield /bg_perc table->table[{i}]}}}
	SD_diff_perc_bg = SD_diff_perc_bg + (diff_perc_val - mean_diff_perc_bg)**2
  end
  SD_diff_perc_bg = (SD_diff_perc_bg / (bg_perc_xmax - 1))**0.5
  echo SD_diff_perc_bg = {SD_diff_perc_bg}
end


//---------------------------------------------------------------------------------------------------

// functions to set pyr AMPA hebbsynchan params

function set_post_thresh_hi(postval)
	float postval
	int i
	for(i = 0; i < {NPYR}; i = i + 1)
	   setfield /pyr_layer/pyr[{i}]/apical/AMPA_channel \
	   post_thresh_hi {postval}
	end
end
// set_post_thresh_hi -45e-3

function set_weight_change_rate(val)
	float val
	int i
	for(i = 0; i < {NPYR}; i = i + 1)
	   setfield /pyr_layer/pyr[{i}]/apical/AMPA_channel \
	   weight_change_rate {val}
	end
end

//---------------------------------------------------------------------------------------

// following functions set AMPA_channel parameters to allow learning only when either pre or postsynaptic
// cell is bursting. c_start and (c_end - 1) specify the range of cells over which the action is taken.

// post_tau set to normal value for non-ILN cells and high value for ILN cells
function set_post_tau(val,c_start,c_end)
	float val
	int   c_start, c_end, i

	if ((c_start < 0) || (c_start >= {NPYR}))
	   echo c_start out of range. Not changing post_tau value.
	   return
	end

	if ((c_end < 0) || (c_end > {NPYR}))
	   echo c_end out of range.
	   if (c_start < {NPYR-1})
		echo Changing post_tau for cells from {c_start} to NPYR {NPYR}
	   	c_end = {NPYR}
	   else
		echo Not changing value of post_tau
		return
	   end
	end

	for (i = c_start; i < c_end; i=i+1)
	   setfield /pyr_layer/pyr[{i}]/apical/AMPA_channel \
	   post_tau {val}
	end
	
end
// set_post_tau 3e-3 {num_ILN} {NPYR} // non-ILN cells
set_post_tau 1e-3 0 {num_ILN} // ILN cells

// pre_thresh_hi set to large value for non-ILN cells and 0.1 for ILN cells.
function set_pre_thresh_hi(val,c_start,c_end) 
	float val
	int   c_start, c_end, i

	if ((c_start < 0) || (c_start >= {NPYR}))
	   echo c_start out of range. Not changing pre_thresh_hi value.
	   return
	end

	if ((c_end < 0) || (c_end > {NPYR}))
	   echo c_end out of range.
	   if (c_start < {NPYR-1})
		echo Changing pre_thresh for cells from {c_start} to NPYR {NPYR}
	   	c_end = {NPYR}
	   else
		echo Not changing value of pre_thresh_hi
		return
	   end
	end

	for (i = c_start; i < c_end; i=i+1)
	   setfield /pyr_layer/pyr[{i}]/apical/AMPA_channel \
	   pre_thresh_hi {val}
	end
	
end
// set_pre_thresh_hi 1.0 {num_ILN} {NPYR}  // non-ILN cells
// set_pre_thresh_hi 0.1 0 {num_ILN}  // ILN cells

// weight_change_rate set to large value for non-ILN cells and lower for ILN cells to prevent saturation of
// within pattern learning

function set_wt_change_rate(val,c_start,c_end) 
	float val
	int   c_start, c_end, i

	if ((c_start < 0) || (c_start >= {NPYR}))
	   echo c_start out of range. Not changing weight_change_rate.
	   return
	end

	if ((c_end < 0) || (c_end > {NPYR}))
	   echo c_end out of range.
	   if (c_start < {NPYR-1})
		echo Changing weight_change_rate for cells from {c_start} to NPYR {NPYR}
	   	c_end = {NPYR}
	   else
		echo Not changing value of weight_change_rate
		return
	   end
	end

	for (i = c_start; i < c_end; i=i+1)
	   setfield /pyr_layer/pyr[{i}]/apical/AMPA_channel \
	   weight_change_rate {val}
	end
	
end
// set_wt_change_rate {non_ILN_learn_rate} {num_ILN} {NPYR}  // non-ILN cells
// set_wt_change_rate {ILN_learn_rate} 0 {num_ILN}  // ILN cells


function set_pyr2pyr_learnrate
	int i

	for (i=0; i < {NPYR}; i=i+1)
	    setfield /pyr_layer/pyr[{i}]/apical/AMPA_channel weight_change_rate {learn_rate}
	end
end

set_pyr2pyr_learnrate

function set_EC2pyr_learnrate
	int i

	for (i=0; i < {NPYR}; i=i+1)
	    setfield /pyr_layer/pyr[{i}]/apical_2/AMPA_channel weight_change_rate {EC2pyr_learnrate}
	end
end

set_EC2pyr_learnrate

function set_EC2pyr_maxweight
	int i

	for (i=0; i < {NPYR}; i=i+1)
	    setfield /pyr_layer/pyr[{i}]/apical_2/AMPA_channel max_weight {EC2pyr_maxweight}
	end
end

set_EC2pyr_maxweight

//---------------------------------------------------------------------------------------

// function to print out weights of AMPA synapses onto ILN cells from other ILN cells

function print_weights
  int cpost, cpre, syn_num
  for(cpost = 0; cpost < {num_ILN}; cpost = cpost + 1)
    for (cpre = 0; cpre <= cpost; cpre = cpre + 1)
	if ({getsyncount /pyr_layer/pyr[{cpre}]/soma/spike /pyr_layer/pyr[{cpost}]/apical/AMPA_channel})
	  syn_num = {getsynindex /pyr_layer/pyr[{cpre}]/soma/spike \
			 	/pyr_layer/pyr[{cpost}]/apical/AMPA_channel}
	  echo Weight of synapse from pyr[{cpre}] to pyr[{cpost}]: \ 
	     {getfield /pyr_layer/pyr[{cpost}]/apical/AMPA_channel \
	               synapse[{syn_num}].weight}
	else
	  echo No synapse between pyr[{cpre}] and pyr[{cpost}]
	end
    end
  end
end

//---------------------------------------------------------------------------------------

// function to print out weights of AMPA synapses that have changed by more than a certain amount.
// Prints out name of presynaptic cell, name of postsynaptic cell and wight of synapse for all p->p
// synapses whose weights are larger than {LARGE_WT}

function print_changed_synapses(w_file)
  int   cpost, imes, i
  str   src
  float weight
  float LARGE_WT = 0.5 // weights greater than this amount are printed to file

  echo Synapses that learn (initial weight = {pyr_ampa_wt}, final weight >= {LARGE_WT})  > {w_file}
  echo learning rate = {learn_rate} >> {w_file}
  echo >> {w_file} // newline

  for(cpost = 0; cpost < {NPYR}; cpost = cpost + 1)
    imes = {getfield /pyr_layer/pyr[{cpost}]/apical/AMPA_channel nsynapses}
    for(i = 0; i < imes; i = i + 1)
	weight = {getfield /pyr_layer/pyr[{cpost}]/apical/AMPA_channel synapse[{i}].weight}
	if (weight > {LARGE_WT})
		src = {getsynsrc /pyr_layer/pyr[{cpost}]/apical/AMPA_channel {i}}
		echo pre: {src} >> {w_file}  
		echo post: /pyr_layer/pyr[{cpost}]/apical/AMPA_channel >> {w_file}
		echo {weight} >> {w_file}  
		echo >> {w_file} // newline
	end
    end
  end
end

//---------------------------------------------------------------------------------------

// For each ILN cell gets (1) sum of weights of AMPA synapses (2) max wt. of AMPA synapse
// Averages sum over ILN cells belonging to same pattern

function get_sum_max_weights(fname)
  int   cpost, nsyn, i, pat, c_start, c_end
  float weight, sum_weight, max_weight

  for (pat = 0; pat < num_aff_pat; pat=pat+1)
    c_start = pat*5
    c_end   = c_start + 4
    sum_weight = 0
    max_weight = -1
    for(cpost = c_start; cpost <= c_end; cpost = cpost + 1)
      nsyn = {getfield /pyr_layer/pyr[{cpost}]/apical/AMPA_channel nsynapses}
      for (i=0; i < nsyn; i=i+1)
	weight = {getfield /pyr_layer/pyr[{cpost}]/apical/AMPA_channel synapse[{i}].weight}
	sum_weight = sum_weight + weight
	if (weight > max_weight)
		max_weight = weight
	end
      end
    end
    echo pat{pat+1}: mean total weight = {sum_weight/5}, max weight = {max_weight} >> {fname}
  end
end

//-----------------------------------------------------------------------------------------


// set_aff_input presents burst of 4 spikes at 100Hz to each cell in each pattern when width1 
// of gategen is 45e-3. To present only single spike instead of a burst set width1 to 15ms.
// Number of patterns is specified in argument num_aff_pat

function set_aff_input(num_aff_pat)
  float num_aff_pat
  int i

  set_gate_delays 0 {aff_input_interval}
  for (i = 0; i < {GATE_NY}; i = i+1)
    setfield /gate_layer/gate[{i}]/gategen width1 15e-3 // (for single spike use 15ms width)
  end
  connect_aff {num_aff_pat}
end

//----------------------------------------------------------------------------------------

// set_recall_cue presents 1 spike (aff. input) to each cell of pattern 1
function set_recall_cue(start_time)
  float start_time //  latency (sec) of recall cue 
  int i

  set_gate_delays {start_time} 1000 // interval very long to prevent activation of succ. patterns
  for (i = 0; i < {GATE_NY}; i = i+1)
    setfield /gate_layer/gate[{i}]/gategen width1 15e-3 // only 1 spike. 10ms latency before 1st spike 
  end
  connect_aff 1
end

//------------------------------------------------------------------------------------------

// set_cyclic_input causes 5th (last)  pattern to be the same as the 1st pattern, so that the
// 4th pattern starts the sequence again

function set_cyclic_input
  int gatenum = 4 // 5th gategen will be reconnected to appropriate pulsegens
  int startaff = 0 // 5th gategen will connect to affs 0 - 4 (which connect to pyr[0-4] resp.)

  change_gate_conn {gatenum} {startaff}
end

//------------------------------------------------------------------------------------------

// set_second_input provides aff input to pyr cells for a second sequence
// Seq 1 = 1->2->3->4->5. 	Seq 2 = 1->6->7->8->9.
// set_second_input connects gategens 1,2,3 and 4 (for patterns 6-9) to
// appropriate pulsegens, and also sets gate delays (gategen 0 -> pat 1).

function set_second_input(second_seq_start_time)
  float second_seq_start_time

  // set gate delays for second input
  set_gate_delays {second_seq_start_time} {aff_input_interval}

  // connect gategens for last 4 patterns to appropriate pulsegens
  change_gate_conn 1 25 // gategen 3 connects to pulsegens 25-29
  change_gate_conn 2 30 // gategen 3 connects to pulsegens 30-34
  change_gate_conn 3 35 // gategen 3 connects to pulsegens 35-39
  change_gate_conn 4 40 // gategen 4 connects to pulsegens 40-44

  // connect pulsegens to pyr cells
  connect_aff 9 // although patterns 2-5 receive connections from pulsegens, those pulsegens
		// are not activated by any gategens.

end

//------------------------------------------------------------------------------------------

// current_inject injects current specified by inject_val into pyr, int_a or int_b cells (specified
// by cell_type)

function inject_current(inject_val, cell_type)
  int   i
  float inject_val
  str   cell_type

  if (cell_type == "pyr")
    for (i=0; i < {NPYR}; i=i+1)
	setfield /pyr_layer/pyr[{i}]/soma inject {inject_val}
    end
  elif (cell_type == "int_a")
    for (i=0; i < {NINTA}; i=i+1)
	setfield /int_a_layer/int_a[{i}]/soma inject {inject_val}
    end
  elif (cell_type == "int_b")
    for (i=0; i < {NINTB}; i=i+1)
	setfield /int_b_layer/int_b[{i}]/soma inject {inject_val}
    end
  else
    echo cell_type should be one of pyr, int_a or int_b.
    echo No action taken.
  end
end

//------------------------------------------------------------------------------------------

// depol_ILN_bg injects current (specified by inject_val) into ILN pyr cells of seq 1 
// (except those of pat 1 (pyr[0-4]) and all bg cells

function depol_ILN_bg(inject_val)
  float inject_val
  int   i
  float delay, width

  if (num_seq == 1)
	delay = learn_time + 10e-3
  elif (num_seq == 2)
	delay = learn_time + gap_time + learn_time + 10e-3
  end
  width = 100e-3 // recall_time - 50e-3

  echo Adding depolarizing pulse of {inject_val} Amps to ILN pyr cells of seq. 1 
  for (i=5; i < {25}; i = i+1) // Note: no current to pat 1
	setfield /pyr_layer/pyr[{i}]/soma/pulse delay1 {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse width1 {width}
	setfield /pyr_layer/pyr[{i}]/soma/pulse level1 {inject_val}
  end
/*
  echo Adding depolarizing pulse of {inject_val} Amps to all bg cells 
  for (i={num_ILN}; i < {NPYR}; i = i+1) // Note: no current to pat 1
	setfield /pyr_layer/pyr[{i}]/soma/pulse delay1 {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse width1 {width}
	setfield /pyr_layer/pyr[{i}]/soma/pulse level1 {inject_val}
  end
*/
end

// depol_ILN_bg {EC_ILN_depol}

//------------------------------------------------------------------------------------------

// hyperpol_ILN injects hyper. current pulse into ILN cells to prevent them from firing once they
// have fired in sequence.

function hyperpol_ILN
  int i
  float inject_val = -100e-9 // -1 nA
  float delay
  float width

  echo Adding hyperpol. pulse of {inject_val} Amps to ILN cells...

  for (i=0; i <= 4 ; i=i+1)
	delay = 80e-3  // (15ms burst latency + 45 ms burst duration)
	width = {learn_time} - {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse delay1 {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse width1 {width}
	setfield /pyr_layer/pyr[{i}]/soma/pulse level1 {inject_val}
  end

  for (i=5; i <= 9 ; i=i+1)
	delay = 180e-3  // (115ms burst latency + 45 ms burst duration)
	width = {learn_time} - {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse delay1 {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse width1 {width}
	setfield /pyr_layer/pyr[{i}]/soma/pulse level1 {inject_val}
  end

  for (i=10; i <= 14 ; i=i+1)
	delay = 280e-3  // (215ms burst latency + 45 ms burst duration)
	width = {learn_time} - {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse delay1 {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse width1 {width}
	setfield /pyr_layer/pyr[{i}]/soma/pulse level1 {inject_val}
  end

  for (i=15; i <= 19 ; i=i+1)
	delay = 380e-3  // (315ms burst latency + 45 ms burst duration)
	width = {learn_time} - {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse delay1 {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse width1 {width}
	setfield /pyr_layer/pyr[{i}]/soma/pulse level1 {inject_val}
  end

  for (i=20; i <= 24 ; i=i+1)
	delay = 480e-3  // (415ms burst latency + 45 ms burst duration)
	width = {learn_time} - {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse delay1 {delay}
	setfield /pyr_layer/pyr[{i}]/soma/pulse width1 {width}
	setfield /pyr_layer/pyr[{i}]/soma/pulse level1 {inject_val}
  end

end

// hyperpol_ILN

//------------------------------------------------------------------------------------------

// quench_activity resets all elements in all cells of the network

function quench_activity
  int i
  str elt_name

  echo resetting all elements in all cells of the network (quench activity) ....

  for (i=0; i < NPYR; i=i+1)
    foreach elt_name ({el /pyr_layer/pyr[{i}]/##})
	call {elt_name} RESET
    end
  end

  for (i=0; i < NINTA; i=i+1)
    foreach elt_name ({el /int_a_layer/int_a[{i}]/##})
	call {elt_name} RESET
    end
  end

  for (i=0; i < NINTB; i=i+1)
    foreach elt_name ({el /int_b_layer/int_b[{i}]/##})
	call {elt_name} RESET
    end
  end

  for (i=0; i < NPYR; i=i+1)
    foreach elt_name ({el /bkgnd_pulse_layer/bkgnd_pulse[{i}]/##})
	call {elt_name} RESET
    end
  end

  for (i=0; i < NEC; i=i+1)
    foreach elt_name ({el /EC_layer/EC[{i}]/##})
	call {elt_name} RESET
    end
  end

end

//----------------------------------------------------------------------------------------

// activate_EC injects a depol. pulse into the appropriate (ecnum) EC cell so that
// it fires for the duration of that sequence and has its connections to all ILN cells
// in that sequence be strengthened.

function activate_EC(ecnum, delay, width, ival)
	int 	ecnum
	float	delay // = 1e-3 // start prior to first pattern of sequence (~15 ms)
	float	width // = {learn_time} - 50e-3 // end little before learning phase ends
	float	ival //= 0.4e-9 // Amps

	setfield /EC_layer/EC[{ecnum}]/soma/pulse delay1 {delay}
	setfield /EC_layer/EC[{ecnum}]/soma/pulse width1 {width}
	setfield /EC_layer/EC[{ecnum}]/soma/pulse level1 {ival}
end

//----------------------------------------------------------------------------------------

// set_pyr_rand_rate sets the value of 'rate' of the randomspike elements of all pyr cells.

function set_pyr_rand_rate(rs_rate)
  int i
  float rs_rate

  for (i=0; i < {NPYR}; i=i+1)
	setfield /bkgnd_rs_layer/bkgnd_rs[{i}]/randspike rate {rs_rate}
  end
end

//------------------------------------------------------------------------------------------

// Following functions save Vm of pyr cells (c_start to c_end)

function do_disk_out(diskpath,srcpath,field)
  str diskpath, srcpath, field

  create asc_file /output/{diskpath}
  setfield /output/{diskpath} leave_open 1  flush 1 notime 1 \
	   filename data/Vm/{diskpath}
  addmsg {srcpath} /output/{diskpath} SAVE {field} 
  useclock /output/{diskpath} 0
  echo {diskpath}
end

function Save_Vm(c_start, c_end)
  int c_start, c_end
  int i

  for (i=c_start; i < c_end; i=i+1)
	do_disk_out pyr_{i}_Vm.dat  /pyr_layer/pyr[{i}]/soma  Vm
  end

end

// Save_Vm 30 31 

//------------------------------------------------------------------------------------------

// create bin table to hold binned activity of pyr cells and int_a and int_b cells
// reset needs to be called after these functions
make_pyr_stabs {tmax}
make_inta_stabs {tmax}
make_intb_stabs {tmax}

reset

//------------------------------------------------------------------------------------------

//*******************************************************************************************
// lrs: main function
// 1. presents afferent input with weight changes allowable to learn sequence of patterns
// 2. Fixes weights and cues retrieval with single shocks to pyr cells of first pattern
// 3. Calls performance measure function, save and plot spike times and/or %active functions.
//*******************************************************************************************

function lrs

  date // call unix command to see start time

  int save_perc_active = 0 // 1 => save perc. activity to file, 0 => dont save to file
  int   learn_nsteps  = {learn_time / {dt}} 
  int   recall_nsteps = {recall_time / {dt}}
  int	gap_nsteps = {gap_time / {dt}} // gap between successive learning phases in 2 seq. task
  float	second_seq_start_time 
  float recall_start_time
  float tot_simtime = {tmax} // total simulation time in sec.
  int   i

  //--------------INITIALIZATIONS----------------------------------------------------------

 // initialize non-hebb weights
  init_nonchg_weights

  // initialize hebb weights
  init_hebb_weights

  // add oscillating depol to bg cells
  add_sawtooth_depol {bg_depol_current} {bg_depol_plateau_period} {bg_depol_ramp_period}

  // reset initVm of bg cells to {bg_steady_depol} and adjust their Em accordingly
//  reset_Em_pyr {num_ILN} {NPYR} {bg_steady_depol}

  // reset initVm of ILN cells to -70 mV and adjust their Em accordingly 
  reset_Em_pyr 0 {num_ILN} -0.07

  // set params for depol. (current pulse) into appropriate EC cell
//  activate_EC 0 1e-3 {learn_time - 50e-3} 0.4e-9

  // initialize all cells in network
  quench_activity

  // free AMPA hebbsynchans on all pyr cells 
  free_hebb_weights

  // Connect afferent input (pulse generators) to fire pyr cells.
  set_aff_input {num_aff_pat} // each pattern is burst of 4 spikes at 100Hz to 5 pyr cells 
  
  // set background firing rate of pyr cells (input from /bkgnd_pulse_layer/bkgnd_pulse[{i}]/pulse) 
  // (make sure set_bkgnd_conn is not commented out (in background_layer.g) for this to work)
//  set_pyr_rand_rate {background_rate}

  //--------------LEARNING PHASE----------------------------------------------------------

  // present afferent input sequence to network and let network learn
  // simulate for duration of learning_time
  echo presenting sequence to network and changing weights...
  echo learning time == {learn_time} sec
  step {learn_nsteps}


  //--------------LEARNING PHASE FOR SECOND SEQUENCE--------------------------------------

  if (num_seq == 2)
    // quench activity of network
    echo Quenching activity of network prior to presentation of second sequence ...
    quench_activity

    step {gap_nsteps}

    // present afferent input sequence to network and let network learn
    second_seq_start_time = {getstat -t}
    set_second_input {second_seq_start_time}

    // set params for depol. (current pulse) into appropriate EC cell
//    activate_EC 1 {second_seq_start_time + 1e-3} {learn_time - 50e-3} 0.4e-9

    // simulate for duration of learning_time
    echo presenting second sequence to network and changing weights...
    echo learning time == {learn_time} sec
    step {learn_nsteps}
  end


  //--------------INITIALIZATIONS BEFORE RECALL PHASE--------------------------------------

  if (do_recall_in_lrs == 1)
    // get current time
    recall_start_time = {getstat -t}

    // print out synapses that have become potentiated (pre, post, weight)
//    print_changed_synapses data/changed_synapses.dat

    // fix weights of all AMPA hebbsynchans for recall
    fix_hebb_weights

    // remove oscillating depol from bg cells
    add_sawtooth_depol 0.0e-9 {bg_depol_plateau_period} {bg_depol_ramp_period} 

    // reduce feedback GABA_A inhibition (reset pyr_gabaa_wt)
    reset_pyr_gaba_weights

    // multiply pyr_ampa_wts of POTENTIATED SYNAPSES ONLY by recall_exc_wt_mult
    reset_pyr_ampa_wt

    // set EC_to_pyr weights of POTENTIATED SYNAPSES ONLY to EC2pyr_recall_maxwt
//    reset_EC2pyr_wt

    // reset_p2p_syndelay(p2pd) changes the delay between p2p connections just prior to recall
    // phase, so that overlapping (jitter) of recalled patterns may be avoided.
//    reset_p2p_syndelay {recall_p2p_syndelay}

    // Reset initVm of all cells to -70 mV and adjust their Em accordingly 
//    reset_Em_pyr 0 {NPYR} -0.07

    // set params for depol. (current pulse) into appropriate EC cell. ecnum, delay, width, val
    activate_EC 0 {recall_start_time} 1e-3 1e-7 // brief, large pulse for 1 spike

    // set weights of EC2pyr connections (from EC[0] to ILN cells of seq. 1)
    set_EC2pyr_weights

    // reset membrane potential of all cells in the network and their subelements (channels etc)
    // so that activity is zero before presentation of recall cue.
    echo Quenching activity of network prior to recall ...
    quench_activity

    // print out post-learning weights of AMPA synapses of ILN cells
//    echo ------------------------------------------------------------- >> data/weights.dat
//    echo LEARNING RATE = {learn_rate}, p-p probability = {p_to_p}      >> data/weights.dat
//    echo ------------------------------------------------------------- >> data/weights.dat
//    get_sum_max_weights data/weights.dat

//    print_weights

    //--------------RECALL PHASE----------------------------------------------------------

    // test recall with weights fixed for duration of recall_time
    // present first pattern and see if sequence generated
    set_recall_cue {recall_start_time + expand_recall_time} // only 1 spike to each cell of pattern 1
    echo presenting first pattern of sequence and eliciting recall ...
    echo recall time == {recall_time} sec
    step {recall_nsteps}


    //--------------ANALYZE/PLOT/SAVE MEASURES----------------------------------------------------------


    // measure recall performance
    // Start from just after the recall cue is given. Since first spike of recall cue occurs
    // 15 ms after recall_start_time (delays set in aff_layer.g), (recall_start_time + 17e-3)
    // will start from first bin after the recall cue.
    gpm3  {recall_start_time + expand_recall_time + 17e-3} {num_aff_pat} 

    // Measure variability of background firing across sequence
    compute_bg_perc

    // SAVE SPIKE TIME DATA
//    save_pyr_stabs data/pyr_stab.dat 25
//    save_inta_stabs data/inta_stab.dat
//    save_intb_stabs data/intb_stab.dat

    // PLOT PERCENT ACTIVE PLOTS AND SAVE IF NECESSARY
//    plot_perc_active {NSTAB_PYR} 0 {tot_simtime} {save_perc_active} 
//    plot_perc_inta_active {NSTAB_INTA} 0 {tot_simtime} {save_perc_active} 
//    plot_perc_intb_active {NSTAB_INTB} 0 {tot_simtime} {save_perc_active} 
 
    // PLOT RASTER PLOTS OF SPIKE ACTIVITY OF CELLS
//    rplot {NSTAB_PYR} {tot_simtime}
//    rplot_inta {NSTAB_INTA} {tot_simtime}
//    rplot_intb {NSTAB_INTB} {tot_simtime}

  end // if (do_recall_in_lrs)

  date // see end time
end  


//-----------------------------------------------------------------------------------------------

// recall runs only recall phase of lrs. This can be run after learning has already occured
// (by running lrs), to assess effects of various parameters on recall dynamics.
 
function recall

  date // call unix command to see start time

  int save_perc_active = 0 // 1 => save perc. activity to file, 0 => dont save to file
  int   recall_nsteps = {recall_time / {dt}}
  float recall_start_time
  float tot_simtime = {recall_time} // total simulation time in sec.
  int   i

  //--------------INITIALIZATIONS BEFORE RECALL PHASE--------------------------------------

  // get current time
  recall_start_time = {getstat -t}
  echo recall start time is {recall_start_time}

  // fix weights of all AMPA hebbsynchans for recall
  fix_hebb_weights

  // remove oscillating depol from bg cells
  add_sawtooth_depol 0.0e-9 {bg_depol_plateau_period} {bg_depol_ramp_period} 

  // set params for depol. (current pulse) into appropriate EC cell  
//  activate_EC 0 {recall_start_time + 10e-3} {recall_time - 50e-3} 0.4e-9

  // Reset initVm of all cells to -70 mV and adjust their Em accordingly 
//  reset_Em_pyr 0 {NPYR} -0.07

  // reset membrane potential of all cells in the network and their subelements (channels etc)
  // so that activity is zero before presentation of recall cue.
  echo Quenching activity of network prior to recall ...
  quench_activity


  // create bin table to hold binned activity of pyr cells and int_a and int_b cells
  // reset needs to be called after these functions
  make_pyr_stabs {tot_simtime}
  make_inta_stabs {tot_simtime}
  make_intb_stabs {tot_simtime}
  reset


  //--------------RECALL PHASE----------------------------------------------------------

  // test recall with weights fixed for duration of recall_time
  // present first pattern and see if sequence generated
  set_recall_cue {recall_start_time + expand_recall_time} // only 1 spike to each cell of pattern 1
  echo presenting first pattern of sequence and eliciting recall ...
  echo recall time == {recall_time} sec
  step {recall_nsteps}

  //--------------ANALYZE/PLOT/SAVE MEASURES----------------------------------------------------------


  // measure recall performance
  // Start from just after the recall cue is given. Since first spike of recall cue occurs
  // 15 ms after recall_start_time (delays set in aff_layer.g), (recall_start_time + 17e-3)
  // will start from first bin after the recall cue.
  gpm3  {recall_start_time + expand_recall_time + 17e-3} {num_aff_pat} 

  // SAVE SPIKE TIME DATA
//  save_pyr_stabs data/pyr_stab.dat 25
//  save_inta_stabs data/inta_stab.dat
//  save_intb_stabs data/intb_stab.dat

  // PLOT PERCENT ACTIVE PLOTS AND SAVE IF NECESSARY
//  plot_perc_active {NSTAB_PYR} 0 {tot_simtime} {save_perc_active} 
//  plot_perc_inta_active {NSTAB_INTA} 0 {tot_simtime} {save_perc_active} 
//  plot_perc_intb_active {NSTAB_INTB} 0 {tot_simtime} {save_perc_active} 
 
  // PLOT RASTER PLOTS OF SPIKE ACTIVITY OF CELLS
//  rplot {NSTAB_PYR} {tot_simtime}
//  rplot_inta {NSTAB_INTA} {tot_simtime}
//  rplot_intb {NSTAB_INTB} {tot_simtime}

  date // see end time
end  


//-----------------------------------------------------------------------------------------------


// set_activity_level (sal) is used to explore how network activity varies as a function of
// synaptic strengths, when no learning is included in the model. The strengths and probability
// of connections for AMPA, NMDA, GABAa and GABAb synapses on pyr, int_a and int_b cells
// are varied and percenage of pyr cells active in every bin is calculated.

function sal
  float tot_simtime = {tmax} // duration of simulation
  date // call unix command to see start time
  int save_perc_active = 1 // 1 => save perc. activity to file, 0 => dont save to file
  float background_rate = 20 // 20Hz rate should make 10% of all pyr cells fire at any time

  // create bin table to hold binned activity of all pyr, int_a and int_b cells
  make_pyr_stabs  {tot_simtime} // (calls reset)
//  make_inta_stabs {tot_simtime}
//  make_intb_stabs {tot_simtime}

  // initialize non-hebb weights
  init_nonchg_weights

  // initialize hebb weights and dont allow them to change
  init_hebb_weights
  fix_hebb_weights
 
  // Connect afferent input (pulse generators) to fire 5 pyr cells once only.
  // This is similar to presenting recall cue.
//  set_recall_cue 0 
//  set_aff_input 3

  // set background firing rate
  set_pyr_rand_rate {background_rate}

  // run simulation for duration of tot_simtime
  echo running simulation for duration {tot_simtime} sec ...
  step {tot_simtime} -t

  // plot percent active of each cell type
  plot_perc_active {NSTAB_PYR} 0 {tot_simtime} {save_perc_active}
//  plot_perc_inta_active {NSTAB_INTA} 0 {tot_simtime} 
//  plot_perc_intb_active {NSTAB_INTB} 0 {tot_simtime} 
  
  date // see end time

end

//---------------------------------------------------------------------------------------------