binaries/src/disembl/Tisean_3.0.1/source_c/nrlazy.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: Nov 30, 2000 */
  21 /*Changes:
  22   12/11/05: Going multivariate
  23 */
  24
  25 #include <stdio.h>
  26 #include <stdlib.h>
  27 #include <string.h>
  28 #include <math.h>
  29 #include <limits.h>
  30 #include "routines/tsa.h"
  31
  32 #define WID_STR "Performs simple noise reduction."
  33
  34 #define BOX (unsigned int)512
  35
  36 unsigned long length=ULONG_MAX,exclude=0;
  37 unsigned int comp=1,embed=5,delay=1,iterations=1,alldim;
  38 unsigned int verbosity=0x3;
  39 char *column=NULL;
  40 double eps=1.0e-3,epsvar;
  41
  42 char *outfile=NULL,epsset=0,stdo=1,epsvarset=0;
  43 char *infile=NULL;
  44 double **series,**corr,*interval,*min,*hcor;
  45 long **box,*list,**nf;
  46 unsigned int **indexes;
  47 char dimset=0;
  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 use [default: whole file]\n");
  59   fprintf(stderr,"\t-x # of lines to be ignored [default: 0]\n");
  60   fprintf(stderr,"\t-c column to read [default: 1]\n");
  61   fprintf(stderr,"\t-m no. of comp.,embedding dim. [default: %u,%u]\n",
  62           comp,embed);
  63   fprintf(stderr,"\t-d delay [default: 1]\n");
  64   fprintf(stderr,"\t-i iterations [default: 1]\n");
  65   fprintf(stderr,"\t-r neighborhoud size [default: (interval of data)/1000]\n");
  66   fprintf(stderr,"\t-v neighborhoud size (in units of the std. dev. of the "
  67           "data \n\t\t(overwrites -r) [default: not set]\n");
  68   fprintf(stderr,"\t-o output file name [Default: 'datafile'.laz.n,"
  69           "\n\t\twhere n is the number of the last iteration,"
  70           "\n\t\twithout -o the last iteration is written to stdout.]\n");
  71   fprintf(stderr,"\t-V verbosity level [Default: 3]\n\t\t"
  72           "0='only panic messages'\n\t\t"
  73           "1='+ input/output messages'\n\t\t"
  74           "2='+ write output of all iterations to files'\n\t\t"
  75           "4='+ write the number of neighbors found for each point\n");
  76   fprintf(stderr,"\t-h show these options\n");
  77   exit(0);
  78 }
  79
  80 void scan_options(int n,char **in)
  81 {
  82   char *out;
  83
  84   if ((out=check_option(in,n,'l','u')) != NULL)
  85     sscanf(out,"%lu",&length);
  86   if ((out=check_option(in,n,'x','u')) != NULL)
  87     sscanf(out,"%lu",&exclude);
  88   if ((out=check_option(in,n,'c','c')) != NULL) {
  89     column=out;
  90     dimset=1;
  91   }
  92   if ((out=check_option(in,n,'m','2')) != NULL)
  93     sscanf(out,"%u,%u",&comp,&embed);
  94   if ((out=check_option(in,n,'d','u')) != NULL)
  95     sscanf(out,"%u",&delay);
  96   if ((out=check_option(in,n,'i','u')) != NULL)
  97     sscanf(out,"%u",&iterations);
  98   if ((out=check_option(in,n,'r','f')) != NULL) {
  99     epsset=1;
 100     sscanf(out,"%lf",&eps);
 101   }
 102   if ((out=check_option(in,n,'V','u')) != NULL)
 103     sscanf(out,"%u",&verbosity);
 104   if ((out=check_option(in,n,'v','f')) != NULL) {
 105     epsvarset=1;
 106     sscanf(out,"%lf",&epsvar);
 107   }
 108   if ((out=check_option(in,n,'o','o')) != NULL) {
 109     stdo=0;
 110     if (strlen(out) > 0)
 111       outfile=out;
 112   }
 113 }
 114
 115 unsigned int correct(unsigned long n)
 116 {
 117   int i,i1,i2,j,j1,k;
 118   int ibox=BOX-1;
 119   unsigned int hdel,hcomp;
 120   double epsinv,dx;
 121   long element,nfound=0;
 122
 123   epsinv=1./eps;
 124
 125   for (i=0;i<alldim;i++)
 126     hcor[i]=0.0;
 127
 128   i=(int)(series[0][n]*epsinv)&ibox;
 129   j=(int)(series[comp-1][n-(embed-1)*delay]*epsinv)&ibox;
 130
 131   for (i1=i-1;i1<=i+1;i1++) {
 132     i2=i1&ibox;
 133     for (j1=j-1;j1<=j+1;j1++) {
 134       element=box[i2][j1&ibox];
 135       while (element != -1) {
 136         for (k=0;k<alldim;k++) {
 137           hcomp=indexes[0][k];
 138           hdel=indexes[1][k];
 139           dx=fabs(series[hcomp][n-hdel]-series[hcomp][element-hdel]);
 140           if (dx > eps)
 141             break;
 142         }
 143         if (k == alldim) {
 144           nfound++;
 145           for (k=0;k<alldim;k++) {
 146             hcomp=indexes[0][k];
 147             hdel=indexes[1][k];
 148             hcor[k] += series[hcomp][element-hdel];
 149           }
 150         }
 151         element=list[element];
 152       }
 153     }
 154   }
 155   for (k=0;k<alldim;k++) {
 156     hcomp=indexes[0][k];
 157     hdel=indexes[1][k];
 158     corr[hcomp][n-hdel] += hcor[k]/nfound;
 159     nf[hcomp][n-hdel]++;
 160   }
 161
 162   return nfound;
 163 }
 164
 165 int main(int argc,char **argv)
 166 {
 167   char *ofname;
 168   char stdi=0;
 169   int iter;
 170   unsigned int *nmf;
 171   unsigned long n,i;
 172   double dav,dvar,maxinterval,maxdvar;
 173   FILE *file=NULL;
 174
 175   if (scan_help(argc,argv))
 176     show_options(argv[0]);
 177
 178   scan_options(argc,argv);
 179 #ifndef OMIT_WHAT_I_DO
 180   if (verbosity&VER_INPUT)
 181     what_i_do(argv[0],WID_STR);
 182 #endif
 183
 184   infile=search_datafile(argc,argv,NULL,verbosity);
 185   if (infile == NULL)
 186     stdi=1;
 187
 188   if (outfile == NULL) {
 189     if (!stdi) {
 190       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
 191       check_alloc(ofname=(char*)calloc(strlen(infile)+9,(size_t)1));
 192       sprintf(outfile,"%s.laz",infile);
 193     }
 194     else {
 195       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
 196       check_alloc(ofname=(char*)calloc((size_t)14,(size_t)1));
 197       sprintf(outfile,"stdin.laz");
 198     }
 199   }
 200   else
 201     check_alloc(ofname=(char*)calloc(strlen(outfile)+10,(size_t)1));
 202
 203
 204   if (column == NULL)
 205     series=(double**)get_multi_series(infile,&length,exclude,&comp,"",dimset,
 206                                       verbosity);
 207   else
 208     series=(double**)get_multi_series(infile,&length,exclude,&comp,column,
 209                                       dimset,verbosity);
 210
 211   check_alloc(interval=(double*)malloc(sizeof(double)*comp));
 212   check_alloc(min=(double*)malloc(sizeof(double)*comp));
 213
 214   maxinterval=maxdvar=0.0;
 215   for (i=0;i<comp;i++) {
 216     rescale_data(series[i],length,&min[i],&interval[i]);
 217     if (interval[i] > maxinterval) maxinterval=interval[i];
 218     variance(series[i],length,&dav,&dvar);
 219     if (dvar > maxdvar)  maxdvar=dvar;
 220   }
 221   alldim=comp*embed;
 222
 223   check_alloc(nmf=(unsigned int*)malloc(sizeof(int)*length));
 224   check_alloc(list=(long*)malloc(sizeof(long)*length));
 225   check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
 226   for (n=0;n<BOX;n++)
 227     check_alloc(box[n]=(long*)malloc(sizeof(long)*BOX));
 228
 229   check_alloc(nf=(long**)malloc(sizeof(long*)*comp));
 230   check_alloc(corr=(double**)malloc(sizeof(double*)*comp));
 231   for (i=0;i<comp;i++) {
 232     check_alloc(nf[i]=(long*)malloc(sizeof(long)*length));
 233     check_alloc(corr[i]=(double*)malloc(sizeof(double)*length));
 234   }
 235
 236   indexes=make_multi_index(comp,embed,delay);
 237
 238   if (epsset)
 239     eps/=maxinterval;
 240   else
 241     eps=1.0/1000.;
 242
 243   if (epsvarset)
 244     eps=epsvar*maxdvar;
 245
 246   for (iter=1;iter<=iterations;iter++) {
 247     make_multi_box2(series,box,list,length,BOX,comp,embed,delay,eps);
 248     for (n=0;n<length;n++) {
 249       for (i=0;i<comp;i++) {
 250         corr[i][n]=0.0;
 251         nf[i][n]=0;
 252       }
 253       nmf[n]=1;
 254     }
 255
 256     check_alloc(hcor=(double*)malloc(sizeof(double)*alldim));
 257     for (n=(embed-1)*delay;n<length;n++)
 258       nmf[n]=correct(n);
 259     free(hcor);
 260
 261     for (n=0;n<length;n++)
 262       for (i=0;i<comp;i++)
 263         if (nf[i][n])
 264           series[i][n]=corr[i][n]/nf[i][n];
 265
 266     if ((verbosity&VER_USR1) && (iter < iterations)) {
 267       sprintf(ofname,"%s.%d",outfile,iter);
 268       test_outfile(ofname);
 269       file=fopen(ofname,"w");
 270       if (verbosity&VER_INPUT)
 271         fprintf(stderr,"Opened %s for writing\n",ofname);
 272       if (stdo && (iter == iterations)) {
 273         if (verbosity&VER_INPUT)
 274           fprintf(stderr,"Writing to stdout\n");
 275       }
 276       for (n=0;n<length;n++) {
 277         if (stdo && (iter == iterations)) {
 278           if (verbosity&VER_USR2) {
 279             for (i=0;i<comp;i++)
 280               fprintf(stdout,"%e ",series[i][n]*interval[i]+min[i]);
 281             fprintf(stdout,"%u\n",nmf[n]);
 282           }
 283           else {
 284             fprintf(stdout,"%e",series[0][n]*interval[0]+min[0]);
 285             for (i=1;i<comp;i++)
 286               fprintf(stdout,"%e ",series[i][n]*interval[i]+min[i]);
 287             fprintf(stdout,"\n");
 288           }
 289         }
 290         if (verbosity&VER_USR2) {
 291           for (i=0;i<comp;i++)
 292             fprintf(file,"%e ",series[i][n]*interval[i]+min[i]);
 293           fprintf(file,"%u\n",nmf[n]);
 294         }
 295         else {
 296           fprintf(file,"%e",series[0][n]*interval[0]+min[0]);
 297           for (i=1;i<comp;i++)
 298             fprintf(file," %e",series[i][n]*interval[i]+min[i]);
 299           fprintf(file,"\n");
 300         }
 301       }
 302       fclose(file);
 303     }
 304     if (iter == iterations) {
 305       if (!stdo || (verbosity&VER_USR1)) {
 306         sprintf(ofname,"%s.%d",outfile,iter);
 307         test_outfile(ofname);
 308         file=fopen(ofname,"w");
 309         if (verbosity&VER_INPUT)
 310           fprintf(stderr,"Opened %s for writing\n",ofname);
 311         if (stdo && (iter == iterations)) {
 312           if (verbosity&VER_INPUT)
 313             fprintf(stderr,"Writing to stdout\n");
 314         }
 315       }
 316       for (n=0;n<length;n++) {
 317         if (stdo) {
 318           if (verbosity&VER_USR2) {
 319             for (i=0;i<comp;i++)
 320               fprintf(stdout,"%e ",series[i][n]*interval[i]+min[i]);
 321             fprintf(stdout,"%u\n",nmf[n]);
 322           }
 323           else {
 324             fprintf(stdout,"%e",series[0][n]*interval[0]+min[0]);
 325             for (i=1;i<comp;i++)
 326               fprintf(stdout," %e",series[i][n]*interval[i]+min[i]);
 327             fprintf(stdout,"\n");
 328           }
 329         }
 330         if (!stdo || (verbosity&VER_USR1)) {
 331           if (verbosity&VER_USR2) {
 332             for (i=0;i<comp;i++)
 333               fprintf(file,"%e ",series[i][n]*interval[i]+min[i]);
 334             fprintf(file,"%u\n",nmf[n]);
 335           }
 336           else {
 337             fprintf(file,"%e",series[0][n]*interval[0]+min[0]);
 338             for (i=1;i<comp;i++)
 339               fprintf(file," %e",series[i][n]*interval[i]+min[i]);
 340             fprintf(file,"\n");
 341           }
 342         }
 343       }
 344       if (!stdo || (verbosity&VER_USR1))
 345         fclose(file);
 346     }
 347   }
 348
 349   /*cleaning up */
 350   for (i=0;i<comp;i++) {
 351     free(series[i]);
 352     free(nf[i]);
 353     free(corr[i]);
 354   }
 355   free(series);
 356   free(nf);
 357   free(corr);
 358
 359   for (i=0;i<2;i++)
 360     free(indexes[i]);
 361   free(indexes);
 362
 363   free(list);
 364   free(nmf);
 365   free(interval);
 366   free(min);
 367
 368   for (i=0;i<BOX;i++)
 369     free(box[i]);
 370   free(box);
 371
 372   if (outfile != NULL)
 373     free(outfile);
 374   if (ofname != NULL)
 375     free(ofname);
 376   if (infile != NULL)
 377     free(infile);
 378   if (column != NULL)
 379     free(column);
 380   /* end cleaning up */
 381
 382   return 0;
 383 }