/******************************************************************* minfo.c -- Eric R. Weeks -- started 2/28/97 does the mutual information algorithm discussed by Fraser & Swinney (Phys Rev A 33 (1986) p1134-1140) v01: 2-28-97: taken from shell.c (7/1/96) quicksort routine taken from sane.c (original version) v02: 2-28-97: revised sorting for s[] (different than q[]) sped up math v03: 2-28-97: add in tau loop 3-01-97: fix for variable number of input; add -b option v04: 3-01-97: take out chi2 tests for substructure v05: 3-01-97: realize that with chi2 tests taken out, number() function doesn't need to be called very often. remove a[] and b[][] arrays! Much faster now. This program is public domain, although please leave my name and email address attached. email: weeks@physics.emory.edu web: http://www.physics.emory.edu/~weeks/ explanation of how to use this program: http://www.physics.emory.edu/~weeks/software/minfo.html *******************************************************************/ #include #include #include #define PI 3.14159265358979323846264338328 #define EE 2.71828182845904523536 #define MAXNUM 100000 #define KMAX 25 int s[MAXNUM],q[MAXNUM]; int numpts,n0,mmax; float pop[MAXNUM]; int pindex[MAXNUM]; int pow2[KMAX]; int bin; int globalivalue[4]; float logtwo; int number(int karray2[],int m2); float findisq(); float ffunct(int kmarray[],int m); void saneqsort(int p,int r); int qpartition_neurons(int p,int r); main(argc,argv) int argc; char **argv; { int c; double x,y; FILE *fp; extern int optind; extern char *optarg; float array[MAXNUM]; int taumax,i,j,tau; float info; bin = -1; taumax = 20; while ((c = getopt(argc, argv, "ht:b:")) != EOF) switch (c) { case 'h': fprintf(stderr,"Usage: %s [options] [<] [file]\n\n",argv[0]); fprintf(stderr," -h : this help message\n"); fprintf(stderr," -t # : set taumax, max tau to find mut. info [%d]\n",taumax); fprintf(stderr," -b # : bin data into 2^# bins [default: all possible]\n"); exit(1); break; case 't': taumax = atoi(optarg); break; case 'b': bin = atoi(optarg); break; } argc -= (optind-1) ; argv += (optind-1) ; fp = (argc > 1) ? fopen(*++argv, "r") : stdin; numpts = 0; while ((fscanf(fp,"%lf",&x)) == 1) { array[numpts] = x; ++numpts; } /*----- done reading in data -----*/ logtwo = 1.0/log(2.0); n0=1; while ((n0+taumax)<=numpts) n0 *= 2; n0 /= 2; /* n0 = numpts - taumax; */ j = n0; for (i=0;i0) { los = 0;loq = 0; his = n0; hiq = n0; for (i=1;i<=m2;i++) { if (karray2[i]%2==0) his -= pow2[i]; else los += pow2[i]; if (karray2[i]<2) hiq -= pow2[i]; else loq += pow2[i]; } ivalue = 0; for (i=los;i=loq)&&(j x); if (i < j) { /* here's where the in-place swap takes place */ /* fortunately pindex[] stores the *original* location */ temp = pop[i]; pop[i] = pop[j]; pop[j] = temp; tempi = pindex[i]; pindex[i] = pindex[j]; pindex[j] = tempi; } else return j; } }