/* file: makecv.c	G. Moody	11 July 1988
			Last revised:	21 July 1999
-------------------------------------------------------------------------------
makecv: Derive covariance matrix from pattern vectors
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 "defs.h"

FILE *ifile;
struct patternvector pattern_vector;

main(int argc, char **argv)
{
  static double c[PVDIM][PVDIM], mean[PVDIM], z[PVDIM];
  int i, j, read_pattern_vector();
  long npv = 0L;
  
  if (argc < 2) {
    fprintf(stderr, "usage: %s input-file\n", argv[0]);
    exit(1);
  }
  if ((ifile = fopen(argv[1], "r")) == NULL) {
    fprintf(stderr, "%s: can't open %s\n", argv[0], argv[1]);
    exit(2);
  }
  
  /* Pass 1: read the pattern vectors and calculate the mean. */
  while (read_pattern_vector()) {
    for (i = 0; i < PVDIM; i++)
      mean[i] += pattern_vector.v[i];
    npv++;
  }
  if (npv == 0) {
    fprintf(stderr, "%s: input file %s contains no pattern vectors\n",
	    argv[0], argv[1]);
    exit(2);
  }
  for (i = 0; i < PVDIM; i++)
    mean[i] /= npv;
  fclose(ifile);
  
  /* Pass 2: read the pattern vectors again and calculate the covariance
     matrix. */
  
  if ((ifile = fopen(argv[1], "r")) == NULL) {
    fprintf(stderr, "%s: can't reopen %s\n", argv[0], argv[1]);
    exit(2);
  }
  while (read_pattern_vector()) {
    for (i = 0; i < PVDIM; i++)
      z[i] = pattern_vector.v[i] - mean[i];
    for (i = 0; i < PVDIM; i++)
      for (j = 0; j < PVDIM; j++)
	c[i][j] += z[i]*z[j];
  }
  fclose(ifile);
  
  /* Normalize and print the covariance matrix. */
  for (i = 0; i < PVDIM; i++)
    for (j = 0; j < PVDIM; j++)
      printf("%g%c", c[i][j]/npv, (j == PVDIM-1) ? '\n' : ' ');
  exit(0);
}

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