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 }