/* file: hp_scatter.c
/****************************************************************
# This program replaces premature ventricular (V)-beats 
# by normal (N) beats, i.e. it interpolates normal-normal (N-N) time intervals.
# These N-N time intervals are averaged over a number of beats, 
# typically 4 to 10 beats. 
# given by the variable window.
# Then, for each V-beat we print the
#  - coupling interval: 
#    interval between V-beat and previous normal beat, if the previous beat 
#    was not a normal beat, a zero is printed
#  - V-V time interval:
#    interval between consecutive V-beats
#  - Number of Intervening normal Beats (NIB):
#    number of normal beats appearing between consecutive V-beats.
#  - corresponding interpolated and averaged N-N interval
# The output is the file file_root.scatter, wich contains 
# first the histogram of the interpolated N-N intervals in 2 columns,
# then, in five columns, coupling interval, VV interval, NIB, 
# and corresponding N-N interval
#
# The input needs the format:
# column 1: Annotation V=V-beat, N=normal, and others
# column 2: RR-interval
#
# command-line: hp_scatter file_name NN_min NN_max binsize max 
#                      > file_name.scatter
#
# where   
# binsize = 1/sampling frequency
# NN_min  = minimum of NN-interval range under consideration
# NN_max  = maximum of NN-interval range under consideration
# max     = maximum number of data points
# These values have to be taken from the data, e.g. by plotting 
# the RR-intervals.
#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# The interpolation of the V-beats by normal beats can be tricky if there 
# are many V-beats, especially many consecutive V-beats (runs of v-tach) 
# or interpolated beats (i.e. no compensatory pause). 
# Note, that we replace each V-beat by a normal beat since we are only 
# interested in the value of the time interal between the underlying normal
# beats. This might generate more normal beats than were present originally 
# if long runs of ventricular tachycardia are present in the data. 
#
# We consider the following cases:
# - isolated V-beats followed by compensatory pause (N-V-N):
#   V-beat is replaced by normal beat which is placed half way [(NV+VN)/2] 
#   between the surrounding N-beats. 
# - isolated V-beat without compensatory pause:
#   V-beat should be just taken out, while the N-N interval results 
#   from adding NV+VN. But in order to keep the number of beats unchanged, 
#   we replace both, the NV-interval and the VN-interval, by NV+VN.   
# - two consecutive V-beats blocking two normal beats
#   (= with compensatory pause)
# - two consecutive V-beats blocking one normal beat 
#   (= without compensatory pause)
# - three consecutive V-beats blocking two or three normal beats 
#   (= with compensatory pause)
# - three consecutive V-beats blocking one normal beat 
#   (= without compensatory pause)
# The identification of these cases involves thresholds that might have 
# to to be adjusted to the data at hand. Also, other cases might be 
# possible that is not dealt with presently. The result of the interpolation 
# should therefore be checked by printing out RRneu!!!  
**********************************************************/

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

#define NUM_OF_ARGC 6
#define MAX_CHAR_IN_LINE 4096
#define STRL     80

#define IA 16807 
#define IM 2147483647 
#define AM (1.0/IM) 
#define IQ 127773 
#define IR 2836 
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB) 
#define EPS 1.2e-7 
#define RNMX (1.0-EPS) 

 
main ( int argc, char *argv[] )
{
  FILE *fp, *ofp;
  double binsize, NN_min, NN_max;
  double *RR, *RRneu, *RRave, *ttime;
  double Voldbeat, x;
  long *hist_RRave;
  char *type;
  long i,l,k,jj,j, max, min, steps, window, Voldindex, NIB;
  double coupint, comp_pause, Ts, Vtime;
  long round(double );
  char dummy;

  if ( argc!=NUM_OF_ARGC )
    {
      printf( "Usage: %s RRLIST NNmin NNmax BINSIZE LENGTH\n\n",
	      argv[0]);
      printf( "argc = %d\n", argc);
      for (i = 1; i < argc; i++)
	  printf("argv[%d] = %s\n", i, argv[i]);
      printf( "(steps= (NN_max-NN_min)/binsize must be integer \n\n" );
      exit(0);
    }

  NN_min  = atof(argv[2]);  /* minimum of NN-interval range  */
  NN_max  = atof(argv[3]);  /* maximum of NN-interval range  */
  binsize = atof(argv[4]);  /*  1/sampling frequency         */
  max     = atol(argv[5]);  /* maximum number of data points */

  Voldbeat=0;
  min= 1;
  steps=round((NN_max-NN_min)/binsize);
  window=10;


  RR = (double *) malloc( (max+1)*sizeof(double) );
  for (i=0; i<=max; i++)
    RR[i]=0.0;

  RRneu = (double *) malloc( (max+1)*sizeof(double) );
  for (i=0; i<=max; i++)
    RRneu[i]=0.0;

  RRave = (double *) malloc( (max+1)*sizeof(double) );
  for (i=0; i<=max; i++)
    RRave[i]=0.;

  type = (char *) malloc( (max+1)*sizeof(char) );
  for (i=0; i<=max; i++)
    type[i]=0;

  ttime = (double *) malloc( (max+1)*sizeof(double) );
  for (i=0; i<=max; i++)
    ttime[i]=0.0;

  hist_RRave = (long *) malloc( (steps+1)*sizeof(long) );
  for (i=0; i<=steps; i++)
      hist_RRave[i]=1;



  fp = fopen( argv[1], "r" );
  i = 1;
  /*  while ((i<max)&&(fscanf(fp, "%c %lf\n", &type[i], &RR[i]) == 2) )*/
  /* this is how is used to be */
  while ((i<max)&&(fscanf(fp, "%lf %c\n", &RR[i], &type[i]) == 2) ) {
    if (type[i] == 'r' || type[i] == '!')
      type[i] = 'V';   /* count R-on-T PVCs and ventricular flutter waves as
			  V-beats */
    i++;
  }
  max=i;
  fclose( fp );  
 
  for( i=min; i<=max;i++)
    {
      ttime[i] = ttime[i-1] + RR[i];
    }
  

  
  /****************************
   * replace V-beats:
   * this part may have to be adapted to the data at hand!!!!
   * note that the number of beats is not changed
   ****************************/
  for (i=min+1;i<max;i++)
    {
      /*** if not a V-beat, just copy beat ***/
      if (type[i]!='V')                         
	{ 
	  RRneu[i]=RR[i];
	}
      /*** if isolated V-beat check if compensatory pause ***/
      if ((type[i]=='V')&&(type[i+1]!='V'))    
	{
	  x=(RR[i]+RR[i+1])/2;
	  if (x > 2*RRneu[i-1]/3)
	    { 
	      RRneu[i]=RRneu[i+1]=x;
	      i++;
	    }
	  else if (x < 2*RRneu[i-1]/3)
	    {	 
	      RRneu[i]=RRneu[i+1]=RR[i]+RR[i+1];
	      i++;
	    }
	}
      /*** if 2 consecutive V-beats check if compensatory pause ***/
      if ((type[i]=='V')&&(type[i+1]=='V')&&(type[i+2]!='V')) 
	{
	  x=(RR[i]+RR[i+1]+RR[i+2])/3;
	  if (x > 3*RRneu[i-1]/4)
	    { 
	      RRneu[i]=RRneu[i+1]=RRneu[i+2]=x;
	      i+=2;
	    }
	  else
	    {	 
	      RRneu[i]=RRneu[i+1]=RRneu[i+2]=(RR[i]+RR[i+1]+RR[i+2])/2;
	      i+=2;
	    }
	}
      /*** if 3 consecutive V-beats check if compensatory pause ***/
      if ((type[i]=='V')&&(type[i+1]=='V')&&(type[i+2]=='V')&&(type[i+3]!='V'))
	{
	  x=(RR[i]+RR[i+1]+RR[i+2]+RR[i+3])/4;
	  if (x > 3*RRneu[i-1]/4)
	    { 
	      RRneu[i]=RRneu[i+1]=RRneu[i+2]=RRneu[i+3]=x;
	      i+=3;
	    }
	  else
	    {
	      x=(RR[i]+RR[i+1]+RR[i+2]+RR[i+3])/3;
	      if (x > 3*RRneu[i-1]/4)
		{	 
		  RRneu[i]=RRneu[i+1]=RRneu[i+2]=RRneu[i+3]=x;
		  i+=3;
		}
	      else 
		{
		  x=(RR[i]+RR[i+1]+RR[i+2]+RR[i+3])/2;
		  RRneu[i]=RRneu[i+1]=RRneu[i+2]=RRneu[i+3]=x;
		  i+=3;
		}
	    }
	  }
      /*** if more then 3 consecutive v-beats or 2 v-beats 
       *** and then not a normal beat fill with interpolated 
       *** NN-interval of previous beat until normal beat appears ***/
	else if ((type[i]=='v')&&(type[i+1]=='v'))
	  {
	    while (type[i]!='n')
	      {
		RRneu[i]=RRneu[i-1];
		i++;
	      }
	    RRneu[i]=RRneu[i-1];
	  }
    }



  /*** print reconstructed normal sinus beats on file: **/
  ofp = fopen("NN_intervals", "w" );
  for (i=1;i<max;i++)
    fprintf(ofp,"%d %f %f\n", i, RRneu[i], RR[i]);



  /*** calculate average and RR_ave histogram: ***/
  for (i=min;i<min+window/2;i++)
    RRave[i]=RRneu[i];
  for (i=min+1;i<=min+window;i++)
    RRave[window/2+min]+=RRneu[i];
  for (i=window+min+1;i<=max;i++)
    RRave[i-window/2]=RRave[i-window/2-1]-RRneu[i-window]+RRneu[i];

  for (i=window+min;i<=max;i++)
    {
      RRave[i]= RRave[i]/window;
      if ((RRave[i]<=NN_max)&&(RRave[i]>=NN_min))
	{
	hist_RRave[round((RRave[i]-NN_min)/binsize)]++;
	}
    }
  for (i=0;i<=steps;i++)
    printf("%f %d\n", i*binsize+NN_min, hist_RRave[i]); 
  
  
  /* calculate and print calculate and print vv-intervals, NIB and coupling intervals:*/
  Voldindex=1;
  for (i=min;i<max;i++)
    {
      if (type[i]=='V')
	{
	  if (NIB!=0)        /* V-V interval is not a coupling interval! */
	    coupint=RR[Voldindex];
	  else
	    coupint=0;
	  comp_pause=RR[Voldindex+1];
	  NIB=i-Voldindex-1; 
	  Vtime=ttime[i]-ttime[Voldindex];
	  /* for interpolated beats you may want to subtract the NIB by one
	   * to match the NIB triplets in parasystole: */
	   	if (((coupint+comp_pause)/2 < 2*RR[Voldindex-1]/3)&&(NIB>0))
	   	    NIB=NIB-1;
	  /* you can use the averaged RR-intervals or the non-averaged RR-intervals:*/
	  /*        Ts=RRneu[Voldindex];*/
	  Ts=RRave[Voldindex];
	  if (Ts != 0)
	    {  
      	      printf("%f %f %d %f %d\n", coupint, Vtime, NIB, Ts, Voldindex);
	    }
	  Voldindex=i;
	}
    }
  
  free(RR);
  free(type);
  free(RRneu);
  free(RRave);
  free(ttime);
  free(hist_RRave);
}

int round(double x)
{
  if (ceil(x)-floor(x) < 0.5)
    return floor(x);
  else 
    return ceil(x);
}