website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_c/ghkss.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: Jun 10, 2006 */
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include <math.h>
  25 #include <limits.h>
  26 #include "routines/tsa.h"
  27
  28 #define WID_STR "Multivariate noise reduction using the GHKSS algorithm"
  29
  30
  31 #define BOX (unsigned int)1024
  32
  33 unsigned long length=ULONG_MAX,exclude=0;
  34 unsigned int dim,qdim=2,delay=1,minn=50,iterations=1,comp=1,embed=5;
  35 unsigned int verbosity=0xff;
  36 double mineps,epsfac;
  37 char *column=NULL;
  38 char eps_set=0,euclidean=0,dimset=0,resize_eps;
  39 char *outfile=NULL,stdo=1;
  40 char *infile=NULL;
  41
  42 double *d_min,*d_max,d_max_max;
  43 double **series,**delta,**corr;
  44 double *metric,trace;
  45 long **box,*list;
  46 unsigned long *flist;
  47 int emb_offset;
  48 unsigned int ibox=BOX-1;
  49 unsigned int *index_comp,*index_embed;
  50
  51 /*these are global to save time*/
  52 int *sorted;
  53 double *av,**mat,*matarray,*eig;
  54
  55 void show_options(char *progname)
  56 {
  57   what_i_do(progname,WID_STR);
  58   fprintf(stderr,"Usage: %s [options]\n",progname);
  59   fprintf(stderr,"Options:\n");
  60   fprintf(stderr,"Everything not being a valid option will be interpreted"
  61           " as a possible"
  62           " datafile.\nIf no datafile is given stdin is read. Just - also"
  63           " means stdin\n");
  64   fprintf(stderr,"\t-l # of data to use [Default: whole file]\n");
  65   fprintf(stderr,"\t-x # of lines to be ignored [Default: 0]\n");
  66   fprintf(stderr,"\t-c column to read [Default: 1,..,# of components]\n");
  67   fprintf(stderr,"\t-m # of components,embedding dimension [Default: 1,5]\n");
  68   fprintf(stderr,"\t-d delay [Default: 1]\n");
  69   fprintf(stderr,"\t-q dimension to project to [Default: 2]\n");
  70   fprintf(stderr,"\t-k minimal number of neighbours [Default: 50]\n");
  71   fprintf(stderr,"\t-r minimal neighbourhood size \n\t\t"
  72           "[Default: (interval of data)/1000]\n");
  73   fprintf(stderr,"\t-i # of iterations [Default: 1]\n");
  74   fprintf(stderr,"\t-2 use euklidean metric [Default: non euklidean]\n");
  75   fprintf(stderr,"\t-o name of output file \n\t\t"
  76           "[Default: 'datafile'.opt.n, where n is the iteration.\n\t\t"
  77           " If no -o is given, the last iteration is also"
  78           " written to stdout]\n");
  79   fprintf(stderr,"\t-V verbosity level [Default: 7]\n\t\t"
  80           "0='only panic messages'\n\t\t"
  81           "1='+ input/output messages'\n\t\t"
  82           "2='+ average correction and trend'\n\t\t"
  83           "4='+ how many points for which epsilon'\n");
  84   fprintf(stderr,"\t-h show these options\n");
  85   exit(0);
  86 }
  87
  88 void scan_options(int n,char **in)
  89 {
  90   char *out;
  91
  92   if ((out=check_option(in,n,'l','u')) != NULL)
  93     sscanf(out,"%lu",&length);
  94   if ((out=check_option(in,n,'x','u')) != NULL)
  95     sscanf(out,"%lu",&exclude);
  96   if ((out=check_option(in,n,'c','s')) != NULL) {
  97     column=out;
  98     dimset=1;
  99   }
 100   if ((out=check_option(in,n,'m','2')) != NULL)
 101     sscanf(out,"%u,%u",&comp,&embed);
 102   if ((out=check_option(in,n,'d','u')) != NULL)
 103     sscanf(out,"%u",&delay);
 104   if ((out=check_option(in,n,'q','u')) != NULL)
 105     sscanf(out,"%u",&qdim);
 106   if ((out=check_option(in,n,'k','u')) != NULL)
 107     sscanf(out,"%u",&minn);
 108   if ((out=check_option(in,n,'r','f')) != NULL) {
 109     eps_set=1;
 110     sscanf(out,"%lf",&mineps);
 111   }
 112   if ((out=check_option(in,n,'i','u')) != NULL)
 113     sscanf(out,"%u",&iterations);
 114   if ((out=check_option(in,n,'V','u')) != NULL)
 115     sscanf(out,"%u",&verbosity);
 116   if ((out=check_option(in,n,'2','n')) != NULL)
 117     euclidean=1;
 118   if ((out=check_option(in,n,'o','o')) != NULL) {
 119     stdo=0;
 120     if (strlen(out) > 0)
 121       outfile=out;
 122   }
 123 }
 124
 125 void sort(double *x,int *n)
 126 {
 127   long i,j,iswap;
 128   double dswap;
 129
 130   for (i=0;i<dim;i++)
 131     n[i]=i;
 132
 133   for (i=0;i<dim-1;i++)
 134     for (j=i+1;j<dim;j++)
 135       if (x[j] > x[i]) {
 136         dswap=x[i];
 137         x[i]=x[j];
 138         x[j]=dswap;
 139         iswap=n[i];
 140         n[i]=n[j];
 141         n[j]=iswap;
 142       }
 143 }
 144
 145 void mmb(double eps)
 146 {
 147   long i,x,y;
 148   double ieps=1.0/eps;
 149
 150   for (x=0;x<BOX;x++)
 151     for (y=0;y<BOX;y++)
 152       box[x][y] = -1;
 153
 154   for (i=emb_offset;i<length;i++) {
 155     x=(int)(series[0][i]*ieps)&ibox;
 156     y=(int)(series[comp-1][i-emb_offset]*ieps)&ibox;
 157     list[i]=box[x][y];
 158     box[x][y]=i;
 159   }
 160 }
 161
 162 unsigned long fmn(long which,double eps)
 163 {
 164   unsigned long nf=0;
 165   long i,i1,i2,j,j1,k,k1,li;
 166   long element;
 167   double dx=0.0;
 168
 169   i=(int)(series[0][which]/eps)&ibox;
 170   j=(int)(series[comp-1][which-emb_offset]/eps)&ibox;
 171
 172   for (i1=i-1;i1<=i+1;i1++) {
 173     i2=i1&ibox;
 174     for (j1=j-1;j1<=j+1;j1++) {
 175       element=box[i2][j1&ibox];
 176       while (element != -1) {
 177         for (k=0;k<embed;k++) {
 178           k1= -k*(int)delay;
 179           for (li=0;li<comp;li++) {
 180             dx=fabs(series[li][which+k1]-series[li][element+k1]);
 181             if (dx > eps)
 182               break;
 183           }
 184           if (dx > eps)
 185             break;
 186         }
 187         if (dx <= eps)
 188           flist[nf++]=element;
 189         element=list[element];
 190       }
 191     }
 192   }
 193   return nf;
 194 }
 195
 196 void make_correction(unsigned long n,unsigned long nf)
 197 {
 198   long i,i1,i2,j,j1,j2,k,k1,k2,hs;
 199   double help;
 200
 201   for (i=0;i<dim;i++) {
 202     i1=index_comp[i];
 203     i2=index_embed[i];
 204     help=0.0;
 205     for (j=0;j<nf;j++)
 206       help += series[i1][flist[j]-i2];
 207     av[i]=help/nf;
 208   }
 209
 210   for (i=0;i<dim;i++) {
 211     i1=index_comp[i];
 212     i2=index_embed[i];
 213     for (j=i;j<dim;j++) {
 214       help=0.0;
 215       j1=index_comp[j];
 216       j2=index_embed[j];
 217       for (k=0;k<nf;k++) {
 218         hs=flist[k];
 219         help += series[i1][hs-i2]*series[j1][hs-j2];
 220       }
 221       mat[i][j]=(help/nf-av[i]*av[j])*metric[i]*metric[j];
 222       mat[j][i]=mat[i][j];
 223     }
 224   }
 225
 226   eigen(mat,(unsigned long)dim,eig);
 227   sort(eig,sorted);
 228
 229   for (i=0;i<dim;i++) {
 230     help=0.0;
 231     for (j=qdim;j<dim;j++) {
 232       hs=sorted[j];
 233       for (k=0;k<dim;k++) {
 234         k1=index_comp[k];
 235         k2=index_embed[k];
 236         help += (series[k1][n-k2]-av[k])*mat[k][hs]*mat[i][hs]*metric[k];
 237       }
 238     }
 239     corr[n][i]=help/metric[i];
 240   }
 241 }
 242
 243 void handle_trend(unsigned long n,unsigned long nf)
 244 {
 245   long i,i1,i2,j;
 246   double help;
 247
 248   for (i=0;i<dim;i++) {
 249     help=0.0;
 250     for (j=0;j<nf;j++)
 251       help += corr[flist[j]][i];
 252     av[i]=help/nf;
 253   }
 254
 255   for (i=0;i<dim;i++) {
 256     i1=index_comp[i];
 257     i2=index_embed[i];
 258     delta[i1][n-i2] += (corr[n][i]-av[i])/(trace*metric[i]);
 259   }
 260 }
 261
 262 void set_correction(void)
 263 {
 264   long i,j;
 265   double *hav,*hsigma,help;
 266
 267   check_alloc(hav=(double*)malloc(sizeof(double)*comp));
 268   check_alloc(hsigma=(double*)malloc(sizeof(double)*comp));
 269   for (j=0;j<comp;j++)
 270     hav[j]=hsigma[j]=0.0;
 271
 272   for (i=0;i<length;i++)
 273     for (j=0;j<comp;j++) {
 274       hav[j] += (help=delta[j][i]);
 275       hsigma[j] += help*help;
 276     }
 277
 278   for (j=0;j<comp;j++) {
 279     hav[j] /= length;
 280     hsigma[j]=sqrt(fabs(hsigma[j]/length-hav[j]*hav[j]));
 281   }
 282   if (verbosity&(VER_USR1|VER_USR2)) {
 283     for (i=0;i<comp;i++) {
 284       fprintf(stderr,"Average shift of component %ld = %e\n",i+1,
 285               hav[i]*d_max[i]);
 286       fprintf(stderr,"Average rms correction of comp. %ld = %e\n\n",
 287               i+1,hsigma[i]*d_max[i]);
 288     }
 289   }
 290   for (i=0;i<length;i++)
 291     for (j=0;j<comp;j++)
 292       series[j][i] -= delta[j][i];
 293
 294   if (resize_eps) {
 295     mineps /= epsfac;
 296     if (verbosity&VER_USR2)
 297       fprintf(stderr,"Reset minimal neighbourhood size to %e\n",
 298               mineps*d_max_max);
 299   }
 300
 301   resize_eps=0;
 302   free(hav);
 303   free(hsigma);
 304 }
 305
 306 int main(int argc,char **argv)
 307 {
 308   char stdi=0;
 309   int iter,epscount,*ok;
 310   long i,j;
 311   char all_done;
 312   char *ofname;
 313   unsigned long nfound,n,allfound;
 314   double epsilon;
 315   double **hser;
 316   FILE *file;
 317
 318   if (scan_help(argc,argv))
 319     show_options(argv[0]);
 320
 321   scan_options(argc,argv);
 322 #ifndef OMIT_WHAT_I_DO
 323   if (verbosity&VER_INPUT)
 324     what_i_do(argv[0],WID_STR);
 325 #endif
 326
 327   dim=comp*embed;
 328   emb_offset=(embed-1)*delay;
 329
 330   infile=search_datafile(argc,argv,NULL,verbosity);
 331   if (infile == NULL)
 332     stdi=1;
 333
 334   if (outfile == NULL) {
 335     if (!stdi) {
 336       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
 337       check_alloc(ofname=(char*)calloc(strlen(infile)+9,(size_t)1));
 338       sprintf(outfile,"%s.opt",infile);
 339     }
 340     else {
 341       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
 342       check_alloc(ofname=(char*)calloc((size_t)14,(size_t)1));
 343       sprintf(outfile,"stdin.opt");
 344     }
 345   }
 346   else
 347     check_alloc(ofname=(char*)calloc(strlen(outfile)+10,(size_t)1));
 348
 349   if (column == NULL)
 350     series=(double**)get_multi_series(infile,&length,exclude,&comp,"",
 351                                      dimset,verbosity);
 352   else
 353     series=(double**)get_multi_series(infile,&length,exclude,&comp,column,
 354                                       dimset,verbosity);
 355
 356   if (length < minn) {
 357     fprintf(stderr,"With %lu data you will never find %u neighbors."
 358             " Exiting!\n",length,minn);
 359     exit(GHKSS__TOO_MANY_NEIGHBORS);
 360   }
 361
 362   check_alloc(d_min=(double*)malloc(sizeof(double)*comp));
 363   check_alloc(d_max=(double*)malloc(sizeof(double)*comp));
 364   d_max_max=0.0;
 365   for (i=0;i<comp;i++) {
 366     rescale_data(series[i],length,&d_min[i],&d_max[i]);
 367     if (d_max[i] > d_max_max)
 368       d_max_max=d_max[i];
 369   }
 370
 371   if (!eps_set)
 372     mineps=1./1000.;
 373   else
 374     mineps /= d_max_max;
 375   epsfac=sqrt(2.0);
 376
 377   check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
 378   for (i=0;i<BOX;i++)
 379     check_alloc(box[i]=(long*)malloc(sizeof(long)*BOX));
 380
 381   check_alloc(list=(long*)malloc(sizeof(long)*length));
 382   check_alloc(flist=(unsigned long*)malloc(sizeof(long)*length));
 383
 384   check_alloc(metric=(double*)malloc(sizeof(double)*dim));
 385   trace=0.0;
 386   if (euclidean) {
 387     for (i=0;i<dim;i++) {
 388       metric[i]=1.0;
 389       trace += 1./metric[i];
 390     }
 391   }
 392   else {
 393     for (i=0;i<dim;i++) {
 394       if ((i >= comp) && (i < ((long)dim-(long)comp)))
 395         metric[i]=1.0;
 396       else
 397         metric[i]=1.0e3;
 398       trace += 1./metric[i];
 399     }
 400   }
 401
 402   check_alloc(corr=(double**)malloc(sizeof(double*)*length));
 403   for (i=0;i<length;i++)
 404     check_alloc(corr[i]=(double*)malloc(sizeof(double)*dim));
 405   check_alloc(ok=(int*)malloc(sizeof(int)*length));
 406   check_alloc(delta=(double**)malloc(sizeof(double*)*comp));
 407   for (i=0;i<comp;i++)
 408     check_alloc(delta[i]=(double*)malloc(sizeof(double)*length));
 409   check_alloc(index_comp=(unsigned int*)malloc(sizeof(int)*dim));
 410   check_alloc(index_embed=(unsigned int*)malloc(sizeof(int)*dim));
 411   check_alloc(av=(double*)malloc(sizeof(double)*dim));
 412   check_alloc(sorted=(int*)malloc(sizeof(int)*dim));
 413   check_alloc(eig=(double*)malloc(sizeof(double)*dim));
 414   check_alloc(matarray=(double*)malloc(sizeof(double)*dim*dim));
 415   check_alloc(mat=(double**)malloc(sizeof(double*)*dim));
 416   for (i=0;i<dim;i++)
 417     mat[i]=(double*)(matarray+dim*i);
 418   check_alloc(hser=(double**)malloc(sizeof(double*)*comp));
 419
 420   for (i=0;i<dim;i++) {
 421     index_comp[i]=i%comp;
 422     index_embed[i]=(i/comp)*delay;
 423   }
 424
 425   resize_eps=0;
 426   for (iter=1;iter<=iterations;iter++) {
 427     for (i=0;i<length;i++) {
 428       ok[i]=0;
 429       for (j=0;j<dim;j++)
 430         corr[i][j]=0.0;
 431       for (j=0;j<comp;j++)
 432         delta[j][i]=0.0;
 433     }
 434     epsilon=mineps;
 435     all_done=0;
 436     epscount=1;
 437     allfound=0;
 438     if (verbosity&(VER_USR1|VER_USR2))
 439       fprintf(stderr,"Starting iteration %d\n",iter);
 440     while(!all_done) {
 441       mmb(epsilon);
 442       all_done=1;
 443       for (n=emb_offset;n<length;n++)
 444         if (!ok[n]) {
 445           nfound=fmn(n,epsilon);
 446           if (nfound >= minn) {
 447             make_correction(n,nfound);
 448             ok[n]=epscount;
 449             if (epscount == 1)
 450               resize_eps=1;
 451             allfound++;
 452           }
 453           else
 454             all_done=0;
 455         }
 456       if (verbosity&VER_USR2)
 457         fprintf(stderr,"Corrected %ld points with epsilon= %e\n",allfound,
 458                 epsilon*d_max_max);
 459       epsilon *= epsfac;
 460       epscount++;
 461     }
 462     if (verbosity&VER_USR2)
 463       fprintf(stderr,"Start evaluating the trend\n");
 464
 465     epsilon=mineps;
 466     allfound=0;
 467     for (i=1;i<epscount;i++) {
 468       mmb(epsilon);
 469       for (n=emb_offset;n<length;n++)
 470         if (ok[n] == i) {
 471           nfound=fmn(n,epsilon);
 472           handle_trend(n,nfound);
 473           allfound++;
 474         }
 475       if (verbosity&VER_USR2)
 476         fprintf(stderr,"Trend subtracted for %ld points with epsilon= %e\n",
 477                 allfound,epsilon*d_max_max);
 478       epsilon *= epsfac;
 479     }
 480     set_correction();
 481
 482     sprintf(ofname,"%s.%d",outfile,iter);
 483     test_outfile(ofname);
 484
 485     file=fopen(ofname,"w");
 486     if (verbosity&VER_INPUT)
 487       fprintf(stderr,"Opened %s for writing\n\n",ofname);
 488     for (i=0;i<length;i++) {
 489       for (j=0;j<comp;j++) {
 490         fprintf(file,"%e ",series[j][i]*d_max[j]+d_min[j]);
 491       }
 492       fprintf(file,"\n");
 493       if (stdo && (iter == iterations)) {
 494         for (j=0;j<comp;j++)
 495           fprintf(stdout,"%e ",series[j][i]*d_max[j]+d_min[j]);
 496         fprintf(stdout,"\n");
 497       }
 498     }
 499     fclose(file);
 500   }
 501
 502   return 0;
 503 }