/*
  l1d2.c ... generates scaling data for calculating the largest
             Lyapunov exponent and the correlation dimension

             Copyright (c) 1999. Michael T. Rosenstein.

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 (mtr@cs.umass.edu) or postal mail
(c/o Department of Computer Science, University of Massachusetts, 140 Governors
Drive, Amherst, MA 01003-4610 USA).  For updates to this software, please visit
PhysioNet (http://www.physionet.org/).

	     reference: M.T. Rosenstein, J.J. Collins, C.J. De Luca,
                        A practical method for calculating largest
                        Lyapunov exponents from small data sets,
                        Physica D 65:117-134, 1993.

             email contact: mtr@cs.umass.edu	     
*/

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

#define IOTA			10e-15
#define MAX_LN_R		12
#define MIN_LN_R		-12
#define N_LN_R			600

typedef struct
{
  char  fileName[16];
  long  startIndex, stopIndex;
  int   seriesN, m, J, W, divergeT;
} test;
	
double	**AllocateDMatrix(long nRows, long nCols);
void	ComputeSlopes(void);
void	FreeDMatrix(double **mat, long nRows);
void	GenerateTemplateFile(void);
long	GetData(char *fileName, int tsN, long start, long stop);
void	LineFit(double *data, double n, double *m, double *b, double *rr);
void	PercentDone(int percentDone);
void	ProcessTest(int testN);
void	ReadTestFile(char *fileName);
void	SaveL1Results(char *fileRoot);
void	SaveD2Results(char *fileRoot);

int     gNTests, gMaxDivergeT;
double  *gData, *gNDivergence;
double  **gDivergence, **gCSum;
test    *gTest;

int    main(void)
{
  int   i;
  char  str[256];
	
  printf("\n");
  printf("*** L1D2: generates scaling information for L1 and D2 ***\n");
  printf("          v1.0 Copyright (c) 1999 M.T. Rosenstein\n");
  printf("                                  (mtr@cs.umass.edu)\n\n");
  printf("          reference: M.T. Rosenstein, J.J. Collins, C.J. De Luca,\n");
  printf("                     A practical method for calculating largest\n");
  printf("                     Lyapunov exponents from small data sets,\n");
  printf("                     Physica D 65:117-134, 1993.\n\n");
	
  GenerateTemplateFile();
  
  printf("Enter test file name: ");
  scanf("%s", str);
  ReadTestFile(str);
	
  printf("\nEnter output file root (no extension): ");
  scanf("%s", str);

  /* allocate the divergence and correlation sum arrays */
  for(i=0; i<gNTests; i++)
    if(gTest[i].divergeT > gMaxDivergeT)
      gMaxDivergeT = 1+gTest[i].divergeT;
  gNDivergence = (double *) calloc(gMaxDivergeT, sizeof(double));
  gDivergence = AllocateDMatrix(gMaxDivergeT, 2*gNTests);
  gCSum = AllocateDMatrix(N_LN_R, 2*gNTests);
	
  for(i=0; i<gNTests; i++)
    ProcessTest(i);

  ComputeSlopes();
  printf("\n");
  SaveL1Results(str);
  SaveD2Results(str);
	
  printf("\nSuccess!\n");
  return(0);
}

double **AllocateDMatrix(long nRows, long nCols)
{
  long   i;
  double **mat;

  /* allocate space for the row pointers */
  mat = (double **) calloc(nRows, sizeof(double *));
  if(mat == NULL)
    {
      printf("OUT OF MEMORY: AllocateDMatrix(%ld, %ld)\n\n", nRows, nCols);
      exit(1);
    };

  /* allocate space for each row pointer */
  for(i=0; i<nRows; i++)
    mat[i] = (double *) calloc(nCols, sizeof(double));
  if(mat[i-1] == NULL)
    {
      printf("OUT OF MEMORY: AllocateDMatrix(%ld, %ld)\n\n", nRows, nCols);
      exit(1);
    };

  return(mat);
}

void   ComputeSlopes(void)
{
  int    i, i2, j;
  double k, m, b, rr;
  double *data;
	
  data = (double *) calloc(N_LN_R>gMaxDivergeT ? N_LN_R : gMaxDivergeT, 
			   sizeof(double));
  if(data == NULL)
    {
      printf("OUT OF MEMORY: ComputeSlopes\n\n");
      exit(1);
    };
	
  for(i=0; i<gNTests; i++)
    {
      i2 = i+gNTests;
      
      /*** work on correlation dimension first ***/
      
      k = (double) N_LN_R/(MAX_LN_R-MIN_LN_R);
      
      /* pack the array column into the local data array */
      for(j=0; j<N_LN_R; j++)
	data[j] = gCSum[j][i];		
      /* compute the 7-point slopes */	
      for(j=3; j<N_LN_R-3; j++)
	{
	  LineFit(data+j-3, 7, &m, &b, &rr);
	  gCSum[j][i2] = k*m;
	};
      /* handle the edges */
      LineFit(data, 5, &m, &b, &rr); gCSum[2][i2] = k*m;
      LineFit(data+N_LN_R-5, 5, &m, &b, &rr); gCSum[N_LN_R-3][i2] = k*m;
      LineFit(data, 3, &m, &b, &rr); gCSum[1][i2] = k*m;
      LineFit(data+N_LN_R-3, 3, &m, &b, &rr); gCSum[N_LN_R-2][i2] = k*m;
      gCSum[0][i2] = k*(data[1]-data[0]);
      gCSum[N_LN_R-1][i2] = k*(data[N_LN_R-1]-data[N_LN_R-2]);
      
      /*** now work on divergence data ***/
		
      /* pack the array column into the local data array */
      for(j=0; j<gMaxDivergeT; j++)
	data[j] = gDivergence[j][i];		
      /* compute the 7-point slopes */	
      for(j=3; j<gMaxDivergeT-3; j++)
	{
	  LineFit(data+j-3, 7, &m, &b, &rr);
	  gDivergence[j][i2] = m;
	};
      /* handle the edges */
      LineFit(data, 5, &m, &b, &rr); gDivergence[2][i2] = m;
      LineFit(data+gMaxDivergeT-5, 5, &m, &b, &rr);
      gDivergence[gMaxDivergeT-3][i2] = m;
      LineFit(data, 3, &m, &b, &rr);
      gDivergence[1][i2] = m;
      LineFit(data+gMaxDivergeT-3, 3, &m, &b, &rr);
      gDivergence[gMaxDivergeT-2][i2] = m;
      gDivergence[0][i2] = data[1]-data[0];
      gDivergence[gMaxDivergeT-1][i2] = data[gMaxDivergeT-1]-
	data[gMaxDivergeT-2];
    };
}

void   FreeDMatrix(double **mat, long nRows)
{
  long   i;
	
  /* free space for each row pointer */
  for(i=nRows-1; i>=0; i--)
    free(mat[i]);
  
  /* free space for the row pointers */
  free(mat);
}

void   GenerateTemplateFile(void)
{
  FILE   *outFile;

  outFile = fopen("sample.l1d2", "w");

  fprintf(outFile, "* Header info starts with an asterisk.\n*\n");
  fprintf(outFile, "* Each line of the test file contains the name of a data file followed\n");
  fprintf(outFile, "* by the parameters for the test:\n");
  fprintf(outFile, "*   file_name series# startIndex stopIndex m J W divergeT\n*\n");
  fprintf(outFile, "*      file_name  = name of the data file\n");
  fprintf(outFile, "*      series#    = time series number to use for delay reconstruction\n");
  fprintf(outFile, "*      startIndex = index of first data point to read (usually 1)\n");
  fprintf(outFile, "*      stopIndex  = index of last data point to read\n");
  fprintf(outFile, "*                      (enter 0 for maximum)\n");
  fprintf(outFile, "*      m          = embedding dimension\n");
  fprintf(outFile, "*      J          = delay in samples\n");
  fprintf(outFile, "*      W          = window size for skipping temporally close nearest neighbors\n");
  fprintf(outFile, "*      divergT    = total divergence time in samples\n");
  fprintf(outFile, "*   example: lorenz.dat 1 1 0 5 7 100 300\n");
	
  fclose(outFile);
}

long   GetData(char *fileName, int tsN, long start, long stop)
{
  long   i, j, len;
  long   nHead, nCols, nRows, nPts;
  char   str[1024];
  double dummy;
  FILE   *inFile;

  /* try to open the data file */
  inFile = fopen(fileName, "r");
  if(inFile == NULL)
    {
      printf("FILE ERROR: GetData(%s)\n\n", fileName);
      exit(1);
    };

  /* skip the header */
  nHead = 0;
  fgets(str, 1024, inFile);
  while(str[0]=='*')
    {
      nHead++;
      fgets(str, 1024, inFile);
    };
		
  /* figure out the number of columns */
  len = strlen(str);
  i = 0;
  nCols = 0;
  while(i<len)
    {
      if(!isspace(str[i]))
	{
	  nCols++;
	  while(!isspace(str[i]) && i<len)
	    i++;
	  while(isspace(str[i]) && i<len)
	    i++;
	}
      else
	i++;
    };

  /* figure out the number of rows; assume there's at least one */
  nRows = 1;
  while(fgets(str, 1024, inFile) != NULL)
    nRows++;
  
  if(stop<start)
    stop = nRows;
  else if(stop>nRows)
    stop = nRows;
  if(start<1 || start>stop-3)
    start = 1;
  nPts = stop-start+1;
  gData = calloc(nPts, sizeof(double));
	
  /* now read the time series data */
  rewind(inFile);
  for(i=0; i<nHead; i++)
    fgets(str, 1024, inFile);
  
  for(i=1; i<start; i++)
    for(j=0; j<nCols; j++)
      fscanf(inFile, "%lf", &dummy);
  for(i=0; i<nPts; i++)
    {
      for(j=0; j<tsN; j++)
	fscanf(inFile, "%lf", &dummy);
      gData[i] = dummy;
      for(; j<nCols; j++)
	fscanf(inFile, "%lf", &dummy);
    };			
  fclose(inFile);

  return(nPts);
}

void   LineFit(double *data, double n, double *m, double *b, double *rr)
{
  int    i;
  double sx, sy, sxy, sx2, sy2;
  double x, y, k, mTemp, bTemp, rrTemp;

  sx = sy = sxy = sx2 = sy2 = 0;
  for(i=0; i<n; i++)
    {
      x = i;
      y = data[i];
      sx += x; sy += y;
      sx2 += x*x; sy2 += y*y;
      sxy += x*y;
    };
  k = n*sx2-sx*sx;
  mTemp = (n*sxy-sx*sy)/k;
  bTemp = (sx2*sy-sx*sxy)/k;
  k = sy*sy/n;
  if(k==sy2)
    rrTemp = 1.0;
  else
    {
      rrTemp = (bTemp*sy+mTemp*sxy-k)/(sy2-k);
      rrTemp = 1.0 - (1.0-rrTemp)*(n-1.0)/(n-2.0);
    };
  *m = mTemp;
  *b = bTemp;
  *rr = rrTemp;
}

void   PercentDone(int percentDone)
{
  static last=100;
  
  if(percentDone<last)
    {
      last = 0;
      printf("0");
      fflush(stdout);
    }
  else if(percentDone>last && percentDone%2==0)
    {
      last = percentDone;
      if(percentDone%10==0)
	printf("%d", percentDone/10);
      else
	printf(".");
      fflush(stdout);
    };
}

void   ProcessTest(int testN)
{
  long   m, J, W, divergeT, neighborIndex, maxIndex;
  long   i, j, k, T, CSumIndex, percentDone;
  long   nPts, nCompletedPairs=0, nVectors;
  char   *isNeighbor;
  double distance, d;
  double k1, k2, temp, kNorm;

  printf("\nProcessing test %d of %d:  ", testN+1, gNTests);
  fflush(stdout);
	
  m = gTest[testN].m;
  J = gTest[testN].J;
  W = gTest[testN].W;
  divergeT = gTest[testN].divergeT;
  nPts = GetData(gTest[testN].fileName, gTest[testN].seriesN, 
		 gTest[testN].startIndex, gTest[testN].stopIndex);
	
  k1 = (double) N_LN_R/(MAX_LN_R-MIN_LN_R);
  k1 *= 0.5; /* accounts for the SQUARED distances later on */
  k2 = N_LN_R/2;

  nVectors = nPts-J*(m-1);
	
  isNeighbor = (char *) calloc(nVectors, sizeof(char));
  if(isNeighbor==NULL)
    {
      printf("\nOUT OF MEMORY: ProcessTest\n\n");
      exit(1);
    };	
	
  /* clear the divergence arrays */
  for(i=0; i<gMaxDivergeT; i++)
    gNDivergence[i] = gDivergence[i][testN] = 0;

  /* loop through vectors */
  i = 0;
  while(i<nVectors)
    {
      percentDone = 100.0*nCompletedPairs/nVectors*2+0.5;
      percentDone = 100.0*i/nVectors+0.5;
      PercentDone(percentDone);
      
      if(!isNeighbor[i])
	{
	  distance = 10e10;
	  
	  /* find the nearest neighbor for the vector at i */
	  /* ignore temporally close neighbors using W */
	  if(i>W)
	    for(j=0; j<i-W; j++)
	      {
		/* calculate distance squared */
		d=0;
		for(k=0; k<m; k++)
		  {
		    T = k*J;
		    temp = gData[i+T]-gData[j+T];
		    temp *= temp;
		    d += temp;
		  };
		d += IOTA;
	
		/* map the squared distance to array position */
		CSumIndex = k1*log(d)+k2;
		if(CSumIndex<0)
		  CSumIndex = 0;
		if(CSumIndex>=N_LN_R)
		  CSumIndex = N_LN_R-1;
		
		/* increment the correlation sum array */
		gCSum[CSumIndex][testN]++;

		/* now compare to current nearest neighbor */
		/* use IOTA above to ignore nbrs that are too close */
		if(d<distance)
		  {
		    distance=d;
		    neighborIndex=j;
		  };
	      };
	  if(i<nVectors-W)
	    for(j=i+W; j<nVectors; j++)
	      {
		d=0;
		for(k=0; k<m; k++)
		  {
		    T = k*J;
		    temp = gData[i+T]-gData[j+T];
		    temp *= temp;
		    d += temp;
		  };
		d += IOTA;

		CSumIndex = k1*log(d)+k2;
		if(CSumIndex<0)
		  CSumIndex = 0;
		if(CSumIndex>=N_LN_R)
		  CSumIndex = N_LN_R-1;

		gCSum[CSumIndex][testN]++;
		
		if(d<distance)
		  {
		    distance=d;
		    neighborIndex=j;
		  };
	      };
	  isNeighbor[neighborIndex] = 1;

	  /* track divergence */
	  for(j=0; j<=divergeT; j++)
	    {
	      maxIndex = nPts-m*J-j-1;
	      if(i<maxIndex && neighborIndex<maxIndex)
		{
		  d=0;
		  for(k=0; k<m; k++)
		    {
		      T = k*J+j;
		      temp = gData[i+T]-gData[neighborIndex+T];
		      temp *= temp;
		      d += temp;
		    };
		  d += IOTA;
		  gNDivergence[j]++;
		  temp = 0.5*log(d);
		  gDivergence[j][testN] += temp;
		};
	    };
	  nCompletedPairs++;
	};
      i++;
    };

  /* integrate the correlation sum array to get the correlation sum curve */
  for(i=1; i<N_LN_R; i++)
    gCSum[i][testN] += gCSum[i-1][testN];

  /* next normalize values */
  kNorm = 1.0/gCSum[N_LN_R-1][testN];
  for(i=0; i<N_LN_R; i++)
    gCSum[i][testN] *= kNorm;

  /* now take natural log of the values */
  for(i=0; i<N_LN_R; i++)
    {
      temp = gCSum[i][testN];
      if( (temp<0.000045) || (temp>0.990050) )
	gCSum[i][testN] = 0;
      else
	gCSum[i][testN] = log(temp);
    };

  /* now take care of Lyapunovv average */
  for(i=0; i<=divergeT; i++)
    if(gNDivergence[i]>0)
      gDivergence[i][testN] /= gNDivergence[i];

  free(isNeighbor);
  free(gData);
}

void   ReadTestFile(char *fileName)
{
  int    i;
  int    nHead, nRows;
  char   str[1024];
  FILE   *inFile;

  printf("\nReading Test File...\n");

  /* try to open the data file */
  inFile = fopen(fileName, "r");
  if(inFile == NULL)
    {
      printf("FILE ERROR: ReadTestFile(%s)\n\n", fileName);
      exit(1);
    };

  /* skip the header */
  nHead = 0;
  fgets(str, 1024, inFile);
  while(str[0]=='*')
    {
      nHead++;
      fgets(str, 1024, inFile);
    };

  /* figure out the number of rows; assume there's at least one */
  nRows = 1;
  while(fgets(str, 1024, inFile) != NULL && !isspace(str[0]))
    nRows++;
  gNTests = nRows;
  	
  /* allocate the test array */
  gTest = (test *) calloc(gNTests, sizeof(test));
  if(gTest == NULL)
    {
      printf("OUT OF MEMORY: ReadTestFile(%d)\n\n", gNTests);
      exit(1);
    };
  	
  printf("detected %d %s\n", gNTests, gNTests==1 ? "test" : "tests");
	
  /* rewind the file and skip the header */
  rewind(inFile);
  for(i=0; i<nHead; i++)
    fgets(str, 1024, inFile);
  
  for(i=0; i<gNTests; i++)
    fscanf(inFile, "%s %d %ld %ld %d %d %d %d\n",
	   gTest[i].fileName, &gTest[i].seriesN, &gTest[i].startIndex,
	   &gTest[i].stopIndex, &gTest[i].m, &gTest[i].J, &gTest[i].W,
	   &gTest[i].divergeT);

  fclose(inFile);
}

void   SaveD2Results(char *fileRoot)
{
  int    i, i1, i2, testN, keepGoing;
  char   str[256];
  double k1, k2;
  FILE   *outFile;
	
  printf("\nSaving data for correlation dimension...\n");

  sprintf(str, "%s.d2", fileRoot);
  outFile = fopen(str, "w");
	
  k1 = (double) (MAX_LN_R-MIN_LN_R)/N_LN_R;
  k2 = MIN_LN_R;

  /* don't save rows of just zeros */
  keepGoing = 1;
  i1 = 0;
  while(keepGoing)
    {
      for(testN=0; testN<gNTests; testN++)
	if(gCSum[i1][testN]!=0)
	  {
	    keepGoing = 0;
	    break;
	  };
      i1 += keepGoing;
    };
  i1--;
  if(i1<0 || i1>=N_LN_R)
    i1 = 0;
  keepGoing = 1;
  i2 = N_LN_R-1;
  while(keepGoing)
    {
      for(testN=0; testN<gNTests; testN++)
	if(gCSum[i2][testN]!=0)
	  {
	    keepGoing = 0;
	    break;
	  };
      i2 -= keepGoing;
    };
  i2++;
  if(i2<0 || i2>=N_LN_R)
    i2 = N_LN_R-1;

  /* write the data */
  for(i=i1; i<i2+1; i++)
    {
      fprintf(outFile, "%lf\t", k1*i+k2);
      for(testN=0; testN<gNTests; testN++)
	fprintf(outFile, "%lf\t", gCSum[i][testN]);
		
      /* write slope data */
      fprintf(outFile, "%lf\t", k1*i+k2);
      for(; testN<2*gNTests-1; testN++)
	fprintf(outFile, "%lf\t", gCSum[i][testN]);
      fprintf(outFile, "%lf\n", gCSum[i][testN]);
    };

  fclose(outFile);
}

void    SaveL1Results(char *fileRoot)
{
  int   i, testN;
  char  str[256];
  FILE  *outFile;
	
  printf("\nSaving data for largest Lyapunov exponent...\n");

  sprintf(str, "%s.l1", fileRoot);
  outFile = fopen(str, "w");
	
  for(i=0; i<gMaxDivergeT; i++)
    {
      fprintf(outFile, "%d\t", i);
      for(testN=0; testN<gNTests; testN++)
	if(i<=gTest[testN].divergeT)
	  fprintf(outFile, "%lf\t", gDivergence[i][testN]);
	else
	  fprintf(outFile, "\t");
      
      /* write slope data */
      fprintf(outFile, "%d\t", i);
      for(; testN<2*gNTests-1; testN++)
	if(i<=gTest[testN-gNTests].divergeT)
	  fprintf(outFile, "%lf\t", gDivergence[i][testN]);
	else
	  fprintf(outFile, "\t");
      if(i<=gTest[testN-gNTests].divergeT)
	fprintf(outFile, "%lf\n", gDivergence[i][testN]);
      else
	fprintf(outFile, "\n");
    };

  fclose(outFile);
}