/* trackOverlap - Correlate a track with a series of tracks specified in specFile. */

/* Copyright (C) 2011 The Regents of the University of California 
 * See kent/LICENSE or http://genome.ucsc.edu/license/ for licensing information. */
#include "common.h"
#include "linefile.h"
#include "hash.h"
#include "cheapcgi.h"
#include "hdb.h"
#include "jksql.h"
#include "bits.h"
#include "featureBits.h"


void usage()
/* Explain usage and exit. */
{
errAbort(
  "trackOverlap- Overlap how much of a track is overlapped by\n"
  "   other tracks and vice versa. This is done by correlating\n"
  "   series of bitmap projections (i.e. featureBits multiple times).\n"
  "usage:\n"
  "   trackOverlap database chromosome track listFile\n"
  "\n"
  "Where listFile is a file with one featureBits specification per\n"
  "line. For example: 'intronEst:exon:10', see featureBits for\n"
  "more examples of syntax.\n"
//  "Use 'all' in the chromosome to cover the whole genome\n"
  );
}

void explainSome(char *database, Bits *homo, Bits *once, Bits *bits, char *chrom, int chromSize, 
	struct sqlConnection *conn, char *trackSpec, char *homologyTrack)
/* Explain some of homology. */
{
int trackSize = 0, homoSize = 0, andSize = 0, cumSize = 0, newSize = 0;

homoSize = bitCountRange(homo, 0, chromSize);
bitClear(bits, chromSize);
if (trackSpec != NULL)
    {
    fbOrTableBits(database, bits, trackSpec, chrom, chromSize, conn);
    trackSize = bitCountRange(bits, 0, chromSize);
    bitAnd(bits, homo, chromSize);
    andSize = bitCountRange(bits, 0, chromSize);
    bitAnd(bits, once, chromSize);
    newSize = bitCountRange(bits, 0, chromSize);
    bitNot(bits, chromSize);
    bitAnd(once, bits, chromSize);
    cumSize = homoSize - bitCountRange(once, 0, chromSize);
    }
else
    {
    trackSpec = homologyTrack;
    trackSize = andSize = homoSize;
    cumSize = newSize = 0;
    }

printf("%-21s %8d %8d %5.2f%% %6.2f%% %6.2f%% %5.2f%% %5.2f%%\n",
	trackSpec, trackSize, andSize, 
	100.0*trackSize/chromSize,
	100.0*andSize/trackSize,
	100.0*andSize/homoSize,
	100.0*newSize/homoSize,
	100.0*cumSize/homoSize);
}

void trackOverlap(char *database, char *chrom, char *homologyTrack, char *specFile)
/* trackOverlap - Correlate a track with a series of tracks specified in specFile. */
{
struct lineFile *lf = NULL;
char *line = NULL;
struct sqlConnection *conn;
int chromSize;
Bits *homo = NULL;
Bits *bits = NULL;
Bits *once = NULL;

lf = lineFileOpen(specFile, TRUE);
conn = hAllocConn(database);
chromSize = hChromSize(database, chrom);

homo = bitAlloc(chromSize);
bits = bitAlloc(chromSize);
once = bitAlloc(chromSize);

/* Get homology bitmap and set once mask to be the same. */
fbOrTableBits(database, homo, homologyTrack, chrom, chromSize, conn);
bitOr(once, homo, chromSize);

/* printHeader */
printf("%-21s %8s %8s %5s %6s  %6s %5s  %5s \n",
	"Track Specification", "track", "overlap", 
	"track", "cov",  "track", "new",  "cum");
printf("%-21s %8s %8s %5s %6s  %6s %5s  %5s \n",
	"", "size", "size", 
	"geno", "track", "cov",  "cov",  "cov");
printf("-----------------------------------------------------------------------------\n");

/* Whittle awway at homology... */
explainSome(database, homo, once, bits, chrom, chromSize, conn, NULL, homologyTrack);
while(lineFileNextReal(lf, &line))
    {
    explainSome(database, homo, once, bits, chrom, chromSize, conn, line, NULL);
    }
lineFileClose(&lf);
hFreeConn(&conn);
}



int main(int argc, char *argv[])
/* Process command line. */
{
cgiSpoof(&argc, argv);
if (argc != 5)
    usage();
trackOverlap(argv[1], argv[2], argv[3], argv[4]);
return 0;
}