binaries/src/disembl/Tisean_3.0.1/source_c/boxcount.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: Feb 22, 2006 */
  21 /* Changes:
  22    02/22/06: Remove this strange else in start_box that
  23              did not compile anyways
  24 */
  25
  26 #include <stdio.h>
  27 #include <stdlib.h>
  28 #include <math.h>
  29 #include <string.h>
  30 #include <limits.h>
  31 #include "routines/tsa.h"
  32
  33 #define WID_STR "Estimates the Renyi entropy of Qth order\n\t\
  34 using a partition instead of a covering."
  35
  36 typedef struct {
  37   double *hist;
  38   void *ptr;
  39 } hliste;
  40
  41 unsigned long LENGTH=ULONG_MAX,exclude=0;
  42 unsigned int maxembed=10,dimension=1,DELAY=1,EPSCOUNT=20;
  43 unsigned int verbosity=0xff;
  44 double Q=2.0,EPSMIN=1.e-3,EPSMAX=1.0;
  45 char dimset=0,epsminset=0,epsmaxset=0;
  46 char *outfile=NULL;
  47 char *column=NULL;
  48
  49 int epsi;
  50 unsigned long length;
  51 double EPSFAKTOR;
  52 unsigned int **which_dims;
  53 double *histo;
  54 double **series;
  55
  56 void show_options(char *progname)
  57 {
  58   what_i_do(progname,WID_STR);
  59   fprintf(stderr,"Usage: %s [Options]\n",progname);
  60   fprintf(stderr,"Options:\n");
  61   fprintf(stderr,"\t-l # of datapoints [Default: whole file]\n");
  62   fprintf(stderr,"\t-x # of lines to ignore [Default: %lu]\n",exclude);
  63   fprintf(stderr,"\t-M # of columns,maximal embedding dimension "
  64           "[Default: %u,%u]\n",dimension,maxembed);
  65   fprintf(stderr,"\t-c columns to read  [Default: 1,...,#of compon.]\n");
  66   fprintf(stderr,"\t-d delay [Default: %u]\n",DELAY);
  67   fprintf(stderr,"\t-Q order of the Renyi entropy [Default: %.1f]\n",Q);
  68   fprintf(stderr,"\t-r minimal epsilon [Default: (data interval)/1000]\n");
  69   fprintf(stderr,"\t-R maximal epsilon [Default: data interval]\n");
  70   fprintf(stderr,"\t-# # of epsilons to use [Default: %u]\n",EPSCOUNT);
  71   fprintf(stderr,"\t-o output file name [Default: 'datafile'.box]\n");
  72   fprintf(stderr,"\t-V verbosity level [Default: 1]\n\t\t"
  73           "0='only panic messages'\n\t\t"
  74           "1='+ input/output messages'\n");
  75   fprintf(stderr,"\t-h show these options\n\n");
  76   exit(0);
  77 }
  78
  79 void scan_options(int n,char **in)
  80 {
  81   char *out;
  82
  83   if ((out=check_option(in,n,'l','u')) != NULL)
  84     sscanf(out,"%lu",&LENGTH);
  85   if ((out=check_option(in,n,'c','s')) != NULL)
  86     column=out;
  87   if ((out=check_option(in,n,'x','u')) != NULL)
  88     sscanf(out,"%lu",&exclude);
  89   if ((out=check_option(in,n,'M','2')) != NULL) {
  90     sscanf(out,"%u,%u",&dimension,&maxembed);
  91     dimset=1;
  92   }
  93   if ((out=check_option(in,n,'d','u')) != NULL)
  94     sscanf(out,"%u",&DELAY);
  95   if ((out=check_option(in,n,'Q','f')) != NULL)
  96     sscanf(out,"%lf",&Q);
  97   if ((out=check_option(in,n,'r','f')) != NULL) {
  98     sscanf(out,"%lf",&EPSMIN);
  99     epsminset=1;
 100   }
 101   if ((out=check_option(in,n,'R','f')) != NULL) {
 102     sscanf(out,"%lf",&EPSMAX);
 103     epsmaxset=1;
 104   }
 105   if ((out=check_option(in,n,'#','u')) != NULL)
 106     sscanf(out,"%u",&EPSCOUNT);
 107   if ((out=check_option(in,n,'V','u')) != NULL)
 108     sscanf(out,"%u",&verbosity);
 109   if ((out=check_option(in,n,'o','s')) != NULL)
 110     outfile=out;
 111 }
 112
 113 hliste *make_histo(void)
 114 {
 115   int i;
 116   hliste *element;
 117
 118   check_alloc(element=(hliste*)malloc(sizeof(hliste)));
 119   element->ptr=NULL;
 120   check_alloc(element->hist=(double*)malloc(sizeof(double)*maxembed*dimension));
 121   for (i=0;i<maxembed*dimension;i++)
 122     element->hist[i]=0.0;
 123
 124   return element;
 125 }
 126
 127 void next_dim(int wd,int n,unsigned int *first)
 128 {
 129   int i,which,d1,comp;
 130   double epsinv,norm,p;
 131   unsigned int **act;
 132   int *found,hf;
 133
 134   comp=which_dims[wd][0];
 135   d1=which_dims[wd][1]*DELAY;
 136
 137   epsinv=(double)epsi;
 138   norm=(double)length;
 139
 140   check_alloc(act=(unsigned int**)malloc(epsi*sizeof(int*)));
 141   check_alloc(found=(int*)malloc(epsi*sizeof(int)));
 142
 143   for (i=0;i<epsi;i++) {
 144     found[i]=0;
 145     act[i]=NULL;
 146   }
 147
 148   for (i=0;i<n;i++) {
 149     which=(int)(series[comp][first[i]+d1]*epsinv);
 150     hf= ++found[which];
 151     check_alloc(act[which]=
 152                 realloc((unsigned int*)act[which],hf*sizeof(unsigned int)));
 153     act[which][hf-1]=first[i];
 154   }
 155
 156   for (i=0;i<epsi;i++)
 157     if (found[i]) {
 158       p=(double)(found[i])/(norm);
 159       if (Q == 1.0)
 160         histo[wd] -= p*log(p);
 161       else
 162         histo[wd] += pow(p,Q);
 163     }
 164
 165   if (wd<(maxembed*dimension-1))
 166     for (i=0;i<epsi;i++)
 167       if (found[i])
 168         next_dim(wd+1,found[i],act[i]);
 169
 170   for (i=0;i<epsi;i++)
 171     if (found[i])
 172       free(act[i]);
 173
 174   free(act);
 175   free(found);
 176 }
 177
 178 void start_box(void)
 179 {
 180   int i,which;
 181   double epsinv,norm,p;
 182   unsigned int **act;
 183   int *found,hf;
 184   void next_dim();
 185
 186   epsinv=(double)epsi;
 187   norm=(double)length;
 188
 189   check_alloc(act=(unsigned int**)malloc(epsi*sizeof(int*)));
 190   check_alloc(found=(int*)malloc(epsi*sizeof(int)));
 191
 192   for (i=0;i<epsi;i++) {
 193     found[i]=0;
 194     act[i]=NULL;
 195   }
 196
 197   for (i=0;i<length;i++) {
 198     which=(int)(series[0][i]*epsinv);
 199     hf= ++found[which];
 200     check_alloc(act[which]=
 201                 realloc((unsigned int*)act[which],hf*sizeof(unsigned int)));
 202     act[which][hf-1]=i;
 203   }
 204
 205   for (i=0;i<epsi;i++)
 206     if (found[i]) {
 207       p=(double)(found[i])/(norm);
 208       if (Q == 1.0)
 209         histo[0] -= p*log(p);
 210       else
 211         histo[0] += pow(p,Q);
 212     }
 213
 214   if (1<dimension*maxembed) {
 215     for (i=0;i<epsi;i++) {
 216       if (found[i])
 217         next_dim(1,found[i],act[i]);
 218     }
 219   }
 220   /*
 221   else {
 222     if (1<maxembed)
 223       for (i=0;i<epsi;i++) {
 224         if (found[i])
 225           next_dim(1,found[i],act[i]);
 226       }
 227   }
 228   */
 229
 230   for (i=0;i<epsi;i++)
 231     if (found[i])
 232       free(act[i]);
 233
 234   free(act);
 235   free(found);
 236 }
 237
 238 int main(int argc,char **argv)
 239 {
 240   int i,j,k,count,epsi_old=0,epsi_test;
 241   void *root;
 242   hliste *histo_el;
 243   double *deps,heps;
 244   double min,interval,maxinterval;
 245   char *infile=NULL,stdi=0;
 246   FILE *fHq;
 247
 248
 249   if (scan_help(argc,argv))
 250     show_options(argv[0]);
 251
 252   scan_options(argc,argv);
 253 #ifndef OMIT_WHAT_I_DO
 254   if (verbosity&VER_INPUT)
 255     what_i_do(argv[0],WID_STR);
 256 #endif
 257
 258   infile=search_datafile(argc,argv,NULL,verbosity);
 259   if (infile == NULL)
 260     stdi=1;
 261
 262   if (outfile == NULL) {
 263     if (!stdi) {
 264       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
 265       sprintf(outfile,"%s.box",infile);
 266     }
 267     else {
 268       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
 269       sprintf(outfile,"stdin.box");
 270     }
 271   }
 272   test_outfile(outfile);
 273
 274   if (column == NULL)
 275     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dimension,"",
 276                                       dimset,verbosity);
 277   else
 278     series=(double**)get_multi_series(infile,&LENGTH,exclude,&dimension,
 279                                       column,dimset,verbosity);
 280   maxinterval=0.0;
 281   for (i=0;i<dimension;i++) {
 282     rescale_data(series[i],LENGTH,&min,&interval);
 283     if (interval > maxinterval)
 284       maxinterval=interval;
 285   }
 286   if (epsminset)
 287     EPSMIN /= maxinterval;
 288   if (epsmaxset)
 289     EPSMAX /= maxinterval;
 290   for (i=0;i<dimension;i++) {
 291     for (j=0;j<LENGTH;j++)
 292       if (series[i][j] >= 1.0)
 293         series[i][j] -= EPSMIN/2.0;
 294   }
 295
 296   check_alloc(histo=(double*)malloc(sizeof(double)*maxembed*dimension));
 297   check_alloc(deps=(double*)malloc(sizeof(double)*EPSCOUNT));
 298   check_alloc(which_dims=(unsigned int**)malloc(sizeof(int*)*
 299                                                 maxembed*dimension));
 300   for (i=0;i<maxembed*dimension;i++)
 301     check_alloc(which_dims[i]=(unsigned int*)malloc(sizeof(int)*2));
 302   for (i=0;i<maxembed;i++)
 303     for (j=0;j<dimension;j++) {
 304       which_dims[i*dimension+j][0]=j;
 305       which_dims[i*dimension+j][1]=i;
 306     }
 307
 308   histo_el=make_histo();
 309   root=histo_el;
 310
 311   if (EPSCOUNT >1)
 312     EPSFAKTOR=pow(EPSMAX/EPSMIN,1.0/(double)(EPSCOUNT-1));
 313   else
 314     EPSFAKTOR=1.0;
 315
 316   length=LENGTH-(maxembed-1)*DELAY;
 317
 318   heps=EPSMAX*EPSFAKTOR;
 319
 320   for (k=0;k<EPSCOUNT;k++) {
 321     count++;
 322     for (i=0;i<maxembed*dimension;i++)
 323       histo[i]=0.0;
 324     do {
 325       heps /= EPSFAKTOR;
 326       epsi_test=(int)(1./heps);
 327     } while (epsi_test <= epsi_old);
 328
 329     epsi=epsi_test;
 330     epsi_old=epsi;
 331     deps[k]=heps;
 332
 333     start_box();
 334     histo_el=root;
 335     while (histo_el->ptr != NULL)
 336       histo_el=histo_el->ptr;
 337
 338     for (i=0;i<maxembed*dimension;i++)
 339       if (Q == 1.0)
 340         histo_el->hist[i]=histo[i];
 341       else
 342         histo_el->hist[i]=log(histo[i])/(1.0-Q);
 343
 344     histo_el->ptr=make_histo();
 345     histo_el=histo_el->ptr;
 346     fHq=fopen(outfile,"w");
 347     if (verbosity&VER_INPUT)
 348       fprintf(stderr,"Opened %s for writing\n",outfile);
 349
 350     for (i=0;i<maxembed*dimension;i++) {
 351       fprintf(fHq,"#component = %d embedding = %d\n",which_dims[i][0]+1,
 352               which_dims[i][1]+1);
 353       histo_el=root;
 354       for (j=0;j<=k;j++) {
 355         if (i == 0)
 356           fprintf(fHq,"%e %e %e\n",deps[j]*maxinterval,
 357                   histo_el->hist[i],histo_el->hist[i]);
 358         else
 359           fprintf(fHq,"%e %e %e\n",deps[j]*maxinterval,
 360                   histo_el->hist[i],histo_el->hist[i]-histo_el->hist[i-1]);
 361         histo_el=histo_el->ptr;
 362       }
 363       fprintf(fHq,"\n");
 364     }
 365     fclose(fHq);
 366   }
 367
 368   return 0;
 369 }