IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / 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 }