/*
  corr.c
    print correlation coeff's for 2 input files of 'short' values

	Copyright (C) 1999 Stanley J. Brooks <stabro@megsinet.net> 

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 
*/

#include <string.h>
// for open,read,write:
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <errno.h>
#include <unistd.h>
// for clock/time:
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define BLEN 0x20000
#define HLEN 0x1000

static u_int32_t dist[0x10000];
#define dis0 (dist+0x8000)

static short buff1[BLEN];
static short buff2[BLEN];
static short hist[HLEN];

inline static double square(double x) {return (x*x);}

static void readin(int fd, short *p, int bs)
{
	int r;
 	do {
		r=read(fd,(char*)p,2*bs);
	}while (r==-1 && errno==EINTR);
	if (r==-1 || r!=2*bs) {
		perror("Input error");
		exit(1);
	}
}

int main(int argct, char **argv)
{
	int fd1,fd2;
	int avgn=0,assay=0,code=0;
	int hx=0,hsum=0;
	int len,len1,len2,x,r;
	char *fnam1,*fnam2;
	double v1=0,v2=0,s1=0;

	/*
	 * Parse the command line arguments
	 */
	while (*++argv && **argv == '-')
		while (*(++*argv))
			switch (**argv) {
			case 'a':
				assay=1;
				break;
			case 'c':
				code=1;
				break;
			case 'h':
				fprintf(stderr,"./stat [-n avgn] <file>\n");
				break;
			case 'n':
				if (!*++*argv) { // end of this param string
					if (!argv[1]) break;
					++argv;
				}
				avgn=atoi(*argv);
				*argv += strlen(*argv)-1;	// skip to end-1 of this param string
				break;
			}

	if (avgn<=0) avgn=1; else if (avgn>HLEN) avgn=HLEN;
	bzero(hist,sizeof(hist));

	fnam1=*argv++;
	fd1=open(fnam1,O_RDWR);
	if (fd1<0) {
		fprintf(stderr,"Open: %s %s\n",fnam1,strerror(errno)); return(1);
	}
	len1=lseek(fd1,0,SEEK_END);
	r=lseek(fd1,0,SEEK_SET);

	fnam2=*argv++;
	fd2=open(fnam2,O_RDWR);
	if (fd2<0) {
		fprintf(stderr,"Open: %s %s\n",fnam2,strerror(errno)); return(1);
	}
	len2=lseek(fd2,0,SEEK_END);
	r=lseek(fd2,0,SEEK_SET);

	bzero(dist,sizeof(dist));
	len = (len1<len2)? len1:len2;
	len /= 2; /* now len is number of shorts data to read */
	for (x=0; x<len; ){
		int bs;
		int j;
		double c11=0, c12=0, c22=0;
		
		bs=len-x; if (bs>BLEN) bs=BLEN; /* number of shorts to read */

		readin(fd1,buff1,bs);
		readin(fd2,buff2,bs);
		
		for (j=0; j<bs; j++) {
			c11 += buff1[j]*buff1[j];
			c12 += buff1[j]*buff2[j];
			c22 += buff2[j]*buff2[j];
		}
		c11 /= bs;
		c12 /= bs;
		c22 /= bs;

		{
			double d11=(c11+2*c12+c22)/4;
			double d22=(c11-2*c12+c22)/4;
			double d12=(c11-c22)/4;

			printf("%8.1f%8.1f cf=%f",sqrt(c11),sqrt(c22),c12/sqrt(c11*c22));
			printf(" | %8.1f%8.1f cf=%f\n",sqrt(d11),sqrt(d22),d12/sqrt(d11*d22));
		}

		for (j=0; j<bs; j++) {
			int y;
			//y=(abs(buff1[j])<abs(buff2[j]))? buff1[j]:buff2[j];
			y=(buff1[j]+buff2[j]+1)/2;
			hsum -= hist[hx];
			hsum += y;
			hist[hx] = y;
			if (++hx == avgn) hx=0;
			s1 += abs(y);
			v1 += hsum;
			v2 += square(hsum);
			y = (hsum+avgn/2)/avgn;
			dis0[y]++;
		}
		
		x += bs;
	}
	v1 /= len;
	v2 /= len;
	printf("%8d %5d=n %10.2f=avg %10.2f=rms\n",len,avgn,v1/avgn,sqrt(v2-v1*v1)/avgn);

	if (code) {
		int32_t j,tot,cut;
		cut = len/2;
		for (tot=0,j=0; j<0x8000; j++) {
			int32_t tot1;
			tot1 = tot+dis0[j];
			if (j!=0) tot1 += dis0[-j];
			if (tot1>cut && tot<=cut) {
				printf(" |%d| %d of %d\n",j-1,tot,cut);
				cut += (len-cut)/2;
			}
			tot=tot1;
		}
	}

	if (assay)
		for (x=-0x8000; x<=0x7fff; x++) {
			printf("%6d %6d\n",x,dis0[x]);
		}
	
	return 0;
}