IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / d2.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 May 10, 2000 */
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include <math.h>
25 #include <limits.h>
26 #include <time.h>
27 #include "routines/tsa.h"
28
29 #define WID_STR "Estimates the correlation sum, -dimension and -entropy"
30
31 /* output is written every WHEN seconds */
32 #define WHEN 120
33 /* Size of the field for box assisted neighbour searching 
34    (has to be a power of 2)*/
35 #define NMAX 256
36 /* Size of the box for the scramble routine */
37 #define SCBOX 4096
38
39 double **series;
40 long *scr;
41 char dimset=0,rescale_set=0,eps_min_set=0,eps_max_set=0;
42 char *FOUT=NULL;
43 double epsfactor,epsinv,lneps,lnfac;
44 double EPSMAX=1.0,EPSMIN=1.e-3;
45 double min,interval;
46 int imax=NMAX-1,howoften1,imin;
47 long box[NMAX][NMAX],*list,boxc1[NMAX],*listc1;
48 unsigned long nmax;
49 double **found,*norm;
50 unsigned long MINDIST=0,MAXFOUND=1000;
51 unsigned long length=ULONG_MAX,exclude=0;
52 unsigned int DIM=1,EMBED=10,HOWOFTEN=100,DELAY=1;
53 unsigned int verbosity=0x1;
54 char *column=NULL;
55 char *infile=NULL;
56
57 void show_options(char *progname)
58 {
59   what_i_do(progname,WID_STR);
60   fprintf(stderr,"  Usage: %s [options]\n",progname);
61   fprintf(stderr,"  Options:\n");
62   fprintf(stderr,"Everything not being a valid option will be interpreted"
63           " as a possible"
64           " datafile.\nIf no datafile is given stdin is read. Just - also"
65           " means stdin\n");
66   fprintf(stderr,"\t-l datapoints [default is whole file]\n");
67   fprintf(stderr,"\t-x exclude # points [default 0]\n");
68   fprintf(stderr,"\t-d delay  [default 1]\n");
69   fprintf(stderr,"\t-M # of components, max. embedding dim. [default 1,10]\n");
70   fprintf(stderr,"\t-c columns [default 1,...,# of components]\n");
71   fprintf(stderr,"\t-t theiler-window [default 0]\n");
72   fprintf(stderr,"\t-R max-epsilon "
73           "[default: max data interval]\n");
74   fprintf(stderr,"\t-r min-epsilon [default: (max data interval)/1000]\n");
75   fprintf(stderr,"\t-# #-of-epsilons [default 100]\n");
76   fprintf(stderr,"\t-N max-#-of-pairs (0 means all) [default 1000]\n");
77   fprintf(stderr,"\t-E use rescaled data [default: not rescaled]\n");
78   fprintf(stderr," \t-o outfiles"
79           " [without exts.! default datafile[.d2][.h2][.stat][.c2]]\n");
80   fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
81           "0='only panic messages'\n\t\t"
82           "1='+ input/output messages'\n\t\t"
83           "2='+ output message each time output is done\n");
84   
85   fprintf(stderr,"\t-h show these options\n");
86   fprintf(stderr,"\n");
87   exit(0);
88 }
89
90 void scan_options(int n,char **argv)
91 {
92   char *out;
93   
94   if ((out=check_option(argv,n,'l','u')) != NULL)
95     sscanf(out,"%lu",&length);
96   if ((out=check_option(argv,n,'x','u')) != NULL)
97     sscanf(out,"%lu",&exclude);
98   if ((out=check_option(argv,n,'c','s')) != NULL)
99     column=out;
100   if ((out=check_option(argv,n,'d','u')) != NULL)
101     sscanf(out,"%u",&DELAY);
102   if ((out=check_option(argv,n,'M','2')) != NULL) {
103     sscanf(out,"%u,%u",&DIM,&EMBED);
104     dimset=1;
105   }
106   if ((out=check_option(argv,n,'t','u')) != NULL)
107     sscanf(out,"%lu",&MINDIST);
108   if ((out=check_option(argv,n,'R','f')) != NULL) {
109     sscanf(out,"%lf",&EPSMAX);
110     eps_max_set=1;
111   }
112   if ((out=check_option(argv,n,'r','f')) != NULL) {
113     sscanf(out,"%lf",&EPSMIN);
114     eps_min_set=1;
115   }
116   if ((out=check_option(argv,n,'#','u')) != NULL)
117     sscanf(out,"%u",&HOWOFTEN);
118   if ((out=check_option(argv,n,'N','u')) != NULL) {
119     sscanf(out,"%lu",&MAXFOUND);
120     if (MAXFOUND == 0)
121       MAXFOUND=ULONG_MAX;
122   }
123   if ((out=check_option(argv,n,'E','n')) != NULL)
124     rescale_set=1;
125   if ((out=check_option(argv,n,'V','u')) != NULL)
126     sscanf(out,"%u",&verbosity);
127   if ((out=check_option(argv,n,'o','o')) != NULL)
128     if (strlen(out) > 0)
129       FOUT=out;
130 }
131       
132 void scramble(void)
133 {
134   long i,j,k,m;
135   unsigned long rnd,rndf,hlength,allscr=0;
136   long *scfound,*scnhelp,scnfound;
137   long scbox[SCBOX],lswap,element,scbox1=SCBOX-1;
138   double *rz,*schelp,sceps=(double)SCBOX-0.001,swap;
139   
140   hlength=length-(EMBED-1)*DELAY;
141
142   if (sizeof(long) == 8) {
143     rndf=13*13*13*13;
144     rndf=rndf*rndf*rndf*13;
145     rnd=0x849178L;
146   }
147   else {
148     rndf=69069;
149     rnd=0x234571L;
150   }
151   for (i=0;i<1000;i++)
152     rnd=rnd*rndf+1;
153
154   check_alloc(rz=(double*)malloc(sizeof(double)*hlength));
155   check_alloc(scfound=(long*)malloc(sizeof(long)*hlength));
156   check_alloc(scnhelp=(long*)malloc(sizeof(long)*hlength));
157   check_alloc(schelp=(double*)malloc(sizeof(double)*hlength));
158
159   for (i=0;i<hlength;i++)
160     rz[i]=(double)(rnd=rnd*rndf+1)/ULONG_MAX;
161   
162   for (i=0;i<SCBOX;i++)
163     scbox[i]= -1;
164   for (i=0;i<hlength;i++) {
165     m=(int)(rz[i]*sceps)&scbox1;
166     scfound[i]=scbox[m];
167     scbox[m]=i;
168   }
169   for (i=0;i<SCBOX;i++) {
170     scnfound=0;
171     element=scbox[i];
172     while(element != -1) {
173       scnhelp[scnfound]=element;
174       schelp[scnfound++]=rz[element];
175       element=scfound[element];
176     }
177     
178     for (j=0;j<scnfound-1;j++)
179       for (k=j+1;k<scnfound;k++)
180         if (schelp[k] < schelp[j]) {
181           swap=schelp[k];
182           schelp[k]=schelp[j];
183           schelp[j]=swap;
184           lswap=scnhelp[k];
185           scnhelp[k]=scnhelp[j];
186           scnhelp[j]=lswap;
187         }
188     for (j=0;j<scnfound;j++)
189       scr[allscr+j]=scnhelp[j];
190     allscr += scnfound;
191   }
192
193   free(rz);
194   free(scfound);
195   free(schelp);
196 }
197
198 void make_c2_dim(int n)
199 {
200   char small;
201   long i,j,k,x,y,i1,i2,j1,element,n1,maxi,count,hi;
202   double *hs,max,dx;
203   
204   check_alloc(hs=(double*)malloc(sizeof(double)*EMBED*DIM));
205   n1=scr[n];
206
207   count=0;
208   for (i1=0;i1<EMBED;i1++) {
209     i2=i1*DELAY;
210     for (j=0;j<DIM;j++)
211       hs[count++]=series[j][n1+i2];
212   }
213
214   x=(int)(hs[0]*epsinv)&imax;
215   y=(int)(hs[1]*epsinv)&imax;
216   
217   for (i1=x-1;i1<=x+1;i1++) {
218     i2=i1&imax;
219     for (j1=y-1;j1<=y+1;j1++) {
220       element=box[i2][j1&imax];
221       while (element != -1) {
222         if (labs((long)(element-n1)) > MINDIST) {
223           count=0;
224           max=0.0;
225           maxi=howoften1;
226           small=0;
227           for (i=0;i<EMBED;i++) {
228             hi=i*DELAY;
229             for (j=0;j<DIM;j++) {
230               dx=fabs(hs[count]-series[j][element+hi]);
231               if (dx <= EPSMAX) {
232                 if (dx > max) {
233                   max=dx;
234                   if (max < EPSMIN) {
235                     maxi=howoften1;
236                   }
237                   else {
238                     maxi=(lneps-log(max))/lnfac;
239                   }
240                 }
241                 if (count > 0)
242                   for (k=imin;k<=maxi;k++)
243                     found[count][k] += 1.0;
244               }
245               else {
246                 small=1;
247                 break;
248               }
249               count++;
250             }
251             if (small)
252               break;
253           }
254         }
255         element=list[element];
256       }
257     }
258   }
259
260   free(hs);
261 }
262
263 void make_c2_1(int n)
264 {
265   int i,x,i1,maxi;
266   long element,n1;
267   double hs,max;
268   
269   n1=scr[n];
270   hs=series[0][n1];
271   
272   x=(int)(hs*epsinv)&imax;
273   
274   for (i1=x-1;i1<=x+1;i1++) {
275     element=boxc1[i1&imax];
276     while (element != -1) {
277       if (labs(element-n1) > MINDIST) {
278         max=fabs(hs-series[0][element]);
279         if (max <= EPSMAX) {
280           if (max < EPSMIN)
281             maxi=howoften1;
282           else
283             maxi=(lneps-log(max))/lnfac;
284           for (i=imin;i<=maxi;i++)
285             found[0][i] += 1.0;
286         }
287       }
288       element=listc1[element];
289     }
290   }
291 }
292
293 int main(int argc,char **argv)
294 {
295   char smaller,stdi=0;
296   FILE *fout,*fstat;
297   char *outd1,*outc1,*outh1,*outstat;
298   int maxembed;
299   long i1,j1,x,y,sn,n,i,j,n1,n2;
300   long *oscr;
301   long lnorm;
302   double eps,*epsm,EPSMAX1,maxinterval;
303   time_t mytime,lasttime;
304   
305   if (scan_help(argc,argv)) 
306     show_options(argv[0]);
307   
308   scan_options(argc,argv);
309 #ifndef OMIT_WHAT_I_DO
310   if (verbosity&VER_INPUT)
311     what_i_do(argv[0],WID_STR);
312 #endif
313
314   infile=search_datafile(argc,argv,NULL,verbosity);
315   if (infile == NULL)
316     stdi=1;
317   
318   if (FOUT == NULL) {
319     if (!stdi) {
320       check_alloc(FOUT=calloc(strlen(infile)+1,(size_t)1));
321       strcpy(FOUT,infile);
322     }
323     else {
324       check_alloc(FOUT=calloc((size_t)6,(size_t)1));
325       strcpy(FOUT,"stdin");
326     }
327   }
328   if (column == NULL)
329     series=(double**)get_multi_series(infile,&length,exclude,&DIM,"",dimset,
330                                       verbosity);
331   else
332     series=(double**)get_multi_series(infile,&length,exclude,&DIM,column,
333                                       dimset,verbosity);
334
335   if (rescale_set) {
336     for (i=0;i<DIM;i++)
337       rescale_data(series[i],length,&min,&interval);
338     maxinterval=1.0;
339   }
340   else {
341     maxinterval=0.0;
342     for (i=0;i<DIM;i++) {
343       min=interval=series[i][0];
344       for (j=1;j<length;j++) {
345         if (min > series[i][j])
346           min=series[i][j];
347         if (interval < series[i][j])
348           interval=series[i][j];
349       }
350       interval -= min;
351       if (interval > maxinterval)
352         maxinterval=interval;
353     }
354   }
355   if (!eps_max_set)
356     EPSMAX *= maxinterval;
357   if (!eps_min_set)
358     EPSMIN *= maxinterval;
359   EPSMAX=(fabs(EPSMAX)<maxinterval) ? fabs(EPSMAX) : maxinterval;
360   EPSMIN=(fabs(EPSMIN)<EPSMAX) ? fabs(EPSMIN) : EPSMAX/2.;
361   EPSMAX1=EPSMAX;
362
363   howoften1=HOWOFTEN-1;
364   maxembed=DIM*EMBED-1;
365
366   check_alloc(outd1=(char*)calloc(strlen(FOUT)+4,(size_t)1));
367   check_alloc(outc1=(char*)calloc(strlen(FOUT)+4,(size_t)1));
368   check_alloc(outh1=(char*)calloc(strlen(FOUT)+4,(size_t)1));
369   check_alloc(outstat=(char*)calloc(strlen(FOUT)+6,(size_t)1));
370   strcpy(outd1,FOUT);
371   strcpy(outc1,FOUT);
372   strcpy(outh1,FOUT);
373   strcpy(outstat,FOUT);
374   strcat(outd1,".d2");
375   strcat(outc1,".c2");
376   strcat(outh1,".h2");
377   strcat(outstat,".stat");
378   test_outfile(outd1);
379   test_outfile(outc1);
380   test_outfile(outh1);
381   test_outfile(outstat);
382
383   check_alloc(list=(long*)malloc(length*sizeof(long)));
384   check_alloc(listc1=(long*)malloc(length*sizeof(long)));
385   if ((long)(length-(EMBED-1)*DELAY) <= 0) {
386     fprintf(stderr,"Embedding dimension and delay are too large.\n"
387             "The delay vector would be longer than the whole series."
388             " Exiting\n");
389     exit(VECTOR_TOO_LARGE_FOR_LENGTH);
390   }
391   check_alloc(scr=(long*)malloc(sizeof(long)*(length-(EMBED-1)*DELAY)));
392   check_alloc(oscr=(long*)malloc(sizeof(long)*(length-(EMBED-1)*DELAY)));
393   check_alloc(found=(double**)malloc(DIM*EMBED*sizeof(double*)));
394   for (i=0;i<EMBED*DIM;i++)
395     check_alloc(found[i]=(double*)malloc(HOWOFTEN*sizeof(double)));
396   check_alloc(norm=(double*)malloc(HOWOFTEN*sizeof(double)));
397   check_alloc(epsm=(double*)malloc(HOWOFTEN*sizeof(double)));
398   
399   epsinv=1.0/EPSMAX;
400   epsfactor=pow(EPSMAX/EPSMIN,1.0/(double)howoften1);
401   lneps=log(EPSMAX);
402   lnfac=log(epsfactor);
403
404   epsm[0]=EPSMAX;
405   norm[0]=0.0;
406   for (i=1;i<HOWOFTEN;i++) {
407     norm[i]=0.0;
408     epsm[i]=epsm[i-1]/epsfactor;
409   }
410   imin=0;
411
412   scramble();
413   for (i=0;i<(length-(EMBED-1)*DELAY);i++)
414     oscr[scr[i]]=i;
415
416   for (i=0;i<DIM*EMBED;i++)
417     for (j=0;j<HOWOFTEN;j++)
418       found[i][j]=0.0;
419   
420   nmax=length-DELAY*(EMBED-1);
421
422   for (i1=0;i1<NMAX;i1++) {
423     boxc1[i1]= -1;
424     for (j1=0;j1<NMAX;j1++)
425       box[i1][j1]= -1;
426   }
427   time(&lasttime);
428   lnorm=0;
429   
430   for (n=1;n<nmax;n++) {
431     smaller=0;
432     sn=scr[n-1];
433     if (DIM > 1) {
434       x=(long)(series[0][sn]*epsinv)&imax;
435       y=(long)(series[1][sn]*epsinv)&imax;
436     }
437     else {
438       x=(long)(series[0][sn]*epsinv)&imax;
439       y=(long)(series[0][sn+DELAY]*epsinv)&imax;
440     }
441     list[sn]=box[x][y];
442     box[x][y]=sn;
443     listc1[sn]=boxc1[x];
444     boxc1[x]=sn;
445
446     i=imin;
447     while (found[maxembed][i] >= MAXFOUND) {
448       smaller=1;
449       if (++i > howoften1)
450         break;
451     }
452     if (smaller) {
453       imin=i;
454       if (imin <= howoften1) {
455         EPSMAX=epsm[imin];
456         epsinv=1.0/EPSMAX;
457         for (i1=0;i1<NMAX;i1++) {
458           boxc1[i1]= -1;
459           for (j1=0;j1<NMAX;j1++)
460             box[i1][j1]= -1;
461         }
462         for (i1=0;i1<n;i1++) {
463           sn=scr[i1];
464           if (DIM > 1) {
465             x=(long)(series[0][sn]*epsinv)&imax;
466             y=(long)(series[1][sn]*epsinv)&imax;
467           }
468           else {
469             x=(long)(series[0][sn]*epsinv)&imax;
470             y=(long)(series[0][sn+DELAY]*epsinv)&imax;
471           }
472           list[sn]=box[x][y];
473           box[x][y]=sn;
474           listc1[sn]=boxc1[x];
475           boxc1[x]=sn;
476         }
477       }
478     }
479
480     if (imin <= howoften1) {
481       lnorm=n;
482       if (MINDIST > 0) {
483         sn=scr[n];
484         n1=(sn-(long)MINDIST>=0)?sn-(long)MINDIST:0;
485         n2=(sn+MINDIST<length-(EMBED-1)*DELAY)?sn+MINDIST:
486           length-(EMBED-1)*DELAY-1;
487         for (i1=n1;i1<=n2;i1++)
488           if ((oscr[i1] < n))
489             lnorm--;
490       }
491       
492       if (EMBED*DIM > 1)
493         make_c2_dim(n);
494       make_c2_1(n);
495       for (i=imin;i<HOWOFTEN;i++)
496         norm[i] += (double)(lnorm);
497     }
498     
499     if (((time(&mytime)-lasttime) > WHEN) || (n == (nmax-1)) || 
500         (imin > howoften1)) {
501       time(&lasttime);
502       fstat=fopen(outstat,"w");
503       if (verbosity&VER_USR1)
504         fprintf(stderr,"Opened %s for writing\n",outstat);
505       fprintf(fstat,"Center points treated so far= %ld\n",n);
506       fprintf(fstat,"Maximal epsilon in the moment= %e\n",epsm[imin]);
507       fclose(fstat);
508       fout=fopen(outc1,"w");
509       if (verbosity&VER_USR1)
510         fprintf(stderr,"Opened %s for writing\n",outc1);
511       fprintf(fout,"#center= %ld\n",n);
512       for (i=0;i<EMBED*DIM;i++) {
513         fprintf(fout,"#dim= %ld\n",i+1);
514         eps=EPSMAX1*epsfactor;
515         for (j=0;j<HOWOFTEN;j++) {
516           eps /= epsfactor;
517           if (norm[j] > 0.0)
518             fprintf(fout,"%e %e\n",eps,found[i][j]/norm[j]);
519         }
520         fprintf(fout,"\n\n");
521       }
522       fclose(fout);
523       fout=fopen(outh1,"w");
524       if (verbosity&VER_USR1)
525         fprintf(stderr,"Opened %s for writing\n",outh1);
526       fprintf(fout,"#center= %ld\n",n);
527       fprintf(fout,"#dim= 1\n");
528       eps=EPSMAX1*epsfactor;
529       for (j=0;j<HOWOFTEN;j++) {
530         eps /= epsfactor;
531         if (found[0][j] > 0.0)
532           fprintf(fout,"%e %e\n",eps,-log(found[0][j]/norm[j]));
533       }
534       fprintf(fout,"\n\n");
535       for (i=1;i<DIM*EMBED;i++) {
536         fprintf(fout,"#dim= %ld\n",i+1);
537         eps=EPSMAX1*epsfactor;
538         for (j=0;j<HOWOFTEN;j++) {
539           eps /= epsfactor;
540           if ((found[i-1][j] > 0.0) && (found[i][j] > 0.0))
541             fprintf(fout,"%e %e\n",eps,log(found[i-1][j]/found[i][j]));
542         }
543         fprintf(fout,"\n\n");
544       }
545       fclose(fout);
546       fout=fopen(outd1,"w");
547       if (verbosity&VER_USR1)
548         fprintf(stderr,"Opened %s for writing\n",outd1);
549       fprintf(fout,"#center= %ld\n",n);
550       for (i=0;i<DIM*EMBED;i++) {
551         fprintf(fout,"#dim= %ld\n",i+1);
552         eps=EPSMAX1;
553         for (j=1;j<HOWOFTEN;j++) {
554           eps /= epsfactor;
555           if ((found[i][j] > 0.0) && (found[i][j-1] > 0.0))
556             fprintf(fout,"%e %e\n",eps,log(found[i][j-1]/found[i][j]
557                                            /norm[j-1]*norm[j])/lnfac);
558         }
559         fprintf(fout,"\n\n");
560       }
561       fclose(fout);
562       if (imin > howoften1)
563         exit(0);
564     }
565   }
566
567   if (infile != NULL)
568     free(infile);
569   free(outd1);
570   free(outh1);
571   free(outc1);
572   free(outstat);
573   free(list);
574   free(listc1);
575   free(scr);
576   free(oscr);
577   free(norm);
578   free(epsm);
579   for (i=0;i<EMBED*DIM;i++)
580     free(found[i]);
581   free(found);
582   for (i=0;i<DIM;i++)
583     free(series[i]);
584   free(series);
585
586   return 0;
587 }