lyap_k.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: Sep 3, 1999*/
  21 #include <math.h>
  22 #include <limits.h>
  23 #include <stdio.h>
  24 #include <stdlib.h>
  25 #include <string.h>
  26 #include "routines/tsa.h"
  27
  28 #define WID_STR "Estimates the maximal Lyapunov exponent using the Kantz\n\t\
  29 algorithm"
  30
  31 #define BOX 128
  32 const unsigned int ibox=BOX-1;
  33
  34 unsigned long length=ULONG_MAX;
  35 unsigned long exclude=0;
  36 unsigned long reference=ULONG_MAX;
  37 unsigned int maxdim=2;
  38 unsigned int mindim=2;
  39 unsigned int delay=1;
  40 unsigned int column=1;
  41 unsigned int epscount=5;
  42 unsigned int maxiter=50;
  43 unsigned int window=0;
  44 unsigned int verbosity=0xff;
  45 double epsmin=1.e-3,epsmax=1.e-2;
  46 char eps0set=0,eps1set=0;
  47 char *outfile=NULL;
  48 char *infile=NULL;
  49
  50 double *series,**lyap;
  51 long box[BOX][BOX],*liste,**lfound,*found,**count;
  52 double max,min;
  53
  54 void show_options(char *progname)
  55 {
  56   what_i_do(progname,WID_STR);
  57
  58   fprintf(stderr," Usage: %s [options]\n",progname);
  59   fprintf(stderr," Options:\n");
  60   fprintf(stderr,"Everything not being a valid option will be "
  61           "interpreted as a possible datafile.\nIf no datafile "
  62           "is given stdin is read. Just - also means stdin\n");
  63   fprintf(stderr,"\t-l # of data [default: whole file]\n");
  64   fprintf(stderr,"\t-x # of lines to be ignored [default: 0]\n");
  65   fprintf(stderr,"\t-c column to read [default: 1]\n");
  66   fprintf(stderr,"\t-M maxdim [default: 2]\n");
  67   fprintf(stderr,"\t-m mindim [default: 2]\n");
  68   fprintf(stderr,"\t-d delay [default: 1]\n");
  69   fprintf(stderr,"\t-r mineps [default: (data interval)/1000]\n");
  70   fprintf(stderr,"\t-R maxeps [default: (data interval)/100]\n");
  71   fprintf(stderr,"\t-# # of eps [default: 5]\n");
  72   fprintf(stderr,"\t-n # of reference points [default: # of data]\n");
  73   fprintf(stderr,"\t-s # of iterations [default: 50]\n");
  74   fprintf(stderr,"\t-t time window [default: 0]\n");
  75   fprintf(stderr,"\t-o outfile [default: 'datafile'.lyap]\n");
  76   fprintf(stderr,"\t-V verbosity level [default: 3]\n\t\t"
  77           "0='only panic messages'\n\t\t"
  78           "1='+ input/output messages'\n\t\t"
  79           "2='+ plus statistics'\n");
  80   fprintf(stderr,"\t-h show these options\n");
  81   exit(0);
  82 }
  83
  84 void scan_options(int n,char **str)
  85 {
  86   char *out;
  87
  88   if ((out=check_option(str,n,'l','u')) != NULL)
  89     sscanf(out,"%lu",&length);
  90   if ((out=check_option(str,n,'x','u')) != NULL)
  91     sscanf(out,"%lu",&exclude);
  92   if ((out=check_option(str,n,'c','u')) != NULL)
  93     sscanf(out,"%u",&column);
  94   if ((out=check_option(str,n,'M','u')) != NULL)
  95     sscanf(out,"%u",&maxdim);
  96   if ((out=check_option(str,n,'m','u')) != NULL)
  97     sscanf(out,"%u",&mindim);
  98   if ((out=check_option(str,n,'d','u')) != NULL)
  99     sscanf(out,"%u",&delay);
 100   if ((out=check_option(str,n,'r','f')) != NULL) {
 101     eps0set=1;
 102     sscanf(out,"%lf",&epsmin);
 103   }
 104   if ((out=check_option(str,n,'R','f')) != NULL) {
 105     eps1set=1;
 106     sscanf(out,"%lf",&epsmax);
 107   }
 108   if ((out=check_option(str,n,'#','u')) != NULL)
 109     sscanf(out,"%u",&epscount);
 110   if ((out=check_option(str,n,'n','u')) != NULL)
 111     sscanf(out,"%lu",&reference);
 112   if ((out=check_option(str,n,'s','u')) != NULL)
 113     sscanf(out,"%u",&maxiter);
 114   if ((out=check_option(str,n,'t','u')) != NULL)
 115     sscanf(out,"%u",&window);
 116   if ((out=check_option(str,n,'V','u')) != NULL)
 117     sscanf(out,"%u",&verbosity);
 118   if ((out=check_option(str,n,'o','o')) != NULL)
 119     if (strlen(out) > 0)
 120       outfile=out;
 121 }
 122
 123 void put_in_boxes(double eps)
 124 {
 125   unsigned long i;
 126   long j,k;
 127   static unsigned long blength;
 128
 129   blength=length-(maxdim-1)*delay-maxiter;
 130
 131   for (i=0;i<BOX;i++)
 132     for (j=0;j<BOX;j++)
 133       box[i][j]= -1;
 134
 135   for (i=0;i<blength;i++) {
 136     j=(long)(series[i]/eps)&ibox;
 137     k=(long)(series[i+delay]/eps)&ibox;
 138     liste[i]=box[j][k];
 139     box[j][k]=i;
 140   }
 141 }
 142
 143 void lfind_neighbors(long act,double eps)
 144 {
 145   unsigned int hi,k,k1;
 146   long i,j,i1,i2,j1,element;
 147   static long lwindow;
 148   double dx,eps2=sqr(eps);
 149
 150   lwindow=(long)window;
 151   for (hi=0;hi<maxdim-1;hi++)
 152     found[hi]=0;
 153   i=(long)(series[act]/eps)&ibox;
 154   j=(long)(series[act+delay]/eps)&ibox;
 155   for (i1=i-1;i1<=i+1;i1++) {
 156     i2=i1&ibox;
 157     for (j1=j-1;j1<=j+1;j1++) {
 158       element=box[i2][j1&ibox];
 159       while (element != -1) {
 160         if ((element < (act-lwindow)) || (element > (act+lwindow))) {
 161           dx=sqr(series[act]-series[element]);
 162           if (dx <= eps2) {
 163             for (k=1;k<maxdim;k++) {
 164               k1=k*delay;
 165               dx += sqr(series[act+k1]-series[element+k1]);
 166               if (dx <= eps2) {
 167                 k1=k-1;
 168                 lfound[k1][found[k1]]=element;
 169                 found[k1]++;
 170               }
 171               else
 172                 break;
 173             }
 174           }
 175         }
 176         element=liste[element];
 177       }
 178     }
 179   }
 180 }
 181
 182 void iterate_points(long act)
 183 {
 184   double **lfactor;
 185   double *dx;
 186   unsigned int i,j,l,l1;
 187   long k,element,**lcount;
 188
 189   check_alloc(lfactor=(double**)malloc(sizeof(double*)*(maxdim-1)));
 190   check_alloc(lcount=(long**)malloc(sizeof(long*)*(maxdim-1)));
 191   for (i=0;i<maxdim-1;i++) {
 192     check_alloc(lfactor[i]=(double*)malloc(sizeof(double)*(maxiter+1)));
 193     check_alloc(lcount[i]=(long*)malloc(sizeof(long)*(maxiter+1)));
 194   }
 195   check_alloc(dx=(double*)malloc(sizeof(double)*(maxiter+1)));
 196
 197   for (i=0;i<=maxiter;i++)
 198     for (j=0;j<maxdim-1;j++) {
 199       lfactor[j][i]=0.0;
 200       lcount[j][i]=0;
 201     }
 202
 203   for (j=mindim-2;j<maxdim-1;j++) {
 204     for (k=0;k<found[j];k++) {
 205       element=lfound[j][k];
 206       for (i=0;i<=maxiter;i++)
 207         dx[i]=sqr(series[act+i]-series[element+i]);
 208       for (l=1;l<j+2;l++) {
 209         l1=l*delay;
 210         for (i=0;i<=maxiter;i++)
 211           dx[i] += sqr(series[act+i+l1]-series[element+l1+i]);
 212       }
 213       for (i=0;i<=maxiter;i++)
 214         if (dx[i] > 0.0){
 215           lcount[j][i]++;
 216           lfactor[j][i] += dx[i];
 217         }
 218     }
 219   }
 220   for (i=mindim-2;i<maxdim-1;i++)
 221     for (j=0;j<=maxiter;j++)
 222       if (lcount[i][j]) {
 223         count[i][j]++;
 224         lyap[i][j] += log(lfactor[i][j]/lcount[i][j])/2.0;
 225       }
 226
 227   for (i=0;i<maxdim-1;i++){
 228     free(lfactor[i]);
 229     free(lcount[i]);
 230   }
 231   free(lcount);
 232   free(lfactor);
 233   free(dx);
 234 }
 235
 236 int main(int argc,char **argv)
 237 {
 238   char stdi=0;
 239   double eps_fak;
 240   double epsilon;
 241   unsigned int i,j,l;
 242   FILE *fout;
 243
 244   if (scan_help(argc,argv))
 245     show_options(argv[0]);
 246
 247   scan_options(argc,argv);
 248 #ifndef OMIT_WHAT_I_DO
 249   if (verbosity&VER_INPUT)
 250     what_i_do(argv[0],WID_STR);
 251 #endif
 252
 253   infile=search_datafile(argc,argv,&column,verbosity);
 254   if (infile == NULL)
 255     stdi=1;
 256
 257   if (outfile == NULL) {
 258     if (!stdi) {
 259       check_alloc(outfile=(char*)calloc(strlen(infile)+6,1));
 260       sprintf(outfile,"%s.lyap",infile);
 261     }
 262     else {
 263       check_alloc(outfile=(char*)calloc(11,1));
 264       sprintf(outfile,"stdin.lyap");
 265     }
 266   }
 267   test_outfile(outfile);
 268
 269   series=get_series(infile,&length,exclude,column,verbosity);
 270   rescale_data(series,length,&min,&max);
 271
 272   if (eps0set)
 273     epsmin /= max;
 274   if (eps1set)
 275     epsmax /= max;
 276
 277   if (epsmin >= epsmax) {
 278     epsmax=epsmin;
 279     epscount=1;
 280   }
 281
 282   if (reference > (length-maxiter-(maxdim-1)*delay))
 283     reference=length-maxiter-(maxdim-1)*delay;
 284   if ((maxiter+(maxdim-1)*delay) >= length) {
 285     fprintf(stderr,"Too few points to handle these parameters!\n");
 286     exit(LYAP_K__MAXITER_TOO_LARGE);
 287   }
 288
 289   if (maxdim < 2)
 290     maxdim=2;
 291   if (mindim < 2)
 292     mindim=2;
 293   if (mindim > maxdim)
 294     maxdim=mindim;
 295
 296   check_alloc(liste=(long*)malloc(sizeof(long)*(length)));
 297   check_alloc(found=(long*)malloc(sizeof(long)*(maxdim-1)));
 298   check_alloc(lfound=(long**)malloc(sizeof(long*)*(maxdim-1)));
 299   for (i=0;i<maxdim-1;i++)
 300     check_alloc(lfound[i]=(long*)malloc(sizeof(long)*(length)));
 301   check_alloc(count=(long**)malloc(sizeof(long*)*(maxdim-1)));
 302   for (i=0;i<maxdim-1;i++)
 303     check_alloc(count[i]=(long*)malloc(sizeof(long)*(maxiter+1)));
 304   check_alloc(lyap=(double**)malloc(sizeof(double*)*(maxdim-1)));
 305   for (i=0;i<maxdim-1;i++)
 306     check_alloc(lyap[i]=(double*)malloc(sizeof(double)*(maxiter+1)));
 307
 308   if (epscount == 1)
 309     eps_fak=1.0;
 310   else
 311     eps_fak=pow(epsmax/epsmin,1.0/(double)(epscount-1));
 312
 313   fout=fopen(outfile,"w");
 314   if (verbosity&VER_INPUT)
 315     fprintf(stderr,"Opened %s for writing\n",outfile);
 316   for (l=0;l<epscount;l++) {
 317     epsilon=epsmin*pow(eps_fak,(double)l);
 318     for (i=0;i<maxdim-1;i++)
 319       for (j=0;j<=maxiter;j++) {
 320         count[i][j]=0;
 321         lyap[i][j]=0.0;
 322       }
 323     put_in_boxes(epsilon);
 324     for (i=0;i<reference;i++) {
 325       lfind_neighbors(i,epsilon);
 326       iterate_points(i);
 327     }
 328     if (verbosity&VER_USR1)
 329       fprintf(stderr,"epsilon= %e\n",epsilon*max);
 330     for (i=mindim-2;i<maxdim-1;i++) {
 331       fprintf(fout,"#epsilon= %e  dim= %d\n",epsilon*max,i+2);
 332       for (j=0;j<=maxiter;j++)
 333         if (count[i][j])
 334           fprintf(fout,"%d %e %ld\n",j,lyap[i][j]/count[i][j],count[i][j]);
 335       fprintf(fout,"\n");
 336     }
 337     fflush(fout);
 338   }
 339   fclose(fout);
 340   return 0;
 341 }