/* file: sampen.c	Doug Lake	2 August 2002
			Last revised:	6 January 2004 (by george@mit.edu)
-------------------------------------------------------------------------------
sampen: calculate Sample Entropy
Copyright (C) 2002 Doug Lake

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 also view the agreement
at http://www.fsf.org/copyleft/gpl.html.

You may contact the author via electronic mail (dlake@virginia.edu).  For
updates to this software, please visit PhysioNet (http://www.physionet.org/).

_______________________________________________________________________________

Compile this program using any standard C compiler, linking with the standard C
math library.  For example, if your compiler is gcc, use:
    gcc -o sampen -O sampen.c -lm

For brief instructions, use the '-h' option:
    sampen -h

Additional information is available at:
    http://www.physionet.org/physiotools/sampen/.

*/

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

float *readdata(char *filenm, int *filelen);
void sampen(float *y, int mm, float r, int n);
void sampen2(float *y, int mm, float r, int n);
void normalize(float *data, int n);
void help(void);

int main(int argc, char* argv[])
{
  float *data = NULL;
  float r = .2;
  int n, j, mm = 2;
  char *filenm=NULL;
  int nflag = 0;
  int vflag = 0;
  int i;

  for (i = 1; i < argc; i++)
    /* Read input and optional arguments. */
    {
      if(argv[i][0] == '-')
	{
	  switch(argv[i][1])
	    {
	    case 'h':
		help();
		exit(0);
		break;
	    case 'r':
	      if(argv[i][2] != '\0')
		{
		  fprintf(stderr, "sampen:  Unrecognized option '%s'\n", argv[i]);
		  help();
		  exit(1);
		}
	      if(i+1 < argc)
		{
		  i++;
		  r = atof(argv[i]);
		  break;
		}
	      fprintf(stderr, "sampen:  No value specified for option '-t'\n");
	      help();
	      exit(1);
	      break;
	    case 'm':
	      if(argv[i][2] != '\0')
                {
                  fprintf(stderr, "sampen:  Unrecognized option '%s'\n", argv[i]);
		  help();
                  exit(1);
                }
	      if(i+1 < argc)
		{
		  i++;
		  mm = atoi(argv[i]);
		  break;
		}
	      fprintf(stderr, "sampen:  No value specified for option '-t'\n");
	      help();
	      exit(1);
	      break;
	    case 'n':
	      if(argv[i][2] != '\0')
                {
                  fprintf(stderr, "sampen:  Unrecognized option '%s'\n", argv[i]);
		  help();
                  exit(1);
                }
	      nflag = 1;
	      break;
	    case 'v':
	      if(argv[i][2] != '\0')
         	{
		  fprintf(stderr, "sampen:  Unrecognized option '%s'\n", argv[i]);
		  help();
		  exit(1);
		}
	      vflag = 1;
	      break;
	    default:
	      fprintf(stderr, "sampen: Unrecognized option '%s'\n", argv[i]);
	      help();
	      exit(1);
	    }
	}
      else
	{
	  if(!filenm)
	    filenm = argv[i];
	  else	{ fprintf(stderr, "sampen:  Too many inputs\n"); exit(1); }
	}
    }
  if(!filenm)
	filenm = "-";	/* use the standard input if no input file specified */
  data=readdata(filenm, &n);
  if(mm > n/2)
  	{fprintf(stderr, "sampen:  m too large for time series of length %d\n", n); exit(1);}

  if(nflag)
    normalize(data, n);
  if(vflag)
    sampen2(data,mm,r,n);
  else
    sampen(data, mm, r, n);

  free(data);
  return 0;
}

float *readdata(char *filenm, int *filelen)
{
    FILE *ifile;
    long maxdat = 0L, npts = 0L;
    float *data = NULL, y, yp = 0.0;

    if (strcmp(filenm, "-") == 0) {
	filenm = "standard input";
	ifile = stdin;
    }
    else if ((ifile = fopen(filenm, "rt")) == NULL) {
      fprintf(stderr, "sampen:  Could not open %s \n", filenm);
      exit(1);
    }
	
    while (fscanf(ifile, "%f", &y) == 1) {
        if (++npts >= maxdat) {
	    float *s;

	    maxdat += 5000;	/* allow the input buffer to grow (the
				   increment is arbitrary) */
	    if ((s = realloc(data, maxdat * sizeof(float))) == NULL) {
		fprintf(stderr,
		   "sampen: insufficient memory, truncating input at row %d\n",
			npts);
	        break;
	    }
	    data = s;
	}
	data[npts] = y;
    }

    fclose(ifile);

    if (npts < 1) {
	fprintf(stderr, "sampen: %s contains no data\n", filenm);
	help();
	exit(1);
    }

    *filelen = npts;
    return (data);
}

void normalize(float *data, int n)
     /* subtracts the mean from data, then divides data by
	its variance. */
{
  int i;
  float mean = 0;
  float var = 0;
  for (i = 0; i < n; i++)
    mean += data[i];
  mean = mean/n;
  for (i = 0; i < n; i++)
    data[i] = data[i]-mean;
  for (i = 0; i <n; i++)
    var += data[i]*data[i];
  var = var/n;
  for (i = 0; i < n; i++)
    data[i] = data[i]/var;
}

void sampen2(float *y, int mm, float r, int n)
     /* calculates an estimate of sample entropy
	and the variance of the estimate. */
{
  float *p=NULL;
  float *v1=NULL, *v2=NULL, *s1=NULL, dv;
  int *R1=NULL, *R2=NULL,*F2=NULL, *F1=NULL, *F=NULL, FF;
  int *run=NULL, *run1=NULL;
  float  *A=NULL, *B=NULL;
  float *K=NULL, *n1=NULL, *n2=NULL;
  int MM;
  int m,m1,i,j,nj,jj,d,d2,i1,i2,dd;
  int  nm1, nm2, nm3, nm4;
  float y1;

  MM=2*mm;

  if ((run=(int*)calloc(n, sizeof(int)))==NULL)	exit(1);
  if ((run1=(int*)calloc(n, sizeof(int)))==NULL)	exit(1);
  if ((R1=(int*)calloc(n*MM, sizeof(int)))==NULL)	exit(1);
  if ((R2=(int*)calloc(n*MM, sizeof(int)))==NULL)	exit(1);
  if ((F=(int*)calloc(n*MM, sizeof(int)))==NULL)	exit(1);
  if ((F1=(int*)calloc(n*mm, sizeof(int)))==NULL)	exit(1);
  if ((F2=(int*)calloc(n*mm, sizeof(int)))==NULL)	exit(1);
  if ((K=(float*)calloc((mm+1)*mm, sizeof(float)))==NULL) exit(1);
  if ((A=(float*)calloc(mm, sizeof(float)))==NULL)	exit(1);
  if ((B=(float*)calloc(mm, sizeof(float)))==NULL)	exit(1);
  if ((p=(float*)calloc(mm, sizeof(float)))==NULL)	exit(1);
  if ((v1=(float*)calloc(mm, sizeof(float)))==NULL)	exit(1);
  if ((v2=(float*)calloc(mm, sizeof(float)))==NULL)	exit(1);
  if ((s1=(float*)calloc(mm, sizeof(float)))==NULL)	exit(1);
  if ((n1=(float*)calloc(mm, sizeof(float)))==NULL)	exit(1);
  if ((n2=(float*)calloc(mm, sizeof(float)))==NULL)	exit(1);

  for ( i=0; i<n-1; i++)
    {
      nj = n-i-1;
      y1=y[i];
      for ( jj=0; jj<nj; jj++)
	{
	  j=jj+i+1;
	  if (((y[j]-y1)<r)&& ((y1-y[j])<r))
	    {
	      run[jj]=run1[jj]+1;
	      m1= (mm<run[jj]) ? mm: run[jj];
	      for ( m=0; m<m1; m++)
		{
		  A[m]++;
		  if (j<n-1) B[m]++;
		  F1[i+m*n]++;   F[i+n*m]++;   F[j+n*m]++;
		}
	    }
	  else  run[jj] = 0;
	} /* for jj*/

      for (j=0; j<MM; j++)
	{
	  run1[j]=run[j];  R1[i+n*j]=run[j];

	}
      if (nj>MM-1)
	for (j=MM; j<nj; j++)   run1[j]=run[j];
    } /* for i*/

  for (i=1; i<MM; i++)
    for ( j=0; j<i-1; j++)
      R2[i+n*j] = R1[i-j-1+n*j];
  for (i=MM; i<n; i++)
    for (j=0; j<MM; j++)
      R2[i+n*j]=R1[i-j-1+n*j];
  for (i=0; i<n; i++)
    for (m=0; m<mm; m++)
      {
	FF = F[i+n*m];
	F2[i+n*m]=FF-F1[i+n*m];
	K[(mm+1)*m]+=FF*(FF-1);
      }

  for (m=mm-1; m > 0; m--)  B[m]=B[m-1];
  B[0]= (float)n*(n-1)/2;
  for (m=0; m<mm; m++)
    {
      p[m]=(float)A[m]/B[m];
      v2[m] = p[m]*(1-p[m])/B[m];
    }
  dd=1;
  for (m=0; m<mm; m++)
    {
      d2 = m+1<mm-1? m+1: mm-1;
      for (d=0; d<d2+1; d++)
	{
	  for (i1=d+1; i1<n; i1++)
	    {
	      i2 =i1-d-1;
	      nm1 = F1[i1+n*m];   nm3 = F1[i2+n*m];
	      nm2 = F2[i1+n*m];   nm4 = F2[i2+n*m];
	      for (j=0; j<(dd-1); j++)
		{
		  if (R1[i1+n*j] >=m+1)    nm1--;
		  if (R2[i1+n*j] >=m+1)    nm4--;
		}
	      for (j=0; j<2*(d+1); j++)
		if (R2[i1+n*j]>=m+1)   nm2--;
	      for (j=0; j<(2*d+1); j++)
		if (R1[i2+n*j] >=m+1)  nm3--;
	      K[d+1+(mm+1)*m]+=(float)2*(nm1+nm2)*(nm3+nm4);
	    }
	}
    }

  n1[0] = (float)n*(n-1)*(n-2);
  for (m=0; m<mm-1; m++)
    for (j=0; j<m+2; j++)
      n1[m+1]+=K[j+(mm+1)*m];
  for (m=0; m<mm; m++)
    {
      for (j=0; j<m+1; j++)  n2[m]+=K[j+(mm+1)*m];
    }

  for (m=0; m<mm; m++)
    {
      v1[m] =v2[m];
      dv = (n2[m]-n1[m]*p[m]*p[m])/(B[m]*B[m]);
      if (dv>0)   v1[m]+=dv;
      s1[m]=(float)sqrt((double)(v1[m]));
    }

  for (m=0; m<mm; m++)
    printf("SampEn(%d,%g,%d) = %f (standard deviation = %f)\n",
	   m+1, r, n, -log(p[m]), s1[m]);

  free(A);
  free(B);
  free(p);
  free(run);
  free(run1);
  free(s1);
  free(K);
  free(n1);
  free(R1);
  free(R2);
  free(v1);
  free(v2);
  free(F);
  free(F1);
  free(F2);

}

void sampen(float *y, int M, float r, int n)
     /* calculates an estimate of sample entropy but does NOT
   caclulate the variance of the estimate */
{
  float *p=NULL;
  float *e=NULL;
  long *run=NULL, *lastrun=NULL, N;
  float *A=NULL, *B=NULL;
  int M1, j, nj, jj, m;
  int i;
  float y1;

  if ((run=(long*)calloc(n, sizeof(long)))==NULL)	exit(1);
  if ((lastrun=(long*)calloc(n, sizeof(long)))==NULL)	exit(1);
  if ((A=(float*)calloc(M, sizeof(float)))==NULL)	exit(1);
  if ((B=(float*)calloc(M, sizeof(float)))==NULL)	exit(1);
  if ((p=(float*)calloc(M, sizeof(float)))==NULL)	exit(1);

  /*start running*/
  for ( i=0; i<n-1; i++)
    {
      nj = n-i-1;
      y1=y[i];
      for ( jj=0; jj<nj; jj++)
	{
	  j=jj+i+1;
	  if (((y[j]-y1)<r)&&((y1-y[j])<r))
	    {
	      run[jj]=lastrun[jj]+1;
	      M1= M<run[jj]? M: run[jj];
	      for ( m=0; m<M1; m++)
		{
		  A[m]++;
		  if (j<n-1) B[m]++;
		}
	    }
	  else  run[jj] = 0;
	} /* for jj*/
      for (j=0; j<nj; j++)  lastrun[j] = run[j];
    } /* for i*/

  N = (long)(n*(n-1)/2);
  p[0]=A[0]/N;

  printf("SampEn(1,%g,%d) = %f\n", r, n, -log(p[0]));
  for (m=1; m<M; m++) {
      p[m] = A[m] / B[m-1];
      printf("SampEn(%d,%g,%d) = %f\n", m+1, r, n, -log(p[m]));
  }

  free(A);
  free(B);
  free(p);
  free(run);
  free(lastrun);
}

static char *help_strings[] = {
 "usage: %s [OPTIONS ...] [TEXT-FILE]\n",
 "where OPTIONS may include:",
 " -h          print this usage summary",
 " -m M        set the maximum epoch length to M (default: 2)",
 " -n          normalize such that the mean of the input is 0 and the sample",
 "             variance is 1",
 " -r R        set the tolerance to R (default: 0.2)",
 " -v          output an estimate of the standard deviation of each SampEn",
 "TEXT-FILE should contain the time series (a column of decimal numbers);",
 "if omitted, sampen reads its standard input.  The output contains values of",
 "SampEn(k,r,N) where k is the epoch length, r is the tolerance, and N is the",
 "length of the input series.",
NULL
};

void help()
{
    int i;

    (void)fprintf(stderr, help_strings[0], "sampen");
    for (i = 1; help_strings[i] != NULL; i++)
	(void)fprintf(stderr, "%s\n", help_strings[i]);
}