esl_gencode: compact block for translated ORFs

esl_gencode.c

@@ -726,7 +726,7 @@ esl_gencode_WorkstateDestroy(ESL_GENCODE_WORKSTATE *wrk)
 
       if(wrk->orf_block != NULL)
       {
-         esl_sq_DestroyBlock(wrk->orf_block);
+         esl_gencode_OrfBlockDestroy(wrk->orf_block);
          wrk->orf_block = NULL;
       }
 
@@ -779,11 +779,135 @@ esl_gencode_WorkstateCreate(ESL_GETOPTS *go, ESL_GENCODE *gcode)
  *  6. Functions for processing ORFs
  *****************************************************************/
 
+/* Residue storage of an ORF block comes in chunks of at least this many bytes.
+ * Kept small: when a block is created and freed for every query, 1 MB chunks
+ * raised peak memory about 30%, and 512 kB chunks sometimes did; 64 kB to
+ * 256 kB did not. */
+#define eslORF_CHUNK (1 << 16)
+
+/* Function:  esl_gencode_OrfBlockCreate()
+ * Synopsis:  Create an empty block of ORFs.
+ *
+ * Purpose:   Create a block with room for <count> ORFs; it grows as needed.
+ *
+ * Returns:   the new block; <NULL> on allocation failure.
+ */
+ESL_ORF_BLOCK *
+esl_gencode_OrfBlockCreate(int count)
+{
+  ESL_ORF_BLOCK *block = NULL;
+  int            status;
+
+  if (count < 1) count = 1;
+  ESL_ALLOC(block, sizeof(ESL_ORF_BLOCK));
+  block->list       = NULL;
+  block->chunk      = NULL;
+  block->chunkalloc = NULL;
+  block->count      = 0;
+  block->listSize   = count;
+  block->nchunks    = 0;
+  block->cur        = 0;
+  block->used       = 0;
+
+  ESL_ALLOC(block->list,       sizeof(ESL_ORF)   * count);
+  ESL_ALLOC(block->chunk,      sizeof(ESL_DSQ *));
+  ESL_ALLOC(block->chunkalloc, sizeof(int64_t));
+  block->chunk[0] = NULL;
+  ESL_ALLOC(block->chunk[0],   sizeof(ESL_DSQ) * eslORF_CHUNK);
+  block->chunkalloc[0] = eslORF_CHUNK;
+  block->nchunks       = 1;
+  return block;
+
+ ERROR:
+  esl_gencode_OrfBlockDestroy(block);
+  return NULL;
+}
+
+/* Function:  esl_gencode_OrfBlockReuse()
+ * Synopsis:  Empty a block of ORFs, keeping its memory.
+ */
+void
+esl_gencode_OrfBlockReuse(ESL_ORF_BLOCK *block)
+{
+  block->count = 0;
+  block->cur   = 0;
+  block->used  = 0;
+}
+
+/* Function:  esl_gencode_OrfBlockDestroy()
+ * Synopsis:  Free a block of ORFs.
+ */
+void
+esl_gencode_OrfBlockDestroy(ESL_ORF_BLOCK *block)
+{
+  int c;
+
+  if (block == NULL) return;
+  if (block->chunk != NULL)
+    for (c = 0; c < block->nchunks; c++) free(block->chunk[c]);
+  free(block->chunk);
+  free(block->chunkalloc);
+  free(block->list);
+  free(block);
+}
+
+/* orf_block_add()
+ * Add the finished ORF <psq> to <block>.
+ */
+static int
+orf_block_add(ESL_ORF_BLOCK *block, const ESL_SQ *psq)
+{
+  ESL_ORF *orf;
+  int64_t  nbytes = psq->n + 2;   /* residues and the two sentinels */
+  int      status;
+
+  if (block->count == block->listSize)
+    {
+      ESL_REALLOC(block->list, sizeof(ESL_ORF) * block->listSize * 2);
+      block->listSize *= 2;
+    }
+
+  if (block->used + nbytes > block->chunkalloc[block->cur])
+    { /* move to the next chunk. Chunks after <cur> hold no ORFs, so one of them can be created or resized */
+      int next = (block->used > 0 ? block->cur + 1 : block->cur);
+      if (next == block->nchunks)
+        {
+          ESL_REALLOC(block->chunk,      sizeof(ESL_DSQ *) * (block->nchunks + 1));
+          ESL_REALLOC(block->chunkalloc, sizeof(int64_t)   * (block->nchunks + 1));
+          block->chunk[next]      = NULL;
+          block->chunkalloc[next] = 0;
+          block->nchunks++;
+        }
+      if (block->chunkalloc[next] < nbytes)
+        {
+          int64_t newalloc = ESL_MAX(nbytes, eslORF_CHUNK);
+          ESL_REALLOC(block->chunk[next], sizeof(ESL_DSQ) * newalloc);
+          block->chunkalloc[next] = newalloc;
+        }
+      block->cur  = next;
+      block->used = 0;
+    }
+
+  orf        = block->list + block->count;
+  orf->dsq   = block->chunk[block->cur] + block->used;
+  memcpy(orf->dsq, psq->dsq, sizeof(ESL_DSQ) * nbytes);
+  orf->n     = psq->n;
+  orf->start = psq->start;
+  orf->end   = psq->end;
+  block->used += nbytes;
+  block->count++;
+  return eslOK;
+
+ ERROR:
+  return status;
+}
+
 int
 esl_gencode_ProcessOrf(ESL_GENCODE_WORKSTATE *wrk, ESL_SQ *sq)
 {
 
   int              status   = eslOK;
+  int              named    = FALSE;   /* TRUE if <psq> was given a name and description for output */
   ESL_SQ *psq = wrk->psq[wrk->frame];
   psq->end = (wrk->is_revcomp ? wrk->apos+1 : wrk->apos-1);
   if (wrk->in_orf[wrk->frame] && psq->n >= wrk->minlen)
@@ -796,6 +920,7 @@ esl_gencode_ProcessOrf(ESL_GENCODE_WORKSTATE *wrk, ESL_SQ *sq)
       if (wrk->orf_block == NULL) {  /* esl-translate: name and description go to FASTA output */
         esl_sq_FormatName(psq, "orf%d", wrk->orfcount);
         esl_sq_FormatDesc(psq, "source=%s coords=%" PRId64 "..%" PRId64 " length=%" PRId64 " frame=%d desc=%s", psq->source, psq->start, psq->end, psq->n, wrk->frame + 1 + (wrk->is_revcomp ? 3 : 0), sq->desc);
+        named = TRUE;
       }
       /* if we do not have a block to write ORFs to then write ORFs to file */
       if (wrk->orf_block == NULL)
@@ -804,20 +929,26 @@ esl_gencode_ProcessOrf(ESL_GENCODE_WORKSTATE *wrk, ESL_SQ *sq)
       }
       else
       {
-        if (wrk->orf_block->count == wrk->orf_block->listSize)
-        {
-          status = esl_sq_BlockGrowTo(wrk->orf_block, wrk->orf_block->listSize + 128, TRUE, psq->abc);
-          if (status != eslOK) ESL_XEXCEPTION(eslEMEM, "Cannot increase size of ORF sequence block");
-        }
-
-        esl_sq_Copy(psq, &(wrk->orf_block->list[wrk->orf_block->count]));
-
-        wrk->orf_block->count++;
+        if ((status = orf_block_add(wrk->orf_block, psq)) != eslOK) ESL_XEXCEPTION(eslEMEM, "Cannot add to the ORF block");
       }
     }
 
-  esl_sq_Reuse(psq);
-  esl_sq_SetSource(psq, sq->name);
+  /* Ready <psq> for the next ORF. Unless it was named for output above, only
+   * its residues and coordinates have changed since it was last reset, and
+   * its source is still the one esl_gencode_ProcessStart() set.
+   */
+  if (named)
+    {
+      esl_sq_Reuse(psq);
+      esl_sq_SetSource(psq, sq->name);
+    }
+  else
+    {
+      psq->dsq[1] = eslDSQ_SENTINEL;
+      psq->n      = 0;
+      psq->start  = 0;
+      psq->end    = 0;
+    }
   wrk->in_orf[wrk->frame] = FALSE;
 
  ERROR:

esl_gencode.h

@@ -22,6 +22,35 @@ typedef struct {
   const ESL_ALPHABET *aa_abc;  // A reference to amino alphabet that caller is maintaining 
 } ESL_GENCODE;
 
+/* A translated ORF in an <ESL_ORF_BLOCK>: its residues, and where it lies on
+ * the nucleotide sequence it came from.
+ */
+typedef struct {
+  ESL_DSQ *dsq;     // residues 1..n, with sentinel bytes at 0 and n+1; points into the block's residue storage
+  int64_t  n;       // length in residues
+  int64_t  start;   // nucleotide coord of the ORF's first base
+  int64_t  end;     // nucleotide coord of its last base; < start on the reverse strand
+} ESL_ORF;
+
+/* The ORFs translated from one or more nucleotide sequences. Their residues
+ * are kept end to end in a few large chunks, so an ORF costs its residues
+ * and one <ESL_ORF>.
+ */
+typedef struct {
+  ESL_ORF  *list;        // ORFs 0..count-1
+  int       count;       // number of ORFs
+  int       listSize;    // number of ORFs <list> is allocated for
+  ESL_DSQ **chunk;       // residue storage, chunks 0..nchunks-1; a chunk doesn't move while it holds ORFs
+  int64_t  *chunkalloc;  // allocated size of each chunk, in bytes
+  int       nchunks;     // number of chunks
+  int       cur;         // chunk being filled
+  int64_t   used;        // bytes used in chunk <cur>
+} ESL_ORF_BLOCK;
+
+extern ESL_ORF_BLOCK *esl_gencode_OrfBlockCreate (int count);
+extern void           esl_gencode_OrfBlockReuse  (ESL_ORF_BLOCK *block);
+extern void           esl_gencode_OrfBlockDestroy(ESL_ORF_BLOCK *block);
+
 /* struct esl_gencode_workstate_s
  *   keeps state in DNA sequence <sq>, allowing us to process a sequence
  *   either in a single gulp (using ReadSeq) or in overlapping windows
@@ -42,7 +71,7 @@ typedef struct esl_gencode_workstate_s {
   int     is_revcomp;   // TRUE|FALSE: TRUE if we're doing reverse complement strand
   int     orfcount;     // >=0:   How many ORFs we've processed so far
 
-  ESL_SQ_BLOCK  *orf_block; // sequences block to hold ORFs
+  ESL_ORF_BLOCK *orf_block; // block to hold ORFs
 
   /* one-time configuration information (from options) */
   int     do_watson;         // TRUE|FALSE:  TRUE if we translate the top strand