IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / polynomp.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 5, 2004 */
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include <math.h>
25 #include <limits.h>
26 #include <time.h>
27 #include "routines/tsa.h"
28
29 #define WID_STR "Fits a polynomial to the data."
30
31 char *outfile=NULL,stdo=1;
32 char *parin=NULL,*infile=NULL;
33 unsigned long length=ULONG_MAX,insample=ULONG_MAX,exclude=0;
34 unsigned long plength=UINT_MAX;
35 unsigned long step=1000;
36 unsigned int column=1,dim=2,delay=1,down_to=1;
37 unsigned int **order;
38 unsigned int verbosity=0xff;
39 double *series,*param;
40
41 void show_options(char *progname)
42 {
43   what_i_do(progname,WID_STR);
44   fprintf(stderr,"Usage: %s [Options]\n",progname);
45   fprintf(stderr,"Options:\n");
46   fprintf(stderr,"Everything not being a valid option will be interpreted"
47           " as a possible"
48           " datafile.\nIf no datafile is given stdin is read. Just - also"
49           " means stdin\n");
50   fprintf(stderr,"\t-l # of data to use [default: whole file]\n");
51   fprintf(stderr,"\t-x # of lines to ignore [default: %lu]\n",exclude);
52   fprintf(stderr,"\t-c column to read [default: %u]\n",column);
53   fprintf(stderr,"\t-m embedding dimension [default: %u]\n",dim);
54   fprintf(stderr,"\t-d delay [default: %u]\n",delay);
55   fprintf(stderr,"\t-n insample data [default: all]\n");
56   fprintf(stderr,"\t-L length of forecasted series [default: %lu]\n",step);
57   fprintf(stderr,"\t-p name of parameter file [default: parameter.pol]\n");
58   fprintf(stderr,"\t-o output file name [default: 'datafile'.pbf]\n");
59   fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
60           "0='only panic messages'\n\t\t"
61           "1='+ input/output messages'\n");
62   fprintf(stderr,"\t-h show these options\n");
63   exit(0);
64 }
65
66 void scan_options(int n,char **in)
67 {
68   char *out;
69
70   if ((out=check_option(in,n,'l','u')) != NULL)
71     sscanf(out,"%lu",&length);
72   if ((out=check_option(in,n,'x','u')) != NULL)
73     sscanf(out,"%lu",&exclude);
74   if ((out=check_option(in,n,'c','u')) != NULL)
75     sscanf(out,"%u",&column);
76   if ((out=check_option(in,n,'m','u')) != NULL)
77     sscanf(out,"%u",&dim);
78   if ((out=check_option(in,n,'d','u')) != NULL)
79     sscanf(out,"%u",&delay);
80   if ((out=check_option(in,n,'V','u')) != NULL)
81     sscanf(out,"%u",&verbosity);
82   if ((out=check_option(in,n,'n','u')) != NULL)
83     sscanf(out,"%lu",&insample);
84   if ((out=check_option(in,n,'L','u')) != NULL)
85     sscanf(out,"%lu",&step);
86   if ((out=check_option(in,n,'p','s')) != NULL)
87     parin=out;
88   if ((out=check_option(in,n,'o','o')) != NULL) {
89     stdo=0;
90     if (strlen(out) > 0)
91       outfile=out;
92   } 
93 }
94
95 double polynom(unsigned long act,unsigned int which)
96 {
97   unsigned int i,j;
98   double ret=1.0,h;
99   
100   for (i=0;i<dim;i++) {
101     h=series[act-i*delay];
102     for (j=0;j<order[which][i];j++)
103       ret *= h;
104   }
105   
106   return ret;
107 }
108
109 void make_fit(void)
110 {
111   double **mat,*vec;
112   double h;
113   unsigned long n,hn;
114   unsigned int i,j;
115
116   check_alloc(vec=(double*)malloc(sizeof(double)*plength));
117   check_alloc(mat=(double**)malloc(sizeof(double*)*plength));
118   for (i=0;i<plength;i++)
119     check_alloc(mat[i]=(double*)malloc(sizeof(double)*plength));
120
121   for (i=0;i<plength;i++) {
122     vec[i]=0.0;
123     for (j=0;j<plength;j++)
124       mat[i][j]=0.0;
125   }
126   
127   for (n=(dim-1)*delay;n<insample-1;n++) {
128     hn=n+1;
129     for (i=0;i<plength;i++) {
130       vec[i] += series[hn]*(h=polynom(n,i));
131       for (j=i;j<plength;j++)
132         mat[i][j] += polynom(n,j)*h;
133     }
134   }
135   for (i=0;i<plength;i++) {
136     vec[i] /= (insample-(dim-1)*delay-1);
137     for (j=i;j<plength;j++)
138       mat[j][i]=(mat[i][j]/=(insample-(dim-1)*delay)-1);
139   }
140   
141   solvele(mat,vec,plength);
142
143   for (i=0;i<plength;i++)
144     param[i]=vec[i];
145
146   free(vec);
147   for (i=0;i<plength;i++)
148     free(mat[i]);
149   free(mat);
150 }
151
152 double forecast_error(unsigned long i0,unsigned long i1)
153 {
154   unsigned int i;
155   unsigned long n;
156   double h,error=0.0;
157
158   for (n=i0+(dim-1)*delay;n<i1-1;n++) {
159     h=0.0;
160     for (i=0;i<plength;i++)
161       h += param[i]*polynom(n,i);
162     error += (series[n+1]-h)*(series[n+1]-h);
163   }
164
165   return sqrt(error/(i1-i0-(dim-1)*delay-1));
166 }
167
168 void make_cast(FILE *fcast)
169 {
170   int i,j,hi;
171   unsigned int k;
172   double casted;
173   
174   for (i=0;i<=(dim-1)*delay;i++)
175     series[i]=series[length-(dim-1)*delay+i-1];
176
177   hi=(dim-1)*delay;
178   for (i=1;i<=step;i++) {
179     casted=0.0;
180     for (k=0;k<plength;k++)
181       casted += param[k]*polynom((unsigned long)((dim-1)*delay),k);
182     if (!stdo) {
183       fprintf(fcast,"%e\n",casted);
184       fflush(fcast);
185     }
186     else {
187       fprintf(stdout,"%e\n",casted);
188       fflush(stdout);
189     }
190     for (j=0;j<(dim-1)*delay;j++)
191       series[j]=series[j+1];
192     series[hi]=casted;
193   }
194 }
195
196 int main(int argc,char **argv)
197 {
198   int i,j;
199   char stdi=0,oose=1;
200   double **dummy,withalli,withallo;
201   double av,varianz;
202   FILE *file;
203
204   if (scan_help(argc,argv))
205     show_options(argv[0]);
206
207   scan_options(argc,argv);
208 #ifndef OMIT_WHAT_I_DO
209   if (verbosity&VER_INPUT)
210     what_i_do(argv[0],WID_STR);
211 #endif
212
213   infile=search_datafile(argc,argv,&column,verbosity);
214   if (infile == NULL)
215     stdi=1;
216   
217   if (outfile == NULL) {
218     if (!stdi) {
219       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
220       sprintf(outfile,"%s.pbf",infile);
221     }
222     else {
223       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
224       sprintf(outfile,"stdin.pbf");
225     }
226   }
227   if (!stdo)
228     test_outfile(outfile);
229
230   if (parin == NULL) {
231     check_alloc(parin=(char*)calloc((size_t)14,(size_t)1));
232     sprintf(parin,"parameter.pol");
233   }
234   file=fopen(parin,"r");
235   if (file == NULL) {
236     fprintf(stderr,"File %s does not exist. Exiting!\n",parin);
237     exit(POLYNOMP__WRONG_PARAMETER_FILE);
238   }
239   fclose(file);
240
241   dummy=(double**)get_multi_series(parin,&plength,0LU,
242                                    &dim,"",(char)"1",verbosity);
243   
244   check_alloc(order=(unsigned int**)malloc(sizeof(int*)*plength));
245   for (i=0;i<plength;i++) {
246     check_alloc(order[i]=(unsigned int*)malloc(sizeof(int)*dim));
247     for (j=0;j<dim;j++)
248       order[i][j]=(unsigned int)dummy[j][i];
249   }
250
251   series=(double*)get_series(infile,&length,exclude,column,verbosity);
252   variance(series,length,&av,&varianz);
253
254   if (insample >= length) {
255     insample=length;
256     oose=0;
257   }
258
259   check_alloc(param=(double*)malloc(sizeof(double)*plength));
260
261   make_fit();
262   withalli=forecast_error(0LU,insample);
263   withallo=0.0;
264   if (oose)
265     withallo=forecast_error(insample+1,length);
266
267   if (stdo) {
268     if (verbosity&VER_INPUT)
269       fprintf(stderr,"Writing to stdout\n");
270     fprintf(stdout,"#FCE: %e %e\n",withalli/varianz,withallo/varianz);
271     for (i=0;i<plength;i++) {
272       fprintf(stdout,"# ");
273       for (j=0;j<dim;j++)
274         fprintf(stdout,"%u ",order[i][j]);
275       fprintf(stdout,"%e\n",param[i]);
276     }
277     fflush(stdout);
278   }
279   else {
280     file=fopen(outfile,"w");
281     if (verbosity&VER_INPUT)
282       fprintf(stderr,"Opened %s for writing\n",outfile);
283     fprintf(file,"#FCE: %e %e\n",withalli/varianz,withallo/varianz);
284     for (i=0;i<plength;i++) {
285       fprintf(file,"# ");
286       for (j=0;j<dim;j++)
287         fprintf(file,"%u ",order[i][j]);
288       fprintf(file,"%e\n",param[i]);
289     }
290     fflush(file);
291   }
292   
293   make_cast(file);
294
295   if (!stdo)
296     fclose(file);
297
298   return 0;
299 }