/*_PROGRAM__ack_______________________________________________________________*/

/*
         1         2         3         4         5         6         7         8
12345678901234567890123456789012345678901234567890123456789012345678901234567890


DESCRIPTION:

This program calculates figures of merit for the correlation of
two waveforms as defined in the I/O Buffer Accuracy Handbook.

To compile, type the following command from the unix prompt:
cc ack.c -lm -o ack

Don't let the disclaimer scare you; it's just there to protect IBM from
get-rich-quick schemes.  I got approval from our lawyers for anyone to
use this source code.  It's "freeware," so hack away!

Program flow chart:
1) parse command line options
2) read waveform files
3) find threshold voltage crossing points
4) shift one waveform to overlay the other
5) truncate waveforms
6) interpolate both waveforms onto the same x-axis
7) compute figures of merit
8) print two interpolated waveforms to stdio


REVISION HISTORY:

06/12/00  Release 1.0
          Named in honor of Bill the Cat.
          Greg Edlund
10/23/01  Release 1.1
          Fixed bug in read_file.  I was resetting x[1] to zero for some
          reason that I don't remember.
          Variables in function linterp were indexed to 0 instead of 1.
          I was looking at junk memory.
          Added a trap for non-monotonically increasing x data, which
          prevents divide by zero errors on linterp.
          Greg Edlund
02/25/02  Release 1.2
          Matt Flora found a bug in regrid that was causing bus errors
          on solaris but not aix.  The variables xint and yint were
          misreferenced.
          Greg Edlund


LICENSE AND DISCLAIMER

This material contains IBM copyrighted sample programming source code
("Sample Code").  IBM grants you a nonexclusive license to compile, link,
execute, display, reproduce, distribute and prepare derivative works of this
Sample Code.  The Sample Code has not been thoroughly tested under all
conditions.  IBM, therefore, does not  guarantee or imply its reliability,
serviceability, or function.  IBM provides no program services for the Sample
Code.

All Sample Code contained herein is provided to you "AS IS" without any
warranties of any kind.  THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NON-INFRINGMENT ARE EXPRESSLY DISCLAIMED.  SOME
JURISDICTIONS DO NOT ALLOW THE EXCLUSION OF IMPLIED WARRANTIES, SO THE ABOVE
EXCLUSIONS MAY NOT APPLY TO YOU.  IN NO EVENT WILL IBM BE LIABLE TO ANY PARTY
FOR ANY DIRECT, INDIRECT, SPECIAL OR OTHER CONSEQUENTIAL DAMAGES FOR ANY USE OF
THE SAMPLE CODE INCLUDING, WITHOUT LIMITATION, ANY LOST PROFITS, BUSINESS
INTERRUPTION, LOSS OF PROGRAMS OR OTHER DATA ON YOUR INFORMATION HANDLING
SYSTEM OR OTHERWISE, EVEN IF WE ARE EXPRESSLY ADVISED OF THE POSSIBILITY OF
SUCH DAMAGES.

COPYRIGHT

(C) Copyright IBM CORP. 2000
All rights reserved.
US Government Users Restricted Rights -
Use, duplication or disclosure restricted
by GSA ADP Schedule Contract with IBM Corp.
Licensed Material - Property of IBM

*/

/*_IMPORT_____________________________________________________________________*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*_FUNCTIONS__________________________________________________________________*/

void usage(void);
int read_file(char*,float**,float**,double);
float find_crossing(float*,float*,double,int,int,char*);
void regrid(float*,float*,int,float*,float**,int);
void compute_fom(float*,float*,int,double,float*,float*,float*);
float *vector(int,int);
void free_vector(float*,int,int);
void polint(float[],float[],int,float,float*,float*);
float linterp(float,float*,float*,int);

/*_MACROS_____________________________________________________________________*/

/*
** Note:  WINDOW must be an even number for Lagrange interpolation to work.
*/

#define MAXBUF 120
#define TRUE 1
#define FALSE 0
#define RISE 1
#define FALL 0
#define WINDOW 4

/*_DECLARATIONS_______________________________________________________________*/


main(int argc,char **argv)
{

   char fname1[MAXBUF];      /*file name*/
   char fname2[MAXBUF];      /*file name*/
   int ndata1;               /*number of data points in file*/
   int ndata2;               /*number of data points in file*/
   float *x1, *y1;           /*waveform 1*/
   float *x2, *y2;           /*waveform 2*/
   int ndatap;               /*size of new (prime) vectors*/
   float *xp, *y1p, *y2p;    /*new vectors*/
   float xx1;                /*x value of threshold crossing 1*/
   float xx2;                /*x value of threshold crossing 2*/
   float dx;                 /*x increment*/
   float xa;                 /*minimum value of new x axis*/
   float xz;                 /*maximum value of new x axis*/
   float xu = 1e-9f;         /*units for scaling x input data*/
   float yt = 0.5;           /*y threshold used to shift waveform*/
   float yw = 1.0;           /*weighing factor used to compute average*/
   int edge = RISE;          /*rising or falling edge*/
   int i = 1;                /*data vector index*/
   int converted = 0;        /*sucessfully converted command line option count*/
   int counting = TRUE;      /*flag for counting ndatap*/
   float xcount;             /*test variable for size of new vector*/
   float fom1,fom2,fom3;     /*figures of merit*/

/*____________________________________________________________________________*/


/*
** parse command line options.
*/
   fprintf(stderr,"Accuracy Figure-of-Merit Calculator\n\n");

   for (i=1; i<argc; ++i) {
      if (strcmp(argv[i],"-h") == 0) {
         usage();
         exit(1);
      }
      if (strcmp(argv[i],"-f1") == 0) {
         ++i;
         strcpy(fname1,argv[i]);
         ++converted;
      }
      if (strcmp(argv[i],"-f2") == 0) {
         ++i;
         strcpy(fname2,argv[i]);
         ++converted;
      }
      if (strcmp(argv[i],"-e") == 0) {
         ++i;
         if ((strcmp(argv[i],"r") == 0) || (strcmp(argv[i],"f") == 0)) {
            if (strcmp(argv[i],"r") == 0) {
               edge = RISE;
            }
            if (strcmp(argv[i],"f") == 0) {
               edge = FALL;
            }
            ++converted;
         }
         else {
            printf("-e argument must be r or f.\n");
            exit(1);
         }
      }
      if (strcmp(argv[i],"-dx") == 0) {
         ++i;
         if (sscanf(argv[i],"%f",&dx) != 1) {
            printf("-dx argument must be floating point data.\n");
            exit(1);
         }
         ++converted;
      }
      if (strcmp(argv[i],"-xu") == 0) {
         ++i;
         if (sscanf(argv[i],"%f",&xu) != 1) {
            printf("-xu argument must be floating point data.\n");
            exit(1);
         }
         ++converted;
         if (xu == 0.0) {
            printf("xu = 0 will cause divide by zero error.\n");
            exit(1);
         }
      }
      if (strcmp(argv[i],"-yt") == 0) {
         ++i;
         if (sscanf(argv[i],"%f",&yt) != 1) {
            printf("-yt argument must be floating point data.\n");
            exit(1);
         }
         ++converted;
      }
      if (strcmp(argv[i],"-yw") == 0) {
         ++i;
         if (sscanf(argv[i],"%f",&yw) != 1) {
            printf("-yw argument must be floating point data.\n");
            exit(1);
         }
         ++converted;
      }
   }

   if (argc != 15 || converted != 7) {
      usage();
      exit(1);
   }

/*
** read files and return number of data points in each file.
*/
   ndata1 = read_file(fname1,&x1,&y1,xu);
   ndata2 = read_file(fname2,&x2,&y2,xu);

   fprintf(stderr,"Found %d datapoints in %s.\n",ndata1,fname1);
   fprintf(stderr,"Found %d datapoints in %s.\n",ndata2,fname2);

/*
** find point at which waveform crosses y = yth.
*/
   xx1 = find_crossing(x1,y1,yt,ndata1,edge,fname1);
   xx2 = find_crossing(x2,y2,yt,ndata2,edge,fname2);

   fprintf(stderr,"Waveform 1 crosses y = %11.4e at x = %11.4e.\n",yt,xx1);
   fprintf(stderr,"Waveform 2 crosses y = %11.4e at x = %11.4e.\n",yt,xx2);

/*
** if waveform 2 crosses right of waveform 1, shift waveform 1 right.
** if waveform 1 crosses right of waveform 2, shift waveform 2 right.
*/
   if ((xx2-xx1) > 0) {
      fprintf(stderr,"Shifting waveform 1 right by %11.4e.\n",(xx2-xx1));
      for (i=1; i<=ndata1; ++i) x1[i] = x1[i]+(xx2-xx1);
   }
   else {
      fprintf(stderr,"Shifting waveform 2 right by %11.4e.\n",(xx1-xx2));
      for (i=1; i<=ndata2; ++i) x2[i] = x2[i]+(xx1-xx2);
   }

/*
** find bounding interval of new composite waveforms.
** bounding region is defined by max x[1] and min x[ndata].
*/
   if (x1[1] > x2[1]) xa = x1[1]; else xa = x2[1];
   if (x1[ndata1] < x2[ndata2]) xz = x1[ndata1]; else xz = x2[ndata2];
   fprintf(stderr,"Bounding interval = ( %11.4e, %11.4e )\n",xa,xz);
   xcount = xa;
   ndatap = 0;
   while (counting) {
      xcount = xcount + dx;
      ++ndatap;
      if (xcount > xz) counting = FALSE;
   }
   xp  = vector(1,ndatap);

/*
** put both waveforms on the same x-axis using interpolation.
*/
   for (i=1; i<=ndatap; ++i) xp[i] = xa + (i-1)*dx;
   regrid(x1,y1,ndata1,xp,&y1p,ndatap);
   regrid(x2,y2,ndata2,xp,&y2p,ndatap);

   for (i=1; i<=ndatap; ++i) {
      xp[i] = (i-1)*dx;
      printf("%11.4e   %11.4e   %11.4e\n",xp[i],y1p[i],y2p[i]);
   }

   compute_fom(y1p,y2p,ndatap,yw,&fom1,&fom2,&fom3);
   fprintf(stderr,"\nAverage Relative Error = ");
   fprintf(stderr,"%5.2f pct.\n",fom1);
   fprintf(stderr,"Maximum Relative Error = ");
   fprintf(stderr,"%5.2f pct.\n",fom2);
   fprintf(stderr,"Maximum Absolute Error = ");
   fprintf(stderr,"%6.3f\n",fom3);

/*
** free space.
*/
   free_vector(x1,1,ndata1);
   free_vector(y1,1,ndata1);
   free_vector(x2,1,ndata2);
   free_vector(y2,1,ndata2);
   free_vector(xp,1,ndatap);
   free_vector(y1p,1,ndatap);
   free_vector(y2p,1,ndatap);

   return(0);

}

/*_FUNCTIONS__________________________________________________________________*/


void usage(void)
{
   printf("usage:  ack -f1 file1 -f2 file2 -e edge -dx dx ");
   printf("-xu xu -yt yt -yw yw\n\n");
   printf("where:  file1 = two columns of floating-point data, ");
   printf("space-delimited\n");
   printf("        file2 = two columns of floating-point data, ");
   printf("space-delimited\n");
   printf("        edge  = r for rising and f for falling\n");
   printf("        dx    = desired x increment in units of xu\n");
   printf("        xu    = units for scaling x input data\n");
   printf("        yt    = y threshold used to shift waveform\n");
   printf("        yw    = weighing factor used to compute average\n");
}

/*____________________________________________________________________________*/

int read_file(char *fname,float **px,float **py,double xu)
{

   char buf[MAXBUF];         /*character buffer*/
   FILE *fptr;               /*file stream pointer*/
   int line = 0;             /*current line number*/
   float dum1,dum2;          /*dummy for reading real data with sscanf*/
   int dataerror = FALSE;    /*data error flag*/
   float *x, *y;             /*data vectors*/
   int j,k;                  /*data vector indecies*/

/*
** open file and count lines.
*/
   fptr = fopen(fname,"r");
   if (fptr == NULL) {
      fprintf(stderr,"Cannot open file %s.\n",fname);
      exit(1);
   }
   while (fgets(buf,MAXBUF,fptr) != NULL) {
      ++line;
   }
   if (line == 0) {
      fprintf(stderr,"Empty data file %s.\n",fname);
      exit(1);
   }

/*
** allocate space for data storage.
*/
   x = vector(1,line);
   y = vector(1,line);

/*
** return pointer to beginning of file.
*/
   fseek(fptr,0L,0);
   line = 0;

/*
** read two columns of floating point data and close file.
** scale x values and reset x vector to x[1] = 0.
** check for identical x values within WINDOW.
*/
   while (fgets(buf,MAXBUF,fptr) != NULL) {
      ++line;
      if (sscanf(buf,"%f %f",&dum1,&dum2) != 2) {
         fprintf(stderr,"Data in file %s ",fname);
         fprintf(stderr,"is not in floating point format");
         fprintf(stderr,"at line %d.\n",line);
         dataerror = TRUE;
      }
      else {
         x[line] = dum1/xu;
         y[line] = dum2;
         if (line > WINDOW) {
            for (k=1;k<=WINDOW;++k) {
               for (j=k;j<=WINDOW;++j) {
                  if ((j!=k) && (x[line-j]==x[line-k])) {
                     fprintf(stderr,"Data in file %s ",fname);
                     fprintf(stderr,"is not monotonically increasing ");
                     fprintf(stderr,"between lines %d ",(line-WINDOW));
                     fprintf(stderr,"and %d.\n",line);
                     dataerror = TRUE;
                  }
               }
            }
         }
      }
      if (dataerror == TRUE) exit(1);
   }

/*
** close file.
** return pointers to new vectors and number of lines.
*/
   fclose(fptr);
   *px = x;
   *py = y;
   return(line);

}

/*____________________________________________________________________________*/


float find_crossing(float *x,float *y,double yth,int ndata,int edge,char *fname)
{

   int i;                    /*data vector index*/
   float m;                  /*slope*/
   float b;                  /*intercept*/
   float xx;                 /*x value of threshold crossing*/
   int found_xx = FALSE;     /*threshold flag*/

/*
** loop through data vectors and look for threshold crossing.
** if crossing is found, use linear interpolation to find y value.
** y = mx + b, where y is unknown.
*/
   for (i=1; i<ndata; ++i) {
      if (y[i] == yth) {
         xx = x[i];
         found_xx = TRUE;
      }
      else {
         if (edge == RISE) {
            if ((y[i]<=yth)&&(y[i+1]>yth)) {
               m = (y[i+1]-y[i])/(x[i+1]-x[i]);
               b = y[i]-(m*x[i]);
               xx = (yth-b)/m;
               found_xx = TRUE;
            }
         }
         if (edge == FALL) {
            if ((y[i]>=yth)&&(y[i+1]<yth)) {
               m = (y[i+1]-y[i])/(x[i+1]-x[i]);
               b = y[i]-(m*x[i]);
               xx = (yth-b)/m;
               found_xx = TRUE;
            }
         }
      }
   }
   if (found_xx == FALSE) {
      fprintf(stderr,"No threshold crossing found in %s.\n",fname);
      exit(1);
   }
   return(xx);
}

/*____________________________________________________________________________*/


void regrid(float *x,float *y,int ndata,float *xp,float **pyp,int ndatap)
{

   int i,j,k;                /*loop counters*/
   float ytmp;               /*temporary storage for yp*/
   float dy;                 /*uncertainty*/
   float *yp;                /*yp (y-prime) vector whose pointer is passed*/
   float xint[WINDOW + 1];   /*x interpolation vector*/
   float yint[WINDOW + 1];   /*y interpolation vector*/

   yp = vector(1,ndatap);

/*
** old waveform = x,y
** new waveform = xp,yp
** loop through new x-axis, defined in calling function.
** loop through old x-axis, looking for points that lie BETWEEN points on
** new x-axis.
** grab a window of data from old waveform and interpolate.
** you can't interpolate the whole thing at once; the interpolator will puke.
** three cases:
**   1. window bumps up against left end of x-axis
**   2. window lies safely between left and right ends of x-axis
**   3. window bumps up against right end of x-axis
*/

   for (i=1; i<=ndatap; ++i) {
      for (j=1; j<ndata; ++j) {
         if ((xp[i] >= x[j]) && (xp[i] < x[j+1])) {
            if (j <= WINDOW) {
               for (k=1; k<=WINDOW; ++k) {
                  xint[k] = x[k];
                  yint[k] = y[k];
               }
               ytmp = linterp(xp[i],xint,yint,WINDOW);
            }
            if ((j > WINDOW) && (j < (ndata-WINDOW))) {
               for (k=1; k<=WINDOW; ++k) {
                  xint[k] = x[j-(WINDOW/2)+k];
                  yint[k] = y[j-(WINDOW/2)+k];
               }
               ytmp = linterp(xp[i],xint,yint,WINDOW);
            }
            if (j >= (ndata-WINDOW)) {
               for (k=1; k<=WINDOW; ++k) {
                  xint[k] = x[(ndata-WINDOW)+k];
                  yint[k] = y[(ndata-WINDOW)+k];
               }
               ytmp = linterp(xp[i],xint,yint,WINDOW);
            }
         }
      }
      yp[i] = ytmp;
   }
   *pyp = yp;

/*
** To use the polint function from "Numerical Recipes in C:"
** polint(xint,yint,WINDOW,xp[i],&ytmp,&dy);
*/

}

/*____________________________________________________________________________*/

void compute_fom(float *y1,float *y2,int ndata,double yw,float *p_fom1,float *p_fom2,float *p_fom3)
{

   int i;                    /*data vector index*/
   float sum = 0.0;          /*sum of absolute differences*/
   float max_abs = 0.0;      /*maximum absolute error*/

   for (i=1; i<=ndata; ++i) {
      sum = sum + fabs(y1[i]-y2[i]);
      if ((fabs(y1[i]-y2[i])) > max_abs) max_abs = fabs(y1[i]-y2[i]);
   }
   sum = 100*(1-(sum/(ndata*yw)));
   *p_fom1 = sum;
   *p_fom2 = 100*(max_abs/yw);
   *p_fom3 = max_abs;

}

/*____________________________________________________________________________*/

/*
** Allocate memory for an array of floats whose index is bounded by nl and nh.
** Let's say I want an array my_floats[1,50].  I simply say:
**     my_floats = vector(1,50);
*/
float *vector(int nl,int nh)
{
   float *v;
   v = (float *)malloc((unsigned) (nh-nl+1)*sizeof(float));
   if (!v) {
      fprintf(stderr,"allocation failure in vector().\n");
      exit(1);
   }
   return v-nl;
}

/*____________________________________________________________________________*/


void free_vector(float *v,int nl,int nh)
{
   free((char*) (v+nl));
}

/*____________________________________________________________________________*/

float linterp(float xx, float *x, float *y, int n)
{
   int j,k;
   float xprod,lk,*c,pnx;

/*
** Lagrange interpolation.  See "Elementary Numerical Computing with
** Mathematica," Skeel and Keiper, 1993, McGraw-Hill, pp. 185-189.
** trap known data points before they cause a divide by zero.
*/
   for (j=1;j<=n;++j) if (xx==x[j]) return(y[j]);

   c = vector(1,n);
   xprod = 1.0;
   for (j=1;j<=n;++j) xprod=xprod*(xx-x[j]);
   for (k=1;k<=n;++k) {
      lk=1.0;
      for (j=1;j<=n;++j) if (j!=k) lk=lk*(x[k]-x[j]);
      if (lk==0.0) {
         fprintf(stderr,"Attempted divide by zero in function linterp.\n");
         exit(1);
      }
      c[k]=y[k]/lk;
   }
   pnx=0;
   for (k=1;k<=n;++k) pnx=pnx+(c[k]/(xx-x[k]));
   pnx=xprod*pnx;
   free_vector(c,1,n);
   return(pnx);
}

