/*
 * xpvcvclass.c
 *
 * classify feature files from x-probed chromosomes
 */
#include <stdio.h>
#include <wstruct.h>
#include <chromanal.h>
#include <otypes.h>
#include "xprobe.h"
#include <vcv.h>


/* logs of prior probabilities */
static float priorp[MAXCLASS] = {
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
1.0, 1.0, 1.0, 1.0};
static int nclass,ndim,nsfv,sfv[MAXDIM];

extern int POOL;

main(argc,argv)
char **argv;
{
	FILE *fclass;
	int i,j,trueclass,class,vclass,ivclass,xrank;
	int TPM,FPM,FNM,TPF,FPF,FNF,TNF,TNM;
	struct xfeats xf;
	FILE *in;
	
	TPM=FPM=FNM=TPF=FPF=FNF=TNF=TNM=0;
	/* pooled classifier ? */
	if (argc > 1 && argv[1][0] == '-' && argv[1][1] == 'p') {
		POOL = 1;
		argc--;
		argv++;
	}
	/* read class mean vectors and variance-covariance matrix */
	fclass = fopen(argv[1],"r");
	if (fclass == NULL) {
		fprintf(stderr,"Cannot open classifier parameter file\n");
		exit(2);
	}
	readclassparams(fclass,&nclass,&ndim,sfv,&nsfv);

	if (argc > 2) {
		in = fopen(argv[2],"r");
		if (in == NULL) {
			fprintf(stderr,"Cannot open file %s\n",argv[2]);
			exit(2);
		}
	}
	else
		in = stdin;
	printf("rm -f TPMs TPFs FPMs FPFs FNMs FNFs\n");

	/*
	 * read feature data file
	 */
	while (readxpboxfeat(in,&xf) != EOF) {
		/*
		 * reset rank - data from new cell
		 */
		if (xf.xrank == 1)
			xrank = 0;
		vcvclassxfeat(&xf);
		if (xf.vcvclass == 1) {
#ifdef LIKLTEST
printf("True %d Rank %d likl %E  : ",xf.Cbtype,xrank+1,xf.vcvscore);
#endif LIKLTEST
			xf.xrank = ++xrank;	/* revise xf.xrank */
			xf.classifmale = xf.classiffemale = 1;
			if (xrank > 2)
				xf.classiffemale = 2;
			if (xrank > 1)
				xf.classifmale = 2;
		}
		else
			xf.classifmale = xf.classiffemale = 2;
		/*
		 * Accumulate TP, FP, FN
		 */
		if (xf.Cbtype == 1) {
			if (xf.classifmale == 1) {
				TPM++;
				printf("nsel %d <%s >>TPMs\n",xf.number,xf.fname);
			}
			else {
				FNM++;
				printf("nsel %d <%s >>FNMs\n",xf.number,xf.fname);
			}
			if (xf.classiffemale == 1) {
				TPF++;
				printf("nsel %d <%s >>TPFs\n",xf.number,xf.fname);
			}
			else {
				FNF++;
				printf("nsel %d <%s >>FNFs\n",xf.number,xf.fname);
			}
		}
		else {
			if (xf.classifmale == 1) {
				FPM++;
				printf("nsel %d <%s >>FPMs\n",xf.number,xf.fname);
			}
			else
				TNM++;
			if (xf.classiffemale == 1) {
				FPF++;
				printf("nsel %d <%s >>FPFs\n",xf.number,xf.fname);
			}
			else
				TNF++;
		}
/*
		printf("%1d %1d %d\n",xf.Cbtype,xf.classiffemale,xf.classifmale);
*/
	}
	printf("If male:   TP%3d FN%3d FP%3d TN%3d\n",TPM,FNM,FPM,TNM);
	printf("If female: TP%3d FN%3d FP%3d TN%3d\n",TPF,FNF,FPF,TNF);
}

vcvclassxfeat(xf)
register struct xfeats *xf;
{
	float fvec[MAXDIM],svec[MAXDIM];
	struct chromcl chrcl;
	int c, i, class, vclass;

	getxfeat(fvec,xf);
	for (i=0; i<nsfv; i++)
		svec[i] = fvec[sfv[i]];
	/* likelihood classifier */
	xf->vcvclass = classdist(svec,priorp,chrcl.cllik,nclass,nsfv);
	xf->vcvscore = chrcl.cllik[0].lval;
}