/*eyehand.c - correlate saccades with button responses in a continuous EEG data
  file set up with fastshft.c and scored with score.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 <setjmp.h>
#include "fastshft.h"
#include "median.h"
#include "fisher.h"
#include "intel.h"
#include "gnuro.h"

#define MAX_NUM_POINTS 1024

extern void mrqmin(float x[], float y[], float sig[], int ndata, float a[],
	int ia[], int ma, float**, float**, float*,
	void (*funcs)(double, float[], float*, float[], int), float*);
jmp_buf env;

/*y(x) = p1*exp(-x**2/p2**2)+p3*exp(-x**2/p4**2) is the sum of two Gaussians
  whose amplitudes are the odd elements of 'param' and whose widths are the
  even elements of 'param'.  Each element of 'deriv' is the partial derivative
  of the corresponding element of 'param' w.r.t. x.*/
void saccade_funcs(x, param, y, deriv, four)
double x;
/*must be double instead of float to prevent stack allignment problems due to
  default argument promotion - corresponding changes must be made to the
  Numerical Recipes files mrqmin.c and mrqcof.c.*/
float *param, *y, *deriv;
int four;
  {
  float chain_factor, scaled_x;
  scaled_x = x/param[2];
  deriv[1] = exp(-scaled_x*scaled_x);
  chain_factor = 2.0*param[1]*deriv[1]*scaled_x;
  *y = param[1]*deriv[1];
  deriv[2] = chain_factor*scaled_x/param[2];
  scaled_x = x/param[4];
  deriv[3] = exp(-scaled_x*scaled_x);
  chain_factor = param[3]*deriv[3]*2.0*scaled_x;
  *y += param[3]*deriv[3];
  deriv[4] = chain_factor*scaled_x/param[4];
  }

/*Find the t-value associated with two-tailed probability 'p' for 'df' degrees
  of freedom.*/
double tcrit(df, p)
int df;
double p;
  {
  register double t, ub, lb, cur_p;
  ub = 1.0;
  lb = 0.0;
  while(1.0-student(df, ub) > p)
    {
    lb = ub;
    ub *= 2.0;
    }
  t = (lb+ub)/2.0;
  while((t != lb) && (t != ub))
    {
    cur_p = 1.0-student(df, t);
    if(cur_p < p)
      ub = t;
    else
      lb = t;
    t = (lb+ub)/2.0;
    }
  return(t);
  }

void nrerror(mesg)
char *mesg;
  {
  fprintf(stderr, "%s\n", mesg);
  longjmp(env, 1);
  }

void main(argc, argv)
int argc;
char *argv[];
  {
/*counters*/
  register int sample, channel, num_points, num_iterations, num_latencies;
  register long event;
/*registers*/
  register int stimtype, heog, prev_heog, deriv, v0, step_interval;
  register int last_nonzero_sample;
  register long start_of_block;
  float chi_squared, lambda;
  register float prev_lambda;
  register int heog_threshold;
  register double saccade_latency, button_latency;
/*constants*/
  register long response_window_start, response_window_end;
  register int num_channels, srate;
  register double heog_scale;
  int heog_channel;
  long evtpos, num_events;
  static int var_flags[] = {1, 1, 1, 1};
/*files*/
  FILE *cnt, *evt;
/*structure sizes*/
  int heog_filter_len;
/*statistics*/
  register double button_sum, button_ss, saccade_sum, saccade_ss, product_sum;
  double button_variance_times_n, saccade_variance_times_n, covariance_times_n, coeff_of_determination, f_ratio, slope, intercept;
/*buffers*/
  float parameters[4], *(covar[4]), *(alpha[4]), covar_storage[16], alpha_storage[16], x[MAX_NUM_POINTS], y[MAX_NUM_POINTS], dummy_sd[MAX_NUM_POINTS];
  btree *heog_filter;
  char	erp[packed_sizeof_SETUP],
	channel_hdr[packed_sizeof_ELECTLOC],
	evt_tbl_hdr[packed_sizeof_TEEG],
	evtbuf2[packed_sizeof_EVENT2],
	prev_event[packed_sizeof_EVENT2];
  unsigned short *histogram, hist[0x1fffd];
  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
/*set up data areas for density estimation*/
  histogram = hist+0xfffe;
  for(sample = -0xfffe; sample != 0xffff; sample++)
    histogram[sample] = 0;
  for(sample = 0; sample != 4; sample++)
    {
    covar[sample] = covar_storage+sample*4-1;
    alpha[sample] = alpha_storage+sample*4-1;
    }
/*check parameters*/
  if(argc < 2)
    {
    fprintf(stderr, "Usage: %s <scored_file.cnt>\n", *argv);
    exit(1);
    }
/*open data file*/
  if((evt = fopen(argv[1],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[1]);
    exit(errno);
    }
/*read main header*/
  if(fread(erp, packed_sizeof_SETUP, 1, evt) != 1)
    {
    perror(argv[1]);
    fclose(evt);
    exit(errno);
    }
  num_channels = S_nchannels(erp);
  srate = S_rate(erp);
/*determine channel assignments*/
  heog_channel = -1;
  for(channel = 0; channel != num_channels; channel++)
    {
    fread(channel_hdr, packed_sizeof_ELECTLOC, 1, evt);
    if(strcasecmp("heog", channel_hdr) == 0)
      {
      heog_channel = channel;
      heog_scale = ((double)(EL_sensitivity(channel_hdr)))*((double)(EL_calib(channel_hdr)))/204.8;
      }
    }
  if(heog_channel == -1)
    {
    fprintf(stderr, "Couldn't find HEOG channel.\n");
    fclose(evt);
    exit(1);
    }
/*read event table header*/
  fseek(evt, S_EventTablePos(erp), SEEK_SET);
  if(fread(evt_tbl_hdr, packed_sizeof_TEEG, 1, evt) != 1)
    {
    perror(argv[1]);
    fclose(evt);
    exit(errno);
    }
/*verify that this is the scored version of the event table structure*/
  if(T_Teeg(evt_tbl_hdr) != TEEG_EVENT_TAB2)
    {
    fprintf(stderr, "You must first score %s.\n", argv[1]);
    fclose(evt);
    exit(1);
    }
/*compute window lengths*/
  heog_filter_len = (int)(((long)RE_SACC_FILTER_LEN)*((long)srate)/1000L);
  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));
  num_events = T_Size(evt_tbl_hdr)/(long)packed_sizeof_EVENT2;
/*position to first event*/
  fseek(evt, T_Offset(evt_tbl_hdr), SEEK_CUR);
  evtpos = ftell(evt);
/*process events*/
  event = 0L;
/*open a duplicate file pointer to obviate seeking between EEG and event table*/
  cnt = fopen(argv[1],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    );

/*compute the density function of the sample of saccades made within blocks*/
  while(event != num_events)
  /*inv: 'event' events have been processed*/
    {
  /*find the beginning of a block of trials*/
    do
    /*inv: 'event' events have been processed, and 'stimtype'
      is the event code from the last event processed.*/
      {
      fread(evtbuf2, packed_sizeof_EVENT2, 1, evt);
      stimtype = E_StimType(evtbuf2);
      event++;
      } while((stimtype != EV_START) && (stimtype != EV_START+1) && (event != num_events));
    if(event != num_events)
      {
    /*get the block's first target event*/
      do
	{
	fread(evtbuf2, packed_sizeof_EVENT2, 1, evt);
	event++;
	stimtype = E_StimType(evtbuf2);
	}
      while(((stimtype != 2) && (stimtype != 4)) && (stimtype != EV_FINISH) && (event != num_events));
      if(((stimtype == 2) || (stimtype == 4)) && (event != num_events))
	{
      /*remember offset of start of block*/
	start_of_block = E_Offset(evtbuf2);
      /*find end of block*/
	do
	  {
	  fread(evtbuf2, packed_sizeof_EVENT2, 1, evt);
	  event++;
	  stimtype = E_StimType(evtbuf2);
	  }
	while((stimtype != EV_FINISH) && (event != num_events));
      /*insert boundary points into HEOG filter*/
	fseek(cnt, start_of_block+(heog_channel-(heog_filter_len/2L)*num_channels)*sizeof(short), SEEK_SET);
	heog_filter = create_btree(heog_filter_len);
	for(sample = 0; sample != heog_filter_len; sample++)
	  {
	  insert_newest(read_Intel_short(cnt), heog_filter);
	  fseek(cnt, (num_channels-1)*sizeof(short), SEEK_CUR);
	  sample++;
	  }
	heog = extract_median(heog_filter);
	sample = (E_Offset(evtbuf2)-start_of_block)/(num_channels*sizeof(short));
	do
	/*inv: 'heog_filter' contains the last 'heog_filter_len' samples,
	  'heog' is the median of these samples, and 'cnt' is positioned at the
	  next sample.*/
	  {
	  prev_heog = heog;
	  delete_oldest(heog_filter);
	  insert_newest(read_Intel_short(cnt), heog_filter);
	  heog = extract_median(heog_filter);
	  fseek(cnt, (num_channels-1)*sizeof(short), SEEK_CUR);
	  sample--;
	  }
	while((prev_heog == heog) && (sample != 1));
      /*the instant between 'prev_heog' and 'heog' is the first increase or
	decrease during this block - perhaps the tail end of a variation that
	began before the start of the block*/
	deriv = heog-prev_heog;
	do
	/*inv: 'heog_filter' contains the last 'heog_filter_len' samples,
	  'heog' is the median of these samples, and 'cnt' is positioned at the
	  next sample.*/
	  {
	  prev_heog = heog;
	  delete_oldest(heog_filter);
	  insert_newest(read_Intel_short(cnt), heog_filter);
	  heog = extract_median(heog_filter);
	  fseek(cnt, (num_channels-1)*sizeof(short), SEEK_CUR);
	  sample--;
	  }
	while((deriv*(heog-prev_heog) >= 0) && (sample != 0));
      /*the instant between 'prev_heog' and 'heog' is the beginning of the first
	decrease or increase whose full extent is contained within this block*/
	while(sample != 0)
	  {
	/*inv: all continuously nonincreasing or nondecreasing intervals up to
	  the current sample have been included in the histogram, 'prev_heog' is
	  the point at the end of the previous interval and at the beginning of
	  the current interval, and 'heog' is the point immediately following
	  'prev_heog'.*/
	  deriv = heog-prev_heog;
	  v0 = prev_heog;
	  do
	  /*inv: 'heog_filter' contains the last 'heog_filter_len' samples,
	    'heog' is the median of these samples, 'cnt' is positioned at the
	    next sample, and 'v0' is the HEOG value at the beginning of the
	    current nondecreasing or nonincreasing interval, according to
	    whether 'deriv' is positive or negative, respectively.*/
	    {
	    prev_heog = heog;
	    delete_oldest(heog_filter);
	    insert_newest(read_Intel_short(cnt), heog_filter);
	    heog = extract_median(heog_filter);
	    fseek(cnt, (num_channels-1)*sizeof(short), SEEK_CUR);
	    sample--;
	    }
	  while((deriv*(heog-prev_heog) >= 0) && (sample != 0));
	  histogram[prev_heog-v0]++;
	  }
	destroy_btree(heog_filter);
	}
      }
    }

/*eliminate density bins that are empty because of discretisation error*/
  num_points = 0;
  for(step_interval = 1; histogram[step_interval] == 0; step_interval++)
    ;
  last_nonzero_sample = 1;
  for(sample = 1; (sample != 0xffff) && (num_points != MAX_NUM_POINTS/2); sample++)
    {
    if(histogram[sample] != 0)
      {
      if((step_interval != 1) && (sample-last_nonzero_sample == 1))
	{
	x[num_points-1] = (sample*(float)(histogram[sample]) + x[num_points-1]*(float)(histogram[sample-1]))/(float)(histogram[sample-1]+histogram[sample]);
	y[num_points-1] += histogram[sample];
	}
      else
	{
	x[num_points] = sample;
	y[num_points] = histogram[sample];
	dummy_sd[num_points++] = 1.0;
	}
      last_nonzero_sample = sample;
      }
    else if(sample-last_nonzero_sample == 3*step_interval/2)
      {
      x[num_points] = last_nonzero_sample + step_interval;
      y[num_points] = 0.0;
      dummy_sd[num_points++] = 1.0;
      last_nonzero_sample += step_interval;
      }
    }
  last_nonzero_sample = -1;
  for(sample = -1; (sample != -0xffff) && (num_points != MAX_NUM_POINTS); sample--)
    {
    if(histogram[sample] != 0)
      {
      if((step_interval != 1) && (last_nonzero_sample-sample == 1))
	{
	x[num_points-1] = (sample*(float)(histogram[sample]) + x[num_points-1]*(float)(histogram[sample-1]))/(float)(histogram[sample-1]+histogram[sample]);
	y[num_points-1] += histogram[sample];
	}
      else
	{
	x[num_points] = sample;
	y[num_points] = histogram[sample];
	dummy_sd[num_points++] = 1.0;
	}
      last_nonzero_sample = sample;
      }
    else if(last_nonzero_sample-sample == 3*step_interval/2)
      {
      x[num_points] = last_nonzero_sample - step_interval;
      y[num_points] = 0.0;
      dummy_sd[num_points++] = 1.0;
      last_nonzero_sample -= step_interval;
      }
    }

/*initialise the Levenberg-Marquardt minimiser*/
  parameters[0] = (y[0]+y[1]+y[num_points/2]+y[num_points/2+1])/4.0;
  parameters[1] = 1.0/heog_scale;
  parameters[2] = (y[num_points/12]+y[num_points/12+1]+y[num_points/2+num_points/12]+y[num_points/2+num_points/12+1])/4.0;
  parameters[3] = 8.0/heog_scale;
  lambda = -1.0;
  num_iterations = 20;
/*run it to convergence*/
  if(setjmp(env))
  /*'heog_threshold' receives a default value if a singular matrix arose during
    the Levenberg-Marquardt algorithm*/
    heog_threshold = 0.8/heog_scale;
  else
    {
    do
      {
      prev_lambda = lambda;
      mrqmin(x-1, y-1, dummy_sd-1, num_points, parameters-1, var_flags-1, 4, covar-1, alpha-1, &chi_squared, saccade_funcs, &lambda);
      }
    while((--num_iterations > 0) || (lambda < prev_lambda));
  /*clean up*/
    lambda = 0.0;
    mrqmin(x-1, y-1, dummy_sd-1, num_points, parameters-1, var_flags-1, 4, covar-1, alpha-1, &chi_squared, saccade_funcs, &lambda);
/*********FOR DEBUGGING ONLY*********
    for(sample = num_points-1; sample != num_points/2-1; sample--)
      {
      float t0, t1;
      t0 = x[sample]/parameters[1];
      t1 = x[sample]/parameters[3];
      printf("%f\t%f\t%f\n", x[sample]*heog_scale, y[sample], parameters[0]*exp(-t0*t0)+parameters[2]*exp(-t1*t1));
      }
    for(sample = 0; sample != num_points/2; sample++)
      {
      float t0, t1;
      t0 = x[sample]/parameters[1];
      t1 = x[sample]/parameters[3];
      printf("%f\t%f\t%f\n", x[sample]*heog_scale, y[sample], parameters[0]*exp(-t0*t0)+parameters[2]*exp(-t1*t1));
      }
    fprintf(stderr, "Plot file using 1:2 with lines.\n");
************************************/
    printf("HEOG scale: %f\nParameter values:\n", heog_scale);
    for(sample = 0; sample != 4; sample++)
      printf("\t%f\n", parameters[sample]);
    heog_threshold = (int)(0.5+parameters[1]*parameters[3]*sqrt(log(parameters[0]/parameters[2])/(parameters[3]*parameters[3]-parameters[1]*parameters[1])));
   if(heog_threshold < 0.4/heog_scale)
      {
      fprintf(stderr, "Ridiculous parameter values - using default model instead.\n");
      heog_threshold = 0.8/heog_scale;
      }
    }
  printf("Crossover occurs at %fuV.\n", heog_scale*heog_threshold);

/*useful for debugging curve-fitting*/
  if(argc == 3)
    heog_threshold *= atoi(argv[2]);

/*rewind event table and initialise counters and accumulators*/
  fseek(evt, evtpos, SEEK_SET);
  button_sum = button_ss = saccade_sum = saccade_ss = product_sum = 0.0;
  num_latencies = 0;
  event = 0L;
  put_E_StimType(evtbuf2, 0);
  while(event < num_events)
  /*inv: The first 'event' events, up through the current value of 'prev_event',
    have been processed.*/
    {
    memcpy(prev_event, evtbuf2, packed_sizeof_EVENT2);
    stimtype = E_StimType(prev_event);
  /*find the next target that elicited a button response*/
    while(!(((stimtype == 2) || (stimtype == 4)) && E_EpochEvent(prev_event) && E_Accuracy(prev_event) && (E_Latency(prev_event) > 0.0)) && (event != num_events))
    /*inv: There are no accurate responses to targets, with valid latencies,
      between the most recently processed target and 'prev_event'.*/
      {
      fread(prev_event, packed_sizeof_EVENT2, 1, evt);
      stimtype = E_StimType(prev_event);
      event++;
      }
  /*find the (contralateral) target following that target*/
    /*inv: There are no targets between 'prev_event' and 'evtbuf2'.*/
    while(((E_StimType(evtbuf2) != 2+stimtype%4) && (E_StimType(evtbuf2) != EV_FINISH)) && (event != num_events))
      {
      fread(evtbuf2, packed_sizeof_EVENT2, 1, evt);
      event++;
      }
  /*position to beginning of response window for 'prev_event'*/
    fseek(cnt, E_Offset(prev_event)+response_window_start+(heog_channel-(heog_filter_len/2L)*num_channels)*sizeof(short), SEEK_SET);
  /*insert boundary points into HEOG filter*/
    heog_filter = create_btree(heog_filter_len);
    for(sample = 0; sample != heog_filter_len; sample++)
      {
      insert_newest(read_Intel_short(cnt), heog_filter);
      fseek(cnt, (num_channels-1)*sizeof(short), SEEK_CUR);
      sample++;
      }
    saccade_latency = 0.0;
    v0 = heog = extract_median(heog_filter);
    for(sample = 0;
	(sample < (E_Offset(evtbuf2)-E_Offset(prev_event))/(num_channels*sizeof(short)))
	 && (sample < (response_window_end-response_window_start)/(num_channels*sizeof(short)))
	 && (saccade_latency == 0.0);
	sample++)
    /*inv: 'v0' is the HEOG at the beginning of the most recent nonincreasing
      interval if this is a left hit (stimtype==2), or the HEOG at the
      beginning of the most recent nondecreasing interval if this is a right
      hit (stimtype==4).  'heog' is the most recent sample of HEOG.
      'heog_filter' is full.*/
      {
      prev_heog = heog;
      delete_oldest(heog_filter);
      insert_newest(read_Intel_short(cnt), heog_filter);
      heog = extract_median(heog_filter);
      fseek(cnt, (num_channels-1)*sizeof(short), SEEK_CUR);
      if((stimtype-3)*(heog-prev_heog) < 0)
	v0 = heog;
      else if((stimtype-3)*(heog-v0) >= heog_threshold)
	saccade_latency = ST_FASTEST_RESPONSE/1000.0+((double)sample)/srate;
      }
    destroy_btree(heog_filter);
    if(saccade_latency != 0.0)
      {
      button_latency = E_Latency(prev_event);
      button_sum += button_latency;
      button_ss += button_latency*button_latency;
      saccade_sum += saccade_latency;
      saccade_ss += saccade_latency*saccade_latency;
      product_sum += button_latency*saccade_latency;
      num_latencies++;
      }
    }
  fclose(cnt);
  fclose(evt);
  button_variance_times_n = button_ss - button_sum*button_sum/num_latencies;
  saccade_variance_times_n = saccade_ss - saccade_sum*saccade_sum/num_latencies;
  covariance_times_n = product_sum - button_sum*saccade_sum/num_latencies;
  coeff_of_determination = covariance_times_n*covariance_times_n/(button_variance_times_n*saccade_variance_times_n);
  f_ratio = coeff_of_determination*(num_latencies-2) / (1.0-coeff_of_determination);
  slope = covariance_times_n/saccade_variance_times_n;
  intercept = (button_sum-slope*saccade_sum)/num_latencies;
  printf("button mean=%.0fms SD=%.0fms, saccade mean=%.0fms SD=%.0fms\nr=%f, t(%d)=%f, p=%f\nbutton press follows saccade by %.0fms, regression slope is %f+-%f\n",
    1000.0*button_sum/num_latencies,
    1000.0*sqrt(button_variance_times_n/(num_latencies-1)),
    1000.0*saccade_sum/num_latencies,
    1000.0*sqrt(saccade_variance_times_n/(num_latencies-1)),
    sqrt(coeff_of_determination),
    num_latencies-2,
    sqrt(f_ratio),
    1.0-fisher(1, num_latencies-2, f_ratio),
    1000.0*intercept,
    slope,
    tcrit(num_latencies-2, 0.05)
      *sqrt(
	(saccade_ss+button_ss
	  -2.0*slope*product_sum
	  +intercept*(2.0*(slope*saccade_sum-button_sum)+num_latencies*intercept))
	/((num_latencies-2)*(saccade_ss-(saccade_sum*saccade_sum)/num_latencies))));
  exit(0);
  }
