binaries/src/disembl/Tisean_3.0.1/source_c/ar-model.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*/
  21 /*Changes:
  22   Jun 24, 2005: Output average error for multivariate data
  23   Nov 25, 2005: Handle model order = 0
  24   Jan 31, 2006: Add verbosity 4 to print data+residuals
  25  */
  26 #include <stdio.h>
  27 #include <stdlib.h>
  28 #include <string.h>
  29 #include <limits.h>
  30 #include <math.h>
  31 #include "routines/tsa.h"
  32
  33 #define WID_STR "Fits an multivariate AR model to the data and gives\
  34  the coefficients\n\tand the residues (or an iterated model)"
  35
  36 unsigned long length=ULONG_MAX,exclude=0;
  37 unsigned int dim=1,poles=1,ilength;
  38 unsigned int verbosity=1;
  39 char *outfile=NULL,*column=NULL,stdo=1,dimset=0,run_model=0;
  40 char *infile=NULL;
  41 double **series,*my_average;
  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 length of file [default is whole file]\n");
  53   fprintf(stderr,"\t-x # of lines to be ignored [default is 0]\n");
  54   fprintf(stderr,"\t-m dimension [default is 1]\n");
  55   fprintf(stderr,"\t-c columns to read [default is 1,...,dimension]\n");
  56   fprintf(stderr,"\t-p #order of AR-Fit [default is 1]\n");
  57   fprintf(stderr,"\t-s length of iterated model [default no iteration]\n");
  58   fprintf(stderr,"\t-o output file name [default is 'datafile'.ar]\n");
  59   fprintf(stderr,"\t-V verbosity level [default is 1]\n\t\t"
  60           "0='only panic messages'\n\t\t"
  61           "1='+ input/output messages'\n\t\t"
  62           "2='+ print residuals though iterating a model'\n\t\t"
  63           "4='+ print original data plus residuals'\n");
  64   fprintf(stderr,"\t-h show these options\n\n");
  65   exit(0);
  66 }
  67
  68 void scan_options(int argc,char **argv)
  69 {
  70   char *out;
  71
  72   if ((out=check_option(argv,argc,'p','u')) != NULL) {
  73     sscanf(out,"%u",&poles);
  74     if (poles < 1) {
  75       fprintf(stderr,"The order should at least be one!\n");
  76       exit(127);
  77     }
  78   }
  79   if ((out=check_option(argv,argc,'l','u')) != NULL)
  80     sscanf(out,"%lu",&length);
  81   if ((out=check_option(argv,argc,'x','u')) != NULL)
  82     sscanf(out,"%lu",&exclude);
  83   if ((out=check_option(argv,argc,'m','u')) != NULL) {
  84     sscanf(out,"%u",&dim);
  85     dimset=1;
  86   }
  87   if ((out=check_option(argv,argc,'c','u')) != NULL)
  88     column=out;
  89   if ((out=check_option(argv,argc,'V','u')) != NULL)
  90     sscanf(out,"%u",&verbosity);
  91   if ((out=check_option(argv,argc,'s','u')) != NULL) {
  92     sscanf(out,"%u",&ilength);
  93     run_model=1;
  94   }
  95   if ((out=check_option(argv,argc,'o','o')) != NULL) {
  96     stdo=0;
  97     if (strlen(out) > 0)
  98       outfile=out;
  99   }
 100 }
 101
 102 void set_averages_to_zero(void)
 103 {
 104   double var;
 105   long i,j;
 106
 107   for (i=0;i<dim;i++) {
 108     variance(series[i],length,&my_average[i],&var);
 109     for (j=0;j<length;j++)
 110       series[i][j] -= my_average[i];
 111   }
 112 }
 113
 114 double** build_matrix(double **mat)
 115 {
 116   long n,i1,j1,i2,j2,hi,hj;
 117   double norm;
 118
 119   norm=1./((double)length-(double)poles);
 120
 121   for (i1=0;i1<dim;i1++)
 122     for (i2=0;i2<poles;i2++) {
 123       hi=i1*poles+i2;
 124       for (j1=0;j1<dim;j1++)
 125         for (j2=0;j2<poles;j2++) {
 126           hj=j1*poles+j2;
 127           mat[hi][hj]=0.0;
 128           for (n=poles-1;n<length-1;n++)
 129             mat[hi][hj] += series[i1][n-i2]*series[j1][n-j2];
 130           mat[hi][hj] *= norm;
 131         }
 132     }
 133
 134   return invert_matrix(mat,(unsigned int)(dim*poles));
 135 }
 136
 137 void build_vector(double *vec,long comp)
 138 {
 139   long i1,i2,hi,n;
 140   double norm;
 141
 142   norm=1./((double)length-(double)poles);
 143
 144   for (i1=0;i1<poles*dim;i1++)
 145     vec[i1]=0.0;
 146
 147   for (i1=0;i1<dim;i1++)
 148     for (i2=0;i2<poles;i2++) {
 149       hi=i1*poles+i2;
 150       for (n=poles-1;n<length-1;n++)
 151         vec[hi] += series[comp][n+1]*series[i1][n-i2];
 152       vec[hi] *= norm;
 153     }
 154 }
 155
 156 double* multiply_matrix_vector(double **mat,double *vec)
 157 {
 158   long i,j;
 159   double *new_vec;
 160
 161   check_alloc(new_vec=(double*)malloc(sizeof(double)*poles*dim));
 162
 163   for (i=0;i<poles*dim;i++) {
 164     new_vec[i]=0.0;
 165     for (j=0;j<poles*dim;j++)
 166       new_vec[i] += mat[i][j]*vec[j];
 167   }
 168   return new_vec;
 169 }
 170
 171 double* make_residuals(double **diff,double **coeff)
 172 {
 173   long n,d,i,j;
 174   double *resi;
 175
 176   check_alloc(resi=(double*)malloc(sizeof(double)*dim));
 177   for (i=0;i<dim;i++)
 178     resi[i]=0.0;
 179
 180   for (n=poles-1;n<length-1;n++) {
 181     for (d=0;d<dim;d++) {
 182       diff[d][n+1]=series[d][n+1];
 183       for (i=0;i<dim;i++)
 184         for (j=0;j<poles;j++)
 185           diff[d][n+1] -= coeff[d][i*poles+j]*series[i][n-j];
 186       resi[d] += sqr(diff[d][n+1]);
 187     }
 188   }
 189   for (i=0;i<dim;i++)
 190     resi[i]=sqrt(resi[i]/((double)length-(double)poles));
 191
 192   return resi;
 193 }
 194
 195 void iterate_model(double **coeff,double *sigma,FILE *file)
 196 {
 197   long i,j,i1,i2,n,d;
 198   double **iterate,*swap;
 199
 200   check_alloc(iterate=(double**)malloc(sizeof(double*)*(poles+1)));
 201   for (i=0;i<=poles;i++)
 202     check_alloc(iterate[i]=(double*)malloc(sizeof(double)*dim));
 203   rnd_init(0x44325);
 204   for (i=0;i<1000;i++)
 205     gaussian(1.0);
 206   for (i=0;i<dim;i++)
 207     for (j=0;j<poles;j++)
 208       iterate[j][i]=gaussian(sigma[i]);
 209
 210   for (n=0;n<ilength;n++) {
 211     for (d=0;d<dim;d++) {
 212       iterate[poles][d]=gaussian(sigma[d]);
 213       for (i1=0;i1<dim;i1++)
 214         for (i2=0;i2<poles;i2++)
 215           iterate[poles][d] += coeff[d][i1*poles+i2]*iterate[poles-1-i2][i1];
 216     }
 217     if (file != NULL) {
 218       for (d=0;d<dim;d++)
 219         fprintf(file,"%e ",iterate[poles][d]);
 220       fprintf(file,"\n");
 221     }
 222     else {
 223       for (d=0;d<dim;d++)
 224         printf("%e ",iterate[poles][d]);
 225       printf("\n");
 226     }
 227
 228     swap=iterate[0];
 229     for (i=0;i<poles;i++)
 230       iterate[i]=iterate[i+1];
 231     iterate[poles]=swap;
 232   }
 233
 234   for (i=0;i<=poles;i++)
 235     free(iterate[i]);
 236   free(iterate);
 237 }
 238
 239 int main(int argc,char **argv)
 240 {
 241   char stdi=0;
 242   double *pm;
 243   long i,j;
 244   FILE *file;
 245   double **mat,**inverse,*vec,**coeff,**diff,avpm;
 246
 247   if (scan_help(argc,argv))
 248     show_options(argv[0]);
 249
 250   scan_options(argc,argv);
 251 #ifndef OMIT_WHAT_I_DO
 252   if (verbosity&VER_INPUT)
 253     what_i_do(argv[0],WID_STR);
 254 #endif
 255
 256   infile=search_datafile(argc,argv,NULL,verbosity);
 257   if (infile == NULL)
 258     stdi=1;
 259
 260   if (outfile == NULL) {
 261     if (!stdi) {
 262       check_alloc(outfile=(char*)calloc(strlen(infile)+4,(size_t)1));
 263       strcpy(outfile,infile);
 264       strcat(outfile,".ar");
 265     }
 266     else {
 267       check_alloc(outfile=(char*)calloc((size_t)9,(size_t)1));
 268       strcpy(outfile,"stdin.ar");
 269     }
 270   }
 271   if (!stdo)
 272     test_outfile(outfile);
 273
 274   if (column == NULL)
 275     series=(double**)get_multi_series(infile,&length,exclude,&dim,"",dimset,
 276                                       verbosity);
 277   else
 278     series=(double**)get_multi_series(infile,&length,exclude,&dim,column,
 279                                       dimset,verbosity);
 280
 281   check_alloc(my_average=(double*)malloc(sizeof(double)*dim));
 282   set_averages_to_zero();
 283
 284   if (poles >= length) {
 285     fprintf(stderr,"It makes no sense to have more poles than data! Exiting\n");
 286     exit(AR_MODEL_TOO_MANY_POLES);
 287   }
 288
 289
 290   check_alloc(vec=(double*)malloc(sizeof(double)*poles*dim));
 291   check_alloc(mat=(double**)malloc(sizeof(double*)*poles*dim));
 292   for (i=0;i<poles*dim;i++)
 293     check_alloc(mat[i]=(double*)malloc(sizeof(double)*poles*dim));
 294
 295   check_alloc(coeff=(double**)malloc(sizeof(double*)*dim));
 296   inverse=build_matrix(mat);
 297   for (i=0;i<dim;i++) {
 298     build_vector(vec,i);
 299     coeff[i]=multiply_matrix_vector(inverse,vec);
 300   }
 301
 302   check_alloc(diff=(double**)malloc(sizeof(double*)*dim));
 303   for (i=0;i<dim;i++)
 304     check_alloc(diff[i]=(double*)malloc(sizeof(double)*length));
 305
 306   pm=make_residuals(diff,coeff);
 307
 308   if (stdo) {
 309     avpm=pm[0]*pm[0];
 310     for (i=1;i<dim;i++)
 311       avpm += pm[i]*pm[i];
 312     avpm=sqrt(avpm/dim);
 313     printf("#average forcast error= %e\n",avpm);
 314     printf("#individual forecast errors: ");
 315     for (i=0;i<dim;i++)
 316       printf("%e ",pm[i]);
 317     printf("\n");
 318     for (i=0;i<dim*poles;i++) {
 319       printf("# ");
 320       for (j=0;j<dim;j++)
 321         printf("%e ",coeff[j][i]);
 322       printf("\n");
 323     }
 324     if (!run_model || (verbosity&VER_USR1)) {
 325       for (i=poles;i<length;i++) {
 326         if (run_model)
 327           printf("#");
 328         for (j=0;j<dim;j++)
 329           if (verbosity&VER_USR2)
 330             printf("%e %e ",series[j][i]+my_average[j],diff[j][i]);
 331           else
 332             printf("%e ",diff[j][i]);
 333         printf("\n");
 334       }
 335     }
 336     if (run_model && (ilength > 0))
 337       iterate_model(coeff,pm,NULL);
 338   }
 339   else {
 340     file=fopen(outfile,"w");
 341     if (verbosity&VER_INPUT)
 342       fprintf(stderr,"Opened %s for output\n",outfile);
 343     avpm=pm[0]*pm[0];
 344     for (i=1;i<dim;i++)
 345       avpm += pm[i]*pm[i];
 346     avpm=sqrt(avpm/dim);
 347     fprintf(file,"#average forcast error= %e\n",avpm);
 348     fprintf(file,"#individual forecast errors: ");
 349     for (i=0;i<dim;i++)
 350       fprintf(file,"%e ",pm[i]);
 351     fprintf(file,"\n");
 352     for (i=0;i<dim*poles;i++) {
 353       fprintf(file,"# ");
 354       for (j=0;j<dim;j++)
 355         fprintf(file,"%e ",coeff[j][i]);
 356       fprintf(file,"\n");
 357     }
 358     if (!run_model || (verbosity&VER_USR1)) {
 359       for (i=poles;i<length;i++) {
 360         if (run_model)
 361           fprintf(file,"#");
 362         for (j=0;j<dim;j++)
 363           if (verbosity&VER_USR2)
 364             fprintf(file,"%e %e ",series[j][i]+my_average[j],diff[j][i]);
 365           else
 366             fprintf(file,"%e ",diff[j][i]);
 367         fprintf(file,"\n");
 368       }
 369     }
 370     if (run_model && (ilength > 0))
 371       iterate_model(coeff,pm,file);
 372     fclose(file);
 373   }
 374
 375   if (outfile != NULL)
 376     free(outfile);
 377   if (infile != NULL)
 378     free(infile);
 379   free(vec);
 380   for (i=0;i<poles*dim;i++) {
 381     free(mat[i]);
 382     free(inverse[i]);
 383   }
 384   free(mat);
 385   free(inverse);
 386   for (i=0;i<dim;i++) {
 387     free(coeff[i]);
 388     free(diff[i]);
 389   }
 390   free(coeff);
 391   free(diff);
 392   free(pm);
 393
 394   return 0;
 395 }