/*chisquare.c - do a Chi-square test on the z-transformed amplitudes output by
  ztrans.
  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"
extern double pochisq(double, int);

/*write the output file 'name', whose file header is 'erp', whose channel
  headers are ch_hdr[0..num_channels-1], and whose data is in
  data[0..num_channels-1][0..num_samples-1]*/
void write_output_file(name, erp, ch_hdr, data, num_channels, num_samples)
char *name, *erp, **ch_hdr;
double **data;
int num_channels, num_samples;
  {
  register int sample, channel;
  FILE *fp;
  if((fp = fopen(name,
#ifdef MSDOS
    "wb"
#else
    "w"
#endif
    )) == NULL)
    {
    perror(name);
    exit(errno);
    }
/*write the file header*/
  fwrite(erp, 1, packed_sizeof_SETUP, fp);
/*write the channel headers*/
  for(channel = 0; channel != num_channels; channel++)
    fwrite(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, fp);
/*write the channel data*/
  for(channel = 0; channel != num_channels; channel++)
    {
  /*5-byte dummy sweep header*/
    for(sample = 0; sample != 5; sample++)
      putc(0, fp);
  /*vector of floats*/
    for(sample = 0; sample != num_samples; sample++)
      write_Intel_float(data[channel][sample], fp);
    }
  fclose(fp);
  }

void main(argc, argv)
int argc;
char **argv;
  {
  register int sample, channel, file;
  register double z;
  register float calibration;
  int num_samples, num_channels;
  double **chisq, **zsum;
  FILE *avg;
  char **ch_hdr, erp[packed_sizeof_SETUP];
  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  if(argc < 5)
    {
    fprintf(stderr, "usage: %s <zscores.avg> ... <zscores.avg> <avg_zscores.avg> <probabilities.avg>\n", *argv);
    exit(1);
    }
/*open the last z-scores file*/
  if((avg = fopen(argv[argc-3],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[argc-3]);
    exit(1);
    }
/*read the file header*/
  if(fread(erp, 1, packed_sizeof_SETUP, avg) != packed_sizeof_SETUP)
    {
    perror(argv[argc-3]);
    fclose(avg);
    exit(errno);
    }
  num_channels = S_nchannels(erp);
  num_samples = S_pnts(erp);
/*read the channel headers*/
  ch_hdr = (char **)calloc(num_channels, sizeof(char *));
  for(channel = 0; channel != num_channels; channel++)
    {
    ch_hdr[channel] = (char *)malloc(packed_sizeof_ELECTLOC);
    if(fread(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, avg) != packed_sizeof_ELECTLOC)
      {
      perror(argv[1]);
      fclose(avg);
      exit(errno);
      }
    }
/*create storage for summed z-scores and for summed squared z-scores*/
  zsum = (double **)calloc(num_channels, sizeof(double *));
  chisq = (double **)calloc(num_channels, sizeof(double *));
  for(channel = 0; channel != num_channels; channel++)
    {
  /*get calibration from header*/
    calibration = EL_calib(ch_hdr[channel]);
  /*skip 5-byte sweep header*/
    fseek(avg, 5L, SEEK_CUR);
    zsum[channel] = (double *)calloc(num_samples, sizeof(double));
    chisq[channel] = (double *)calloc(num_samples, sizeof(double));
  /*read all z-scores for the current channel*/
    for(sample = 0; sample != num_samples; sample++)
      {
      z = calibration*read_Intel_float(avg);
      zsum[channel][sample] = z;
      chisq[channel][sample] = z*z;
      }
    }
/*done with the current file*/
  fclose(avg);
/*do the same for all other input files*/
  for(file = argc-4; file; file--)
    {
    if((avg = fopen(argv[file],
#ifdef MSDOS
      "rb"
#else
      "r"
#endif
      )) == NULL)
      {
      perror(argv[file]);
      exit(errno);
      }
  /*read the file header*/
    if(fread(erp, 1, packed_sizeof_SETUP, avg) != packed_sizeof_SETUP)
      {
      perror(argv[file]);
      fclose(avg);
      exit(errno);
      }
  /*make sure the number of channels and number of samples match*/
    if((S_nchannels(erp) != num_channels) || (S_pnts(erp) != num_samples))
      {
      fprintf(stderr, "%s (%d channels, %d samples) and %s (%d channels, %d samples) are incompatible.\n", argv[file], S_nchannels(erp), S_pnts(erp), argv[argc-4], num_channels, num_samples);
      fclose(avg);
      exit(1);
      }
  /*read the channel headers*/
    for(channel = 0; channel != num_channels; channel++)
      {
      if(fread(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, avg) != packed_sizeof_ELECTLOC)
	{
	perror(argv[file]);
	fclose(avg);
	exit(errno);
	}
      }
  /*read the channel data*/
    for(channel = 0; channel != num_channels; channel++)
      {
    /*get calibration factor from channel header*/
      calibration = EL_calib(ch_hdr[channel]);
    /*skip 5-byte sweep header*/
      fseek(avg, 5L, SEEK_CUR);
    /*read all z-scores for the current file*/
      for(sample = 0; sample != num_samples; sample++)
	{
	z = calibration*read_Intel_float(avg);
	zsum[channel][sample] += z;
	chisq[channel][sample] += z*z;
	}
      }
    fclose(avg);
    }
/*compute average z-scores, and Chi-square probabilities*/
  for(channel = 0; channel != num_channels; channel++)
    for(sample = 0; sample != num_samples; sample++)
      {
      zsum[channel][sample] /= argc-3;
      chisq[channel][sample] = pochisq(chisq[channel][sample], argc-4);
      }
/*variance data will not be saved, so clear the variance flag*/
  S_variance(erp) = 0;
  for(channel = 0; channel != num_channels; channel++)
  /*calibration has been multiplied in, so set the calibration factor to 1*/
    put_EL_calib(ch_hdr[channel], 1.0);
/*write average z-scores*/
  write_output_file(argv[argc-2], erp, ch_hdr, zsum, num_channels, num_samples);
/*write Chi-square probabilities*/
  write_output_file(argv[argc-1], erp, ch_hdr, chisq, num_channels, num_samples);
  for(channel = num_channels-1; channel >= 0; channel--)
    {
    free(chisq[channel]);
    free(zsum[channel]);
    free(ch_hdr[channel]);
    }
  free(chisq);
  free(zsum);
  free(ch_hdr);
  exit(0);
  }
