/* Function:  esl_workqueue_Reset()
 * Synopsis:  Reset the queue for another run.
 * Incept:    MSF, Thu Jun 18 11:51:39 2009
 *
 * Purpose:   Reset the queue for another run.  This is done by moving
 *            all the queued object to the reader's list (i.e. producer).
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslESYS> if thread synchronization fails somewhere.
 *            <eslEINVAL> if something's wrong with <queue>.
 */
int 
esl_workqueue_Reset(ESL_WORK_QUEUE *queue)
{
  int inx;
  int queueSize;

  if (queue == NULL)                                ESL_EXCEPTION(eslEINVAL, "Invalid queue object");
  if (pthread_mutex_lock (&queue->queueMutex) != 0) ESL_EXCEPTION(eslESYS,   "mutex lock failed");

  queueSize = queue->queueSize;

  /* move all buffers back to the reader queue */
  while (queue->workerQueueCnt > 0) 
    {
      inx = (queue->readerQueueHead + queue->readerQueueCnt) % queueSize;
      queue->readerQueue[inx] = queue->workerQueue[queue->workerQueueHead];
      ++queue->readerQueueCnt;

      queue->workerQueue[queue->workerQueueHead] = NULL;
      queue->workerQueueHead = (queue->workerQueueHead + 1) % queueSize;
      --queue->workerQueueCnt;
    }

  queue->pendingWorkers = 0;

  if (pthread_mutex_unlock (&queue->queueMutex) != 0) ESL_EXCEPTION(eslESYS, "mutex unlock failed");

  return eslOK;
}

/* Function:  esl_rmx_ValidateQ()
 * Incept:    SRE, Sun Mar 11 10:30:50 2007 [Janelia]
 *
 * Purpose:   Validates an instantaneous rate matrix <Q> for a
 *            continuous-time Markov process, whose elements $q_{ij}$
 *            represent instantaneous transition rates $i \rightarrow
 *            j$. 
 *            
 *            Rows satisfy the condition that
 *            $q_{ii} = -\sum_{i \neq j} q_{ij}$, and also
 *            that $q_{ij} \geq 0$ for all $j \neq i$. 
 *            
 *            <tol> specifies the floating-point tolerance to which
 *            that condition must hold: <fabs(sum-q_ii) <= tol>.
 *            
 *            <errbuf> is an optional error message buffer. The caller
 *            may pass <NULL> or a pointer to a buffer of at least
 *            <eslERRBUFSIZE> characters.
 *            
 * Args:      Q      - rate matrix to validate
 *            tol    - floating-point tolerance (0.00001, for example)      
 *            errbuf - OPTIONAL: ptr to an error buffer of at least
 *                     <eslERRBUFSIZE> characters.
 *
 * Returns:   <eslOK> on successful validation. 
 *            <eslFAIL> on failure, and if a non-<NULL> <errbuf> was
 *            provided by the caller, a message describing
 *            the reason for the failure is put there.
 *
 * Throws:    (no abnormal error conditions)
 */
int
esl_rmx_ValidateQ(ESL_DMATRIX *Q, double tol, char *errbuf)
{
  int    i,j;
  double qi;

  if (Q->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "Q must be type eslGENERAL to be validated");
  if (Q->n    != Q->m)       ESL_EXCEPTION(eslEINVAL, "a rate matrix Q must be square");

  for (i = 0; i < Q->n; i++)
    {
      qi = 0.;
      for (j = 0; j < Q->m; j++)
	{
	  if (i != j) {
	    if (Q->mx[i][j] < 0.)       ESL_FAIL(eslFAIL, errbuf, "offdiag elem %d,%d < 0",i,j);
	    qi += Q->mx[i][j];
	  } else {
	    if (Q->mx[i][j] > 0.)       ESL_FAIL(eslFAIL, errbuf, "diag elem %d,%d < 0", i,j);
	  }
	}
      if (fabs(qi + Q->mx[i][i]) > tol) ESL_FAIL(eslFAIL, errbuf, "row %d does not sum to 0.0", i);
    }
  return eslOK;
}

/* Function:  esl_dmatrix_Copy()
 *
 * Purpose:   Copies <src> matrix into <dest> matrix. <dest> must
 *            be allocated already by the caller.
 * 
 *            You may copy to a matrix of a different type, so long as
 *            the copy makes sense. If <dest> matrix is a packed type
 *            and <src> is not, the values that should be zeros must
 *            be zero in <src>, else the routine throws
 *            <eslEINCOMPAT>. If the <src> matrix is a packed type and
 *            <dest> is not, the values that are implicitly zeros are
 *            set to zeros in the <dest> matrix.
 *            
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEINCOMPAT> if <src>, <dest> are different sizes,
 *            or if their types differ and <dest> cannot represent
 *            <src>.
 */
int
esl_dmatrix_Copy(const ESL_DMATRIX *src, ESL_DMATRIX *dest)
{
  int i,j;

  if (dest->n != src->n || dest->m != src->m)
    ESL_EXCEPTION(eslEINCOMPAT, "matrices of different size");

  if (src->type == dest->type)   /* simple case. */
    memcpy(dest->mx[0], src->mx[0], src->ncells * sizeof(double));

  else if (src->type == eslGENERAL && dest->type == eslUPPER)		
    {
      for (i = 1; i < src->n; i++)
	for (j = 0; j < i; j++)
	  if (src->mx[i][j] != 0.) 
	    ESL_EXCEPTION(eslEINCOMPAT, "general matrix isn't upper triangular, can't be copied/packed");
      for (i = 0; i < src->n; i++)
	for (j = i; j < src->m; j++)
	  dest->mx[i][j] = src->mx[i][j];
    }
  
  else if (src->type == eslUPPER && dest->type == eslGENERAL)		
    {
      for (i = 1; i < src->n; i++)
	for (j = 0; j < i; j++)
	  dest->mx[i][j] = 0.;
      for (i = 0; i < src->n; i++)
	for (j = i; j < src->m; j++)
	  dest->mx[i][j] = src->mx[i][j];      
    }

  return eslOK;
}

/* Using FChoose() here would mean allocating tmp space for 2M-1 paths;
 * instead we use the fact that E(i) is itself the necessary normalization
 * factor, and implement FChoose's algorithm here for an on-the-fly 
 * calculation.
 */
static inline int
select_e(ESL_RANDOMNESS *rng, const P7_OPROFILE *om, const P7_OMX *ox, int i, int *ret_k)
{
  int    Q     = p7O_NQF(ox->M);
  float  sum   = 0.0;
  float  roll  = esl_random(rng);
  float  norm  = 1.0 / ox->xmx[i*p7X_NXCELLS+p7X_E]; /* all M, D already scaled exactly the same */
  __m128 xEv   = _mm_set1_ps(norm);
  union { __m128 v; float p[4]; } u;
  int    q,r;

  while (1) {
    for (q = 0; q < Q; q++)
      {
	u.v = _mm_mul_ps(ox->dpf[i][q*3 + p7X_M], xEv);
	for (r = 0; r < 4; r++) {
	  sum += u.p[r];
	  if (roll < sum) { *ret_k = r*Q + q + 1; return p7T_M;}
	}

	u.v = _mm_mul_ps(ox->dpf[i][q*3 + p7X_D], xEv);
	for (r = 0; r < 4; r++) {
	  sum += u.p[r];
	  if (roll < sum) { *ret_k = r*Q + q + 1; return p7T_D;}
	}
      }
    if (sum < 0.99) 
      ESL_EXCEPTION(-1, "Probabilities weren't normalized - failed to trace back from an E");
  }
  /*UNREACHED*/
  ESL_EXCEPTION(-1, "unreached code was reached. universe collapses.");
} 

/* Function:  esl_recorder_ResizeTo()
 * Synopsis:  Reallocate an <ESL_RECORDER> for a new <maxlines>
 * Incept:    SRE, Fri Dec 25 17:02:46 2009 [Casa de Gatos]
 *
 * Purpose:   Reallocate the <ESL_RECORDER> <rc> to have a new
 *            window size <maxlines>. 
 *            
 *            The new <maxlines> may be more or less than the previous
 *            window size for <rc>.
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEMEM> if (re-)allocation fails.
 *
 *            <eslEINVAL> if the recorder has a marked line (for start
 *            of a block) and you try to shrink it so much that that
 *            marked line would be lost.
 *
 *            <eslEINCONCEIVABLE> on any baseline resetting problem;
 *            this would have to be an internal error in the module.
 * 
 * Note:      We may have to repermute the line array, and reset its
 *            baseline, as follows.
 *            
 *            In the growth case: if the line array is out of order
 *            (circularly permuted) we must straighten it out, which
 *            means resetting the baseline.  
 *            i.e. to grow 3 1 2 to nalloc=6, we need 1 2 3 x x x; 
 *            simple reallocation to 3 1 2 x x x doesn't work,
 *            next read would make 3 4 2 x x x.
 *          
 *            In the shrinkage case: if the line array is in use beyond the
 *            new array size, we set a new baseline to keep as much of the
 *            old array as possible.
 * 
 *            i.e. for 6->3  
 *            1 2 3 x x x -> 1 2 3
 *            1 2 3 4 x x -> 2 3 4 with new baseline=2.
 *            4 5 0 1 2 3 -> 3 4 5 with new baseline=3
 */
int
esl_recorder_ResizeTo(ESL_RECORDER *rc, int new_maxlines)
{
  int   idx;
  int   newbase;
  void *tmp;
  int   minlines;
  int   status;

  if (new_maxlines == rc->nalloc) return eslOK;

  if (new_maxlines > rc->nalloc) /* growth case */
    {
      if ((rc->nread - rc->baseline) / rc->nalloc != 0)	/* array is permuted; reorder it */
	{
	  newbase = ESL_MAX(rc->baseline, rc->nread - rc->nalloc);
	  status = recorder_new_baseline(rc, newbase);
	  if (status) ESL_EXCEPTION(eslEINCONCEIVABLE, "baseline reset failed unexpectedly");
	}
    }
  else 				/* shrinkage case */
    {
      /* check that the marked line (if any) will stay in window */
      if (rc->markline >= 0)	
	{
	  minlines = rc->nread - rc->markline;
	  if (new_maxlines < minlines)
	    ESL_EXCEPTION(eslEINVAL, "can't shrink that far without losing marked line");
	}
      /* check that current line will stay in window */
      minlines = rc->nread - rc->ncurr + 1;
      if (new_maxlines < minlines)
	ESL_EXCEPTION(eslEINVAL, "can't shrink that far without losing current line");

      if (rc->nread - rc->baseline > new_maxlines) /* baseline needs to move up */
	{
	  newbase = rc->nread - new_maxlines;
	  status = recorder_new_baseline(rc, newbase);
	  if (status) ESL_EXCEPTION(eslEINCONCEIVABLE, "baseline reset failed unexpectedly");
	}

      for (idx = new_maxlines; idx < rc->nalloc; idx++)
	if (rc->line[idx]) free(rc->line[idx]);
    }

  ESL_RALLOC(rc->line,   tmp, sizeof(char *) * new_maxlines);
  ESL_RALLOC(rc->lalloc, tmp, sizeof(int)    * new_maxlines);
  ESL_RALLOC(rc->offset, tmp, sizeof(off_t)  * new_maxlines);
  for (idx = rc->nalloc; idx < new_maxlines; idx++) /* no-op in shrinkage case */
    { 
      rc->line[idx]   = NULL;
      rc->lalloc[idx] = 0;
      rc->offset[idx] = 0;
    }
  rc->nalloc = new_maxlines;
  return eslOK;

 ERROR:
  return status;
}

/* Function:  p7_hmm_mpi_Send()
 * Synopsis:  Send an HMM as an MPI work unit.
 *
 * Purpose:   Sends an HMM <hmm> as a work unit to MPI process
 *            <dest> (where <dest> ranges from 0..<nproc-1>), tagged
 *            with MPI tag <tag>, for MPI communicator <comm>, as 
 *            the sole workunit or result. 
 *            
 *            Work units are prefixed by a status code indicating the
 *            number of HMMs sent. If <hmm> is <NULL>, this code is 0,
 *            and <_Recv()> interprets such a unit as an EOD
 *            (end-of-data) signal, a signal to cleanly shut down
 *            worker processes.
 *            
 *            In order to minimize alloc/free cycles in this routine,
 *            caller passes a pointer to a working buffer <*buf> of
 *            size <*nalloc> characters. If necessary (i.e. if <hmm> is
 *            too big to fit), <*buf> will be reallocated and <*nalloc>
 *            increased to the new size. As a special case, if <*buf>
 *            is <NULL> and <*nalloc> is 0, the buffer will be
 *            allocated appropriately, but the caller is still
 *            responsible for free'ing it.
 *            
 * Returns:   <eslOK> on success; <*buf> may have been reallocated and
 *            <*nalloc> may have been increased.
 * 
 * Throws:    <eslESYS> if an MPI call fails; <eslEMEM> if a malloc/realloc
 *            fails. In either case, <*buf> and <*nalloc> remain valid and useful
 *            memory (though the contents of <*buf> are undefined). 
 * 
 * Note:      Compare to p7_hmmfile_WriteBinary(). The two operations (sending
 *            an HMM via MPI, or saving it as a binary file to disk) are
 *            similar.
 */
int
p7_hmm_mpi_Send(const P7_HMM *hmm, int dest, int tag, MPI_Comm comm, char **buf, int *nalloc)
{
  int   n = 0;
  int   code;
  int   sz, pos;
  int   status;

  /* Figure out size. We always send at least a status code (0=EOD=nothing sent) */
  if ( MPI_Pack_size(1, MPI_INT, comm, &sz)          != MPI_SUCCESS) ESL_EXCEPTION(eslESYS, "mpi pack size failed");  n += sz;
  if ((status = p7_hmm_mpi_PackSize(hmm, comm, &sz)) != eslOK)       return status;                                   n += sz;

  /* Make sure the buffer is allocated appropriately */
  if (*buf == NULL || n > *nalloc) 
    {
      ESL_REALLOC(*buf, sizeof(char) * n);
      *nalloc = n; 
    }

  /* Pack the status code and HMM into the buffer */
  /* The status code is the # of HMMs being sent as one MPI message; here 1 or 0 */
  pos  = 0;
  code = (hmm ? 1 : 0);
  if (MPI_Pack(&code, 1, MPI_INT,           *buf, n, &pos, comm)  != MPI_SUCCESS) ESL_EXCEPTION(eslESYS, "mpi pack failed");
  if (hmm && (status = p7_hmm_mpi_Pack(hmm, *buf, n, &pos, comm)) != eslOK)       return status;
  
  /* Send the packed HMM to the destination. */
  if (MPI_Send(*buf, n, MPI_PACKED, dest, tag, comm) != MPI_SUCCESS)  ESL_EXCEPTION(eslESYS, "mpi send failed");
  return eslOK;

 ERROR:
  return status;
}

/* Function:  p7_hmm_mpi_Recv()
 * Synopsis:  Receives an HMM as a work unit from an MPI sender.
 *
 * Purpose:   Receive a work unit that consists of a single HMM
 *            sent by MPI <source> (<0..nproc-1>, or
 *            <MPI_ANY_SOURCE>) tagged as <tag> for MPI communicator <comm>.
 *            
 *            Work units are prefixed by a status code that gives the
 *            number of HMMs to follow; here, 0 or 1 (but in the future,
 *            we could easily extend to sending several HMMs in one 
 *            packed buffer). If we receive a 1 code and we successfully
 *            unpack an HMM, this routine will return <eslOK> and a non-<NULL> <*ret_hmm>.
 *            If we receive a 0 code (a shutdown signal), 
 *            this routine returns <eslEOD> and <*ret_hmm> is <NULL>.
 *   
 *            Caller provides a working buffer <*buf> of size
 *            <*nalloc> characters. These are passed by reference, so
 *            that <*buf> can be reallocated and <*nalloc> increased
 *            if necessary. As a special case, if <*buf> is <NULL> and
 *            <*nalloc> is 0, the buffer will be allocated
 *            appropriately, but the caller is still responsible for
 *            free'ing it.
 *            
 *            Caller may or may not already know what alphabet the HMM
 *            is expected to be in.  A reference to the current
 *            alphabet is passed in <byp_abc>. If the alphabet is unknown,
 *            pass <*byp_abc = NULL>, and when the HMM is received, an
 *            appropriate new alphabet object is allocated and passed
 *            back to the caller via <*abc>.  If the alphabet is
 *            already known, <*byp_abc> is that alphabet, and the new
 *            HMM's alphabet type is verified to agree with it. This
 *            mechanism allows an application to let the first HMM
 *            determine the alphabet type for the application, while
 *            still keeping the alphabet under the application's scope
 *            of control.
 *
 * Args:      source  - index of MPI sender, 0..nproc-1 (0=master), or MPI_ANY_SOURCE
 *            tag     - MPI message tag;  MPI_ANY_TAG, or a specific message tag (0..32767 will work on any MPI)
 *            comm    - MPI communicator; MPI_COMM_WORLD, or a specific MPI communicator
 *            buf     - working buffer (for receiving packed message);
 *                      if <*buf> == NULL, a <*buf> is allocated and returned;
 *                      if <*buf> != NULL, it is used (and may be reallocated)
 *            nalloc  - allocation size of <*buf> in bytes; pass 0 if <*buf==NULL>.           
 *            byp_abc - BYPASS: <*byp_abc> == ESL_ALPHABET *> if known;
 *                              <*byp_abc> == NULL> if alphabet unknown.
 *            ret_hmm  - RETURN: newly allocated/received profile
 *
 * Returns:   <eslOK> on success. <*ret_hmm> contains the received HMM;
 *            it is allocated here, and the caller is responsible for
 *            free'ing it.  <*buf> may have been reallocated to a
 *            larger size, and <*nalloc> may have been increased.  If
 *            <*abc> was passed as <NULL>, it now points to an
 *            <ESL_ALPHABET> object that was allocated here; caller is
 *            responsible for free'ing this.
 *            
 *            Returns <eslEOD> if an end-of-data signal was received.
 *            In this case, <*buf>, <*nalloc>, and <*abc> are left unchanged,
 *            and <*ret_hmm> is <NULL>.
 *            
 *            Returns <eslEINCOMPAT> if the HMM is in a different alphabet
 *            than <*abc> said to expect. In this case, <*abc> is unchanged,
 *            <*buf> and <*nalloc> may have been changed, and <*ret_hmm> is
 *            <NULL>.
 *            
 * Throws:    <eslEMEM> on allocation error, and <eslESYS> on MPI communication
 *            errors; in either case <*ret_hmm> is <NULL>.           
 */
int
p7_hmm_mpi_Recv(int source, int tag, MPI_Comm comm, char **buf, int *nalloc, ESL_ALPHABET **byp_abc, P7_HMM **ret_hmm)
{
  int         pos = 0;
  int         code;
  int         n;
  MPI_Status  mpistatus;
  int         status;

  /* Probe first, because we need to know if our buffer is big enough. */
  if ( MPI_Probe(source, tag, comm, &mpistatus)  != MPI_SUCCESS) ESL_EXCEPTION(eslESYS, "mpi probe failed");
  if ( MPI_Get_count(&mpistatus, MPI_PACKED, &n) != MPI_SUCCESS) ESL_EXCEPTION(eslESYS, "mpi get count failed");

  /* Make sure the buffer is allocated appropriately */
  if (*buf == NULL || n > *nalloc) 
    {
      ESL_REALLOC(*buf, sizeof(char) * n);
      *nalloc = n; 
    }

  /* Receive the entire packed work unit */
  if (MPI_Recv(*buf, n, MPI_PACKED, source, tag, comm, &mpistatus) != MPI_SUCCESS) ESL_EXCEPTION(eslESYS, "mpi recv failed");

  /* Unpack the status code prefix */
  if (MPI_Unpack(*buf, n, &pos, &code, 1, MPI_INT, comm) != MPI_SUCCESS) ESL_EXCEPTION(eslESYS, "mpi unpack failed");

  if      (code == 0) { status = eslEOD; *ret_hmm = NULL; }
  else if (code == 1)   status = p7_hmm_mpi_Unpack(*buf, *nalloc, &pos, comm, byp_abc, ret_hmm);
  else                  ESL_EXCEPTION(eslESYS, "bad mpi buffer transmission code");
  return status;

 ERROR: /* from ESL_REALLOC only */
  *ret_hmm = NULL;
//.........这里部分代码省略.........

/* tau_by_moments()
 * 
 * Obtain an initial estimate for tau by 
 * matching moments. Also returns mean and
 * logsum, which we need for ML fitting.
 * To obtain a lambda estimate, use
 * lambda = tau / mean.
 */
static int
tau_by_moments(double *x, int n, double mu, double *ret_tau, double *ret_mean, double *ret_logsum)
{
  int    i;
  double mean, var, logsum;

  mean = var = logsum = 0.;
  for (i = 0; i < n; i++)
    {
      if (x[i] < mu) ESL_EXCEPTION(eslEINVAL, "No x[i] can be < mu in gamma data");
      mean   += x[i] - mu;	   /* mean is temporarily just the sum */
      logsum += log(x[i] - mu);
      var  += (x[i]-mu)*(x[i]-mu); /* var is temporarily the sum of squares */
    }
  var     = (var - mean*mean/(double)n) / ((double)n-1); /* now var is the variance */
  mean   /= (double) n;		/* and now mean is the mean */
  logsum /= (double) n;

  if (var == 0.)		/* and if mean = 0, var = 0 anyway. */
    ESL_EXCEPTION(eslEINVAL, "Zero variance in allegedly gamma-distributed dataset");
  
  if (ret_tau    != NULL) *ret_tau    = mean * mean / var;
  if (ret_mean   != NULL) *ret_mean   = mean;
  if (ret_logsum != NULL) *ret_logsum = logsum;
  return eslOK;
}

/* Function:  esl_workqueue_Remove()
 * Synopsis:  Removes a queued object from the producers list.
 * Incept:    MSF, Thu Jun 18 11:51:39 2009
 *
 * Purpose:   Removes a queued object from the producers list.
 *
 *            A object <void> that has already been consumed by a worker
 *            is removed the the producers list.  If there are no empty
 *            objects, a <obj> is set to NULL.
 *
 *            The pointer to the object is returned in the obj arguement.
 *
 * Returns:   <eslOK>  on success.
 *            <eslEOD> if no objects are in the queue.
 *
 * Throws:    <eslESYS> if thread synchronization fails somewhere.
 *            <eslEINVAL> if something's wrong with <queue>.
 */
int 
esl_workqueue_Remove(ESL_WORK_QUEUE *queue, void **obj)
{
  int inx;
  int status = eslEOD;

  if (obj == NULL)   ESL_EXCEPTION(eslEINVAL, "Invalid object pointer");
  if (queue == NULL) ESL_EXCEPTION(eslEINVAL, "Invalid queue object");

  if (pthread_mutex_lock (&queue->queueMutex) != 0) ESL_EXCEPTION(eslESYS, "mutex lock failed");

  /* check if there are any items on the readers list */
  *obj = NULL;
  if (queue->readerQueueCnt > 0)
    {
      inx = (queue->readerQueueHead + queue->readerQueueCnt) % queue->queueSize;
      *obj = queue->readerQueue[inx];
      queue->readerQueue[inx] = NULL;
      --queue->readerQueueCnt;
      status = eslOK;
    }

  if (pthread_mutex_unlock (&queue->queueMutex) != 0) ESL_EXCEPTION(eslESYS, "mutex unlock failed");

  return status;
}

/* Function:  esl_msashuffle_Shuffle()
 * Synopsis:  Shuffle an alignment's columns.
 *
 * Purpose:   Returns a column-shuffled version of <msa> in <shuf>,
 *            using random generator <r>. Shuffling by columns
 *            preserves the \% identity of the original
 *            alignment. <msa> and <shuf> can be identical, to shuffle
 *            in place.
 *            
 *            The caller sets up the rest of the data (everything but
 *            the alignment itself) in <shuf> the way it wants,
 *            including sequence names, MSA name, and other
 *            annotation. The easy thing to do is to make <shuf>
 *            a copy of <msa>: the caller might create <shuf> by
 *            a call to <esl_msa_Clone()>.
 *            
 *            The alignments <msa> and <shuf> can both be in digital
 *            mode, or can both be in text mode; you cannot mix
 *            digital and text modes.
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEINVAL> if <msa>,<shuf> aren't in the same mode (digital vs. text).
 */
int
esl_msashuffle_Shuffle(ESL_RANDOMNESS *r, ESL_MSA *msa, ESL_MSA *shuf)
{
  int i, pos, alen;

  if (! (msa->flags & eslMSA_DIGITAL))
    {
      char c;
      if (shuf->flags & eslMSA_DIGITAL) ESL_EXCEPTION(eslEINVAL, "<shuf> must be in text mode if <msa> is");
      if (msa != shuf) {
	for (i = 0; i < msa->nseq; i++)
	  strcpy(shuf->aseq[i], msa->aseq[i]);
      }

      for (i = 0; i < msa->nseq; i++)
	shuf->aseq[i][msa->alen] = '\0';

      for (alen = msa->alen; alen > 1; alen--)
	{
	  pos = esl_rnd_Roll(r, alen);
	  for (i = 0; i < msa->nseq; i++)
	    {
	      c                     = msa->aseq[i][pos];
	      shuf->aseq[i][pos]    = shuf->aseq[i][alen-1];
	      shuf->aseq[i][alen-1] = c;
	    }
	}
    }
#ifdef eslAUGMENT_ALPHABET
  else 
    {
      ESL_DSQ x;
      if (! (shuf->flags & eslMSA_DIGITAL)) ESL_EXCEPTION(eslEINVAL, "<shuf> must be in digital mode if <msa> is");

      if (msa != shuf) {
	for (i = 0; i < msa->nseq; i++)
	  memcpy(shuf->ax[i], msa->ax[i], (msa->alen + 2) * sizeof(ESL_DSQ));
      }

      for (i = 0; i < msa->nseq; i++)
	shuf->ax[i][msa->alen+1] = eslDSQ_SENTINEL;

      for (alen = msa->alen; alen > 1; alen--)
	{
	  pos = esl_rnd_Roll(r, alen) + 1;
	  for (i = 0; i < msa->nseq; i++)
	    {
	      x                 = msa->ax[i][pos];
	      shuf->ax[i][pos]  = shuf->ax[i][alen];
	      shuf->ax[i][alen] = x;
	    }
	}
    }
#endif /*eslAUGMENT_ALPHABET*/

  return eslOK;
}

/* Function:  esl_hyperexp_Write()
 *
 * Purpose:   Write hyperexponential parameters from <hxp> to an open <fp>.
 *            
 *            The output format is suitable for input by <esl_hyperexp_Read()>.
 *
 * Returns:   <eslOK> on success.
 * 
 * Throws:    <eslEWRITE> on any write error.
 */
int
esl_hyperexp_Write(FILE *fp, ESL_HYPEREXP *hxp)
{
  int k;

  if (fprintf(fp, "%8d     # number of components\n", hxp->K)     < 0) ESL_EXCEPTION(eslEWRITE, "hyperexp write failed");
  if (fprintf(fp, "%8.2f   # mu (for all components)\n", hxp->mu) < 0) ESL_EXCEPTION(eslEWRITE, "hyperexp write failed");
  for (k = 0; k < hxp->K; k++)
    if (fprintf(fp, "%8.6f %12.6f  # q[%d], lambda[%d]\n",
		hxp->q[k], hxp->lambda[k], k, k)                  < 0) ESL_EXCEPTION(eslEWRITE, "hyperexp write failed");
  return eslOK;
}

/* Function:  p7_Forward()
 * Synopsis:  The Forward algorithm, full matrix fill version.
 * Incept:    SRE, Fri Aug 15 18:59:43 2008 [Casa de Gatos]
 *
 * Purpose:   Calculates the Forward algorithm for sequence <dsq> of
 *            length <L> residues, using optimized profile <om>, and a
 *            preallocated DP matrix <ox>. Upon successful return, <ox>
 *            contains the filled Forward matrix, and <*opt_sc>
 *            optionally contains the raw Forward score in nats.
 *            
 *            This calculation requires $O(ML)$ memory and time.
 *            The caller must provide a suitably allocated full <ox>
 *            by calling <ox = p7_omx_Create(M, L, L)> or
 *            <p7_omx_GrowTo(ox, M, L, L)>.
 *
 *            The model <om> must be configured in local alignment
 *            mode. The sparse rescaling method used to keep
 *            probability values within single-precision floating
 *            point dynamic range cannot be safely applied to models in
 *            glocal or global mode.
 *
 * Args:      dsq     - digital target sequence, 1..L
 *            L       - length of dsq in residues          
 *            om      - optimized profile
 *            ox      - RETURN: Forward DP matrix
 *            opt_sc  - RETURN: Forward score (in nats)          
 *
 * Returns:   <eslOK> on success. 
 *
 * Throws:    <eslEINVAL> if <ox> allocation is too small, or if the profile
 *            isn't in local alignment mode.
 *            <eslERANGE> if the score exceeds the limited range of
 *            a probability-space odds ratio.
 *            In either case, <*opt_sc> is undefined.
 */
int
p7_Forward(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, P7_OMX *ox, float *opt_sc)
{
#ifdef p7_DEBUGGING		
  if (om->M >  ox->allocQ4*4)    ESL_EXCEPTION(eslEINVAL, "DP matrix allocated too small (too few columns)");
  if (L     >= ox->validR)       ESL_EXCEPTION(eslEINVAL, "DP matrix allocated too small (too few MDI rows)");
  if (L     >= ox->allocXR)      ESL_EXCEPTION(eslEINVAL, "DP matrix allocated too small (too few X rows)");
  if (! p7_oprofile_IsLocal(om)) ESL_EXCEPTION(eslEINVAL, "Forward implementation makes assumptions that only work for local alignment");
#endif

  return forward_engine(TRUE, dsq, L, om, ox, opt_sc);
}

/* Function:  p7_oprofile_ReadRest()
 * Synopsis:  Read the rest of an optimized profile.
 *
 * Purpose:   Read the rest of an optimized profile <om> from
 *            the <.h3p> file associated with an open HMM
 *            file <hfp>. 
 *            
 *            This is the second part of a two-part calling sequence.
 *            The <om> here must be the result of a previous
 *            successful <p7_oprofile_ReadMSV()> call on the same
 *            open <hfp>.
 *
 * Args:      hfp - open HMM file, from which we've previously
 *                  called <p7_oprofile_ReadMSV()>.
 *            om  - optimized profile that was successfully
 *                  returned by  <p7_oprofile_ReadMSV()>.
 *
 * Returns:   <eslOK> on success, and <om> is now a complete
 *            optimized profile.
 *            
 *            Returns <eslEFORMAT> if <hfp> has no <.h3p> file open,
 *            or on any parsing error, and set <hfp->errbuf> to
 *            an informative error message.
 *
 * Throws:    <eslESYS> if an <fseek()> fails to reposition the
 *            binary <.h3p> file.
 *            
 *            <eslEMEM> on allocation error.
 */
int
p7_oprofile_ReadRest(P7_HMMFILE *hfp, P7_OPROFILE *om)
{
  uint32_t       magic;
  int            M;
  int            alphatype;
  int            status = eslOK;

#ifdef HMMER_THREADS
  /* lock the mutex to prevent other threads from reading from the optimized
   * profile at the same time.
   */
  if (hfp->syncRead)
    {
      if (pthread_mutex_lock (&hfp->readMutex) != 0) ESL_EXCEPTION(eslESYS, "mutex lock failed");
    }
#endif

  if (! fread( (char *) &magic,          sizeof(uint32_t), 1,           hfp->pfp)) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "failed to read magic");
  if (magic == v3a_pmagic) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "binary auxfiles are in an outdated HMMER format (3/a); please hmmpress your HMM file again");
  if (magic == v3b_pmagic) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "binary auxfiles are in an outdated HMMER format (3/b); please hmmpress your HMM file again");
  if (magic == v3c_pmagic) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "binary auxfiles are in an outdated HMMER format (3/c); please hmmpress your HMM file again");
  if (magic == v3d_pmagic) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "binary auxfiles are in an outdated HMMER format (3/d); please hmmpress your HMM file again");
  if (magic == v3e_pmagic) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "binary auxfiles are in an outdated HMMER format (3/e); please hmmpress your HMM file again");
  if (magic != v3f_pmagic) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "bad magic; not an HMM database file?");

  if (! fread( (char *) &M,              sizeof(int),      1,           hfp->pfp)) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "failed to read model size M");
  if (! fread( (char *) &alphatype,      sizeof(int),      1,           hfp->pfp)) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "failed to read alphabet type");  
  if (! fread( (char *) &magic,          sizeof(uint32_t), 1,           hfp->pfp)) ESL_XFAIL(eslEFORMAT, hfp->errbuf, "no sentinel magic: .h3p file corrupted?");

  if (magic     != v3f_pmagic)     ESL_XFAIL(eslEFORMAT, hfp->errbuf, "bad sentinel magic: .h3p file corrupted?");
  if (M         != om->M)          ESL_XFAIL(eslEFORMAT, hfp->errbuf, "p/f model length mismatch");
  if (alphatype != om->abc->type)  ESL_XFAIL(eslEFORMAT, hfp->errbuf, "p/f alphabet type mismatch");

#ifdef HMMER_THREADS
  if (hfp->syncRead)
    {
      if (pthread_mutex_unlock (&hfp->readMutex) != 0) ESL_EXCEPTION(eslESYS, "mutex unlock failed");
    }
#endif

  return eslOK;

 ERROR:
#ifdef HMMER_THREADS
  if (hfp->syncRead)
    {
      if (pthread_mutex_unlock (&hfp->readMutex) != 0) ESL_EXCEPTION(eslESYS, "mutex unlock failed");
    }
#endif

  return status;
}

/* Function:  p7_oprofile_Position()
 * Synopsis:  Reposition an open hmm file to an offset.
 * Incept:    MSF, Thu Oct 15, 2009 [Janelia]
 *
 * Purpose:   Reposition an open <hfp> to offset <offset>.
 *            <offset> would usually be the first byte of a
 *            desired hmm record.
 *            
 * Returns:   <eslOK>     on success;
 *            <eslEOF>    if no data can be read from this position.
 *
 * Throws:    <eslEINVAL>  if the <sqfp> is not positionable.
 *            <eslEFORMAT> if no msv profile opened.
 *            <eslESYS>    if the fseeko() call fails.
 */
int
p7_oprofile_Position(P7_HMMFILE *hfp, off_t offset)
{
  if (hfp->ffp == NULL)  ESL_EXCEPTION(eslEFORMAT, hfp->errbuf, "no MSV profile file; hmmpress probably wasn't run");
  if (hfp->do_stdin)     ESL_EXCEPTION(eslEINVAL, "can't Position() in standard input");
  if (hfp->do_gzip)      ESL_EXCEPTION(eslEINVAL, "can't Position() in a gzipped file");
  if (offset < 0)        ESL_EXCEPTION(eslEINVAL, "bad offset");

  if (fseeko(hfp->ffp, offset, SEEK_SET) != 0) ESL_EXCEPTION(eslESYS, "fseeko() failed");

  return eslOK;
}

/* Function:  p7_BackwardParser()
 * Synopsis:  The Backward algorithm, linear memory parsing version.
 * Incept:    SRE, Sat Aug 16 08:34:13 2008 [Janelia]
 *
 * Purpose:   Same as <p7_Backward()> except that the full matrix isn't
 *            kept. Instead, a linear $O(M+L)$ memory algorithm is
 *            used, keeping only the DP matrix values for the special
 *            (BENCJ) states. These are sufficient to do posterior
 *            decoding to identify high-probability regions where
 *            domains are.
 *       
 *            The caller must provide a suitably allocated "parsing"
 *            <bck> by calling <bck = p7_omx_Create(M, 0, L)> or
 *            <p7_omx_GrowTo(bck, M, 0, L)>.
 *
 * Args:      dsq     - digital target sequence, 1..L
 *            L       - length of dsq in residues          
 *            om      - optimized profile
 *            fwd     - filled Forward DP matrix, for scale factors
 *            bck     - RETURN: filled Backward matrix
 *            opt_sc  - optRETURN: Backward score (in nats)          
 *
 * Returns:   <eslOK> on success. 
 *
 * Throws:    <eslEINVAL> if <ox> allocation is too small, or if the profile
 *            isn't in local alignment mode.
 *            <eslERANGE> if the score exceeds the limited range of
 *            a probability-space odds ratio.
 *            In either case, <*opt_sc> is undefined.
 */
int 
p7_BackwardParser(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, const P7_OMX *fwd, P7_OMX *bck, float *opt_sc)
{
#ifdef p7_DEBUGGING		
  if (om->M >  bck->allocQ4*4)    ESL_EXCEPTION(eslEINVAL, "DP matrix allocated too small (too few columns)");
  if (bck->validR < 1)            ESL_EXCEPTION(eslEINVAL, "DP matrix allocated too small (too few MDI rows)");
  if (L     >= bck->allocXR)      ESL_EXCEPTION(eslEINVAL, "DP matrix allocated too small (too few X rows)");
  if (L     != fwd->L)            ESL_EXCEPTION(eslEINVAL, "fwd matrix size doesn't agree with length L");
  if (! p7_oprofile_IsLocal(om))  ESL_EXCEPTION(eslEINVAL, "Forward implementation makes assumptions that only work for local alignment");
#endif

  return backward_engine(FALSE, dsq, L, om, fwd, bck, opt_sc);
}

/* Function:  esl_workqueue_Complete()
 * Synopsis:  Signals the end of the queue.
 * Incept:    MSF, Thu Jun 18 11:51:39 2009
 *
 * Purpose:   Signal the end of the queue.  If there are any threads
 *            waiting on an object, signal them to wake up and complete
 *            their processing.
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslESYS> if thread synchronization fails somewhere.
 *            <eslEINVAL> if something's wrong with <queue>.
 */
int 
esl_workqueue_Complete(ESL_WORK_QUEUE *queue)
{
  if (queue == NULL)                                ESL_EXCEPTION(eslEINVAL, "Invalid queue object");
  if (pthread_mutex_lock (&queue->queueMutex) != 0) ESL_EXCEPTION(eslESYS,   "mutex lock failed");

  if (queue->pendingWorkers != 0)
    {
      if (pthread_cond_broadcast (&queue->workerQueueCond) != 0) ESL_EXCEPTION(eslESYS, "broadcast failed");
    }

  if (pthread_mutex_unlock (&queue->queueMutex) != 0) ESL_EXCEPTION(eslESYS, "mutex unlock failed");

  return eslOK;
}

/* Function:  p7_emit_FancyConsensus()
 * Synopsis:  Emit a fancier consensus with upper/lower case and N/X's.
 * Incept:    SRE, Fri May 14 09:33:10 2010 [Janelia]
 *
 * Purpose:   Generate a consensus sequence for model <hmm>, consisting
 *            of the maximum probability residue in each match state;
 *            store this sequence in text-mode <sq> provided by the caller.
 *            
 *            If the probability of the consensus residue is less than
 *            <min_lower>, show an ``any'' residue (N or X) instead.
 *            If the probability of the consensus residue is $\geq$
 *            <min_lower>  and less than <min_upper>, show the residue
 *            as lower case; if it is $\geq$ <min_upper>, show it as
 *            upper case.
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEINVAL> if the <sq> isn't in text mode.
 *
 * Xref:      SRE:J6/59.
 */
int
p7_emit_FancyConsensus(const P7_HMM *hmm, float min_lower, float min_upper, ESL_SQ *sq)
{
  int   k, x;
  float p;
  char  c;
  int   status;

  if (! esl_sq_IsText(sq)) ESL_EXCEPTION(eslEINVAL, "p7_emit_FancyConsensus() expects a text-mode <sq>");

  if ((status = esl_sq_GrowTo(sq, hmm->M)) != eslOK) return status;

  for (k = 1; k <= hmm->M; k++)
  {

    if (hmm->mm && hmm->mm[k] == 'm') { //masked position, spit out the degenerate code
      if ((status = esl_sq_CAddResidue(sq, tolower(esl_abc_CGetUnknown(hmm->abc))) ) != eslOK) return status;
    } else {
      p = esl_vec_FMax(   hmm->mat[k], hmm->abc->K);
      x = esl_vec_FArgMax(hmm->mat[k], hmm->abc->K);
  
      if      (p <  min_lower)  c = tolower(esl_abc_CGetUnknown(hmm->abc));
      else if (p >= min_upper)  c = toupper(hmm->abc->sym[x]);
      else                      c = tolower(hmm->abc->sym[x]);

      if ((status = esl_sq_CAddResidue(sq, c)) != eslOK) return status;
    }
  }
  if ((status = esl_sq_CAddResidue(sq, '\0')) != eslOK) return status;
  return eslOK;
}