677b513b231891e30d32024de7ded7766ef72058
braney
  Wed May 24 14:09:42 2017 -0700
add mathWig track type which is an arithmetic combination of bigWi.

diff --git src/hg/hgTracks/bigWigTrack.c src/hg/hgTracks/bigWigTrack.c
index a6298e7..8c9ae03 100644
--- src/hg/hgTracks/bigWigTrack.c
+++ src/hg/hgTracks/bigWigTrack.c
@@ -4,30 +4,31 @@
  * See README in this or parent directory for licensing information. */
 
 #include "common.h"
 #include "hash.h"
 #include "linefile.h"
 #include "jksql.h"
 #include "hdb.h"
 #include "hgTracks.h"
 #include "localmem.h"
 #include "wigCommon.h"
 #include "bbiFile.h"
 #include "bigWig.h"
 #include "errCatch.h"
 #include "container.h"
 #include "bigWarn.h"
+#include "mathWig.h"
 
 struct preDrawContainer *bigWigLoadPreDraw(struct track *tg, int seqStart, int seqEnd, int width)
 /* Do bits that load the predraw buffer tg->preDrawContainer. */
 {
 /* Just need to do this once... */
 if (tg->preDrawContainer != NULL)
     return tg->preDrawContainer;
 
 struct preDrawContainer *preDrawList = NULL;
 /* protect against temporary network error */
 struct errCatch *errCatch = errCatchNew();
 if (errCatchStart(errCatch))
     {
     /* Allocate predraw area. */
 
@@ -113,30 +114,136 @@
     {
     struct bbiFile *bbiFile = bigWigFileOpen(fileName);
     slAddHead(&tg->bbiFile, bbiFile);
     }
 errCatchEnd(errCatch);
 if (errCatch->gotError)
     {
     tg->networkErrMsg = cloneString(errCatch->message->string);
     tg->drawItems = bigDrawWarning;
     if (!sameOk(parentContainerType(tg), "multiWig"))
 	tg->totalHeight = bigWarnTotalHeight;
     }
 errCatchFree(&errCatch);
 }
 
+static void dataToPixelsUp(double *data, struct preDrawContainer *pre)
+/* Up sample data into pixels. */
+{
+int preDrawZero = pre->preDrawZero;
+unsigned size = winEnd - winStart;
+double dataPerPixel = size / (double) insideWidth;
+int pixel;
+
+for (pixel=0; pixel<insideWidth; ++pixel)
+    {
+    struct preDrawElement *pe = &pre->preDraw[pixel + preDrawZero];
+    unsigned index = pixel * dataPerPixel;
+    pe->count = 1;
+    pe->min = data[index];
+    pe->max = data[index];
+    pe->sumData = data[index] ;
+    pe->sumSquares = data[index] * data[index];
+    }
+}
+
+static void dataToPixelsDown(double *data, struct preDrawContainer *pre)
+/* Down sample data into pixels. */
+{
+int preDrawZero = pre->preDrawZero;
+unsigned size = winEnd - winStart;
+double dataPerPixel = size / (double) insideWidth;
+int pixel;
+
+for (pixel=0; pixel<insideWidth; ++pixel)
+    {
+    struct preDrawElement *pe = &pre->preDraw[pixel + preDrawZero];
+    double startReal = pixel * dataPerPixel;
+    double endReal = (pixel + 1) * dataPerPixel;
+    unsigned startUns = startReal;
+    unsigned endUns = endReal;
+    double realCount, realSum, realSumSquares, max, min;
+
+    realCount = realSum = realSumSquares = 0.0;
+    max = min = data[startUns];
+
+    assert(startUns != endUns);
+    unsigned ceilUns = ceil(startReal);
+
+    if (ceilUns != startUns)
+	{
+	/* need a fraction of the first count */
+	double frac = (double)ceilUns - startReal;
+	realCount = frac;
+	realSum = frac * data[startUns];
+	realSumSquares = realSum * realSum;
+	startUns++;
+	}
+
+    // add in all the data that are totally in this pixel
+    for(; startUns < endUns; startUns++)
+	{
+	realCount += 1.0;
+	realSum += data[startUns];
+	realSumSquares += data[startUns] * data[startUns];
+	if (max < data[startUns])
+	    max = data[startUns];
+	if (min > data[startUns])
+	    min = data[startUns];
+	}
+
+    // add any fraction of the count that's only partially in this pixel
+    double lastFrac = endReal - endUns;
+    double lastSum = lastFrac * data[endUns];
+    if ((lastFrac > 0.0) && (endUns < size))
+	{
+	if (max < data[endUns])
+	    max = data[endUns];
+	if (min > data[endUns])
+	    min = data[endUns];
+	realCount += lastFrac;
+	realSum += lastSum;
+	realSumSquares += lastSum * lastSum;
+	}
+
+    pe->count = normalizeCount(pe, realCount, min, max, realSum, realSumSquares);
+    }
+}
+
+static void dataToPixels(double *data, struct preDrawContainer *pre)
+/* Sample data into pixels. */
+{
+unsigned size = winEnd - winStart;
+double dataPerPixel = size / (double) insideWidth;
+
+if (dataPerPixel <= 1.0)
+    dataToPixelsUp(data, pre);
+else
+    dataToPixelsDown(data, pre);
+}
+
+static void mathWigLoadItems(struct track *tg)
+/* Fill up tg->items with bedGraphItems derived from a bigWig file */
+{
+char *equation = cloneString(trackDbSetting(tg->tdb, "mathDataUrl"));
+
+struct preDrawContainer *pre = tg->preDrawContainer = initPreDrawContainer(insideWidth);
+double *data = mathWigGetValues(equation, chromName, winStart, winEnd);
+dataToPixels(data, pre);
+
+free(data);
+}
 
 static void bigWigLoadItems(struct track *tg)
 /* Fill up tg->items with bedGraphItems derived from a bigWig file */
 {
 char *extTableString = trackDbSetting(tg->tdb, "extTable");
 
 if (tg->bbiFile == NULL)
     {
     /* Figure out bigWig file name. */
     if (tg->parallelLoading) // do not use mysql during parallel-fetch
 	{
 	char *fileName = cloneString(trackDbSetting(tg->tdb, "bigDataUrl"));
 	bigWigOpenCatch(tg, fileName);
 	}
     else
@@ -144,26 +251,34 @@
 	struct sqlConnection *conn = hAllocConnTrack(database, tg->tdb);
 	char *fileName = bbiNameFromSettingOrTable(tg->tdb, conn, tg->table);
 	bigWigOpenCatch(tg, fileName);
 	// if there's an extra table, read this one in too
 	if (extTableString != NULL)
 	    {
 	    fileName = bbiNameFromSettingOrTable(tg->tdb, conn, extTableString);
 	    bigWigOpenCatch(tg, fileName);
 	    }
 	hFreeConn(&conn);
 	}
     }
 bigWigLoadPreDraw(tg, winStart, winEnd, insideWidth);
 }
 
+void mathWigMethods(struct track *track, struct trackDb *tdb, 
+	int wordCount, char *words[])
+/* Set up bigWig methods. */
+{
+bigWigMethods(track, tdb, wordCount, words);
+track->loadItems = mathWigLoadItems;
+}
+
 void bigWigMethods(struct track *track, struct trackDb *tdb, 
 	int wordCount, char *words[])
 /* Set up bigWig methods. */
 {
 bedGraphMethods(track, tdb, wordCount, words);
 track->loadItems = bigWigLoadItems;
 track->preDrawItems = bigWigPreDrawItems;
 track->preDrawMultiRegion = wigMultiRegionGraphLimits;
 track->drawItems = bigWigDrawItems;
 track->loadPreDraw = bigWigLoadPreDraw;
 }