/* file: makefv.c	G. Moody	12 January 1992
			Last revised:	21 July 1999
-------------------------------------------------------------------------------
makefv: Derive feature vectors from pattern vectors and basis functions
Copyright (C) 1999 George B. Moody

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.

You may contact the author by e-mail (george@mit.edu) or postal mail
(MIT Room E25-505A, Cambridge, MA 02139 USA).  For updates to this software,
please visit PhysioNet (http://physionet.org/).
_______________________________________________________________________________

*/

#include <stdio.h>
#include <wfdb/wfdb.h>
#include "defs.h"

FILE *ifile;
struct patternvector pattern_vector;
struct eigenvector ev[FVDIM];
struct featurevector feature_vector;

main(int argc, char **argv)
{
  double fv_energy, pv_energy, vv;
  int i, j;
  void write_feature_vector();
  
  if (argc != 2) {
    fprintf(stderr, "usage: %s eigenvector-file <pattern-vector-file\n",
	    argv[0]);
    exit(1);
  }
  if ((ifile = fopen(argv[1], "r")) == NULL) {
    fprintf(stderr, "%s: can't read %s\n", argv[0], argv[1]);
    exit(2);
  }
  
  /* Read eigenvectors. */
  for (i = 0; i < FVDIM; i++)
    if (read_eigenvector(i) < 0) {
      fprintf(stderr, "%s: error while reading %s\n", argv[0],
	      argv[1]);
      exit(2);
    }
  fclose(ifile);
  
  /* Main processing loop. */
  while (read_pattern_vector()) {
    pv_energy = fv_energy = 0.;
    for (i = 0; i < FVDIM; i++)
      feature_vector.v[i] = 0.;
    
    /* Form dot products of pattern vector and each eigenvector. */
    for (j = 0; j < PVDIM; j++) {
      vv = pattern_vector.v[j];
      for (i = 0; i < FVDIM; i++)
	feature_vector.v[i] += vv*ev[i].e[j];
      /* Calculate squared magnitude of pattern vector. */
      pv_energy += vv * vv;
    }
    
    /* Calculate squared magnitude of feature vector. */
    for (i = 0; i < FVDIM; i++)
      fv_energy += feature_vector.v[i]*feature_vector.v[i];
    
    /* Calculate squared residual error and normalize it. */
    if (pv_energy > 0.)
      feature_vector.rsq = 1. - fv_energy/pv_energy;
    else
      feature_vector.rsq = 0.;
    
    /* Copy tag fields. */
    feature_vector.time = pattern_vector.time;
    feature_vector.anntyp = pattern_vector.anntyp;
    write_feature_vector();
  }
}

int read_eigenvector(int i)
{
  int j;
  
  if (fscanf(ifile, "%lf", &ev[i].lambda) < 1)
    return (0);
  for (j = 0; j < PVDIM; j++)
    if (fscanf(ifile, "%lf", &ev[i].e[j]) < 1)
      return (0);
  return (1);
}

int read_pattern_vector()
{
  char s[5];
  int i;
  long t;
  
  if (scanf("%ld%s", &t, s) < 2)
    return (0);
  pattern_vector.time = t;
  pattern_vector.anntyp = strann(s);
  for (i = 0; i < PVDIM; i++)
    if (scanf("%d", &pattern_vector.v[i]) < 1)
      return (0);
  return (1);
}

void write_feature_vector()
{
  int i;
  
  printf("%8ld %s", feature_vector.time, annstr(feature_vector.anntyp));
  for (i = 0; i < FVDIM; i++)
    printf(" %g", feature_vector.v[i]);
  printf(" %g\n", feature_vector.rsq);
}