/*fplotavg.c - translate into PostScript code a GnuroScan averaged EEG file that
  has been generated by the fastshft experiment.  Based on plotavg.c.
  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>
#ifdef NOFLOATS
#  define fceil(x) ((float)(ceil((double)x)))
#  define ffloor(x) ((float)(floor((double)x)))
#endif
#include <errno.h>
#include "intel.h"
#include "gnuro.h"
#include "fastshft.h"

#define FI_MAX 4
#define CH_MAX 16
#define NUM_DIFFERENCE_SAMPLES 1000

int fontsize, num_samples;
float lmargin, rmargin, tmargin, bmargin, page_height, page_width, plot_height, plot_leading, title_height, tick_length;
float xmin, xmax, image_width;

#include <math.h>
#include "fmedian.h"

/*Return the number of items in distrib[0..k-1], which is sorted in
  nondecreasing order, that are strictly less than 'x'.*/
int rank(x, distrib, k)
float x, *distrib;
int k;
  {
  register int lb, ub;
  lb = 0;
  ub = k;
  while(lb < ub)
  /*inv: distrib[0..lb-1] < x <= distrib[ub..k-1]*/
    {
    k = (lb+ub)/2;
    if(distrib[k] < x)
      lb = k+1;
    else
      ub = k;
    }
  return(lb);
  }

void draw_axes(channel, label, ymin, ymax)
int channel;
char *label;
float ymin, ymax;
  {
  float x_axis_elevation, y_axis_elevation, hundreds_of_msecs;
  x_axis_elevation = page_height-tmargin-title_height-channel*plot_height-(ymax<=0.0? 0.0: ymin>=0.0? 1.0: ymax/(ymax-ymin))*(plot_height-plot_leading)-plot_leading/2.0;
  y_axis_elevation = (xmax<=0.0)? page_width-rmargin: (xmin>=0.0)? lmargin: lmargin-(xmin/(xmax-xmin))*image_width;
  printf("newpath\n2 setlinewidth\n0.75 setgray\n");
  printf("%f %f moveto\n", lmargin, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0);
  printf("(%s) show\n", label);
  printf("%f %f moveto\n", lmargin, x_axis_elevation);
  printf("%f 0 rlineto\n", image_width);
  printf("%f %f moveto\n", y_axis_elevation, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0);
  printf("0 %f rlineto\n", plot_leading-plot_height);
  printf("(%.0fuV) show\n", ymin);
  printf("%f %f moveto\n", y_axis_elevation, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0);
  printf("(%.0fuV) show\n", ymax);
  for(hundreds_of_msecs = fceil(10.0*xmin); hundreds_of_msecs != ffloor(10.0*xmax); hundreds_of_msecs += 1.0)
    printf("%f %f moveto\n0 %f rlineto\n", y_axis_elevation+(hundreds_of_msecs/(10.0*(xmax-xmin)))*image_width, x_axis_elevation-tick_length/2.0, tick_length);
  printf("stroke\n0 setgray\n");
  }

/*Plot points epoch[0..epoch_len-1], at the given scale, comparing the
  difference (diff_epoch1-diff_epoch0)[0..epoch_len-1] to the distribution
  distrb[0..distrb_len-1] in order to ascertain significance.*/
void plotpoints(epoch, diff_epoch0, diff_epoch1, epoch_len, time_zero, distrb, distrb_len, scale)
float *epoch, *diff_epoch0, *diff_epoch1;
int epoch_len;
float *distrb;
int distrb_len;
float scale;
  {
  register int sample, zero_crossing, prev_zero_crossing;
  register float cur_diff, prev_diff;
  register double integral, min_p, max_p, p, grey_level, prev_grey_level;
/*clip significance range s.t. baseline is excluded from significance*/
  min_p = 0.025;
  max_p = 0.975;
  zero_crossing = 1+time_zero/2;
  cur_diff = diff_epoch1[time_zero/2]-diff_epoch0[time_zero/2];
  while(zero_crossing != time_zero)
    {
    prev_zero_crossing = zero_crossing;
    integral = 0.0;
    do
      {
      integral += cur_diff;
      zero_crossing++;
      prev_diff = cur_diff;
      cur_diff = diff_epoch1[zero_crossing]-diff_epoch0[zero_crossing];
      }
    while((zero_crossing != time_zero) && (cur_diff*prev_diff > 0));
    p = rank(integral, distrb, distrb_len)/(float)distrb_len;
    if(p == 0.0)
      {
      distrb[0] = integral;
      min_p = 1.0/distrb_len;
      }
    else if(p == 1.0)
      {
      distrb[distrb_len-1] = integral;
      max_p = (distrb_len-1)/(float)distrb_len;
      }
    else if(p < min_p)
      min_p = p;
    else if(p > max_p)
      max_p = p;
    }
/*now use the clipped significance range for the rest of the computation*/
  prev_grey_level = 0.0;
  zero_crossing = 1;
  cur_diff = diff_epoch1[0]-diff_epoch0[0];
  sample = 0;
  printf("1 setlinewidth\n0 %f rmoveto\n", epoch[0]*scale);
  while(zero_crossing != epoch_len-1)
    {
    prev_zero_crossing = zero_crossing;
    integral = 0.0;
    do
      {
      integral += cur_diff;
      zero_crossing++;
      prev_diff = cur_diff;
      cur_diff = diff_epoch1[zero_crossing]-diff_epoch0[zero_crossing];
      }
    while((zero_crossing != epoch_len-1) && (zero_crossing != time_zero) && (cur_diff*prev_diff > 0));
  /*the added conjunct in the guard above produces an artificial 'zero-crossing'
    at time zero - this separates any spurious deviation immediately preceding
    time zero from any significant deviation immediately following time zero
    and having the same polarity*/
    p = rank(integral, distrb, distrb_len)/(float)distrb_len;
    grey_level = (((p > max_p) || (p < min_p))? 0.0: 0.5);
    if((epoch == diff_epoch0) && (integral != 0.0) && (zero_crossing != time_zero) && (zero_crossing != epoch_len-1))
      fprintf(stderr, "\t%s-going zero-crossing at %.0fms\n", ((integral>0)?"positive":"negative"), 0.5+1000.0*((zero_crossing-time_zero)*(xmax-xmin)/num_samples));
    if(grey_level != prev_grey_level)
      {
      printf("currentpoint\nstroke\nmoveto\n%f setgray\n", grey_level);
      prev_grey_level = grey_level;
      }
    for(sample = prev_zero_crossing; sample != zero_crossing; sample++)
      printf("%f %f rlineto\n", image_width/(epoch_len-1), (epoch[sample]-epoch[sample-1])*scale);
    }
  printf("%f %f rlineto\n", image_width/(epoch_len-1), (epoch[epoch_len-1]-epoch[epoch_len-2])*scale);
  printf("stroke\n0 setgray\n");
  }

/*Read a vector of floats, of length 'num_samples', from file 'fp' into array
  'epoch'.  Return the index of the leftmost nonzero element.*/
int read_channel(fp, epoch, num_samples, calibration, min, max)
FILE *fp;
float *epoch;
int num_samples;
float calibration, *min, *max;
  {
  register int sample, leftmost_nonzero_sample;
  leftmost_nonzero_sample = num_samples;
  for(sample = 0; sample != num_samples; sample++)
    {
    epoch[sample] = read_Intel_float(fp)*calibration;
    if((sample < leftmost_nonzero_sample) && (epoch[sample] != 0.0))
      leftmost_nonzero_sample = sample;
    if(epoch[sample] < *min)
      *min = epoch[sample];
    if(epoch[sample] > *max)
      *max = epoch[sample];
    }
  return(leftmost_nonzero_sample);
  }

/*Read a 'num_channels'x'num_samples' matrix of floats from file 'fp' into the
  two-dimensional array 'epoch'.  Make 'startpos' the maximum index, over all
  the channels, at which the leftmost nonzero sample occurs.*/
void read_channels(fp, epoch, num_samples, num_channels, calibration, min, max, startpos)
FILE *fp;
float **epoch;
int num_samples, num_channels;
float *calibration, *min, *max;
int *startpos;
  {
  register int channel, leftmost_nonzero_sample;
  for(channel = 0; channel != num_channels; channel++)
    {
    fseek(fp, 5L, SEEK_CUR);	/*seek past channel header*/
    epoch[channel] = (float *)calloc(num_samples, sizeof(**epoch));
    leftmost_nonzero_sample = read_channel(fp, epoch[channel], num_samples, calibration[channel], min+channel, max+channel);
    if(leftmost_nonzero_sample > *startpos)
      *startpos = leftmost_nonzero_sample;
    }
  }

/*Find the integrals over the intervals between zero-crossings of
  (ch1%offset-ch0)[0..baseline_length-offset-1], and insert them into binary
  tree 't'.  Produce at most 'max_num_integrals' integrals, where
  max_num_integrals > 0.*/
int find_integrals(ch0, ch1, baseline_length, offset, t, max_num_integrals)
float *ch0, *ch1;
int baseline_length, offset;
btree *t;
int max_num_integrals;
  {
  register int sample, num_integrals;
  register float cur_diff, prev_diff, integral;
  cur_diff = ch1[offset]-*ch0;
  sample = 0;
  do
  /*inv: (ch1%offset-ch0)[0..sample-1] contains no zero-crossings, and
    cur_diff==ch1[offset+sample]-ch0[sample]*/
    {
    sample++;
    prev_diff = cur_diff;
    cur_diff = ch1[offset+sample]-ch0[sample];
    }
  while((cur_diff*prev_diff > 0.0) && (sample != baseline_length-offset-1));
  num_integrals = 0;
/*if a zero-crossing was found, begin computing integrals between crossings*/
  if(sample != baseline_length-offset-1)
    {
    integral = 0.0;
    do
    /*inv: 'integral' is the discrete integral from the point just past the
      most recent zero-crossing to sample-1, and
      cur_diff==((ch1+offset)-(ch0))[sample].*/
      {
      integral += cur_diff;
      sample++;
      prev_diff = cur_diff;
      cur_diff = ch1[offset+sample]-ch0[sample];
      }
    while((sample != baseline_length-offset-1) && (cur_diff*prev_diff > 0.0));
    while((sample != baseline_length-offset-1) && (num_integrals != max_num_integrals))
      {
    /*inv: 'integral' is the most recent integral computed, 't' contains all
      the other integrals computed thus far, 'num_integrals' is the size of 't',
      and cur_diff==ch1[offset+sample]-ch0[sample].*/
      insert_newest(integral, t);
      num_integrals++;
      integral = 0.0;
	do
	{
	integral += cur_diff;
	sample++;
	prev_diff = cur_diff;
	cur_diff = ch1[offset+sample]-ch0[sample];
	}
      while((sample != baseline_length-offset-1) && (cur_diff*prev_diff > 0.0));
      }
    }
  return(num_integrals);
  }

/*Sample the distribution of integrals between zero-crossings of the difference
  between ch1[0..baseline_length-1] and ch0[0..baseline_length-1], with a
  randomly selected offset.  Place the sampled values, in nondecreasing order,
  in distrib[0..num_distrib_samples-1].*/
void sample_difference_distribution(ch0, ch1, baseline_length, distrib, num_distrib_samples)
float *ch0, *ch1;
int baseline_length;
float *distrib;
int num_distrib_samples;
  {
  register btree *t;
  register int distrib_size;
  t = create_btree(num_distrib_samples);
  distrib_size = find_integrals(ch0, ch1, baseline_length, 0, t, num_distrib_samples);
  if(distrib_size == 0)
    {
  /*there are no pairs of zero-crossings in the entire interval; sum it and
    return a constant distribution*/
    int sample;
    float sum;
    sum = 0.0;
    for(sample = 0; sample != baseline_length; sample++)
      sum += ch1[sample]-ch0[sample];
    if(sum == 0.0)
      for(sample = 0; sample != num_distrib_samples; sample++)
	distrib[sample] = 0.0;
    else if(sum < 0.0)
      {
      for(sample = 0; sample != num_distrib_samples/2; sample++)
	distrib[sample] = -1e+19;
      while(sample != num_distrib_samples)
	distrib[sample++] = sum;
      }
    else /*sum > 0.0*/
      {
      for(sample = 0; sample != num_distrib_samples/2; sample++)
	distrib[sample] = sum;
      while(sample != num_distrib_samples)
	distrib[sample++] = 1e+19;
      }
    }
  else
    {
  /*slide ch1 a random distance along ch0*/
    while(distrib_size != num_distrib_samples)
      distrib_size += find_integrals(ch0, ch1, baseline_length, (int)(random()%(baseline_length-3)), t, num_distrib_samples-distrib_size);
    traverse_btree(t, distrib);
    }
  destroy_btree(t);
  }

void main(argc, argv)
int argc;
char **argv;
  {
  register int channel, file, i;
  int num_channels, srate;
  double min_range, round_step;
  float *(epoch[FI_MAX][CH_MAX]);
  FILE *avg;
  char erp[packed_sizeof_SETUP], ch_hdr[CH_MAX][packed_sizeof_ELECTLOC];
  int startpos, endpos;	/*endpoints of interval to be used as baseline*/
  float calibration[FI_MAX][CH_MAX], yscale[CH_MAX], ymin[CH_MAX], ymax[CH_MAX], difference_distribution[FI_MAX/2][CH_MAX][NUM_DIFFERENCE_SAMPLES];

  fprintf(stderr, "Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  fontsize = 12;
  title_height = fontsize;
  lmargin = 72*0.5;
  rmargin = 72*0.5;
  tmargin = 72*0.5;
  bmargin = 72*0.5;
  page_height = 72*11.0;
  page_width = 72*8.5;
  plot_leading = 12.0;
  tick_length = 5.0;
  image_width = page_width-lmargin-rmargin;

  if((argc < 6) || (argc%2 != 0))
    {
    fprintf(stderr, "usage: %s <label> <min. range> <rounding step> <file0.avg> <file0'.avg> ... <fileN.avg> <fileN'.avg>\n", *argv);
    exit(1);
    }
  sscanf(argv[2], "%lf", &min_range);
  sscanf(argv[3], "%lf", &round_step);
  if((min_range <= 0.0) || (round_step <= 0.0))
    {
    fprintf(stderr, "crazy parameter values (%f, %f)\n", min_range, round_step);
    exit(1);
    }
  argv += 4;
  argc -= 4;
  if(argc > FI_MAX)
    {
    fprintf(stderr, "recompile with FI_MAX=%d\n", argc);
    exit(1);
    }
  if((avg = fopen(*argv,
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(*argv);
    exit(errno);
    }
/*read # of channels, # of samples, and real time bounds*/
  fread(erp, packed_sizeof_SETUP, 1, avg);
  num_samples = S_pnts(erp);
  num_channels = S_nchannels(erp);
  srate = S_rate(erp);
  if(num_channels > CH_MAX)
    {
    fprintf(stderr, "Recompile %s with CH_MAX=%d\n", *argv, num_channels);
    fclose(avg);
    exit(1);
    }
  plot_height = (page_height-tmargin-bmargin)/num_channels;
  xmin = S_xmin(erp);
  xmax = S_xmax(erp);
/*read channel labels and calibrations*/
  fseek(avg, (long)packed_sizeof_SETUP, SEEK_SET);
  for(channel = 0; channel != num_channels; channel++)
    {
    fread(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, avg);
    calibration[0][channel] = EL_calib(ch_hdr[channel]);
    ymin[channel] = 1e+19;
    ymax[channel] = -1e+19;
    }
  startpos = 0;
  endpos = (int)(num_samples*xmin/(xmin-xmax));
/*load channel data*/
  fseek(avg, (long)(packed_sizeof_SETUP+num_channels*packed_sizeof_ELECTLOC), SEEK_SET);
  read_channels(avg, epoch[0], num_samples, num_channels, calibration[0], ymin, ymax, &startpos);
  fclose(avg);
  for(file = 1; file < argc; file++)
    {
    if((avg = fopen(argv[file],
#ifdef MSDOS
      "rb"
#else
      "r"
#endif
      )) == NULL)
      {
      perror(argv[file]);
      exit(errno);
      }
    fread(erp, packed_sizeof_SETUP, 1, avg);
    if((num_samples != S_pnts(erp)) || (num_channels != S_nchannels(erp)) || (xmin != S_xmin(erp)) || (xmax != S_xmax(erp)))
      {
      fclose(avg);
      fprintf(stderr, "%s (%d points, %d channels, time bounds %f..%f) doesn't match %s (%d points, %d channels, time bounds %f..%f)\n", argv[file], S_pnts(erp), S_nchannels(erp), S_xmin(erp), S_xmax(erp), *argv, num_samples, num_channels, xmin, xmax);
      exit(1);
      }
    for(channel = 0; channel != num_channels; channel++)
      {
      int c;
      i = 0;
      c = getc(avg);
      while((i != EL_lab_length-1) && (c != '\0') && (c == ch_hdr[channel][EL_lab_offset+i]))
	{
	i++;
	c = getc(avg);
	}
      if((i != EL_lab_length-1) && (c != ch_hdr[channel][EL_lab_offset+i]))
	{
	fclose(avg);
	fprintf(stderr, "Channel label %d in file %s doesn't match channel label %d (%s) in file %s\n", channel, argv[file], channel, EL_lab(ch_hdr[channel]), *argv);
	exit(1);
	}
      fread(ch_hdr[channel]+i+1, packed_sizeof_ELECTLOC-i-1, 1, avg);
      calibration[file][channel] = EL_calib(ch_hdr[channel]);
      }
  /*load channel data*/
    fseek(avg, (long)(packed_sizeof_SETUP+num_channels*packed_sizeof_ELECTLOC), SEEK_SET);
    read_channels(avg, epoch[file], num_samples, num_channels, calibration[file], ymin, ymax, &startpos);
    fclose(avg);
    }
/*round ranges to nearest 10.0uV, and set scaling factors*/
  for(channel = 0; channel != num_channels; channel++)
    if(ymin[channel] != ymax[channel])
      {
      ymin[channel] = (float)(floor(ymin[channel]/round_step)*round_step);
      ymax[channel] = (float)(ceil(ymax[channel]/round_step)*round_step);
      if(ymax[channel]-ymin[channel] < min_range)
	ymax[channel] = ymin[channel]+min_range;
      yscale[channel] = (plot_height-plot_leading)/(ymax[channel]-ymin[channel]);
      }
    else
      yscale[channel] = 0.0;
/*clip boundaries of baseline interval, which may be contaminated by
  neighbouring signals*/
  if(endpos-startpos > 2*((ST_MSECS_PER_TICK*srate+500)/1000))
    {
    startpos += (ST_MSECS_PER_TICK*srate+500)/1000;
    endpos -= (ST_MSECS_PER_TICK*srate+500)/1000;
    }
/*set font and print title banner*/
  printf("%%!PS\n/Helvetica findfont\n%d scalefont\nsetfont\n", fontsize);
  printf("%f %f moveto\n", lmargin, page_height-tmargin);
  printf("(%s) show\n", argv[-3]);
  for(channel = 0; channel != num_channels; channel++)
    draw_axes(channel, EL_lab(ch_hdr[channel]), ymin[channel], ymax[channel]);
/*sample baseline distributions of channel differences*/
  for(file = 0; file < argc; file += 2)
    for(channel = 0; channel != num_channels; channel++)
      sample_difference_distribution(epoch[file][channel]+startpos+(endpos-startpos)/2, epoch[1+file][channel]+startpos+(endpos-startpos)/2, endpos-startpos-(endpos-startpos)/2, difference_distribution[file/2][channel], NUM_DIFFERENCE_SAMPLES);
  for(file = 0; file != argc; file++)
    {
  /*set a new line type for this file*/
    if(file > 0)
      printf("[%d 1] 0 setdash\n", argc-file);
  /*plot the curves*/
    for(channel = 0; channel != num_channels; channel++)
      {
      printf("newpath\n%f %f moveto\n", lmargin, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0-ymax[channel]*yscale[channel]);
    /*send file names and channel name to stderr*/
      if(file%2 == 0)
	{
	int i;
	fprintf(stderr, "%s %s ", argv[file], argv[1+file]);
	for(i = 0; (ch_hdr[channel][EL_lab_offset+i] != '\0') && (i != EL_lab_length); i++)
	  putc(ch_hdr[channel][EL_lab_offset+i], stderr);
	putc('\n', stderr);
	}
      plotpoints(epoch[file][channel], epoch[file&~1][channel], epoch[file|1][channel], num_samples, (int)(num_samples*xmin/(xmin-xmax)), difference_distribution[file/2][channel], NUM_DIFFERENCE_SAMPLES, yscale[channel]);
      free(epoch[file][channel]);
      free(epoch[1+file][channel]);
      }
    }
  printf("showpage\n");
  exit(0);
  }
