// $Id: perf_measure.g,v 16.1 2000/03/02 14:43:39 ajay Exp $

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


// 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 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 + 0.002 // 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

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


function g_p_m2(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   prob_pattern_evaluated = 0 // set to 1 when prob_pattern calculated
  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
  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 xdivs = {num_ILN-1} // used for table /curr_stab_idx.
  float xmin = 0
  float xmax = {xdivs}
  float prev_win_start // start time of previous window
  float curr_win_start // start time of current window
  float curr_win_end   // end time of current window
  int   length_idx = 1
  float win_length = bin_width * length_idx    // duration of window
  float win_separation = win_length // gap between start (or end) times of windows
  float max_rcl_acc = 0 // maximum over various values of win_length

  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

  for (length_idx = 1; length_idx <= 6; length_idx = length_idx + 1)
  call /curr_stab_idx TABCREATE {xdivs} {xmin} {xmax} // resets elements to 0
  stab_idx = 0
  num_ILN_spikes_over = 0
  prob_sequence = 0
  prob_pattern = 0 
  win_length = bin_width * length_idx
  win_separation = win_length
//  echo --------------------------------------------------------------
//  echo win_length = {win_length}
//  echo -------------------------------------------------------------- 

  prev_win_start = recall_start_time - win_separation // so that curr_win_start = 0

  for (pattern_idx=1; pattern_idx < num_pattern; pattern_idx=pattern_idx+1) // start after cue pattern
    if (num_ILN_spikes_over < num_ILN) // some ILNs have yet to fire
      prob_pattern_evaluated = 0
      while (!{prob_pattern_evaluated}) // get prob_pattern

	curr_win_start = prev_win_start + win_separation
	curr_win_end = curr_win_start + win_length
	prev_win_start = curr_win_start

//	echo curr_win_start {curr_win_start}
//	echo curr_win_end {curr_win_end}

	pat_cell_start = pattern_idx * 5 // 5 cells per pattern
	pat_cell_end = pat_cell_start + 4 
	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
		    end
		  end
	      end
	  end // if last spike not reached
	end // for

	if ((num_fired_in_pattern == 0) && (num_fired_out_pattern == 0)) // ignore this window unless
									// no more ILNs left to fire
 	  if (num_ILN_spikes_over == num_ILN) // no more ILNs left to fire
//		echo prob_pattern {pattern_idx} == 0
//		echo prob_sequence = {prob_sequence}
		prob_pattern_evaluated = 1
	  else // ignore window and set prev_win_start so next win starts at curr_win_end
		echo No ILN fired in interval [{curr_win_start}, {curr_win_end})
		prev_win_start = curr_win_end - win_separation
	  end

	else // add prob_pattern to prob_sequence
//	  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)) 
	  prob_sequence = prob_sequence + prob_pattern
//	  echo prob_pattern {pattern_idx} == {prob_pattern}
//	  echo prob_sequence = {prob_sequence}
	  prob_pattern_evaluated = 1
	end
      end // while
    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 

  if (RCL_ACC > max_rcl_acc)
    max_rcl_acc = RCL_ACC
  end

//  echo *********************************
//  echo RECALL ACCURACY = {RCL_ACC}
//  echo *********************************
  end // for length_idx
  echo max recall accuracy = {max_rcl_acc}
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 intervals of duration (5 * win_length) ms, (5 ms is bin_width
// used to store spike data).
// 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 g_p_m(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   num_ILN = num_pattern * 5 // number of pyr cells coding input
  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
  int   num_fired_in_pattern  // num. of ILNs firing that are part of current pattern
  int   num_fired_out_pattern // num. of ILNs firing that are not in current pattern
  int   prob_pattern_evaluated = 0 // set to 1 when prob_pattern calculated
  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
  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 xdivs = {num_ILN-1} // used for table /curr_stab_idx.
  float xmin = 0
  float xmax = {xdivs}
  float prev_win_start // start time of previous window
  float curr_win_start // start time of current window
  float curr_win_end   // end time of current window
  float win_separation = bin_width * 2 // 10 ms gap between start (or end) times of windows
  float win_length = bin_width * 6    // 30 ms duration of window


  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

  prev_win_start = -win_separation // so that curr_win_start = 0

  for (pattern_idx=1; pattern_idx < num_pattern; pattern_idx=pattern_idx+1) // start after cue pattern
    if (num_ILN_spikes_over < num_ILN) // some ILNs have yet to fire
      prob_pattern_evaluated = 0
      while (!{prob_pattern_evaluated}) // get prob_pattern

	curr_win_start = prev_win_start + win_separation + recall_start_time
	curr_win_end = curr_win_start + win_length
	prev_win_start = curr_win_start
	pat_cell_start = pattern_idx * 5 // 5 cells per pattern
	pat_cell_end = pat_cell_start + 4 
	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}
	      else
	//	  echo 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
		    end
		  end
	      end
	  end // if last spike not reached
	end // for

	if ((num_fired_in_pattern == 0) && (num_fired_out_pattern == 0)) // ignore this window unless
									// no more ILNs left to fire
 	  if (num_ILN_spikes_over == num_ILN) // no more ILNs left to fire
		echo prob_pattern {pattern_idx} == 0
		echo prob_sequence = {prob_sequence}
		prob_pattern_evaluated = 1
	  else // ignore window and set prev_win_start so next win starts at curr_win_end
		echo No ILN fired in interval [{curr_win_start}, {curr_win_end})
		prev_win_start = curr_win_end - win_separation
	  end

	else // add prob_pattern to prob_sequence
	  prob_pattern = (num_fired_in_pattern / 5) - (num_fired_out_pattern /(num_ILN-5)) 
	  prob_sequence = prob_sequence + prob_pattern
	  echo prob_pattern {pattern_idx} == {prob_pattern}
	  echo prob_sequence = {prob_sequence}
	  prob_pattern_evaluated = 1
	end
      end // while
    end // if (num_ILN_spikes_over < num_ILN)..
  end // for

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

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

end

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


//*********************OLD METHOD**********************************************************

// 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 5ms intervals. (This is same as refractory period
// of pyr cells (and ints), so a cell cant fire more than once in a bin).
// RA is calculated by looking at activities of input layer neurons (ILNs), 
// which are the pyr cells receiving afferent input (pyr[0-74]), and which
// code the patterns in the sequnce. If even one ILN fires in a bin, the
// number of pyr cells, out of 5 in a pattern, that are supposed to fire
// in that bin are determined. (For e.g., pyr[0-4] are supposed to fire in
// first bin, pyr[5-9] in bin 2, etc.) This number, divided by 5, gives
// 'prob_pattern' value. The sum of prob_pattern over all the patterns 
// (called 'prob_sequence'), divided by the pattern length (typically 15),
// gives the RCL_ACC value.


function get_perf_measure(recall_start_time, num_pattern)
  float recall_start_time // RCL_ACC measured after learning phase completed
  int   num_pattern    // number of patterns in sequence
  float prob_sequence = 0  // sum of prob_pattern over all patterns in sequence
  float prob_pattern = 0 // (# appropriate cells firing) / 5
  int   pattern_idx  // index that loops over patterns in the sequence
  int   cell_idx  // index that loops over appropriate cells in a pattern
  int   cell_start  // start value of cell_idx in a pattern
  int   cell_end  // end value of cell_idx in a pattern
  int   bin_idx = 0 // keeps track of bin number   
  float bin_start // start time of bin. closed end.
  float bin_end // open end
  int   num_fired_in_pattern = 0 // used to calculate prob_pattern
  int   prob_pattern_evaluated = 0 // set to 1 when prob_pattern calculated
  float cell_spike_time // value got from pyr_stab_i
  int   stab_idx = 0 // index into stab_i, got from stab_indices
  int   total_fired // total number of ILNs (out of 0-74) firing in a particular bin
  int   i
  float xdivs = {NSTAB_PYR-1} // used for table /last_spike_reached. This stores 1 if last spike
			  // has been reached for that cell (-1 encountered in stab).
  float xmin = 0
  float xmax = {xdivs}
  int more_pyr_to_fire = 1

  echo Measuring performance on sequence recall from time {recall_start_time} ....

  if (!{exists /last_spike_reached}) // used to stop searches if known that no more spikes 
	create table /last_spike_reached
  end
  call /last_spike_reached TABCREATE {xdivs} {xmin} {xmax} // resets elements to 0
  for (i={xmin}; i<={xmax}; i=i+1) // initialize elements to 0
	setfield /last_spike_reached table->table[{i}] {0}
  end 

  bin_idx = 0
  more_pyr_to_fire = 1
  for (pattern_idx=0; pattern_idx < num_pattern; pattern_idx=pattern_idx+1)
    prob_pattern_evaluated = 0
    while (!{prob_pattern_evaluated}) // get prob_pattern
	bin_start = bin_idx * bin_width + recall_start_time
	bin_end = bin_start + bin_width
	cell_start = pattern_idx * 5 // 5 cells per pattern
	cell_end = cell_start + 4 
	num_fired_in_pattern = 0

//	echo bin_start = {bin_start}
//	echo bin_end = {bin_end}

	for (cell_idx = cell_start; cell_idx <= cell_end; cell_idx=cell_idx+1)
	  // get first spike time >= bin_start from /stables/pyr_stab_{cell_idx}
	  // Do this only if cell's last spike time is not < bin_start
	  if ({getfield /last_spike_reached table->table[{cell_idx}]} == 0) // last spike not reached  
	      stab_idx = 0
	      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 < bin_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 < bin_start. Else, check if < bin_end
	      if ((cell_spike_time == -1) || (stab_idx == (nbins-1)))
		  setfield /last_spike_reached table->table[{cell_idx}] 1
	//	  echo last spike reached: cell {cell_idx}
	      else
	//	  echo stime = {cell_spike_time}
		  if (cell_spike_time < bin_end) // note bin_end is open end of interval
			num_fired_in_pattern = num_fired_in_pattern + 1
		  end
	      end
	  end // if last spike not reached
	end // for

	prob_pattern = num_fired_in_pattern / 5
	echo prob_pattern {pattern_idx} == {prob_pattern}
	if (prob_pattern == 0) // check if atleast one ILN fired unless last_spike_reached
				// condition is true for all stabs (cells)
	  if (more_pyr_to_fire) // all tables not reached -1

	     total_fired = {get_num_fired {NSTAB_PYR} {bin_start} {bin_end}} // num ILNs firing in bin
	     if (total_fired > 0) // prob_pattern taken as 0
		prob_pattern_evaluated = 1
		bin_idx = bin_idx + 1
	     else // No ILN fired. Skip to next bin for calculation of prob_pattern
		echo No ILN fired in interval [{bin_start}, {bin_end})
		bin_idx = bin_idx + 1
	     end

	  else // no more spikes. set prob_pattern_evaluated

	     echo No more spikes fired by any ILN
	     prob_pattern_evaluated = 1
	     bin_idx = bin_idx + 1

	  end

	else // add prob_pattern to prob_sequence

	  prob_sequence = prob_sequence + prob_pattern
	  prob_pattern_evaluated = 1
	  bin_idx = bin_idx + 1

	end
    end // while
    echo prob_sequence = {prob_sequence}
  end // for

  RCL_ACC = (prob_sequence - 1) / num_pattern  // first pattern was cue for retrieval

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

end

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