website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_c/xcor.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 4, 1999 */
  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 crosscorrelations of two data sets\n\t\
  29 given as two columns of one file."
  30
  31 char *columns=NULL,*outfile=NULL,stout=1;
  32 unsigned long length=ULONG_MAX,exclude=0;
  33 long tau=100;
  34 unsigned int verbosity=0xff;
  35 double *array1,*array2;
  36 char *infile=NULL;
  37
  38 void show_options(char *progname)
  39 {
  40   what_i_do(progname,WID_STR);
  41   fprintf(stderr," Usage: %s [Options]\n",progname);
  42   fprintf(stderr," Options:\n");
  43   fprintf(stderr,"Everything not being a valid option will be interpreted"
  44           " as a possible"
  45           " datafile.\nIf no datafile is given stdin is read. Just - also"
  46           " means stdin\n");
  47   fprintf(stderr,"\t-l length [default is whole file]\n");
  48   fprintf(stderr,"\t-x # of lines to be ignored [default 0]\n");
  49   fprintf(stderr,"\t-c which columns (separated by commas) [default is 1,2]\n");
  50   fprintf(stderr,"\t-D corrlength  [default is 100]\n");
  51   fprintf(stderr,"\t-o output_file  [default is 'datafile'.crc; no -o"
  52   " means stdout]\n");
  53   fprintf(stderr,"\t-V verbosity level [default is 1]\n\t\t"
  54           "0='only panic messages'\n\t\t"
  55           "1='+ input/output messages'\n");
  56   fprintf(stderr,"\t-h show these options\n");
  57   fprintf(stderr,"\n");
  58   exit(0);
  59 }
  60
  61 void scan_options(int argc,char **argv)
  62 {
  63   char *out;
  64
  65   if ((out=check_option(argv,argc,'l','u')) != NULL)
  66     sscanf(out,"%lu",&length);
  67   if ((out=check_option(argv,argc,'x','u')) != NULL)
  68     sscanf(out,"%lu",&exclude);
  69   if ((out=check_option(argv,argc,'c','s')) != NULL)
  70     columns=out;
  71   if ((out=check_option(argv,argc,'D','u')) != NULL)
  72     sscanf(out,"%ld",&tau);
  73   if ((out=check_option(argv,argc,'V','u')) != NULL)
  74     sscanf(out,"%u",&verbosity);
  75   if ((out=check_option(argv,argc,'o','o')) != NULL) {
  76     stout=0;
  77     if (strlen(out) > 0)
  78       outfile=out;
  79   }
  80 }
  81
  82 double corr(long i)
  83 {
  84   unsigned long count=0;
  85   long j,hi;
  86   double c=0.0;
  87
  88   for (j=0;j<length;j++) {
  89     hi=j+i;
  90     if ((hi >= 0) && (hi < length)) {
  91       count++;
  92       c += array1[j]*array2[hi];
  93     }
  94   }
  95   return c/(double)count;
  96 }
  97
  98 int main(int argc,char** argv)
  99 {
 100   char stdi=0;
 101   long i;
 102   unsigned int dummy=2;
 103   FILE *fout=NULL;
 104   double **both;
 105   double av1,var1,av2,var2;
 106
 107   if (scan_help(argc,argv))
 108     show_options(argv[0]);
 109
 110   scan_options(argc,argv);
 111 #ifndef OMIT_WHAT_I_DO
 112   if (verbosity&VER_INPUT)
 113     what_i_do(argv[0],WID_STR);
 114 #endif
 115
 116   infile=search_datafile(argc,argv,0L,verbosity);
 117   if (infile == NULL)
 118     stdi=1;
 119
 120   if (outfile == NULL) {
 121     if (!stdi) {
 122       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
 123       strcpy(outfile,infile);
 124       strcat(outfile,".ccr");
 125     }
 126     else {
 127       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
 128       strcpy(outfile,"stdin.ccr");
 129     }
 130   }
 131   if (!stout)
 132     test_outfile(outfile);
 133
 134   if (columns == NULL)
 135     both=(double**)get_multi_series(infile,&length,exclude,&dummy,"",(char)1,
 136                                     verbosity);
 137   else
 138     both=(double**)get_multi_series(infile,&length,exclude,&dummy,columns,
 139                                     (char)1,verbosity);
 140
 141   array1=both[0];
 142   array2=both[1];
 143
 144   if (tau >= length)
 145     tau=length-1;
 146
 147   variance(array1,length,&av1,&var1);
 148   variance(array2,length,&av2,&var2);
 149
 150   for (i=0;i<length;i++) {
 151     array1[i] -= av1;
 152     array2[i] -= av2;
 153   }
 154
 155   if (!stout) {
 156     fout=fopen(outfile,"w");
 157     if (verbosity&VER_INPUT)
 158       fprintf(stderr,"Opened %s for writing\n",outfile);
 159     fprintf(fout,"# average of first comp.=%e\n",av1);
 160     fprintf(fout,"# standard deviation of first comp.=%e\n",var1);
 161     fprintf(fout,"# average of sec. comp.=%e\n",av2);
 162     fprintf(fout,"# standard deviation of sec. comp.=%e\n",var2);
 163   }
 164   else {
 165     if (verbosity&VER_INPUT)
 166       fprintf(stderr,"Writing to stdout\n");
 167     fprintf(stdout,"# average of first comp.=%e\n",av1);
 168     fprintf(stdout,"# standard deviation of first comp.=%e\n",var1);
 169     fprintf(stdout,"# average of sec. comp.=%e\n",av2);
 170     fprintf(stdout,"# standard deviation of sec. comp.=%e\n",var2);
 171   }
 172
 173   for (i= -tau;i<=tau;i++)
 174     if (!stout) {
 175       fprintf(fout,"%ld %e\n",i,corr(i)/var1/var2);
 176       fflush(fout);
 177     }
 178     else {
 179       fprintf(stdout,"%ld %e\n",i,corr(i)/var1/var2);
 180       fflush(stdout);
 181     }
 182   if (!stout)
 183     fclose(fout);
 184
 185   return 0;
 186 }