/*eyecomp.c - Principal Components transformation of EOG channels in a
  Scan386 continuous multiplexed (non-SynAmps) data (.CNT) file.
  Copyright (c) 1995 by Matthew Belmonte.

  Licence is hereby granted to any and all individuals or organisations to
  use and to modify this software for the purpose of nonprofit research only.
  Use of this software in conjunction or association with any for-profit
  endeavour requires specific, written permission from Matthew Belmonte.

  Licence is hereby granted to copy and to distribute this software, provided
  that such copies and distributions include this copyright and licence
  notice in its entirety, and provided that the software is not altered,
  edited, amended, embellished, or otherwise modified in any way.
  Distribution of this software in other than its original form requires
  specific, written permission from Matthew Belmonte.

  This software is provided 'as is', with no warranty, express or implied,
  including but not limited to warranties of merchantability or fitness for a
  particular purpose.  The user of this software assumes liability for any and
  all damages, whether direct or consequential, arising from its use.  The
  author of this software will not be liable for any such damages.

  This software makes use of the following copyrighted programs and header
  file from _Numerical Recipes in C_, second edition, by Press, Teukolsky,
  Vetterling, and Flannery:
	tred2.c
	tqli.c
	pythag.c
	nrutil.h
  These programs and header are not included in this distribution, and are
  available from Cambridge University Press.  The above licence for use,
  copying, and distribution DOES NOT apply to these programs.

  This software makes use of the following proprietary header file from the
  Scan386 EEG data processing system:
	sethead.h
  This file is not included in this distribution, and is available from
  NeuroScan, Inc.  The above licence for use, copying, and distribution
  DOES NOT apply to this file.
*/

/*If you've found this software useful, please tell me so.  Hopefully you'll
  be able to contact me at one of the following addresses:
  mkb4@cornell.edu
  mbelmonte@ucsd.edu
  matthew@salk.edu
  Neurosciences Department, UCSD, La Jolla, CA 92093-0608, USA.
  Neuropsychology 5056, Children's Hospital, San Diego, CA 92123-4282, USA.
  Check http://salk.edu/~matthew/ for more information.
  Also, I would appreciate acknowledgement in, and a copy of, any publication
  based on data processed by this program or any derivative of it.*/

/*Define MSDOS if you're compiling for MS-DOS.  Define M_I86 if you're
  running on an Intel processor (or any processor with backwards
  byte-order).  (If you're using Microsoft QuickC, these preprocessor
  constants will be automatically defined for you.)  This code should port
  to any system but has been verified only on an 80486 MS-DOS system and an
  R4000 IRIX system.  With Microsoft QuickC 2.00, the _Numerical Recipes_
  module tqli.c compiles incorrectly unless optimisations are turned off.
  Make sure your C compiler packs structures - the program includes a
  dynamic check for this.  If you want to preserve your original data,
  make a copy of it before running this program.  The program will overwrite
  your original eye channels after computing the principal component vectors
  and receiving confirmation from the user.  The program will not overwrite
  anything unless you tell it explictly to do so.*/

/*Makefile for Microsoft QuickC:
OPTIONS=/c /G2 /FPi87 /W0 /Zp

eyecomp.exe: eyecomp.obj tred2.obj tqli.obj pythag.obj
	link /BATCH /STACK:0x4000 /EXEPACK eyecomp.obj tred2.obj tqli.obj pythag.obj
eyecomp.obj: eyecomp.c sethead.h
	qcl $(OPTIONS) /Ox eyecomp.c
tred2.obj: tred2.c
	qcl $(OPTIONS) /Ox tred2.c
tqli.obj: tqli.c nrutil.h
	qcl $(OPTIONS) /Gs tqli.c
pythag.obj: pythag.c nrutil.h
	qcl $(OPTIONS) /Ox pythag.c
*/

#define NCHANS 6 /*number of eye channels (must be lowest-numbered channels)*/
#define MAX_FLATLINE_LENGTH 21.0 /*length of amp. saturation test (ms)*/
#define EYE_SCALE 4.0 /*coarse enough to prevent clipping*/
#define VOLTAGE_MULTIPLIER 204.8 /*see page Headers-5 of the Scan 3.0 manual*/

#include <stdio.h>
#include <errno.h>
#ifndef MSDOS
# include <unistd.h>
# include <values.h>
# include <sys/param.h>
#else
# define MAXPATHLEN 80
#endif

FILE *cnt; /*continuous EEG file*/

/*routines from _Numerical Recipes in C_:*/
extern void tred2(float**, int, float*, float*);
extern void tqli(float*, float*, int, float**);
/*a simpler substitute for the NRC error handler:*/
void nrerror(char *message)
  {
  fprintf(stderr, "%s\n", message);
  fclose(cnt);
  exit(1);
  }

typedef int Boolean;
#define false 0
#define true 1

#define N_ELECT 0
#define int long                /*Scan386 code denotes longs by "int"*/
#include "sethead.h"            /*Scan386 header file*/
#undef int
#undef N_ELECT

#ifdef M_I86                    /*Scan386 data is in LSB:MSB byte order*/
#  define swap_bytes(word) word
#  define reverse_long(double_word) double_word
#  define reverse_float(double_word) double_word
#  define reverse_double(quad_word) quad_word
#else

short swap_bytes(word)
short word;
  {
  return(((word >> 8) & 0xff)|(word << 8));
  }

void reverse_bytes(src, dest, len)
char *src, *dest;
int len;
  {
  register int i;
  for(i = 0; i != len; i++)
    dest[i] = src[len-1-i];
  }

long reverse_long(double_word)
long double_word;
  {
  long result;
  reverse_bytes((char *)&double_word, (char *)&result, sizeof(long));
  return(result);
  }

float reverse_float(double_word)
float double_word;
  {
  float result;
  reverse_bytes((char *)&double_word, (char *)&result, sizeof(long));
  return(result);
  }

double reverse_double(quad_word)
double quad_word;
  {
  double result;
  reverse_bytes((char *)&quad_word, (char *)&result, sizeof(long));
  return(result);
  }
#endif

short read_swapped_word(stream)
FILE *stream;
  {
  short result;
  result = getc(stream);
  result |= getc(stream) << 8;
  return(result);
  }

float dotproduct(x, y, n)
float *x, *y;
int n;
  {
  register int j;
  float sum;
  sum = 0.0;
  for(j=1; j<=n; j++)
    sum += x[j]*y[j];
  return(sum);
  }

void main(argc, argv)
int argc;
char *argv[];
  {
  register int column, channel;
  int num_channels, max_flatline_length, *flatline_length, *new_baselines;
  register long sample;
  long num_samples, num_useable_samples, num_overflows;
  Boolean compute_covariance;
  char filename[MAXPATHLEN];
  SETUP erp;
/*vectors and matrices*/
  ELECTLOC *channel_header;
  double **summed_products, *summed_eog;
  float **covariance_matrix,
	**eigenvectors,
	*eigenvalues,
	*subdiagonals,
	*calibration,
	*eog;
/*data areas*/
  int f[NCHANS];
  ELECTLOC derivations[NCHANS];
  double doubles[(NCHANS+1)*NCHANS], *(dbl_pointers[NCHANS]);
  float data[(NCHANS+4)*NCHANS], *(pointers[NCHANS]);
/*data pointers - NRC routines use 1 rather than 0 as lowest subscript*/
  flatline_length = f-1;
  new_baselines = f-1;
  channel_header = derivations-1;
  summed_eog = doubles-1;
  summed_products = dbl_pointers-1;
  covariance_matrix = pointers-1;
  eigenvectors = pointers-1;
  eigenvalues = data-1+NCHANS*NCHANS;
  subdiagonals = data-1+(NCHANS+1)*NCHANS;
  calibration = data-1+(NCHANS+2)*NCHANS;
  eog = data-1+(NCHANS+3)*NCHANS;
  if((sizeof(SETUP) != 900) || (sizeof(ELECTLOC) != 75))
    {
    fprintf(stderr, "Your C compiler doesn't pack structures - recompile with packed structures.\n");
    exit(2);
    }
  if(argc != 2)
    {
    fprintf(stderr, "Usage: %s <continuous file prefix>\n", argv[0]);
    exit(2);
    }
  sprintf(filename, "%s.cnt", argv[1]);
  cnt = fopen(filename,
#ifdef MSDOS
	"rb+"
#else
	"r+"
#endif
	);
  if(cnt == NULL)
    {
    perror(filename);
    exit(errno);
    }
  if(fread(&erp, sizeof(erp), 1, cnt) != 1) /*read main header*/
    {
    perror(filename);
    fclose(cnt);
    exit(errno);
    }
  num_channels = swap_bytes(erp.nchannels);
  if(num_channels < NCHANS)
    {
    fprintf(stderr, "Looking for %d eye channels, found only %d total channels!\n", NCHANS, num_channels);
    fclose(cnt);
    exit(1);
    }
  num_samples = (reverse_long(erp.EventTablePos)-sizeof(SETUP)-num_channels*sizeof(ELECTLOC))/(num_channels*sizeof(short));
  max_flatline_length = (int)(((float)(erp.rate*MAX_FLATLINE_LENGTH))/1000.0);
  for(channel = 1; channel <= NCHANS; channel++)
    {
    covariance_matrix[channel] /*==eigenvectors[channel]*/ = data-1+(channel-1)*NCHANS;
    summed_products[channel] = doubles-1+channel*NCHANS;
    eog[channel] = 32767.0; /*to initialise for flatline_length computation*/
    covariance_matrix[channel][channel] = 1.0; /*covariance of channel w/ self is always 1*/
    for(column = 1; column <= channel; column++) /*zero upper right triangle*/
      summed_products[channel][column] = 0.0;
    summed_eog[channel] = 0.0;
    fread(channel_header+channel, sizeof(*channel_header), 1, cnt);
    calibration[channel] = reverse_float(channel_header[channel].sensitivity)*reverse_float(channel_header[channel].calib)/VOLTAGE_MULTIPLIER;
    }
  num_useable_samples = 0;
  for(sample = 0; sample != num_samples; sample++)
    {
    compute_covariance = true;
    for(channel = 1; channel <= NCHANS; channel++)
      {
      float temp;
      temp = calibration[channel]*read_swapped_word(cnt);
      flatline_length[channel] = (eog[channel]==temp? 1+flatline_length[channel]: 0);
      eog[channel] = temp;
      compute_covariance = (compute_covariance && (flatline_length[channel] < max_flatline_length));
      }
    if(compute_covariance)
      {
      for(channel = 1; channel <= NCHANS; channel++)
	{
	for(column = 1; column <= channel; column++)
	  summed_products[channel][column] += eog[channel]*eog[column];
	summed_eog[channel] += eog[channel];
	}
      num_useable_samples++;
      }
    }
/*compute covariance from sums*/
  for(channel = 1; channel <= NCHANS; channel++)
    for(column = 1; column < channel; column++)
      {
      covariance_matrix[channel][column] = (summed_products[channel][column] - summed_eog[channel]*summed_eog[column]/num_useable_samples)/num_useable_samples;
      covariance_matrix[column][channel] = covariance_matrix[channel][column];
      }
/*for the diagonal, this reduces to simple variance:*/
  for(channel = 1; channel <= NCHANS; channel++)
    covariance_matrix[channel][channel] = (summed_products[channel][channel] - summed_eog[channel]*summed_eog[channel]/num_useable_samples)/num_useable_samples;
/*compute principal components*/
  tred2(covariance_matrix, NCHANS, eigenvalues, subdiagonals);
  tqli(eigenvalues, subdiagonals, NCHANS, eigenvectors);
/*sort in order of decreasing variance
 - not many values to sort, so we can use a simple O(N**2) routine*/
  for(channel = 1; channel <= NCHANS; channel++)
  /*inv: eigenvalues[1..channel-1] >= eigenvalues[channel..NCHANS], and
    eigenvalues[i] >= eigenvalues[i+1] for 1<=i<=channel-1.  eigenvectors[]
    has been permuted identically to eigenvalues[].*/
    {
    float max_eigenvalue, *temp;
    int max_eigenvalue_index;
    max_eigenvalue = eigenvalues[channel];
    for(column = channel+1; column <= NCHANS; column++)
      /*inv: max_eigenvalue is the maximum of eigenvalues[channel..column-1].*/
      if(eigenvalues[column] > max_eigenvalue)
	{
	max_eigenvalue = eigenvalues[column];
	max_eigenvalue_index = column;
	}
    eigenvalues[max_eigenvalue_index] = eigenvalues[channel];
    eigenvalues[channel] = max_eigenvalue;
    temp = eigenvectors[max_eigenvalue_index];
    eigenvectors[max_eigenvalue_index] = eigenvectors[channel];
    eigenvectors[channel] = temp;
    printf("PRINCIPAL COMPONENT %d (variance %f)\n", channel, max_eigenvalue);
    for(column = 1; column <= NCHANS; column++)
      printf("(%s %.3f) ", channel_header[column].lab, eigenvectors[channel][column]);
    putchar('\n');
    }
/*update*/
  printf("Permanently replace eye channels with these components (Y/N)? ");
  if((getchar() & 0xdf) != 'Y')
    {
    fclose(cnt);
    printf("File was not altered.\n");
    exit(1);
    }
  fseek(cnt, (long)(sizeof(SETUP)), SEEK_SET);
  for(channel = 1; channel <= NCHANS; channel++)
    {
    float new_baseline;
    for(channel = 1; channel <= NCHANS; channel++)
      {
      new_baseline = 0.0;
      for(column = 1; column <= NCHANS; column++)
	new_baseline += eigenvectors[channel][column]*swap_bytes(channel_header[column].baseline)*reverse_float(channel_header[column].sensitivity)*reverse_float(channel_header[column].calib);
      new_baseline *= VOLTAGE_MULTIPLIER/EYE_SCALE;
      }
    new_baselines[channel] = (int)new_baseline;
    }
  for(channel = 1; channel <= NCHANS; channel++)
    {
    channel_header[channel].sensitivity = reverse_float(1.0);
    channel_header[channel].calib = reverse_float(EYE_SCALE);
    channel_header[channel].baseline = swap_bytes(new_baselines[channel]);
    channel_header[channel].lab[0] = 'E';
    channel_header[channel].lab[1] = 'O';
    channel_header[channel].lab[2] = 'G';
    channel_header[channel].lab[3] = '0'+channel;
    channel_header[channel].lab[4] = '\0';
    fwrite(channel_header+channel, sizeof(*channel_header), 1, cnt);
    }
  fseek(cnt, (long)(sizeof(SETUP)+num_channels*sizeof(ELECTLOC)), SEEK_SET);
  num_overflows = 0L;
  for(sample = 0; sample != num_samples; sample++)
    {
    for(channel = 1; channel <= NCHANS; channel++)
      eog[channel] = calibration[channel]*read_swapped_word(cnt);
    fseek(cnt, (long)(-NCHANS*sizeof(short)), SEEK_CUR);
    for(channel = 1; channel <= NCHANS; channel++)
      {
      float temp;
      short buf;
      temp = VOLTAGE_MULTIPLIER/EYE_SCALE*dotproduct(eog, eigenvectors[channel], NCHANS);
      if(temp > 32767.0)
	{
	temp = 32767.0;
	num_overflows++;
	}
      else if(temp < -32768.0)
	{
	temp = -32768.0;
	num_overflows++;
	}
      buf = swap_bytes((short)temp);
      fwrite(&buf, sizeof(short), 1, cnt);
      }
    fseek(cnt, (long)((num_channels-NCHANS)*sizeof(short)), SEEK_CUR);
    }
  fclose(cnt);
  if(num_overflows != 0L)
    {
    fprintf(stderr, "%ld overflows (%ld%% of samples)\n", num_overflows, 100L*num_overflows/num_samples);
    fprintf(stderr, "If this is too many, recompile with a larger EYE_SCALE.\n");
    }
  exit(0);
  }