/*
 * dcsummary.c		Jim Piper		02-02-88
 *
 * Summarise results of dicentric classifier
 *
 * MODIFICATIONS
 * -------------
 * 12-04-88	jimp	Exclude overlapped FP frpm FP rate computation.
 * 28-03-88	jimp	Find overlap contribution to FN and FP
 * 25-02-88	jimp	Divide input data into buckets; compute FP and FN
 *			variance between buckets.
 * 24-02-88	jimp	Do not permit dicentrics smaller than cell median area
 */
#include <stdio.h>
#include <dicen.h>

main(argc,argv)
char ** argv;
{
	static struct dcresults r,tr = {0};
	static int fpd[5] = {0};
	int bucketsize = 90;
	int i,nbuckets,count;
	FILE *pf = stdin;
	int FPa,FNa,TPa,NTPa;
	int BS[100],TTP[100];
	double FP[100],FN[100],NFN,FNsigstat,FNsigsyst,sqrt();
	double TNFN,TNFP,NFP,FPsigstat,FPsigsyst;
	double FPm,FPsd,FNm,FNsd;
	nbuckets = count = FPa = TPa = FNa = NTPa = 0;
	NFN = FNsigsyst = FNsigstat = 0.0;
	while (fscanf(pf,"%d %d %d [%d %d]",&r.TP,&r.FP,&r.FN,&r.UTP,&r.UFN) != EOF) {
		fscanf(pf," (%1d %1d %1d) (%1d %1d %1d %1d)",&r.compoFP,&r.twistedFP,&r.ovFP,&r.acroFN,&r.onecandFN,&r.areaFN,&r.ovFN);
		fscanf(pf," <%2d %2d>\n",&r.nchrom,&r.nwithcands);
		tr.TP += r.TP;
		tr.FP += r.FP;
		tr.FN += r.FN;
		tr.UTP += r.UTP;
		tr.UFN += r.UFN;
		tr.compoFP += r.compoFP;
		tr.twistedFP += r.twistedFP;
		tr.acroFN += r.acroFN;
		tr.onecandFN += r.onecandFN;
		tr.areaFN += r.areaFN;
		tr.ovFN += r.ovFN;
		tr.ovFP += r.ovFP;
		tr.ncells++;
		/*
		 * fill the bucket
		 */
		FPa += r.FP - r.ovFP;
		FNa += r.FN - r.acroFN - r.UFN;
		TPa += r.TP - r.UTP;
		if (++count == bucketsize) {
			TPa += FNa;
			NTPa += TPa;
			NFN += FNa;
			NFP += FPa;
			TTP[nbuckets] = TPa;
			if (TPa > 0)
				FN[nbuckets] = (double) FNa / TPa;
			FP[nbuckets] = (double) FPa / count;
			BS[nbuckets] = bucketsize;
			nbuckets++;
			count = FPa = FNa = TPa = 0;
		}
		if (r.FP > 4)
			r.FP = 4;
		fpd[r.FP]++;
	}
	printf("%d cells analysed\n",tr.ncells);
	printf("TP %3d (unsure    %2d)\n",tr.TP,tr.UTP);
	printf("FN %3d (too small %2d, acrocentric  %2d, unsure %2d, overlapped %2d, <2 candidates %2d)\n",tr.FN,tr.areaFN,tr.acroFN,tr.UFN,tr.ovFN,tr.onecandFN);
	printf("FP %3d (overlap   %2d, twisted arms %2d, missed compo %2d)\n",tr.FP,tr.ovFP,tr.twistedFP,tr.compoFP);
	printf("       (FP distribution 0:%d 1:%d 2:%d 3:%d >=4:%d)\n",
		fpd[0],fpd[1],fpd[2],fpd[3],fpd[4]);
	/*
	 * complete the last bucket
	 */
	if (count > 0) {
		TPa += FNa;
		NTPa += TPa;
		NFN += FNa;
		NFP += FPa;
		TTP[nbuckets] = TPa;
		if (TPa > 0)
			FN[nbuckets] = (double) FNa / TPa;
		FP[nbuckets] = (double) FPa / count;
		BS[nbuckets] = count;
		nbuckets++;
	}
#ifdef BUCKETS
	printf("#cells:#dicen     ");
	for (i=0; i<nbuckets; i++)
		printf("%3d:%2d ",BS[i],TTP[i]);
	printf("\n");
	printf("FN                ");
	for (i=0; i<nbuckets; i++)
		printf("%6.3f ",FN[i]);
	printf("\n");
	printf("FP per cell       ");
	for (i=0; i<nbuckets; i++)
		printf("%6.3f ",FP[i]);
	printf("\n");
#endif BUCKETS

	NFN /= NTPa;
	TNFN = (double) tr.FN / (tr.FN + tr.TP);
	TNFP = (double) tr.FP / tr.ncells;
	for (i=0; i<nbuckets; i++) {
		FNsigstat += TTP[i] * FN[i] * (1.0 - FN[i]);
		FNsigsyst += TTP[i] * (FN[i] - NFN) * (FN[i] - NFN);
	}
	FNsigstat = sqrt(FNsigstat) / NTPa;
	FNsigsyst = sqrt(FNsigsyst) / NTPa;
	printf("FN: rate %5.3f (%5.3f) sigma(stat) %6.3f sigma(syst) %6.3f\n",TNFN,NFN,FNsigstat,FNsigsyst);

	NFP /= tr.ncells;
	for (i=0; i<nbuckets; i++) {
		FPsigstat += BS[i] * FP[i] * (1.0 - FP[i]);
		FPsigsyst += BS[i] * (FP[i] - NFP) * (FP[i] - NFP);
	}
	FPsigstat = sqrt(FPsigstat) / tr.ncells;
	FPsigsyst = sqrt(FPsigsyst) / tr.ncells;
	printf("FP: rate %5.3f (%5.3f) sigma(stat) %6.3f sigma(syst) %6.3f\n",TNFP,NFP,FPsigstat,FPsigsyst);
	printf("%5.3f %5.3f\n",TNFN,TNFP);	/* to allow pick-up / stuff */
}