binaries/src/disembl/Tisean_3.0.1/source_c/d2.c

   1 /*
   2  *   This file is part of TISEAN
   3  *
   4  *   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   5  *
   6  *   TISEAN is free software; you can redistribute it and/or modify
   7  *   it under the terms of the GNU General Public License as published by
   8  *   the Free Software Foundation; either version 2 of the License, or
   9  *   (at your option) any later version.
  10  *
  11  *   TISEAN is distributed in the hope that it will be useful,
  12  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  *   GNU General Public License for more details.
  15  *
  16  *   You should have received a copy of the GNU General Public License
  17  *   along with TISEAN; if not, write to the Free Software
  18  *   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  19  */
  20 /*Author: Rainer Hegger. Last modified May 10, 2000 */
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include <math.h>
  25 #include <limits.h>
  26 #include <time.h>
  27 #include "routines/tsa.h"
  28
  29 #define WID_STR "Estimates the correlation sum, -dimension and -entropy"
  30
  31 /* output is written every WHEN seconds */
  32 #define WHEN 120
  33 /* Size of the field for box assisted neighbour searching
  34    (has to be a power of 2)*/
  35 #define NMAX 256
  36 /* Size of the box for the scramble routine */
  37 #define SCBOX 4096
  38
  39 double **series;
  40 long *scr;
  41 char dimset=0,rescale_set=0,eps_min_set=0,eps_max_set=0;
  42 char *FOUT=NULL;
  43 double epsfactor,epsinv,lneps,lnfac;
  44 double EPSMAX=1.0,EPSMIN=1.e-3;
  45 double min,interval;
  46 int imax=NMAX-1,howoften1,imin;
  47 long box[NMAX][NMAX],*list,boxc1[NMAX],*listc1;
  48 unsigned long nmax;
  49 double **found,*norm;
  50 unsigned long MINDIST=0,MAXFOUND=1000;
  51 unsigned long length=ULONG_MAX,exclude=0;
  52 unsigned int DIM=1,EMBED=10,HOWOFTEN=100,DELAY=1;
  53 unsigned int verbosity=0x1;
  54 char *column=NULL;
  55 char *infile=NULL;
  56
  57 void show_options(char *progname)
  58 {
  59   what_i_do(progname,WID_STR);
  60   fprintf(stderr,"  Usage: %s [options]\n",progname);
  61   fprintf(stderr,"  Options:\n");
  62   fprintf(stderr,"Everything not being a valid option will be interpreted"
  63           " as a possible"
  64           " datafile.\nIf no datafile is given stdin is read. Just - also"
  65           " means stdin\n");
  66   fprintf(stderr,"\t-l datapoints [default is whole file]\n");
  67   fprintf(stderr,"\t-x exclude # points [default 0]\n");
  68   fprintf(stderr,"\t-d delay  [default 1]\n");
  69   fprintf(stderr,"\t-M # of components, max. embedding dim. [default 1,10]\n");
  70   fprintf(stderr,"\t-c columns [default 1,...,# of components]\n");
  71   fprintf(stderr,"\t-t theiler-window [default 0]\n");
  72   fprintf(stderr,"\t-R max-epsilon "
  73           "[default: max data interval]\n");
  74   fprintf(stderr,"\t-r min-epsilon [default: (max data interval)/1000]\n");
  75   fprintf(stderr,"\t-# #-of-epsilons [default 100]\n");
  76   fprintf(stderr,"\t-N max-#-of-pairs (0 means all) [default 1000]\n");
  77   fprintf(stderr,"\t-E use rescaled data [default: not rescaled]\n");
  78   fprintf(stderr," \t-o outfiles"
  79           " [without exts.! default datafile[.d2][.h2][.stat][.c2]]\n");
  80   fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
  81           "0='only panic messages'\n\t\t"
  82           "1='+ input/output messages'\n\t\t"
  83           "2='+ output message each time output is done\n");
  84
  85   fprintf(stderr,"\t-h show these options\n");
  86   fprintf(stderr,"\n");
  87   exit(0);
  88 }
  89
  90 void scan_options(int n,char **argv)
  91 {
  92   char *out;
  93
  94   if ((out=check_option(argv,n,'l','u')) != NULL)
  95     sscanf(out,"%lu",&length);
  96   if ((out=check_option(argv,n,'x','u')) != NULL)
  97     sscanf(out,"%lu",&exclude);
  98   if ((out=check_option(argv,n,'c','s')) != NULL)
  99     column=out;
 100   if ((out=check_option(argv,n,'d','u')) != NULL)
 101     sscanf(out,"%u",&DELAY);
 102   if ((out=check_option(argv,n,'M','2')) != NULL) {
 103     sscanf(out,"%u,%u",&DIM,&EMBED);
 104     dimset=1;
 105   }
 106   if ((out=check_option(argv,n,'t','u')) != NULL)
 107     sscanf(out,"%lu",&MINDIST);
 108   if ((out=check_option(argv,n,'R','f')) != NULL) {
 109     sscanf(out,"%lf",&EPSMAX);
 110     eps_max_set=1;
 111   }
 112   if ((out=check_option(argv,n,'r','f')) != NULL) {
 113     sscanf(out,"%lf",&EPSMIN);
 114     eps_min_set=1;
 115   }
 116   if ((out=check_option(argv,n,'#','u')) != NULL)
 117     sscanf(out,"%u",&HOWOFTEN);
 118   if ((out=check_option(argv,n,'N','u')) != NULL) {
 119     sscanf(out,"%lu",&MAXFOUND);
 120     if (MAXFOUND == 0)
 121       MAXFOUND=ULONG_MAX;
 122   }
 123   if ((out=check_option(argv,n,'E','n')) != NULL)
 124     rescale_set=1;
 125   if ((out=check_option(argv,n,'V','u')) != NULL)
 126     sscanf(out,"%u",&verbosity);
 127   if ((out=check_option(argv,n,'o','o')) != NULL)
 128     if (strlen(out) > 0)
 129       FOUT=out;
 130 }
 131
 132 void scramble(void)
 133 {
 134   long i,j,k,m;
 135   unsigned long rnd,rndf,hlength,allscr=0;
 136   long *scfound,*scnhelp,scnfound;
 137   long scbox[SCBOX],lswap,element,scbox1=SCBOX-1;
 138   double *rz,*schelp,sceps=(double)SCBOX-0.001,swap;
 139
 140   hlength=length-(EMBED-1)*DELAY;
 141
 142   if (sizeof(long) == 8) {
 143     rndf=13*13*13*13;
 144     rndf=rndf*rndf*rndf*13;
 145     rnd=0x849178L;
 146   }
 147   else {
 148     rndf=69069;
 149     rnd=0x234571L;
 150   }
 151   for (i=0;i<1000;i++)
 152     rnd=rnd*rndf+1;
 153
 154   check_alloc(rz=(double*)malloc(sizeof(double)*hlength));
 155   check_alloc(scfound=(long*)malloc(sizeof(long)*hlength));
 156   check_alloc(scnhelp=(long*)malloc(sizeof(long)*hlength));
 157   check_alloc(schelp=(double*)malloc(sizeof(double)*hlength));
 158
 159   for (i=0;i<hlength;i++)
 160     rz[i]=(double)(rnd=rnd*rndf+1)/ULONG_MAX;
 161
 162   for (i=0;i<SCBOX;i++)
 163     scbox[i]= -1;
 164   for (i=0;i<hlength;i++) {
 165     m=(int)(rz[i]*sceps)&scbox1;
 166     scfound[i]=scbox[m];
 167     scbox[m]=i;
 168   }
 169   for (i=0;i<SCBOX;i++) {
 170     scnfound=0;
 171     element=scbox[i];
 172     while(element != -1) {
 173       scnhelp[scnfound]=element;
 174       schelp[scnfound++]=rz[element];
 175       element=scfound[element];
 176     }
 177
 178     for (j=0;j<scnfound-1;j++)
 179       for (k=j+1;k<scnfound;k++)
 180         if (schelp[k] < schelp[j]) {
 181           swap=schelp[k];
 182           schelp[k]=schelp[j];
 183           schelp[j]=swap;
 184           lswap=scnhelp[k];
 185           scnhelp[k]=scnhelp[j];
 186           scnhelp[j]=lswap;
 187         }
 188     for (j=0;j<scnfound;j++)
 189       scr[allscr+j]=scnhelp[j];
 190     allscr += scnfound;
 191   }
 192
 193   free(rz);
 194   free(scfound);
 195   free(schelp);
 196 }
 197
 198 void make_c2_dim(int n)
 199 {
 200   char small;
 201   long i,j,k,x,y,i1,i2,j1,element,n1,maxi,count,hi;
 202   double *hs,max,dx;
 203
 204   check_alloc(hs=(double*)malloc(sizeof(double)*EMBED*DIM));
 205   n1=scr[n];
 206
 207   count=0;
 208   for (i1=0;i1<EMBED;i1++) {
 209     i2=i1*DELAY;
 210     for (j=0;j<DIM;j++)
 211       hs[count++]=series[j][n1+i2];
 212   }
 213
 214   x=(int)(hs[0]*epsinv)&imax;
 215   y=(int)(hs[1]*epsinv)&imax;
 216
 217   for (i1=x-1;i1<=x+1;i1++) {
 218     i2=i1&imax;
 219     for (j1=y-1;j1<=y+1;j1++) {
 220       element=box[i2][j1&imax];
 221       while (element != -1) {
 222         if (labs((long)(element-n1)) > MINDIST) {
 223           count=0;
 224           max=0.0;
 225           maxi=howoften1;
 226           small=0;
 227           for (i=0;i<EMBED;i++) {
 228             hi=i*DELAY;
 229             for (j=0;j<DIM;j++) {
 230               dx=fabs(hs[count]-series[j][element+hi]);
 231               if (dx <= EPSMAX) {
 232                 if (dx > max) {
 233                   max=dx;
 234                   if (max < EPSMIN) {
 235                     maxi=howoften1;
 236                   }
 237                   else {
 238                     maxi=(lneps-log(max))/lnfac;
 239                   }
 240                 }
 241                 if (count > 0)
 242                   for (k=imin;k<=maxi;k++)
 243                     found[count][k] += 1.0;
 244               }
 245               else {
 246                 small=1;
 247                 break;
 248               }
 249               count++;
 250             }
 251             if (small)
 252               break;
 253           }
 254         }
 255         element=list[element];
 256       }
 257     }
 258   }
 259
 260   free(hs);
 261 }
 262
 263 void make_c2_1(int n)
 264 {
 265   int i,x,i1,maxi;
 266   long element,n1;
 267   double hs,max;
 268
 269   n1=scr[n];
 270   hs=series[0][n1];
 271
 272   x=(int)(hs*epsinv)&imax;
 273
 274   for (i1=x-1;i1<=x+1;i1++) {
 275     element=boxc1[i1&imax];
 276     while (element != -1) {
 277       if (labs(element-n1) > MINDIST) {
 278         max=fabs(hs-series[0][element]);
 279         if (max <= EPSMAX) {
 280           if (max < EPSMIN)
 281             maxi=howoften1;
 282           else
 283             maxi=(lneps-log(max))/lnfac;
 284           for (i=imin;i<=maxi;i++)
 285             found[0][i] += 1.0;
 286         }
 287       }
 288       element=listc1[element];
 289     }
 290   }
 291 }
 292
 293 int main(int argc,char **argv)
 294 {
 295   char smaller,stdi=0;
 296   FILE *fout,*fstat;
 297   char *outd1,*outc1,*outh1,*outstat;
 298   int maxembed;
 299   long i1,j1,x,y,sn,n,i,j,n1,n2;
 300   long *oscr;
 301   long lnorm;
 302   double eps,*epsm,EPSMAX1,maxinterval;
 303   time_t mytime,lasttime;
 304
 305   if (scan_help(argc,argv))
 306     show_options(argv[0]);
 307
 308   scan_options(argc,argv);
 309 #ifndef OMIT_WHAT_I_DO
 310   if (verbosity&VER_INPUT)
 311     what_i_do(argv[0],WID_STR);
 312 #endif
 313
 314   infile=search_datafile(argc,argv,NULL,verbosity);
 315   if (infile == NULL)
 316     stdi=1;
 317
 318   if (FOUT == NULL) {
 319     if (!stdi) {
 320       check_alloc(FOUT=calloc(strlen(infile)+1,(size_t)1));
 321       strcpy(FOUT,infile);
 322     }
 323     else {
 324       check_alloc(FOUT=calloc((size_t)6,(size_t)1));
 325       strcpy(FOUT,"stdin");
 326     }
 327   }
 328   if (column == NULL)
 329     series=(double**)get_multi_series(infile,&length,exclude,&DIM,"",dimset,
 330                                       verbosity);
 331   else
 332     series=(double**)get_multi_series(infile,&length,exclude,&DIM,column,
 333                                       dimset,verbosity);
 334
 335   if (rescale_set) {
 336     for (i=0;i<DIM;i++)
 337       rescale_data(series[i],length,&min,&interval);
 338     maxinterval=1.0;
 339   }
 340   else {
 341     maxinterval=0.0;
 342     for (i=0;i<DIM;i++) {
 343       min=interval=series[i][0];
 344       for (j=1;j<length;j++) {
 345         if (min > series[i][j])
 346           min=series[i][j];
 347         if (interval < series[i][j])
 348           interval=series[i][j];
 349       }
 350       interval -= min;
 351       if (interval > maxinterval)
 352         maxinterval=interval;
 353     }
 354   }
 355   if (!eps_max_set)
 356     EPSMAX *= maxinterval;
 357   if (!eps_min_set)
 358     EPSMIN *= maxinterval;
 359   EPSMAX=(fabs(EPSMAX)<maxinterval) ? fabs(EPSMAX) : maxinterval;
 360   EPSMIN=(fabs(EPSMIN)<EPSMAX) ? fabs(EPSMIN) : EPSMAX/2.;
 361   EPSMAX1=EPSMAX;
 362
 363   howoften1=HOWOFTEN-1;
 364   maxembed=DIM*EMBED-1;
 365
 366   check_alloc(outd1=(char*)calloc(strlen(FOUT)+4,(size_t)1));
 367   check_alloc(outc1=(char*)calloc(strlen(FOUT)+4,(size_t)1));
 368   check_alloc(outh1=(char*)calloc(strlen(FOUT)+4,(size_t)1));
 369   check_alloc(outstat=(char*)calloc(strlen(FOUT)+6,(size_t)1));
 370   strcpy(outd1,FOUT);
 371   strcpy(outc1,FOUT);
 372   strcpy(outh1,FOUT);
 373   strcpy(outstat,FOUT);
 374   strcat(outd1,".d2");
 375   strcat(outc1,".c2");
 376   strcat(outh1,".h2");
 377   strcat(outstat,".stat");
 378   test_outfile(outd1);
 379   test_outfile(outc1);
 380   test_outfile(outh1);
 381   test_outfile(outstat);
 382
 383   check_alloc(list=(long*)malloc(length*sizeof(long)));
 384   check_alloc(listc1=(long*)malloc(length*sizeof(long)));
 385   if ((long)(length-(EMBED-1)*DELAY) <= 0) {
 386     fprintf(stderr,"Embedding dimension and delay are too large.\n"
 387             "The delay vector would be longer than the whole series."
 388             " Exiting\n");
 389     exit(VECTOR_TOO_LARGE_FOR_LENGTH);
 390   }
 391   check_alloc(scr=(long*)malloc(sizeof(long)*(length-(EMBED-1)*DELAY)));
 392   check_alloc(oscr=(long*)malloc(sizeof(long)*(length-(EMBED-1)*DELAY)));
 393   check_alloc(found=(double**)malloc(DIM*EMBED*sizeof(double*)));
 394   for (i=0;i<EMBED*DIM;i++)
 395     check_alloc(found[i]=(double*)malloc(HOWOFTEN*sizeof(double)));
 396   check_alloc(norm=(double*)malloc(HOWOFTEN*sizeof(double)));
 397   check_alloc(epsm=(double*)malloc(HOWOFTEN*sizeof(double)));
 398
 399   epsinv=1.0/EPSMAX;
 400   epsfactor=pow(EPSMAX/EPSMIN,1.0/(double)howoften1);
 401   lneps=log(EPSMAX);
 402   lnfac=log(epsfactor);
 403
 404   epsm[0]=EPSMAX;
 405   norm[0]=0.0;
 406   for (i=1;i<HOWOFTEN;i++) {
 407     norm[i]=0.0;
 408     epsm[i]=epsm[i-1]/epsfactor;
 409   }
 410   imin=0;
 411
 412   scramble();
 413   for (i=0;i<(length-(EMBED-1)*DELAY);i++)
 414     oscr[scr[i]]=i;
 415
 416   for (i=0;i<DIM*EMBED;i++)
 417     for (j=0;j<HOWOFTEN;j++)
 418       found[i][j]=0.0;
 419
 420   nmax=length-DELAY*(EMBED-1);
 421
 422   for (i1=0;i1<NMAX;i1++) {
 423     boxc1[i1]= -1;
 424     for (j1=0;j1<NMAX;j1++)
 425       box[i1][j1]= -1;
 426   }
 427   time(&lasttime);
 428   lnorm=0;
 429
 430   for (n=1;n<nmax;n++) {
 431     smaller=0;
 432     sn=scr[n-1];
 433     if (DIM > 1) {
 434       x=(long)(series[0][sn]*epsinv)&imax;
 435       y=(long)(series[1][sn]*epsinv)&imax;
 436     }
 437     else {
 438       x=(long)(series[0][sn]*epsinv)&imax;
 439       y=(long)(series[0][sn+DELAY]*epsinv)&imax;
 440     }
 441     list[sn]=box[x][y];
 442     box[x][y]=sn;
 443     listc1[sn]=boxc1[x];
 444     boxc1[x]=sn;
 445
 446     i=imin;
 447     while (found[maxembed][i] >= MAXFOUND) {
 448       smaller=1;
 449       if (++i > howoften1)
 450         break;
 451     }
 452     if (smaller) {
 453       imin=i;
 454       if (imin <= howoften1) {
 455         EPSMAX=epsm[imin];
 456         epsinv=1.0/EPSMAX;
 457         for (i1=0;i1<NMAX;i1++) {
 458           boxc1[i1]= -1;
 459           for (j1=0;j1<NMAX;j1++)
 460             box[i1][j1]= -1;
 461         }
 462         for (i1=0;i1<n;i1++) {
 463           sn=scr[i1];
 464           if (DIM > 1) {
 465             x=(long)(series[0][sn]*epsinv)&imax;
 466             y=(long)(series[1][sn]*epsinv)&imax;
 467           }
 468           else {
 469             x=(long)(series[0][sn]*epsinv)&imax;
 470             y=(long)(series[0][sn+DELAY]*epsinv)&imax;
 471           }
 472           list[sn]=box[x][y];
 473           box[x][y]=sn;
 474           listc1[sn]=boxc1[x];
 475           boxc1[x]=sn;
 476         }
 477       }
 478     }
 479
 480     if (imin <= howoften1) {
 481       lnorm=n;
 482       if (MINDIST > 0) {
 483         sn=scr[n];
 484         n1=(sn-(long)MINDIST>=0)?sn-(long)MINDIST:0;
 485         n2=(sn+MINDIST<length-(EMBED-1)*DELAY)?sn+MINDIST:
 486           length-(EMBED-1)*DELAY-1;
 487         for (i1=n1;i1<=n2;i1++)
 488           if ((oscr[i1] < n))
 489             lnorm--;
 490       }
 491
 492       if (EMBED*DIM > 1)
 493         make_c2_dim(n);
 494       make_c2_1(n);
 495       for (i=imin;i<HOWOFTEN;i++)
 496         norm[i] += (double)(lnorm);
 497     }
 498
 499     if (((time(&mytime)-lasttime) > WHEN) || (n == (nmax-1)) ||
 500         (imin > howoften1)) {
 501       time(&lasttime);
 502       fstat=fopen(outstat,"w");
 503       if (verbosity&VER_USR1)
 504         fprintf(stderr,"Opened %s for writing\n",outstat);
 505       fprintf(fstat,"Center points treated so far= %ld\n",n);
 506       fprintf(fstat,"Maximal epsilon in the moment= %e\n",epsm[imin]);
 507       fclose(fstat);
 508       fout=fopen(outc1,"w");
 509       if (verbosity&VER_USR1)
 510         fprintf(stderr,"Opened %s for writing\n",outc1);
 511       fprintf(fout,"#center= %ld\n",n);
 512       for (i=0;i<EMBED*DIM;i++) {
 513         fprintf(fout,"#dim= %ld\n",i+1);
 514         eps=EPSMAX1*epsfactor;
 515         for (j=0;j<HOWOFTEN;j++) {
 516           eps /= epsfactor;
 517           if (norm[j] > 0.0)
 518             fprintf(fout,"%e %e\n",eps,found[i][j]/norm[j]);
 519         }
 520         fprintf(fout,"\n\n");
 521       }
 522       fclose(fout);
 523       fout=fopen(outh1,"w");
 524       if (verbosity&VER_USR1)
 525         fprintf(stderr,"Opened %s for writing\n",outh1);
 526       fprintf(fout,"#center= %ld\n",n);
 527       fprintf(fout,"#dim= 1\n");
 528       eps=EPSMAX1*epsfactor;
 529       for (j=0;j<HOWOFTEN;j++) {
 530         eps /= epsfactor;
 531         if (found[0][j] > 0.0)
 532           fprintf(fout,"%e %e\n",eps,-log(found[0][j]/norm[j]));
 533       }
 534       fprintf(fout,"\n\n");
 535       for (i=1;i<DIM*EMBED;i++) {
 536         fprintf(fout,"#dim= %ld\n",i+1);
 537         eps=EPSMAX1*epsfactor;
 538         for (j=0;j<HOWOFTEN;j++) {
 539           eps /= epsfactor;
 540           if ((found[i-1][j] > 0.0) && (found[i][j] > 0.0))
 541             fprintf(fout,"%e %e\n",eps,log(found[i-1][j]/found[i][j]));
 542         }
 543         fprintf(fout,"\n\n");
 544       }
 545       fclose(fout);
 546       fout=fopen(outd1,"w");
 547       if (verbosity&VER_USR1)
 548         fprintf(stderr,"Opened %s for writing\n",outd1);
 549       fprintf(fout,"#center= %ld\n",n);
 550       for (i=0;i<DIM*EMBED;i++) {
 551         fprintf(fout,"#dim= %ld\n",i+1);
 552         eps=EPSMAX1;
 553         for (j=1;j<HOWOFTEN;j++) {
 554           eps /= epsfactor;
 555           if ((found[i][j] > 0.0) && (found[i][j-1] > 0.0))
 556             fprintf(fout,"%e %e\n",eps,log(found[i][j-1]/found[i][j]
 557                                            /norm[j-1]*norm[j])/lnfac);
 558         }
 559         fprintf(fout,"\n\n");
 560       }
 561       fclose(fout);
 562       if (imin > howoften1)
 563         exit(0);
 564     }
 565   }
 566
 567   if (infile != NULL)
 568     free(infile);
 569   free(outd1);
 570   free(outh1);
 571   free(outc1);
 572   free(outstat);
 573   free(list);
 574   free(listc1);
 575   free(scr);
 576   free(oscr);
 577   free(norm);
 578   free(epsm);
 579   for (i=0;i<EMBED*DIM;i++)
 580     free(found[i]);
 581   free(found);
 582   for (i=0;i<DIM;i++)
 583     free(series[i]);
 584   free(series);
 585
 586   return 0;
 587 }