/*signtest.c - sign test on a group of GnuroScan averaged EEG files.
  Based on ttest.c by Matthew Belmonte.
  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 <math.h>
#include <errno.h>
#include "intel.h"
#include "gnuro.h"

/*compute two-tailed cumulative binomial probabilities for n trials at even
  chances*/
void compute_binomial_probs(b, n)
double *b;
int n;
  {
  register double term;
  register int k;
  term = 1.0;
/*base probability is 1/2**N, times two for two-tailed*/
  for(k = 0; k < n-1; k++)
    term /= 2.0;
  b[0] = term;
  k = 0;
  printf("BINOMIAL PROBABILITIES FOR n=%d:\n\t k\tindiv.\tcum.\n\t 0\t%f\t%f\n", n, term, term);
  while(k < n/2-1)
  /*inv: b[0..k] is the two-tailed binomial probability of the sign count of n
    signs totalling k or less*/
    {
    term *= (n-k)/(double)(k+1);
    k++;
    b[k] = term+b[k-1];
    printf("\t%2d\t%f\t%f\n", k, term, b[k]);
    }
  b[k+1] = 1.0;
  printf("\t%2d\t%f\t1.0\n", k+1, term*(n-k)/(double)(k+1));
  }

/*open the named file, and read its main header into 'erp'*/
FILE *open_file(name, erp)
char *name, *erp;
  {
  FILE *avg;
  if((avg = fopen(name,
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(name);
    exit(errno);
    }
  fread(erp, 1, packed_sizeof_SETUP, avg);
  return(avg);
  }

/*read the channel headers from 'avg' into ch_hdr[0..num_channels-1]*/
void read_headers(avg, ch_hdr, num_channels)
FILE *avg;
char **ch_hdr;
int num_channels;
  {
  register int channel;
  for(channel = 0; channel != num_channels; channel++)
    fread(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, avg);
  }

/*sum the data from 'avg' into s[0..num_channels-1][0..num_samples-1]
  (sign count)*/
void read_data(avg, s, num_channels, num_samples)
FILE *avg;
int **s;
int num_channels, num_samples;
  {
  register int sample, channel;
  for(channel = 0; channel != num_channels; channel++)
    {
  /*seek past 5-byte sweep header*/
    fseek(avg, 5L, SEEK_CUR);
  /*read one channel of data*/
    for(sample = 0; sample != num_samples; sample++)
      s[channel][sample] += (read_Intel_float(avg) > 0.0);
    }
  }

void main(argc, argv)
int argc;
char **argv;
  {
  register int sample, channel, file;
  int num_samples, num_channels, num_files;
  FILE *avg;
  char **ch_hdr, erp[packed_sizeof_SETUP];
  double *binomial_table;
  int **s;
  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  if(argc < 4)
    {
    fprintf(stderr, "usage: %s <file.avg> ... <file.avg> <probabilities.avg>\n", *argv);
    exit(1);
    }
  num_files = argc-2;
  avg = open_file(argv[1], erp);
  num_channels = S_nchannels(erp);
  num_samples = S_pnts(erp);
  ch_hdr = (char **)calloc(num_channels, sizeof(char *));
  s = (int **)calloc(num_channels, sizeof(double *));
  for(channel = 0; channel != num_channels; channel++)
    {
    ch_hdr[channel] = (char *)malloc(packed_sizeof_ELECTLOC);
    s[channel] = (int *)calloc(num_samples, sizeof(double));
    for(sample = 0; sample != num_samples; sample++)
       s[channel][sample] = 0;
    }
  read_headers(avg, ch_hdr, num_channels);
  read_data(avg, s, num_channels, num_samples);
  fclose(avg);
  for(file = 2; file != argc-1; file++)
    {
    avg = open_file(argv[file], erp);
    if((S_nchannels(erp) != num_channels) || (S_pnts(erp) != num_samples))
      {
      fprintf(stderr, "%s has incorrect dimensions (%d x %d)\n", argv[file], S_nchannels(erp), S_pnts(erp));
      exit(1);
      }
    read_headers(avg, ch_hdr, num_channels);
    read_data(avg, s, num_channels, num_samples);
    fclose(avg);
    }
/*compute binomial probabilities*/
  binomial_table = (double *)(calloc(num_files/2+1, sizeof(double)));
  compute_binomial_probs(binomial_table, num_files);
/*create output file*/
  if((avg = fopen(argv[argc-1],
#ifdef MSDOS
    "wb"
#else
    "w"
#endif
    )) == NULL)
    {
    perror(argv[argc-1]);
    exit(errno);
    }
/*no individual variance data included*/
  S_variance(erp) = 0;
/*write the file header*/
  fwrite(erp, 1, packed_sizeof_SETUP, avg);
/*write the channel headers, with all calibrations set to unity*/
  for(channel = 0; channel != num_channels; channel++)
    {
    put_EL_calib(ch_hdr[channel], 1.0);
    fwrite(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, avg);
    }
  for(channel = 0; channel != num_channels; channel++)
    {
  /*write dummy 5-byte sweep header*/
    for(sample = 0; sample != 5; sample++)
      putc(0, avg);
  /*print channel number*/
    printf("\nCHANNEL %d:\n", channel);
  /*write probabilities for this channel*/
    for(sample = 0; sample != num_samples; sample++)
      {
    /*print raw score*/
      printf("%4d %2d ", sample, s[channel][sample]);
    /*adjust for symmetry of the binomial distribution*/
      if(s[channel][sample] > num_files/2)
	s[channel][sample] = num_files - s[channel][sample];
    /*look up binomial probability and write it to the file*/
      write_Intel_float(binomial_table[s[channel][sample]], avg);
    /*print binomial probability*/
      printf("%f\n", binomial_table[s[channel][sample]]);
      }
    }
  fclose(avg);
  exit(0);
  }
