IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / lzo-test.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: Aug 27, 2004 */
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include <limits.h>
25 #include "routines/tsa.h"
26
27 #define WID_STR "Estimates the average forecast error for a zeroth\n\t\
28 order fit from a multidimensional time series"
29
30
31 #ifndef _MATH_H
32 #include <math.h>
33 #endif
34
35 /*number of boxes for the neighbor search algorithm*/
36 #define NMAX 512
37
38 unsigned int nmax=(NMAX-1);
39 long **box,*list;
40 unsigned long *found;
41 double **series,**diffs;
42 double interval,min,epsilon;
43
44 char epsset=0,dimset=0,clengthset=0,causalset=0;
45 char *infile=NULL;
46 char *outfile=NULL,stdo=1;
47 char *COLUMNS=NULL;
48 unsigned int embed=2,dim=1,DELAY=1,MINN=30;
49 unsigned long STEP=1,causal;
50 unsigned int verbosity=0x1;
51 double EPS0=1.e-3,EPSF=1.2;
52 unsigned long refstep=1;
53 unsigned long LENGTH=ULONG_MAX,exclude=0,CLENGTH=ULONG_MAX;
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 columns to read [default: 1,...,X]\n");
67   fprintf(stderr,"\t-m dimension and embedding dimension"
68           " [default: %d,%d]\n",dim,embed);
69   fprintf(stderr,"\t-d delay [default: %d]\n",DELAY);
70   fprintf(stderr,"\t-n # of reference points [default: length]\n");
71   fprintf(stderr,"\t-S temporal distance between the reference points"
72           " [default: %lu]\n",refstep);
73   fprintf(stderr,"\t-k minimal number of neighbors for the fit "
74           "[default: %d]\n",MINN);
75   fprintf(stderr,"\t-r neighborhoud size to start with "
76           "[default: (data interval)/1000]\n");
77   fprintf(stderr,"\t-f factor to increase size [default: 1.2]\n");
78   fprintf(stderr,"\t-s steps to forecast [default: 1]\n");
79   fprintf(stderr,"\t-C width of causality window [default: steps]\n");
80   fprintf(stderr,"\t-o output file [default: 'datafile.zer',"
81           " without -o: stdout]\n");
82   fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
83           "0='only panic messages'\n\t\t"
84           "1='+ input/output messages'\n\t\t"
85           "2='give individual forecast errors for the max step'\n");
86   fprintf(stderr,"\t-h show these options\n");
87   exit(0);
88 }
89
90 void scan_options(int n,char **in)
91 {
92   char *out;
93
94   if ((out=check_option(in,n,'l','u')) != NULL)
95     sscanf(out,"%lu",&LENGTH);
96   if ((out=check_option(in,n,'x','u')) != NULL)
97     sscanf(out,"%lu",&exclude);
98   if ((out=check_option(in,n,'c','s')) != NULL)
99     COLUMNS=out;
100   if ((out=check_option(in,n,'m','2')) != NULL) {
101     dimset=1;
102     sscanf(out,"%u%*c%u",&dim,&embed);
103     if (embed == 0)
104       embed=1;
105   }
106   if ((out=check_option(in,n,'d','u')) != NULL)
107     sscanf(out,"%u",&DELAY);
108   if ((out=check_option(in,n,'n','u')) != NULL) {
109     sscanf(out,"%lu",&CLENGTH);
110     clengthset=1;
111   }
112   if ((out=check_option(in,n,'S','u')) != NULL)
113     sscanf(out,"%lu",&refstep);
114   if ((out=check_option(in,n,'k','u')) != NULL)
115     sscanf(out,"%u",&MINN);
116   if ((out=check_option(in,n,'r','f')) != NULL) {
117     epsset=1;
118     sscanf(out,"%lf",&EPS0);
119   }
120   if ((out=check_option(in,n,'f','f')) != NULL)
121     sscanf(out,"%lf",&EPSF);
122   if ((out=check_option(in,n,'s','u')) != NULL)
123     sscanf(out,"%lu",&STEP);
124   if ((out=check_option(in,n,'C','u')) != NULL) {
125     sscanf(out,"%lu",&causal);
126     causalset=1;
127   }
128   if ((out=check_option(in,n,'V','u')) != NULL)
129     sscanf(out,"%u",&verbosity);
130   if ((out=check_option(in,n,'o','o')) != NULL) {
131     stdo=0;
132     if (strlen(out) > 0)
133       outfile=out;
134   }
135 }
136
137 void make_fit(long act,unsigned long number,long istep,double **error)
138 {
139   double casted,*help;
140   long i,j,h;
141   
142   h=istep-1;
143   for (j=0;j<dim;j++) {
144     casted=0.0;
145     help=series[j]+istep;
146     for (i=0;i<number;i++)
147       casted += help[found[i]];
148     casted /= number;
149     diffs[j][act]=casted-help[act];
150     error[j][h] += sqr(casted-help[act]);
151   }
152 }
153
154 int main(int argc,char **argv)
155 {
156   char stdi=0;
157   char alldone,*done;
158   long i,j,hi;
159   unsigned long *hfound;
160   unsigned long actfound;
161   unsigned long clength;
162   double *rms,*av,**error,**hser,*hinter;
163   FILE *file;
164
165   if (scan_help(argc,argv))
166     show_options(argv[0]);
167   
168   scan_options(argc,argv);
169
170   if ((2*STEP+causal) >= ((long)LENGTH-(long)(embed*DELAY)-(long)MINN)) {
171     fprintf(stderr,"steps to forecast (-s) too large. Exiting!\n");
172     exit(ZEROTH__STEP_TOO_LARGE);
173   }
174   if (!causalset)
175     causal=STEP;
176
177 #ifndef OMIT_WHAT_I_DO
178   if (verbosity&VER_INPUT)
179     what_i_do(argv[0],WID_STR);
180 #endif
181
182   infile=search_datafile(argc,argv,NULL,verbosity);
183   if (infile == NULL)
184     stdi=1;
185
186   if (outfile == NULL) {
187     if (!stdi) {
188       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
189       sprintf(outfile,"%s.zer",infile);
190     }
191     else {
192       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
193       sprintf(outfile,"stdin.zer");
194     }
195   }
196   if (!stdo)
197     test_outfile(outfile);
198   
199   if (COLUMNS == NULL)
200     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,"",dimset,
201                                       verbosity);
202   else
203     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,COLUMNS,
204                                       dimset,verbosity);
205
206   check_alloc(hser=(double**)malloc(sizeof(double*)*dim));
207   check_alloc(av=(double*)malloc(sizeof(double)*dim));
208   check_alloc(rms=(double*)malloc(sizeof(double)*dim));
209   check_alloc(hinter=(double*)malloc(sizeof(double)*dim));
210   interval=0.0;
211   for (i=0;i<dim;i++) {
212     rescale_data(series[i],LENGTH,&min,&hinter[i]);
213     variance(series[i],LENGTH,&av[i],&rms[i]);
214     interval += hinter[i];
215   }
216   interval /= (double)dim;
217
218   check_alloc(list=(long*)malloc(sizeof(long)*LENGTH));
219   check_alloc(found=(unsigned long*)malloc(sizeof(long)*LENGTH));
220   check_alloc(hfound=(unsigned long*)malloc(sizeof(long)*LENGTH));
221   check_alloc(done=(char*)malloc(sizeof(char)*LENGTH));
222   check_alloc(box=(long**)malloc(sizeof(long*)*NMAX));
223   check_alloc(error=(double**)malloc(sizeof(double*)*dim));
224   check_alloc(diffs=(double**)malloc(sizeof(double*)*dim));
225   for (j=0;j<dim;j++) {
226     check_alloc(diffs[j]=(double*)malloc(sizeof(double)*LENGTH));
227     check_alloc(error[j]=(double*)malloc(sizeof(double)*STEP));
228     for (i=0;i<STEP;i++)
229       error[j][i]=0.0;
230   }
231   
232   for (i=0;i<NMAX;i++)
233     check_alloc(box[i]=(long*)malloc(sizeof(long)*NMAX));
234     
235   for (i=0;i<LENGTH;i++)
236     done[i]=0;
237
238   alldone=0;
239   if (epsset)
240     EPS0 /= interval;
241
242   epsilon=EPS0/EPSF;
243
244   if (!clengthset)
245     CLENGTH=LENGTH;
246   clength=((CLENGTH*refstep+STEP) <= LENGTH) ? CLENGTH : 
247     (LENGTH-(long)STEP)/refstep;
248
249   while (!alldone) {
250     alldone=1;
251     epsilon*=EPSF;
252     make_multi_box(series,box,list,LENGTH-(long)STEP,NMAX,(unsigned int)dim,
253                    (unsigned int)embed,(unsigned int)DELAY,epsilon);
254     for (i=(embed-1)*DELAY;i<clength;i++)
255       if (!done[i]) {
256         hi=i*refstep;
257         for (j=0;j<dim;j++)
258           hser[j]=series[j]+hi;
259         actfound=find_multi_neighbors(series,box,list,hser,LENGTH,NMAX,
260                                        (unsigned int)dim,(unsigned int)embed,
261                                        (unsigned int)DELAY,epsilon,hfound);
262         actfound=exclude_interval(actfound,hi-(long)causal+1,
263                                   hi+causal+(embed-1)*DELAY-1,hfound,found);    
264         if (actfound >= MINN) {
265           for (j=1;j<=STEP;j++) {
266             make_fit(hi,actfound,j,error);
267           }
268           done[i]=1;
269         }
270         alldone &= done[i];
271       }
272   }
273   if (stdo) {
274     if (verbosity&VER_INPUT)
275       fprintf(stderr,"Writing to stdout\n");
276     for (i=0;i<STEP;i++) {
277       if (verbosity&VER_USR1)
278         fprintf(stdout,"# %lu ",i+1);
279       else
280         fprintf(stdout,"%lu ",i+1);
281       for (j=0;j<dim;j++) 
282         fprintf(stdout,"%e ",
283                 sqrt(error[j][i]/(clength-(embed-1)*DELAY))/rms[j]);
284       fprintf(stdout,"\n");
285     }
286     if (verbosity&VER_USR1) {
287       for (i=(embed-1)*DELAY;i<clength;i++) {
288         hi=i*refstep;
289         for (j=0;j<dim;j++)
290           fprintf(stdout,"%e ",diffs[j][hi]*hinter[j]);
291         fprintf(stdout,"\n");
292       }
293     }
294   }
295   else {
296     file=fopen(outfile,"w");
297     if (verbosity&VER_INPUT)
298       fprintf(stderr,"Opened %s for writing\n",outfile);
299     for (i=0;i<STEP;i++) {
300       if (verbosity&VER_USR1)
301         fprintf(file,"# %lu ",i+1);
302       else
303         fprintf(file,"%lu ",i+1);
304       for (j=0;j<dim;j++) 
305         fprintf(file,"%e ",sqrt(error[j][i]/(clength-(embed-1)*DELAY))/rms[j]);
306       fprintf(file,"\n");
307     }
308     if (verbosity&VER_USR1) {
309       for (i=(embed-1)*DELAY;i<clength;i++) {
310         hi=i*refstep;
311         for (j=0;j<dim;j++)
312           fprintf(file,"%e ",diffs[j][hi]*hinter[j]);
313         fprintf(file,"\n");
314       }
315     }
316     fclose(file);
317   }
318
319   if (outfile != NULL)
320     free(outfile);
321   if (infile != NULL)
322     free(infile);
323   if (COLUMNS != NULL)
324     free(COLUMNS);
325   for (i=0;i<dim;i++) {
326     free(series[i]);
327     free(diffs[i]);
328     free(error[i]);
329   }
330   free(series);
331   free(diffs);
332   free(hser);
333   free(error);
334   free(av);
335   free(rms);
336   free(list);
337   free(found);
338   free(hfound);
339   free(done);
340   for (i=0;i<NMAX;i++)
341     free(box[i]);
342   free(box);
343
344   return 0;
345 }