/*zplotavg.c - translate into PostScript code a pair consisting of a GnuroScan
  averaged z-scores and a GnuroScan probabilities file.
  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 <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"

#define CH_MAX 16

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;

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("(%.0f) show\ngsave\n/Symbol findfont\n%d scalefont\nsetfont\n(s) show\ngrestore\n", ymin, fontsize);
  printf("%f %f moveto\n", y_axis_elevation, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0);
  printf("(%.0f) show\ngsave\n/Symbol findfont\n%d scalefont\nsetfont\n(s) show\ngrestore\n", ymax, fontsize);
  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");
  }

void plotpoints(fp, pfile, n, scale, bonferroni_factor)
FILE *fp, *pfile;
int n;
float scale, bonferroni_factor;
  {
  register int sample;
  register double cur_y, prev_y;
  register float cur_p, prev_p, grey_level, prev_grey_level;
  register int significance_flag;
  significance_flag = 1;
  printf("1 setlinewidth\n");
  cur_y = read_Intel_float(fp);
  cur_p = read_Intel_float(pfile)*bonferroni_factor;
  grey_level = 0.0;
  printf("0 %f rmoveto\n", scale*cur_y);
  for(sample = 1; sample != n; sample++)
    {
    prev_y = cur_y;
    cur_y = read_Intel_float(fp);
    prev_p = cur_p;
    cur_p = read_Intel_float(pfile);
    prev_grey_level = grey_level;
    if(significance_flag && !((prev_p <= 0.05) && (cur_p <= 0.05)))
      {
      significance_flag = 0;
      grey_level = 0.5;
      }
    else if((prev_p <= 0.05) && (cur_p <= 0.05))
      {
      significance_flag = 1;
      grey_level = prev_p+cur_p; /*0.0 .. 0.1*/
      }
    if(grey_level != prev_grey_level)
      printf("currentpoint\nstroke\nmoveto\n%f setgray\n", grey_level);
    printf("%f %f rlineto\n", image_width/(n-1), (cur_y-prev_y)*scale);
    }
  printf("stroke\n0 setgray\n");
  }

void find_range(fp, n, minp, maxp, calib)
FILE *fp;
int n;
float *minp, *maxp, calib;
  {
  register int sample;
  register double cur_y;
  for(sample = 0; sample != n; sample++)
    {
    cur_y = read_Intel_float(fp);
    if(cur_y*calib < *minp)
      *minp = cur_y*calib;
    else if(cur_y*calib > *maxp)
      *maxp = cur_y*calib;
    }
  }

void main(argc, argv)
int argc;
char **argv;
  {
  register int channel;
  int num_channels;
  double min_range, round_step;
  FILE *avg, *pfile;
  char erp[packed_sizeof_SETUP], ch_hdr[CH_MAX][packed_sizeof_ELECTLOC];
  float zcalibration[CH_MAX], pcalibration[CH_MAX], yscale[CH_MAX], ymin[CH_MAX], ymax[CH_MAX];

  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 < 5)
    {
    fprintf(stderr, "usage: %s <label> <min. range> <rounding step> <zscores.avg> <probabilities.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);
    }
  if((avg = fopen(argv[4],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[4]);
    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);
  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*/
  for(channel = 0; channel != num_channels; channel++)
    {
    fread(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, avg);
    zcalibration[channel] = EL_calib(ch_hdr[channel]);
    }
/*find ranges of all channels*/
  fseek(avg, (long)(packed_sizeof_SETUP+num_channels*packed_sizeof_ELECTLOC), SEEK_SET);
  for(channel = 0; channel != num_channels; channel++)
    {
    fseek(avg, 5L, SEEK_CUR);
    ymin[channel] = 1e+38;
    ymax[channel] = -1e+38;
    find_range(avg, num_samples, ymin+channel, ymax+channel, zcalibration[channel]);
    }
/*round ranges, 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;
/*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[1]);
  for(channel = 0; channel != num_channels; channel++)
    draw_axes(channel, EL_lab(ch_hdr[channel]), ymin[channel], ymax[channel]);

  /*plot the curve*/
  fseek(avg, (long)(packed_sizeof_SETUP+num_channels*packed_sizeof_ELECTLOC), SEEK_SET);
  if((pfile = fopen(argv[5],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[5]);
    fclose(avg);
    exit(errno);
    }
  for(channel = 0; channel != num_channels; channel++)
    {
    fseek(avg, (long)(packed_sizeof_SETUP+channel*packed_sizeof_ELECTLOC+EL_calib_offset), SEEK_SET);
    zcalibration[channel] = read_Intel_float(avg);
    fseek(pfile, (long)(packed_sizeof_SETUP+channel*packed_sizeof_ELECTLOC+EL_calib_offset), SEEK_SET);
    pcalibration[channel] = read_Intel_float(pfile);
    }
  for(channel = 0; channel != num_channels; channel++)
    {
    fseek(avg, 5L, SEEK_CUR);
    fseek(pfile, 5L, SEEK_CUR);
    printf("newpath\n%f %f moveto\n", lmargin, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0-ymax[channel]*yscale[channel]);
    plotpoints(avg, pfile, num_samples, yscale[channel]*zcalibration[channel], pcalibration[channel]);
    }
  fclose(pfile);
  fclose(avg);
  printf("showpage\n");
  exit(0);
  }
