website/archive/binaries/mac/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