/*mannwhitney.c - The Mann-Whitney 'U' test, applied pointwise to two groups 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 "gnuro.h"
#include "intel.h"
/*#define free(x) DEBUG*/

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 sample, channel, file;
  int	m, n,			/*group sizes*/
	m_ranksum, n_ranksum,	/*rank sums for each group*/
	*group;			/*ptr to group flags - 0 or 1 for each datum*/
  float **calibration, *eeg;
  FILE **avg;
  int num_samples, num_channels;
  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  if(argc < 4)
    {
    fprintf(stderr, "usage: %s <file.avg> ... - <file.avg> ...\n", *argv);
    exit(1);
    }
  for(m = 2; (m != argc-1) && !((argv[m][0] == '-') && (argv[m][1] == '\0')); m++)
    ;
  if(m == argc-1)
    {
    fprintf(stderr, "usage: %s <file.avg> ... - <file.avg> ...\n", *argv);
    exit(1);
    }
  n = argc-1-m--;
  avg = (FILE **)calloc(m+n, sizeof(FILE *));
  calibration = (float **)calloc(m+n, sizeof(float *));
  if((avg[0] = fopen(argv[1], "r")) == NULL)
    {
    perror(argv[1]);
    exit(errno);
    }
  fseek(avg[0], (long)S_pnts_offset, SEEK_SET);
  num_samples = read_Intel_short(avg[0]);
  num_channels = read_Intel_short(avg[0]);
  for(file = 0; file != m+n; file++)
    calibration[file] = (float *)calloc(num_channels, sizeof(float));
  read_calibrations(avg[0], calibration[0], num_channels);
  for(file = 1; file != m; file++)
    {
    avg[file] = open_file(argv[1+file], num_channels, num_samples);
    read_calibrations(avg[file], calibration[file], num_channels);
    }
  for(file = 0; file != n; file++)
    {
    avg[m+file] = open_file(argv[2+m+file], num_channels, num_samples);
    read_calibrations(avg[m+file], calibration[m+file], num_channels);
    }
  eeg = (float *)calloc(m+n, sizeof(float));
  group = (int *)calloc(m+n, sizeof(int));
  printf("m=%d n=%d\n", m, n);
  for(channel = 0; channel != num_channels; channel++)
    {
    printf("CHANNEL %d:\nsample Sm Sn   Um Un\n", channel);
    for(file = 0; file != m+n; file++)
      fseek(avg[file], 5L, SEEK_CUR);
    for(sample = 0; sample != num_samples; sample++)
      {
      for(file = 0; file != m+n; file++)
	eeg[file] = calibration[file][channel]*read_Intel_float(avg[file]);
    /*group[file]=0 if file is in 1st group, 1 if file is in 2nd group*/
      for(file = 0; file != m; file++)
	group[file] = 0;
      for(file = 0; file != n; file++)
	group[m+file] = 1;
      sortrecords(eeg, group, m+n);
      m_ranksum = 0;
      n_ranksum = 0;
      for(file = 1; file <= m+n; file++)
	{
	if(group[file-1])
	  n_ranksum += file;
	else
	  m_ranksum += file;
	}
      printf("%4d   %2d %2d   %2d %2d\n", sample, m_ranksum, n_ranksum, m*n+m*(m+1)/2-m_ranksum, m*n+n*(n+1)/2-n_ranksum);
      }
    }
  free(group);
  free(eeg);
  for(channel = num_channels-1; channel >= 0; channel--)
    free(calibration[channel]);
  free(calibration);
  free(avg);
  exit(0);
  }
