IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / nrlazy.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: Nov 30, 2000 */
21 /*Changes:
22   12/11/05: Going multivariate
23 */
24
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <string.h>
28 #include <math.h>
29 #include <limits.h>
30 #include "routines/tsa.h"
31
32 #define WID_STR "Performs simple noise reduction."
33
34 #define BOX (unsigned int)512
35
36 unsigned long length=ULONG_MAX,exclude=0;
37 unsigned int comp=1,embed=5,delay=1,iterations=1,alldim;
38 unsigned int verbosity=0x3;
39 char *column=NULL;
40 double eps=1.0e-3,epsvar;
41
42 char *outfile=NULL,epsset=0,stdo=1,epsvarset=0;
43 char *infile=NULL;
44 double **series,**corr,*interval,*min,*hcor;
45 long **box,*list,**nf;
46 unsigned int **indexes;
47 char dimset=0;
48
49 void show_options(char *progname)
50 {
51   what_i_do(progname,WID_STR);
52   fprintf(stderr," Usage: %s [Options]\n",progname);
53   fprintf(stderr," Options:\n");
54   fprintf(stderr,"Everything not being a valid option will be interpreted"
55           " as a possible"
56           " datafile.\nIf no datafile is given stdin is read. Just - also"
57           " means stdin\n");
58   fprintf(stderr,"\t-l # of data to use [default: whole file]\n");
59   fprintf(stderr,"\t-x # of lines to be ignored [default: 0]\n");
60   fprintf(stderr,"\t-c column to read [default: 1]\n");
61   fprintf(stderr,"\t-m no. of comp.,embedding dim. [default: %u,%u]\n",
62           comp,embed);
63   fprintf(stderr,"\t-d delay [default: 1]\n");
64   fprintf(stderr,"\t-i iterations [default: 1]\n");
65   fprintf(stderr,"\t-r neighborhoud size [default: (interval of data)/1000]\n");
66   fprintf(stderr,"\t-v neighborhoud size (in units of the std. dev. of the "
67           "data \n\t\t(overwrites -r) [default: not set]\n");
68   fprintf(stderr,"\t-o output file name [Default: 'datafile'.laz.n,"
69           "\n\t\twhere n is the number of the last iteration,"
70           "\n\t\twithout -o the last iteration is written to stdout.]\n");
71   fprintf(stderr,"\t-V verbosity level [Default: 3]\n\t\t"
72           "0='only panic messages'\n\t\t"
73           "1='+ input/output messages'\n\t\t"
74           "2='+ write output of all iterations to files'\n\t\t"
75           "4='+ write the number of neighbors found for each point\n");
76   fprintf(stderr,"\t-h show these options\n");
77   exit(0);
78 }
79
80 void scan_options(int n,char **in)
81 {
82   char *out;
83
84   if ((out=check_option(in,n,'l','u')) != NULL)
85     sscanf(out,"%lu",&length);
86   if ((out=check_option(in,n,'x','u')) != NULL)
87     sscanf(out,"%lu",&exclude);
88   if ((out=check_option(in,n,'c','c')) != NULL) {
89     column=out;
90     dimset=1;
91   }
92   if ((out=check_option(in,n,'m','2')) != NULL)
93     sscanf(out,"%u,%u",&comp,&embed);
94   if ((out=check_option(in,n,'d','u')) != NULL)
95     sscanf(out,"%u",&delay);
96   if ((out=check_option(in,n,'i','u')) != NULL)
97     sscanf(out,"%u",&iterations);
98   if ((out=check_option(in,n,'r','f')) != NULL) {
99     epsset=1;
100     sscanf(out,"%lf",&eps);
101   }
102   if ((out=check_option(in,n,'V','u')) != NULL)
103     sscanf(out,"%u",&verbosity);
104   if ((out=check_option(in,n,'v','f')) != NULL) {
105     epsvarset=1;
106     sscanf(out,"%lf",&epsvar);
107   }
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 unsigned int correct(unsigned long n)
116 {
117   int i,i1,i2,j,j1,k;
118   int ibox=BOX-1;
119   unsigned int hdel,hcomp;
120   double epsinv,dx;
121   long element,nfound=0;
122
123   epsinv=1./eps;
124
125   for (i=0;i<alldim;i++)
126     hcor[i]=0.0;
127
128   i=(int)(series[0][n]*epsinv)&ibox;
129   j=(int)(series[comp-1][n-(embed-1)*delay]*epsinv)&ibox;
130   
131   for (i1=i-1;i1<=i+1;i1++) {
132     i2=i1&ibox;
133     for (j1=j-1;j1<=j+1;j1++) {
134       element=box[i2][j1&ibox];
135       while (element != -1) {
136         for (k=0;k<alldim;k++) {
137           hcomp=indexes[0][k];
138           hdel=indexes[1][k];
139           dx=fabs(series[hcomp][n-hdel]-series[hcomp][element-hdel]);
140           if (dx > eps)
141             break;
142         }
143         if (k == alldim) {
144           nfound++;
145           for (k=0;k<alldim;k++) {
146             hcomp=indexes[0][k];
147             hdel=indexes[1][k];
148             hcor[k] += series[hcomp][element-hdel];
149           }
150         }
151         element=list[element];
152       }
153     }
154   }
155   for (k=0;k<alldim;k++) {
156     hcomp=indexes[0][k];
157     hdel=indexes[1][k];
158     corr[hcomp][n-hdel] += hcor[k]/nfound;
159     nf[hcomp][n-hdel]++;
160   }
161
162   return nfound;
163 }
164
165 int main(int argc,char **argv)
166 {
167   char *ofname;
168   char stdi=0;
169   int iter;
170   unsigned int *nmf;
171   unsigned long n,i;
172   double dav,dvar,maxinterval,maxdvar;
173   FILE *file=NULL;
174
175   if (scan_help(argc,argv))
176     show_options(argv[0]);
177   
178   scan_options(argc,argv);
179 #ifndef OMIT_WHAT_I_DO
180   if (verbosity&VER_INPUT)
181     what_i_do(argv[0],WID_STR);
182 #endif
183
184   infile=search_datafile(argc,argv,NULL,verbosity);
185   if (infile == NULL)
186     stdi=1;
187
188   if (outfile == NULL) {
189     if (!stdi) {
190       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
191       check_alloc(ofname=(char*)calloc(strlen(infile)+9,(size_t)1));
192       sprintf(outfile,"%s.laz",infile);
193     }
194     else {
195       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
196       check_alloc(ofname=(char*)calloc((size_t)14,(size_t)1));
197       sprintf(outfile,"stdin.laz");
198     }
199   }
200   else
201     check_alloc(ofname=(char*)calloc(strlen(outfile)+10,(size_t)1));
202
203
204   if (column == NULL)
205     series=(double**)get_multi_series(infile,&length,exclude,&comp,"",dimset,
206                                       verbosity);
207   else
208     series=(double**)get_multi_series(infile,&length,exclude,&comp,column,
209                                       dimset,verbosity);
210
211   check_alloc(interval=(double*)malloc(sizeof(double)*comp));
212   check_alloc(min=(double*)malloc(sizeof(double)*comp));
213
214   maxinterval=maxdvar=0.0;
215   for (i=0;i<comp;i++) {
216     rescale_data(series[i],length,&min[i],&interval[i]);
217     if (interval[i] > maxinterval) maxinterval=interval[i];
218     variance(series[i],length,&dav,&dvar);
219     if (dvar > maxdvar)  maxdvar=dvar;
220   }
221   alldim=comp*embed;
222
223   check_alloc(nmf=(unsigned int*)malloc(sizeof(int)*length));
224   check_alloc(list=(long*)malloc(sizeof(long)*length));
225   check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
226   for (n=0;n<BOX;n++)
227     check_alloc(box[n]=(long*)malloc(sizeof(long)*BOX));
228
229   check_alloc(nf=(long**)malloc(sizeof(long*)*comp));
230   check_alloc(corr=(double**)malloc(sizeof(double*)*comp));
231   for (i=0;i<comp;i++) {
232     check_alloc(nf[i]=(long*)malloc(sizeof(long)*length));
233     check_alloc(corr[i]=(double*)malloc(sizeof(double)*length));
234   }
235
236   indexes=make_multi_index(comp,embed,delay);
237
238   if (epsset)
239     eps/=maxinterval;
240   else
241     eps=1.0/1000.;
242
243   if (epsvarset)
244     eps=epsvar*maxdvar;
245
246   for (iter=1;iter<=iterations;iter++) {
247     make_multi_box2(series,box,list,length,BOX,comp,embed,delay,eps);
248     for (n=0;n<length;n++) {
249       for (i=0;i<comp;i++) {
250         corr[i][n]=0.0;
251         nf[i][n]=0;
252       }
253       nmf[n]=1;
254     }
255     
256     check_alloc(hcor=(double*)malloc(sizeof(double)*alldim));
257     for (n=(embed-1)*delay;n<length;n++)
258       nmf[n]=correct(n);
259     free(hcor);
260     
261     for (n=0;n<length;n++)
262       for (i=0;i<comp;i++)
263         if (nf[i][n])
264           series[i][n]=corr[i][n]/nf[i][n];
265
266     if ((verbosity&VER_USR1) && (iter < iterations)) {
267       sprintf(ofname,"%s.%d",outfile,iter);
268       test_outfile(ofname);
269       file=fopen(ofname,"w");
270       if (verbosity&VER_INPUT)
271         fprintf(stderr,"Opened %s for writing\n",ofname);
272       if (stdo && (iter == iterations)) {
273         if (verbosity&VER_INPUT)
274           fprintf(stderr,"Writing to stdout\n");
275       }
276       for (n=0;n<length;n++) {
277         if (stdo && (iter == iterations)) {
278           if (verbosity&VER_USR2) {
279             for (i=0;i<comp;i++) 
280               fprintf(stdout,"%e ",series[i][n]*interval[i]+min[i]);
281             fprintf(stdout,"%u\n",nmf[n]);
282           }
283           else {
284             fprintf(stdout,"%e",series[0][n]*interval[0]+min[0]);
285             for (i=1;i<comp;i++)
286               fprintf(stdout,"%e ",series[i][n]*interval[i]+min[i]);
287             fprintf(stdout,"\n");
288           }
289         }
290         if (verbosity&VER_USR2) {
291           for (i=0;i<comp;i++) 
292             fprintf(file,"%e ",series[i][n]*interval[i]+min[i]);
293           fprintf(file,"%u\n",nmf[n]);
294         }
295         else {
296           fprintf(file,"%e",series[0][n]*interval[0]+min[0]);
297           for (i=1;i<comp;i++)
298             fprintf(file," %e",series[i][n]*interval[i]+min[i]);
299           fprintf(file,"\n");
300         }
301       }
302       fclose(file);
303     }
304     if (iter == iterations) {
305       if (!stdo || (verbosity&VER_USR1)) {
306         sprintf(ofname,"%s.%d",outfile,iter);
307         test_outfile(ofname);
308         file=fopen(ofname,"w");
309         if (verbosity&VER_INPUT)
310           fprintf(stderr,"Opened %s for writing\n",ofname);
311         if (stdo && (iter == iterations)) {
312           if (verbosity&VER_INPUT)
313             fprintf(stderr,"Writing to stdout\n");
314         }
315       }
316       for (n=0;n<length;n++) {
317         if (stdo) {
318           if (verbosity&VER_USR2) {
319             for (i=0;i<comp;i++) 
320               fprintf(stdout,"%e ",series[i][n]*interval[i]+min[i]);
321             fprintf(stdout,"%u\n",nmf[n]);
322           }
323           else {
324             fprintf(stdout,"%e",series[0][n]*interval[0]+min[0]);
325             for (i=1;i<comp;i++)
326               fprintf(stdout," %e",series[i][n]*interval[i]+min[i]);
327             fprintf(stdout,"\n");
328           }
329         }
330         if (!stdo || (verbosity&VER_USR1)) {
331           if (verbosity&VER_USR2) {
332             for (i=0;i<comp;i++) 
333               fprintf(file,"%e ",series[i][n]*interval[i]+min[i]);
334             fprintf(file,"%u\n",nmf[n]);
335           }
336           else {
337             fprintf(file,"%e",series[0][n]*interval[0]+min[0]);
338             for (i=1;i<comp;i++)
339               fprintf(file," %e",series[i][n]*interval[i]+min[i]);
340             fprintf(file,"\n");
341           }
342         }
343       }
344       if (!stdo || (verbosity&VER_USR1))
345         fclose(file);
346     }
347   }
348
349   /*cleaning up */
350   for (i=0;i<comp;i++) {
351     free(series[i]);
352     free(nf[i]);
353     free(corr[i]);
354   }
355   free(series);
356   free(nf);
357   free(corr);
358
359   for (i=0;i<2;i++)
360     free(indexes[i]);
361   free(indexes);
362
363   free(list);
364   free(nmf);
365   free(interval);
366   free(min);
367
368   for (i=0;i<BOX;i++)
369     free(box[i]);
370   free(box);
371
372   if (outfile != NULL)
373     free(outfile);
374   if (ofname != NULL)
375     free(ofname);
376   if (infile != NULL)
377     free(infile);
378   if (column != NULL)
379     free(column);
380   /* end cleaning up */
381
382   return 0;
383 }