最大lyapunov指数与相关维数c程序.txt

最大lyapunov指数与相关维数C程序

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);
}