IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / rbf.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: Mar 11, 2002 */
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 "Fits a RBF-model to the data"
29
30 char *outfile=NULL,stdo=1,MAKECAST=0;
31 char *infile=NULL;
32 char setdrift=1;
33 int DIM=2,DELAY=1,CENTER=10,STEP=1;
34 unsigned int COLUMN=1;
35 unsigned int verbosity=0xff;
36 long CLENGTH=1000;
37 unsigned long LENGTH=ULONG_MAX,INSAMPLE=ULONG_MAX,exclude=0;
38
39 double *series,*coefs;
40 double varianz,interval,min;
41 double **center;
42
43 void show_options(char *progname)
44 {
45   what_i_do(progname,WID_STR);
46   fprintf(stderr," Usage: %s [options]\n",progname);
47   fprintf(stderr," Options:\n");
48   fprintf(stderr,"Everything not being a valid option will be interpreted"
49           " as a possible"
50           " datafile.\nIf no datafile is given stdin is read. Just - also"
51           " means stdin\n");
52   fprintf(stderr,"\t-l # of data to use [default: all from file]\n");
53   fprintf(stderr,"\t-x # of lines to be ignored [default: 0]\n");
54   fprintf(stderr,"\t-c column to read [default: %u]\n",COLUMN);
55   fprintf(stderr,"\t-m embedding dimension [default: %d]\n",DIM);
56   fprintf(stderr,"\t-d delay [default: %d]\n",DELAY);
57   fprintf(stderr,"\t-p number of centers [default: %d]\n",CENTER);
58   fprintf(stderr,"\t-X deactivate drift [default: activated]\n");
59   fprintf(stderr,"\t-s steps to forecast [default: %d]\n",STEP);
60   fprintf(stderr,"\t-n # of points for insample [default: # of data]\n");
61   fprintf(stderr,"\t-L steps to cast [default: none]\n");
62   fprintf(stderr,"\t-o output file name [default: 'datafile'.rbf]\n");
63   fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
64           "0='only panic messages'\n\t\t"
65           "1='+ input/output messages'\n");
66   fprintf(stderr,"\t-h show these options\n");
67   exit(0);
68 }
69
70 void scan_options(int n,char **in)
71 {
72   char *out;
73
74   if ((out=check_option(in,n,'l','u')) != NULL)
75     sscanf(out,"%lu",&LENGTH);
76   if ((out=check_option(in,n,'x','u')) != NULL)
77     sscanf(out,"%lu",&exclude);
78   if ((out=check_option(in,n,'c','u')) != NULL)
79     sscanf(out,"%u",&COLUMN);
80   if ((out=check_option(in,n,'m','u')) != NULL)
81     sscanf(out,"%u",&DIM);
82   if ((out=check_option(in,n,'d','u')) != NULL)
83     sscanf(out,"%u",&DELAY);
84   if ((out=check_option(in,n,'p','u')) != NULL)
85     sscanf(out,"%u",&CENTER);
86   if ((out=check_option(in,n,'X','n')) != NULL)
87     setdrift=0;
88   if ((out=check_option(in,n,'s','u')) != NULL)
89     sscanf(out,"%u",&STEP);
90   if ((out=check_option(in,n,'V','u')) != NULL)
91     sscanf(out,"%u",&verbosity);
92   if ((out=check_option(in,n,'n','u')) != NULL)
93     sscanf(out,"%lu",&INSAMPLE);
94   if ((out=check_option(in,n,'L','u')) != NULL) {
95     MAKECAST=1;
96     sscanf(out,"%lu",&CLENGTH);
97   }
98   if ((out=check_option(in,n,'o','o')) != NULL) {
99     stdo=0;
100     if (strlen(out) > 0)
101       outfile=out;
102   }
103 }
104
105 double avdistance(void)
106 {
107   int i,j,k;
108   double dist=0.0;
109   
110   for (i=0;i<CENTER;i++)
111     for (j=0;j<CENTER;j++)
112       if (i != j)
113         for (k=0;k<DIM;k++)
114           dist += sqr(center[i][k]-center[j][k]);
115
116   return sqrt(dist/(CENTER-1)/CENTER/DIM);
117 }
118
119 double rbf(double *act,double *cen)
120 {
121   static double denum;
122   double r=0;
123   int i;
124
125   denum=2.0*varianz*varianz;
126
127   for (i=0;i<DIM;i++)
128     r += sqr(*(act-i*DELAY)-cen[i]);
129   
130   return exp(-r/denum);
131 }
132
133 void drift(void) 
134 {
135   double *force,h,h1,step=1e-2,step1;
136   int i,j,k,l,d2=DIM;
137
138   check_alloc(force=(double*)malloc(sizeof(double)*d2));
139   for (l=0;l<20;l++) {
140     for (i=0;i<CENTER;i++) {
141       for (j=0;j<d2;j++) {
142         force[j]=0.0;
143         for (k=0;k<CENTER;k++) {
144           if (k != i) {
145             h=center[i][j]-center[k][j];
146             force[j] += h/sqr(h)/fabs(h);
147           }
148         }
149       }
150       h=0.0;
151       for (j=0;j<d2;j++) 
152         h += sqr(force[j]);
153       step1=step/sqrt(h);
154       for (j=0;j<d2;j++) {
155         h1 = step1*force[j];
156         if (((center[i][j]+h1) > -0.1) && ((center[i][j]+h1) < 1.1))
157           center[i][j] += h1;
158       }
159     }
160   }
161   free(force);
162 }
163
164 void make_fit(void)
165 {
166   double **mat,*hcen;
167   double h;
168   int i,j,n,nst;
169
170   check_alloc(mat=(double**)malloc(sizeof(double*)*(CENTER+1)));
171   for (i=0;i<=CENTER;i++)
172     check_alloc(mat[i]=(double*)malloc(sizeof(double)*(CENTER+1)));
173   check_alloc(hcen=(double*)malloc(sizeof(double)*CENTER));
174
175   for (i=0;i<=CENTER;i++) {
176     coefs[i]=0.0;
177     for (j=0;j<=CENTER;j++)
178       mat[i][j]=0.0;
179   }
180
181   for (n=(DIM-1)*DELAY;n<INSAMPLE-STEP;n++) {
182     nst=n+STEP;
183     for (i=0;i<CENTER;i++)
184       hcen[i]=rbf(&series[n],center[i]);
185     coefs[0] += series[nst];
186     mat[0][0] += 1.0;
187     for (i=1;i<=CENTER;i++)
188       mat[i][0] += hcen[i-1];
189     for (i=1;i<=CENTER;i++) {
190       coefs[i] += series[nst]*(h=hcen[i-1]);
191       for (j=1;j<=i;j++)
192         mat[i][j] += h*hcen[j-1];
193     }
194   }
195   
196   h=(double)(INSAMPLE-STEP-(DIM-1)*DELAY);
197   for (i=0;i<=CENTER;i++) {
198     coefs[i] /= h;
199     for (j=0;j<=i;j++) {
200       mat[i][j] /= h;
201       mat[j][i]=mat[i][j];
202     }
203   }
204
205   solvele(mat,coefs,(unsigned int)(CENTER+1));
206
207   for (i=0;i<=CENTER;i++)
208     free(mat[i]);
209   free(mat);
210   free(hcen);
211 }
212
213 double forecast_error(unsigned long i0,unsigned long i1)
214 {
215   int i,n;
216   double h,error=0.0;
217
218   for (n=i0+(DIM-1)*DELAY;n<i1-STEP;n++) {
219     h=coefs[0];
220     for (i=1;i<=CENTER;i++)
221       h += coefs[i]*rbf(&series[n],center[i-1]);
222     error += (series[n+STEP]-h)*(series[n+STEP]-h);
223   }
224   
225   return sqrt(error/(i1-i0-STEP-(DIM-1)*DELAY));
226 }
227
228 void make_cast(FILE *out)
229 {
230   double *cast,new_el;
231   int i,n,dim;
232   
233   dim=(DIM-1)*DELAY;
234   check_alloc(cast=(double*)malloc(sizeof(double)*(dim+1)));
235   for (i=0;i<=dim;i++)
236     cast[i]=series[LENGTH-1-dim+i];
237   
238   for (n=0;n<CLENGTH;n++) {
239     new_el=coefs[0];
240     for (i=1;i<=CENTER;i++)
241       new_el += coefs[i]*rbf(&cast[dim],center[i-1]);
242     fprintf(out,"%e\n",new_el*interval+min);
243     for (i=0;i<dim;i++)
244       cast[i]=cast[i+1];
245     cast[dim]=new_el;
246   }
247 }
248
249 int main(int argc,char **argv)
250 {
251   char stdi=0;
252   int i,j,cstep;
253   double sigma,av;
254   FILE *file=NULL;
255
256   if (scan_help(argc,argv))
257     show_options(argv[0]);
258
259   scan_options(argc,argv);
260 #ifndef OMIT_WHAT_I_DO
261   if (verbosity&VER_INPUT)
262     what_i_do(argv[0],WID_STR);
263 #endif
264
265   infile=search_datafile(argc,argv,&COLUMN,verbosity);
266   if (infile == NULL) 
267     stdi=1;
268
269   if (outfile == NULL) {
270     if (!stdi) {
271       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
272       strcpy(outfile,infile);
273       strcat(outfile,".rbf");
274     }
275     else {
276       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
277       strcpy(outfile,"stdin.rbf");
278     }
279   }
280   if (!stdo)
281     test_outfile(outfile);
282
283   series=(double*)get_series(infile,&LENGTH,exclude,COLUMN,verbosity);
284   rescale_data(series,LENGTH,&min,&interval);
285   variance(series,LENGTH,&av,&varianz);
286
287   if (INSAMPLE > LENGTH)
288     INSAMPLE=LENGTH;
289   
290   if (CENTER > LENGTH) 
291     CENTER = LENGTH;
292   
293   if (MAKECAST)
294     STEP=1;
295   
296   check_alloc(coefs=(double*)malloc(sizeof(double)*(CENTER+1)));
297   check_alloc(center=(double**)malloc(sizeof(double*)*CENTER));
298   for (i=0;i<CENTER;i++)
299     check_alloc(center[i]=(double*)malloc(sizeof(double)*DIM));
300   
301   cstep=LENGTH-1-(DIM-1)*DELAY;
302   for (i=0;i<CENTER;i++)
303     for (j=0;j<DIM;j++)
304       center[i][j]=series[(DIM-1)*DELAY-j*DELAY+(i*cstep)/(CENTER-1)];
305
306   if (setdrift)
307     drift();
308   varianz=avdistance();
309   make_fit();
310
311   if (!stdo) {
312     file=fopen(outfile,"w");
313     if (verbosity&VER_INPUT)
314       fprintf(stderr,"Opened %s for writing\n",outfile);
315     fprintf(file,"#Center points used:\n");
316     for (i=0;i<CENTER;i++) {
317       fprintf(file,"#");
318       for (j=0;j<DIM;j++)
319         fprintf(file," %e",center[i][j]*interval+min);
320       fprintf(file,"\n");
321     }
322     fprintf(file,"#variance= %e\n",varianz*interval);
323     fprintf(file,"#Coefficients:\n");
324     fprintf(file,"#%e\n",coefs[0]*interval+min);
325     for (i=1;i<=CENTER;i++)
326       fprintf(file,"#%e\n",coefs[i]*interval);
327   }
328   else {
329     if (verbosity&VER_INPUT)
330       fprintf(stderr,"Writing to stdout\n");
331     fprintf(stdout,"#Center points used:\n");
332     for (i=0;i<CENTER;i++) {
333       fprintf(stdout,"#");
334       for (j=0;j<DIM;j++)
335         fprintf(stdout," %e",center[i][j]*interval+min);
336       fprintf(stdout,"\n");
337     }
338     fprintf(stdout,"#variance= %e\n",varianz*interval);
339     fprintf(stdout,"#Coefficients:\n");
340     fprintf(stdout,"#%e\n",coefs[0]*interval+min);
341     for (i=1;i<=CENTER;i++)
342       fprintf(stdout,"#%e\n",coefs[i]*interval);
343   }
344   av=sigma=0.0;
345   for (i=0;i<INSAMPLE;i++) {
346     av += series[i];
347     sigma += series[i]*series[i];
348   }
349   av /= INSAMPLE;
350   sigma=sqrt(fabs(sigma/INSAMPLE-av*av));
351   if (!stdo)
352     fprintf(file,"#insample error= %e\n",forecast_error(0LU,INSAMPLE)/sigma);
353   else
354     fprintf(stdout,"#insample error= %e\n",forecast_error(0LU,INSAMPLE)/sigma);
355
356   if (INSAMPLE < LENGTH) {
357     av=sigma=0.0;
358     for (i=INSAMPLE;i<LENGTH;i++) {
359       av += series[i];
360       sigma += series[i]*series[i];
361     }
362     av /= (LENGTH-INSAMPLE);
363     sigma=sqrt(fabs(sigma/(LENGTH-INSAMPLE)-av*av));
364     if (!stdout)
365       fprintf(file,"#out of sample error= %e\n",
366               forecast_error(INSAMPLE,LENGTH)/sigma);
367     else
368       fprintf(stdout,"#out of sample error= %e\n",
369               forecast_error(INSAMPLE,LENGTH)/sigma);
370   }
371
372   if (MAKECAST) {
373     if (!stdo)
374       make_cast(file);
375     else
376       make_cast(stdout);
377   }
378
379   if (!stdo)
380     fclose(file);
381
382   return 0;
383 }