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