binaries/src/disembl/Tisean_3.0.1/source_c/lyap_r.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: Apr 25, 2002 */
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <math.h>
  24 #include <limits.h>
  25 #include <string.h>
  26 #include "routines/tsa.h"
  27
  28 #define WID_STR "Estimates the maximal Lyapunov exponent; Rosenstein et al."
  29
  30 #define NMAX 256
  31
  32 char *outfile=NULL;
  33 char *infile=NULL;
  34 char epsset=0;
  35 double *series,*lyap;
  36 long box[NMAX][NMAX],*list;
  37 unsigned int dim=2,delay=1,steps=10,mindist=0;
  38 unsigned int column=1;
  39 unsigned int verbosity=0xff;
  40 const unsigned int nmax=NMAX-1;
  41 unsigned long length=ULONG_MAX,exclude=0;
  42 long *found;
  43 double eps0=1.e-3,eps,epsinv;
  44
  45 void show_options(char *progname)
  46 {
  47   what_i_do(progname,WID_STR);
  48   fprintf(stderr," Usage: %s [options]\n",progname);
  49   fprintf(stderr," Options:\n");
  50   fprintf(stderr,"Everything not being a valid option will be interpreted"
  51           " as a possible"
  52           " datafile.\nIf no datafile is given stdin is read. Just - also"
  53           " means stdin\n");
  54   fprintf(stderr,"\t-l # of datapoints [default is whole file]\n");
  55   fprintf(stderr,"\t-x # of lines to be ignored [default is 0]\n");
  56   fprintf(stderr,"\t-c column to read[default 1]\n");
  57   fprintf(stderr,"\t-m embedding dimension [default 2]\n");
  58   fprintf(stderr,"\t-d delay  [default 1]\n");
  59   fprintf(stderr,"\t-t time window to omit [default 0]\n");
  60   fprintf(stderr,"\t-r epsilon size to start with [default "
  61           "(data interval)/1000]\n");
  62   fprintf(stderr,"\t-s # of iterations [default 10]\n");
  63   fprintf(stderr,"\t-o name of output file [default 'datafile'.ros]\n");
  64   fprintf(stderr,"\t-V verbosity level [default 3]\n\t\t"
  65           "0='only panic messages'\n\t\t"
  66           "1='+ input/output messages'\n\t\t"
  67           "2='+ give more detailed information about the length scales\n");
  68   fprintf(stderr,"\t-h show these options\n");
  69   fprintf(stderr,"\n");
  70   exit(0);
  71 }
  72
  73 void scan_options(int n,char **argv)
  74 {
  75   char *out;
  76
  77   if ((out=check_option(argv,n,'l','u')) != NULL)
  78     sscanf(out,"%lu",&length);
  79   if ((out=check_option(argv,n,'x','u')) != NULL)
  80     sscanf(out,"%lu",&exclude);
  81   if ((out=check_option(argv,n,'c','u')) != NULL)
  82     sscanf(out,"%u",&column);
  83   if ((out=check_option(argv,n,'m','u')) != NULL)
  84     sscanf(out,"%u",&dim);
  85   if ((out=check_option(argv,n,'d','u')) != NULL)
  86     sscanf(out,"%u",&delay);
  87   if ((out=check_option(argv,n,'t','u')) != NULL)
  88     sscanf(out,"%u",&mindist);
  89   if ((out=check_option(argv,n,'r','f')) != NULL) {
  90     epsset=1;
  91     sscanf(out,"%lf",&eps0);
  92   }
  93   if ((out=check_option(argv,n,'s','u')) != NULL)
  94     sscanf(out,"%u",&steps);
  95   if ((out=check_option(argv,n,'V','u')) != NULL)
  96     sscanf(out,"%u",&verbosity);
  97   if ((out=check_option(argv,n,'o','o')) != NULL)
  98     if (strlen(out) > 0)
  99       outfile=out;
 100 }
 101
 102 void put_in_boxes(void)
 103 {
 104   int i,j,x,y,del;
 105
 106   for (i=0;i<NMAX;i++)
 107     for (j=0;j<NMAX;j++)
 108       box[i][j]= -1;
 109
 110   del=delay*(dim-1);
 111   for (i=0;i<length-del-steps;i++) {
 112     x=(int)(series[i]*epsinv)&nmax;
 113     y=(int)(series[i+del]*epsinv)&nmax;
 114     list[i]=box[x][y];
 115     box[x][y]=i;
 116   }
 117 }
 118
 119 char make_iterate(long act)
 120 {
 121   char ok=0;
 122   int x,y,i,j,i1,k,del1=dim*delay;
 123   long element,minelement= -1;
 124   double dx,mindx=1.0;
 125
 126   x=(int)(series[act]*epsinv)&nmax;
 127   y=(int)(series[act+delay*(dim-1)]*epsinv)&nmax;
 128   for (i=x-1;i<=x+1;i++) {
 129     i1=i&nmax;
 130     for (j=y-1;j<=y+1;j++) {
 131       element=box[i1][j&nmax];
 132       while (element != -1) {
 133         if (labs(act-element) > mindist) {
 134           dx=0.0;
 135           for (k=0;k<del1;k+=delay) {
 136             dx += (series[act+k]-series[element+k])*
 137               (series[act+k]-series[element+k]);
 138             if (dx > eps*eps)
 139               break;
 140           }
 141           if (k==del1) {
 142             if (dx < mindx) {
 143               ok=1;
 144               if (dx > 0.0) {
 145                 mindx=dx;
 146                 minelement=element;
 147               }
 148             }
 149           }
 150         }
 151         element=list[element];
 152       }
 153     }
 154   }
 155   if ((minelement != -1) ) {
 156     act--;
 157     minelement--;
 158     for (i=0;i<=steps;i++) {
 159       act++;
 160       minelement++;
 161       dx=0.0;
 162       for (j=0;j<del1;j+=delay) {
 163         dx += (series[act+j]-series[minelement+j])*
 164           (series[act+j]-series[minelement+j]);
 165       }
 166       if (dx > 0.0) {
 167         found[i]++;
 168         lyap[i] += log(dx);
 169       }
 170     }
 171   }
 172   return ok;
 173 }
 174
 175 int main(int argc,char **argv)
 176 {
 177   char stdi=0,*done,alldone;
 178   int i;
 179   long n;
 180   long maxlength;
 181   double min,max;
 182   FILE *file;
 183
 184   if (scan_help(argc,argv))
 185     show_options(argv[0]);
 186
 187   scan_options(argc,argv);
 188 #ifndef OMIT_WHAT_I_DO
 189   if (verbosity&VER_INPUT)
 190     what_i_do(argv[0],WID_STR);
 191 #endif
 192
 193   infile=search_datafile(argc,argv,&column,verbosity);
 194   if (infile == NULL)
 195     stdi=1;
 196
 197   if (outfile == NULL) {
 198     if (!stdi) {
 199       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
 200       strcpy(outfile,infile);
 201       strcat(outfile,".ros");
 202     }
 203     else {
 204       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
 205       strcpy(outfile,"stdin.ros");
 206     }
 207   }
 208   test_outfile(outfile);
 209
 210   series=(double*)get_series(infile,&length,exclude,column,verbosity);
 211   rescale_data(series,length,&min,&max);
 212
 213   if (epsset)
 214     eps0 /= max;
 215
 216   check_alloc(list=(long*)malloc(length*sizeof(long)));
 217   check_alloc(lyap=(double*)malloc((steps+1)*sizeof(double)));
 218   check_alloc(found=(long*)malloc((steps+1)*sizeof(long)));
 219   check_alloc(done=(char*)malloc(length));
 220
 221   for (i=0;i<=steps;i++) {
 222     lyap[i]=0.0;
 223     found[i]=0;
 224   }
 225   for (i=0;i<length;i++)
 226     done[i]=0;
 227
 228   maxlength=length-delay*(dim-1)-steps-1-mindist;
 229   alldone=0;
 230   file=fopen(outfile,"w");
 231   if (verbosity&VER_INPUT)
 232     fprintf(stderr,"Opened %s for writing\n",outfile);
 233   for (eps=eps0;!alldone;eps*=1.1) {
 234     epsinv=1.0/eps;
 235     put_in_boxes();
 236     alldone=1;
 237     for (n=0;n<=maxlength;n++) {
 238       if (!done[n])
 239         done[n]=make_iterate(n);
 240       alldone &= done[n];
 241     }
 242     if (verbosity&VER_USR1)
 243       fprintf(stderr,"epsilon: %e already found: %ld\n",eps*max,found[0]);
 244   }
 245   for (i=0;i<=steps;i++)
 246     if (found[i])
 247       fprintf(file,"%d %e\n",i,lyap[i]/found[i]/2.0);
 248   fclose(file);
 249
 250   return 0;
 251 }