// $Id: bin_intb_spikes.g,v 16.1 2000/03/02 14:43:39 ajay Exp $
// functions to:
// 1. Bin spiking activity of int_b (GABAb) cells
// 2. Plot the percent of active intb cells in a bin and 
// 3. Plot (raster) spikes of int_b cells whose activity has been binned. 
// The spike time of each int_b cell is stored in a table (intb_stab) that operates in
// TAB_SPIKE mode (mode 4).

//----------------------------------------------------------------------------------------
// global variables

int	NSTAB_INTB = {NINTB} // activity of int_b[0 - {NSTAB_INTB - 1}] binned into tables

// bin_width (sec) gives min. interval between spikes. 
if (!{exists bin_width})
	float   bin_width = {getfield /int_b_layer/int_b[0]/soma/spike abs_refract}
end

// nbins gives #bins into which spike activity is binned. It is (tot_simtime / bin_width)
// and set in make_intb_stabs
if (!{exists nbins})
  int	nbins  
end

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

function make_intb_stabs(tot_simtime) // create spike time tables for intb cells

  float tot_simtime  // number in sec giving total simulation time
  int xelts  = {tot_simtime}/{bin_width}  // enough bins for total sim
  int xdivs  = {xelts - 1} 
  float xmin = 0
  float xmax = {xdivs}
  int i, j

  nbins = {xelts}
  if (!{exists /stables})
    create neutral /stables
  end
  for (i=0; i<{NSTAB_INTB}; i=i+1)
    if (!{exists /stables/intb_stab_{i}}) 
      create table /stables/intb_stab_{i}
    end
    call /stables/intb_stab_{i} TABCREATE {xdivs} {xmin} {xmax} // resets elements to 0
    for (j={xmin}; j<={xmax}; j=j+1) // initialize elements to -1
	setfield /stables/intb_stab_{i} table->table[{j}] {-1}
    end
    setfield /stables/intb_stab_{i} step_mode 4 stepsize 0.5 // stepsize is spike threshold
					        // output_amp of spikegen == 1

    if (!{getmsg /stables/intb_stab_{i} -in -count})
	addmsg /int_b_layer/int_b[{i}]/soma/spike   /stables/intb_stab_{i} INPUT state
    end

  end
  echo created {NSTAB_INTB} bin tables for int_b cells. Number of bins per cell = {xelts}  
//  reset  
end

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

function save_intb_stabs(fname) // save spike time tables to a file in dir data.
  int i,j

  echo Saving spike time info to {fname}...
  echo spike time info about input layer intb cells > {fname}
  echo number of intb cells == {NSTAB_INTB} >> {fname}
  echo number of bins for each cell == {nbins} >> {fname}
  echo *************************************************** >> {fname}
  echo >> {fname} // newline

  for (i=0; i<{NSTAB_INTB}; i=i+1)
    for (j=0; j<{nbins}; j=j+1)
	echo  {getfield /stables/intb_stab_{i} table->table[{j}]} -n >> {fname}
	echo "  " -n >> {fname} 
    end
    echo >> {fname} // newline
    echo >> {fname} // newline
  end
end

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

// get_num_intb_fired returns the number of ILNs (see below) firing between time
// bin_start and bin_end (including bin_start and excluding bin_end).
// search is done over intb[0] - intb[numstab-1]. A spike time table (/stables/intb_stab_i)
// should exist for each of these cells. Also, bin interval (end - start) should
// be <= abs_refract of intb spikegen, so that only one spike per cell per bin.

function get_num_intb_fired(numstab, bin_start, bin_end)
  int    numstab // looks for activity of intb[0] to intb[numstab-1]
  float  bin_start // start of closed interval (sec)
  float  bin_end  // end of open interval (sec)
  int    cell_idx // loops over ILNs
  int    num_fired = 0  // number of ILNs that fired in that bin interval
  int    stab_idx // index into stab
  int	 more_intb_to_fire
  float  spike_time // retrieved from stab

  more_intb_to_fire = 0 // assume no cells fire after bin_start. set to 1 if later spike_time found

  for (cell_idx = 0; cell_idx < numstab; cell_idx=cell_idx+1)
    if ({exists /stables/intb_stab_{cell_idx}})
      // get first spike time >= bin_start from /stables/intb_stab_{cell_idx}
      stab_idx = 0
      spike_time = {getfield /stables/intb_stab_{cell_idx} table->table[{stab_idx}]}
      while ((stab_idx < (nbins-1)) && (spike_time != -1) && (spike_time < bin_start))
	stab_idx = stab_idx + 1
	spike_time = {getfield /stables/intb_stab_{cell_idx} table->table[{stab_idx}]}
      end
    
      // Check if spike_time < bin_end only if not reached -1 or end of table
      if ((spike_time != -1) && (stab_idx < (nbins-1)))
	  // spike time that is later than bin_start found. Thus, set more_intb_to_fire to 1
	  more_intb_to_fire = 1
	  if (spike_time < bin_end) // note bin_end is open end of interval
	    num_fired = num_fired + 1
	  end
      end
    else // intb_stab_{cell_idx} does not exist
      echo /stables/intb_stab_{cell_idx} does not exist
    end // if exists
  end // for
//  echo total num fired in {bin_start} to {bin_end} = {num_fired}
  return {num_fired}
end

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

// plot_perc_intb_active stores the percentage of numstab intb neurons that are active in each
// bin from start_time to end_time in a table, 'perc_active_intb_tab', and then plots the
// contents of the table vs time. 
// bin_width is the global var bin_width (= abs_refract of intb spikegens).
// max number bins is 'nbins' (= total_simtime / bin_width). These had been created in make_intb_stabs 
// Activity of int_b[0] to int_br[numstab-1] is determined, and perc. of active cells
// is (#active/numstab)*100 for each bin.

function plot_perc_intb_active(numstab, start_time, end_time, save_perc_active)

  int    numstab
  float  start_time 
  float  end_time
  float  num_fired // number of intb cells firing in a bin
  float  perc_fired // percentage of cells firing in a bin
  float  bin_start // passed to get_num_intb_fired
  float  bin_end   // passed to get_num_intb_fired
  float  xdivs = {nbins - 1}
  float  xmin = 0
  float  xmax = {xdivs}
  float  spike_time // retrieved from stab
  int    start_time_idx = 0 // loops over time
  int    end_time_idx = {end_time} / {bin_width}
  int    tab_idx  // index in 'perc_active_intb_tab'
  int    cell_idx // loops over all cells, 0 to numstab-1
  int    stab_idx // index in stables/stab_i
  int    i
  int    num_cells_spikes_over = 0 // number of cells that have finished firing
  float  graph_min, graph_max // to set graph axes
  float  xval, yval // used to set points in xplot

  if (numstab > {NSTAB_INTB})
	echo Error: numstab cannot be larger than {NSTAB_INTB}
	return
  end
  if (end_time_idx > nbins)
	echo Error: end_time cannot be larger than {tmax}
	return
  end
  echo calculating percent active intb cells out of intb[0 - {numstab-1}] ....

  if ({save_perc_active} == 1)
    echo Saving percent active data to data/perc_intb_active.dat ...
    echo # Percent pyr cells active in each bin > data/perc_intb_active.dat
    echo # number of pyr cells == {NSTAB_PYR} >> data/perc_intb_active.dat
    echo # number of bins for each cell == {nbins} >> data/perc_intb_active.dat
    echo # pyr_ampa_wt = pyr_nmda_wt = {pyr_ampa_wt} >> data/perc_intb_active.dat
    echo # *************************************************** >> data/perc_intb_active.dat
    echo # >> data/perc_intb_active.dat // newline
  end

  // create table to store percentages
  if (!{exists /perc_active_intb_tab})
	create table /perc_active_intb_tab
  end
  call /perc_active_intb_tab TABCREATE {xdivs} {xmin} {xmax} // resets elements to 0
  for (i={xmin}; i<={xmax}; i=i+1) // initialize elements to -1
	setfield /perc_active_intb_tab table->table[{i}] {-1}
  end

  // table '/curr_stab_idx' stores for each cell (0 to numstab-1), the current position 
  // in the stab table. If that cell 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 {numstab-1} {0} {numstab-1} // resets elements to 0   

  // get num fired in each bin, calculate perecntage, and store in table
  for (tab_idx = start_time_idx; tab_idx < end_time_idx; tab_idx = tab_idx + 1)
    if (num_cells_spikes_over < numstab) // some cells yet to fire
	bin_start = start_time + tab_idx * bin_width
	bin_end   = bin_start + bin_width
	num_fired = 0
	for (cell_idx = 0; cell_idx < numstab; cell_idx=cell_idx+1)
	  if ({exists /stables/intb_stab_{cell_idx}})
	    // get first spike time >= bin_start from /stables/intb_stab_{cell_idx}
	    // Do this only if cell has not already stopped firing
	    stab_idx = {getfield /curr_stab_idx table->table[{cell_idx}]} 
	    if (stab_idx != -1) // cell still to fire
		spike_time = {getfield /stables/intb_stab_{cell_idx} table->table[{stab_idx}]}
		while ((stab_idx < (nbins-1)) && (spike_time != -1) && (spike_time < bin_start))
		  stab_idx = stab_idx + 1
		  spike_time = {getfield /stables/intb_stab_{cell_idx} table->table[{stab_idx}]}
		end

		// Check if spike_time < bin_end only if not reached -1 or end of table
		if ((spike_time == -1) || (stab_idx == (nbins-1)))
		  setfield /curr_stab_idx table->table[{cell_idx}] -1
		  num_cells_spikes_over = num_cells_spikes_over + 1
		else
		  if (spike_time < bin_end) // increment num_fired and stab_idx
			num_fired = num_fired + 1
			setfield /curr_stab_idx table->table[{cell_idx}] {stab_idx+1}
		  else // read the same spike time for next bin
			setfield /curr_stab_idx table->table[{cell_idx}] {stab_idx}
		  end
		end
	    end // if (stab_idx != -1) ..
	  else // intb_stab_{cell_idx} does not exist
	    echo /stables/intb_stab_{cell_idx} does not exist
	  end // if exists
	end // for (cell_idx ...)
    else // if (num_cell_spikes_over < numstab)...)
    num_fired = 0
    end
    perc_fired = {{num_fired} / {numstab} * 100}
    setfield /perc_active_intb_tab table->table[{tab_idx}] {perc_fired}
  end

  // create xgraph and plot percentages
  graph_min = {start_time}
  graph_max = {end_time}

  if ({exists /perc_active_intb_form})
	delete /perc_active_intb_form
  end

  create xform /perc_active_intb_form [225, 90, 350, 350]

  create xgraph /perc_active_intb_form/g  -hgeom 100% -title "Percent Intb. Cells Active"
  setfield ^ XUnits sec YUnits "% Active"
  setfield ^ xmin {graph_min} xmax {graph_max}
  setfield ^ ymin {-5} ymax {100}

  create xplot /perc_active_intb_form/g/p
  setfield ^ fg blue
  setfield /perc_active_intb_form/g/x_axis fg blue
  setfield /perc_active_intb_form/g/x_axis textcolor blue
  setfield /perc_active_intb_form/g/y_axis fg blue
  setfield /perc_active_intb_form/g/y_axis textcolor blue

  for (i = start_time_idx; i < end_time_idx; i = i + 1)
	xval = start_time + (i * bin_width)
	yval = {getfield /perc_active_intb_tab table->table[{i}]}
	call /perc_active_intb_form/g/p ADDPTS {xval} {yval} 
	if ({save_perc_active} == 1)
  		echo {xval}	{yval} >> data/perc_intb_active.dat
	end
  end

  xshow /perc_active_intb_form

end

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


// rplot_intb plots the spike times of intb cells [0-{numstab-1}] from time 0 to {end_time}

function rplot_intb(numstab, end_time)
  int   numstab
  float end_time
  int   cell_idx, stab_idx, xpts_idx
  int   last_bin
  float next_spike_time
  float bin_end_time
  float xval, yval
  float graph_height

  if (end_time > (bin_width * nbins + bin_width))
	echo end_time should be less than or equal to {bin_width*nbins}
	return
  else // compute last bin number to be plotted
	last_bin = (end_time / bin_width) - 1
  end 

  if ({exists /raster_form_intb})
	delete /raster_form_intb
  end

  create xform /raster_form_intb [225, 50, 700, 700]
  create xgraph /raster_form_intb/g -hgeom 100%
  setfield ^ title "spiking of intb[0 - {numstab-1}]"  
  setfield ^ XUnits sec YUnits "intb[#]"
  setfield ^ xmin 0  xmax {end_time}
  setfield ^ ymin 0  ymax {numstab}

  for (cell_idx = 0; cell_idx < {numstab}; cell_idx = cell_idx + 1)
    setfield /raster_form_intb/g yoffset 1
    create xplot /raster_form_intb/g/p_{cell_idx}
    setfield ^ fg blue
    setfield /raster_form_intb/g/x_axis fg blue
    setfield /raster_form_intb/g/x_axis textcolor blue
    setfield /raster_form_intb/g/y_axis fg blue
    setfield /raster_form_intb/g/y_axis textcolor blue

    // fill interpol-structs xpts and ypts of xplot
    stab_idx = 0
    next_spike_time = {getfield /stables/intb_stab_{cell_idx} table->table[{stab_idx}]}
    for (xpts_idx = 0; xpts_idx <= last_bin; xpts_idx = xpts_idx + 1)
	xval = xpts_idx * bin_width
	bin_end_time = xval + bin_width
	if ((next_spike_time == -1) || (next_spike_time >= bin_end_time)) // open end interval
		yval = 0
	else
		yval = 1
		stab_idx = stab_idx + 1
		next_spike_time = {getfield /stables/intb_stab_{cell_idx} table->table[{stab_idx}]}
	end
	call /raster_form_intb/g/p_{cell_idx} ADDPTS {xval} {yval} 
    end // for xpts_idx
  end // for cell_idx
  xshow /raster_form_intb

end

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