IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / false_nearest.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: Dec 10, 2005 */
21 /*Changes:
22   12/10/05: It's multivariate now
23   12/16/05: Scaled <eps> and sigma(eps)
24 */
25
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <string.h>
29 #include <limits.h>
30 #include <math.h>
31 #include "routines/tsa.h"
32
33 #define WID_STR "Determines the fraction of false nearest neighbors."
34
35 char *outfile=NULL;
36 char *infile=NULL;
37 char stdo=1,dimset=0;
38 char *column=NULL;
39 unsigned long length=ULONG_MAX,exclude=0,theiler=0;
40 unsigned int delay=1,maxdim=5,minemb=1;
41 unsigned int comp=1,maxemb=5;
42 unsigned int verbosity=0xff;
43 double rt=2.0;
44 double eps0=1.0e-5;
45 double **series;
46 double aveps,vareps;
47 double varianz;
48
49 #define BOX 1024
50 int ibox=BOX-1;
51 long **box,*list;
52 unsigned int *vcomp,*vemb;
53 unsigned long toolarge;
54
55 void show_options(char *progname)
56 {
57   what_i_do(progname,WID_STR);
58   fprintf(stderr," Usage: %s [options]\n",progname);
59   fprintf(stderr," Options:\n");
60   fprintf(stderr,"Everything not being a valid option will be interpreted"
61           " as a possible"
62           " datafile.\nIf no datafile is given stdin is read. Just - also"
63           " means stdin\n");
64   fprintf(stderr,"\t-l # of data [default: whole file]\n");
65   fprintf(stderr,"\t-x # of lines to ignore [default: 0]\n");
66   fprintf(stderr,"\t-c columns to read [default: 1]\n");
67   fprintf(stderr,"\t-m min. test embedding dimension [default: %u]\n",minemb);
68   fprintf(stderr,"\t-M # of components,max. emb. dim. [default: %u,%u]\n",
69           comp,maxemb);
70   fprintf(stderr,"\t-d delay [default: 1]\n");
71   fprintf(stderr,"\t-f escape factor [default: %.2lf]\n",rt);
72   fprintf(stderr,"\t-t theiler window [default: 0]\n");
73   fprintf(stderr,"\t-o output file [default: 'datafile'.fnn; without -o"
74           " stdout]\n");
75   fprintf(stderr,"\t-V verbosity level [default: 3]\n\t\t"
76           "0='only panic messages'\n\t\t"
77           "1='+ input/output messages'\n\t\t"
78           "2='+ information about the current state\n");
79   fprintf(stderr,"\t-h show these options\n");
80   exit(0);
81 }
82
83 void scan_options(int n,char **in)
84 {
85   char *out;
86
87   if ((out=check_option(in,n,'l','u')) != NULL)
88     sscanf(out,"%lu",&length);
89   if ((out=check_option(in,n,'x','u')) != NULL)
90     sscanf(out,"%lu",&exclude);
91   if ((out=check_option(in,n,'c','s')) != NULL)
92     column=out;
93   if ((out=check_option(in,n,'m','u')) != NULL)
94     sscanf(out,"%u",&minemb);
95   if ((out=check_option(in,n,'M','2')) != NULL) {
96     sscanf(out,"%u,%u",&comp,&maxemb);
97     maxdim=comp*(maxemb+1);
98     dimset=1;
99   }
100   if ((out=check_option(in,n,'d','u')) != NULL)
101     sscanf(out,"%u",&delay);
102   if ((out=check_option(in,n,'f','f')) != NULL)
103     sscanf(out,"%lf",&rt);
104   if ((out=check_option(in,n,'t','u')) != NULL)
105     sscanf(out,"%lu",&theiler);
106   if ((out=check_option(in,n,'V','u')) != NULL)
107     sscanf(out,"%u",&verbosity);
108   if ((out=check_option(in,n,'o','o')) != NULL) {
109     stdo=0;
110     if (strlen(out) > 0)
111       outfile=out;
112   }
113 }
114
115 void mmb(unsigned int hdim,unsigned int hemb,double eps)
116 {
117   unsigned long i;
118   long x,y;
119
120   for (x=0;x<BOX;x++)
121     for (y=0;y<BOX;y++)
122       box[x][y] = -1;
123
124   for (i=0;i<length-(maxemb+1)*delay;i++) {
125     x=(long)(series[0][i]/eps)&ibox;
126     y=(long)(series[hdim][i+hemb]/eps)&ibox;
127     list[i]=box[x][y];
128     box[x][y]=i;
129   }
130 }
131
132 char find_nearest(long n,unsigned int dim,double eps)
133 {
134   long x,y,x1,x2,y1,i,i1,ic,ie;
135   long element,which= -1;
136   double dx,maxdx,mindx=1.1,hfactor,factor;
137
138   ic=vcomp[dim];
139   ie=vemb[dim];
140   x=(long)(series[0][n]/eps)&ibox;
141   y=(long)(series[ic][n+ie]/eps)&ibox;
142   
143   for (x1=x-1;x1<=x+1;x1++) {
144     x2=x1&ibox;
145     for (y1=y-1;y1<=y+1;y1++) {
146       element=box[x2][y1&ibox];
147       while (element != -1) {
148         if (labs(element-n) > theiler) {
149           maxdx=fabs(series[0][n]-series[0][element]);
150           for (i=1;i<=dim;i++) {
151             ic=vcomp[i];
152             i1=vemb[i];
153             dx=fabs(series[ic][n+i1]-series[ic][element+i1]);
154             if (dx > maxdx)
155               maxdx=dx;
156           }
157           if ((maxdx < mindx) && (maxdx > 0.0)) {
158             which=element;
159             mindx=maxdx;
160           }
161         }
162         element=list[element];
163       }
164     }
165   }
166
167   if ((which != -1) && (mindx <= eps) && (mindx <= varianz/rt)) {
168     aveps += mindx;
169     vareps += mindx*mindx;
170     factor=0.0;
171     for (i=1;i<=comp;i++) {
172       ic=vcomp[dim+i];
173       ie=vemb[dim+i];
174       hfactor=fabs(series[ic][n+ie]-series[ic][which+ie])/mindx;
175       if (hfactor > factor) 
176         factor=hfactor;
177     }
178     if (factor > rt)
179       toolarge++;
180     return 1;
181   }
182   return 0;
183 }
184
185 int main(int argc,char **argv)
186 {
187   char stdi=0;
188   FILE *file=NULL;
189   double min,inter=0.0,ind_inter,epsilon,av,ind_var;
190   char *nearest,alldone;
191   long i;
192   unsigned int dim,emb;
193   unsigned long donesofar;
194
195   if (scan_help(argc,argv))
196     show_options(argv[0]);
197   
198   scan_options(argc,argv);
199 #ifndef OMIT_WHAT_I_DO
200   if (verbosity&VER_INPUT)
201     what_i_do(argv[0],WID_STR);
202 #endif
203
204   infile=search_datafile(argc,argv,NULL,verbosity);
205   if (infile == NULL)
206     stdi=1;
207
208   if (outfile == NULL) {
209     if (!stdi) {
210       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
211       strcpy(outfile,infile);
212       strcat(outfile,".fnn");
213     }
214     else {
215       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
216       strcpy(outfile,"stdin.fnn");
217     }
218   }
219   if (!stdo)
220     test_outfile(outfile);
221
222   if (column == NULL)
223     series=(double**)get_multi_series(infile,&length,exclude,&comp,"",dimset,
224                                       verbosity);
225   else
226     series=(double**)get_multi_series(infile,&length,exclude,&comp,column,
227                                       dimset,verbosity);
228
229   for (i=0;i<comp;i++) {
230     rescale_data(series[i],length,&min,&ind_inter);
231     variance(series[i],length,&av,&ind_var);
232     if (i == 0) {
233       varianz=ind_var;
234       inter=ind_inter;
235     }
236     else {
237       varianz=(varianz>ind_var)?ind_var:varianz;
238       inter=(inter<ind_inter)?ind_inter:inter;
239     }
240   }
241
242   check_alloc(list=(long*)malloc(sizeof(long)*length));
243   check_alloc(nearest=(char*)malloc(length));
244   check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
245   for (i=0;i<BOX;i++)
246     check_alloc(box[i]=(long*)malloc(sizeof(long)*BOX));
247
248   if (!stdo) {
249     file=fopen(outfile,"w");
250     if (verbosity&VER_INPUT)
251       fprintf(stderr,"Opened %s for writing\n",outfile);
252   }
253   else {
254     if (verbosity&VER_INPUT)
255       fprintf(stderr,"Writing to stdout\n");
256   }
257   check_alloc(vcomp=(unsigned int*)malloc(sizeof(int)*(maxdim)));
258   check_alloc(vemb=(unsigned int*)malloc(sizeof(int)*(maxdim)));
259   for (i=0;i<maxdim;i++) {
260     if (comp == 1) {
261       vcomp[i]=0;
262       vemb[i]=i;
263     }
264     else {
265       vcomp[i]=i%comp;
266       vemb[i]=(i/comp)*delay;
267     }
268   }
269   for (emb=minemb;emb<=maxemb;emb++) {
270     dim=emb*comp-1;
271     epsilon=eps0;
272     toolarge=0;
273     alldone=0;
274     donesofar=0;
275     aveps=0.0;
276     vareps=0.0;
277     for (i=0;i<length;i++)
278       nearest[i]=0;
279     if (verbosity&VER_USR1)
280       fprintf(stderr,"Start for dimension=%u\n",dim+1);
281     while (!alldone && (epsilon < 2.*varianz/rt)) {
282       alldone=1;
283       mmb(vcomp[dim],vemb[dim],epsilon);
284       for (i=0;i<length-maxemb*delay;i++)
285         if (!nearest[i]) {
286           nearest[i]=find_nearest(i,dim,epsilon);
287           alldone &= nearest[i];
288           donesofar += (unsigned long)nearest[i];
289         }
290       if (verbosity&VER_USR1)
291         fprintf(stderr,"Found %lu up to epsilon=%e\n",donesofar,epsilon*inter);
292       epsilon*=sqrt(2.0);
293       if (!donesofar)
294         eps0=epsilon;
295     }
296     if (donesofar == 0) {
297       fprintf(stderr,"Not enough points found!\n");
298       exit(FALSE_NEAREST_NOT_ENOUGH_POINTS);
299     }
300     aveps *= (1./(double)donesofar);
301     vareps *= (1./(double)donesofar);
302     if (stdo) {
303       fprintf(stdout,"%u %e %e %e\n",dim+1,(double)toolarge/(double)donesofar,
304               aveps*inter,sqrt(vareps)*inter);
305       fflush(stdout);
306     }
307     else {
308       fprintf(file,"%u %e %e %e\n",dim+1,(double)toolarge/(double)donesofar,
309               aveps*inter,sqrt(vareps)*inter);
310       fflush(file);
311     }
312   }
313   if (!stdo)
314     fclose(file);
315
316   if (infile != NULL)
317     free(infile);
318   if (outfile != NULL)
319     free(outfile);
320   free(series);
321   free(list);
322   free(nearest);
323   for (i=0;i<BOX;i++)
324     free(box[i]);
325   free(box);
326
327   return 0;
328 }