Skip to content

Commit

Permalink
Simplify numOverlapsPerRead(). Issue marbl#619.
Browse files Browse the repository at this point in the history
  • Loading branch information
brianwalenz committed Sep 11, 2017
1 parent c8f3151 commit 651b37b
Show file tree
Hide file tree
Showing 3 changed files with 21 additions and 39 deletions.
53 changes: 20 additions & 33 deletions src/stores/ovStore.C
Original file line number Diff line number Diff line change
Expand Up @@ -487,51 +487,38 @@ ovStore::numOverlapsInRange(void) {



// Fills olapsPerRead with the number of overlaps per read.
// olapsPerRead must be at least 'lastRead+1' in size, not just lastRead-firstRead+1.
// olapsPerRead is not initialized; must be cleared before calling.
// Returns array with the number of overlaps per read. If numReads is more than zero,
// only those reads will be loaded.
//
void
ovStore::numOverlapsPerRead(uint32 *olapsPerRead,
uint32 firstRead,
uint32 lastRead) {

off_t originalPosition = AS_UTL_ftell(_offtFile);

AS_UTL_fseek(_offtFile, (off_t)firstRead * sizeof(ovStoreOfft), SEEK_SET);

uint64 len = lastRead - firstRead + 1;

ovStoreOfft *offsets = new ovStoreOfft [len];

uint64 act = AS_UTL_safeRead(_offtFile, offsets, "ovStore::numOverlapsPerRead", sizeof(ovStoreOfft), len);

if (len != act)
fprintf(stderr, "ovStore::numOverlapsPerRead()-- short read on offsets! Expected len=" F_U64 " read act=" F_U64 "\n", len, act), exit(1);
uint32 *
ovStore::numOverlapsPerRead(uint32 numReads) {

for (uint64 i=0; i<len; i++)
olapsPerRead[firstRead+i] = offsets[i]._numOlaps;
if ((numReads == 0) && (_gkp != NULL))
numReads = _gkp->gkStore_getNumReads();

delete [] offsets;
assert(numReads > 0);

AS_UTL_fseek(_offtFile, originalPosition, SEEK_SET);
}
numReads++; // Because read 0 doesn't exist.

uint32 *olapsPerRead = new uint32 [numReads];
ovStoreOfft *offsets = new ovStoreOfft [numReads];

off_t originalPosition = AS_UTL_ftell(_offtFile);

uint32 *
ovStore::numOverlapsPerRead(uint32 numReads) {
AS_UTL_fseek(_offtFile, 0, SEEK_SET);

if ((numReads == 0) && (_gkp != NULL))
numReads = _gkp->gkStore_getNumReads();
uint64 act = AS_UTL_safeRead(_offtFile, offsets, "ovStore::numOverlapsPerRead", sizeof(ovStoreOfft), numReads);

assert(numReads > 0);
if (numReads != act)
fprintf(stderr, "ovStore::numOverlapsPerRead()-- short read on offsets! Expected len=" F_U64 " read act=" F_U64 "\n",
numReads, act), exit(1);

uint32 *olapsPerRead = new uint32 [numReads+1];
for (uint64 i=0; i<numReads; i++)
olapsPerRead[i] = offsets[i]._numOlaps;

memset(olapsPerRead, 0, (numReads+1) * sizeof(uint32));
delete [] offsets;

numOverlapsPerRead(olapsPerRead, 1, numReads);
AS_UTL_fseek(_offtFile, originalPosition, SEEK_SET);

return(olapsPerRead);
}
Expand Down
3 changes: 0 additions & 3 deletions src/stores/ovStore.H
Original file line number Diff line number Diff line change
Expand Up @@ -317,9 +317,6 @@ public:

uint64 numOverlapsInRange(void);

void numOverlapsPerRead(uint32 *olapsPerRead,
uint32 firstRead,
uint32 lastRead);
uint32 *numOverlapsPerRead(uint32 numReads=0);

// Add new evalues for reads between bgnID and endID. No checking of IDs is done, but the number
Expand Down
4 changes: 1 addition & 3 deletions src/stores/ovStoreDump.C
Original file line number Diff line number Diff line change
Expand Up @@ -276,9 +276,7 @@ dumpStore(ovStore *ovlStore,
// bgnID, so we need to rewrite everything.

if ((asCounts) && (dumpType == 0)) {
counts = new uint32 [endID - bgnID + 1];

ovlStore->numOverlapsPerRead(counts, bgnID, endID);
counts = ovlStore->numOverlapsPerRead(endID);
ovlStore->setRange(1, 0);

for (uint32 ii=bgnID; ii<=endID; ii++)
Expand Down

0 comments on commit 651b37b

Please sign in to comment.