IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / lzo-gm.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: Sep 7, 2004 */
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include <limits.h>
25 #include "routines/tsa.h"
26 #include <math.h>
27
28 #define WID_STR "Estimates the average forecast error for a local\n\t\
29 constant fit as a function of the neighborhood size."
30
31
32 /*number of boxes for the neighbor search algorithm*/
33 #define NMAX 256
34
35 unsigned int nmax=(NMAX-1);
36 long **box,*list;
37 unsigned long *found;
38 double *error;
39 double **series;
40
41 char eps0set=0,eps1set=0,causalset=0,dimset=0;
42 char *outfile=NULL,stdo=1;
43 char *column=NULL;
44 unsigned int dim=1,embed=2,delay=1;
45 unsigned int verbosity=0xff;
46 int STEP=1;
47 double EPS0=1.e-3,EPS1=1.0,EPSF=1.2;
48 unsigned long LENGTH=ULONG_MAX,exclude=0,CLENGTH=ULONG_MAX,causal;
49 char *infile=NULL;
50
51 void show_options(char *progname)
52 {
53   what_i_do(progname,WID_STR);
54   fprintf(stderr," Usage: %s [options]\n",progname);
55   fprintf(stderr," Options:\n");
56   fprintf(stderr,"Everything not being a valid option will be interpreted"
57           " as a possible"
58           " datafile.\nIf no datafile is given stdin is read. Just - also"
59           " means stdin\n");
60   fprintf(stderr,"\t-l # of data to use [default: whole file]\n");
61   fprintf(stderr,"\t-x # of lines to be ignored [default: 0]\n");
62   fprintf(stderr,"\t-c columns to read [default: 1,...,# of components]\n");
63   fprintf(stderr,"\t-m # of components,embedding dimension [default: 1,2]\n");
64   fprintf(stderr,"\t-d delay [default: 1]\n");
65   fprintf(stderr,"\t-i iterations [default: length]\n");
66   fprintf(stderr,"\t-r neighborhood size to start with [default:"
67           " (interval of data)/1000)]\n");
68   fprintf(stderr,"\t-R neighborhood size to end with [default:"
69           " interval of data]\n");
70   fprintf(stderr,"\t-f factor to increase size [default: 1.2]\n");
71   fprintf(stderr,"\t-s steps to forecast [default: 1]\n");
72   fprintf(stderr,"\t-C width of causality window [default: steps]\n");
73   fprintf(stderr,"\t-o output file name [default: 'datafile.lm']\n");
74   fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
75           "0='only panic messages'\n\t\t"
76           "1='+ input/output messages'\n");
77   fprintf(stderr,"\t-h show these options\n");
78   exit(0);
79 }
80
81 void scan_options(int n,char **in)
82 {
83   char *out;
84
85   if ((out=check_option(in,n,'l','u')) != NULL)
86     sscanf(out,"%lu",&LENGTH);
87   if ((out=check_option(in,n,'x','u')) != NULL)
88     sscanf(out,"%lu",&exclude);
89   if ((out=check_option(in,n,'c','s')) != NULL) {
90     column=out;
91     dimset=1;
92   }
93   if ((out=check_option(in,n,'m','2')) != NULL)
94     sscanf(out,"%u,%u",&dim,&embed);
95   if ((out=check_option(in,n,'d','u')) != NULL)
96     sscanf(out,"%u",&delay);
97   if ((out=check_option(in,n,'i','u')) != NULL)
98     sscanf(out,"%lu",&CLENGTH);
99   if ((out=check_option(in,n,'r','f')) != NULL) {
100     eps0set=1;
101     sscanf(out,"%lf",&EPS0);
102   }
103   if ((out=check_option(in,n,'R','f')) != NULL) {
104     eps1set=1;
105     sscanf(out,"%lf",&EPS1);
106   }
107   if ((out=check_option(in,n,'f','f')) != NULL)
108     sscanf(out,"%lf",&EPSF);
109   if ((out=check_option(in,n,'s','u')) != NULL)
110     sscanf(out,"%u",&STEP);
111   if ((out=check_option(in,n,'C','u')) != NULL) {
112     sscanf(out,"%lu",&causal);
113     causalset=1;
114   }
115   if ((out=check_option(in,n,'V','u')) != NULL)
116     sscanf(out,"%u",&verbosity);
117   if ((out=check_option(in,n,'o','o')) != NULL) {
118     stdo=0;
119     if (strlen(out) > 0)
120       outfile=out;
121   }
122 }
123
124 void make_fit(long act,unsigned long number)
125 {
126   double *si,cast;
127   long i,j;
128   
129   for (i=0;i<dim;i++) {
130     si=series[i];
131     cast=si[found[0]+STEP];
132     for (j=1;j<number;j++)
133       cast += si[found[j]+STEP];
134     cast /= (double)number;
135     error[i] += sqr(cast-series[i][act+STEP]);
136   }
137 }
138
139 int main(int argc,char **argv)
140 {
141   char stdi=0;
142   unsigned long actfound;
143   unsigned long *hfound;
144   long pfound,i,j;
145   unsigned long clength;
146   double interval,min,maxinterval;
147   double epsilon;
148   double **hser;
149   double avfound,*hrms,*hav,sumerror=0.0;
150   FILE *file=NULL;
151
152   if (scan_help(argc,argv))
153     show_options(argv[0]);
154   
155   scan_options(argc,argv);
156 #ifndef OMIT_WHAT_I_DO
157   if (verbosity&VER_INPUT)
158     what_i_do(argv[0],WID_STR);
159 #endif
160
161   if (!causalset)
162     causal=STEP;
163
164   infile=search_datafile(argc,argv,NULL,verbosity);
165   if (infile == NULL)
166     stdi=1;
167
168   if (outfile == NULL) {
169     if (!stdi) {
170       check_alloc(outfile=(char*)calloc(strlen(infile)+4,(size_t)1));
171       sprintf(outfile,"%s.lm",infile);
172     }
173     else {
174       check_alloc(outfile=(char*)calloc((size_t)9,(size_t)1));
175       sprintf(outfile,"stdin.lm");
176     }
177   }
178   if (!stdo)
179     test_outfile(outfile);
180
181   if (column == NULL)
182     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,"",dimset,
183                                       verbosity);
184   else
185     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,column,
186                                       dimset,verbosity);
187   maxinterval=0.0;
188   for (i=0;i<dim;i++) {
189     rescale_data(series[i],LENGTH,&min,&interval);
190     if (interval > maxinterval)
191       maxinterval=interval;
192   }
193   interval=maxinterval;
194
195   check_alloc(list=(long*)malloc(sizeof(long)*LENGTH));
196   check_alloc(found=(unsigned long*)malloc(sizeof(long)*LENGTH));
197   check_alloc(hfound=(unsigned long*)malloc(sizeof(long)*LENGTH));
198   check_alloc(box=(long**)malloc(sizeof(long*)*NMAX));
199   for (i=0;i<NMAX;i++)
200     check_alloc(box[i]=(long*)malloc(sizeof(long)*NMAX));
201   check_alloc(error=(double*)malloc(sizeof(double)*dim));
202   check_alloc(hrms=(double*)malloc(sizeof(double)*dim));
203   check_alloc(hav=(double*)malloc(sizeof(double)*dim));
204   check_alloc(hser=(double**)malloc(sizeof(double*)*dim));
205   
206   if (eps0set)
207     EPS0 /= interval;
208   if (eps1set)
209     EPS1 /= interval;
210
211   clength=(CLENGTH <= LENGTH) ? CLENGTH-STEP : LENGTH-STEP;
212
213   if (!stdo) {
214     file=fopen(outfile,"w");
215     if (verbosity&VER_INPUT)
216       fprintf(stderr,"Opened %s for writing\n",outfile);
217     fprintf(file,"#1. size 2. relative forecast error 3. fraction of points\n"
218             "#4. av neighbors found 5. absolute variance of the points\n");
219   }
220   else {
221     if (verbosity&VER_INPUT)
222       fprintf(stderr,"Writing to stdout\n");
223   }
224
225   for (epsilon=EPS0;epsilon<EPS1*EPSF;epsilon*=EPSF) {
226     pfound=0;
227     for (i=0;i<dim;i++)
228       error[i]=hrms[i]=hav[i]=0.0;
229     avfound=0.0;
230     make_multi_box(series,box,list,LENGTH-STEP,NMAX,dim,
231                    embed,delay,epsilon);
232     for (i=(embed-1)*delay;i<clength;i++) {
233       for (j=0;j<dim;j++)
234         hser[j]=series[j]+i;
235       actfound=find_multi_neighbors(series,box,list,hser,LENGTH,
236                                     NMAX,dim,embed,delay,epsilon,hfound);
237       actfound=exclude_interval(actfound,i-causal+1,i+causal+(embed-1)*delay-1,
238                                 hfound,found);
239       if (actfound > 2*(dim*embed+1)) {
240         make_fit(i,actfound);
241         pfound++;
242         avfound += (double)(actfound-1);
243         for (j=0;j<dim;j++) {
244           hrms[j] += series[j][i+STEP]*series[j][i+STEP];
245           hav[j] += series[j][i+STEP];
246         }
247       }
248     }
249     if (pfound > 1) {
250       sumerror=0.0;
251       for (j=0;j<dim;j++) {
252         hav[j] /= pfound;
253         hrms[j]=sqrt(fabs(hrms[j]/(pfound-1)-hav[j]*hav[j]*pfound/(pfound-1)));
254         error[j]=sqrt(error[j]/pfound)/hrms[j];
255         sumerror += error[j];
256       }
257     }
258     if (stdo) {
259       if (pfound > 1) {
260         fprintf(stdout,"%e %e ",epsilon*interval,sumerror/(double)dim);
261         for (j=0;j<dim;j++)
262           fprintf(stdout,"%e ",error[j]);
263         fprintf(stdout,"%e %e\n",(double)pfound/(clength-(embed-1)*delay),
264                 avfound/pfound);
265         fflush(stdout);
266       }
267     }
268     else {
269       if (pfound > 1) {
270         fprintf(file,"%e %e ",epsilon*interval,sumerror/(double)dim);
271         for (j=0;j<dim;j++)
272           fprintf(file,"%e ",error[j]);
273         fprintf(file,"%e %e\n",(double)pfound/(clength-(embed-1)*delay),
274                 avfound/pfound);
275         fflush(file);
276       }
277     }
278   }
279   if (!stdo)
280     fclose(file);
281
282   free(list);
283   free(hfound);
284   free(error);
285   free(hrms);
286   free(hav);
287   free(hser);
288   for (i=0;i<NMAX;i++)
289     free(box[i]);
290   free(box);
291
292   return 0;
293 }