binaries/src/disembl/Tisean_3.0.1/source_c/lfo-run.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 29, 2000 */
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include <limits.h>
  25 #include <math.h>
  26 #include "routines/tsa.h"
  27
  28 #define WID_STR "Makes a local linear fit for multivariate data\n\
  29 and iterates a trajectory"
  30
  31 #define NMAX 128
  32
  33 char onscreen=1,epsset=0,*outfile=NULL;
  34 char *infile=NULL;
  35 unsigned int nmax=(NMAX-1);
  36 unsigned int verbosity=0xff;
  37 long **box,*list,*found;
  38 double **series,**cast;
  39 double *interval,*min,epsilon;
  40
  41 unsigned int embed=2,dim=1,dim1,DELAY=1;
  42 char *column=NULL,dimset=0,do_zeroth=0;
  43 int MINN=30;
  44 unsigned long LENGTH=ULONG_MAX,FLENGTH=1000,exclude=0;
  45 double EPS0=1.e-3,EPSF=1.2;
  46
  47 double **mat,**imat,*vec,*localav,*foreav;
  48
  49 void show_options(char *progname)
  50 {
  51   what_i_do(progname,WID_STR);
  52   fprintf(stderr," Usage: %s [Options]\n",progname);
  53   fprintf(stderr," Options:\n");
  54   fprintf(stderr,"Everything not being a valid option will be interpreted"
  55           " as a possible"
  56           " datafile.\nIf no datafile is given stdin is read. Just - also"
  57           " means stdin\n");
  58   fprintf(stderr,"\t-l # of data to be used [default whole file]\n");
  59   fprintf(stderr,"\t-x # of lines to be ignored [default 0]\n");
  60   fprintf(stderr,"\t-c column [default 1,...,# of components]\n");
  61   fprintf(stderr,"\t-m #of components,embedding dimension [default 1,2]\n");
  62   fprintf(stderr,"\t-d delay for the embedding [default 1]\n");
  63   fprintf(stderr,"\t-L # of iterations [default 1000]\n");
  64   fprintf(stderr,"\t-k # of neighbors  [default 30]\n");
  65   fprintf(stderr,"\t-r size of initial neighborhood ["
  66           " default (data interval)/1000]\n");
  67   fprintf(stderr,"\t-f factor to increase size [default 1.2]\n");
  68   fprintf(stderr,"\t-0 perfom a zeroth order fit [default not set]\n");
  69   fprintf(stderr,"\t-o output file [default 'datafile'.cast;"
  70           " no -o means write to stdout]\n");
  71   fprintf(stderr,"\t-V verbosity level [default 1]\n\t\t"
  72           "0='only panic messages'\n\t\t"
  73           "1='+ input/output messages'\n");
  74   fprintf(stderr,"\t-h  show these options\n");
  75   exit(0);
  76 }
  77
  78 void scan_options(int n,char **in)
  79 {
  80   char *out;
  81
  82   if ((out=check_option(in,n,'l','u')) != NULL)
  83     sscanf(out,"%lu",&LENGTH);
  84   if ((out=check_option(in,n,'x','u')) != NULL)
  85     sscanf(out,"%lu",&exclude);
  86   if ((out=check_option(in,n,'c','s')) != NULL) {
  87     column=out;
  88     dimset=1;
  89   }
  90   if ((out=check_option(in,n,'m','2')) != NULL)
  91     sscanf(out,"%u,%u",&dim,&embed);
  92   if ((out=check_option(in,n,'d','u')) != NULL)
  93     sscanf(out,"%u",&DELAY);
  94   if ((out=check_option(in,n,'L','u')) != NULL)
  95     sscanf(out,"%lu",&FLENGTH);
  96   if ((out=check_option(in,n,'k','u')) != NULL)
  97     sscanf(out,"%u",&MINN);
  98   if ((out=check_option(in,n,'0','n')) != NULL)
  99     do_zeroth=1;
 100   if ((out=check_option(in,n,'V','u')) != NULL)
 101     sscanf(out,"%u",&verbosity);
 102   if ((out=check_option(in,n,'r','f')) != NULL) {
 103     epsset=1;
 104     sscanf(out,"%lf",&EPS0);
 105   }
 106   if ((out=check_option(in,n,'f','f')) != NULL)
 107     sscanf(out,"%lf",&EPSF);
 108   if ((out=check_option(in,n,'o','o')) != NULL) {
 109     onscreen=0;
 110     if (strlen(out) > 0)
 111       outfile=out;
 112   }
 113 }
 114
 115 void put_in_boxes(void)
 116 {
 117   int i,j,n;
 118   static int hdim;
 119   double epsinv;
 120
 121   hdim=(embed-1)*DELAY;
 122   epsinv=1.0/epsilon;
 123   for (i=0;i<NMAX;i++)
 124     for (j=0;j<NMAX;j++)
 125       box[i][j]= -1;
 126
 127   for (n=hdim;n<LENGTH-1;n++) {
 128     i=(int)(series[0][n]*epsinv)&nmax;
 129     j=(int)(series[dim1][n-hdim]*epsinv)&nmax;
 130     list[n]=box[i][j];
 131     box[i][j]=n;
 132   }
 133 }
 134
 135 unsigned int hfind_neighbors(void)
 136 {
 137   char toolarge;
 138   int i,j,i1,i2,j1,k,l,element;
 139   static int hdim;
 140   unsigned nfound=0;
 141   double max,dx,epsinv;
 142
 143   hdim=(embed-1)*DELAY;
 144   epsinv=1.0/epsilon;
 145   i=(int)(cast[hdim][0]*epsinv)&nmax;
 146   j=(int)(cast[0][dim1]*epsinv)&nmax;
 147
 148   for (i1=i-1;i1<=i+1;i1++) {
 149     i2=i1&nmax;
 150     for (j1=j-1;j1<=j+1;j1++) {
 151       element=box[i2][j1&nmax];
 152       while (element != -1) {
 153         max=0.0;
 154         toolarge=0;
 155         for (l=0;l<dim;l++) {
 156           for (k=0;k<=hdim;k += DELAY) {
 157             dx=fabs(series[l][element-k]-cast[hdim-k][l]);
 158             max=(dx>max) ? dx : max;
 159             if (max > epsilon) {
 160               toolarge=1;
 161               break;
 162             }
 163           }
 164           if (toolarge)
 165             break;
 166         }
 167         if (max <= epsilon)
 168           found[nfound++]=element;
 169         element=list[element];
 170       }
 171     }
 172   }
 173   return nfound;
 174 }
 175
 176 void multiply_matrix(double **mat,double *vec)
 177 {
 178   double *hvec;
 179   long i,j;
 180
 181   check_alloc(hvec=(double*)malloc(sizeof(double)*dim*embed));
 182   for (i=0;i<dim*embed;i++) {
 183     hvec[i]=0.0;
 184     for (j=0;j<dim*embed;j++)
 185       hvec[i] += mat[i][j]*vec[j];
 186   }
 187   for (i=0;i<dim*embed;i++)
 188     vec[i]=hvec[i];
 189   free(hvec);
 190 }
 191
 192 void make_fit(int number,double *newcast)
 193 {
 194   double *sj,*si,lavi,lavj,fav;
 195   long i,i1,j,j1,hi,hj,hi1,hj1,n,which;
 196   static int hdim;
 197
 198   hdim=(embed-1)*DELAY;
 199
 200   for (i=0;i<dim*embed;i++)
 201     localav[i]=0.0;
 202   for (i=0;i<dim;i++)
 203     foreav[i]=0.0;
 204
 205   for (n=0;n<number;n++) {
 206     which=found[n];
 207     for (j=0;j<dim;j++) {
 208       sj=series[j];
 209       foreav[j] += sj[which+1];
 210       for (j1=0;j1<embed;j1++) {
 211         hj=j*embed+j1;
 212         localav[hj] += sj[which-j1*DELAY];
 213       }
 214     }
 215   }
 216
 217   for (i=0;i<dim*embed;i++)
 218     localav[i] /= number;
 219   for (i=0;i<dim;i++)
 220     foreav[i] /= number;
 221
 222   for (i=0;i<dim;i++) {
 223     si=series[i];
 224     for (i1=0;i1<embed;i1++) {
 225       hi=i*embed+i1;
 226       lavi=localav[hi];
 227       hi1=i1*DELAY;
 228       for (j=0;j<dim;j++) {
 229         sj=series[j];
 230         for (j1=0;j1<embed;j1++) {
 231           hj=j*embed+j1;
 232           lavj=localav[hj];
 233           hj1=j1*DELAY;
 234           mat[hi][hj]=0.0;
 235           if (hj >= hi) {
 236             for (n=0;n<number;n++) {
 237               which=found[n];
 238               mat[hi][hj] += (si[which-hi1]-lavi)*(sj[which-hj1]-lavj);
 239             }
 240           }
 241         }
 242       }
 243     }
 244   }
 245
 246   for (i=0;i<dim*embed;i++)
 247     for (j=i;j<dim*embed;j++) {
 248       mat[i][j] /= number;
 249       mat[j][i]=mat[i][j];
 250     }
 251
 252   imat=invert_matrix(mat,dim*embed);
 253
 254   for (i=0;i<dim;i++) {
 255     si=series[i];
 256     fav=foreav[i];
 257     for (j=0;j<dim;j++) {
 258       sj=series[j];
 259       for (j1=0;j1<embed;j1++) {
 260         hj=j*embed+j1;
 261         lavj=localav[hj];
 262         hj1=j1*DELAY;
 263         vec[hj]=0.0;
 264         for (n=0;n<number;n++) {
 265           which=found[n];
 266           vec[hj] += (si[which+1]-fav)*(sj[which-hj1]-lavj);
 267         }
 268         vec[hj] /= number;
 269       }
 270     }
 271
 272     multiply_matrix(imat,vec);
 273
 274     newcast[i]=foreav[i];
 275     for (j=0;j<dim;j++) {
 276       for (j1=0;j1<embed;j1++) {
 277         hj=j*embed+j1;
 278         newcast[i] += vec[hj]*(cast[hdim-j1*DELAY][j]-localav[hj]);
 279       }
 280     }
 281   }
 282
 283   for (i=0;i<dim*embed;i++)
 284     free(imat[i]);
 285   free(imat);
 286 }
 287
 288 void make_zeroth(int number,double *newcast)
 289 {
 290   unsigned long i,d;
 291   double *sj;
 292
 293   for (d=0;d<dim;d++) {
 294     newcast[d]=0.0;
 295     sj=series[d]+1;
 296     for (i=0;i<number;i++)
 297       newcast[d] += sj[found[i]];
 298     newcast[d] /= number;
 299   }
 300 }
 301
 302 int main(int argc,char **argv)
 303 {
 304   char stdi=0,done;
 305   long i,j,hdim,actfound;
 306   double maxinterval,*swap,*newcast;
 307   FILE *file=NULL;
 308
 309   if (scan_help(argc,argv))
 310     show_options(argv[0]);
 311
 312   scan_options(argc,argv);
 313 #ifndef OMIT_WHAT_I_DO
 314   if (verbosity&VER_INPUT)
 315     what_i_do(argv[0],WID_STR);
 316 #endif
 317
 318   infile=search_datafile(argc,argv,NULL,verbosity);
 319   if (infile == NULL)
 320     stdi=1;
 321
 322   if (outfile == NULL) {
 323     if (!stdi) {
 324       check_alloc(outfile=(char*)calloc(strlen(infile)+6,(size_t)1));
 325       strcpy(outfile,infile);
 326       strcat(outfile,".cast");
 327     }
 328     else {
 329       check_alloc(outfile=(char*)calloc((size_t)11,(size_t)1));
 330       strcpy(outfile,"stdin.cast");
 331     }
 332   }
 333   if (!onscreen)
 334     test_outfile(outfile);
 335
 336   hdim=(embed-1)*DELAY+1;
 337   if (column == NULL)
 338     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,"",dimset,
 339                                       verbosity);
 340   else
 341     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dim,column,
 342                                       dimset,verbosity);
 343   check_alloc(min=(double*)malloc(sizeof(double)*dim));
 344   check_alloc(interval=(double*)malloc(sizeof(double)*dim));
 345   dim1=dim-1;
 346   maxinterval=0.0;
 347   for (i=0;i<dim;i++) {
 348     rescale_data(series[i],LENGTH,&min[i],&interval[i]);
 349     if (interval[i] > maxinterval)
 350       maxinterval=interval[i];
 351   }
 352
 353   check_alloc(cast=(double**)malloc(sizeof(double*)*hdim));
 354   for (i=0;i<hdim;i++)
 355     check_alloc(cast[i]=(double*)malloc(sizeof(double)*dim));
 356   check_alloc(newcast=(double*)malloc(sizeof(double)*dim));
 357
 358   check_alloc(list=(long*)malloc(sizeof(long)*LENGTH));
 359   check_alloc(found=(long*)malloc(sizeof(long)*LENGTH));
 360   check_alloc(box=(long**)malloc(sizeof(long*)*NMAX));
 361   for (i=0;i<NMAX;i++)
 362     check_alloc(box[i]=(long*)malloc(sizeof(long)*NMAX));
 363
 364   check_alloc(localav=(double*)malloc(sizeof(double)*dim*embed));
 365   check_alloc(foreav=(double*)malloc(sizeof(double)*dim));
 366   check_alloc(vec=(double*)malloc(sizeof(double)*dim*embed));
 367   check_alloc(mat=(double**)malloc(sizeof(double*)*dim*embed));
 368   for (i=0;i<dim*embed;i++)
 369     check_alloc(mat[i]=(double*)malloc(sizeof(double)*dim*embed));
 370
 371   if (epsset)
 372     EPS0 /= maxinterval;
 373
 374   for (j=0;j<dim;j++)
 375     for (i=0;i<hdim;i++)
 376       cast[i][j]=series[j][LENGTH-hdim+i];
 377
 378   if (!onscreen) {
 379     file=fopen(outfile,"w");
 380     if (verbosity&VER_INPUT)
 381       fprintf(stderr,"Opened %s for writing\n",outfile);
 382   }
 383   else {
 384     if (verbosity&VER_INPUT)
 385       fprintf(stderr,"Writing to stdout\n");
 386   }
 387
 388   for (i=0;i<FLENGTH;i++) {
 389     done=0;
 390     epsilon=EPS0/EPSF;
 391     while (!done) {
 392       epsilon*=EPSF;
 393       put_in_boxes();
 394       actfound=hfind_neighbors();
 395       if (actfound >= MINN) {
 396         if (!do_zeroth)
 397           make_fit(actfound,newcast);
 398         else
 399           make_zeroth(actfound,newcast);
 400         if (onscreen) {
 401           for (j=0;j<dim-1;j++)
 402             printf("%e ",newcast[j]*interval[j]+min[j]);
 403           printf("%e\n",newcast[dim-1]*interval[dim-1]+min[dim-1]);
 404           fflush(stdout);
 405         }
 406         else {
 407           for (j=0;j<dim-1;j++)
 408             fprintf(file,"%e ",newcast[j]*interval[j]+min[j]);
 409           fprintf(file,"%e\n",newcast[dim-1]*interval[dim-1]+min[dim-1]);
 410           fflush(file);
 411         }
 412         done=1;
 413         for (j=0;j<dim;j++) {
 414           if ((newcast[j] > 2.0) || (newcast[j] < -1.0)) {
 415             fprintf(stderr,"Forecast failed. Escaping data region!\n");
 416             exit(NSTEP__ESCAPE_REGION);
 417           }
 418         }
 419
 420         swap=cast[0];
 421         for (j=0;j<hdim-1;j++)
 422           cast[j]=cast[j+1];
 423         cast[hdim-1]=swap;
 424         for (j=0;j<dim;j++)
 425           cast[hdim-1][j]=newcast[j];
 426       }
 427     }
 428   }
 429   if (!onscreen)
 430     fclose(file);
 431
 432   if (outfile != NULL)
 433     free(outfile);
 434   for (i=0;i<embed*dim;i++)
 435     free(mat[i]);
 436   free(mat);
 437   for (i=0;i<hdim;i++)
 438     free(cast[i]);
 439   free(cast);
 440   free(newcast);
 441   free(found);
 442   free(list);
 443   for (i=0;i<NMAX;i++)
 444     free(box[i]);
 445   free(box);
 446   free(vec);
 447   free(localav);
 448   free(foreav);
 449   free(min);
 450   free(interval);
 451   for (i=0;i<dim;i++)
 452     free(series[i]);
 453   free(series);
 454
 455   return 0;
 456 }