lya_section.c

混沌分析的C语言程序的

/*

                 LYAPUNOV SPECTRUM CALCULATION
	         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This program calculates the lyapunov spectrum for time series data.
The options available are ;

 (a) Time delay embedding
 (b) Global SVD reduction
 (c) Previously embedded (or reduced) data

 (d) Diagnostics on/off (runs quicker if set to off)
Notes:

The program prompts the user for all required parameters (as above).
The output is suitable for display using Gnuplot.
The parameters are all written as comments onto the output file.


Ref : M.Banbrook, G.Ushaw, S.McLaughlin, "Lyapunov exponents from a
      time series : a noise robust algorithm", submitted to Physica D,
      1995.
*/


#include "Recipes/recipes_pseudo.c"
#include "Recipes/qr.c"
#include "Recipes/recipes_svdred.c"
#include "Recipes/n_dim_traversal.c"
#include "Recipes/city_block.c"

int    MATRIX_NUM=50;
double MATRIX_NORM;

struct linked_list{
          int                d;
	  double             distance;
	  struct linked_list *next;
	  struct linked_list *prior;
	};
typedef struct linked_list ELEMENT;
typedef ELEMENT            *LINK;

int insert_into_list(info,start,last,num,total_num,re_init)
     struct linked_list **info,**start,**last;
     int num,re_init,total_num;
{
  struct linked_list *p,*q,*old,*i,*tmp;
  int count=0,put_in_list=0;

  old=NULL;
  q= *start;
  p= *last;
  i= *info;
  while (count <num && count <= total_num ) {
    count++;
    if ((i->d >= q->d-re_init) && (i->d <= q->d+re_init)){
      free(*info);
      return(0);
    }
    if (i->distance > q->distance){
      q=q->next;
    }
    else {
      if (put_in_list==0)
	tmp=q;
      put_in_list=1;
      q=q->next;
    }
  }
  if (put_in_list==0)
    tmp=q;
  if (put_in_list==1 || num<total_num){
    if (tmp->prior){
      tmp->prior->next = i;
      i->prior= tmp->prior;
      i->next = tmp;
      tmp->prior=i;
      p->prior->next = NULL;
      *last=p->prior;
      free(p);
      return(1);
    }
    else {
      i->next = tmp;
      i->prior=NULL;
      tmp->prior=i;
      *start=i;
      p->prior->next = NULL;
      *last=p->prior;
      free(p);
      return(1);
    }
  }
  else {
    free(*info);
    return(0);
  }
}

fill_R(start,last,M)
     struct linked_list **start,**last;
     int M;
{
  struct linked_list *info;
  int count=0;

  
  info =malloc(sizeof(ELEMENT));
  info->d=0;
  info->distance=0;
  info->next=NULL;
  info->prior=NULL;
  *last=info;
  *start=info;

  while (count <=M){
    info =malloc(sizeof(ELEMENT));
    info->d=0;
    info->distance=0;
    info->prior=*last;
    (*last)->next=info;
    *last=info;
    count++;
  }
}

/**************************************************************************
  Obtain_R  : locate points that are close to test point and that
              are on a similar trajectory
**************************************************************************/
    
void lya_obtain_R(X,R_start,R_end,num,time_pos,row,col,radius,k,
	      re_init,matrix)
     struct linked_list   **R_start,**R_end;
     double **X,radius; 
     int    num,*time_pos,row,col,k,re_init;
     void   **matrix;

{
  LEAFTYPE out_list;
  LEAFTYPE tmp_list;
  int i=0,j,count=0,test=0,p,q,loop=0,temp_count,x_int[10],repeat=0;
  double distance,ev_distance,start_radius,ev_ratio,max_dist,max_evdist;
  struct linked_list *info;

  (*time_pos)+=re_init;
  for (j=0;j<col;j++)
      x_int[j]=(int)((X[*time_pos][j]+MAX_VAL)/(MATRIX_NORM*MAX_VAL));
  while (count<num || repeat==0 && i<=MATRIX_NUM/2){
    if (count>=num) repeat=1;
    out_list=NULL;
    city_block(matrix,i,col,&out_list,x_int);
    while (out_list != NULL) {
      if (out_list->val != *time_pos){
    	max_evdist=0;
	max_dist=0;      
	for (j=1;j<=col;j++){
	  if ((distance=fabs(X[*time_pos][j]-
			     X[out_list->val][j])) > max_dist)
	    max_dist=distance;
	  if ((ev_distance=fabs(X[*time_pos+re_init][j]-
				X[out_list->val+re_init][j]))>max_evdist)
	  max_evdist=ev_distance;
	}
	if (((max_evdist < 2*max_dist) || (max_evdist < radius))){
	  info =malloc(sizeof(ELEMENT));
	  info->d=out_list->val;
	  info->distance=max_dist;
	  if (insert_into_list(&info,R_start,R_end,
			       count,num,re_init)==1){
	    count++;
	  }
	}
      }
      tmp_list = out_list;
      out_list=out_list->next;
      free((char *)tmp_list);
    }
    i++;
  }
}

/**************************************************************************
  Obtain_B  : Construct neighbourhood set and the evolved neighbourhood
              set from the neighbourhood matrix
**************************************************************************/

void lya_obtain_B(X,R_start,R_end,B,B_k,l,row,col,k,B_num,M_lya,centre,re_init)
     struct linked_list   *R_start,*R_end;
     double **B,**B_k,**X;
     int   l,k,row,col,B_num,re_init,M_lya;
     char  centre;
     
{
  int i,j,count,indexes[1000];
  double distance,test;
  struct linked_list *head;
  
  for (i=1;i<=B_num;i++){         /* initialise B matrices */
    for (j=1;j<=col;j++){
      B[i][j]=0;
      B_k[i][j]=0;
    }
  }
  count=0;
  head=R_start;
  while (head && count < M_lya){
    count++;
    indexes[count]=head->d;
    head = head -> next;
  }
  if (centre!='c'){
    for (i=1;i<=count/2;i++){
      for (j=1;j<=col;j++){         /* average B vectors */
	B[i%B_num + 1][j]+=X[indexes[i+(count/2)]][j]-X[indexes[i]][j];
	B_k[i%B_num +1][j]+=X[indexes[i+(count/2)]+k][j]-X[indexes[i]+k][j];
      }
    }
  }
  else
    for (i=1;i<=count;i++){
      for (j=1;j<=col;j++){         /* average B vectors */
	B[i%B_num + 1][j]+=X[indexes[i]][j]-X[l][j];
	B_k[i%B_num +1][j]+=X[indexes[i]+k][j]-X[l+k][j];
      }
    }
}



/**************************************************************************
  readredfile : reads in a previously embedded time series
                
**************************************************************************/

void lyareadredfile(filename,dataarray,window,num_of_points,traj_check)
     char    *filename;
     double  **dataarray;
     int     window,num_of_points,*traj_check;

{
  double databit;
  FILE   *ioptr;
  int    i,l;


  ioptr = fopen(filename,"r");
  if (ioptr == (FILE *)NULL) {
    printf("\nCannot open file : %s\n",filename);
    exit(1);
  }
 
  for (i=1;i<=num_of_points;i++){
    for (l=0;l< window;l++){
      if (read_double(ioptr,&databit)==1)
	traj_check[i-1]=1;
      dataarray[i][l+1]=databit;
    }
  }

  printf("\nFinished reading data ... \n");
  fclose(ioptr);
} 



/*************************************************************************
  local_reducer : reduce B matrix into local dimensional space

*************************************************************************/
void local_reducer(B,B_evolve,corr,num_row,glob_dim,dim)
     double **B,**B_evolve,**corr;
     int    dim,glob_dim,num_row;

{
  double **s;
  double *E,**tmp_B,**tmp_B_evolve;
  int   i,j,k;

  s=dmatrix(1,num_row,1,glob_dim);
  E=dvector(1,glob_dim);
  tmp_B=dmatrix(1,num_row,1,glob_dim);
  tmp_B_evolve=dmatrix(1,num_row,1,glob_dim);
  
  for (i=1;i<=num_row;i++)
    for(j=1;j<=glob_dim;j++){
      s[i][j]=B[i][j];
      tmp_B[i][j]=B[i][j];
      tmp_B_evolve[i][j]=B_evolve[i][j];
    }
  
  multiply_matrix(B,tmp_B,corr,num_row,glob_dim,glob_dim,dim);

  svdcmp(s,num_row,glob_dim,E,corr);

  multiply_matrix(B_evolve,tmp_B_evolve,corr,num_row,glob_dim,glob_dim,dim);
  
  free_dmatrix(s,1,num_row,1,glob_dim);
  free_dvector(E,1,glob_dim);
  free_dmatrix(tmp_B,1,num_row,1,glob_dim);
  free_dmatrix(tmp_B_evolve,1,num_row,1,glob_dim);   
}


/**************************************************************************
  main : main control loop for the program
                
**************************************************************************/



calc_lyapunov(X)
     double **X;

{
  char   embed=EMBED,         /* type of embedding (s/t/a)      */
         centre='c',        /* difference or constant centre  */
         local='y';         /* local tangent reduction ?      */
  double **B,           /* neighbourhood matrix           */
         **B_inv,       /* pseudo inverse of **B          */
         **B_evolve,    /* evolved neighbourhood matrix   */
         **corr,        /* Covariance matrix (lsvd calc)  */
         **T,           /* Trajectory matrix              */
         **T_tran,      /* transpose of traj matrix, T    */
         **T_Q,         /* T.Q used in qr technique       */
         **Q_qrd,       /* Q used in qr technique         */
         **R_qrd,       /* R used in qr technique         */
         **Q_last,      /* previous loop's Q matrix       */   
         radius=MAX_VAL/10.0,        /* radius of nieghbourhood        */
         in_rad,        /* annulus radius (inner)         */
         tau=TAU;           /* 1/sampling rate                */ 
  int    *traj_check,   /* check of discontinuity         */
         N=NO_POINTS,             /* Number of vectors to form X    */
         k=EVOLVE,             /* total evolve period            */
         time_pos,      /* current evolve time position   */
         evolve,        /* current evolve step            */
         B_num=B_VECT,         /* number of vectors in B         */
         M_lya=B_AVE,             /* number of vectors in R         */
         re_init=REINIT,       /* steps between reinitialisation */
         glob_dim=GLOBAL_DIM,      /* global embedding dimension     */
         dim=LOC_DIM,           /* embedding dimension            */
         svd_win=10,     /* svd window length              */
         m=5,           /* time delay embedding delay     */
         j,i,p,q,
         test,           /* general counters               */
         y_int,          /*used for fast neighbourhood algorithm*/
         x_int[10];
  FILE   *optr;
  struct linked_list *R_start,*R_end,*info;
  Xv_Window window=canvas_paint_window(canvas_lya);
  void *matrix;
  matrix = (void *)0;


  MATRIX_NUM=(int)(pow(10,(12.0/(double)glob_dim)));
  if (MATRIX_NUM>1000)
    MATRIX_NUM=1000;
  dimsize=MATRIX_NUM;
  MATRIX_NORM = 0.001 + 2/(double)MATRIX_NUM;
  
  M_lya=M_lya * B_num;
 
/************* initialise matrices, vectors and variables *****************/

  
  B=dmatrix(1,400,1,glob_dim);
  B_inv=dmatrix(1,glob_dim,1,400);
  B_evolve=dmatrix(1,400,1,glob_dim);
  Q_qrd=dmatrix(1,glob_dim,1,glob_dim);
  Q_last=dmatrix(1,glob_dim,1,glob_dim);
  R_qrd=dmatrix(1,glob_dim,1,glob_dim);
  T=dmatrix(1,glob_dim,1,glob_dim);
  T_tran=dmatrix(1,glob_dim,1,glob_dim);
  T_Q=dmatrix(1,glob_dim,1,glob_dim);
  corr=dmatrix(1,glob_dim,1,glob_dim);
  traj_check=ivector(1,N);
  R_start=NULL;
  R_end=NULL;

  fill_R(&R_start,&R_end,400); /* initialise R list */
 
  time_pos=100;                 /* start calculation at start of file */
  for (i=1;i<=dim;i++){
    LYA_EXP[0][i]=0.0;         /* initialise lya_exp */
    for (j=1;j<=dim;j++)
      if (i==j)
	Q_qrd[i][j]=1.0;       /* initialise Q_qrd matrix */
      else
	Q_qrd[i][j]=0.0;
  }

  for (i=re_init+1;i<N-re_init;i++){
    for (j=0;j<glob_dim;j++)
      x_int[j]=(int)((X[i][j]+MAX_VAL)/(MATRIX_NORM*MAX_VAL));
    enter_value(i,&matrix,glob_dim,x_int);
  }

  lya_obtain_R(X,&R_start,&R_end,
	       M_lya,&time_pos,N,glob_dim,radius,
	       k,re_init,&matrix);
  lya_obtain_B(X,R_start,R_end,B,B_evolve,time_pos,N,glob_dim,
	   re_init,B_num,M_lya,centre,re_init);
  local_reducer(B,B_evolve,corr,B_num,glob_dim,dim);
  time_pos-=re_init;
  
  /****** main loop for evolving around the attractor *******************/
  evolve=0;
  while(evolve<=k && LYA_STOP==-1){
    evolve++;
    lya_obtain_R(X,&R_start,&R_end,
		 M_lya,&time_pos,N,glob_dim,radius,k,
		 re_init,&matrix);
    lya_obtain_B(X,R_start,R_end,B,B_evolve,time_pos,N,glob_dim,
	     re_init,B_num,M_lya,centre,re_init);
    local_reducer(B,B_evolve,corr,B_num,glob_dim,dim);
    pseudo(B_inv,B,B_num,dim);
    for (q=1;q<=dim;q++){
      for (p=1;p<=dim;p++)
	T_Q[q][p]=0.0;
      for (j=1;j<=dim;j++){
	T[q][j]=0.0;
	for (i=1;i<=B_num;i++)
	  T[q][j]+=B_inv[j][i] * B_evolve[i][q];
	for (p=1;p<=dim;p++)
	  T_Q[q][p]+=T[q][j]*Q_qrd[j][p];
      }
    }
    qrd(T_Q,dim,Q_qrd,R_qrd);  
    
    for (i=1;i<=dim;i++){
      LYA_EXP[evolve][i]=LYA_EXP[evolve-1][i]+(1/(tau* re_init))* log(R_qrd[i][i]);
    }
    notify_dispatch();
    XFlush(dpy);
    EV_POS=evolve;
    lya_paint_proc((Canvas)NULL,window,
		     dpy,
		     canvas_win_lya, (Xv_xrectlist *) NULL); 
  }
}