IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / binaries / 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