recipes_svdred.c

混沌分析的C语言程序的

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

/*  This program takes a time series and produces a trajectory matrix */
/*  no_of_columns wide by pushing a no_of_columns window through it.  */
/*  The resulting matrix is no_of_rows long. Maxrows is just the      */
/*  maximum array length set aside in memory.                         */
/*  It then calculates the covariance matrix, which is a square       */
/*  matrix of size no_of_columns.                                     */
/*  The eigenvalues of this matrix are then calculated and          */
/*  normalised to the highest one. This is singular value decomposition*/

#define max_rows 51000


void trajectory(time_series,no_of_rows,no_of_columns,ifp,trajmat,no_of_points)
     double *time_series,**trajmat;
     int *no_of_rows,no_of_columns,no_of_points;
     FILE *ifp;
     
{
  int count=0;
  int row,column;
  double x0,normalise;
  
  row=-1*no_of_columns+1;                     /* Ignores non-full rows */
  while(fscanf(ifp,"%lf",&x0)==1 && count<=no_of_points+no_of_columns){
    count++;
    time_series[row+no_of_columns]=x0;
    for (column=1;column<=no_of_columns;++column)
      if ((row+column)>0) trajmat[row+column][no_of_columns-column+1]=x0;
		  
    ++row;
  }
  *no_of_rows=row;

  /* Normalise */
  normalise = pow((double)*no_of_rows,-0.5);
  for (row=1;row<=*no_of_rows;++row)
    for (column=1;column<=no_of_columns;++column)
      trajmat[row][column]=trajmat[row][column]*normalise;
      
}

void covariance(time_series,no_of_rows,no_of_columns,covmat)
     double *time_series,**covmat;
     int no_of_rows,no_of_columns;
{
  int row,column,i;
  double x0;
      
  /* Covariance matrix calculation  */
  for (row=1;row<=no_of_columns;++row)
    for (column=1;column<=no_of_columns;++column){
      x0=0.0;
      for(i=0;i<no_of_rows;++i)
	x0+=time_series[i+row]*time_series[i+column];
      x0=x0/(double)no_of_rows;
      covmat[row][column]=x0;
    }
}




void eignormalise(d,n)
double d[];
int n;
{
      int i;
      double highest=0.0;
      for (i=1;i<=n;++i) d[i]=sqrt(d[i]);    /* square root for svd*/
      for (i=1;i<=n;++i) highest+=d[i];
      for (i=1;i<=n;++i) d[i]=log(d[i]/highest);
}

void reduced_trajectory(x,c,N,n,d,r)
double **x,**c,**r;
int N,n,d;
{
  int i,j,k;
  double entry;
  
  for (i=1;i<=N;++i){
    for(j=1;j<=d;++j){
      entry=0.0;
      for(k=1;k<=n;++k)
	entry+=x[i][k]*pow((double)N,0.5)*c[k][j];
      /* lose normalisation factor 1/(root N)*/
      r[i][j]=entry;
    }
  }
  
}

svdred(reduced,reduced_columns,file_in,no_of_points,no_of_columns)
     double **reduced;
     int reduced_columns,*no_of_points,no_of_columns;
     char file_in[100];
{
  FILE *ifp;
  double *time_series;
  double **trajmat;
  double **covmat;
  double *eigvals;
  double **eigvecs;


  int column,row,i,no;
  double normalise;
  int no_of_rows;
 
 
  if(!(ifp = fopen (file_in,"r"))) {
    printf("Couldn't open in file : %s\n",file_in);
    exit(1);
  }
  
  time_series=dvector(1,max_rows);
  trajmat = dmatrix(1,max_rows,1,no_of_columns);
  covmat = dmatrix(1,no_of_columns,1,no_of_columns);
  eigvecs = dmatrix(1,no_of_columns,1,no_of_columns);
  eigvals = dvector(1,no_of_columns);
  
  trajectory(time_series,&no_of_rows,no_of_columns,ifp,trajmat,*no_of_points);
  covariance(time_series,no_of_rows,no_of_columns,covmat);

  svdcmp(covmat,no_of_columns,no_of_columns,eigvals,eigvecs);

  
  reduced_trajectory(trajmat,eigvecs,no_of_rows,no_of_columns,reduced_columns,reduced);
  if (*no_of_points>no_of_rows-no_of_columns)
    *no_of_points=no_of_rows-no_of_columns;
  free_dvector(time_series,1,max_rows);
  free_dmatrix(trajmat,1,max_rows,1,no_of_columns); 
  free_dmatrix(covmat,1,no_of_columns,1,no_of_columns);
  free_dmatrix(eigvecs,1,no_of_columns,1,no_of_columns);
  free_dvector(eigvals,1,no_of_columns);
}