website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_c/mutual.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 20, 2000 */
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <math.h>
  24 #include <limits.h>
  25 #include <string.h>
  26 #include "routines/tsa.h"
  27
  28 #define WID_STR "Estimates the time delayed mutual information\n\t\
  29 of the data set"
  30
  31
  32 char *file_out=NULL,stout=1;
  33 char *infile=NULL;
  34 unsigned long length=ULONG_MAX,exclude=0;
  35 unsigned int column=1;
  36 unsigned int verbosity=0xff;
  37 long partitions=16,corrlength=20;
  38 long *array,*h1,*h11,**h2;
  39
  40 void show_options(char *progname)
  41 {
  42   what_i_do(progname,WID_STR);
  43   fprintf(stderr," Usage: %s [Options]\n\n",progname);
  44   fprintf(stderr," Options:\n");
  45   fprintf(stderr,"Everything not being a valid option will be interpreted"
  46           " as a possible"
  47           " datafile.\nIf no datafile is given stdin is read. Just - also"
  48           " means stdin\n");
  49   fprintf(stderr,"\t-l # of points to be used [Default is all]\n");
  50   fprintf(stderr,"\t-x # of lines to be ignored [Default is 0]\n");
  51   fprintf(stderr,"\t-c column to read  [Default is 1]\n");
  52   fprintf(stderr,"\t-b # of boxes [Default is 16]\n");
  53   fprintf(stderr,"\t-D max. time delay [Default is 20]\n");
  54   fprintf(stderr,"\t-o output file [-o without name means 'datafile'.mut;"
  55           "\n\t\tNo -o means write to stdout]\n");
  56   fprintf(stderr,"\t-V verbosity level [Default is 1]\n\t\t"
  57           "0='only panic messages'\n\t\t"
  58           "1='+ input/output messages'\n");
  59   fprintf(stderr,"\t-h  show these options\n");
  60   fprintf(stderr,"\n");
  61   exit(0);
  62 }
  63
  64 void scan_options(int n,char** in)
  65 {
  66   char *out;
  67
  68   if ((out=check_option(in,n,'l','u')) != NULL)
  69     sscanf(out,"%lu",&length);
  70   if ((out=check_option(in,n,'x','u')) != NULL)
  71     sscanf(out,"%lu",&exclude);
  72   if ((out=check_option(in,n,'c','u')) != NULL)
  73     sscanf(out,"%u",&column);
  74   if ((out=check_option(in,n,'b','u')) != NULL)
  75     sscanf(out,"%lu",&partitions);
  76   if ((out=check_option(in,n,'D','u')) != NULL)
  77     sscanf(out,"%lu",&corrlength);
  78   if ((out=check_option(in,n,'V','u')) != NULL)
  79     sscanf(out,"%u",&verbosity);
  80   if ((out=check_option(in,n,'o','o')) != NULL) {
  81     stout=0;
  82     if (strlen(out) > 0)
  83       file_out=out;
  84   }
  85 }
  86
  87 double make_cond_entropy(long t)
  88 {
  89   long i,j,hi,hii,count=0;
  90   double hpi,hpj,pij,cond_ent=0.0,norm;
  91
  92   for (i=0;i<partitions;i++) {
  93     h1[i]=h11[i]=0;
  94     for (j=0;j<partitions;j++)
  95       h2[i][j]=0;
  96   }
  97   for (i=0;i<length;i++)
  98     if (i >= t) {
  99       hii=array[i];
 100       hi=array[i-t];
 101       h1[hi]++;
 102       h11[hii]++;
 103       h2[hi][hii]++;
 104       count++;
 105     }
 106
 107   norm=1.0/(double)count;
 108   cond_ent=0.0;
 109
 110   for (i=0;i<partitions;i++) {
 111     hpi=(double)(h1[i])*norm;
 112     if (hpi > 0.0) {
 113       for (j=0;j<partitions;j++) {
 114         hpj=(double)(h11[j])*norm;
 115         if (hpj > 0.0) {
 116           pij=(double)h2[i][j]*norm;
 117           if (pij > 0.0)
 118             cond_ent += pij*log(pij/hpj/hpi);
 119         }
 120       }
 121     }
 122   }
 123
 124   return cond_ent;
 125 }
 126
 127 int main(int argc,char** argv)
 128 {
 129   char stdi=0;
 130   long tau,i;
 131   double *series,min,interval,shannon;
 132   FILE *file;
 133
 134   if (scan_help(argc,argv))
 135     show_options(argv[0]);
 136
 137   scan_options(argc,argv);
 138 #ifndef OMIT_WHAT_I_DO
 139   if (verbosity&VER_INPUT)
 140     what_i_do(argv[0],WID_STR);
 141 #endif
 142
 143   infile=search_datafile(argc,argv,&column,verbosity);
 144   if (infile == NULL)
 145     stdi=1;
 146
 147   if (file_out == NULL) {
 148     if (!stdi) {
 149       check_alloc(file_out=(char*)calloc(strlen(infile)+5,(size_t)1));
 150       strcpy(file_out,infile);
 151       strcat(file_out,".mut");
 152     }
 153     else {
 154       check_alloc(file_out=(char*)calloc((size_t)10,(size_t)1));
 155       strcpy(file_out,"stdin.mut");
 156     }
 157   }
 158   if (!stout)
 159     test_outfile(file_out);
 160
 161   series=(double*)get_series(infile,&length,exclude,column,verbosity);
 162   rescale_data(series,length,&min,&interval);
 163
 164   check_alloc(h1=(long *)malloc(sizeof(long)*partitions));
 165   check_alloc(h11=(long *)malloc(sizeof(long)*partitions));
 166   check_alloc(h2=(long **)malloc(sizeof(long *)*partitions));
 167   for (i=0;i<partitions;i++)
 168     check_alloc(h2[i]=(long *)malloc(sizeof(long)*partitions));
 169   check_alloc(array=(long *)malloc(sizeof(long)*length));
 170   for (i=0;i<length;i++)
 171     if (series[i] < 1.0)
 172       array[i]=(long)(series[i]*(double)partitions);
 173     else
 174       array[i]=partitions-1;
 175   free(series);
 176
 177   shannon=make_cond_entropy(0);
 178   if (corrlength >= length)
 179     corrlength=length-1;
 180
 181   if (!stout) {
 182     file=fopen(file_out,"w");
 183     if (verbosity&VER_INPUT)
 184       fprintf(stderr,"Opened %s for writing\n",file_out);
 185     fprintf(file,"#shannon= %e\n",shannon);
 186     fprintf(file,"%d %e\n",0,shannon);
 187     for (tau=1;tau<=corrlength;tau++) {
 188       fprintf(file,"%ld %e\n",tau,make_cond_entropy(tau));
 189       fflush(file);
 190     }
 191     fclose(file);
 192   }
 193   else {
 194     if (verbosity&VER_INPUT)
 195       fprintf(stderr,"Writing to stdout\n");
 196     fprintf(stdout,"#shannon= %e\n",shannon);
 197     fprintf(stdout,"%d %e\n",0,shannon);
 198     for (tau=1;tau<=corrlength;tau++) {
 199       fprintf(stdout,"%ld %e\n",tau,make_cond_entropy(tau));
 200       fflush(stdout);
 201     }
 202   }
 203
 204   return 0;
 205 }
 206