/*score.c - Do scoring and rejection for the fastshft.c experiment, by
  replacing the event table in a GnuroScan file that has been collected using
  fastshft.c.
  Copyright (c) 1996 Matthew Belmonte

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2
  of the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

  If you find this program useful, please send mail to Matthew Belmonte.
  <mkb4@Cornell.edu>.  If you base a publication on data processed by this
  program, please notify Matthew Belmonte and include the following citation
  in your publication:

	Matthew Belmonte, `A Software System for Analysis of
	Steady-State Evoked Potentials', Association for Computing
	Machinery SIGBIO Newsletter 17:1:9-14 (April 1997).
*/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <errno.h>
#include <string.h>
#ifdef MSDOS
#  define strcasecmp(s, t) strcmpi(s, t)
#  include <memory.h>
#endif
#include <math.h>
#include "fastshft.h"
#include "median.h"
#include "queue.h"
#include "intel.h"
#include "gnuro.h"

static char EVENT_FILE[] = "events.tmp";

static void swap(b, i, j)
short *b;
int i, j;
  {
  register short temp;
  temp = b[i];
  b[i] = b[j];
  b[j] = temp;
  }

/*Find the median of an array 'b' of short integers.  This is accomplished using
  partial quicksort; thus the contents of 'b' end up as a permutation of the
  original sequence.
  This runs in linear expected time (cN + cN/2 +cN/4 + cN/8 + ... = 2cN).  The
  key is deciding what parts of the quicksort really need to be done (see the
  invariant below).
*/
short median(b, size)
short *b;
int size;
  {
  register int h, k, m, n;
  m = 0;
  n = size-1;
  /*inv: m <= size/2 < n and b[m..n] is a permutation of the subsequence
    b[m..n] of sorted sequence b[0..size-1]*/
  do
    {
    h = m+1;
    k = n;
    /*inv: b[m+1..h-1] <= b[m] and b[k+1..n] > b[m]*/
    while(h <= k)
      {
      while((h <= k) && (b[h] <= b[m]))
        h++;
      while((h <= k) && (b[k] > b[m]))
        k--;
      if(h < k)
        swap(b, h++, k--);
      } /*elihw*/
    /*k=h-1*/
    if(k==m)
      m++;
    else if(h==n+1)
      swap(b, m, n--);
    else
      {
      swap(b, m, k);
      if(k <= size/2)
        m=k;
      else
        n=k;
      } /*fi*/
    /*b[k] is now in sorted position*/
    } while(k != size/2);
  return(b[k]);
  }

void main(argc, argv)
int argc;
char *argv[];
  {
/*counters*/
  register int sample, channel, cal;
  register long event;
/*registers*/
  register short stimtype;
  register int attended_target;
  register long prev_offset, offset_of_last_scored_response;
  register long window_start, window_end, next_window_end;
/*constants*/
  register long response_window_start, response_window_end;
  register int num_channels, srate;
  register double heog_scale;
  double heog_threshold_angle;
/*files*/
  FILE *cnt, *cnt2, *evt;
/*eye channel assignments*/
  int heog_channel, blink_channel;
/*rejection thresholds and centre values*/
  long blink_threshold, heog_threshold;
/*structure sizes*/
  int blink_filter_len, blink_deriv_len, heog_filter_len, heog_deriv_len;
/*statistics*/
  long num_events, num_blink_rejects, num_saccade_rejects;
  int num_hits, num_misses, num_false_alarms, num_glitches, num_dropped_responses, num_heog_calibrations, num_good_calibrations;
  float summed_latencies, summed_squared_latencies;
/*temporaries*/
  long temp;
/*buffers*/
  long heog_calibration_offset[RE_NUM_HEOG_CALIBS];
  short heog_calibration_value[RE_NUM_HEOG_CALIBS];
  btree *blink_filter, *heog_filter, *blink_window;
  QUEUE *blinkq, *heog_window;
  char	erp[packed_sizeof_SETUP],
	channel_hdr[packed_sizeof_ELECTLOC],
	evt_tbl_hdr[packed_sizeof_TEEG],
	evtbuf[packed_sizeof_EVENT1],
	evtbuf2[packed_sizeof_EVENT2];

  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  if((argc != 2) && (argc != 3))
    {
    fprintf(stderr, "Usage: %s <file.cnt> [HEOG rejection threshold (degrees)]\n", *argv);
    exit(1);
    }
  heog_threshold_angle = ((argc==3)? atof(argv[2]): RE_SACC_THRESH);
  if((cnt = fopen(argv[1],
#ifdef MSDOS
    "rb+"
#else
    "r+"
#endif
    )) == NULL)
    {
    perror(argv[1]);
    exit(errno);
    }
/*read the main header*/
  if(fread(erp, packed_sizeof_SETUP, 1, cnt) != 1)
    {
    perror(argv[1]);
    fclose(cnt);
    exit(errno);
    }
  num_channels = S_nchannels(erp);
  srate = S_rate(erp);

/*determine channel assignments, and set blink_threshold*/
  blink_threshold = -32768;
  heog_channel = blink_channel = -1;
  for(channel = 0; channel != num_channels; channel++)
    {
    fread(channel_hdr, packed_sizeof_ELECTLOC, 1, cnt);
    if(strcasecmp("heog", channel_hdr) == 0)
      {
      heog_channel = channel;
      heog_scale = ((double)(EL_sensitivity(channel_hdr)))*((double)(EL_calib(channel_hdr)))/204.8;
      }
    else if((strcasecmp("blink", channel_hdr) == 0) || ((blink_channel == -1) && ((strcasecmp("veog", channel_hdr) == 0) || (strcasecmp("upe", channel_hdr) == 0))))
      {
      blink_channel = channel;
      blink_threshold = (long)(204.8*RE_BLINK_THRESH/(EL_sensitivity(channel_hdr)*EL_calib(channel_hdr)));
      }
    }
  if(blink_channel == -1)
    {
    fprintf(stderr, "Couldn't find BLINK or UPE channel.\n");
    fclose(cnt);
    exit(1);
    }
  if(heog_channel == -1)
    {
    fprintf(stderr, "Couldn't find HEOG channel.\n");
    fclose(cnt);
    exit(1);
    }

/*read the event table header*/
  fseek(cnt, S_EventTablePos(erp), SEEK_SET);
  if(fread(evt_tbl_hdr, packed_sizeof_TEEG, 1, cnt) != 1)
    {
    perror(argv[1]);
    fclose(cnt);
    exit(errno);
    }
/*verify that this is the raw version of the event table structure*/
  if(T_Teeg(evt_tbl_hdr) != TEEG_EVENT_TAB1)
    {
    fprintf(stderr, "Task data have already been merged.\n");
    fclose(cnt);
    exit(1);
    }
  fseek(cnt, T_Offset(evt_tbl_hdr), SEEK_CUR);

/*compute window lengths*/
  blink_deriv_len = (int)(((long)RE_BLINK_DERIV_LEN)*((long)srate)/1000L);
  blink_filter_len = (int)(((long)RE_BLINK_FILTER_LEN)*((long)srate)/1000L);
  heog_deriv_len = (int)(((long)RE_SACC_DERIV_LEN)*((long)srate)/1000L);
  heog_filter_len = (int)(((long)RE_SACC_FILTER_LEN)*((long)srate)/1000L);
  if((heog_filter_len-heog_deriv_len+1 < blink_filter_len-blink_deriv_len+1)
   ||(blink_deriv_len+blink_filter_len < heog_filter_len-heog_deriv_len+1)
   ||(heog_deriv_len+heog_filter_len < blink_deriv_len+blink_filter_len))
    {
    fclose(cnt);
    fprintf(stderr, "Incompatible window lengths %d,%d,%d,%d at %dHz\n", blink_deriv_len, blink_filter_len, heog_deriv_len, heog_filter_len, srate);
    exit(1);
    }
  num_events = T_Size(evt_tbl_hdr)/(long)packed_sizeof_EVENT1;
  response_window_start = ((long)(0.5+(srate*((double)ST_FASTEST_RESPONSE+ST_LATENCY))/1000.0))*(num_channels*sizeof(short));
  response_window_end = ((long)(0.5+(srate*((double)ST_SLOWEST_RESPONSE+ST_LATENCY))/1000.0))*(num_channels*sizeof(short));

/*use the HEOG calibrations to find resting HEOG & rejection threshold:*/
  num_heog_calibrations = 0;
/*count the calibrations*/
  for(event = 0L; (event != num_events) && (num_heog_calibrations != RE_NUM_HEOG_CALIBS) && (fread(evtbuf, packed_sizeof_EVENT1, 1, cnt)); event++)
    {
    if(E_StimType(evtbuf) == EV_HEOG_CALIB)
      heog_calibration_offset[num_heog_calibrations++] = E_Offset(evtbuf);
    }
  if((num_heog_calibrations == RE_NUM_HEOG_CALIBS) && (event != num_events))
    fprintf(stderr, "Warning: %s: more than %d calibrations in '%s' - excess not processed\n", *argv, RE_NUM_HEOG_CALIBS, argv[1]);
/*average the interval [-25ms,+25ms] around each calibration mark*/
  for(cal = 0; cal != num_heog_calibrations; cal++)
    {
    temp = 0L;
    fseek(cnt, heog_calibration_offset[cal]+(heog_channel-(srate/40)*num_channels)*sizeof(short), SEEK_SET);
    for(sample = 0; sample != srate/20; sample++)
      {
      temp += read_Intel_short(cnt);
      fseek(cnt, (long)((num_channels-1)*sizeof(short)), SEEK_CUR);
      }
    heog_calibration_value[cal] = (short)(temp/((long)(srate/20)));
    }
/*examine each pair of calibrations*/
  num_good_calibrations = 0;
  for(cal = 0; cal != num_heog_calibrations-1; cal++)
    {
  /*verify that the pair comes reasonably close together*/
    if((heog_calibration_offset[cal+1]-heog_calibration_offset[cal])/(num_channels*sizeof(short)) <= ((long)RE_MAX_CALIB_INTERVAL)*((long)srate))
      {
      heog_calibration_value[num_good_calibrations] = heog_calibration_value[cal]-heog_calibration_value[cal+1];
    /*verify that the leftward-gazing and rightward-gazing marks are far apart*/
      if(heog_calibration_value[num_good_calibrations]*heog_scale >= RE_MIN_CALIB_LEVEL*ST_ANGLE)
	num_good_calibrations++;
      }
    }
  if(num_good_calibrations == 0)
    {
    fprintf(stderr, "%s: warning: no good HEOG calibrations - using default level\n", argv[1]);
    heog_threshold = (long)(0.5+RE_DFL_CALIB_LEVEL*heog_threshold_angle/heog_scale);
    }
  else
    heog_threshold = (long)(0.5+(heog_threshold_angle/ST_ANGLE)*median(heog_calibration_value, num_good_calibrations));
  printf("Saccade threshold is %ld (%.2fuV)\n", heog_threshold, heog_scale*heog_threshold);

/*set up a Type 2 event table in a temporary file*/
  T_Teeg(evt_tbl_hdr) = TEEG_EVENT_TAB2;
  put_T_Size(evt_tbl_hdr, num_events*(long)packed_sizeof_EVENT2);
  if((evt = fopen(EVENT_FILE,
#ifdef MSDOS
    "wb+"
#else
    "w+"
#endif
    )) == NULL)
    {
    fclose(cnt);
    fprintf(stderr, "Couldn't open temporary file.\n");
    exit(errno);
    }
  event = T_Offset(evt_tbl_hdr);
  put_T_Offset(evt_tbl_hdr, 0L);
  if(fwrite(evt_tbl_hdr, packed_sizeof_TEEG, 1, evt) != 1)
    {
    fclose(cnt);
    fprintf(stderr, "Couldn't write to temporary file.\n");
    exit(errno);
    }
    
/*return to the first event*/
  fseek(cnt, S_EventTablePos(erp)+packed_sizeof_TEEG+event, SEEK_SET);
  num_blink_rejects = num_saccade_rejects = 0L;
  num_hits = num_misses = num_false_alarms = num_glitches = num_dropped_responses = 0;
  summed_latencies = summed_squared_latencies = 0.0;
  offset_of_last_scored_response = -1L;
/*process events*/
  event = 0L;
  offset_of_last_scored_response = 0L;
/*open a duplicate file pointer to obviate seeking between EEG and event table*/
  cnt2 = fopen(argv[1],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    );
  while(event != num_events)
  /*inv: 'event' events have been processed*/
    {
    E_EpochEvent(evtbuf2) = 0;
    E_Accuracy(evtbuf2) = 0;
    put_E_Latency(evtbuf2, 0.0);
    E_Accept(evtbuf2) = 1;
    put_E_Code(evtbuf2, 0);
    do
    /*inv: 'event' events have been copied to the temporary file, and 'stimtype'
      is the event code from the last event copied.*/
      {
      fread(evtbuf, packed_sizeof_EVENT1, 1, cnt);
      memcpy(evtbuf2, evtbuf, packed_sizeof_EVENT1);
      stimtype = E_StimType(evtbuf);
      put_E_Type(evtbuf2, stimtype);
      fwrite(evtbuf2, packed_sizeof_EVENT2, 1, evt);
      event++;
      } while((stimtype != EV_START) && (stimtype != EV_START+1) && (event != num_events));
    if(event != num_events)
      {
    /*this is the beginning of a block of stimuli*/
      attended_target = 2*(stimtype-EV_START+1);
    /*get its first stimulus event; transcribe any intervening joystick events*/
      fread(evtbuf, packed_sizeof_EVENT1, 1, cnt);
      event++;
      stimtype = E_StimType(evtbuf);
      while(((stimtype < 1) || (stimtype > 4)) && (stimtype != EV_FINISH) && (event != num_events))
	{
	memcpy(evtbuf2, evtbuf, packed_sizeof_EVENT1);
	put_E_Type(evtbuf2, stimtype);
	E_EpochEvent(evtbuf2) = 0;
	E_Accuracy(evtbuf2) = 1;
	put_E_Latency(evtbuf2, 0.0);
	E_Accept(evtbuf2) = 1;
	put_E_Code(evtbuf2, 0);
	fwrite(evtbuf2, packed_sizeof_EVENT2, 1, evt);
	fread(evtbuf, packed_sizeof_EVENT1, 1, cnt);
	event++;
	stimtype = E_StimType(evtbuf);
	}
    /*set up offset to an imaginary preceding stimulus*/
      prev_offset = E_Offset(evtbuf)-(((long)ST_MSECS_PER_TICK)*srate/1000L/2L)*num_channels*sizeof(short);
    /*insert boundary points into HEOG filter*/
      fseek(cnt2, E_Offset(evtbuf)-(((long)ST_MSECS_PER_TICK)*srate/1000L/4L+heog_deriv_len/2L+heog_filter_len/2L+1L)*num_channels*sizeof(short)+heog_channel*sizeof(short), SEEK_SET);
      heog_filter = create_btree(heog_filter_len);
      for(sample = -heog_deriv_len/2-heog_filter_len/2+1; sample != -blink_deriv_len/2-blink_filter_len/2+1; sample++)
	{
	insert_newest(read_Intel_short(cnt2), heog_filter);
	fseek(cnt2, (num_channels-1)*sizeof(short), SEEK_CUR);
	sample++;
	}
    /*HEOG filter size is HD/2+HF/2-BD/2-BF/2*/
    /*insert boundary points into HEOG & blink filters*/
      blink_filter = create_btree(blink_filter_len);
      while(sample != -blink_deriv_len/2+(blink_filter_len+1)/2)
	{
	insert_newest(read_Intel_short(cnt2), heog_filter);
	fseek(cnt2, (long)((blink_channel-heog_channel-1)*sizeof(short)), SEEK_CUR);
	insert_newest(read_Intel_short(cnt2), blink_filter);
	fseek(cnt2, (long)((num_channels+heog_channel-blink_channel-1)*sizeof(short)), SEEK_CUR);
	sample++;
	}
    /*HEOG filter size is HD/2+HF/2-BD/2+(BF+1)/2-1, blink filter size is BF-1*/
    /*insert into HEOG filter, and window the output of blink filter*/
      blink_window = create_btree(blink_deriv_len);
      while(sample != -heog_deriv_len/2+(heog_filter_len+1)/2)
	{
	insert_newest(read_Intel_short(cnt2), heog_filter);
	fseek(cnt2, (long)((blink_channel-heog_channel-1)*sizeof(short)), SEEK_CUR);
	insert_newest(read_Intel_short(cnt2), blink_filter);
	fseek(cnt2, (long)((num_channels+heog_channel-blink_channel-1)*sizeof(short)), SEEK_CUR);
	insert_newest(extract_median(blink_filter), blink_window);
	delete_oldest(blink_filter);
	sample++;
	}
    /*HEOG filter size is HF-1, blink filter size is BF-1,
      blink window size is -HD/2+(HF+1)/2+BD/2-(BF+1)/2*/
    /*window the outputs of blink and HEOG filters*/
      heog_window = create_queue(heog_deriv_len-1);
      while(sample != (blink_deriv_len+1)/2+(blink_filter_len+1)/2-1)
	{
	insert_newest(read_Intel_short(cnt2), heog_filter);
	fseek(cnt2, (long)((blink_channel-heog_channel-1)*sizeof(short)), SEEK_CUR);
	insert_newest(read_Intel_short(cnt2), blink_filter);
	fseek(cnt2, (long)((num_channels+heog_channel-blink_channel-1)*sizeof(short)), SEEK_CUR);
	insert_newest(extract_median(blink_filter), blink_window);
	delete_oldest(blink_filter);
	queue_insert(extract_median(heog_filter), heog_window);
	delete_oldest(heog_filter);
	sample++;
	}
    /*HEOG filter size is HF-1, blink filter size is BF-1, blink window size is
      BD-1, HEOG window size is (BD+1)/2+(BF+1)/2+HD/2-(HF+1)/2-1*/
    /*window output of the HEOG filter, and queue output of the blink filter*/
      blinkq = create_queue((heog_deriv_len+1)/2-(blink_deriv_len+1)/2-(blink_filter_len+1)/2+(heog_filter_len+1)/2);
      while(sample != (heog_deriv_len+1)/2+(heog_filter_len+1)/2-1)
	{
	insert_newest(read_Intel_short(cnt2), heog_filter);
	fseek(cnt2, (long)((blink_channel-heog_channel-1)*sizeof(short)), SEEK_CUR);
	insert_newest(read_Intel_short(cnt2), blink_filter);
	fseek(cnt2, (long)((num_channels+heog_channel-blink_channel-1)*sizeof(short)), SEEK_CUR);
	queue_insert(extract_median(blink_filter), blinkq);
	delete_oldest(blink_filter);
	queue_insert(extract_median(heog_filter), heog_window);
	delete_oldest(heog_filter);
	sample++;
	}
      stimtype = E_StimType(evtbuf);
      while((stimtype != EV_FINISH) && (event != num_events))
      /*inv: 'evtbuf' contains the current stimulus event, 'stimtype' contains
	its type, 'prev_offset' contains offset of the previous stimulus event,
	HEOG filter size is HF-1, blink filter size is BF-1,
	HEOG window size is HD-1, blink window size is BD-1,
	blink queue size is (HD+1)/2+(HF+1)/2-(BD+1)/2-(BF+1)/2*/
	{
        long savedpos;
	memcpy(evtbuf2, evtbuf, packed_sizeof_EVENT1);
	put_E_Type(evtbuf2, stimtype);
	E_EpochEvent(evtbuf2) = (stimtype == attended_target);
	E_Accuracy(evtbuf2) = 0;
	put_E_Latency(evtbuf2, 0.0);
	E_Accept(evtbuf2) = 1;
	put_E_Code(evtbuf2, 0);

      /*compute accuracy and latency*/
	savedpos = ftell(cnt);
	if(E_EpochEvent(evtbuf2))
	  {
	/*skip forward to beginning of response window*/
	  do
	    fread(evtbuf, packed_sizeof_EVENT1, 1, cnt);
	  while(E_Offset(evtbuf) - E_Offset(evtbuf2) < response_window_start);
	/*scan forward until end-of-trial, end of response window, or response*/
	  while((E_StimType(evtbuf) != EV_FINISH)
	   && (E_Offset(evtbuf) - E_Offset(evtbuf2) <= response_window_end)
	   && (E_StimType(evtbuf) != EV_RESPONSE+(4-E_StimType(evtbuf2))/2)
	   && !((E_Offset(evtbuf) > offset_of_last_scored_response) && (E_StimType(evtbuf) == EV_RESPONSE-1+E_StimType(evtbuf2)/2)))
	    fread(evtbuf, packed_sizeof_EVENT1, 1, cnt);
	/*Sometimes the flaky event recorder confuses left and right.
	  Fortunately, the data are redundant enough to allow the identification
	  and reconstruction of cases in which this occurs.  If we do find an
	  event within the response window that hasn't already been scored, we
	  assume that it's a response to the most recent target, even if its
	  laterality is incorrect.  We increment 'num_glitches' a count of such
	  events, and then rewrite the event code so that it has the correct
	  laterality.*/
	  if((E_Offset(evtbuf) - E_Offset(evtbuf2) <= response_window_end)
	   && (E_Offset(evtbuf) > offset_of_last_scored_response)
	   && (E_StimType(evtbuf) == EV_RESPONSE-1+E_StimType(evtbuf2)/2))
	    {
	    num_glitches++;
	    put_E_StimType(evtbuf, EV_RESPONSE+(4-E_StimType(evtbuf2))/2);
	    }
        /*fill Accuracy and Latency fields from response data*/
	  E_Accuracy(evtbuf2) = (E_Offset(evtbuf) - E_Offset(evtbuf2) <= response_window_end) && (E_StimType(evtbuf) == EV_RESPONSE+(4-E_StimType(evtbuf2))/2);
	  if(E_Accuracy(evtbuf2))
	    {
	  /*this trial was a hit; record latency and toggle the attended side*/
	    put_E_Latency(evtbuf2, (E_Offset(evtbuf) - E_Offset(evtbuf2))/((float)(sizeof(short)*num_channels*srate))-ST_LATENCY/1000.0);
	    attended_target = 2+attended_target%4;
	  /*Half a tick is added to 'offset_of_last_scored_response', in case
	    debouncing on the stimulus delivery machine didn't work*/
	    offset_of_last_scored_response = E_Offset(evtbuf)+(srate*num_channels*(long)(sizeof(short)*(ST_MSECS_PER_TICK/2)))/1000L;
	    }
	  else
	    {
	  /*this trial seems a miss - but check for a dropped response code*/
	    long offset_of_next_target;
	  /*Search for any responses up until the next target on this same side,
	    or until the end of the current trial, whichever comes first.*/
	    while((E_StimType(evtbuf) != EV_RESPONSE)
	     && (E_StimType(evtbuf) != EV_RESPONSE+1)
	     && (E_StimType(evtbuf) != E_StimType(evtbuf2))
	     && (E_StimType(evtbuf) != EV_FINISH))
	      fread(evtbuf, packed_sizeof_EVENT1, 1, cnt);
	  /*If the search above was terminated by the occurrence of the next
	    target, continue the search until the response window for that
	    target begins.*/
	    offset_of_next_target = E_Offset(evtbuf);
	    while((E_StimType(evtbuf) != EV_RESPONSE)
	     && (E_StimType(evtbuf) != EV_RESPONSE+1)
	     && (E_Offset(evtbuf) - offset_of_next_target < response_window_start)
	     && (E_StimType(evtbuf) != EV_FINISH))
	      fread(evtbuf, packed_sizeof_EVENT1, 1, cnt);
	  /*If a response to a target on the opposite side occurred at any time
	    from the beginning of the response window for the current target up
	    until the beginning of the response window for the next target on
	    the same side, with no intervening response to the current target,
	    then it must be the case that a response to the current target did
	    occur but was not recorded.  Otherwise the joystick couldn't have
	    got over to the other side.  So set the Accuracy flag for this
	    response, but don't process any latency information, because the
	    exact time of the response is unknown.*/
	    if((E_Offset(evtbuf) - offset_of_next_target < response_window_start)
	     && (E_StimType(evtbuf) == EV_RESPONSE-1+attended_target/2))
	      {
	      E_Accuracy(evtbuf2) = 1;
	      attended_target = 2+attended_target%4;
	      }
	    put_E_Latency(evtbuf2, 0.0);
	    }
	  }

      /*test for eye artefact*/
	for(sample = 0; sample != (E_Offset(evtbuf)-prev_offset)/(num_channels*sizeof(short)); sample++)
	/*inv: 'sample' EOG points from the interval surrounding the current
	  stimulus have been tested for ocular artefact*/
	  {
	  short temp;
	  insert_newest(read_Intel_short(cnt2), heog_filter);
	  fseek(cnt2, (long)((blink_channel-heog_channel-1)*sizeof(short)), SEEK_CUR);
	  insert_newest(read_Intel_short(cnt2), blink_filter);
	  fseek(cnt2, (long)((num_channels+heog_channel-blink_channel-1)*sizeof(short)), SEEK_CUR);
	  temp = extract_median(heog_filter);
	  delete_oldest(heog_filter);
	  if(E_Accept(evtbuf2) && (abs(temp - queue_extract(heog_window)) > heog_threshold))
	    {
	    E_Accept(evtbuf2) = 0;
	    put_E_Code(evtbuf2, CC_SACCADE);
	    }
	  queue_insert(temp, heog_window);
	  insert_newest(queue_extract(blinkq), blink_window);
	  if(extract_maximum(blink_window)-extract_minimum(blink_window) > blink_threshold)
	    {
	    E_Accept(evtbuf2) = 0;
	    put_E_Code(evtbuf2, CC_BLINK);
	    }
	  delete_oldest(blink_window);
	  queue_insert(extract_median(blink_filter), blinkq);
	  delete_oldest(blink_filter);
	  }
      /*write the current event record*/
	fwrite(evtbuf2, packed_sizeof_EVENT2, 1, evt);
	prev_offset = E_Offset(evtbuf2);
      /*update statistics*/
	if(!E_Accept(evtbuf2))
	  {
	  if(E_Code(evtbuf2) == CC_BLINK)
	    num_blink_rejects++;
	  else
	    num_saccade_rejects++;
	  }
	if(E_EpochEvent(evtbuf2) && E_Accept(evtbuf2))
	  {
	  if(E_Accuracy(evtbuf2))
	    {
	    if(E_Latency(evtbuf2) == 0.0)
	      num_dropped_responses++;
	    else
	      {
	      num_hits++;
	      summed_latencies += E_Latency(evtbuf2);
	      summed_squared_latencies += E_Latency(evtbuf2) * E_Latency(evtbuf2);
	      }
	    }
	  else
	    num_misses++;
	  }
	fseek(cnt, savedpos, SEEK_SET);
      /*transcribe any intervening joystick events, & get next stimulus event*/
	fread(evtbuf, packed_sizeof_EVENT1, 1, cnt);
	event++;
	stimtype = E_StimType(evtbuf);
	while(((stimtype < 1) || (stimtype > 4)) && (stimtype != EV_FINISH) && (event != num_events))
	  {
	  memcpy(evtbuf2, evtbuf, packed_sizeof_EVENT1);
	  put_E_Type(evtbuf2, stimtype);
	  E_EpochEvent(evtbuf2) = 0;
	  E_Accuracy(evtbuf2) = 1;
	  put_E_Latency(evtbuf2, 0.0);
	  E_Accept(evtbuf2) = 1;
	  put_E_Code(evtbuf2, 0);
	  if((E_Offset(evtbuf) > offset_of_last_scored_response) && ((stimtype == EV_RESPONSE) || (stimtype == EV_RESPONSE+1)))
	    {
	  /*This response has not been associated with any prior stimulus.
	    Flag it s.t. any and all unattended targets in whose response
	    window it falls can be marked during a subsequent pass.  Increment
	    the false alarm count, and toggle the attended side.*/
	    E_EpochEvent(evtbuf2) = 2;
	    num_false_alarms++;
	    attended_target = 2*(E_StimType(evtbuf)-EV_RESPONSE+1);
	    }
	  fwrite(evtbuf2, packed_sizeof_EVENT2, 1, evt);
	  fread(evtbuf, packed_sizeof_EVENT1, 1, cnt);
	  event++;
	  stimtype = E_StimType(evtbuf);
	  } /*end of stimulus*/
        } /*end of block*/
      destroy_queue(blinkq);
      destroy_queue(heog_window);
      destroy_btree(blink_window);
      destroy_btree(blink_filter);
      destroy_btree(heog_filter);
    /*write the end-of-block marker event*/
      memcpy(evtbuf2, evtbuf, packed_sizeof_EVENT1);
      put_E_Type(evtbuf2, EV_FINISH);
      E_EpochEvent(evtbuf2) = 0;
      E_Accuracy(evtbuf2) = 1;
      put_E_Latency(evtbuf2, 0.0);
      E_Accept(evtbuf2) = 1;
      put_E_Code(evtbuf2, 0);
      fwrite(evtbuf2, packed_sizeof_EVENT2, 1, evt);
      }
    }
  fclose(cnt2);

/*Scan backward through the newly created event table, looking for responses
  that have not been associated with stimuli.  For each of these 'false alarm'
  responses, flag all the ignored targets whose response windows include the
  time of occurrence of the false alarm.*/
  window_start = window_end = next_window_end = 0L;
  for(event = 0L; event != num_events; event++)
    {
  /*inv: 'event' events have been processed, 'evtbuf2' contains the next event
    to be processed, 'evt' points just past this event in the event file, and
    one of the following four cases is true:
     (0) Ignored targets are not currently being flagged, no onset of a response
      window is pending, window_end=0, and window_start=0.
     (1) Ignored targets are not currently being flagged, window_start is the
      event table offset of the pending onset of a response window, window_end
      is the event table offset of that window or of the most recently
      encountered window that overlaps it, and next_window_end=0.
     (2) Ignored targets are being flagged, window_end is the event table offset
      beyond which flagging will cease, no onset of a subsequent,
      non-overlapping response window is pending, and window_start=0.
     (3) Ignored targets are being flagged, window_end is the event table offset
      beyond which flagging will cease, window_start is the event table offset
      of the pending onset of a non-overlapping response window, and
      next_window_end is the offset beyond which flagging of that pending
      window will cease.*/
    if(next_window_end != 0L)
    /*condition 3: we're in the middle of one window and have another pending*/
      {
      if(E_EpochEvent(evtbuf2) == 2)
      /*if a flagged response is seen, extend the range of the pending window*/
	{
	next_window_end = E_Offset(evtbuf2)-response_window_end;
	E_EpochEvent(evtbuf2) = 0;
	}
      else if((E_EpochEvent(evtbuf2) == 0) && ((E_StimType(evtbuf2) == 2) || (E_StimType(evtbuf2) == 4)))
      /*if a target in the ignored channel is seen, flag it as the possible
	source of the false alarm*/
	E_EpochEvent(evtbuf2) = 2;
      if(E_Offset(evtbuf2) <= window_end)
      /*if we're at the end of the current window, go to condition (1)*/
	{
	window_end = next_window_end;
	next_window_end = 0L;
	}
      }
    else if(window_start != 0L)
    /*condition 1: a window is pending but has not yet opened*/
      {
      if(E_EpochEvent(evtbuf2) == 2)
      /*if a flagged response is seen, extend the range of the pending window*/
	{
	window_end = E_Offset(evtbuf2)-response_window_end;
	E_EpochEvent(evtbuf2) = 0;
	}
      if(E_Offset(evtbuf2) <= window_start)
      /*if we're at the beginning of the pending window, go to condition (2)*/
	window_start = 0L;
      }
    else if(window_end != 0L)
    /*condition 2: we're in the middle of a window with no window pending*/
      {
      if(E_EpochEvent(evtbuf2) == 2)
	{
	if(E_Offset(evtbuf2)-window_end < response_window_start)
	/*if a flagged response is seen and we're within 200ms of the end of the
	  current window, create a new, pending window by going to condition
	  (3)*/
	  {
	  window_start = E_Offset(evtbuf2)-response_window_start;
	  next_window_end = E_Offset(evtbuf2)-response_window_end;
	  }
	else
	/*if a flagged response is seen and we're not within 200ms of the end of
	  the current window, extend the range of the current window and remain
	  in condition (2)*/
	  window_end = E_Offset(evtbuf2)-response_window_end;
	E_EpochEvent(evtbuf2) = 0;
	}
      else if((E_EpochEvent(evtbuf2) == 0) && ((E_StimType(evtbuf2) == 2) || (E_StimType(evtbuf2) == 4)))
      /*if a target in the ignored channel is seen, flag it as the possible
	source of the false alarm*/
	E_EpochEvent(evtbuf2) = 2;
      if(E_Offset(evtbuf2) <= window_end)
      /*if we've reached the end of the current window, go to condition (0)*/
	window_end = 0L;
      }
    else
    /*condition 0: no window is current and no window is pending*/
      if(E_EpochEvent(evtbuf2) == 2)
      /*if a flagged response is seen, go to condition (1)*/
	{
	window_start = E_Offset(evtbuf2)-response_window_start;
	window_end = E_Offset(evtbuf2)-response_window_end;
	E_EpochEvent(evtbuf2) = 0;
	}
    fseek(evt, -(long)packed_sizeof_EVENT2, SEEK_CUR);
    fwrite(evtbuf2, packed_sizeof_EVENT2, 1, evt);
    fseek(evt, -2L*packed_sizeof_EVENT2, SEEK_CUR);
    fread(evtbuf2, packed_sizeof_EVENT2, 1, evt);
    }
  rewind(evt);

/*overwrite the old event table with the New & Improved event table*/
  fseek(cnt, S_EventTablePos(erp), SEEK_SET);
  while((sample = getc(evt)) != EOF)
    putc(sample, cnt);
  fclose(evt);
  fclose(cnt);
  unlink(EVENT_FILE);
  printf("%ld events: %ld blink rejects (%ld%%), %ld saccade rejects (%ld%%).\n%d hits, %d misses (%d%%), %d false alarms\nresponse time avg = %f SD = %f\n",
    num_events, num_blink_rejects, 100L*num_blink_rejects/num_events, num_saccade_rejects, 100L*num_saccade_rejects/num_events,
    num_hits+num_dropped_responses, num_misses, (int)(100L*(num_hits+num_dropped_responses)/(num_hits+num_dropped_responses+num_misses)), num_false_alarms,
    1000.0*summed_latencies/num_hits, 1000.0*sqrt((summed_squared_latencies - summed_latencies*summed_latencies/num_hits)/num_hits));
  if(num_glitches)
    fprintf(stderr, "Warning: %d responses were incorrectly lateralised by the event recorder.\n", num_glitches);
  if(num_dropped_responses)
    fprintf(stderr, "Warning: %d responses were dropped by the event recorder.  These were assumed\nto be hits rather than spurious responses, and were included in the totals but\nnot in the computation of average response time.\n", num_dropped_responses);
  exit(0);
  }
