083833e1f84528f923fde137a041657635f4564e hiram Thu Jul 20 22:29:00 2017 -0700 remove the CPU cache experiment and change memLimit option to threads refs #18969 diff --git src/utils/crisprKmers/crisprKmers.c src/utils/crisprKmers/crisprKmers.c index 0f6a72c..a69e768 100644 --- src/utils/crisprKmers/crisprKmers.c +++ src/utils/crisprKmers/crisprKmers.c @@ -66,60 +66,54 @@ /* Explain usage and exit. */ { errAbort( "crisprKmers - annotate crispr guide sequences\n" "usage:\n" " crisprKmers <sequence>\n" "options:\n" " where <sequence> is a .2bit file or .fa fasta sequence\n" " -verbose=N - for debugging control processing steps with level of verbose:\n" " - verbose < 2 - scan sequence only, no comparisons\n" " - verbose > 1 - run comparisons\n" " -bed=<file> - output results to given bed9+ file\n" " -offTargets=<file> - output off target data to given file\n" " -ranges=<file> - use specified bed3 file to limit which guides are\n" " - measured, only those with any overlap to these bed items.\n" - " -memLimit=N - N number of gigabytes for each thread, default: no threading\n" - " when memory size goes beyond N gB program will thread gB/N times\n" - " -cacheCount=N - number of guides to check at one time vs. targets,\n" - " - to take advantage of CPU cache, e.g. 1Mb of cache\n" - " - could hold 16,384 64bit words.\n" + " -threads=N - N number of threads to run, default: no threading\n" + " - use when ku is going to allocate more CPUs for big mem jobs.\n" " -dumpKmers=<file> - NOT VALID after scan of sequence, output kmers to file\n" " - process will exit after this, use -loadKmers to continue\n" " -loadKmers=<file> - NOT VALID load kmers from previous scan of sequence from -dumpKmers" ); } static char *bedFileOut = NULL; /* set with -bed=<file> argument, kmer output */ static char *offTargets = NULL; /* write offTargets to <file> */ static FILE *offFile = NULL; /* file handle to write off targets */ static char *ranges = NULL; /* use ranges <file> to limit scanning */ static struct hash *rangesHash = NULL; /* ranges into hash + binkeeper */ static char *dumpKmers = NULL; /* file name to write out kmers from scan */ static char *loadKmers = NULL; /* file name to read in kmers from previous scan */ -static int memLimit = 0; /* gB limit before going into thread mode */ -static int cacheCount = 0; /* number of query guides to run through - * for each target to utilize CPU cache */ +static int threads = 0; /* number of threads to run */ static int threadCount = 1; /* will be adjusted depending upon vmPeak */ /* Command line validation table. */ static struct optionSpec options[] = { {"bed", OPTION_STRING}, {"offTargets", OPTION_STRING}, {"ranges", OPTION_STRING}, - {"memLimit", OPTION_INT}, - {"cacheCount", OPTION_INT}, + {"threads", OPTION_INT}, {"dumpKmers", OPTION_STRING}, {"loadKmers", OPTION_STRING}, {NULL, 0}, }; #define guideSize 20 // 20 bases #define pamSize 3 // 3 bases #define negativeStrand 0x0000800000000000 #define fortyEightBits 0xffffffffffff #define fortySixBits 0x3fffffffffff #define fortyBits 0xffffffffff #define high32bits 0xffffffff00000000 #define low32bits 0x00000000ffffffff #define high16bits 0xffff0000ffff0000 @@ -518,31 +512,34 @@ char pamString[33]; kmerPAMString(pamString, list->sequence[c]); if (bedFH) { long mitScore = 1.0; if (mitScoreCount > 0) /* note: interger <- float conversion */ { mitScore = roundf((list->mitSum[c] / (float) mitScoreCount)); } else if (list->offBy[0][c] > 0) { mitScore = 0.0; } char *color = scoreToColor(mitScore); - fprintf(bedFH, "%s\t%lld\t%lld\t\t%ld\t%c\t%lld\t%lld\t%s\t%s\t%s\t%ld\t%d,%d,%d,%d,%d\n", list->chrom, list->start[c], list->start[c]+pamSize+guideSize, mitScore, strand, txStart, txEnd, color, kmerString, pamString, mitScore, list->offBy[0][c], list->offBy[1][c], list->offBy[2][c], list->offBy[3][c], list->offBy[4][c]); + /* the end zero will be filled in by post-processing, it will be + * the _offset into the crisprDetails.tab file + */ + fprintf(bedFH, "%s\t%lld\t%lld\t\t%ld\t%c\t%lld\t%lld\t%s\t%s\t%s\t%d,%d,%d,%d,%d\t0\n", list->chrom, list->start[c], list->start[c]+pamSize+guideSize, mitScore, strand, txStart, txEnd, color, kmerString, pamString, list->offBy[0][c], list->offBy[1][c], list->offBy[2][c], list->offBy[3][c], list->offBy[4][c]); } ++totalOut; } } if (bedFH) carefulClose(&bedFH); timingMessage("countsOutput:", itemsOutput, "items output", startTime, "items/sec", "seconds/item"); } // static void countsOutput(struct crisprList *all, FILE *bedFH) static struct crisprList *scanSequence(char *inFile) /* scan the given file, return list of crisprs */ { @@ -976,33 +973,31 @@ query->chrom, query->start[qIndex], negativeStrand & query->sequence[qIndex] ? '-' : '+', queryString, queryPAM, targetString, targetPAM, target->chrom, target->start[tIndex], negativeStrand & target->sequence[tIndex] ? '-' : '+', cfdInt, misMatch, bitsOn, mitScore, cfdScore, target->chrom, target->start[tIndex], negativeStrand & target->sequence[tIndex] ? '-' : '+'); } } } // static void recordOffTargets(struct crisprList *query, // struct crisprList *target, int bitsOn, long long qIndex, // long long tIndex) -/* this queryVsTarget can be used by threads, it should be safe, - * remains to be proven. - */ +/* this queryVsTarget can be used by threads, appears to be safe */ static void queryVsTarget(struct crisprList *query, struct crisprList *target, int threadCount, int threadId) /* run the query crisprs list against the target list in the array structures */ { struct crisprList *qList; long long totalCrisprsQuery = 0; long long totalCrisprsTarget = 0; long long totalCompares = 0; struct loopControl *control = NULL; AllocVar(control); long startTime = clock1000(); for (qList = query; qList; qList = qList->next) { @@ -1011,77 +1006,65 @@ { setLoopEnds(control, qList->crisprCount, threadCount, threadId); if (control->listStart >= qList->crisprCount) continue; /* next chrom, no part to be done for this thread */ } else { control->listStart = 0; control->listEnd = qList->crisprCount; } verbose(1, "# thread %d of %d running %s %lld items [ %lld : %lld )\n", 1 + threadId, threadCount, qList->chrom, control->listEnd - control->listStart, control->listStart, control->listEnd); totalCrisprsTarget += control->listEnd - control->listStart; - long long qStart = control->listStart; - long long qEnd = control->listEnd; long long qCount; - for (qCount = control->listStart; qCount < control->listEnd; qCount = qEnd) - { - qStart = qCount; - qEnd = qStart + cacheCount; - if (cacheCount < 1) - qEnd = qStart + 1; - if (qEnd > control->listEnd) - qEnd = control->listEnd; + for (qCount = control->listStart; qCount < control->listEnd; ++qCount) + { struct crisprList *tList; for (tList = target; tList; tList = tList->next) { long long tCount; totalCompares += tList->crisprCount; for (tCount = 0; tCount < tList->crisprCount; ++tCount) { - long long q; - for (q = qStart; q < qEnd; ++q) - { /* the XOR determine differences in two sequences, the * shift right 6 removes the PAM sequence and * the 'fortyBits &' eliminates the negativeStrand bit */ long long misMatch = fortyBits & - ((qList->sequence[q] ^ tList->sequence[tCount]) >> 6); + ((qList->sequence[qCount] ^ tList->sequence[tCount]) >> 6); if (misMatch) { /* possible misMatch bit values: 01 10 11 * turn those three values into just: 01 */ misMatch = (misMatch | (misMatch >> 1)) & 0x5555555555; int bitsOn = _mm_popcnt_u64(misMatch); if (bitsOn < 5) { - recordOffTargets(qList, tList, bitsOn, q, tCount, misMatch); - qList->offBy[bitsOn][q] += 1; + recordOffTargets(qList, tList, bitsOn, qCount, tCount, misMatch); + qList->offBy[bitsOn][qCount] += 1; // tList->offBy[bitsOn][tCount] += 1; not needed } } else { /* no misMatch, identical crisprs */ - qList->offBy[0][q] += 1; + qList->offBy[0][qCount] += 1; // tList->offBy[0][tCount] += 1; not needed } - } // for (q = qStart; q < qEnd; ++q) } // for (tCount = 0; tCount < tList->crisprCount; ++tCount) } // for (tList = target; tList; tList = tList->next) } // for (qCount = 0; qCount < qList->crisprCount; ++qCount) } // for (qList = query; qList; qList = qList->next) timingMessage("queryVsTarget", totalCrisprsQuery, "query crisprs processed", startTime, "crisprs/sec", "seconds/crispr"); timingMessage("queryVsTarget", totalCrisprsTarget, "vs target crisprs", startTime, "crisprs/sec", "seconds/crispr"); timingMessage("queryVsTarget", totalCompares, "total comparisons", startTime, "compares/sec", "seconds/compare"); } /* static struct crisprList *queryVsTarget(struct crisprList *query, struct crisprList *target) */ @@ -1304,34 +1287,35 @@ offFile = mustOpen(offTargets, "w"); /* if ranges have been specified, construct queryList of kmers to measure */ if (rangesHash) { /* result here is two exclusive sets: query, and allGuides * the query crisprs have been extracted from the allGuides */ queryGuides = rangeExtraction(& allGuides); } if (queryGuides) copyToArray(queryGuides); if (allGuides) copyToArray(allGuides); vmPeak = currentVmPeak(); verbose(1, "# vmPeak after copyToArray: %lld kB\n", vmPeak); /* larger example: 62646196 kB */ - if ((memLimit > 0) && ((vmPeak >> 20) > memLimit)) // the >> 20 converts kB to gB + if (threads > 1) { - threadCount = 1 + ((vmPeak >> 20) / memLimit); - verbose(1, "# over %d Gb at %lld kB, threadCount: %d\n", memLimit, vmPeak, threadCount); + int gB = vmPeak >> 20; /* convert kB to Gb */ + threadCount = threads; + verbose(1, "# at %d Gb (%lld kB), running %d threads\n", gB, vmPeak, threadCount); } if (queryGuides) // when range selected some query sequences { if (allGuides) // if there are any left on the all list { if (threadCount > 1) runThreads(threadCount, queryGuides, allGuides); else queryVsTarget(queryGuides, allGuides, 1, 0); } /* then run up the query vs. itself avoiding self vs. self */ queryVsSelf(queryGuides); countsOutput(queryGuides, bedFH); } else @@ -1345,40 +1329,35 @@ } // static void crisprKmers(char *sequence, FILE *bedFH) int main(int argc, char *argv[]) /* Process command line, initialize translation arrays and call the process */ { optionInit(&argc, argv, options); if (argc != 2) usage(); verbose(0, "# running verboseLevel: %d\n", verboseLevel()); bedFileOut = optionVal("bed", bedFileOut); dumpKmers = optionVal("dumpKmers", dumpKmers); loadKmers = optionVal("loadKmers", loadKmers); offTargets = optionVal("offTargets", offTargets); -memLimit = optionInt("memLimit", memLimit); -if (memLimit < 0) - errAbort("specified memLimi is less than 0: %d\n", memLimit); -if (memLimit > 0) - verbose(1, "# memLimit %d gB\n", memLimit); -cacheCount = optionInt("cacheCount", cacheCount); -if (cacheCount < 0) - errAbort("specified cacheCount is less than 0: %d\n", cacheCount); -if (cacheCount > 0) - verbose(1, "# cacheCount %d items\n", cacheCount); +threads = optionInt("threads", threads); +if (threads < 0) + errAbort("specified threads is less than 0: %d\n", threads); +if (threads > 0) + verbose(1, "# requesting threads: %d\n", threads); ranges = optionVal("ranges", ranges); if (ranges) rangesHash = readRanges(ranges); FILE *bedFH = NULL; if (bedFileOut) { bedFH = mustOpen(bedFileOut, "w"); verbose(1, "# output to bed file: %s\n", bedFileOut); } initOrderedNtVal(); /* set up orderedNtVal[] */ crisprKmers(argv[1], bedFH); if (verboseLevel() > 1)