binaries/src/disembl/Tisean_3.0.1/source_c/false_nearest.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: Dec 10, 2005 */
  21 /*Changes:
  22   12/10/05: It's multivariate now
  23   12/16/05: Scaled <eps> and sigma(eps)
  24 */
  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 "Determines the fraction of false nearest neighbors."
  34
  35 char *outfile=NULL;
  36 char *infile=NULL;
  37 char stdo=1,dimset=0;
  38 char *column=NULL;
  39 unsigned long length=ULONG_MAX,exclude=0,theiler=0;
  40 unsigned int delay=1,maxdim=5,minemb=1;
  41 unsigned int comp=1,maxemb=5;
  42 unsigned int verbosity=0xff;
  43 double rt=2.0;
  44 double eps0=1.0e-5;
  45 double **series;
  46 double aveps,vareps;
  47 double varianz;
  48
  49 #define BOX 1024
  50 int ibox=BOX-1;
  51 long **box,*list;
  52 unsigned int *vcomp,*vemb;
  53 unsigned long toolarge;
  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 [default: whole file]\n");
  65   fprintf(stderr,"\t-x # of lines to ignore [default: 0]\n");
  66   fprintf(stderr,"\t-c columns to read [default: 1]\n");
  67   fprintf(stderr,"\t-m min. test embedding dimension [default: %u]\n",minemb);
  68   fprintf(stderr,"\t-M # of components,max. emb. dim. [default: %u,%u]\n",
  69           comp,maxemb);
  70   fprintf(stderr,"\t-d delay [default: 1]\n");
  71   fprintf(stderr,"\t-f escape factor [default: %.2lf]\n",rt);
  72   fprintf(stderr,"\t-t theiler window [default: 0]\n");
  73   fprintf(stderr,"\t-o output file [default: 'datafile'.fnn; without -o"
  74           " stdout]\n");
  75   fprintf(stderr,"\t-V verbosity level [default: 3]\n\t\t"
  76           "0='only panic messages'\n\t\t"
  77           "1='+ input/output messages'\n\t\t"
  78           "2='+ information about the current state\n");
  79   fprintf(stderr,"\t-h show these options\n");
  80   exit(0);
  81 }
  82
  83 void scan_options(int n,char **in)
  84 {
  85   char *out;
  86
  87   if ((out=check_option(in,n,'l','u')) != NULL)
  88     sscanf(out,"%lu",&length);
  89   if ((out=check_option(in,n,'x','u')) != NULL)
  90     sscanf(out,"%lu",&exclude);
  91   if ((out=check_option(in,n,'c','s')) != NULL)
  92     column=out;
  93   if ((out=check_option(in,n,'m','u')) != NULL)
  94     sscanf(out,"%u",&minemb);
  95   if ((out=check_option(in,n,'M','2')) != NULL) {
  96     sscanf(out,"%u,%u",&comp,&maxemb);
  97     maxdim=comp*(maxemb+1);
  98     dimset=1;
  99   }
 100   if ((out=check_option(in,n,'d','u')) != NULL)
 101     sscanf(out,"%u",&delay);
 102   if ((out=check_option(in,n,'f','f')) != NULL)
 103     sscanf(out,"%lf",&rt);
 104   if ((out=check_option(in,n,'t','u')) != NULL)
 105     sscanf(out,"%lu",&theiler);
 106   if ((out=check_option(in,n,'V','u')) != NULL)
 107     sscanf(out,"%u",&verbosity);
 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 void mmb(unsigned int hdim,unsigned int hemb,double eps)
 116 {
 117   unsigned long i;
 118   long x,y;
 119
 120   for (x=0;x<BOX;x++)
 121     for (y=0;y<BOX;y++)
 122       box[x][y] = -1;
 123
 124   for (i=0;i<length-(maxemb+1)*delay;i++) {
 125     x=(long)(series[0][i]/eps)&ibox;
 126     y=(long)(series[hdim][i+hemb]/eps)&ibox;
 127     list[i]=box[x][y];
 128     box[x][y]=i;
 129   }
 130 }
 131
 132 char find_nearest(long n,unsigned int dim,double eps)
 133 {
 134   long x,y,x1,x2,y1,i,i1,ic,ie;
 135   long element,which= -1;
 136   double dx,maxdx,mindx=1.1,hfactor,factor;
 137
 138   ic=vcomp[dim];
 139   ie=vemb[dim];
 140   x=(long)(series[0][n]/eps)&ibox;
 141   y=(long)(series[ic][n+ie]/eps)&ibox;
 142
 143   for (x1=x-1;x1<=x+1;x1++) {
 144     x2=x1&ibox;
 145     for (y1=y-1;y1<=y+1;y1++) {
 146       element=box[x2][y1&ibox];
 147       while (element != -1) {
 148         if (labs(element-n) > theiler) {
 149           maxdx=fabs(series[0][n]-series[0][element]);
 150           for (i=1;i<=dim;i++) {
 151             ic=vcomp[i];
 152             i1=vemb[i];
 153             dx=fabs(series[ic][n+i1]-series[ic][element+i1]);
 154             if (dx > maxdx)
 155               maxdx=dx;
 156           }
 157           if ((maxdx < mindx) && (maxdx > 0.0)) {
 158             which=element;
 159             mindx=maxdx;
 160           }
 161         }
 162         element=list[element];
 163       }
 164     }
 165   }
 166
 167   if ((which != -1) && (mindx <= eps) && (mindx <= varianz/rt)) {
 168     aveps += mindx;
 169     vareps += mindx*mindx;
 170     factor=0.0;
 171     for (i=1;i<=comp;i++) {
 172       ic=vcomp[dim+i];
 173       ie=vemb[dim+i];
 174       hfactor=fabs(series[ic][n+ie]-series[ic][which+ie])/mindx;
 175       if (hfactor > factor)
 176         factor=hfactor;
 177     }
 178     if (factor > rt)
 179       toolarge++;
 180     return 1;
 181   }
 182   return 0;
 183 }
 184
 185 int main(int argc,char **argv)
 186 {
 187   char stdi=0;
 188   FILE *file=NULL;
 189   double min,inter=0.0,ind_inter,epsilon,av,ind_var;
 190   char *nearest,alldone;
 191   long i;
 192   unsigned int dim,emb;
 193   unsigned long donesofar;
 194
 195   if (scan_help(argc,argv))
 196     show_options(argv[0]);
 197
 198   scan_options(argc,argv);
 199 #ifndef OMIT_WHAT_I_DO
 200   if (verbosity&VER_INPUT)
 201     what_i_do(argv[0],WID_STR);
 202 #endif
 203
 204   infile=search_datafile(argc,argv,NULL,verbosity);
 205   if (infile == NULL)
 206     stdi=1;
 207
 208   if (outfile == NULL) {
 209     if (!stdi) {
 210       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
 211       strcpy(outfile,infile);
 212       strcat(outfile,".fnn");
 213     }
 214     else {
 215       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
 216       strcpy(outfile,"stdin.fnn");
 217     }
 218   }
 219   if (!stdo)
 220     test_outfile(outfile);
 221
 222   if (column == NULL)
 223     series=(double**)get_multi_series(infile,&length,exclude,&comp,"",dimset,
 224                                       verbosity);
 225   else
 226     series=(double**)get_multi_series(infile,&length,exclude,&comp,column,
 227                                       dimset,verbosity);
 228
 229   for (i=0;i<comp;i++) {
 230     rescale_data(series[i],length,&min,&ind_inter);
 231     variance(series[i],length,&av,&ind_var);
 232     if (i == 0) {
 233       varianz=ind_var;
 234       inter=ind_inter;
 235     }
 236     else {
 237       varianz=(varianz>ind_var)?ind_var:varianz;
 238       inter=(inter<ind_inter)?ind_inter:inter;
 239     }
 240   }
 241
 242   check_alloc(list=(long*)malloc(sizeof(long)*length));
 243   check_alloc(nearest=(char*)malloc(length));
 244   check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
 245   for (i=0;i<BOX;i++)
 246     check_alloc(box[i]=(long*)malloc(sizeof(long)*BOX));
 247
 248   if (!stdo) {
 249     file=fopen(outfile,"w");
 250     if (verbosity&VER_INPUT)
 251       fprintf(stderr,"Opened %s for writing\n",outfile);
 252   }
 253   else {
 254     if (verbosity&VER_INPUT)
 255       fprintf(stderr,"Writing to stdout\n");
 256   }
 257   check_alloc(vcomp=(unsigned int*)malloc(sizeof(int)*(maxdim)));
 258   check_alloc(vemb=(unsigned int*)malloc(sizeof(int)*(maxdim)));
 259   for (i=0;i<maxdim;i++) {
 260     if (comp == 1) {
 261       vcomp[i]=0;
 262       vemb[i]=i;
 263     }
 264     else {
 265       vcomp[i]=i%comp;
 266       vemb[i]=(i/comp)*delay;
 267     }
 268   }
 269   for (emb=minemb;emb<=maxemb;emb++) {
 270     dim=emb*comp-1;
 271     epsilon=eps0;
 272     toolarge=0;
 273     alldone=0;
 274     donesofar=0;
 275     aveps=0.0;
 276     vareps=0.0;
 277     for (i=0;i<length;i++)
 278       nearest[i]=0;
 279     if (verbosity&VER_USR1)
 280       fprintf(stderr,"Start for dimension=%u\n",dim+1);
 281     while (!alldone && (epsilon < 2.*varianz/rt)) {
 282       alldone=1;
 283       mmb(vcomp[dim],vemb[dim],epsilon);
 284       for (i=0;i<length-maxemb*delay;i++)
 285         if (!nearest[i]) {
 286           nearest[i]=find_nearest(i,dim,epsilon);
 287           alldone &= nearest[i];
 288           donesofar += (unsigned long)nearest[i];
 289         }
 290       if (verbosity&VER_USR1)
 291         fprintf(stderr,"Found %lu up to epsilon=%e\n",donesofar,epsilon*inter);
 292       epsilon*=sqrt(2.0);
 293       if (!donesofar)
 294         eps0=epsilon;
 295     }
 296     if (donesofar == 0) {
 297       fprintf(stderr,"Not enough points found!\n");
 298       exit(FALSE_NEAREST_NOT_ENOUGH_POINTS);
 299     }
 300     aveps *= (1./(double)donesofar);
 301     vareps *= (1./(double)donesofar);
 302     if (stdo) {
 303       fprintf(stdout,"%u %e %e %e\n",dim+1,(double)toolarge/(double)donesofar,
 304               aveps*inter,sqrt(vareps)*inter);
 305       fflush(stdout);
 306     }
 307     else {
 308       fprintf(file,"%u %e %e %e\n",dim+1,(double)toolarge/(double)donesofar,
 309               aveps*inter,sqrt(vareps)*inter);
 310       fflush(file);
 311     }
 312   }
 313   if (!stdo)
 314     fclose(file);
 315
 316   if (infile != NULL)
 317     free(infile);
 318   if (outfile != NULL)
 319     free(outfile);
 320   free(series);
 321   free(list);
 322   free(nearest);
 323   for (i=0;i<BOX;i++)
 324     free(box[i]);
 325   free(box);
 326
 327   return 0;
 328 }