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 }