/*friedman.c - Friedman's rank test for a randomised block design, applied
  pointwise to a set of Gnuroscan averaged EEG files.
  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 <sys/param.h>
#include "gnuro.h"
#include "intel.h"
extern double pochisq(double, int);

void process_args(argc, argv, num_subjs, num_treatments)
int argc;
char **argv;
int	*num_subjs,		/*size of experimental group*/
	*num_treatments;	/*number of treatments*/
  {
  register int arg;		/*argument counter*/
  for(arg = 1; (arg != argc) && !((argv[arg][0] == '-')&&(argv[arg][1] == '\0')); arg++)
    ;
  if((arg == 1) || (arg >= argc-1))
    {
    fprintf(stderr, "usage: %s <SubjDir> ... <SubjDir> - <TreatmentName.avg> ... <TreatmentName.avg>\nThe dash ( '-' ) is a separator.  Each SubjDir is the name of a directory\nwhere the data set for one subject can be found.  Each TreatmentName is the\nname of an averaged EEG file, present in each directory, that corresponds to a\nparticular treatment (for example, different rates of presentation of stimuli).\n", *argv);
    exit(1);
    }
  *num_subjs = arg-1;
  *num_treatments = argc-arg-1;
  }

/*Return a pathname consisting of the given directory name and file name.  The
  pathname is in static storage and thus will be useable until the next
  invocation of make_pathname().*/
char *make_pathname(dir, file)
char *dir, *file;
  {
  static char pathname[MAXPATHLEN];
  register char *s;
  s = pathname;
  while(*dir != '\0')
    *(s++) = *(dir++);
  *(s++) = '/';
  do
    *(s++) = *file;
  while(*(file++) != '\0');
  return(pathname);
  }

FILE *open_file(filename, num_channels, num_samples)
char *filename;
int num_channels, num_samples;
  {
  register FILE *fp;
  int c, s;
  if((fp = fopen(filename,
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(filename);
    exit(errno);
    }
  fseek(fp, (long)S_pnts_offset, SEEK_SET);
  s = read_Intel_short(fp);
  c = read_Intel_short(fp);
  if((c != num_channels) || (s != num_samples))
    {
    fclose(fp);
    fprintf(stderr, "%s: bad dimensions %dx%d\n", filename, c, s);
    exit(1);
    }
  return(fp);
  }

/*Read calibration factors for all channels from the given file, and leave the
  file's input pointer at the end of the last channel header.*/
void read_calibrations(fp, c, num_channels)
FILE *fp;
float *c;
int num_channels;
  {
  register int channel;
  fseek(fp, (long)packed_sizeof_SETUP, SEEK_SET);
  for(channel = 0; channel != num_channels; channel++)
    {
    fseek(fp, (long)EL_calib_offset, SEEK_CUR);
    c[channel] = read_Intel_float(fp);
    }
  }

void swap(b, c, i, j)
float *b;
int *c;
int i, j;
  {
  register float temp;
  register int tempi;
  temp = b[i];
  tempi = c[i];
  b[i] = b[j];
  c[i] = c[j];
  b[j] = temp;
  c[j] = tempi;
  }

void sortrecords(b, c, count)
float *b;	/*key record*/
int *c;		/*companion record*/
int count;	/*number of records*/
  {
  register int h, k, m, n;
  m = 0;
  n = count-1;
  /*inv: b[m..n] is a permutation of the subsequence
    b[m..n] of sorted sequence b[0..count-1]*/
  while(m < n)
    {
    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, c, h++, k--);
      } /*elihw*/
    /*k=h-1*/
    if(k==m)
      m++;
    else if(h==n+1)
      swap(b, c, m, n--);
    else
      {
      swap(b, c, m, k);
      sortrecords(b+m, c+m, k-m);
      m = k+1;
      } /*fi*/
    /*b[k] is now in sorted position*/
    }
  }

void main(argc, argv)
int argc;
char **argv;
  {
  register int	n,			/*temporary*/
		treatment, subject,	/*counters*/
		sample, channel;
  register double chi_square;		/*Chi-square value*/
  float		*cell_obs;		/*subject x treatment observations*/
  int		*group,			/*treatment group of each datum*/
		*ranksum;		/*rank sum for each treatment*/
  int num_subjs, num_treatments, num_samples, num_channels;
  float ***calibration;			/*subj x treatment x channel calibs.*/
  FILE ***avg;				/*subj x treatment input pointers*/

  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
/*process command-line parameters and check them for consistency*/
  process_args(argc, argv, &num_subjs, &num_treatments);
  if((num_subjs <= 5) && (num_treatments <= 5))
    fprintf(stderr, "warning: Chi-square may be a poor approximation\n");
/*allocate data areas for each subject*/
  avg = (FILE ***)calloc(num_subjs, sizeof(FILE **));
  for(subject = 0; subject != num_subjs; subject++)
    avg[subject] = (FILE **)calloc(num_treatments, sizeof(FILE *));
  calibration = (float ***)calloc(num_subjs, sizeof(float **));
/*get input dimensions*/
  if((avg[0][0] = fopen(make_pathname(argv[1], argv[1+num_subjs+1]), "r")) == NULL)
    {
    perror(argv[1]);
    exit(errno);
    }
  fseek(avg[0][0], (long)S_pnts_offset, SEEK_SET);
  num_samples = read_Intel_short(avg[0][0]);
  num_channels = read_Intel_short(avg[0][0]);
  fclose(avg[0][0]);
  for(subject = 0; subject != num_subjs; subject++)
    {
    calibration[subject] = (float **)calloc(num_treatments, sizeof(float*));
    for(treatment = 0; treatment != num_treatments; treatment++)
      calibration[subject][treatment] = (float *)calloc(num_channels, sizeof(float));
    }
/*compute statistics*/
  for(treatment = 0; treatment != num_treatments; treatment++)
    {
    for(subject = 0; subject != num_subjs; subject++)
      {
      avg[subject][treatment] = open_file(make_pathname(argv[1+subject], argv[1+num_subjs+1+treatment]), num_channels, num_samples);
      read_calibrations(avg[subject][treatment], calibration[subject][treatment], num_channels);
      }
    }
  cell_obs = (float *)calloc(num_treatments, sizeof(float));
  group = (int *)calloc(num_treatments, sizeof(int));
  ranksum = (int *)calloc(num_treatments, sizeof(int));
  for(channel = 0; channel != num_channels; channel++)
    {
  /*print banner*/
    printf("CHANNEL %d:\nsmpl Chi2(%d) p\n", channel, num_treatments-1);
  /*skip 5-byte data header in all files*/
    for(subject = 0; subject != num_subjs; subject++)
      for(treatment = 0; treatment != num_treatments; treatment++)
	fseek(avg[subject][treatment], 5L, SEEK_CUR);
  /*separate, pointwise computations - one for each sample*/
    for(sample = 0; sample != num_samples; sample++)
      {
    /*read observations from data files*/
      for(treatment = 0; treatment != num_treatments; treatment++)
	ranksum[treatment] = 0;
      for(subject = 0; subject != num_subjs; subject++)
	{
	for(treatment = 0; treatment != num_treatments; treatment++)
	  {
	  group[treatment] = treatment;
	  cell_obs[treatment] = calibration[subject][treatment][channel]*read_Intel_float(avg[subject][treatment]);
	  }
	sortrecords(cell_obs, group, num_treatments);
	for(treatment = 0; treatment != num_treatments; treatment++)
	  ranksum[group[treatment]] += 1+treatment;
	  }
      chi_square = 0.0;
      for(treatment = 0; treatment != num_treatments; treatment++)
	chi_square += ((double)(ranksum[treatment]))*(double)(ranksum[treatment]);
      chi_square *= 12.0/(num_subjs*num_treatments*(num_treatments+1));
      chi_square -= 3*num_subjs*(num_treatments+1);
      printf("%4d %7.3f %5.3f\n", sample, chi_square, 1.0-pochisq(chi_square, num_treatments-1));
      } /*sample*/
    } /*channel*/
  free(ranksum);
  free(group);
  free(cell_obs);
  for(subject = num_subjs-1; subject >= 0; subject--)
    {
    for(treatment = num_treatments-1; treatment >= 0; treatment--)
      free(calibration[subject][treatment]);
    free(calibration[subject]);
    free(avg[subject]);
    }
  free(calibration);
  free(avg);
  exit(0);
  }
