#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 51000void 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);}