IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / src / disembl / Tisean_3.0.1 / source_c / nstat_z.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 <limits.h>
25 #include <ctype.h>
26 #include "routines/tsa.h"
27
28 #define WID_STR "Tests for nonstationarity by means of the average\n\t\
29 forecast error for a zeroth order fit"
30
31
32 #ifndef _MATH_H
33 #include <math.h>
34 #endif
35
36 /*number of boxes for the neighbor search algorithm*/
37 #define NMAX 128
38
39 unsigned int nmax=(NMAX-1);
40 long **box,*list;
41 unsigned long *found;
42 double *series,*series1,*series2;
43 double interval,min,epsilon;
44
45 char epsset=0,causalset=0;
46 char *infile=NULL;
47 char *outfile=NULL,stdo=1,centerset=0;
48 char *firstwindow,*secondwindow,**window;
49 unsigned int COLUMN=1,pieces;
50 unsigned int verbosity=0xff;
51 int DIM=3,DELAY=1,MINN=30,STEP=1;
52 int firstoffset= -1,secondoffset= -1;
53 double EPS0=1.e-3,EPSF=1.2;
54 unsigned long LENGTH=ULONG_MAX,exclude=0,center,causal;
55
56 void show_options(char *progname)
57 {
58   what_i_do(progname,WID_STR);
59   fprintf(stderr," Usage: %s -# [other 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 column to read [default: 1]\n");
68   fprintf(stderr,"\t-m embedding dimension [default: 3]\n");
69   fprintf(stderr,"\t-d delay [default: 1]\n");
70   fprintf(stderr,"\t-# # of pieces [no default]\n");
71   fprintf(stderr,"\t-1 which pieces for the first window "
72           "[default: 1-pieces]\n");
73   fprintf(stderr,"\t-2 which pieces for the second window "
74           "[default: 1-pieces]\n");
75   fprintf(stderr,"\t-n # of reference points in the window [default: all]\n");
76   fprintf(stderr,"\t-k minimal number of neighbors for the fit "
77           "[default: 30]\n");
78   fprintf(stderr,"\t-r neighborhoud size to start with "
79           "[default: (data interval)/1000]\n");
80   fprintf(stderr,"\t-f factor to increase size [default: 1.2]\n");
81   fprintf(stderr,"\t-s steps to forecast [default: 1]\n");
82   fprintf(stderr,"\t-C width of causality window [default: steps]\n");
83   fprintf(stderr,"\t-o output file [default: 'datafile.nsz',"
84           " without -o: stdout]\n");
85   fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
86           "0='only panic messages'\n\t\t"
87           "1='+ input/output messages'\n");
88   fprintf(stderr,"\t-h show these options\n");
89   fprintf(stderr,"\n\t The -# option has to be set\n");
90   exit(0);
91 }
92
93 void parse_minus(char *str,char *array,char *wopt)
94 {
95   int cm=0,i,strl,n1,n2;
96   
97   strl=strlen(str);
98   for (i=0;i<strl;i++)
99     if (str[i] == '-')
100       cm++;
101   if (cm > 1) {
102     fprintf(stderr,"Invalid string for the %s option! "
103             "Please consult the help-page\n",wopt);
104     exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
105   }
106   if (cm == 0) {
107     sscanf(str,"%d",&n1);
108     n1--;
109     if (n1 < 0) {
110       fprintf(stderr,"Numbers in %s option must be larger than 0!\n",wopt);
111       exit(NSTAT_Z__NOT_UNSIGNED_FOR_OPTION);
112     }
113     if (n1 >= pieces) {
114       fprintf(stderr,"Numbers in %s option must be smaller than %u!\n",wopt,
115               pieces);
116       exit(NSTAT_Z__TOO_LARGE_FOR_OPTION);
117     }
118     array[n1]=1;
119   }
120   else {
121     sscanf(str,"%d-%d",&n1,&n2);
122     n1--;
123     n2--;
124     if ((n1 < 0) || (n2 < 0)) {
125       fprintf(stderr,"Numbers in %s option must be larger than 0!\n",wopt);
126       exit(NSTAT_Z__NOT_UNSIGNED_FOR_OPTION);
127     }
128     if ((n1 >= pieces) || (n2 >= pieces)) {
129       fprintf(stderr,"Numbers in %s option must be smaller than %u!\n",wopt,
130               pieces+1);
131       exit(NSTAT_Z__TOO_LARGE_FOR_OPTION);
132     }
133     if (n2 < n1) {
134       i=n1;
135       n1=n2;
136       n2=i;
137     }
138     for (i=n1;i<=n2;i++)
139       array[i]=1;
140   }
141 }
142
143 void parse_comma(char *str,char *array,char *wopt)
144 {
145   unsigned int strl,i,cp=1,which,iwhich;
146   char **hstr;
147
148   strl=strlen(str);
149   for (i=0;i<strl;i++)
150     if (str[i] == ',')
151       cp++;
152
153   if (cp == 1) {
154     parse_minus(str,array,wopt);
155     return ;
156   }
157   
158   check_alloc(hstr=(char**)malloc(sizeof(char*)*cp));
159   for (i=0;i<cp;i++)
160     check_alloc(hstr[i]=(char*)calloc(strl,1));
161   
162   which=iwhich=0;
163   for (i=0;i<strl;i++) {
164     if (str[i] != ',')
165       hstr[which][iwhich++]=str[i];
166     else {
167       which++;
168       iwhich=0;
169     }
170   }
171   for (i=0;i<cp;i++) {
172     if (hstr[i][0] == '\0') {
173       fprintf(stderr,"Invalid string for the %s option! "
174               "Please consult the help-page\n",wopt);
175       exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
176     }
177     if (!isdigit(hstr[i][strlen(hstr[i])-1])) {
178       fprintf(stderr,"Invalid string for the %s option! "
179               "Please consult the help-page\n",wopt);
180       exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
181     }
182     parse_minus(hstr[i],array,wopt);
183   }
184   for (i=0;i<cp;i++)
185     free(hstr[i]);
186   free(hstr);
187 }
188
189 void parse_out(char *str,char *array,char *which)
190 {
191   unsigned int i;
192   char test;
193
194   for (i=0;i<pieces;i++)
195     array[i]=0;
196   
197   for (i=0;i<strlen(str);i++) {
198     test= (str[i] == '-') || (str[i] == ',') || isdigit(str[i]);
199     if (!test) {
200       fprintf(stderr,"Invalid string for the %s option! "
201               "Please consult the help-page\n",which);
202       exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
203     }
204   }
205   if (!isdigit(str[strlen(str)-1])) {
206     fprintf(stderr,"Invalid string for the %s option! "
207             "Please consult the help-page\n",which);
208     exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
209   }
210   parse_comma(str,array,which);
211 }
212
213 void parse_offset(char *str,int *iwhich,char *array,char *which) 
214 {
215   int i,strl;
216   
217   if (str[0] != '+')
218     return;
219   strl=strlen(str);
220   for (i=1;i<strl;i++)
221     if (!isdigit(str[i])) {
222       fprintf(stderr,"Invalid string for the %s option! "
223               "Please consult the help-page\n",which);
224       exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
225     }
226   sscanf(str,"+%d",iwhich);
227   for (i=0;i<pieces;i++)
228     array[i]=0;
229 }
230       
231 void scan_options(int n,char **in)
232 {
233   unsigned int i;
234   char *out,piecesset=0;
235
236   if ((out=check_option(in,n,'l','u')) != NULL)
237     sscanf(out,"%lu",&LENGTH);
238   if ((out=check_option(in,n,'x','u')) != NULL)
239     sscanf(out,"%lu",&exclude);
240   if ((out=check_option(in,n,'c','u')) != NULL)
241     sscanf(out,"%u",&COLUMN);
242   if ((out=check_option(in,n,'m','u')) != NULL)
243     sscanf(out,"%u",&DIM);
244   if ((out=check_option(in,n,'d','u')) != NULL)
245     sscanf(out,"%u",&DELAY);
246   if ((out=check_option(in,n,'V','u')) != NULL)
247     sscanf(out,"%u",&verbosity);
248   if ((out=check_option(in,n,'#','u')) != NULL) {
249     sscanf(out,"%u",&pieces);
250     if (pieces < 1)
251       pieces=1;
252     piecesset=1;
253     check_alloc(firstwindow=(char*)malloc(pieces));
254     check_alloc(secondwindow=(char*)malloc(pieces));
255     for (i=0;i<pieces;i++)
256       firstwindow[i]=secondwindow[i]=1;
257     check_alloc(window=(char**)malloc(sizeof(char*)*pieces));
258     for (i=0;i<pieces;i++)
259       check_alloc(window[i]=(char*)malloc(pieces));
260   }
261   if (!piecesset) {
262     fprintf(stderr,"\tThe -# option wasn't set. Please add it!\n");
263     exit(NSTAT_Z__OPTION_NOT_SET);
264   }
265   if ((out=check_option(in,n,'1','s')) != NULL) {
266     parse_offset(out,&firstoffset,firstwindow,"-1");
267     if (firstoffset == -1)
268       parse_out(out,firstwindow,"-1");
269   }
270   if ((out=check_option(in,n,'2','s')) != NULL) {
271     parse_offset(out,&secondoffset,secondwindow,"-2");
272     if (secondoffset == -1)
273       parse_out(out,secondwindow,"-2");
274   }
275   if ((out=check_option(in,n,'n','u')) != NULL) {
276     sscanf(out,"%lu",&center);
277     centerset=1;
278   }
279   if ((out=check_option(in,n,'k','u')) != NULL)
280     sscanf(out,"%u",&MINN);
281   if ((out=check_option(in,n,'r','f')) != NULL) {
282     epsset=1;
283     sscanf(out,"%lf",&EPS0);
284   }
285   if ((out=check_option(in,n,'f','f')) != NULL)
286     sscanf(out,"%lf",&EPSF);
287   if ((out=check_option(in,n,'s','u')) != NULL)
288     sscanf(out,"%u",&STEP);
289   if ((out=check_option(in,n,'C','u')) != NULL) {
290     sscanf(out,"%lu",&causal);
291     causalset=1;
292   }
293   if ((out=check_option(in,n,'o','o')) != NULL) {
294     stdo=0;
295     if (strlen(out) > 0)
296       outfile=out;
297   }
298 }
299
300 double make_fit(long act,unsigned long number)
301 {
302   double casted=0.0,*help;
303   int i;
304   
305   help=series1+STEP;
306   for (i=0;i<number;i++) {
307     casted += help[found[i]];
308   }
309   casted /= number;
310
311   return sqr(casted-series2[act+STEP]);
312 }
313
314 int main(int argc,char **argv)
315 {
316   char stdi=0;
317   char alldone,*done,sdone;
318   long i,first,second,pstart;
319   unsigned long *hfound;
320   unsigned long actfound;
321   unsigned long clength;
322   double *rms,av,error;
323   FILE *file=NULL;
324
325   if (scan_help(argc,argv))
326     show_options(argv[0]);
327   
328   scan_options(argc,argv);
329
330   if (!causalset)
331     causal=STEP;
332
333 #ifndef OMIT_WHAT_I_DO
334   if (verbosity&VER_INPUT)
335     what_i_do(argv[0],WID_STR);
336 #endif
337
338   infile=search_datafile(argc,argv,&COLUMN,verbosity);
339   if (infile == NULL)
340     stdi=1;
341
342   if (outfile == NULL) {
343     if (!stdi) {
344       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
345       sprintf(outfile,"%s.nsz",infile);
346     }
347     else {
348       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
349       sprintf(outfile,"stdin.nsz");
350     }
351   }
352   if (!stdo)
353     test_outfile(outfile);
354
355   series=(double*)get_series(infile,&LENGTH,exclude,COLUMN,verbosity);
356
357   rescale_data(series,LENGTH,&min,&interval);
358   
359   check_alloc(list=(long*)malloc(sizeof(long)*LENGTH));
360   check_alloc(found=(unsigned long*)malloc(sizeof(long)*LENGTH));
361   check_alloc(hfound=(unsigned long*)malloc(sizeof(long)*LENGTH));
362   check_alloc(done=(char*)malloc(sizeof(char)*LENGTH));
363   check_alloc(box=(long**)malloc(sizeof(long*)*NMAX));
364
365   for (i=0;i<NMAX;i++)
366     check_alloc(box[i]=(long*)malloc(sizeof(long)*NMAX));
367     
368   if (epsset)
369     EPS0 /= interval;
370
371   clength=(LENGTH-(DIM-1)*DELAY)/pieces;
372   if ((clength-(DIM-1)*DELAY-STEP) < MINN) {
373     fprintf(stderr,"You chose too many pieces and will never find enough"
374             " neighbors!\n");
375     exit(NSTAT_Z__TOO_MANY_PIECES);
376   }
377   check_alloc(rms=(double*)malloc(sizeof(double)*pieces));
378   for (i=0;i<pieces;i++) {
379     series1=series+i*clength;
380     variance(series1,clength,&av,&rms[i]);
381   }
382   
383   pstart=(DIM-1)*DELAY;
384   if (!centerset)
385     center=clength-STEP;
386   else
387     center=(center < (clength-STEP-pstart)) ? center : clength-STEP-pstart;
388
389   if (stdo) {
390     if (verbosity&VER_INPUT)
391       fprintf(stderr,"Writing to stdout\n");
392   }
393   else {
394     file=fopen(outfile,"w");
395     if (verbosity&VER_INPUT)
396       fprintf(stderr,"Opened %s for writing\n",outfile);
397   }
398   for (first=0;first<pieces;first++)
399     for (second=0;second<pieces;second++)
400       window[first][second]=firstwindow[first]&&secondwindow[second];
401   if (firstoffset != -1) {
402     for (second=0;second<pieces;second++)
403       for (first=second-firstoffset;first<=second+firstoffset;first++)
404         if ((first >= 0) && (first < pieces))
405           window[first][second]=secondwindow[second];
406   }
407   if (secondoffset != -1) {
408     for (first=0;first<pieces;first++)
409       for (second=first-secondoffset;second<=first+secondoffset;second++)
410         if ((second >= 0) && (second < pieces))
411           window[first][second]=firstwindow[first];
412   }
413
414   free(firstwindow);
415   free(secondwindow);
416
417   for (first=0;first<pieces;first++) {
418     sdone=0;
419     for (second=0;second<pieces;second++) {
420       if (window[first][second]) {
421         sdone=1;
422         series1=series+first*clength;
423         series2=series+second*clength;
424         for (i=0;i<LENGTH;i++)
425           done[i]=0;
426         alldone=0;
427         epsilon=EPS0/EPSF;
428         error=0.0;
429         while (!alldone) {
430           alldone=1;
431           epsilon*=EPSF;
432           make_box(series1,box,list,clength-STEP,NMAX,(unsigned int)DIM,
433                    (unsigned int)DELAY,epsilon);
434           for (i=pstart;i<pstart+center;i++)
435             if (!done[i]) {
436               actfound=find_neighbors(series1,box,list,series2+i,clength,NMAX,
437                                       (unsigned int)DIM,(unsigned int)DELAY,
438                                       epsilon,hfound);
439               actfound=exclude_interval(actfound,i-causal+1,
440                                         i+causal+pstart-1,hfound,found);
441               if (actfound >= MINN) {
442                 error += make_fit(i,actfound);
443                 done[i]=1;
444               }
445               alldone &= done[i];
446             }
447         }
448         if (stdo)
449           fprintf(stdout,"%ld %ld %e\n",first+1,second+1,
450                   sqrt(error/center)/rms[second]);
451         else {
452           fprintf(file,"%ld %ld %e\n",first+1,second+1,
453                   sqrt(error/center)/rms[second]);
454           fflush(file);
455         }
456       }
457     }
458     if (sdone) {
459       if (stdo)
460         fprintf(stdout,"\n");
461       else
462         fprintf(file,"\n");
463     }
464   }
465   
466   if (!stdo)
467     fclose(file);
468
469   if (outfile != NULL)
470     free(outfile);
471   free(list);
472   free(found);
473   free(hfound);
474   free(done);
475   for (i=0;i<NMAX;i++)
476     free(box[i]);
477   free(box);
478   for (i=0;i<pieces;i++)
479     free(window[i]);
480   free(window);
481   free(rms);
482   free(series);
483
484   return 0;
485 }
486