binaries/src/disembl/Tisean_3.0.1/source_c/lyap_spec.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 Dec 4, 2005  */
  21 /*Changes:
  22   7/14/05: Changed borders of the sort routine to speed things up
  23   11/25/05: Show also absolute forecast errors
  24   12/04/05: Some more changes in sort
  25   12/20/05: Change in increase neighborhood size loop
  26   12/28/05: Found bug in memory allocation (index)
  27 */
  28
  29 #include <stdio.h>
  30 #include <stdlib.h>
  31 #include <math.h>
  32 #include <limits.h>
  33 #include <time.h>
  34 #include <string.h>
  35 #include "routines/tsa.h"
  36
  37 #define WID_STR "Estimates the spectrum of Lyapunov exponents using the\n\t\
  38 method of Sano and Sawada."
  39
  40 #define OUT 10
  41
  42 #define BOX 512
  43 #define EPSMAX 1.0
  44 #define DELAY 1
  45
  46 char epsset=0,stdo=1;
  47 char INVERSE,*outfile=NULL;
  48 char *infile=NULL;
  49 char dimset=0;
  50 char *COLUMNS=NULL;
  51 unsigned long LENGTH=ULONG_MAX,ITERATIONS,exclude=0;
  52 unsigned int EMBED=2,DIMENSION=1/*,DELAY=1*/,MINNEIGHBORS=30;
  53 unsigned int verbosity=0xff;
  54 double EPSSTEP=1.2;
  55
  56 double **series,*averr,avneig=0.0,aveps=0.0;
  57 double **mat,*vec,*abstand;
  58 double epsmin;
  59 long imax=BOX-1,count=0;
  60 long **box,*list;
  61 unsigned long *found;
  62 unsigned int alldim,**indexes;
  63
  64 void show_options(char *progname)
  65 {
  66   what_i_do(progname,WID_STR);
  67   fprintf(stderr," Usage: %s [options]\n",progname);
  68   fprintf(stderr," Options:\n");
  69   fprintf(stderr,"Everything not being a valid option will be interpreted"
  70           " as a possible"
  71           " datafile.\nIf no datafile is given stdin is read. Just - also"
  72           " means stdin\n");
  73   fprintf(stderr,"\t-l # of datapoints [default is whole file]\n");
  74   fprintf(stderr,"\t-x # of lines to be ignored [default is 0]\n");
  75   fprintf(stderr,"\t-c column to read[default 1]\n");
  76   fprintf(stderr,"\t-m # of components,embedding dimension [default %d,%d]\n",
  77           DIMENSION,EMBED);
  78   //  fprintf(stderr,"\t-d delay  [default %d]\n",DELAY);
  79   fprintf(stderr,"\t-r epsilon size to start with [default "
  80   "(data interval)/1000]\n");
  81   fprintf(stderr,"\t-f factor to increase epsilon [default: 1.2]\n");
  82   fprintf(stderr,"\t-k # of neighbors to use [default: 30]\n");
  83   fprintf(stderr,"\t-n # of iterations [default: length]\n");
  84   fprintf(stderr,"\t-I invert the time series [default: no]\n");
  85   fprintf(stderr,"\t-o name of output file [default 'datafile'.lyaps]\n");
  86   fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
  87           "0='only panic messages'\n\t\t"
  88           "1='+ input/output messages'\n");
  89   fprintf(stderr,"\t-h show these options\n");
  90   fprintf(stderr,"\n");
  91   exit(0);
  92 }
  93
  94 void scan_options(int n,char **argv)
  95 {
  96   char *out;
  97
  98   if ((out=check_option(argv,n,'l','u')) != NULL)
  99     sscanf(out,"%lu",&LENGTH);
 100   if ((out=check_option(argv,n,'x','u')) != NULL)
 101     sscanf(out,"%lu",&exclude);
 102   if ((out=check_option(argv,n,'c','s')) != NULL)
 103     COLUMNS=out;
 104   /*  if ((out=check_option(argv,n,'d','u')) != NULL)
 105       sscanf(out,"%u",&DELAY);*/
 106   if ((out=check_option(argv,n,'m','2')) != NULL) {
 107     sscanf(out,"%u,%u",&DIMENSION,&EMBED);
 108     dimset=1;
 109   }
 110   if ((out=check_option(argv,n,'n','u')) != NULL)
 111     sscanf(out,"%lu",&ITERATIONS);
 112   if ((out=check_option(argv,n,'r','f')) != NULL) {
 113     epsset=1;
 114     sscanf(out,"%lf",&epsmin);
 115   }
 116   if ((out=check_option(argv,n,'f','f')) != NULL)
 117     sscanf(out,"%lf",&EPSSTEP);
 118   if ((out=check_option(argv,n,'k','u')) != NULL)
 119     sscanf(out,"%u",&MINNEIGHBORS);
 120   if ((out=check_option(argv,n,'V','u')) != NULL)
 121     sscanf(out,"%u",&verbosity);
 122   if ((out=check_option(argv,n,'I','n')) != NULL)
 123     INVERSE=1;
 124   if ((out=check_option(argv,n,'o','o')) != NULL) {
 125     stdo=0;
 126     if (strlen(out) > 0)
 127       outfile=out;
 128   }
 129 }
 130
 131 double sort(long act,unsigned long* nfound,char *enough)
 132 {
 133   double maxeps=0.0,dx,dswap,maxdx;
 134   long self=0,i,j,del,hf,iswap,n1;
 135   unsigned long imax=*nfound;
 136
 137   *enough=0;
 138
 139   for (i=0;i<imax;i++) {
 140     hf=found[i];
 141     if (hf != act) {
 142       maxdx=fabs(series[0][act]-series[0][hf]);
 143       for (j=1;j<alldim;j++) {
 144         n1=indexes[0][j];
 145         del=indexes[1][j];
 146         dx=fabs(series[n1][act-del]-series[n1][hf-del]);
 147         if (dx > maxdx) maxdx=dx;
 148       }
 149       abstand[i]=maxdx;
 150     }
 151     else {
 152       self=i;
 153     }
 154   }
 155
 156   if (self != (imax-1)) {
 157     abstand[self]=abstand[imax-1];
 158     found[self]=found[imax-1];
 159   }
 160
 161   for (i=0;i<MINNEIGHBORS;i++) {
 162     for (j=i+1;j<imax-1;j++) {
 163       if (abstand[j]<abstand[i]) {
 164         dswap=abstand[i];
 165         abstand[i]=abstand[j];
 166         abstand[j]=dswap;
 167         iswap=found[i];
 168         found[i]=found[j];
 169         found[j]=iswap;
 170       }
 171     }
 172   }
 173
 174   if (!epsset || (abstand[MINNEIGHBORS-1] >= epsmin)) {
 175     *nfound=MINNEIGHBORS;
 176     *enough=1;
 177     maxeps=abstand[MINNEIGHBORS-1];
 178
 179     return maxeps;
 180   }
 181
 182   for (i=MINNEIGHBORS;i<imax-2;i++) {
 183     for (j=i+1;j<imax-1;j++) {
 184       if (abstand[j]<abstand[i]) {
 185         dswap=abstand[i];
 186         abstand[i]=abstand[j];
 187         abstand[j]=dswap;
 188         iswap=found[i];
 189         found[i]=found[j];
 190         found[j]=iswap;
 191       }
 192     }
 193     if (abstand[i] > epsmin) {
 194       (*nfound)=i+1;
 195       *enough=1;
 196       maxeps=abstand[i];
 197
 198       return maxeps;
 199     }
 200   }
 201
 202   maxeps=abstand[imax-2];
 203
 204   return maxeps;
 205 }
 206
 207 void make_dynamics(double **dynamics,long act)
 208 {
 209   long i,hi,j,hj,k,t=act,d;
 210   unsigned long nfound=0;
 211   double **hser,**imat;
 212   double foundeps=0.0,epsilon,hv,hv1;
 213   double new_vec;
 214   char got_enough;
 215
 216   check_alloc(hser=(double**)malloc(sizeof(double*)*DIMENSION));
 217   for (i=0;i<DIMENSION;i++)
 218     hser[i]=series[i]+act;
 219
 220   epsilon=epsmin/EPSSTEP;
 221   do {
 222     epsilon *= EPSSTEP;
 223     if (epsilon > EPSMAX)
 224       epsilon=EPSMAX;
 225     make_multi_box(series,box,list,LENGTH-DELAY,BOX,DIMENSION,EMBED,
 226                    DELAY,epsilon);
 227     nfound=find_multi_neighbors(series,box,list,hser,LENGTH-DELAY,BOX,
 228                                 DIMENSION,EMBED,DELAY,epsilon,found);
 229     if (nfound > MINNEIGHBORS) {
 230       foundeps=sort(act,&nfound,&got_enough);
 231       if (got_enough)
 232         break;
 233     }
 234   } while (epsilon < EPSMAX);
 235
 236   free(hser);
 237
 238   avneig += nfound;
 239   aveps += foundeps;
 240   if (!epsset)
 241     epsmin=aveps/count;
 242   if (nfound < MINNEIGHBORS) {
 243     fprintf(stderr,"#Not enough neighbors found. Exiting\n");
 244     exit(LYAP_SPEC_NOT_ENOUGH_NEIGHBORS);
 245   }
 246
 247   for (i=0;i<=alldim;i++) {
 248     vec[i]=0.0;
 249     for (j=0;j<=alldim;j++)
 250       mat[i][j]=0.0;
 251   }
 252
 253   for (i=0;i<nfound;i++) {
 254     act=found[i];
 255     mat[0][0] += 1.0;
 256     for (j=0;j<alldim;j++)
 257       mat[0][j+1] += series[indexes[0][j]][act-indexes[1][j]];
 258     for (j=0;j<alldim;j++) {
 259       hv1=series[indexes[0][j]][act-indexes[1][j]];
 260       hj=j+1;
 261       for (k=j;k<alldim;k++)
 262         mat[hj][k+1] += series[indexes[0][k]][act-indexes[1][k]]*hv1;
 263     }
 264   }
 265
 266   for (i=0;i<=alldim;i++)
 267     for (j=i;j<=alldim;j++)
 268       mat[j][i]=(mat[i][j]/=(double)nfound);
 269
 270   imat=invert_matrix(mat,alldim+1);
 271
 272   for (d=0;d<DIMENSION;d++) {
 273     for (i=0;i<=alldim;i++)
 274       vec[i]=0.0;
 275     for (i=0;i<nfound;i++) {
 276       act=found[i];
 277       hv=series[d][act+DELAY];
 278       vec[0] += hv;
 279       for (j=0;j<alldim;j++)
 280         vec[j+1] += hv*series[indexes[0][j]][act-indexes[1][j]];
 281     }
 282     for (i=0;i<=alldim;i++)
 283       vec[i] /= (double)nfound;
 284
 285     new_vec=0.0;
 286     for (i=0;i<=alldim;i++)
 287       new_vec += imat[0][i]*vec[i];
 288     for (i=1;i<=alldim;i++) {
 289       hi=i-1;
 290       dynamics[d][hi]=0.0;
 291       for (j=0;j<=alldim;j++)
 292         dynamics[d][hi] += imat[i][j]*vec[j];
 293     }
 294     for (i=0;i<alldim;i++)
 295       new_vec += dynamics[d][i]*series[indexes[0][i]][t-indexes[1][i]];
 296     averr[d] += (new_vec-series[d][t+DELAY])*(new_vec-series[d][t+DELAY]);
 297   }
 298
 299   for (i=0;i<=alldim;i++)
 300     free(imat[i]);
 301   free(imat);
 302 }
 303
 304 void gram_schmidt(double **delta,
 305                   double *stretch)
 306 {
 307   double **dnew,norm,*diff;
 308   long i,j,k;
 309
 310   check_alloc(diff=(double*)malloc(sizeof(double)*alldim));
 311   check_alloc(dnew=(double**)malloc(sizeof(double*)*alldim));
 312   for (i=0;i<alldim;i++)
 313     check_alloc(dnew[i]=(double*)malloc(sizeof(double)*alldim));
 314
 315   for (i=0;i<alldim;i++) {
 316     for (j=0;j<alldim;j++)
 317       diff[j]=0.0;
 318     for (j=0;j<i;j++) {
 319       norm=0.0;
 320       for (k=0;k<alldim;k++)
 321         norm += delta[i][k]*dnew[j][k];
 322       for (k=0;k<alldim;k++)
 323         diff[k] -= norm*dnew[j][k];
 324     }
 325     norm=0.0;
 326     for (j=0;j<alldim;j++)
 327       norm += sqr(delta[i][j]+diff[j]);
 328     stretch[i]=(norm=sqrt(norm));
 329     for (j=0;j<alldim;j++)
 330       dnew[i][j]=(delta[i][j]+diff[j])/norm;
 331   }
 332   for (i=0;i<alldim;i++)
 333     for (j=0;j<alldim;j++)
 334       delta[i][j]=dnew[i][j];
 335
 336   free(diff);
 337   for (i=0;i<alldim;i++)
 338     free(dnew[i]);
 339   free(dnew);
 340 }
 341
 342 void make_iteration(double **dynamics,
 343                     double **delta)
 344 {
 345   double **dnew;
 346   long i,j,k;
 347
 348   check_alloc(dnew=(double**)malloc(sizeof(double*)*alldim));
 349   for (i=0;i<alldim;i++)
 350     check_alloc(dnew[i]=(double*)malloc(sizeof(double)*alldim));
 351
 352   for (i=0;i<alldim;i++) {
 353     for (j=0;j<DIMENSION;j++) {
 354       dnew[i][j]=dynamics[j][0]*delta[i][0];
 355       for (k=1;k<alldim;k++)
 356         dnew[i][j] += dynamics[j][k]*delta[i][k];
 357     }
 358     for (j=DIMENSION;j<alldim;j++)
 359       dnew[i][j]=delta[i][j-1];
 360   }
 361
 362   for (i=0;i<alldim;i++)
 363     for (j=0;j<alldim;j++)
 364       delta[i][j]=dnew[i][j];
 365
 366   for (i=0;i<alldim;i++)
 367     free(dnew[i]);
 368   free(dnew);
 369 }
 370
 371 int main(int argc,char **argv)
 372 {
 373   char stdi=0;
 374   double **delta,**dynamics,*lfactor;
 375   double *factor,dim;
 376   double *hseries;
 377   double *interval,*min,*av,*var,maxinterval;
 378   long start,i,j;
 379   time_t lasttime,newtime;
 380   FILE *file=NULL;
 381
 382   if (scan_help(argc,argv))
 383     show_options(argv[0]);
 384
 385   ITERATIONS=ULONG_MAX;
 386
 387   scan_options(argc,argv);
 388 #ifndef OMIT_WHAT_I_DO
 389   if (verbosity&VER_INPUT)
 390     what_i_do(argv[0],WID_STR);
 391 #endif
 392
 393   infile=search_datafile(argc,argv,NULL,verbosity);
 394   if (infile == NULL)
 395     stdi=1;
 396
 397   if (outfile == NULL) {
 398     if (!stdi) {
 399       check_alloc(outfile=(char*)calloc(strlen(infile)+7,(size_t)1));
 400       strcpy(outfile,infile);
 401       strcat(outfile,".lyaps");
 402     }
 403     else {
 404       check_alloc(outfile=(char*)calloc((size_t)12,(size_t)1));
 405       strcpy(outfile,"stdin.lyaps");
 406     }
 407   }
 408   if (!stdo)
 409     test_outfile(outfile);
 410
 411   alldim=DIMENSION*EMBED;
 412
 413   if (COLUMNS == NULL)
 414     series=(double**)get_multi_series(infile,&LENGTH,exclude,&DIMENSION,"",
 415                                       dimset,verbosity);
 416   else
 417     series=(double**)get_multi_series(infile,&LENGTH,exclude,&DIMENSION,
 418                                       COLUMNS,dimset,verbosity);
 419
 420   if (MINNEIGHBORS > (LENGTH-DELAY*(EMBED-1)-1)) {
 421     fprintf(stderr,"Your time series is not long enough to find %d neighbors!"
 422             " Exiting.\n",MINNEIGHBORS);
 423     exit(LYAP_SPEC_DATA_TOO_SHORT);
 424   }
 425
 426   check_alloc(min=(double*)malloc(sizeof(double)*DIMENSION));
 427   check_alloc(interval=(double*)malloc(sizeof(double)*DIMENSION));
 428   check_alloc(av=(double*)malloc(sizeof(double)*DIMENSION));
 429   check_alloc(var=(double*)malloc(sizeof(double)*DIMENSION));
 430   check_alloc(averr=(double*)malloc(sizeof(double)*DIMENSION));
 431   maxinterval=0.0;
 432   for (i=0;i<DIMENSION;i++) {
 433     averr[i]=0.0;
 434     rescale_data(series[i],LENGTH,&min[i],&interval[i]);
 435     if (interval[i] > maxinterval)
 436       maxinterval=interval[i];
 437     variance(series[i],LENGTH,&av[i],&var[i]);
 438   }
 439
 440   if (INVERSE) {
 441     check_alloc(hseries=(double*)malloc(sizeof(double)*LENGTH));
 442     for (j=0;j<DIMENSION;j++) {
 443       for (i=0;i<LENGTH;i++)
 444         hseries[LENGTH-1-i]=series[j][i];
 445       for (i=0;i<LENGTH;i++)
 446         series[j][i]=hseries[i];
 447     }
 448     free(hseries);
 449   }
 450
 451   if (!epsset)
 452     epsmin=1./1000.;
 453   else
 454     epsmin /= maxinterval;
 455
 456   check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
 457   for (i=0;i<BOX;i++)
 458     check_alloc(box[i]=(long*)malloc(sizeof(long)*BOX));
 459
 460   check_alloc(list=(long*)malloc(sizeof(long)*LENGTH));
 461   check_alloc(found=(unsigned long*)malloc(sizeof(long)*LENGTH));
 462
 463   check_alloc(dynamics=(double**)malloc(sizeof(double*)*DIMENSION));
 464   for (i=0;i<DIMENSION;i++)
 465     check_alloc(dynamics[i]=(double*)malloc(sizeof(double)*alldim));
 466   check_alloc(factor=(double*)malloc(sizeof(double)*alldim));
 467   check_alloc(lfactor=(double*)malloc(sizeof(double)*alldim));
 468   check_alloc(delta=(double**)malloc(sizeof(double*)*alldim));
 469   for (i=0;i<alldim;i++)
 470     check_alloc(delta[i]=(double*)malloc(sizeof(double)*alldim));
 471
 472   check_alloc(vec=(double*)malloc(sizeof(double)*(alldim+1)));
 473   check_alloc(mat=(double**)malloc(sizeof(double*)*(alldim+1)));
 474   for (i=0;i<=alldim;i++)
 475     check_alloc(mat[i]=(double*)malloc(sizeof(double)*(alldim+1)));
 476
 477   indexes=(unsigned int**)make_multi_index(DIMENSION,EMBED,DELAY);
 478
 479   rnd_init(0x098342L);
 480   for (i=0;i<10000;i++)
 481     rnd_long();
 482   for (i=0;i<alldim;i++) {
 483     factor[i]=0.0;
 484     for (j=0;j<alldim;j++)
 485       delta[i][j]=(double)rnd_long()/(double)ULONG_MAX;
 486   }
 487   gram_schmidt(delta,lfactor);
 488
 489   start=ITERATIONS;
 490   if (start>(LENGTH-DELAY))
 491     start=LENGTH-DELAY;
 492
 493   if (!stdo) {
 494     file=fopen(outfile,"w");
 495     if (verbosity&VER_INPUT)
 496       fprintf(stderr,"Opened %s for writing\n",outfile);
 497   }
 498   else {
 499     if (verbosity&VER_INPUT)
 500       fprintf(stderr,"Writing to stdout\n");
 501   }
 502
 503   check_alloc(abstand=(double*)malloc(sizeof(double)*LENGTH));
 504
 505   time(&lasttime);
 506   for (i=(EMBED-1)*DELAY;i<start;i++) {
 507     count++;
 508     make_dynamics(dynamics,i);
 509     make_iteration(dynamics,delta);
 510     gram_schmidt(delta,lfactor);
 511     for (j=0;j<alldim;j++) {
 512       factor[j] += log(lfactor[j])/(double)DELAY;
 513     }
 514     if (((time(&newtime)-lasttime) > OUT) || (i == (start-1))) {
 515       time(&lasttime);
 516       if (!stdo) {
 517         fprintf(file,"%ld ",count);
 518         for (j=0;j<alldim;j++)
 519           fprintf(file,"%e ",factor[j]/count);
 520         fprintf(file,"\n");
 521         fflush(file);
 522       }
 523       else {
 524         fprintf(stdout,"%ld ",count);
 525         for (j=0;j<alldim;j++)
 526           fprintf(stdout,"%e ",factor[j]/count);
 527         fprintf(stdout,"\n");
 528       }
 529     }
 530   }
 531
 532   dim=0.0;
 533   for (i=0;i<alldim;i++) {
 534     dim += factor[i];
 535     if (dim < 0.0)
 536       break;
 537   }
 538   if (i < alldim)
 539     dim=i+(dim-factor[i])/fabs(factor[i]);
 540   else
 541     dim=alldim;
 542   if (!stdo) {
 543     fprintf(file,"#Average relative forecast errors:= ");
 544     for (i=0;i<DIMENSION;i++)
 545       fprintf(file,"%e ",sqrt(averr[i]/count)/var[i]);
 546     fprintf(file,"\n");
 547     fprintf(file,"#Average absolute forecast errors:= ");
 548     for (i=0;i<DIMENSION;i++)
 549       fprintf(file,"%e ",sqrt(averr[i]/count)*interval[i]);
 550     fprintf(file,"\n");
 551     fprintf(file,"#Average Neighborhood Size= %e\n",aveps*maxinterval/count);
 552     fprintf(file,"#Average num. of neighbors= %e\n",avneig/count);
 553     fprintf(file,"#estimated KY-Dimension= %f\n",dim);
 554   }
 555   else {
 556     fprintf(stdout,"#Average relative forecast errors:= ");
 557     for (i=0;i<DIMENSION;i++)
 558       fprintf(stdout,"%e ",sqrt(averr[i]/count)/var[i]);
 559     fprintf(stdout,"\n");
 560     fprintf(stdout,"#Average absolute forecast errors:= ");
 561     for (i=0;i<DIMENSION;i++)
 562       fprintf(stdout,"%e ",sqrt(averr[i]/count)*interval[i]);
 563     fprintf(stdout,"\n");
 564     fprintf(stdout,"#Average Neighborhood Size= %e\n",aveps*maxinterval/count);
 565     fprintf(stdout,"#Average num. of neighbors= %e\n",avneig/count);
 566     fprintf(stdout,"#estimated KY-Dimension= %f\n",dim);
 567   }
 568   if (!stdo)
 569     fclose(file);
 570
 571   free(abstand);
 572
 573   return 0;
 574 }