IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / ghkss.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: Jun 10, 2006 */
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include <math.h>
25 #include <limits.h>
26 #include "routines/tsa.h"
27
28 #define WID_STR "Multivariate noise reduction using the GHKSS algorithm"
29
30
31 #define BOX (unsigned int)1024
32
33 unsigned long length=ULONG_MAX,exclude=0;
34 unsigned int dim,qdim=2,delay=1,minn=50,iterations=1,comp=1,embed=5;
35 unsigned int verbosity=0xff;
36 double mineps,epsfac;
37 char *column=NULL;
38 char eps_set=0,euclidean=0,dimset=0,resize_eps;
39 char *outfile=NULL,stdo=1;
40 char *infile=NULL;
41
42 double *d_min,*d_max,d_max_max;
43 double **series,**delta,**corr;
44 double *metric,trace;
45 long **box,*list;
46 unsigned long *flist;
47 int emb_offset;
48 unsigned int ibox=BOX-1;
49 unsigned int *index_comp,*index_embed;
50
51 /*these are global to save time*/
52 int *sorted;
53 double *av,**mat,*matarray,*eig;
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 to use [Default: whole file]\n");
65   fprintf(stderr,"\t-x # of lines to be ignored [Default: 0]\n");
66   fprintf(stderr,"\t-c column to read [Default: 1,..,# of components]\n");
67   fprintf(stderr,"\t-m # of components,embedding dimension [Default: 1,5]\n");
68   fprintf(stderr,"\t-d delay [Default: 1]\n");
69   fprintf(stderr,"\t-q dimension to project to [Default: 2]\n");
70   fprintf(stderr,"\t-k minimal number of neighbours [Default: 50]\n");
71   fprintf(stderr,"\t-r minimal neighbourhood size \n\t\t"
72           "[Default: (interval of data)/1000]\n");
73   fprintf(stderr,"\t-i # of iterations [Default: 1]\n");
74   fprintf(stderr,"\t-2 use euklidean metric [Default: non euklidean]\n");
75   fprintf(stderr,"\t-o name of output file \n\t\t"
76           "[Default: 'datafile'.opt.n, where n is the iteration.\n\t\t"
77           " If no -o is given, the last iteration is also"
78           " written to stdout]\n");
79   fprintf(stderr,"\t-V verbosity level [Default: 7]\n\t\t"
80           "0='only panic messages'\n\t\t"
81           "1='+ input/output messages'\n\t\t"
82           "2='+ average correction and trend'\n\t\t"
83           "4='+ how many points for which epsilon'\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",&comp,&embed);
102   if ((out=check_option(in,n,'d','u')) != NULL)
103     sscanf(out,"%u",&delay);
104   if ((out=check_option(in,n,'q','u')) != NULL)
105     sscanf(out,"%u",&qdim);
106   if ((out=check_option(in,n,'k','u')) != NULL)
107     sscanf(out,"%u",&minn);
108   if ((out=check_option(in,n,'r','f')) != NULL) {
109     eps_set=1;
110     sscanf(out,"%lf",&mineps);
111   }
112   if ((out=check_option(in,n,'i','u')) != NULL)
113     sscanf(out,"%u",&iterations);
114   if ((out=check_option(in,n,'V','u')) != NULL)
115     sscanf(out,"%u",&verbosity);
116   if ((out=check_option(in,n,'2','n')) != NULL)
117     euclidean=1;
118   if ((out=check_option(in,n,'o','o')) != NULL) {
119     stdo=0;
120     if (strlen(out) > 0)
121       outfile=out;
122   }
123 }
124
125 void sort(double *x,int *n)
126 {
127   long i,j,iswap;
128   double dswap;
129   
130   for (i=0;i<dim;i++)
131     n[i]=i;
132   
133   for (i=0;i<dim-1;i++)
134     for (j=i+1;j<dim;j++)
135       if (x[j] > x[i]) {
136         dswap=x[i];
137         x[i]=x[j];
138         x[j]=dswap;
139         iswap=n[i];
140         n[i]=n[j];
141         n[j]=iswap;
142       }
143 }
144
145 void mmb(double eps)
146 {
147   long i,x,y;
148   double ieps=1.0/eps;
149
150   for (x=0;x<BOX;x++)
151     for (y=0;y<BOX;y++)
152       box[x][y] = -1;
153   
154   for (i=emb_offset;i<length;i++) {
155     x=(int)(series[0][i]*ieps)&ibox;
156     y=(int)(series[comp-1][i-emb_offset]*ieps)&ibox;
157     list[i]=box[x][y];
158     box[x][y]=i;
159   }
160 }
161
162 unsigned long fmn(long which,double eps)
163 {
164   unsigned long nf=0;
165   long i,i1,i2,j,j1,k,k1,li;
166   long element;
167   double dx=0.0;
168   
169   i=(int)(series[0][which]/eps)&ibox;
170   j=(int)(series[comp-1][which-emb_offset]/eps)&ibox;
171   
172   for (i1=i-1;i1<=i+1;i1++) {
173     i2=i1&ibox;
174     for (j1=j-1;j1<=j+1;j1++) {
175       element=box[i2][j1&ibox];
176       while (element != -1) {
177         for (k=0;k<embed;k++) {
178           k1= -k*(int)delay;
179           for (li=0;li<comp;li++) {
180             dx=fabs(series[li][which+k1]-series[li][element+k1]);
181             if (dx > eps)
182               break;
183           }
184           if (dx > eps)
185             break;
186         }
187         if (dx <= eps)
188           flist[nf++]=element;
189         element=list[element];
190       }
191     }
192   }
193   return nf;
194 }
195
196 void make_correction(unsigned long n,unsigned long nf)
197 {
198   long i,i1,i2,j,j1,j2,k,k1,k2,hs;
199   double help;
200   
201   for (i=0;i<dim;i++) {
202     i1=index_comp[i];
203     i2=index_embed[i];
204     help=0.0;
205     for (j=0;j<nf;j++)
206       help += series[i1][flist[j]-i2];
207     av[i]=help/nf;
208   }
209
210   for (i=0;i<dim;i++) {
211     i1=index_comp[i];
212     i2=index_embed[i];
213     for (j=i;j<dim;j++) {
214       help=0.0;
215       j1=index_comp[j];
216       j2=index_embed[j];
217       for (k=0;k<nf;k++) {
218         hs=flist[k];
219         help += series[i1][hs-i2]*series[j1][hs-j2];
220       }
221       mat[i][j]=(help/nf-av[i]*av[j])*metric[i]*metric[j];
222       mat[j][i]=mat[i][j];
223     }
224   }
225
226   eigen(mat,(unsigned long)dim,eig);
227   sort(eig,sorted);
228
229   for (i=0;i<dim;i++) {
230     help=0.0;
231     for (j=qdim;j<dim;j++) {
232       hs=sorted[j];
233       for (k=0;k<dim;k++) {
234         k1=index_comp[k];
235         k2=index_embed[k];
236         help += (series[k1][n-k2]-av[k])*mat[k][hs]*mat[i][hs]*metric[k];
237       }
238     }
239     corr[n][i]=help/metric[i];
240   }
241 }
242
243 void handle_trend(unsigned long n,unsigned long nf)
244 {
245   long i,i1,i2,j;
246   double help;
247   
248   for (i=0;i<dim;i++) {
249     help=0.0;
250     for (j=0;j<nf;j++)
251       help += corr[flist[j]][i];
252     av[i]=help/nf;
253   }
254
255   for (i=0;i<dim;i++) {
256     i1=index_comp[i];
257     i2=index_embed[i];
258     delta[i1][n-i2] += (corr[n][i]-av[i])/(trace*metric[i]);
259   }
260 }
261
262 void set_correction(void)
263 {
264   long i,j;
265   double *hav,*hsigma,help;
266
267   check_alloc(hav=(double*)malloc(sizeof(double)*comp));
268   check_alloc(hsigma=(double*)malloc(sizeof(double)*comp));
269   for (j=0;j<comp;j++)
270     hav[j]=hsigma[j]=0.0;
271
272   for (i=0;i<length;i++)
273     for (j=0;j<comp;j++) {
274       hav[j] += (help=delta[j][i]);
275       hsigma[j] += help*help;
276     }
277
278   for (j=0;j<comp;j++) {
279     hav[j] /= length;
280     hsigma[j]=sqrt(fabs(hsigma[j]/length-hav[j]*hav[j]));
281   }
282   if (verbosity&(VER_USR1|VER_USR2)) {
283     for (i=0;i<comp;i++) {
284       fprintf(stderr,"Average shift of component %ld = %e\n",i+1,
285               hav[i]*d_max[i]);
286       fprintf(stderr,"Average rms correction of comp. %ld = %e\n\n",
287               i+1,hsigma[i]*d_max[i]);
288     }
289   }
290   for (i=0;i<length;i++)
291     for (j=0;j<comp;j++)
292       series[j][i] -= delta[j][i];
293
294   if (resize_eps) {
295     mineps /= epsfac;
296     if (verbosity&VER_USR2)
297       fprintf(stderr,"Reset minimal neighbourhood size to %e\n",
298               mineps*d_max_max);
299   }
300
301   resize_eps=0;
302   free(hav);
303   free(hsigma);
304 }
305
306 int main(int argc,char **argv)
307 {
308   char stdi=0;
309   int iter,epscount,*ok;
310   long i,j;
311   char all_done;
312   char *ofname;
313   unsigned long nfound,n,allfound;
314   double epsilon;
315   double **hser;
316   FILE *file;
317
318   if (scan_help(argc,argv))
319     show_options(argv[0]);
320   
321   scan_options(argc,argv);
322 #ifndef OMIT_WHAT_I_DO
323   if (verbosity&VER_INPUT)
324     what_i_do(argv[0],WID_STR);
325 #endif
326
327   dim=comp*embed;
328   emb_offset=(embed-1)*delay;
329
330   infile=search_datafile(argc,argv,NULL,verbosity);
331   if (infile == NULL)
332     stdi=1;
333
334   if (outfile == NULL) {
335     if (!stdi) {
336       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
337       check_alloc(ofname=(char*)calloc(strlen(infile)+9,(size_t)1));
338       sprintf(outfile,"%s.opt",infile);
339     }
340     else {
341       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
342       check_alloc(ofname=(char*)calloc((size_t)14,(size_t)1));
343       sprintf(outfile,"stdin.opt");
344     }
345   }
346   else
347     check_alloc(ofname=(char*)calloc(strlen(outfile)+10,(size_t)1));
348   
349   if (column == NULL)
350     series=(double**)get_multi_series(infile,&length,exclude,&comp,"",
351                                      dimset,verbosity);
352   else 
353     series=(double**)get_multi_series(infile,&length,exclude,&comp,column,
354                                       dimset,verbosity);
355
356   if (length < minn) {
357     fprintf(stderr,"With %lu data you will never find %u neighbors."
358             " Exiting!\n",length,minn);
359     exit(GHKSS__TOO_MANY_NEIGHBORS);
360   }
361
362   check_alloc(d_min=(double*)malloc(sizeof(double)*comp));
363   check_alloc(d_max=(double*)malloc(sizeof(double)*comp));
364   d_max_max=0.0;
365   for (i=0;i<comp;i++) {
366     rescale_data(series[i],length,&d_min[i],&d_max[i]);
367     if (d_max[i] > d_max_max)
368       d_max_max=d_max[i];
369   }
370
371   if (!eps_set)
372     mineps=1./1000.;
373   else
374     mineps /= d_max_max;
375   epsfac=sqrt(2.0);
376
377   check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
378   for (i=0;i<BOX;i++)
379     check_alloc(box[i]=(long*)malloc(sizeof(long)*BOX));
380
381   check_alloc(list=(long*)malloc(sizeof(long)*length));
382   check_alloc(flist=(unsigned long*)malloc(sizeof(long)*length));
383
384   check_alloc(metric=(double*)malloc(sizeof(double)*dim));
385   trace=0.0;
386   if (euclidean) {
387     for (i=0;i<dim;i++) {
388       metric[i]=1.0;
389       trace += 1./metric[i];
390     }
391   }
392   else {
393     for (i=0;i<dim;i++) {
394       if ((i >= comp) && (i < ((long)dim-(long)comp))) 
395         metric[i]=1.0;
396       else 
397         metric[i]=1.0e3;
398       trace += 1./metric[i];
399     }
400   }
401
402   check_alloc(corr=(double**)malloc(sizeof(double*)*length));
403   for (i=0;i<length;i++)
404     check_alloc(corr[i]=(double*)malloc(sizeof(double)*dim));
405   check_alloc(ok=(int*)malloc(sizeof(int)*length));
406   check_alloc(delta=(double**)malloc(sizeof(double*)*comp));
407   for (i=0;i<comp;i++)
408     check_alloc(delta[i]=(double*)malloc(sizeof(double)*length));
409   check_alloc(index_comp=(unsigned int*)malloc(sizeof(int)*dim));
410   check_alloc(index_embed=(unsigned int*)malloc(sizeof(int)*dim));
411   check_alloc(av=(double*)malloc(sizeof(double)*dim));
412   check_alloc(sorted=(int*)malloc(sizeof(int)*dim));
413   check_alloc(eig=(double*)malloc(sizeof(double)*dim));
414   check_alloc(matarray=(double*)malloc(sizeof(double)*dim*dim));
415   check_alloc(mat=(double**)malloc(sizeof(double*)*dim));
416   for (i=0;i<dim;i++)
417     mat[i]=(double*)(matarray+dim*i);
418   check_alloc(hser=(double**)malloc(sizeof(double*)*comp));
419
420   for (i=0;i<dim;i++) {
421     index_comp[i]=i%comp;
422     index_embed[i]=(i/comp)*delay;
423   }
424
425   resize_eps=0;
426   for (iter=1;iter<=iterations;iter++) {
427     for (i=0;i<length;i++) {
428       ok[i]=0;
429       for (j=0;j<dim;j++)
430         corr[i][j]=0.0;
431       for (j=0;j<comp;j++)
432         delta[j][i]=0.0;
433     }
434     epsilon=mineps;
435     all_done=0;
436     epscount=1;
437     allfound=0;
438     if (verbosity&(VER_USR1|VER_USR2))
439       fprintf(stderr,"Starting iteration %d\n",iter);
440     while(!all_done) {
441       mmb(epsilon);
442       all_done=1;
443       for (n=emb_offset;n<length;n++)
444         if (!ok[n]) {
445           nfound=fmn(n,epsilon);
446           if (nfound >= minn) {
447             make_correction(n,nfound);
448             ok[n]=epscount;
449             if (epscount == 1)
450               resize_eps=1;
451             allfound++;
452           }
453           else
454             all_done=0;
455         }
456       if (verbosity&VER_USR2)
457         fprintf(stderr,"Corrected %ld points with epsilon= %e\n",allfound,
458                 epsilon*d_max_max);
459       epsilon *= epsfac;
460       epscount++;
461     }
462     if (verbosity&VER_USR2)
463       fprintf(stderr,"Start evaluating the trend\n");
464
465     epsilon=mineps;
466     allfound=0;
467     for (i=1;i<epscount;i++) {
468       mmb(epsilon);
469       for (n=emb_offset;n<length;n++)
470         if (ok[n] == i) {
471           nfound=fmn(n,epsilon);
472           handle_trend(n,nfound);
473           allfound++;
474         }
475       if (verbosity&VER_USR2)
476         fprintf(stderr,"Trend subtracted for %ld points with epsilon= %e\n",
477                 allfound,epsilon*d_max_max);
478       epsilon *= epsfac;
479     }
480     set_correction();
481     
482     sprintf(ofname,"%s.%d",outfile,iter);
483     test_outfile(ofname);
484
485     file=fopen(ofname,"w");
486     if (verbosity&VER_INPUT)
487       fprintf(stderr,"Opened %s for writing\n\n",ofname);
488     for (i=0;i<length;i++) {
489       for (j=0;j<comp;j++) {
490         fprintf(file,"%e ",series[j][i]*d_max[j]+d_min[j]);
491       }
492       fprintf(file,"\n");
493       if (stdo && (iter == iterations)) {
494         for (j=0;j<comp;j++)
495           fprintf(stdout,"%e ",series[j][i]*d_max[j]+d_min[j]);
496         fprintf(stdout,"\n");
497       }
498     }
499     fclose(file);
500   }
501
502   return 0;
503 }