44ccfacbe3a3d4b300f80d48651c77837a4b571e galt Tue Apr 26 11:12:02 2022 -0700 SQL INJECTION Prevention Version 2 - this improves our methods by making subclauses of SQL that get passed around be both easy and correct to use. The way that was achieved was by getting rid of the obscure and not well used functions sqlSafefFrag and sqlDyStringPrintfFrag and replacing them with the plain versions of those functions, since these are not needed anymore. The new version checks for NOSQLINJ in unquoted %-s which is used to include SQL clauses, and will give an error the NOSQLINJ clause is not present, and this will automatically require the correct behavior by developers. sqlDyStringPrint is a very useful function, however because it was not enforced, users could use various other dyString functions and they operated without any awareness or checking for SQL correct use. Now those dyString functions are prohibited and it will produce an error if you try to use a dyString function on a SQL string, which is simply detected by the presence of the NOSQLINJ prefix. diff --git src/hg/lib/wigDataStream.c src/hg/lib/wigDataStream.c index ee7a664..66f45b8 100644 --- src/hg/lib/wigDataStream.c +++ src/hg/lib/wigDataStream.c @@ -1,2294 +1,2300 @@ /* wigDataStream - an object to access wiggle data values from * either a DB access or from a .wig text file (==custom track) */ /* Copyright (C) 2014 The Regents of the University of California * See kent/LICENSE or http://genome.ucsc.edu/license/ for licensing information. */ #include "common.h" #include "memalloc.h" #include "wiggle.h" #include "portable.h" #include "hgColors.h" #include "obscure.h" #include "customTrack.h" #include "udc.h" /* size of memory to work with during wiggle data fetch operations, * ~4 Gb is good enough for about 750,000,000 bases */ /* raising this number has a big effect and slows down some tracks ? * namely, psl/bed custom tracks, this doesn't make any sense */ #define HUGE_LIMIT 2100000000 /* Routines that are not strictly part of the wigDataStream object, but they are used to do things with the object. **************************************************************************/ void wigStatsTableHeading(FILE * fh, boolean htmlOut) /* Print the single html (or text) table row for statistics column headings */ { if (htmlOut) { fprintf(fh,"<TR><TH> Chrom </TH><TH> Data <BR> start </TH>"); fprintf(fh,"<TH> Data <BR> end </TH>"); fprintf(fh,"<TH> # of Data <BR> values </TH><TH> " "Each data <BR> value spans <BR> # bases </TH>"); fprintf(fh,"<TH> Bases <BR> covered </TH><TH> Minimum </TH>"); fprintf(fh,"<TH> Maximum </TH><TH> Range </TH><TH> Mean </TH>"); fprintf(fh,"<TH> Variance </TH><TH> Standard <BR> deviation </TH></TR>\n"); } else { fprintf(fh,"# Chrom\tData\tData"); fprintf(fh,"\t# Data\tData"); fprintf(fh,"\tBases\tMinimum"); fprintf(fh,"\tMaximum\tRange\tMean"); fprintf(fh,"\tVariance Standard\n"); fprintf(fh,"#\tstart\tend"); fprintf(fh,"\tvalues\tspan"); fprintf(fh,"\tcovered\t"); fprintf(fh,"\t\t\t"); fprintf(fh,"\t\tdeviation\n"); } } void wigPrintDataConstraint(struct wiggleDataStream *wds, FILE * fh) /* output string to file handle fh indicating current data constraint */ { if (wds->useDataConstraint && wds->dataConstraint) { if (sameWord(wds->dataConstraint,"in range")) fprintf (fh, "#\tdata values in range [%g : %g)\n", wds->limit_0, wds->limit_1); else fprintf (fh, "#\tdata values %s %g\n", wds->dataConstraint, wds->limit_0); } } void statsPreamble(struct wiggleDataStream *wds, char *chrom, int winStart, int winEnd, unsigned span, unsigned long long valuesMatched, char *table2) { char num1Buf[64], num2Buf[64]; /* big enough for 2^64 (and then some) */ sprintLongWithCommas(num1Buf, BASE_1(winStart)); sprintLongWithCommas(num2Buf, winEnd); printf("<P><B> Position: </B> %s:%s-%s</P>\n", chrom, num1Buf, num2Buf ); sprintLongWithCommas(num1Buf, winEnd - winStart); printf("<P><B> Total Bases in view: </B> %s </P>\n", num1Buf); if (wds->useDataConstraint) { if (sameWord(wds->dataConstraint,"in range")) printf("<P><B> Filter: %g <= (data value) < %g </B></P>\n", wds->limit_0, wds->limit_1); else printf("<P><B> Filter: (data value %s %g) </B> </P>\n", wds->dataConstraint, wds->limit_0); } if (table2) printf("<P><B> Intersection with table: %s </B></P>\n", table2); if (table2 && (valuesMatched == 0)) { printf("<P><B> No data found in this intersection. </B></P>\n"); } else { if (valuesMatched == 0) { if ( (span * 3) > (winEnd - winStart)) { printf("<P><B> Viewpoint has too few bases to calculate statistics. </B></P>\n"); printf("<P><B> Zoom out to at least %d bases to see statistics. </B></P>\n", 3 * span); } else printf("<P><B> No data found in this region. </B></P>\n"); } else { unsigned printSpan = wds->stats->span; if ((span > 0) && (span < wds->stats->span)) printSpan = span; sprintLongWithCommas(num1Buf, wds->stats->count * printSpan); printf( "<P><B> Statistics on: </B> %s <B> bases </B> (%% %.4f coverage)</P>\n", num1Buf, 100.0*(wds->stats->count * printSpan)/(winEnd - winStart)); } } } void wigStatsHeader(struct wiggleDataStream *wds, FILE * fh, boolean htmlOut) /* begin wiggle stats table */ { if (htmlOut) { /* For some reason BORDER=1 does not work in our web.c nested table * scheme. So use web.c's trick of using an enclosing table * to provide a border. */ fprintf(fh,"<P><!--outer table is for border purposes-->" "\n" "<TABLE BGCOLOR=\"#" HG_COL_BORDER "\" BORDER=\"0\" CELLSPACING=\"0\" CELLPADDING=\"1\"><TR><TD>"); fprintf(fh, "<TABLE COLS=12 BORDER=1 BGCOLOR='#" HG_COL_INSIDE "' ALIGN=CENTER HSPACE=0><TR>"); if (wds->db) fprintf(fh, "<TH COLSPAN=6 ALIGN=LEFT> Database: %s </TH><TH COLSPAN=6 ALIGN=RIGHT> Table: %s </TH></TR>\n", wds->db, wds->tblName); if (wds->isFile) { if ( (stringIn("trash/ct_",wds->tblName)) || (stringIn("trash/hgtct_",wds->tblName))) fprintf(fh, "<TH COLSPAN=12 ALIGN=LEFT> custom track </TH></TR>\n" ); else fprintf(fh, "<TH COLSPAN=12 ALIGN=LEFT> from file %s </TH></TR>\n", wds->tblName); } wigStatsTableHeading(fh, htmlOut); } else { if (wds->db) fprintf(fh, "#\t Database: %s, Table: %s\n", wds->db, wds->tblName); if (wds->isFile) fprintf(fh, "#\t from file, Table: %s\n", wds->tblName); wigStatsTableHeading(fh, htmlOut); } } /* void wigStatsHeader() */ /* strictly object methods following ************************/ /* PRIVATE METHODS ************************************************/ static void addConstraint(struct wiggleDataStream *wds, char *left, char *op, char *right) { -struct dyString *constrain = dyStringNew(256); +struct dyString *constraint = dyStringNew(256); if (wds->sqlConstraint) - dyStringPrintf(constrain, "%s AND ", wds->sqlConstraint); + sqlDyStringPrintf(constraint, "%-s AND ", wds->sqlConstraint); -sqlDyStringPrintfFrag(constrain, "%s %-s \"%s\"", left, op, right); +sqlDyStringPrintf(constraint, "%s %-s '%s'", left, op, right); freeMem(wds->sqlConstraint); /* potentially previously existing */ -wds->sqlConstraint = cloneString(constrain->string); -dyStringFree(&constrain); +wds->sqlConstraint = dyStringCannibalize(&constraint); } /* *row[] is artifically one too big to allow for a potential bin * column when reading files that may have it. */ static boolean nextRow(struct wiggleDataStream *wds, char *row[], int maxRow) /* read next wig row from sql query or lineFile * FALSE return on no more data */ { int numCols; if (wds->isFile) { /* the file may have a bin column, detect automatically and eliminate */ numCols = lineFileChopNextTab(wds->lf, row, maxRow+1); if (numCols == maxRow+1) { int i; for (i = 1; i < (maxRow+1); ++i) row[i-1] = row[i]; } if (numCols < maxRow) return FALSE; verbose(VERBOSE_PER_VALUE_LEVEL, "#\tnumCols = %d, row[0]: %s, row[1]: %s, row[%d]: %s\n", numCols, row[0], row[1], maxRow-1, row[maxRow-1]); } else { int i; char **sqlRow; sqlRow = sqlNextRow(wds->sr); if (sqlRow == NULL) return FALSE; /* skip the bin column sqlRow[0] */ for (i=1; i <= maxRow; ++i) { row[i-1] = sqlRow[i]; } } return TRUE; } static void findWibFile(struct wiggleDataStream *wds, char *file) /* look for file in full pathname given, or in same directory */ { wds->wibFile = hReplaceGbdb(file); wds->wibFH = udcFileMayOpen(wds->wibFile, NULL); if (wds->wibFH == NULL) { char *baseName = strrchr(wds->wibFile, '/'); if (baseName) wds->wibFH = udcFileMayOpen(baseName+1, NULL); if ((NULL == baseName) || (wds->wibFH == NULL)) errAbort("findWibFile: failed to open %s", wds->wibFile); } } static void openWibFile(struct wiggleDataStream *wds, char *file) { if (wds->wibFile) { /* close and open only if different */ if (differentString(wds->wibFile,file)) { if (wds->wibFH > 0) udcFileClose(&wds->wibFH); freeMem(wds->wibFile); findWibFile(wds, file); } } else findWibFile(wds, file); verbose(VERBOSE_HIGHEST, "#\topened wib file: %s\n", wds->wibFile); } static void setCompareByte(struct wiggleDataStream *wds, double lower, double range) { if (wds->limit_0 < lower) wds->ucLowerLimit = 0; else if (wds->limit_0 > (lower+range)) wds->ucLowerLimit = MAX_WIG_VALUE; else wds->ucLowerLimit = MAX_WIG_VALUE * ((wds->limit_0 - lower)/range); if (wds->limit_1 < lower) wds->ucUpperLimit = 0; else if (wds->limit_1 > (lower+range)) wds->ucUpperLimit = MAX_WIG_VALUE + 1; /* +1 because in range is [) */ else wds->ucUpperLimit = MAX_WIG_VALUE * ((wds->limit_1 - lower)/range); verbose(VERBOSE_HIGHEST, "#\twigSetCompareByte: [%g : %g] becomes [%d : %d]\n", lower, lower+range, wds->ucLowerLimit, wds->ucUpperLimit); } static void resetStats(float *lowerLimit, float *upperLimit, double *sumData, double *sumSquares, unsigned long *statsCount, long int *chromStart, long int *chromEnd) { *lowerLimit = INFINITY; *upperLimit = -1.0 * INFINITY; *sumData = 0.0; *sumSquares = 0.0; *statsCount = 0; *chromStart = -1; *chromEnd = 0; } static void accumStats(struct wiggleDataStream *wds, float lowerLimit, float upperLimit, double sumData, double sumSquares, unsigned statsCount, long int chromStart, long int chromEnd) { if (statsCount > 0) { struct wiggleStats *ws; AllocVar(ws); ws->chrom = cloneString(wds->currentChrom); ws->chromStart = chromStart; ws->chromEnd = chromEnd; ws->span = wds->currentSpan; ws->count = statsCount; ws->lowerLimit = lowerLimit; ws->dataRange = upperLimit - lowerLimit; ws->mean = sumData / (double) statsCount; ws->variance = 0.0; ws->stddev = 0.0; if (statsCount > 1) { ws->variance = (sumSquares - ((sumData * sumData)/(double) statsCount)) / (double) (statsCount - 1); if (ws->variance > 0.0) ws->stddev = sqrt(ws->variance); } slAddHead(&wds->stats, ws); } } static struct bed *bedElement(char *chrom, unsigned start, unsigned end, unsigned lineCount) { struct bed *bed; char name[128]; AllocVar(bed); bed->chrom = cloneString(chrom); bed->chromStart = start; bed->chromEnd = end; safef(name, sizeof(name), "%s.%u", chrom, lineCount); bed->name = cloneString(name); return bed; } static void closeWibFile(struct wiggleDataStream *wds) /* if there is a Wib file open, close it */ { if (wds->wibFH > 0) udcFileClose(&wds->wibFH); wds->wibFH = (struct udcFile*)NULL; if (wds->wibFile) freez(&wds->wibFile); } static void closeWigConn(struct wiggleDataStream *wds) { lineFileClose(&wds->lf); closeWibFile(wds); /* closes only if it is open */ if (wds->conn) { sqlFreeResult(&wds->sr); hFreeConn(&wds->conn); } if (wds->sqlConstraint) freez(&wds->sqlConstraint); /* always reconstructed at open time */ } static void openWigConn(struct wiggleDataStream *wds) /* open connection to db or to a file, prepare SQL result for db */ { if (!wds->tblName) errAbort("openWigConn: table name missing. setDbTable before open."); if (wds->isFile) { struct dyString *fileName = dyStringNew(256); lineFileClose(&wds->lf); /* possibly a previous file */ /* don't add .wig if it is already there, or use whatever filename * was given */ if (fileExists(wds->tblName)) dyStringPrintf(fileName, "%s", wds->tblName); else dyStringPrintf(fileName, "%s.wig", wds->tblName); wds->lf = lineFileOpen(fileName->string, TRUE); dyStringFree(&fileName); } else { struct dyString *query = dyStringNew(256); sqlDyStringPrintf(query, "select * from %s", wds->tblName); + char op[1+NOSQLINJ_SIZE+1]; if (wds->chrName) - addConstraint(wds, "chrom", "=", wds->chrName); + { + sqlSafef(op, sizeof op, "="); + addConstraint(wds, "chrom", op, wds->chrName); + } if (wds->winEnd) { char limits[256]; safef(limits, ArraySize(limits), "%d", wds->winEnd ); - addConstraint(wds, "chromStart", "<", limits); + sqlSafef(op, sizeof op, "<"); + addConstraint(wds, "chromStart", op, limits); safef(limits, ArraySize(limits), "%d", wds->winStart ); - addConstraint(wds, "chromEnd", ">", limits); + sqlSafef(op, sizeof op, ">"); + addConstraint(wds, "chromEnd", op, limits); } if (wds->spanLimit) { struct dyString *dyTmp = dyStringNew(256); dyStringPrintf(dyTmp, "%u", wds->spanLimit); - addConstraint(wds, "span", "=", dyTmp->string); + sqlSafef(op, sizeof op, "="); + addConstraint(wds, "span", op, dyTmp->string); dyStringFree(&dyTmp); } if (wds->sqlConstraint) { - dyStringPrintf(query, " where "); + sqlDyStringPrintf(query, " where "); if (wds->winEnd) { hAddBinToQuery(wds->winStart, wds->winEnd, query); } - dyStringPrintf(query, " (%s)", + sqlDyStringPrintf(query, " (%-s)", wds->sqlConstraint); } - dyStringPrintf(query, " order by "); + sqlDyStringPrintf(query, " order by "); if (!wds->chrName) - dyStringPrintf(query, " chrom ASC,"); + sqlDyStringPrintf(query, " chrom ASC,"); if (!wds->spanLimit) - dyStringPrintf(query, " span ASC,"); - dyStringPrintf(query, " chromStart ASC"); + sqlDyStringPrintf(query, " span ASC,"); + sqlDyStringPrintf(query, " chromStart ASC"); verbose(VERBOSE_SQL_ROW_LEVEL, "#\t%s\n", query->string); if (!wds->conn) wds->conn = hAllocConn(wds->db); wds->sr = sqlGetResult(wds->conn,query->string); } } static void setDbTable(struct wiggleDataStream *wds, char * db, char *table) /* sets the DB and table name, determines if table is a file or not */ { struct dyString *fileName = dyStringNew(256); if (!table) errAbort("setDbTable: table specification missing"); /* Check to see if there is a .wig file, or whatever name was given */ if (fileExists(table)) dyStringPrintf(fileName, "%s", table); else dyStringPrintf(fileName, "%s.wig", table); /* file present ignores db specification */ if ((db == (char *)NULL) && fileExists(fileName->string)) wds->isFile = TRUE; else { if (db) /* db can be NULL */ wds->isFile = FALSE; /* assume it will be in the db */ else errAbort("setDbTable: db is NULL and can not find file %s", fileName->string); } dyStringFree(&fileName); freeMem(wds->db); /* potentially previously existing */ if (!wds->isFile && db) wds->db = cloneString(db); freeMem(wds->tblName); /* potentially previously existing */ wds->tblName = cloneString(table); } static char *dateTimeStamp() { time_t now = time(NULL); char *dateTime; dateTime = sqlUnixTimeToDate(&now,TRUE); /* TRUE == gmTime */ return dateTime; } static void outputIdentification(struct wiggleDataStream *wds, FILE *fh) { char *dateStamp = dateTimeStamp(); if (wds->db) fprintf (fh, "#\tdb: '%s', track: '%s', output date: %s UTC\n", wds->db, wds->tblName, dateStamp); if (wds->isFile) fprintf (fh, "#\tfrom file input, output date: %s UTC\n", dateStamp); freeMem(dateStamp); } static void showResolution(double resolution, FILE *fh) { if (resolution > 0.0) { fprintf (fh, "#\tThis data has been compressed with a minor " "loss in resolution.\n" ); fprintf (fh, "#\t(Worst case: %g) The original source data\n", resolution); fprintf (fh, "#\t(before querying and compression) is available at \n" "#\t\thttp://hgdownload.soe.ucsc.edu/downloads.html\n"); } } static void showResolutionNoDownloads(double resolution, FILE *fh) { if (resolution > 0.0) { fprintf (fh, "#\tThis data has been compressed with a minor " "loss in resolution.\n" ); fprintf (fh, "#\t(Worst case: %g)\n", resolution); } } static void showConstraints(struct wiggleDataStream *wds, FILE *fh) { if (wds->chrName) fprintf (fh, "#\tchrom specified: %s\n", wds->chrName); if (wds->spanLimit) fprintf (fh, "#\tspan specified: %u\n", wds->spanLimit); if (wds->winEnd) fprintf (fh, "#\tposition specified: %d-%d\n", BASE_1(wds->winStart), wds->winEnd); if (wds->bedConstrained && !wds->chrName) fprintf (fh, "#\tconstrained by chr names and coordinates in bed list\n"); else if (wds->bedConstrained) fprintf (fh, "#\tconstrained by coordinates in bed list\n"); wigPrintDataConstraint(wds, fh); } static int arrayDataCmp(const void *va, const void *vb) /* Compare to sort based on chrom,winStart */ { const struct wiggleArray *a = *((struct wiggleArray **)va); const struct wiggleArray *b = *((struct wiggleArray **)vb); int dif; dif = chrStrippedCmp(a->chrom, b->chrom); if (dif == 0) dif = a->winStart - b->winStart; return dif; } static int asciiDataCmp(const void *va, const void *vb) /* Compare to sort based on chrom,span,chromStart. */ { const struct wigAsciiData *a = *((struct wigAsciiData **)va); const struct wigAsciiData *b = *((struct wigAsciiData **)vb); int dif; dif = chrStrippedCmp(a->chrom, b->chrom); if (dif == 0) { dif = a->span - b->span; if (dif == 0) dif = a->data->chromStart - b->data->chromStart; } return dif; } static int statsDataCmp(const void *va, const void *vb) /* Compare to sort based on chrom,span,chromStart. */ { const struct wiggleStats *a = *((struct wiggleStats **)va); const struct wiggleStats *b = *((struct wiggleStats **)vb); int dif; dif = chrStrippedCmp(a->chrom, b->chrom); if (dif == 0) { dif = a->span - b->span; if (dif == 0) dif = a->chromStart - b->chromStart; } return dif; } /* currently this dataArrayOut is private, but it may become public */ static void dataArrayOut(struct wiggleDataStream *wds, char *fileName, boolean rawDataOut, boolean sort) /* print to fileName ascii data values from the first list element * of the data "array" results */ { FILE * fh; fh = mustOpen(fileName, "w"); if (wds->array) { unsigned pointCount; unsigned chromPosition; unsigned inx; float *dataPtr; if (sort) slSort(&wds->array, arrayDataCmp); dataPtr = wds->array->data; /* user visible coordinates are wds->offset relative/based */ chromPosition = wds->array->winStart + wds->offset; pointCount = wds->array->winEnd - wds->array->winStart; for (inx = 0; inx < pointCount; ++inx) { if (!isnan(*dataPtr)) { if (rawDataOut) fprintf (fh, "%g\n", *dataPtr); else fprintf (fh, "%u\t%g\n", chromPosition, *dataPtr); } ++dataPtr; ++chromPosition; } } else { if (!rawDataOut) { showConstraints(wds, fh); fprintf(fh, "#\tdataArray: no data points found\n"); } } carefulClose(&fh); } /* static void dataArrayOut() */ /* PUBLIC METHODS **************************************************/ static void setMaxOutput(struct wiggleDataStream *wds, unsigned long long maxOut) /* set the maximum # of values to return */ { wds->maxOutput = maxOut; verbose(VERBOSE_CHR_LEVEL,"#\tsetMaxOutput: %llu\n", wds->maxOutput); } static void setPositionConstraint(struct wiggleDataStream *wds, int winStart, int winEnd) /* both 0 means no constraint */ { if ((!wds->isFile) && wds->conn) { /* this should not happen */ errAbort("setPositionConstraint: not allowed after openWigConn()"); } /* keep them in proper order */ if (winStart > winEnd) { wds->winStart = winEnd; wds->winEnd = winStart; } else if ((winEnd > 0) && (winStart == winEnd)) errAbort( "setPositionConstraint: can not have winStart == winEnd (%d == %d)", winStart, winEnd); else { wds->winStart = winStart; wds->winEnd = winEnd; } verbose(VERBOSE_SQL_ROW_LEVEL,"#\tsetPosition: %d - %d\n", wds->winStart, wds->winEnd); } static void setChromConstraint(struct wiggleDataStream *wds, char *chr) { freeMem(wds->chrName); /* potentially previously existing */ wds->chrName = cloneString(chr); } static void setSpanConstraint(struct wiggleDataStream *wds, unsigned span) { wds->spanLimit = span; } static void freeConstraints(struct wiggleDataStream *wds) /* free the position, span, chrName and data constraints */ { wds->spanLimit = 0; wds->setPositionConstraint(wds, 0, 0); if (wds->chrName) freez(&wds->chrName); if (wds->dataConstraint) freez(&wds->dataConstraint); wds->wigCmpSwitch = wigNoOp_e; wds->limit_0 = wds->limit_1 = 0.0; wds->ucLowerLimit = wds->ucUpperLimit = 0; if (wds->sqlConstraint) freez(&wds->sqlConstraint); wds->useDataConstraint = FALSE; wds->bedConstrained = FALSE; } static void freeAscii(struct wiggleDataStream *wds) /* free the wiggleAsciiData list */ { if (wds->ascii) { struct wigAsciiData *el, *next; for (el = wds->ascii; el != NULL; el = next) { next = el->next; freeMem(el->chrom); freeMem(el->data); freeMem(el); } } wds->ascii = NULL; } static void freeBed(struct wiggleDataStream *wds) /* free the wiggle bed list */ { bedFreeList(&wds->bed); } static void freeArray(struct wiggleDataStream *wds) /* free the wiggleArray list */ { if (wds->array) { struct wiggleArray *wa; struct wiggleArray *next; for (wa = wds->array; wa; wa = next) { next = wa->next; freeMem(wa->chrom); freeMem(wa->data); freeMem(wa); } } wds->array = NULL; } static void freeStats(struct wiggleDataStream *wds) /* free the wiggleStats list */ { if (wds->stats) { struct wiggleStats *el, *next; for (el = wds->stats; el != NULL; el = next) { next = el->next; freeMem(el->chrom); freeMem(el); } } wds->stats = NULL; } enum wigCompare wigCompareFromString(char *string) /* Given a string return corresponding cmp code. */ { if (sameWord(string,"<")) return wigLessThan_e; else if (sameWord(string,"<=")) return wigLessEqual_e; else if (sameWord(string,"=")) return wigEqual_e; else if (sameWord(string,"!=")) return wigNotEqual_e; else if (sameWord(string,">=")) return wigGreaterEqual_e; else if (sameWord(string,">")) return wigGreaterThan_e; else if (sameWord(string,"in range")) return wigInRange_e; else { errAbort("wigSetCompareFunctions: unknown constraint: '%s'", string); return wigLessThan_e; } } /* the double comparison functions * used to check the wiggle SQL rows which are a bucket of values * between *lower and *upper. Therefore, the value to be checked * which is in wds->limit_0 (and wds->limit_1 in the case of * a range) needs to be compared to the bucket of values. If it * falls within the specified range, then it is considered to be in * that bucket. */ /* InRange means the SQL row begins before the limit_1 (lower<limit_1) * and the row ends after the limit_0 (upper>=limit_0) * i.e. there is at least some overlap of the range * This is a half open inquiry: [ limit_0 : limit_1 ) */ /* LessThan means: the row begins before the limit_0 (value<limit_0) * i.e. there are data values below the specified limit_0 */ /* LessEqual means: the row begins at or before the limit_0 (value<=limit_0) * i.e. there are data values at or below the specified limit_0 */ /* Equal means: similar to InRange, the test value limit_0 can be found * in the SQL row, i.e. lower <= limit_0 <= upper */ /* NotEqual means: the opposite of Equal, the test value limit_0 can not * be found in the SQL row, i.e. (limit_0 < lower) or (upper < limit_0) */ /* GreaterEqual means: the row ends at or after the limit_0 (limit_0<=upper) * i.e. there are data values at or above the specified limit_0 */ /* GreaterEqual means: the row ends after the limit_0 (limit_0<upper) * i.e. there are data values above the specified limit_0 */ /* the unsigned char comparison functions * Unlike the above, these are straighforward, just compare the * byte values */ static void wigSetCompareFunctions(struct wiggleDataStream *wds) { if (!wds->dataConstraint) return; wds->wigCmpSwitch = wigCompareFromString(wds->dataConstraint); } static void setDataConstraint(struct wiggleDataStream *wds, char *dataConstraint, double lowerLimit, double upperLimit) { freeMem(wds->dataConstraint); /* potentially previously existing */ wds->dataConstraint = cloneString(dataConstraint); if (differentWord(wds->dataConstraint, "in range")) { wds->limit_0 = lowerLimit; } else { if (lowerLimit < upperLimit) { wds->limit_0 = lowerLimit; wds->limit_1 = upperLimit; } else if (!(upperLimit < lowerLimit)) errAbort("wigSetDataConstraint: upper and lower limits are equal: %g == %g", lowerLimit, upperLimit); else { wds->limit_0 = upperLimit; wds->limit_1 = lowerLimit; } } wigSetCompareFunctions(wds); wds->useDataConstraint = TRUE; } boolean wigCompareValFilter(double val, enum wigCompare cmp, double ll, double ul) /* Return TRUE if val passes comparison filter. */ { boolean takeIt = FALSE; switch (cmp) { case wigNoOp_e: takeIt = TRUE; break; case wigInRange_e: takeIt = (val < ul) && (val >= ll); break; case wigLessThan_e: takeIt = val < ll; break; case wigLessEqual_e: takeIt = val <= ll; break; case wigEqual_e: takeIt = val == ll; break; case wigNotEqual_e: takeIt = val != ll; break; case wigGreaterEqual_e: takeIt = val >= ll; break; case wigGreaterThan_e: takeIt = val > ll; break; default: errAbort("wigCmpValFilter: illegal wigCmpSwitch ? %#x", cmp); break; } return takeIt; } static unsigned long long getData(struct wiggleDataStream *wds, char *db, char *table, int operations) /* getData - read and return wiggle data */ { char *row[WIGGLE_NUM_COLS+1]; /* potentially a bin column too */ unsigned long long rowCount = 0; unsigned long long validData = 0; unsigned long long valuesMatched = 0; unsigned long long noDataBytes = 0; unsigned long long bytesSkipped = 0; boolean doRawStats = FALSE; boolean doStats = FALSE; boolean doBed = FALSE; boolean doAscii = FALSE; boolean doDataArray = FALSE; boolean doNoOp = FALSE; boolean skipDataRead = FALSE; /* may need this later */ float lowerLimit = INFINITY; float upperLimit = -1.0 * INFINITY; double sumData = 0.0; double sumSquares = 0.0; unsigned long statsCount = 0; long int chromStart = -1; long int chromEnd = 0; boolean summaryOnly = TRUE; unsigned bedElStart = 0; unsigned bedElEnd = 0; unsigned bedElCount = 0; boolean firstSpanDone = FALSE; /* to prevent multiple bed lists */ float *dataArrayPtr = NULL; /* to access the data array values */ unsigned dataArrayPosition = 0; /* marches thru all from beginning to end */ struct wiggle *wiggle; /* one SQL data read results */ boolean maxReached = FALSE; doAscii = operations & wigFetchAscii; doDataArray = operations & wigFetchDataArray; doBed = operations & wigFetchBed; doRawStats = operations & wigFetchRawStats; doStats = (operations & wigFetchStats) || doRawStats; doNoOp = operations & wigFetchNoOp; /* for timing purposes, allow the wigFetchNoOp to go through */ if (doNoOp) doDataArray = doBed = doStats = doAscii = FALSE; if (! (doNoOp || doAscii || doDataArray || doBed || doStats) ) { verbose(VERBOSE_CHR_LEVEL, "getData: no type of data fetch requested ?\n"); return(valuesMatched); /* NOTHING ASKED FOR ? */ } /* not going to do this without a range specified and a chr name */ if (doDataArray && wds->winEnd && wds->chrName) { long startTime = 0; struct wiggleArray *wa; size_t inx, size; size_t winSize; unsigned long long hugeSize; AllocVar(wa); wa->chrom = cloneString(wds->chrName); wa->winStart = wds->winStart; wa->winEnd = wds->winEnd; winSize = wa->winEnd - wa->winStart; hugeSize = (unsigned long long)sizeof(float) * (unsigned long long)winSize; if (sizeof(size_t) == 4) { if (hugeSize > HUGE_LIMIT) { warn("ERROR: Can not perform requested data operation.<BR>"); errAbort("Not enough memory (requested: %llu bytes)", hugeSize); } } size = sizeof(float) * winSize; setMaxAlloc((size_t)HUGE_LIMIT); if (verboseLevel() >= VERBOSE_CHR_LEVEL) startTime = clock1000(); wa->data = (float *)needLargeMem(size); dataArrayPtr = wa->data; /* this init loop turns out to be quite fast, even for huge arrays */ for (inx = 0; inx < winSize; ++inx) *dataArrayPtr++ = NAN; dataArrayPtr = wa->data; dataArrayPosition = wa->winStart; slAddHead(&wds->array, wa); if (verboseLevel() >= VERBOSE_CHR_LEVEL) { long et = clock1000() - startTime; verbose(VERBOSE_CHR_LEVEL, "#\tgetData: created data array for %lld values (%llu b) in %ld ms\n", (unsigned long long)winSize, (unsigned long long)size, et); verbose(VERBOSE_CHR_LEVEL, "#\tdata array begins at %#lx, last at %#lx, size: %#lx\n", (unsigned long) dataArrayPtr, (unsigned long) (((unsigned long) dataArrayPtr) + ((unsigned long) size)), (unsigned long) size); } } else if (doDataArray && !(doNoOp || doAscii || doBed || doStats)) { verbose(VERBOSE_CHR_LEVEL, "getData: must specify a chrom range to FetchDataArray\n"); return(valuesMatched); /* NOTHING ASKED FOR ? */ } setDbTable(wds, db, table); openWigConn(wds); if (doStats && wds->useDataConstraint) summaryOnly = FALSE; if (doDataArray) summaryOnly = FALSE; if (doAscii) summaryOnly = FALSE; if (doBed) summaryOnly = FALSE; if (wds->winEnd) summaryOnly = FALSE; /* nextRow() produces the next SQL row from either DB or file. * * The unusual use of wiggleFree() as the third element of this for() * loop is to take care of the condition when the 'continue;' * statements are used to skip the loop and not have to worry about * what happens at the end of the loop. the wiggleFree() must be * done at the end of each loop even in the case of 'continue;' or * not. This is the single instance of wiggleFree() in the entire * source file, as it should be. */ for ( ; (!maxReached) && nextRow(wds, row, WIGGLE_NUM_COLS); wiggleFree(&wiggle) ) { struct wigAsciiData *wigAscii = NULL; struct asciiDatum *asciiOut = NULL; /* to address data[] in wigAsciiData */ unsigned chromPosition; boolean range0TakesAll = FALSE; if (doRawStats) { accumStats(wds, lowerLimit, upperLimit, sumData, sumSquares, statsCount, chromStart, chromEnd); resetStats(&lowerLimit, &upperLimit, &sumData, &sumSquares, &statsCount, &chromStart, &chromEnd); } ++rowCount; wiggle = wiggleLoad(row); /* constraints have to be done manually for files */ if (wds->isFile) { /* if chrName is not correct, or span is not correct, * or outside of window, go to next SQL row */ if ( ((wds->chrName) && (differentString(wds->chrName,wiggle->chrom))) || ((wds->spanLimit) && (wds->spanLimit != wiggle->span)) || ((wds->winEnd) && ((wiggle->chromStart >= wds->winEnd) || (wiggle->chromEnd < wds->winStart)) ) ) { bytesSkipped += wiggle->count; continue; /* next SQL row */ } } chromPosition = wiggle->chromStart; /* this will be true the very first time for both reasons */ if ( (wds->currentSpan != wiggle->span) || (wds->currentChrom && differentString(wds->currentChrom, wiggle->chrom))) { /* if we have been working on one, then doBed */ if (!firstSpanDone && doBed && wds->currentChrom) { if (bedElEnd > bedElStart) { struct bed *bedEl; bedEl = bedElement(wds->currentChrom, bedElStart, bedElEnd, ++bedElCount); slAddHead(&wds->bed, bedEl); } bedElStart = 0; bedElEnd = 0; bedElCount = 0; } if (wds->currentSpan && (wds->currentSpan != wiggle->span)) firstSpanDone = TRUE; if (wds->currentChrom && differentString(wds->currentChrom, wiggle->chrom)) firstSpanDone = FALSE; /* if we have been working on one, then doStats */ if (doStats && wds->currentChrom) { accumStats(wds, lowerLimit, upperLimit, sumData, sumSquares, statsCount, chromStart, chromEnd); resetStats(&lowerLimit, &upperLimit, &sumData, &sumSquares, &statsCount, &chromStart, &chromEnd); } freeMem(wds->currentChrom); /* set them whether they are changing or not, doesn't matter */ wds->currentChrom = cloneString(wiggle->chrom); wds->currentSpan = wiggle->span; } /* it may seem inefficient to make one of these for each SQL row, * but the alternative would require expanding the data area for * each row and thus a re-alloc - probably more expensive than * just making one of these for each row. */ if (doAscii) { AllocVar(wigAscii); wigAscii->chrom = cloneString(wiggle->chrom); wigAscii->span = wiggle->span; wigAscii->count = 0; /* will count up as values added */ wigAscii->dataRange = wiggle->dataRange; /* for resolution calc */ wigAscii->data = (struct asciiDatum *) needMem((size_t) (sizeof(struct asciiDatum) * wiggle->validCount)); asciiOut = wigAscii->data; slAddHead(&wds->ascii,wigAscii); } verbose(VERBOSE_SQL_ROW_LEVEL, "#\trow: %llu, start: %u, data range: %g: [%g:%g], resolution: %g per bin\n", rowCount, wiggle->chromStart, wiggle->dataRange, wiggle->lowerLimit, wiggle->lowerLimit+wiggle->dataRange, wiggle->dataRange/(double)MAX_WIG_VALUE); if (wds->useDataConstraint) { boolean takeIt = FALSE; switch (wds->wigCmpSwitch) { case wigNoOp_e: takeIt = TRUE; break; case wigInRange_e: takeIt = (wiggle->lowerLimit < wds->limit_1) && ((wiggle->lowerLimit + wiggle->dataRange) >= wds->limit_0); break; case wigLessThan_e: takeIt = wiggle->lowerLimit < wds->limit_0; break; case wigLessEqual_e: takeIt = wiggle->lowerLimit <= wds->limit_0; break; case wigEqual_e: takeIt = (wiggle->lowerLimit <= wds->limit_0) && (wds->limit_0 <= (wiggle->lowerLimit + wiggle->dataRange)); break; case wigNotEqual_e: takeIt = (wds->limit_0 < wiggle->lowerLimit) || ((wiggle->lowerLimit + wiggle->dataRange) < wds->limit_0); break; case wigGreaterEqual_e: takeIt = wds->limit_0 <= (wiggle->lowerLimit + wiggle->dataRange); break; case wigGreaterThan_e: takeIt = wds->limit_0 < (wiggle->lowerLimit + wiggle->dataRange); break; default: errAbort("wig_getData: illegal wigCmpSwitch ? %#x", wds->wigCmpSwitch); break; } /* switch (wds->wigCmpSwitch) */ if (0.0 == wiggle->dataRange) range0TakesAll = takeIt; if (takeIt) setCompareByte(wds, wiggle->lowerLimit, wiggle->dataRange); else { verbose(VERBOSE_HIGHEST, "#\tthis row fails compare, next SQL row\n"); bytesSkipped += wiggle->count; continue; /* next SQL row */ } } if (summaryOnly && doStats) { /* record maximum limits of data values */ if (wiggle->lowerLimit < lowerLimit) lowerLimit = wiggle->lowerLimit; if (wiggle->lowerLimit+wiggle->dataRange > upperLimit) upperLimit = wiggle->lowerLimit+wiggle->dataRange; sumData += wiggle->sumData; sumSquares += wiggle->sumSquares; /* record maximum extents */ if ((chromStart < 0)||(chromStart > wiggle->chromStart)) chromStart = wiggle->chromStart; if (chromEnd < wiggle->chromEnd) chromEnd = wiggle->chromEnd; statsCount += wiggle->validCount; valuesMatched += wiggle->validCount; bytesSkipped += wiggle->count; continue; /* next SQL row */ } if (!skipDataRead) { int j; /* loop counter through readData */ unsigned char *datum; /* to walk through readData bytes */ unsigned char *readData; /* the bytes read in from the file */ size_t bytesToRead; /* to check read validity */ size_t bytesRead; /* to check read validity */ openWibFile(wds, wiggle->file); /* possibly open a new wib file */ readData = (unsigned char *) needMem((size_t) (wiggle->count + 1)); wds->bytesRead += wiggle->count; udcSeek(wds->wibFH, wiggle->offset); bytesToRead = (size_t) wiggle->count * (size_t) sizeof(unsigned char); bytesRead = udcRead(wds->wibFH, readData, bytesToRead); if (bytesToRead != bytesRead) { errAbort("wig_getData: failed to read %llu bytes from %s\n", (unsigned long long)bytesToRead, wiggle->file); } verbose(VERBOSE_PER_VALUE_LEVEL, "#\trow: %llu, reading: %u bytes\n", rowCount, wiggle->count); /* The third element of the for() statement takes care of the end * of loop operations for the case of the 'continue;' * statement as well as the normal end of loop business */ datum = readData; for (j = 0; j < wiggle->count; ++datum, chromPosition += wiggle->span, ++j) { if (wds->maxOutput) { if (doAscii) { /* previous results to this fetch are wds->valuesMatched * this result is currently at valuesMatched */ if (wds->maxOutput <= (wds->valuesMatched + valuesMatched)) { maxReached = TRUE; break; } } if (doBed) { /* previous results to this fetch are * wds->totalBedElements, this result is currently * at bedElCount. The -1 is because there is * one more bed element to go upon exit from * this loop. */ if ( (wds->totalBedElements + bedElCount) >= (wds->maxOutput - 1) ) { maxReached = TRUE; break; } } } if (*datum != WIG_NO_DATA) { boolean takeIt = FALSE; if (wds->winEnd) /* non-zero means a range is in effect */ { /* do not allow item (+span) to run over winEnd */ unsigned span = wiggle->span; /* if only doing data array, span is 1 */ if (doDataArray && !(doStats || doBed || doStats)) span = 1; if ( (chromPosition < wds->winStart) || ((chromPosition+span) > wds->winEnd) ) continue; /* next *datum */ } ++validData; switch (wds->wigCmpSwitch) { case wigNoOp_e: takeIt = TRUE; break; case wigInRange_e: takeIt = (*datum < wds->ucUpperLimit) && (*datum >= wds->ucLowerLimit); break; case wigLessThan_e: takeIt = *datum < wds->ucLowerLimit; break; case wigLessEqual_e: takeIt = *datum <= wds->ucLowerLimit; break; case wigEqual_e: takeIt = *datum == wds->ucLowerLimit; break; case wigNotEqual_e: takeIt = *datum != wds->ucLowerLimit; break; case wigGreaterEqual_e: takeIt = *datum >= wds->ucLowerLimit; break; case wigGreaterThan_e: takeIt = *datum > wds->ucLowerLimit; break; default: errAbort("wig_getData: illegal wigCmpSwitch ? %#x", wds->wigCmpSwitch); break; } /* switch (wds->wigCmpSwitch) */ if (range0TakesAll || takeIt) { float value = 0.0; ++valuesMatched; if (doAscii || doStats || doDataArray || doNoOp) value = BIN_TO_VALUE(*datum, wiggle->lowerLimit,wiggle->dataRange); if (doAscii) { asciiOut->value = value; asciiOut->chromStart = chromPosition; ++asciiOut; ++wigAscii->count; } if (!firstSpanDone && doDataArray) { unsigned span = wiggle->span; unsigned inx; long long gap = (long long) chromPosition - (long long) dataArrayPosition; if (verboseLevel() >= VERBOSE_CHR_LEVEL) { if (chromPosition < dataArrayPosition) verbose(VERBOSE_CHR_LEVEL, "#\tWARNING: coordinates out of order %s %u < %u\n", wiggle->chrom, chromPosition, dataArrayPosition); } dataArrayPosition += gap; /* special case squeezing into the end */ if ((dataArrayPosition + span) > wds->winEnd) span = wds->winEnd - dataArrayPosition; dataArrayPtr += gap; /* move up, leave NaN's */ verbose(VERBOSE_PER_VALUE_LEVEL, "#\t%u-%u <- %g\n", dataArrayPosition, dataArrayPosition + span, value); /* fill all chrom positions for this element */ for (inx = 0; inx < span; ++inx) *dataArrayPtr++ = value; dataArrayPosition += span; } /* beware, coords must come in sequence for this to * work. The original SQL query should be ordering * things properly. We only do the first span * encountered, it should be the lowest due to the * SQL query. */ if (!firstSpanDone && doBed) { /* as long as they are continuous, keep * extending the end */ if (chromPosition == bedElEnd) bedElEnd += wiggle->span; else { /* if one was collected, output it */ if (bedElEnd > bedElStart) { struct bed *bedEl; bedEl = bedElement(wiggle->chrom, bedElStart, bedElEnd, ++bedElCount); slAddHead(&wds->bed, bedEl); } /* start a new element here */ bedElStart = chromPosition; bedElEnd = bedElStart + wiggle->span; } } if (doStats) { /* record maximum data limits */ if (value < lowerLimit) lowerLimit = value; if (value > upperLimit) upperLimit = value; sumData += value; sumSquares += value * value; /* record maximum extents */ if ((chromStart < 0)||(chromStart > chromPosition)) chromStart = chromPosition; if (chromEnd < (chromPosition + wiggle->span)) chromEnd = chromPosition + wiggle->span; ++statsCount; } } } /* if (*datum != WIG_NO_DATA) */ else { ++noDataBytes; } } /* for (j = 0; j < wiggle->count; ++j) */ freeMem(readData); } /* if (!skipDataRead) */ } /* for ( ; nextRow(wds, row, WIGGLE_NUM_COLS); ... ) */ /* there may be one last bed element to output */ if (!firstSpanDone && doBed) { if (bedElEnd > bedElStart) { struct bed *bedEl; bedEl = bedElement(wds->currentChrom, bedElStart, bedElEnd, ++bedElCount); slAddHead(&wds->bed, bedEl); } } /* there are accumulated stats to complete */ if (doStats) { accumStats(wds, lowerLimit, upperLimit, sumData, sumSquares, statsCount, chromStart, chromEnd); resetStats(&lowerLimit, &upperLimit, &sumData, &sumSquares, &statsCount, &chromStart, &chromEnd); } wds->rowsRead += rowCount; wds->validPoints += validData; wds->valuesMatched += valuesMatched; wds->noDataPoints += noDataBytes; wds->bytesSkipped += bytesSkipped; wds->totalBedElements += bedElCount; closeWigConn(wds); return(valuesMatched); } /* unsigned long long getData() */ static float *asciiToDataArray(struct wiggleDataStream *wds, unsigned long long count, size_t *returned) /* convert the AsciiData list to a float array */ { float *floatArray = NULL; float *fptr = NULL; unsigned long long valuesDone = 0; struct wigAsciiData *asciiData = NULL; if (count < 1) return floatArray; setMaxAlloc((size_t)HUGE_LIMIT); AllocArray(floatArray, count); fptr = floatArray; /* the (valuesDone <= count) condition is for safety */ for (asciiData = wds->ascii; asciiData && (valuesDone < count); asciiData = asciiData->next) { if (asciiData->count) { struct asciiDatum *data; unsigned i; data = asciiData->data; for (i = 0; (i < asciiData->count)&&(valuesDone < count); ++i) { *fptr++ = data->value; ++data; ++valuesDone; } } } if (returned) *returned = valuesDone; return floatArray; } static unsigned long long getDataViaBed(struct wiggleDataStream *wds, char *db, char *table, int operations, struct bed **bedList) /* getDataViaBed - constrained by the bedList */ { unsigned long long valuesMatched = 0; unsigned long long prevValuesMatched = 0; if (bedList && *bedList) { struct bed *filteredBed = NULL; boolean chromConstraint = ((char *)NULL != wds->chrName); boolean positionConstraints = (0 != wds->winEnd); int saveWinStart = 0; int saveWinEnd = 0; char * saveChrName = NULL; struct slName *chromList = NULL; struct slName *chr; struct hash *chromSizes = NULL; struct { unsigned long chrStart; unsigned long chrEnd; } *chrStartEnd; unsigned long long valuesFound = 0; boolean maxReached = FALSE; unsigned long long prevMaxOutput = 0; /* remember these constraints so we can reset them afterwards * because we are going to use them for our own internal uses here. */ saveWinStart = wds->winStart; saveWinEnd = wds->winEnd; chromSizes = newHash(0); /* Determine what will control which chrom to do. If a chrom * constraint has been set then that is the one to do. If it has * not been set, then allow the bedList to be the list of chroms to * do. This business finds the maximum extent for each chrom in * the bedList. */ if (chromConstraint) { saveChrName = cloneString(wds->chrName); chr = newSlName(wds->chrName); slAddHead(&chromList, chr); AllocVar(chrStartEnd); chrStartEnd->chrStart = wds->winStart; chrStartEnd->chrEnd = wds->winEnd; hashAdd(chromSizes, wds->chrName, chrStartEnd); } else { struct bed *bed = NULL; for (bed = *bedList; bed; bed = bed->next) { /* out of order bed file ? Maybe already saw this name */ if (hashLookup(chromSizes, bed->chrom) != NULL) { chrStartEnd = hashFindVal(chromSizes, bed->chrom); chrStartEnd->chrStart = min(bed->chromStart,chrStartEnd->chrStart); chrStartEnd->chrEnd = max(bed->chromEnd, chrStartEnd->chrEnd); } else { chr = newSlName(bed->chrom); slAddHead(&chromList, chr); AllocVar(chrStartEnd); chrStartEnd->chrStart = bed->chromStart; chrStartEnd->chrEnd = bed->chromEnd; hashAdd(chromSizes, bed->chrom, chrStartEnd); } } slReverse(&chromList); } /* OK, the chromList and chromSizes hash have been created, run * through that list. This may seem unusual, but if the caller has * set position constraints but no chr constraint, allow those * constraints to override * the bedList coordinates, even though we may be running through a * number of chroms. That's what they wanted, that's what they get. * e.g. they wanted the first 10,000 bases of each chrom, * winStart=0, winEnd = 10,000, no chrName specified */ for (chr = chromList; (!maxReached) && chr; chr = chr->next) { long startTime = 0; float *fptr; boolean doAscii = FALSE; boolean doBed = FALSE; boolean doStats = FALSE; boolean doNoOp = FALSE; boolean doDataArray = FALSE; struct bed *bed; unsigned long bedExtent = 0; /* from first thru last element of bedList*/ unsigned long bedStart = 0; unsigned long bedEnd = 0; unsigned long bedPosition = 0; char *bedArray; /* meaning is boolean but */ char *boolPtr; /* boolean is 4 bytes, char is 1 */ unsigned long bedArraySize = 0; unsigned long winStart = 0; unsigned long winEnd = 0; unsigned long winExtent = 0; unsigned long bedMarked = 0; unsigned long asciiDataSizeLimit = 0; unsigned long chromPosition = 0; unsigned long bedElStart = 0; unsigned long bedElEnd = 0; unsigned long bedElCount = 0; float lowerLimit = INFINITY; /* stats accumulators */ float upperLimit = -1.0 * INFINITY; /* stats accumulators */ double sumData = 0.0; /* stats accumulators */ double sumSquares = 0.0; /* stats accumulators */ unsigned long statsCount = 0; /* stats accumulators */ long int chromStart = -1; /* stats accumulators */ long int chromEnd = 0; /* stats accumulators */ size_t dataArraySize = 0; doDataArray = operations & wigFetchDataArray; /* set chrom name constraint */ wds->setChromConstraint(wds, (char *)&chr->name); chrStartEnd = hashFindVal(chromSizes, wds->chrName); if (NULL == chrStartEnd) errAbort("getDataViaBed: constructed hash of bed coords is broken"); /* This may seem unusual, but if the caller has * set position constraints, allow those constraints to override * the bedList coordinates, even though we may be running through * a number of chroms. */ if (saveWinEnd) { winStart = saveWinStart; winEnd = saveWinEnd; } else { winStart = chrStartEnd->chrStart; winEnd = chrStartEnd->chrEnd; } /* winExtent could still be 0 here */ /* remember this window extent */ winExtent = winEnd - winStart; /* set window position constraint */ wds->setPositionConstraint(wds, winStart, winEnd); /* Going to create a TRUE/FALSE marker array from the bedList. * Each chrom position will be marked TRUE if it belongs to * a bed element. Overlapping bed elements will be * redundant and a waste of time, but no real problem. * Filter by the chr name and the win limits. */ for (bed = *bedList; bed; bed = bed->next) { if (sameString(wds->chrName, bed->chrom)) { /* if (no limits) OR (within limits) use it */ if (!(winEnd) || ( (bed->chromStart < winEnd) && (bed->chromEnd >= winStart)) ) { struct bed *newBed = cloneBed(bed); slAddHead(&filteredBed, newBed); } } } if (NULL == filteredBed) { verbose(VERBOSE_CHR_LEVEL, "#\tfilter found nothing in bed file: %s:%lu-%lu\n", wds->chrName, winStart, winEnd); continue; /* next chrom */ } else slSort(&filteredBed,bedCmp); /* in proper order */ if (verboseLevel() >= VERBOSE_CHR_LEVEL) { int count=slCount(filteredBed); verbose(VERBOSE_CHR_LEVEL, "#\tfiltered by chr,pos constraints, bed el count: %d\n", count); } /* determine actual limits loaded by the filter */ bedStart = min(winStart,filteredBed->chromStart); bedEnd = filteredBed->chromEnd; for (bed = filteredBed; bed; bed = bed->next) { if (bed->chromStart < bedStart) bedStart = bed->chromStart; if (bed->chromEnd > bedEnd) bedEnd = bed->chromEnd; } bedExtent = bedEnd - bedStart; /* must use the larger of the two sizes */ bedArraySize = max(winExtent,bedExtent); bedArraySize *= sizeof(char); /* in bytes */ setMaxAlloc((size_t)HUGE_LIMIT); if (verboseLevel() >= VERBOSE_CHR_LEVEL) startTime = clock1000(); bedArray = (char *)needLargeMem(bedArraySize); verbose(VERBOSE_CHR_LEVEL, "#\tbed array begins at %#lx, last at %#lx, size: %#lx\n", (unsigned long) bedArray, (unsigned long) (((unsigned long) bedArray) + ((unsigned long) bedArraySize)), (unsigned long) bedArraySize ); /* set them all to FALSE to begin with */ memset((void *)bedArray, FALSE, bedArraySize); if (verboseLevel() >= VERBOSE_CHR_LEVEL) { long et = clock1000() - startTime; verbose(VERBOSE_CHR_LEVEL, "#\tgetDataViaBed: created data array for %lu bed indicators in %ld ms\n", bedArraySize, et); } /* use the bed definition to fill this marker array */ /* start at the lesser of the two start positions */ bedPosition = min(winStart, bedStart); bedStart = bedPosition; boolPtr = bedArray; bedEnd = filteredBed->chromEnd; for (bed = filteredBed; bed; bed = bed->next) { unsigned long elSize = bed->chromEnd - bed->chromStart; boolPtr = bedArray + (bed->chromStart - bedStart); memset((void *)boolPtr, TRUE, elSize); boolPtr += elSize; bedPosition = (unsigned long)((void *)boolPtr - (void *)bedArray); bedMarked += elSize; /* number of bases marked */ bedEnd = max(bedEnd, bed->chromEnd); } bedExtent = bedEnd - bedStart; verbose(VERBOSE_CHR_LEVEL, "#\tbed range %lu = %lu - %lu, allocated %lu bytes\n", bedExtent, bedEnd, bedStart, bedArraySize); verbose(VERBOSE_CHR_LEVEL, "#\tbed marked %lu bases, winSize %lu\n", bedMarked, winExtent); /* worst case ascii limit: bedMarked */ asciiDataSizeLimit = bedMarked; /*set window position constraint, this eliminates the 0,0 possibility*/ if (winEnd) { winStart = max(bedStart, winStart); winEnd = min(bedEnd, winEnd); } else { winStart = bedStart; winEnd = bedEnd; } /* must have actual positions for constraints for getData * to do this wigFetchDataArray function. */ wds->setPositionConstraint(wds, winStart, winEnd); /* now fetch all the wiggle data for this position, * constraints have been set at the top of this loop */ /* POTENTIAL SPEEDUP here, if only bed results are * the requested result, then we should just ask for a bed * list from getData and simply intersect that bed list * with the given bed list filter. */ prevValuesMatched = wds->valuesMatched; prevMaxOutput = wds->maxOutput; wds->maxOutput = 0; valuesFound = wds->getData(wds, db, table, wigFetchDataArray); wds->valuesMatched = prevValuesMatched; wds->maxOutput = prevMaxOutput; verbose(VERBOSE_CHR_LEVEL, "#\tback from getData, found %llu valid points\n", valuesFound); doAscii = operations & wigFetchAscii; doBed = operations & wigFetchBed; doStats = operations & wigFetchStats; doNoOp = operations & wigFetchNoOp; /* for timing purposes, allow the wigFetchNoOp to go through */ if (doNoOp) doDataArray = doBed = doStats = doAscii = FALSE; /* if ((something found) and (something to do)) */ if ((valuesFound > 0) && (doAscii || doBed || doStats || doNoOp)) { struct wigAsciiData *wigAscii = NULL; struct asciiDatum *asciiOut = NULL; /* to address data[] in wigAsciiData */ dataArraySize = wds->array->winEnd - wds->array->winStart; /* no bed constraint, then worst case ascii limit: dataArraySize */ if (0 == asciiDataSizeLimit) asciiDataSizeLimit = dataArraySize; if (verboseLevel() > VERBOSE_PER_VALUE_LEVEL) dataArrayOut(wds, "stdout", FALSE, TRUE); /* OK, ready to scan the bedArray in concert with the data * array and pick out those values that are masked by the bed. * The starting position will be the greater of the * two different possible windows. */ if (doAscii) { AllocVar(wigAscii); wigAscii->chrom = cloneString(wds->chrName); wigAscii->span = 1; /* span information has been lost */ wigAscii->count = 0; /* will count up as values added */ wigAscii->dataRange = 0.0; /* to be determined */ setMaxAlloc((size_t)( HUGE_LIMIT * (((sizeof(size_t)/4)*(sizeof(size_t)/4)*(sizeof(size_t)/4))))); /* produces: size_t is 4 ~= 4 Gb * size_t is 8 ~= 32 Gb */ verbose(VERBOSE_CHR_LEVEL, "#\tworst case ascii array needLargeMem (%llu * %llu = %llu)\n", (unsigned long long)sizeof(struct asciiDatum), (unsigned long long)asciiDataSizeLimit, (unsigned long long)(sizeof(struct asciiDatum) * asciiDataSizeLimit)); wigAscii->data = (struct asciiDatum *) needLargeMem((size_t) (sizeof(struct asciiDatum) * asciiDataSizeLimit)); /* maximum area needed, may use less */ asciiOut = wigAscii->data; /* ptr to data area */ slAddHead(&wds->ascii,wigAscii); } chromPosition = max(bedStart, winStart); boolPtr = bedArray + (chromPosition - bedStart); fptr = wds->array->data + (chromPosition - winStart); for ( ; (!maxReached) && (chromPosition < winEnd); ++chromPosition, ++boolPtr, ++fptr) { if (wds->maxOutput) { if (doAscii) { /* previous results to this fetch are wds->valuesMatched * this result is currently at valuesMatched */ if (wds->maxOutput <= (wds->valuesMatched + valuesMatched)) { maxReached = TRUE; break; } } if (doBed) { /* previous results to this fetch are * wds->totalBedElements, this result is currently * at bedElCount. The -1 is because there is * one more bed element to go upon exit from * this loop. */ if ( (wds->totalBedElements + bedElCount) >= (wds->maxOutput - 1) ) { maxReached = TRUE; break; } } } if (*boolPtr && (!isnan(*fptr))) { float value = *fptr; ++valuesMatched; /* construct output listings here */ if (doAscii) { /* record limits when not also doing that * for stats */ if (!doStats) { if (value < lowerLimit) lowerLimit = value; if (value > upperLimit) upperLimit = value; } asciiOut->value = value; asciiOut->chromStart = chromPosition; if (wigAscii->count < asciiDataSizeLimit) { ++asciiOut; ++wigAscii->count; } else maxReached = TRUE; } if (doBed) { /* have we started a bed element yet ? */ if (0 == bedElEnd) bedElEnd = bedElStart = chromPosition; /* That will get it started right here */ /* as long as they are continuous, keep * extending the end */ if (chromPosition == bedElEnd) bedElEnd += 1; /* span information has been lost */ else { /* if one was collected, output it */ if (bedElEnd > bedElStart) { struct bed *bedEl; bedEl = bedElement(wds->chrName, bedElStart, bedElEnd, ++bedElCount); slAddHead(&wds->bed, bedEl); } /* start a new element here */ bedElStart = chromPosition; bedElEnd = bedElStart + 1; /* span = 1 */ } } /* if (doBed) */ if (doStats) { /* record maximum data limits */ if (value < lowerLimit) lowerLimit = value; if (value > upperLimit) upperLimit = value; sumData += value; sumSquares += value * value; /* record maximum extents */ if ((chromStart < 0)||(chromStart > chromPosition)) chromStart = chromPosition; if (chromEnd < (chromPosition + 1)) /* span = 1 */ chromEnd = chromPosition + 1; /* span = 1 */ ++statsCount; } } /* if (*boolPtr && (!isnan(*fptr))) */ } /* if we can save more than 100,000 data locations, move * the results to a smaller area. Could be huge savings on * very sparse result sets. This comes about because at * the beginning we blindly allocated a full * winEnd-winStart data array and we probably did not use * all of that. */ if (doAscii && ((asciiDataSizeLimit - wigAscii->count) > 100000)) { struct asciiDatum *smallerDataArea; size_t newSize; newSize = sizeof(struct asciiDatum) * wigAscii->count; if (newSize > 0) { setMaxAlloc((size_t)( HUGE_LIMIT * (((sizeof(size_t)/4)*(sizeof(size_t)/4)*(sizeof(size_t)/4))))); /* produces: size_t is 4 ~= 4 Gb * size_t is 8 ~= 32 Gb */ verbose(VERBOSE_CHR_LEVEL, "#\tmoving to smaller ascii array needLargeMem( %llu * %llu = %llu)\n", (unsigned long long)sizeof(struct asciiDatum), (unsigned long long)wigAscii->count, (unsigned long long)(sizeof(struct asciiDatum) * wigAscii->count)); smallerDataArea=(struct asciiDatum *)needLargeMem(newSize); memcpy(smallerDataArea, wigAscii->data, newSize); freeMem(wigAscii->data); wigAscii->data = smallerDataArea; } else { freeMem(wigAscii->data); wigAscii->data = NULL; } } /* record dataRange for resolution accounting */ if (doAscii && wigAscii->count) wigAscii->dataRange = (upperLimit - lowerLimit); if (doBed) /* there may be one last element */ { if (bedElEnd > bedElStart) { struct bed *bedEl; bedEl = bedElement(wds->chrName, bedElStart, bedElEnd, ++bedElCount); slAddHead(&wds->bed, bedEl); } } /* there are accumulated stats to complete */ if (doStats) { accumStats(wds, lowerLimit, upperLimit, sumData, sumSquares, statsCount, chromStart, chromEnd); if (wds->stats) wds->stats->span = 1; /* viaBed reduces span to 1 */ resetStats(&lowerLimit, &upperLimit, &sumData, &sumSquares, &statsCount, &chromStart, &chromEnd); } verbose(VERBOSE_CHR_LEVEL, "#\tviaBed found: %llu valid data points\n", valuesMatched); } /* if ((something found) and (something to do)) */ /* did they want the data array returned ? No, then release it */ if (!doDataArray) wds->freeArray(wds); /* we are certainly done with the bedArray */ bedFreeList(&filteredBed); freeMem(bedArray); wds->totalBedElements += bedElCount; } /* for (chr = chromList; chr; chr = chr->next) */ /* restore these constraints we used here */ if (chromConstraint) { wds->setChromConstraint(wds, saveChrName); freeMem(saveChrName); } else { freeMem(wds->chrName); wds->chrName = NULL; } if (positionConstraints) wds->setPositionConstraint(wds, saveWinStart, saveWinEnd); else wds->setPositionConstraint(wds, 0, 0); /* The bedConstrained flag will provide a printout indication of * what happened to the *Out() routines. */ wds->bedConstrained = TRUE; /* to signal the *Out() displays */ freeHashAndVals(&chromSizes); } /* if (bedList && *bedList) */ wds->valuesMatched += valuesMatched; return(valuesMatched); } /* static unsigned long long getDataViaBed() */ static void sortResults(struct wiggleDataStream *wds) /* sort any results that exist */ { if (wds->bed) slSort(&wds->bed, bedCmp); if (wds->ascii) slSort(&wds->ascii, asciiDataCmp); if (wds->stats) slSort(&wds->stats, statsDataCmp); if (wds->array) slSort(&wds->array, arrayDataCmp); } static int bedOut(struct wiggleDataStream *wds, char *fileName, boolean sort) /* print to fileName the bed list */ { int linesOut = 0; FILE * fh; fh = mustOpen(fileName, "w"); if (wds->bed) { struct bed *bedEl, *next; outputIdentification(wds, fh); showConstraints(wds, fh); if (sort) slSort(&wds->bed, bedCmp); for (bedEl = wds->bed; bedEl; bedEl = next ) { fprintf(fh,"%s\t%u\t%u\t%s\n", bedEl->chrom, bedEl->chromStart, bedEl->chromEnd, bedEl->name); next = bedEl->next; ++linesOut; } } else { showConstraints(wds, fh); fprintf(fh, "#\tbed: no data points found for bed format output (maxOutput: %llu)\n", wds->maxOutput); } carefulClose(&fh); return (linesOut); } /* static void bedOut() */ static void statsOut(struct wiggleDataStream *wds, char *db, char *fileName, boolean sort, boolean htmlOut, boolean withHeader, boolean leaveTableOpen) /* print to fileName the statistics */ { FILE * fh; fh = mustOpen(fileName, "w"); if (wds->stats) { struct wiggleStats *stats, *next; int itemsDisplayed = 0; if (sort) slSort(&wds->stats, statsDataCmp); if (withHeader) wigStatsHeader(wds, fh, htmlOut); for (stats = wds->stats; stats; stats = next ) { if (htmlOut) { long long chromSize = 0; if (wds->winEnd) chromSize = wds->winEnd - wds->winStart; else { if (! wds->isFile) chromSize = hChromSize(db, stats->chrom); } fprintf(fh,"<TR><TH ALIGN=LEFT> %s </TH>\n", stats->chrom); fprintf(fh,"\t<TD ALIGN=RIGHT> %u </TD>\n", BASE_1(stats->chromStart)); /* display closed coords */ fprintf(fh,"\t<TD ALIGN=RIGHT> %u </TD>\n", stats->chromEnd); fprintf(fh,"\t<TD ALIGN=RIGHT> "); printLongWithCommas(fh, (long long)stats->count); fprintf(fh,"</TD>\n"); fprintf(fh,"\t<TD ALIGN=RIGHT> %d </TD>\n", stats->span); fprintf(fh,"\t<TD ALIGN=RIGHT> "); printLongWithCommas(fh, (long long) stats->count*stats->span); if (chromSize > 0) fprintf(fh," (%.2f%%) </TD>\n", 100.0*(double)(stats->count*stats->span)/(double)chromSize); else fprintf(fh,"</TD>\n"); fprintf(fh,"\t<TD ALIGN=RIGHT> %g </TD>\n", stats->lowerLimit); fprintf(fh,"\t<TD ALIGN=RIGHT> %g </TD>\n", stats->lowerLimit+stats->dataRange); fprintf(fh,"\t<TD ALIGN=RIGHT> %g </TD>\n", stats->dataRange); fprintf(fh,"\t<TD ALIGN=RIGHT> %g </TD>\n", stats->mean); fprintf(fh,"\t<TD ALIGN=RIGHT> %g </TD>\n", stats->variance); fprintf(fh,"\t<TD ALIGN=RIGHT> %g </TD>\n", stats->stddev); fprintf(fh,"</TR>\n"); } else { fprintf(fh,"%s", stats->chrom); fprintf(fh,"\t%u", BASE_1(stats->chromStart)); /* display closed coords */ fprintf(fh,"\t%u", stats->chromEnd); fprintf(fh,"\t%u", stats->count); fprintf(fh,"\t%d", stats->span); fprintf(fh,"\t%u", stats->count*stats->span); fprintf(fh,"\t%g", stats->lowerLimit); fprintf(fh,"\t%g", stats->lowerLimit+stats->dataRange); fprintf(fh,"\t%g", stats->dataRange); fprintf(fh,"\t%g", stats->mean); fprintf(fh,"\t%g", stats->variance); fprintf(fh,"\t%g", stats->stddev); fprintf(fh,"\n"); } ++itemsDisplayed; next = stats->next; } if (!itemsDisplayed) { if (htmlOut) fprintf(fh,"<TR><TH ALIGN=CENTER COLSPAN=12> No data found matching this request </TH></TR>\n"); else fprintf(fh,"#\tNo data found matching this request\n"); } if (!leaveTableOpen && withHeader && htmlOut) fprintf(fh,"</TABLE></TD></TR></TABLE></P>\n"); } else { if (!htmlOut) { showConstraints(wds, fh); fprintf(fh, "#\tstats: no data points found\n"); } } carefulClose(&fh); } /* static void statsOut() */ static int asciiOut(struct wiggleDataStream *wds, char *db, char *fileName, boolean sort, boolean rawDataOut) /* print to fileName the ascii data values */ { boolean firstLine = TRUE; int linesOut = 0; FILE * fh; fh = mustOpen(fileName, "w"); if (wds->ascii) { struct wigAsciiData *asciiData; char *chrom = NULL; unsigned span = 0; double worstCaseResolution = 0.0; if (sort) slSort(&wds->ascii, asciiDataCmp); /* determine an over-all resolution measurement */ for (asciiData = wds->ascii; asciiData; asciiData = asciiData->next) { double resolution = asciiData->dataRange / (MAX_WIG_VALUE+1); if (resolution > worstCaseResolution) worstCaseResolution = resolution; } for (asciiData = wds->ascii; asciiData; asciiData = asciiData->next) { unsigned i; struct asciiDatum *data; if (asciiData->count) { /* will be true the first time for both reasons */ if ( (span != asciiData->span) || (chrom && differentString(asciiData->chrom, chrom)) ) { freeMem(chrom); chrom = cloneString(asciiData->chrom); span = asciiData->span; if (!rawDataOut) { if (firstLine) { char *dateStamp = dateTimeStamp(); fprintf (fh, "#\toutput date: %s UTC\n", dateStamp); firstLine = FALSE; freeMem(dateStamp); } showConstraints(wds, fh); if (wds->tblName && startsWith("gc5", wds->tblName)) showResolutionNoDownloads(worstCaseResolution, fh); else showResolution(worstCaseResolution, fh); fprintf (fh, "variableStep chrom=%s span=%u\n", chrom, span); } } data = asciiData->data; for (i = 0; i < asciiData->count; ++i) { /* user visible coordinates are wds->offset relative/based */ if (rawDataOut) fprintf (fh, "%g\n", data->value); else fprintf (fh, "%u\t%g\n", data->chromStart + wds->offset, data->value); ++data; ++linesOut; } } } } else { if (!rawDataOut) { showConstraints(wds, fh); fprintf(fh, "#\tasciiOut: no data points found\n"); } } carefulClose(&fh); return (linesOut); } /* static void asciiOut() */ void wiggleDataStreamFree(struct wiggleDataStream **wds) /* free all structures and zero the callers' structure pointer */ { if (wds) { struct wiggleDataStream *wdstream; wdstream=*wds; if (wdstream) { closeWigConn(wdstream); wdstream->freeAscii(wdstream); wdstream->freeBed(wdstream); wdstream->freeStats(wdstream); wdstream->freeArray(wdstream); wdstream->freeConstraints(wdstream); freeMem(wdstream->currentChrom); } freez(wds); } } struct wiggleDataStream *wiggleDataStreamNew() { struct wiggleDataStream *wdstream; AllocVar(wdstream); /* everything is zero which is good since that is NULL for all the * strings and lists. A few items should have some initial values * which are not necessarily NULL */ wdstream->isFile = FALSE; wdstream->useDataConstraint = FALSE; wdstream->wibFH = (struct udcFile*)NULL; wdstream->limit_0 = -1 * INFINITY; wdstream->limit_1 = INFINITY; wdstream->wigCmpSwitch = wigNoOp_e; wdstream->freeConstraints = freeConstraints; wdstream->freeAscii = freeAscii; wdstream->freeBed = freeBed; wdstream->freeStats = freeStats; wdstream->freeArray = freeArray; wdstream->setMaxOutput = setMaxOutput; wdstream->setPositionConstraint = setPositionConstraint; wdstream->setChromConstraint = setChromConstraint; wdstream->setSpanConstraint = setSpanConstraint; wdstream->setDataConstraint = setDataConstraint; wdstream->bedOut = bedOut; wdstream->statsOut = statsOut; wdstream->asciiOut = asciiOut; wdstream->sortResults = sortResults; wdstream->asciiToDataArray = asciiToDataArray; wdstream->getDataViaBed = getDataViaBed; wdstream->getData = getData; wdstream->maxOutput = 0; wdstream->offset = 1; /* default 1-relative/based output */ return wdstream; }