/*
 * fhist.c	Jim Piper	26-09-86
 */

#include <stdio.h>
#include "../h/wstruct.h"
#include "../h/chromanal.h"
#include <math.h>

#define NBUCKETS 20



main(argc,argv)
char **argv;
{
	int clasi,clasj,ndim,g,h,i,j,k,err,sinverr;
	int nsamp,class,nclass,clnsamp[MAXCLASS];
	int bucket[MAXCLASS][MAXDIM][NBUCKETS];
	int TRAIN,STATS,DERIV,OUTLIERS;
	double maxcor;
	double w,v[MAXDIM];
	double ss,t,msd,cc,cd,sqrtdet;
	double cp[MAXDIM][MAXDIM],s[MAXDIM],m[MAXDIM],psd[MAXDIM];
	double clm[MAXCLASS][MAXDIM],clsd[MAXCLASS][MAXDIM];
	double cls[MAXCLASS][MAXDIM];
	double clcp[MAXCLASS][MAXDIM][MAXDIM];
	FILE *proto, *outlf, *infile;
	double diff,outlval;
	char VCVTEMP[100];

	/* mode switches */
	DERIV = 0;
	OUTLIERS = 1;
	proto = NULL;
	while (argc > 1 && argv[1][0] == '-') {
		switch(argv[1][1]) {
		case 'd':
			DERIV++;
			break;
		case 'P':
			proto = fopen(argv[2],"r");
			if (proto == NULL) {
				fprintf(stderr,"-P can't open file %s\n",argv[2]);
				exit(1);
			}
			argc--;
			argv++;
			fscanf(proto,"%d %d\n",&nclass,&ndim);
			fclose(proto);
			break;
		default:
			break;
		}
		argc--;
		argv++;
	}

	/* number of dimensions and classes */
	if (proto == NULL)
		ndim = MAXDIM;
	if (ndim > MAXDIM) {
		printf("Too many dimensions\n");
		exit(1);
	}
	if (proto == NULL)
		nclass = MAXCLASS;
	if (nclass > MAXCLASS) {
		printf("Too many classes\n");
		exit(1);
	}

	/* temporary file if outlier removal in force */
	if (OUTLIERS) {
		sprintf(VCVTEMP,"/tmp/fhist%d",getpid());
		outlf = fopen(VCVTEMP,"w");
		if (outlf == NULL) {
			fprintf(stderr,"Cannot open temp file\n");
			exit(1);
		}
		/* clear arrays */
		nsamp = 0;
		for (i=0; i<ndim; i++) {
			s[i] = 0.0;
			for (j=0; j<ndim; j++)
				cp[j][i] = 0.0;
			for (class=0; class<nclass; class++) {
				clnsamp[class] = 0;
				cls[class][i] = 0.0;
				for (j=0; j<ndim; j++)
					clcp[class][j][i] = 0.0;
			}
		}

		/* get data vector, check class, check for derived objects */
		nsamp = 0;
		while (scanf("%d",&class) > 0) {
			fprintf(outlf,"%d ",class);
			class--;
			for (i=0; i<ndim; i++) {
				scanf("%e",&v[i]);
				fprintf(outlf,"%e ",v[i]);
			}
			fprintf(outlf,"\n");
			if (DERIV != 0 && class >= 100)
				class -= 100;

		/* set up vector sum and matrix of sum of products */
			if (class >= 0 && class < nclass) {
				nsamp++;
				clnsamp[class]++;
				for (i=0; i<ndim; i++) {
					s[i] += v[i];
					cls[class][i] += v[i];
					cp[i][i] += v[i]*v[i];
					clcp[class][i][i] += v[i]*v[i];
				}
			}
		}

		for (i=0; i<ndim; i++) {
			/* compute and print per class means, standard deviations */
			classmeansd(0,i,nclass,nsamp,clnsamp,clm,clsd,m,s,cls,clcp);
		}
		fclose(outlf);
		infile = fopen(VCVTEMP,"r");
	} else
		infile = stdin;

	/* clear histograms */
	nsamp = 0;
	for (class=0; class<nclass; class++)
		for (i=0; i<ndim; i++)
			for (j=0; j<NBUCKETS; j++)
				bucket[class][i][j] = 0;

	/* get data vector, check class, check for derived objects */
	nsamp = 0;
	while (fscanf(infile,"%d",&class) > 0) {
		class--;
		if (DERIV != 0 && class >= 100)
			class -= 100;

		/* increment histogram bucket */
		if (class >= 0 && class < nclass) {
			nsamp++;
			for (i=0; i<ndim; i++) {
				fscanf(infile,"%e",&w);
				w -= clm[class][i];
				w /= clsd[class][i];
				w *= NBUCKETS;
				w /= 6.0;
				w += NBUCKETS/2;
				j = w;
				if (j < 0)
					j = 0;
				if (j >= NBUCKETS)
					j = NBUCKETS - 1;
				++bucket[class][i][j];
		}
		INEND:;
		}
	}
	printf("\n%d samples\n",nsamp);
	printf("class: feature: histogram\n");
	for (class=0; class<nclass; class++) {
		for (i=0; i<ndim; i++) {
			printf("%2d: %2d: ",class,i);
			for (j=0; j<NBUCKETS; j++)
				printf("%5d",bucket[class][i][j]);
			printf("\n");
		}
	}
}

/* compute per class means, standard deviations */
classmeansd(STATS,i,nclass,nsamp,clnsamp,clm,clsd,m,s,cls,clcp)
int *clnsamp;
double clm[MAXCLASS][MAXDIM], clsd[MAXCLASS][MAXDIM];
double m[MAXDIM], s[MAXDIM], cls[MAXCLASS][MAXDIM];
double clcp[MAXCLASS][MAXDIM][MAXDIM];
{
	int class;
	if (STATS > 1)
		printf("\f\nclass   n      mean     standard deviation       for feature %d\n",i);
	m[i] = s[i]/nsamp;
	for (class=0; class<nclass; class++) {
		if (clnsamp[class] > 0)
			clm[class][i] = cls[class][i]/clnsamp[class];
		else
			clm[class][i] = 0.0;
		if (clnsamp[class] > 1)
			clsd[class][i] = sqrt((clcp[class][i][i] - clnsamp[class]*clm[class][i]*clm[class][i])/(clnsamp[class]-1));
		else
			clsd[class][i] = 0.0;
		if (STATS > 1) {
			printf("%2d   %4d  ",class+1,clnsamp[class]);
			printf(" %8.2f",clm[class][i]);
			printf(" %8.2f\n",clsd[class][i]);
		}
	}
}