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