IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / lfo-ar.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 21, 2005 */
21 /*changes:
22   Jun 17, 2005: Comments in the output file updated
23   Jun 21, 2005: free imat in make_fit
24 */
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <string.h>
28 #include <limits.h>
29 #include "routines/tsa.h"
30 #include <math.h>
31
32 #define WID_STR "Estimates the average forecast error for a local\n\t\
33 linear fit as a function of the neighborhood size."
34
35
36 /*number of boxes for the neighbor search algorithm*/
37 #define NMAX 256
38
39 unsigned int nmax=(NMAX-1);
40 long **box,*list;
41 unsigned long *found;
42 double *error;
43 double **series;
44
45 char eps0set=0,eps1set=0,causalset=0,dimset=0;
46 char *outfile=NULL,stdo=1;
47 char *column=NULL;
48 unsigned int dim=1,embed=2,delay=1;
49 unsigned int verbosity=0xff;
50 int STEP=1;
51 double EPS0=1.e-3,EPS1=1.0,EPSF=1.2;
52 unsigned long LENGTH=ULONG_MAX,exclude=0,CLENGTH=ULONG_MAX,causal;
53 char *infile=NULL;
54 double **mat,*vec,*localav,*foreav,*hvec;
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 use [default: whole file]\n");
66   fprintf(stderr,"\t-x # of lines to be ignored [default: 0]\n");
67   fprintf(stderr,"\t-c columns to read [default: 1,...,# of components]\n");
68   fprintf(stderr,"\t-m # of components,embedding dimension [default: 1,2]\n");
69   fprintf(stderr,"\t-d delay [default: 1]\n");
70   fprintf(stderr,"\t-i iterations [default: length]\n");
71   fprintf(stderr,"\t-r neighborhood size to start with [default:"
72           " (interval of data)/1000)]\n");
73   fprintf(stderr,"\t-R neighborhood size to end with [default:"
74           " interval of data]\n");
75   fprintf(stderr,"\t-f factor to increase size [default: 1.2]\n");
76   fprintf(stderr,"\t-s steps to forecast [default: 1]\n");
77   fprintf(stderr,"\t-C width of causality window [default: steps]\n");
78   fprintf(stderr,"\t-o output file name [default: 'datafile.ll']\n");
79   fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
80           "0='only panic messages'\n\t\t"
81           "1='+ input/output messages'\n");
82   fprintf(stderr,"\t-h show these options\n");
83   exit(0);
84 }
85
86 void scan_options(int n,char **in)
87 {
88   char *out;
89
90   if ((out=check_option(in,n,'l','u')) != NULL)
91     sscanf(out,"%lu",&LENGTH);
92   if ((out=check_option(in,n,'x','u')) != NULL)
93     sscanf(out,"%lu",&exclude);
94   if ((out=check_option(in,n,'c','s')) != NULL) {
95     column=out;
96     dimset=1;
97   }
98   if ((out=check_option(in,n,'m','2')) != NULL)
99     sscanf(out,"%u,%u",&dim,&embed);
100   if ((out=check_option(in,n,'d','u')) != NULL)
101     sscanf(out,"%u",&delay);
102   if ((out=check_option(in,n,'i','u')) != NULL)
103     sscanf(out,"%lu",&CLENGTH);
104   if ((out=check_option(in,n,'r','f')) != NULL) {
105     eps0set=1;
106     sscanf(out,"%lf",&EPS0);
107   }
108   if ((out=check_option(in,n,'R','f')) != NULL) {
109     eps1set=1;
110     sscanf(out,"%lf",&EPS1);
111   }
112   if ((out=check_option(in,n,'f','f')) != NULL)
113     sscanf(out,"%lf",&EPSF);
114   if ((out=check_option(in,n,'s','u')) != NULL)
115     sscanf(out,"%u",&STEP);
116   if ((out=check_option(in,n,'C','u')) != NULL) {
117     sscanf(out,"%lu",&causal);
118     causalset=1;
119   }
120   if ((out=check_option(in,n,'V','u')) != NULL)
121     sscanf(out,"%u",&verbosity);
122   if ((out=check_option(in,n,'o','o')) != NULL) {
123     stdo=0;
124     if (strlen(out) > 0)
125       outfile=out;
126   }
127 }
128
129 void multiply_matrix(double **mat,double *vec)
130 {
131   long i,j;
132
133   for (i=0;i<dim*embed;i++) {
134     hvec[i]=0.0;
135     for (j=0;j<dim*embed;j++)
136       hvec[i] += mat[i][j]*vec[j];
137   }
138   for (i=0;i<dim*embed;i++)
139     vec[i]=hvec[i];
140 }
141
142 void make_fit(long act,unsigned long number)
143 {
144   double *si,*sj,lavi,lavj,fav,**imat,cast;
145   long i,i1,hi,hi1,j,j1,hj,hj1,n,which;
146   
147   for (i=0;i<embed*dim;i++)
148     localav[i]=0;
149   for (i=0;i<dim;i++)
150     foreav[i]=0.0;
151   
152   for (n=0;n<number;n++) {
153     which=found[n];
154     for (j=0;j<dim;j++) {
155       sj=series[j];
156       foreav[j] += sj[which+STEP];
157       for (j1=0;j1<embed;j1++) {
158         hj=j*embed+j1;
159         localav[hj] += sj[which-j1*delay];
160       }
161     }
162   }
163
164   for (i=0;i<dim*embed;i++)
165     localav[i] /= number;
166   for (i=0;i<dim;i++)
167     foreav[i] /= number;
168
169   for (i=0;i<dim;i++) {
170     si=series[i];
171     for (i1=0;i1<embed;i1++) {
172       hi=i*embed+i1;
173       lavi=localav[hi];
174       hi1=i1*delay;
175       for (j=0;j<dim;j++) {
176         sj=series[j];
177         for (j1=0;j1<embed;j1++) {
178           hj=j*embed+j1;
179           lavj=localav[hj];
180           hj1=j1*delay;
181           mat[hi][hj]=0.0;
182           if (hj >= hi) {
183             for (n=0;n<number;n++) {
184               which=found[n];
185               mat[hi][hj] += (si[which-hi1]-lavi)*(sj[which-hj1]-lavj);
186             }
187           }
188         }
189       }
190     }
191   }
192   
193   for (i=0;i<dim*embed;i++)
194     for (j=i;j<dim*embed;j++) {
195       mat[i][j] /= number;
196       mat[j][i]=mat[i][j];
197     }
198   
199   imat=invert_matrix(mat,dim*embed);
200
201   for (i=0;i<dim;i++) {
202     si=series[i];
203     fav=foreav[i];
204     for (j=0;j<dim;j++) {
205       sj=series[j];
206       for (j1=0;j1<embed;j1++) {
207         hj=j*embed+j1;
208         lavj=localav[hj];
209         hj1=j1*delay;
210         vec[hj]=0.0;
211         for (n=0;n<number;n++) {
212           which=found[n];
213           vec[hj] += (si[which+STEP]-fav)*(sj[which-hj1]-lavj);
214         }
215         vec[hj] /= number;
216       }
217     }
218
219     multiply_matrix(imat,vec);
220
221     cast=foreav[i];
222     for (j=0;j<dim;j++) {
223       sj=series[j];
224       for (j1=0;j1<embed;j1++) {
225         hj=j*embed+j1;
226         cast += vec[hj]*(sj[act-j1*delay]-localav[hj]);
227       }
228     }
229     error[i] += sqr(cast-series[i][act+STEP]);
230   }
231   for (i=0;i<embed*dim;i++)
232     free(imat[i]);
233   free(imat);
234 }
235
236 int main(int argc,char **argv)
237 {
238   char stdi=0;
239   unsigned long actfound;
240   unsigned long *hfound;
241   long pfound,i,j;
242   unsigned long clength;
243   double interval,min,maxinterval;
244   double epsilon;
245   double **hser;
246   double avfound,*hrms,*hav,sumerror=0.0;
247   FILE *file=NULL;
248
249   if (scan_help(argc,argv))
250     show_options(argv[0]);
251   
252   scan_options(argc,argv);
253 #ifndef OMIT_WHAT_I_DO
254   if (verbosity&VER_INPUT)
255     what_i_do(argv[0],WID_STR);
256 #endif
257
258   if (!causalset)
259     causal=STEP;
260
261   infile=search_datafile(argc,argv,NULL,verbosity);
262   if (infile == NULL)
263     stdi=1;
264
265   if (outfile == NULL) {
266     if (!stdi) {
267       check_alloc(outfile=(char*)calloc(strlen(infile)+4,(size_t)1));
268       sprintf(outfile,"%s.ll",infile);
269     }
270     else {
271       check_alloc(outfile=(char*)calloc((size_t)9,(size_t)1));
272       sprintf(outfile,"stdin.ll");
273     }
274   }
275   if (!stdo)
276     test_outfile(outfile);
277
278   if (column == NULL)
279     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,"",dimset,
280                                       verbosity);
281   else
282     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,column,
283                                       dimset,verbosity);
284   maxinterval=0.0;
285   for (i=0;i<dim;i++) {
286     rescale_data(series[i],LENGTH,&min,&interval);
287     if (interval > maxinterval)
288       maxinterval=interval;
289   }
290   interval=maxinterval;
291
292   check_alloc(list=(long*)malloc(sizeof(long)*LENGTH));
293   check_alloc(found=(unsigned long*)malloc(sizeof(long)*LENGTH));
294   check_alloc(hfound=(unsigned long*)malloc(sizeof(long)*LENGTH));
295   check_alloc(box=(long**)malloc(sizeof(long*)*NMAX));
296   for (i=0;i<NMAX;i++)
297     check_alloc(box[i]=(long*)malloc(sizeof(long)*NMAX));
298   check_alloc(vec=(double*)malloc(sizeof(double)*(embed*dim)));
299   check_alloc(hvec=(double*)malloc(sizeof(double)*(embed*dim)));
300   check_alloc(mat=(double**)malloc(sizeof(double*)*(embed*dim)));
301   for (i=0;i<dim*embed;i++)
302     check_alloc(mat[i]=(double*)malloc(sizeof(double)*(embed*dim)));
303   check_alloc(error=(double*)malloc(sizeof(double)*dim));
304   check_alloc(hrms=(double*)malloc(sizeof(double)*dim));
305   check_alloc(hav=(double*)malloc(sizeof(double)*dim));
306   check_alloc(hser=(double**)malloc(sizeof(double*)*dim));
307   check_alloc(foreav=(double*)malloc(sizeof(double)*dim));
308   check_alloc(localav=(double*)malloc(sizeof(double)*(embed*dim)));
309   
310   if (eps0set)
311     EPS0 /= interval;
312   if (eps1set)
313     EPS1 /= interval;
314
315   clength=(CLENGTH <= LENGTH) ? CLENGTH-STEP : LENGTH-STEP;
316
317   if (!stdo) {
318     file=fopen(outfile,"w");
319     if (verbosity&VER_INPUT)
320       fprintf(stderr,"Opened %s for writing\n",outfile);
321     fprintf(file,"#1.) neighborhood size\n");
322     fprintf(file,"#2.) average relative forecast error\n");
323     fprintf(file,"#next n.) relative forecast error of the n components\n");
324     fprintf(file,"#second last.) fraction of points with enough neighbors\n");
325     fprintf(file,"#last .) average number of neighbors used for the fit\n");
326   }
327   else {
328     if (verbosity&VER_INPUT)
329       fprintf(stderr,"Writing to stdout\n");
330   }
331
332   for (epsilon=EPS0;epsilon<EPS1*EPSF;epsilon*=EPSF) {
333     pfound=0;
334     for (i=0;i<dim;i++)
335       error[i]=hrms[i]=hav[i]=0.0;
336     avfound=0.0;
337     make_multi_box(series,box,list,LENGTH-STEP,NMAX,dim,
338                    embed,delay,epsilon);
339     for (i=(embed-1)*delay;i<clength;i++) {
340       for (j=0;j<dim;j++)
341         hser[j]=series[j]+i;
342       actfound=find_multi_neighbors(series,box,list,hser,LENGTH,
343                                     NMAX,dim,embed,delay,epsilon,hfound);
344       actfound=exclude_interval(actfound,i-causal+1,i+causal+(embed-1)*delay-1,
345                                 hfound,found);
346       if (actfound > 2*(dim*embed+1)) {
347         make_fit(i,actfound);
348         pfound++;
349         avfound += (double)(actfound-1);
350         for (j=0;j<dim;j++) {
351           hrms[j] += series[j][i+STEP]*series[j][i+STEP];
352           hav[j] += series[j][i+STEP];
353         }
354       }
355     }
356     if (pfound > 1) {
357       sumerror=0.0;
358       for (j=0;j<dim;j++) {
359         hav[j] /= pfound;
360         hrms[j]=sqrt(fabs(hrms[j]/(pfound-1)-hav[j]*hav[j]*pfound/(pfound-1)));
361         error[j]=sqrt(error[j]/pfound)/hrms[j];
362         sumerror += error[j];
363       }
364     }
365     if (stdo) {
366       if (pfound > 1) {
367         fprintf(stdout,"%e %e ",epsilon*interval,sumerror/(double)dim);
368         for (j=0;j<dim;j++)
369           fprintf(stdout,"%e ",error[j]);
370         fprintf(stdout,"%e %e\n",(double)pfound/(clength-(embed-1)*delay),
371                 avfound/pfound);
372         fflush(stdout);
373       }
374     }
375     else {
376       if (pfound > 1) {
377         fprintf(file,"%e %e ",epsilon*interval,sumerror/(double)dim);
378         for (j=0;j<dim;j++)
379           fprintf(file,"%e ",error[j]);
380         fprintf(file,"%e %e\n",(double)pfound/(clength-(embed-1)*delay),
381                 avfound/pfound);
382         fflush(file);
383       }
384     }
385   }
386   if (!stdo)
387     fclose(file);
388   
389   return 0;
390 }