chaos_analyser.c

混沌分析工具

    xv_set(LYA_EMBEDDING,PANEL_VALUE,0,NULL);
    xv_set(SVD_M_SLIDER,PANEL_VALUE,M,NULL);
    xv_set(M_SLIDER,PANEL_VALUE,M,NULL);
    xv_set(M_SLIDER,PANEL_LABEL_STRING,"Time Delay",NULL);
  
  }
  else{
    EMBED='s';
    xv_set(SVD_M_SLIDER,PANEL_LABEL_STRING,"SVD window",NULL);
    xv_set(LYA_EMBEDDING,PANEL_VALUE,1,NULL);
    xv_set(SVD_M_SLIDER,PANEL_VALUE,SVD_COLS,NULL);
    xv_set(M_SLIDER,PANEL_VALUE,SVD_COLS,NULL);
    xv_set(M_SLIDER,PANEL_LABEL_STRING,"SVD",NULL);
  }
}



/*************************************************************************

                 Multi Column input
		 
**************************************************************************/

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

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


  ioptr = fopen(filename,"r");

  if (ioptr == (FILE *)NULL) {
    printf("\nCannot open file !");
    exit(1);
  }
 
  for (i=1;i<=num_of_points;i++){
    for (l=0;l< window;l++){
      if (read_double(ioptr,&databit)==1)
	printf("\n.. break in the data !!!");
      dataarray[i][l+1]=databit;
    }
  }

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


/**********************************************************************

                 File input procedures

***********************************************************************/

void readfile(filename,window,num_of_points,m,MAX_VAL,data_array,skip)
     char   *filename;
      int    window,*num_of_points,m,skip;
     double *MAX_VAL,**data_array;

{
  double databit,buffer[100];
  FILE *ioptr;
  int i,l,j,check,sect=1,leave=0;
  char data[256],space;

  ioptr = fopen(filename,"r");
 
  if (ioptr == (FILE *)NULL) {
     xv_set(mes_frame,XV_SHOW,TRUE,NULL); 
     leave=1;
  }
 
  if (leave==0){
    i=0;
    while (read_double(ioptr,&databit)==0 && i<=skip)
      i++;
    for(l=1;l<=m* window;l++){
      read_double(ioptr,&databit);
      buffer[l]=databit;
    }
    i=0;
    while (read_double(ioptr,&databit)==0 && i<=*num_of_points){
      i++;
      for (l=0;l< window;l++)
	data_array[i][l+1]=buffer[l* m +1];
      for (l=1;l<=m* window -1;l++)
	buffer[l]=buffer[l+1];
      buffer[m* window]=databit;
    if (databit<0) databit*=-1;
      if (*MAX_VAL < databit)
	*MAX_VAL=databit;
    }
    *num_of_points=i-1;
    fclose(ioptr);
  }
}


void
newfile()
{
  int num_points=NO_POINTS,skip1,skip2,i,j;
  char filein[50];
  Xv_Window window1=canvas_paint_window(canvas);

  MAX_VAL=0;
  
  COUNTER=1;

  sprintf(filein,"%s",FILENAME);
  
  if (EMBED=='s'){
     svdred(DATA,GLOBAL_DIM,filein,&num_points,SVD_COLS);
     for (i=1;i<=num_points;i++)
       for (j=1;j<=GLOBAL_DIM;j++){
	 if (fabs(DATA[i][j])>MAX_VAL)
	   MAX_VAL=fabs(DATA[i][j]);
       }
   }
  else
     readfile(filein,GLOBAL_DIM,&num_points,M,&MAX_VAL,DATA,0);
		  
  
  NO_POINTS=num_points;
  xv_set(NUM_POINTS,PANEL_VALUE,num_points,NULL);
 
  SIZE=2.5* MULT/MAX_VAL;
  TSIZE=45/MAX_VAL;

  redraw_proc((Canvas)NULL,window1,
		dpy,
	       canvas_win, (Xv_xrectlist *) NULL);
}

void
get_name(item,event)
    Panel_item item;
     Event *event;
{
  strcpy(FILENAME,(char *)xv_get(item,PANEL_VALUE));
   
  newfile();
}


/*********************************************************************

               Error messages

*********************************************************************/

void
message_off(item,value)
     Panel_item item;
     int value;
{
   xv_set(mes_frame,XV_SHOW,FALSE,NULL);
}

/*********************************************************************

               Print messages

*********************************************************************/
void
get_print_name(item,event)
    Panel_item item;
     Event *event;
{
  strcpy(PRINT_NAME,(char *)xv_get(item,PANEL_VALUE));
  xv_set(PRINT_MESSAGE,PANEL_LABEL_STRING,"Printing ....",NULL);
  if (PRINT=='l'){
    save_to_file(PRINT_NAME,canvas_lya,canvas_win_lya);
  }
  if (PRINT=='m'){
    print_attr(PRINT_NAME);
  }
  if (PRINT=='s'){
    print_lsvd(PRINT_NAME);
  }
   if (PRINT=='i'){
    print_mut(PRINT_NAME);
  }
  xv_set(PRINT_MESSAGE,PANEL_LABEL_STRING,NULL,NULL);  
  xv_set(print_frame,XV_SHOW,FALSE,NULL);
}

/********************************************************************

            Time Domain

********************************************************************/


void
tdomain_on_off(item,value)
     Panel_item item;
     int value;
{
  Xv_Window window2=canvas_paint_window(canvas_time); 

  if (value==1)
    xv_set(time_frame,XV_SHOW,FALSE,NULL);
  else
     xv_set(time_frame,XV_SHOW,TRUE,NULL);

  time_paint_proc((Canvas)NULL,window2,
        dpy,
        canvas_win_time, (Xv_xrectlist *) NULL);  
}

void
zoom_out(item,value)
     Panel_item item;
     int value;
{

   Xv_Window window2=canvas_paint_window(canvas_time);

   TMULT*=0.8;
  time_paint_proc((Canvas)NULL,window2,
        dpy,
        canvas_win_time, (Xv_xrectlist *) NULL); 
}

void
zoom_in(item,value)
     Panel_item item;
     int value;
{
  Xv_Window window2=canvas_paint_window(canvas_time);

  TMULT*=1.25;
 time_paint_proc((Canvas)NULL,window2,
        dpy,
        canvas_win_time, (Xv_xrectlist *) NULL); 
}

void
time_paint_proc(canvas_time,window, dpy, xwin, xrects)
Canvas        canvas_time;           /* Ignored */
Xv_Window    window;
Display      *dpy;
Window        xwin;
Xv_xrectlist *xrects;           /* Ignored */

{
 char msg[50]; 
 GC gc2;
 int i,j,start,length;
 
 gc2 = DefaultGC(dpy,DefaultScreen(dpy));

 length=250/TMULT - LENGTH/2;
 XClearWindow(dpy, xwin);
 XSetForeground(dpy, gc, pixel_table[BLACK]);
 sprintf(msg,"%3.0lf",MAX_VAL);
 XDrawString(dpy,xwin,gc,5,10,msg,strlen(msg));
 sprintf(msg,"%3.0lf",-MAX_VAL);
 XDrawString(dpy,xwin,gc,1,95,msg,strlen(msg));
 sprintf(msg,"0",-MAX_VAL);
 XDrawString(dpy,xwin,gc,5,50,msg,strlen(msg));
 XDrawLine(dpy,xwin,gc2,50,5,50,95);
 XDrawLine(dpy,xwin,gc2,55,5,45,5);
 XDrawLine(dpy,xwin,gc2,55,95,45,95);
 XDrawLine(dpy,xwin,gc2,55,50,45,50);
 j=0;
 if ((start=COUNTER-length)<0) start=0;
 XSetForeground(dpy, gc2, pixel_table[GREY]); 
 
 for (i=start;i<=COUNTER;i++){
   j++;
   XDrawLine(dpy,xwin,gc2,j*TMULT+50,50-TSIZE*DATA[i][1],(j+1)*TMULT+50,50-TSIZE*DATA[i+1][1]);
 }

  XSetForeground(dpy, gc2, pixel_table[BLACK]);
 for (i=COUNTER;i<=COUNTER+LENGTH;i++){
   j++;
   XDrawLine(dpy,xwin,gc2,j*TMULT+50,50-TSIZE*DATA[i][1],(j+1)*TMULT+50,50-TSIZE*DATA[i+1][1]);
 } 
 XSetForeground(dpy, gc2, pixel_table[GREY]); 
 
 for (i=COUNTER+LENGTH;i<=COUNTER+LENGTH+length;i++){
   j++;
   XDrawLine(dpy,xwin,gc2,j*TMULT+50,50-TSIZE*DATA[i][1],(j+1)*TMULT+50,50-TSIZE*DATA[i+1][1]);
 }
}



/********************************************************************

            Poincare section procedures

********************************************************************/

void
adjust_p_size(item,value)
     Panel_item item;
     int value;
{
  Xv_Window window=canvas_paint_window(canvas);

  MULT=value;
  P_SIZE=1.5* value/MAX_VAL;
  redraw_proc((Canvas)NULL,window,
        dpy,
        canvas_win, (Xv_xrectlist *) NULL);
}

void
poin_on_off(item,value)
     Panel_item item;
     int value;
{
  Xv_Window window=canvas_paint_window(canvas_poin);
  if (value==1)
    xv_set(poin_frame,XV_SHOW,FALSE,NULL);
  else{

    xv_set(poin_frame,XV_SHOW,TRUE,NULL);
   }
  poin_paint_proc((Canvas)NULL,window,
        dpy,
        canvas_win_poin, (Xv_xrectlist *) NULL);
}

void
poin_paint_proc(canvas_poin,window, dpy, xwin, xrects)
Canvas        canvas_poin;           /* Ignored */
Xv_Window    window;
Display      *dpy;
Window        xwin;
Xv_xrectlist *xrects;           /* Ignored */

{
  int i,start;
  double x,y,x_last,y_last;
  char pos_msg[128];
  
  XClearWindow(dpy, xwin);
  XSetForeground(dpy, gc, pixel_table[BLACK]); 

  P_SIZE=1.5* MULT/MAX_VAL; 
  for (i=2;i<NO_POINTS;i++){
    if (T1*DATA[i-1][1]+T2*DATA[i-1][2]+T3*DATA[i-1][3] <(250-X_SET)/SIZE && T1*DATA[i][1]+T2*DATA[i][2]+T3*DATA[i][3]>(250-X_SET)/SIZE){
      x=T4*DATA[i][1]+T5*DATA[i][2]+T6*DATA[i][3];
      y=T7*DATA[i][1]+T8*DATA[i][2]+T9*DATA[i][3];
      XDrawPoint(dpy,xwin,gc,X_SET_P+x*P_SIZE,Y_SET_P-y*P_SIZE);
    }
  }
}

/********************************************************************

          Mutual Information

*********************************************************************/
void
mut_print(item,value)
     Panel_item item;
     int        value;
{
  xv_set(print_frame,XV_SHOW,TRUE,NULL);
  
  PRINT='i';
}

void
mut_on_off(item,value)
     Panel_item item;
     int value;
{
  Xv_Window window=canvas_paint_window(canvas_mut);
  if (value==1)
    xv_set(mut_frame,XV_SHOW,FALSE,NULL);
  else{
    xv_set(mut_frame,XV_SHOW,TRUE,NULL);
   }
  mut_paint_proc((Canvas)NULL,window,
        dpy,
        canvas_win_mut, (Xv_xrectlist *) NULL);
}

void read_mutfile(filename,num_of_points,data,max,max_delay,radius)
char *filename;
double max,radius;
int *num_of_points,*data,max_delay;

{
  double databit,buffer,how_close,temp;
  FILE *ioptr;
  int i=0,l,j=1,m=0;


  ioptr = fopen(filename,"r");
  if (ioptr == (FILE *)NULL)
    {
      printf("\nCannot open file !");
      exit(1);
    }

  for (i=1;i<COUNTER-*num_of_points/2;i++){
    read_double(ioptr,&databit);
  }
  while (read_double(ioptr,&databit)!=2 && m< *num_of_points + max_delay){
    m++;
    temp=MAX_VAL + databit;
    if (temp<0){
      printf("\nDiameter too small ... increase by at least %lf",temp);
      exit(1);
    }
    data[m]= (int)(temp/radius);
  }
  *num_of_points=m-max_delay-1;
  
  fclose(ioptr);
}


mut_info()
{
  char *file_in;
  double radius,max,Hqs,Hqsj,Iqs,Hq;
  float percentage;
  int num_of_points,*data,num,*s,**sq,delay,max_delay,i,j,n;
   
  file_in = FILENAME;
  num_of_points=MUT_NO_POINTS;
  radius=MUT_RAD;
  max_delay=MUT_DELAY;
  max=MAX_VAL*2;

  data=ivector(1,num_of_points+max_delay);
  sq=imatrix(0,(int)(max/radius),0,(int)(max/radius));
  s=ivector(0,(int)(max/radius));

  read_mutfile(file_in,&num_of_points,data,max,max_delay,radius);
  num=num_of_points;
  
  for (i=0;i<(int)(max/radius);i++)
    s[i]=0;

  for (n=1;n<=num;n++)
    s[data[n]]++;

  Hq=0;
  for (i=0;i<(int)(max/radius);i++){
    if (s[i]!=0)
      Hq+=(s[i]/(float)num)*( log(s[i]/(float)num)/log(2));
  }
  
  for (delay=1;delay<=max_delay;delay++){
    for (i=0;i<(int)(max/radius);i++)
      for (j=0;j<(int)(max/radius);j++)
	sq[i][j]=0;
   
    for (n=1;n<=num;n++)
      sq[data[n]][data[n+delay]]++;
    
    
    Hqs=0;
    for (i=0;i<(int)(max/radius);i++){
      Hqsj=0;
      for (j=0;j<(int)(max/radius);j++){
	if (sq[i][j]!=0 && s[i]!=0){
	  Hqsj += (sq[i][j]/(float)s[i])*(log(sq[i][j]/(float)s[i])/log(2));
	}
      }
      Hqs += Hqsj * (s[i]/(float)num);
    }

    Iqs =-Hq+Hqs;
    mut[delay]=Iqs;
    mut_delay[delay]=delay;
  }
  MUT_NUM=max_delay;
  free_imatrix(sq,0,(int)(max/radius),0,(int)(max/radius));
  free_ivector(s,0,(int)(max/radius));
  free_ivector(data,1,num_of_points+max_delay);

}

mut_num_adjust(item,event)
     Panel_item item;
     Event *event;
{
  MUT_DELAY=xv_get(item,PANEL_VALUE);
}

mut_st_rad(item,event)
     Panel_item item;
     Event      *event;
{
  char *test[50],str[50];
 
  strcpy(str,(char *)xv_get(item,PANEL_VALUE));

  MUT_RAD=strtod(str,test);
}

mut_length(item,event)
     Panel_item item;
     Event *event;
{
  MUT_NO_POINTS=xv_get(item,PANEL_VALUE);
}
  
mut_accept(item,value)
     Panel_item item;
     int value;
{
  Xv_Window window2=canvas_paint_window(canvas_mut);