IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / lzo-run.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: Feb 19, 2007 */
21 /* Changes:
22      2/19/2007:  Changed name and default for noise  
23      10/26/2006: Add seed option
24  */
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <string.h>
28 #include <limits.h>
29 #include <time.h>
30 #include <math.h>
31 #include "routines/tsa.h"
32
33 #define WID_STR "Makes a local zeroth order forecast for multivariate data\n\
34 and iterates a trajectory"
35
36 #define NMAX 128
37
38 char onscreen=1,epsset=0,*outfile=NULL,setsort=1,setnoise=0;
39 char *infile=NULL;
40 unsigned int nmax=(NMAX-1);
41 unsigned int verbosity=0xff;
42 long **box,*list,*found;
43 double **series,**cast,*abstand,*var;
44 double epsilon;
45
46 unsigned int embed=2,dim=1,dim1,DELAY=1;
47 char *column=NULL,dimset=0;
48 unsigned int MINN=50;
49 unsigned int **indexes;
50 unsigned long LENGTH=ULONG_MAX,FLENGTH=1000,exclude=0;
51 unsigned long seed=0x9074325L;
52 double EPS0=1.e-3,EPSF=1.2,Q=10.0;
53
54 double **mat,*vec,*hsum,*newav;
55
56 void show_options(char *progname)
57 {
58   what_i_do(progname,WID_STR);
59   fprintf(stderr," Usage: %s [Options]\n",progname);
60   fprintf(stderr," Options:\n");
61   fprintf(stderr,"Everything not being a valid option will be interpreted"
62           " as a possible"
63           " datafile.\nIf no datafile is given stdin is read. Just - also"
64           " means stdin\n");
65   fprintf(stderr,"\t-l # of data to be used [default whole file]\n");
66   fprintf(stderr,"\t-x # of lines to be ignored [default 0]\n");
67   fprintf(stderr,"\t-c column [default 1,...,# of components]\n");
68   fprintf(stderr,"\t-m #of components,embedding dimension [default 1,2]\n");
69   fprintf(stderr,"\t-d delay for the embedding [default 1]\n");
70   fprintf(stderr,"\t-L # of iterations [default 1000]\n");
71   fprintf(stderr,"\t-k # of neighbors  [default %u]\n",MINN);
72   fprintf(stderr,"\t-K fix # of neighbors  [default no]\n");
73   fprintf(stderr,"\t-%% # variance of noise [default %3.1lf]\n",Q);
74   fprintf(stderr,"\t-I seed for the rnd-generator (If seed=0, the time\n"
75           "\t\tcommand is used to set the seed) [Default: fixed]\n");
76   fprintf(stderr,"\t-r size of initial neighborhood ["
77           " default (data interval)/1000]\n");
78   fprintf(stderr,"\t-f factor to increase size [default 1.2]\n");
79   fprintf(stderr,"\t-o output file [default 'datafile'.lzr;"
80           " no -o means write to stdout]\n");
81   fprintf(stderr,"\t-V verbosity level [default 1]\n\t\t"
82           "0='only panic messages'\n\t\t"
83           "1='+ input/output messages'\n");
84   fprintf(stderr,"\t-h  show these options\n");
85   exit(0);
86 }
87
88 void scan_options(int n,char **in)
89 {
90   char *out;
91
92   if ((out=check_option(in,n,'l','u')) != NULL)
93     sscanf(out,"%lu",&LENGTH);
94   if ((out=check_option(in,n,'x','u')) != NULL)
95     sscanf(out,"%lu",&exclude);
96   if ((out=check_option(in,n,'c','s')) != NULL) {
97     column=out;
98     dimset=1;
99   }
100   if ((out=check_option(in,n,'m','2')) != NULL)
101     sscanf(out,"%u,%u",&dim,&embed);
102   if ((out=check_option(in,n,'d','u')) != NULL)
103     sscanf(out,"%u",&DELAY);
104   if ((out=check_option(in,n,'L','u')) != NULL)
105     sscanf(out,"%lu",&FLENGTH);
106   if ((out=check_option(in,n,'k','u')) != NULL)
107     sscanf(out,"%u",&MINN);
108   if ((out=check_option(in,n,'K','n')) != NULL)
109     setsort=1;
110   if ((out=check_option(in,n,'I','u')) != NULL) {
111     sscanf(out,"%lu",&seed);
112     if (seed == 0)
113       seed=(unsigned long)time((time_t*)&seed);
114   }
115   if ((out=check_option(in,n,'V','u')) != NULL)
116     sscanf(out,"%u",&verbosity);
117   if ((out=check_option(in,n,'r','f')) != NULL) {
118     epsset=1;
119     sscanf(out,"%lf",&EPS0);
120   }
121   if ((out=check_option(in,n,'f','f')) != NULL)
122     sscanf(out,"%lf",&EPSF);
123   if ((out=check_option(in,n,'%','f')) != NULL) {
124     sscanf(out,"%lf",&Q);
125     if (Q>0.0)
126       setnoise=1;
127   }
128   if ((out=check_option(in,n,'o','o')) != NULL) {
129     onscreen=0;
130     if (strlen(out) > 0)
131       outfile=out;
132   }
133 }
134
135 void sort(unsigned long nfound)
136 {
137   double dx,dswap;
138   int i,j,k,hf,iswap,hdim;
139
140   hdim=(embed-1)*DELAY;
141
142   for (i=0;i<nfound;i++) {
143     hf=found[i];
144     abstand[i]=0.0;
145     for (j=0;j<dim;j++) {
146       for (k=0;k<=hdim;k += DELAY) {
147         dx=fabs(series[j][hf-k]-cast[hdim-k][j]);
148         if (dx > abstand[i]) abstand[i]=dx;
149       }
150     }
151   }
152
153   for (i=0;i<MINN;i++)
154     for (j=i+1;j<nfound;j++)
155       if (abstand[j]<abstand[i]) {
156         dswap=abstand[i];
157         abstand[i]=abstand[j];
158         abstand[j]=dswap;
159         iswap=found[i];
160         found[i]=found[j];
161         found[j]=iswap;
162       }
163 }
164
165 void put_in_boxes(void)
166 {
167   int i,j,n;
168   static int hdim;
169   double epsinv;
170
171   hdim=(embed-1)*DELAY;
172   epsinv=1.0/epsilon;
173   for (i=0;i<NMAX;i++)
174     for (j=0;j<NMAX;j++)
175       box[i][j]= -1;
176
177   for (n=hdim;n<LENGTH-1;n++) {
178     i=(int)(series[0][n]*epsinv)&nmax;
179     j=(int)(series[dim1][n-hdim]*epsinv)&nmax;
180     list[n]=box[i][j];
181     box[i][j]=n;
182   }
183 }
184
185 unsigned int hfind_neighbors(void)
186 {
187   char toolarge;
188   int i,j,i1,i2,j1,l,hc,hd,element;
189   static int hdim;
190   unsigned nfound=0;
191   double max,dx,epsinv;
192
193   hdim=(embed-1)*DELAY;
194   epsinv=1.0/epsilon;
195   i=(int)(cast[hdim][0]*epsinv)&nmax;
196   j=(int)(cast[0][dim1]*epsinv)&nmax;
197   
198   for (i1=i-1;i1<=i+1;i1++) {
199     i2=i1&nmax;
200     for (j1=j-1;j1<=j+1;j1++) {
201       element=box[i2][j1&nmax];
202       while (element != -1) {
203         max=0.0;
204         toolarge=0;
205         for (l=0;l<dim*embed;l++) {
206           hc=indexes[0][l];
207           hd=indexes[1][l];
208           dx=fabs(series[hc][element-hd]-cast[hdim-hd][hc]);
209           max=(dx>max) ? dx : max;
210           if (max > epsilon) {
211             toolarge=1;
212             break;
213           }
214         }
215         if (max <= epsilon)
216           found[nfound++]=element;
217         element=list[element];
218       }
219     }
220   }
221   return nfound;
222 }
223
224 void make_zeroth(int number,double *newcast)
225 {
226   long d,i;
227   double *sd;
228   
229   for (d=0;d<dim;d++) {
230     newcast[d]=0.0;
231     sd=series[d]+1;
232     for (i=0;i<number;i++)
233       newcast[d] += sd[found[i]];
234     newcast[d] /= (double)number;
235   }
236
237   if (setnoise) {
238     for (d=0;d<dim;d++)
239       newcast[d] += gaussian(var[d]*Q);
240   }
241 }
242
243 int main(int argc,char **argv)
244 {
245   char stdi=0,done;
246   long i,j,hdim,actfound;
247   unsigned long count=1;
248   double *swap,*newcast,maxinterval,*min,*interval,dummy,epsilon0;
249   FILE *file=NULL;
250   
251   if (scan_help(argc,argv))
252     show_options(argv[0]);
253   
254   scan_options(argc,argv);
255 #ifndef OMIT_WHAT_I_DO
256   if (verbosity&VER_INPUT)
257     what_i_do(argv[0],WID_STR);
258 #endif
259   
260   infile=search_datafile(argc,argv,NULL,verbosity);
261   if (infile == NULL)
262     stdi=1;
263   
264   if (outfile == NULL) {
265     if (!stdi) {
266       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
267       strcpy(outfile,infile);
268       strcat(outfile,".lzr");
269     }
270     else {
271       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
272       strcpy(outfile,"stdin.lzr");
273     }
274   }
275   if (!onscreen)
276     test_outfile(outfile);
277   
278   hdim=(embed-1)*DELAY+1;
279   if (column == NULL)
280     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,"",dimset,
281                                       verbosity);
282   else
283     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,column,
284                                       dimset,verbosity);
285
286   dim1=dim-1;
287
288   check_alloc(min=(double*)malloc(sizeof(double)*dim));
289   check_alloc(interval=(double*)malloc(sizeof(double)*dim));
290   check_alloc(var=(double*)malloc(sizeof(double)*dim));
291
292   maxinterval=0.0;
293
294   for (i=0;i<dim;i++) {
295     rescale_data(series[i],LENGTH,&min[i],&interval[i]);
296     variance(series[i],LENGTH,&dummy,&var[i]);
297     if (interval[i] > maxinterval)
298       maxinterval=interval[i];
299   }
300
301   if (epsset)
302     EPS0 /= maxinterval;
303     
304   check_alloc(cast=(double**)malloc(sizeof(double*)*hdim));
305   for (i=0;i<hdim;i++)
306     check_alloc(cast[i]=(double*)malloc(sizeof(double)*dim));
307   check_alloc(newcast=(double*)malloc(sizeof(double)*dim));
308   check_alloc(newav=(double*)malloc(sizeof(double)*dim));
309     
310   check_alloc(list=(long*)malloc(sizeof(long)*LENGTH));
311   check_alloc(found=(long*)malloc(sizeof(long)*LENGTH));
312   check_alloc(abstand=(double*)malloc(sizeof(double)*LENGTH));
313   check_alloc(box=(long**)malloc(sizeof(long*)*NMAX));
314   for (i=0;i<NMAX;i++)
315     check_alloc(box[i]=(long*)malloc(sizeof(long)*NMAX));
316   
317   check_alloc(vec=(double*)malloc(sizeof(double)*dim));
318   check_alloc(hsum=(double*)malloc(sizeof(double)*dim));
319   check_alloc(mat=(double**)malloc(sizeof(double*)*dim));
320   for (i=0;i<dim;i++) {
321     check_alloc(mat[i]=(double*)malloc(sizeof(double)*dim));
322   }
323
324   for (j=0;j<dim;j++)
325     for (i=0;i<hdim;i++)
326       cast[i][j]=series[j][LENGTH-hdim+i];
327
328   indexes=make_multi_index(dim,embed,DELAY);
329   
330   if (!onscreen) {
331     file=fopen(outfile,"w");
332     if (verbosity&VER_INPUT)
333       fprintf(stderr,"Opened %s for writing\n",outfile);
334   }
335   else {
336     if (verbosity&VER_INPUT)
337       fprintf(stderr,"Writing to stdout\n");
338   }
339
340   rnd_init(seed);
341
342   epsilon0=EPS0/EPSF;
343
344   if (setnoise) 
345     Q /= 100.0;
346
347   for (i=0;i<FLENGTH;i++) {
348     done=0;
349     if (setsort)
350       epsilon= epsilon0/((double)count*EPSF);
351     else
352       epsilon=epsilon0;
353     while (!done) {
354       epsilon*=EPSF;
355       put_in_boxes();
356       actfound=hfind_neighbors();
357       if (actfound >= MINN) {
358         if (setsort) {
359           epsilon0 += epsilon;
360           count++;
361           sort(actfound);
362           actfound=MINN;
363         }
364         make_zeroth(actfound,newcast);
365         if (onscreen) {
366           for (j=0;j<dim-1;j++)
367             printf("%e ",newcast[j]*interval[j]+min[j]);
368           printf("%e\n",newcast[dim-1]*interval[dim-1]+min[dim-1]);
369           fflush(stdout);
370         }
371         else {
372           for (j=0;j<dim-1;j++)
373             fprintf(file,"%e ",newcast[j]*interval[j]+min[j]);
374           fprintf(file,"%e\n",newcast[dim-1]*interval[dim-1]+min[dim-1]);
375           fflush(file);
376         }
377         done=1;
378         swap=cast[0];
379         for (j=0;j<hdim-1;j++)
380           cast[j]=cast[j+1];
381         cast[hdim-1]=swap;
382         for (j=0;j<dim;j++)
383           cast[hdim-1][j]=newcast[j];
384       }
385     }
386   }
387   if (!onscreen)
388     fclose(file);
389   
390   if (outfile != NULL)
391     free(outfile);
392   for (i=0;i<dim;i++)
393     free(mat[i]);
394   free(mat);
395   for (i=0;i<hdim;i++)
396     free(cast[i]);
397   free(cast);
398   free(newcast);
399   free(found);
400   free(list);
401   for (i=0;i<NMAX;i++)
402     free(box[i]);
403   free(box);
404   free(vec);
405   free(newav);
406   for (i=0;i<dim;i++)
407     free(series[i]);
408   free(series);
409
410   return 0;
411 }