/* keyhash_create()
 * 
 * The real creation function, which takes arguments for memory sizes.
 * This is abstracted to a static function because it's used by both
 * Create() and Clone() but slightly differently.
 *
 * Args:  hashsize          - size of hash table; this must be a power of two.
 *        init_key_alloc    - initial allocation for # of keys.
 *        init_string_alloc - initial allocation for total size of key strings.
 *
 * Returns:  An allocated hash table structure; or NULL on failure.
 */
ESL_KEYHASH *
keyhash_create(uint32_t hashsize, int init_key_alloc, int init_string_alloc)
{
  ESL_KEYHASH *kh = NULL;
  int  i;
  int  status;

  ESL_ALLOC(kh, sizeof(ESL_KEYHASH));
  kh->hashtable  = NULL;
  kh->key_offset = NULL;
  kh->nxt        = NULL;
  kh->smem       = NULL;

  kh->hashsize  = hashsize;
  kh->kalloc    = init_key_alloc;
  kh->salloc    = init_string_alloc;

  ESL_ALLOC(kh->hashtable, sizeof(int) * kh->hashsize);
  for (i = 0; i < kh->hashsize; i++)  kh->hashtable[i] = -1;

  ESL_ALLOC(kh->key_offset, sizeof(int) * kh->kalloc);
  ESL_ALLOC(kh->nxt,        sizeof(int) * kh->kalloc);
  for (i = 0; i < kh->kalloc; i++)  kh->nxt[i] = -1;

  ESL_ALLOC(kh->smem,   sizeof(char) * kh->salloc);
  kh->nkeys = 0;
  kh->sn    = 0;
  return kh;

 ERROR:
  esl_keyhash_Destroy(kh);
  return NULL;
}

/* Function:  esl_dmatrix_Create()
 *
 * Purpose:   Creates a general <n> x <m> matrix (<n> rows, <m> 
 *            columns).
 *
 * Args:      <n> - number of rows;    $>= 1$
 *            <m> - number of columns; $>= 1$
 * 
 * Returns:   a pointer to a new <ESL_DMATRIX> object. Caller frees
 *            with <esl_dmatrix_Destroy()>.
 *
 * Throws:    <NULL> if an allocation failed.
 */
ESL_DMATRIX *
esl_dmatrix_Create(int n, int m)
{
  ESL_DMATRIX *A = NULL;
  int r;
  int status;

  ESL_ALLOC(A, sizeof(ESL_DMATRIX));
  A->mx = NULL;
  A->n  = n;
  A->m  = m;

  ESL_ALLOC(A->mx, sizeof(double *) * n);
  A->mx[0] = NULL;

  ESL_ALLOC(A->mx[0], sizeof(double) * n * m);
  for (r = 1; r < n; r++)
    A->mx[r] = A->mx[0] + r*m;

  A->type   = eslGENERAL;
  A->ncells = n * m; 
  return A;
  
 ERROR:
  esl_dmatrix_Destroy(A);
  return NULL;
}

/* Function:  p7_coords2_hash_Create()
 * Synopsis:  Create a <P7_COORDS2_HASH>
 *
 * Purpose:   Allocate and initialize a <P7_COORDS2_HASH> hash table for storing
 *            lots of coord2 arrays (i.e. domain annotations).
 * 
 *            The <init_*> arguments let you set non-default initial
 *            allocation sizes. To use the default for any of these,
 *            pass a 0 value. Defaults are 128 for the initial 
 *            hashtable size <init_hashsize>; 128 for the initial
 *            allocation for number of keys to be stored <init_nkeyalloc>;
 *            and 2048 for the initial allocation for the number
 *            of integers to be stored in key data. 
 *            
 *            In general the initialization defaults should be
 *            fine. All three are grown automatically as needed, as
 *            you add keys to the hash.
 *            
 *            "key data" means <n> <start>/<end> pairs, plus <n>
 *            itself: it takes 2n+1 integers to store a <P7_COORD2>
 *            array of length <n>.
 *            
 *            <hashsize> must be a power of 2; remember that if you
 *            pass a non-default value.
 *            
 * Args:      init_hashsize : initial hashtable size. Power of 2; >0.
 *            init_keyalloc : initial allocation for # keys. >0.
 *            init_calloc   : initial allocation for key data. >0.
 *
 * Returns:   pointer to the new <P7_COORDS2_HASH> object on success.
 *
 * Throws:    <NULL> on allocation failure.
 */
P7_COORDS2_HASH *
p7_coords2_hash_Create(int32_t init_hashsize, int32_t init_nkeyalloc, int32_t init_calloc)
{
  P7_COORDS2_HASH *ch = NULL;
  int32_t          i;
  int              status;

  ESL_DASSERT1(( init_hashsize == 0 || (init_hashsize && ((init_hashsize & (init_hashsize-1)) == 0)))); /* hashsize is a power of 2 (bitshifting trickery) */
  
  ESL_ALLOC(ch, sizeof(P7_COORDS2_HASH));
  ch->hashtable  = NULL;
  ch->key_offset = NULL;
  ch->nxt        = NULL;
  ch->cmem       = NULL;

  ch->nkeys      = 0;
  ch->cn         = 0;

  ch->hashsize   = (init_hashsize  > 0 ? init_hashsize  : 128);
  ch->kalloc     = (init_nkeyalloc > 0 ? init_nkeyalloc : 128);
  ch->calloc     = (init_calloc    > 0 ? init_calloc    : 2048);
  
  ESL_ALLOC(ch->hashtable, sizeof(int32_t) * ch->hashsize);
  for (i = 0; i < ch->hashsize; i++) ch->hashtable[i] = -1;

  ESL_ALLOC(ch->key_offset, sizeof(int32_t) * ch->kalloc);
  ESL_ALLOC(ch->nxt,        sizeof(int32_t) * ch->kalloc);
  ESL_ALLOC(ch->cmem,       sizeof(int32_t) * ch->calloc);
  return ch;
  
 ERROR:
  p7_coords2_hash_Destroy(ch);
  return NULL;
}

P7_GBANDS *
p7_gbands_Create(void)
{
  P7_GBANDS *bnd           = NULL;
  int        init_segalloc = 4;
  int        init_rowalloc = 64;
  int        status;

  ESL_ALLOC(bnd, sizeof(P7_GBANDS));
  bnd->nseg  = 0;
  bnd->nrow  = 0;
  bnd->L     = 0;
  bnd->M     = 0;
  bnd->ncell = 0;
  bnd->imem  = NULL;
  bnd->kmem  = NULL;


  ESL_ALLOC(bnd->imem, sizeof(int) * init_segalloc * 2); /* *2: for ia, ib pairs */
  ESL_ALLOC(bnd->kmem, sizeof(int) * init_rowalloc * p7_GBANDS_NK);
  bnd->segalloc = init_segalloc;
  bnd->rowalloc = init_rowalloc;

  return bnd;
  
 ERROR:
  p7_gbands_Destroy(bnd);
  return NULL;
}

ESL_HMX *
esl_hmx_Create(int allocL, int allocM)
{
  ESL_HMX *mx = NULL;
  int      i;
  int      status;
  
  ESL_ALLOC(mx, sizeof(ESL_HMX));
  mx->dp_mem = NULL;
  mx->dp     = NULL;
  mx->sc     = NULL;

  ESL_ALLOC(mx->dp_mem, sizeof(float) * (allocL+1) * allocM);
  mx->ncells = (allocL+1) * allocM;
  
  ESL_ALLOC(mx->dp, sizeof (float *) * (allocL+1));
  ESL_ALLOC(mx->sc, sizeof (float)   * (allocL+2));
  mx->allocR = allocL+1;

  for (i = 0; i <= allocL; i++)
    mx->dp[i] = mx->dp_mem + i*allocM;
  mx->validR = allocL+1;
  mx->allocM = allocM;

  mx->L = 0;
  mx->M = 0;
  return mx;

 ERROR:
  esl_hmx_Destroy(mx);
  return NULL;
}

/* Function:  esl_recorder_Create()
 * Synopsis:  Create an <ESL_RECORDER>.
 * Incept:    SRE, Fri Dec 25 16:27:40 2009 [Casa de Gatos]
 *
 * Purpose:   Allocate a new <ESL_RECORDER> that will read
 *            line-by-line from input stream <fp>, saving
 *            a history of up to <maxlines> lines.
 *
 * Returns:   pointer to the new <ESL_RECORDER> on success.
 *
 * Throws:    <NULL> on allocation failure.
 */
ESL_RECORDER *
esl_recorder_Create(FILE *fp, int maxlines)
{
  ESL_RECORDER *rc = NULL;
  int           i;
  int           status;

  ESL_ALLOC(rc, sizeof(ESL_RECORDER));
  rc->fp         = fp;
  rc->line       = NULL;
  rc->nalloc     = maxlines;
  rc->lalloc     = NULL;
  rc->offset     = NULL;
  rc->nread      = 0;
  rc->ncurr      = 0;
  rc->baseline   = 0;
  rc->markline   = -1;

  ESL_ALLOC(rc->line,   sizeof(char *) * rc->nalloc);
  for (i = 0; i < rc->nalloc; i++) rc->line[i]   = NULL;

  ESL_ALLOC(rc->lalloc, sizeof(int)    * rc->nalloc);
  for (i = 0; i < rc->nalloc; i++) rc->lalloc[i] = 0;

  ESL_ALLOC(rc->offset, sizeof(off_t)  * rc->nalloc);
  for (i = 0; i < rc->nalloc; i++) rc->offset[i] = 0;

  return rc;

 ERROR:
  esl_recorder_Destroy(rc);
  return NULL;
}

/* Function:  p7_bg_Create()
 * Incept:    SRE, Fri Jan 12 13:32:51 2007 [Janelia]
 *
 * Purpose:   Allocate a <P7_BG> object for digital alphabet <abc>,
 *            initializes it to appropriate default values, and
 *            returns a pointer to it.
 *            
 *            For protein models, default iid background frequencies
 *            are set (by <p7_AminoFrequencies()>) to average
 *            SwissProt residue composition. For DNA, RNA and other
 *            alphabets, default frequencies are set to a uniform
 *            distribution.
 *            
 *            The model composition <bg->mcomp[]> is not initialized
 *            here; neither is the filter null model <bg->fhmm>.  To
 *            use the filter null model, caller will want to
 *            initialize these fields by calling
 *            <p7_bg_SetFilterByHMM()>.
 *
 * Throws:    <NULL> on allocation failure.
 *
 * Xref:      STL11/125.
 */
P7_BG *
p7_bg_Create(const ESL_ALPHABET *abc)
{
  P7_BG *bg = NULL;
  int    status;

  ESL_ALLOC(bg, sizeof(P7_BG));
  bg->f     = NULL;
  bg->fhmm  = NULL;

  ESL_ALLOC(bg->f,     sizeof(float) * abc->K);
  if ((bg->fhmm = esl_hmm_Create(abc, 2)) == NULL) goto ERROR;

  if       (abc->type == eslAMINO)
    {
      if (p7_AminoFrequencies(bg->f) != eslOK) goto ERROR;
    }
  else
    esl_vec_FSet(bg->f, abc->K, 1. / (float) abc->K);

  bg->p1    = 350./351.;
  bg->omega = 1./256.;
  bg->abc   = abc;
  return bg;

 ERROR:
  p7_bg_Destroy(bg);
  return NULL;
}

/* Function:  p7_oprofile_CreateBlock()
 * Synopsis:  Create a new block of empty <P7_OM_BLOCK>.
 * Incept:    
 *
 * Purpose:   Creates a block of empty <P7_OM_BLOCK> profile objects.
 *            
 * Returns:   a pointer to the new <P7_OM_BLOCK>. Caller frees this
 *            with <p7_oprofile_DestroyBlock()>.
 *
 * Throws:    <NULL> if allocation fails.
 */
P7_OM_BLOCK *
p7_oprofile_CreateBlock(int count)
{
  int i = 0;

  P7_OM_BLOCK *block = NULL;
  int status = eslOK;

  ESL_ALLOC(block, sizeof(*block));

  block->count = 0;
  block->listSize = 0;
  block->list  = NULL;

  ESL_ALLOC(block->list, sizeof(P7_OPROFILE *) * count);
  block->listSize = count;

  for (i = 0; i < count; ++i)
    {
      block->list[i] = NULL;
    }

  return block;

 ERROR:
  if (block != NULL)
    {
      if (block->list != NULL)  free(block->list);
      free(block);
    }
  
  return NULL;
}

static int
parse_replace_string(const char *rstring, char **ret_from, char **ret_to)
{
  int    status;
  int    rlen, mid, i;
  int    is_valid = FALSE;
  char  *from = NULL;
  char  *to   = NULL;

  /* Note: we could use ESL_REGEXP but then multiple ':'s in rstring could cause problems */
  rlen = strlen(rstring);
  /* check validity of rstring: must be "<s1>:<s2>" with len(<s1>)==len(<s2>) */
  if((rlen % 2) != 0) { /* odd num chars, good */
    mid = rlen / 2;
    if(rstring[mid] == ':') { /* middle character is ':', good */
      ESL_ALLOC(from, sizeof(char) * (mid+1));
      ESL_ALLOC(to,   sizeof(char) * (mid+1));
      for(i = 0;     i < mid;  i++) from[i]       = rstring[i];
      for(i = mid+1; i < rlen; i++) to[i-(mid+1)] = rstring[i];
      from[mid] = '\0';
      to[mid]   = '\0';
      is_valid = TRUE;
    }
  }
  if(! is_valid) esl_fatal("--replace takes arg of <s1>:<s2> with len(<s1>) == len(<s2>); %s not recognized", rstring);
  *ret_from = from;
  *ret_to   = to;

  return eslOK;

 ERROR: 
  if(from != NULL) free(from);
  if(to   != NULL) free(to);
  return status;
}

/* Function:  esl_msaweight_IDFilter()
 * Synopsis:  Filter by %ID.
 * Incept:    ER, Wed Oct 29 10:06:43 2008 [Janelia]
 * 
 * Purpose:   Constructs a new alignment by removing near-identical 
 *            sequences from a given alignment (where identity is 
 *            calculated *based on the alignment*).
 *            Does not affect the given alignment.
 *            Keeps earlier sequence, discards later one. 
 *           
 *            Usually called as an ad hoc sequence "weighting" mechanism.
 *           
 * Limitations:
 *            Unparsed Stockholm markup is not propagated into the
 *            new alignment.
 *           
 * Return:    <eslOK> on success, and the <newmsa>.
 *
 * Throws:    <eslEMEM> on allocation error. <eslEINVAL> if a pairwise
 *            identity calculation fails because of corrupted sequence 
 *            data. In either case, the <msa> is unmodified.
 *
 * Xref:      squid::weight.c::FilterAlignment().
 */
int
esl_msaweight_IDFilter(const ESL_MSA *msa, double maxid, ESL_MSA **ret_newmsa)
{
  int     *list   = NULL;               /* array of seqs in new msa */
  int     *useme  = NULL;               /* TRUE if seq is kept in new msa */
  int      nnew;			/* number of seqs in new alignment */
  double   ident;                       /* pairwise percentage id */
  int      i,j;                         /* seqs counters*/
  int      remove;                      /* TRUE if sq is to be removed */
  int      status;
  
  /* Contract checks
   */
  ESL_DASSERT1( (msa       != NULL) );
  ESL_DASSERT1( (msa->nseq >= 1)    );
  ESL_DASSERT1( (msa->alen >= 1)    );

  /* allocate */
  ESL_ALLOC(list,  sizeof(int) * msa->nseq);
  ESL_ALLOC(useme, sizeof(int) * msa->nseq);
  esl_vec_ISet(useme, msa->nseq, 0); /* initialize array */

  /* find which seqs to keep (list) */
  nnew = 0;
  for (i = 0; i < msa->nseq; i++)
    {
      remove = FALSE;
      for (j = 0; j < nnew; j++)
	{
	  if (! (msa->flags & eslMSA_DIGITAL)) {
	    if ((status = esl_dst_CPairId(msa->aseq[i], msa->aseq[list[j]], &ident, NULL, NULL))       != eslOK) goto ERROR;
	  } 
#ifdef eslAUGMENT_ALPHABET
	  else {
	    if ((status = esl_dst_XPairId(msa->abc, msa->ax[i], msa->ax[list[j]], &ident, NULL, NULL)) != eslOK) goto ERROR;
	  }
#endif
	  
	  if (ident > maxid)
	    { 
	      remove = TRUE; 
	      break; 
	    }
	}
      if (remove == FALSE) {
	list[nnew++] = i;
	useme[i]     = TRUE;
      }
    }
  if ((status = esl_msa_SequenceSubset(msa, useme, ret_newmsa)) != eslOK) goto ERROR;
 
  free(list);
  free(useme);
  return eslOK;

 ERROR:
  if (list  != NULL) free(list);
  if (useme != NULL) free(useme);
  return status;
}

/* Function:  esl_hyperexp_Create()
 *
 * Purpose:   Creates an object to hold parameters for a <K>-component
 *            hyperexponential. 
 *
 *            Parameters in the object are initialized 
 *            ($q_k = \frac{1}{K}$, $\lambda_k = 1$, $\mu = 0$), but
 *            the caller will want to set these according to its own
 *            purposes.
 *
 * Args:      K  - number of components in the mixture
 *
 * Returns:   ptr to newly allocated/initialized <ESL_HYPEREXP> object.
 *
 * Throws:    NULL on allocation failure.
 */
ESL_HYPEREXP *
esl_hyperexp_Create(int K)
{
  int           status;
  ESL_HYPEREXP *h = NULL;
  int           k;

  ESL_ALLOC(h, sizeof(ESL_HYPEREXP));
  h->q = h->lambda = h->wrk = NULL;
  h->fixlambda = NULL;
  h->K         = K;
  h->fixmix    = FALSE;

  ESL_ALLOC(h->q,         sizeof(double) * K);
  ESL_ALLOC(h->lambda,    sizeof(double) * K);
  ESL_ALLOC(h->wrk,       sizeof(double) * K);
  ESL_ALLOC(h->fixlambda, sizeof(char)   * K);

  for (k = 0; k < K; k++)
    {
      h->q[k]        = 1. / (double) K;
      h->lambda[k]   = 1.;
      h->fixlambda[k]= 0;
    }
  h->mu = 0.;
  return h;
  
 ERROR:
  esl_hyperexp_Destroy(h);
  return NULL;
}

/* Function:  p7_tophits_Create()
 * Synopsis:  Allocate a hit list.
 *
 * Purpose:   Allocates a new <P7_TOPHITS> hit list, for an initial
 *            allocation of <int_hit_alloc> hits (this will be grown
 *            later as needed). Return a pointer to it.
 *            
 * Args:      init_hit_alloc  - initial allocation size, # of hits.
 *                              Often p7_TOPHITS_DEFAULT_INIT_ALLOC.
 *
 * Throws:    <NULL> on allocation failure.
 */
P7_TOPHITS *
p7_tophits_Create(int init_hit_alloc)
{
  P7_TOPHITS *h = NULL;
  int         status;

  ESL_ALLOC(h, sizeof(P7_TOPHITS));
  h->hit    = NULL;
  h->unsrt  = NULL;

  if (( h->unsrt = p7_hit_Create(init_hit_alloc) ) == NULL) goto ERROR;
  ESL_ALLOC(h->hit, sizeof(P7_HIT *) * init_hit_alloc);

  h->Nalloc               = init_hit_alloc;
  h->N                    = 0;
  h->nreported            = 0;
  h->nincluded            = 0;
  h->is_sorted_by_sortkey = TRUE;     /* but only because there's 0 hits */
  h->is_sorted_by_seqidx  = FALSE;
  h->hit[0]               = h->unsrt; /* if you're going to call it "sorted" when it contains just one hit, you need this */
  return h;

 ERROR:
  p7_tophits_Destroy(h);
  return NULL;
}

/* Function: annotate_model()
 * 
 * Purpose:  Transfer rf, cs, and other optional annotation from the alignment
 *           to the new model.
 * 
 * Args:     hmm       - [M] new model to annotate
 *           matassign - which alignment columns are MAT; [1..alen]
 *           msa       - alignment, including annotation to transfer
 *           
 * Return:   <eslOK> on success.
 *
 * Throws:   <eslEMEM> on allocation error.
 */
static int
annotate_model(P7_HMM *hmm, int *matassign, ESL_MSA *msa)
{                      
  int   apos;			/* position in matassign, 1.alen  */
  int   k;			/* position in model, 1.M         */
  int   status;

  /* Reference coord annotation  */
  if (msa->rf != NULL) {
    ESL_ALLOC(hmm->rf, sizeof(char) * (hmm->M+2));
    hmm->rf[0] = ' ';
    for (apos = k = 1; apos <= msa->alen; apos++) 
      if (matassign[apos]) hmm->rf[k++] = msa->rf[apos-1]; /* watch off-by-one in msa's rf */
    hmm->rf[k] = '\0';
    hmm->flags |= p7H_RF;
  }

  /* Model mask annotation  */
  if (msa->mm != NULL) {
    ESL_ALLOC(hmm->mm, sizeof(char) * (hmm->M+2));
    hmm->mm[0] = ' ';
    for (apos = k = 1; apos <= msa->alen; apos++)
      if (matassign[apos]) hmm->mm[k++] = msa->mm[apos-1];
    hmm->mm[k] = '\0';
    hmm->flags |= p7H_MMASK;
  }

  /* Consensus structure annotation */
  if (msa->ss_cons != NULL) {
    ESL_ALLOC(hmm->cs, sizeof(char) * (hmm->M+2));
    hmm->cs[0] = ' ';
    for (apos = k = 1; apos <= msa->alen; apos++)
      if (matassign[apos]) hmm->cs[k++] = msa->ss_cons[apos-1];
    hmm->cs[k] = '\0';
    hmm->flags |= p7H_CS;
  }

  /* Surface accessibility annotation */
  if (msa->sa_cons != NULL) {
    ESL_ALLOC(hmm->ca, sizeof(char) * (hmm->M+2));
    hmm->ca[0] = ' ';
    for (apos = k = 1; apos <= msa->alen; apos++)
      if (matassign[apos]) hmm->ca[k++] = msa->sa_cons[apos-1];
    hmm->ca[k] = '\0';
    hmm->flags |= p7H_CA;
  }

  /* The alignment map (1..M in model, 1..alen in alignment) */
  ESL_ALLOC(hmm->map, sizeof(int) * (hmm->M+1));
  hmm->map[0] = 0;
  for (apos = k = 1; apos <= msa->alen; apos++)
    if (matassign[apos]) hmm->map[k++] = apos;
  hmm->flags |= p7H_MAP;

  return eslOK;

 ERROR:
  return status;
}

/* Function:  cp9_Copy()
 * Synopsis:  Copy a CM plan 9 HMM.
 *
 * Purpose:   Copies cp9 hmm <src> to cp9 hmm <dst>, where <dst>
 *            has already been allocated to be of sufficient size.
 *
 *            <src> should be properly normalized, no check is done to
 *            ensure that. If <src> is logoddsified (src->flags &
 *            CPLAN9_HASBITS) its bit scores will be copied to <dst>,
 *            otherwise they are invalid and won't be copied.
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEMEM> on allocation error; <eslEINVAL> if <dst> is too small
 *            to fit <src>.
 */
int
cp9_Copy(const CP9_t *src, CP9_t *dst)
{
    int status;
    int k;
    int src_has_bits = (src->flags & CPLAN9_HASBITS) ? TRUE : FALSE;

    if (src->M != dst->M) return eslEINVAL;

    dst->abc = src->abc;

    for(k = 0; k <= src->M; k++) {
        esl_vec_FCopy(src->t[k],   cp9_NTRANS,  dst->t[k]);
        esl_vec_FCopy(src->mat[k], src->abc->K, dst->mat[k]);
        esl_vec_FCopy(src->ins[k], src->abc->K, dst->ins[k]);
    }
    esl_vec_FCopy(src->begin, src->M+1, dst->begin);
    esl_vec_FCopy(src->end,   src->M+1, dst->end);
    if(src_has_bits) {
        esl_vec_ICopy(src->bsc_mem, src->M+1, dst->bsc_mem);
        esl_vec_ICopy(src->esc_mem, src->M+1, dst->esc_mem);
    }

    /* exploit linear-memory of these 2d arrays */
    if(src_has_bits) {
        esl_vec_ICopy(src->tsc_mem, cp9_NTRANS   * (src->M+1), dst->tsc_mem);
        esl_vec_ICopy(src->msc_mem, src->abc->Kp * (src->M+1), dst->msc_mem);
        esl_vec_ICopy(src->isc_mem, src->abc->Kp * (src->M+1), dst->isc_mem);
        esl_vec_ICopy(src->otsc,    cp9O_NTRANS  * (src->M+1), dst->otsc);
    }

    /* EL info */
    dst->el_self     = src->el_self;
    dst->el_selfsc   = src->el_selfsc;
    esl_vec_ICopy(src->has_el,     src->M+1,    dst->has_el);
    esl_vec_ICopy(src->el_from_ct, src->M+2,    dst->el_from_ct);
    for(k = 0; k <= src->M+1; k++) {
        if(src->el_from_ct[k] > 0) {
            ESL_ALLOC(dst->el_from_idx[k],  sizeof(int) * src->el_from_ct[k]);
            ESL_ALLOC(dst->el_from_cmnd[k], sizeof(int) * src->el_from_ct[k]);
            esl_vec_ICopy(src->el_from_idx[k],  src->el_from_ct[k], dst->el_from_idx[k]);
            esl_vec_ICopy(src->el_from_cmnd[k], src->el_from_ct[k], dst->el_from_cmnd[k]);
        }
    }

    dst->null2_omega = src->null2_omega;
    dst->null3_omega = src->null3_omega;
    esl_vec_FCopy(src->null, src->abc->K, dst->null);

    dst->p1    = src->p1;
    dst->flags = src->flags;

    return eslOK;

ERROR:
    return status;
}

/* Function:  esl_msacluster_SingleLinkage()
 * Synopsis:  Single linkage clustering by percent identity.
 * Incept:    SRE, Sun Nov  5 10:11:45 2006 [Janelia]
 *
 * Purpose:   Perform single link clustering of the sequences in 
 *            multiple alignment <msa>. Any pair of sequences with
 *            percent identity $\geq$ <maxid> are linked (using
 *            the definition from the \eslmod{distance} module).
 *            
 *            The resulting clustering is optionally returned in one
 *            or more of <opt_c>, <opt_nin>, and <opt_nc>.  The
 *            <opt_c[0..nseq-1]> array assigns a cluster index
 *            <(0..nc-1)> to each sequence. For example, <c[4] = 1>
 *            means that sequence 4 is assigned to cluster 1.  The
 *            <opt_nin[0..nc-1]> array is the number of sequences
 *            in each cluster. <opt_nc> is the number of clusters.
 *
 *            Importantly, this algorithm runs in $O(N)$ memory, and
 *            produces one discrete clustering. Compare to
 *            <esl_tree_SingleLinkage()>, which requires an $O(N^2)$ 
 *            adjacency matrix, and produces a hierarchical clustering
 *            tree.
 *            
 *            The algorithm is worst case $O(LN^2)$ time, for $N$
 *            sequences of length $L$. However, the worst case is no
 *            links at all, and this is unusual. More typically, time
 *            scales as about $LN \log N$. The best case scales as
 *            $LN$, when there is just one cluster in a completely
 *            connected graph.
 *            
 * Args:      msa     - multiple alignment to cluster
 *            maxid   - pairwise identity threshold: cluster if $\geq$ <maxid>
 *            opt_c   - optRETURN: cluster assignments for each sequence, [0..nseq-1]
 *            opt_nin - optRETURN: number of seqs in each cluster, [0..nc-1] 
 *            opt_nc  - optRETURN: number of clusters        
 *
 * Returns:   <eslOK> on success; the <opt_c[0..nseq-1]> array contains
 *            cluster indices <0..nc-1> assigned to each sequence; the
 *            <opt_nin[0..nc-1]> array contains the number of seqs in
 *            each cluster; and <opt_nc> contains the number of
 *            clusters. The <opt_c> array and <opt_nin> arrays will be
 *            allocated here, if non-<NULL>, and must be free'd by the
 *            caller. The input <msa> is unmodified.
 *            
 *            The caller may pass <NULL> for either <opt_c> or
 *            <opt_nc> if it is only interested in one of the two
 *            results.
 *
 * Throws:    <eslEMEM> on allocation failure, and <eslEINVAL> if a pairwise
 *            comparison is invalid (which means the MSA is corrupted, so it
 *            shouldn't happen). In either case, <opt_c> and <opt_nin> are set to <NULL>
 *            and <opt_nc> is set to 0, and the <msa> is unmodified.
 */
int
esl_msacluster_SingleLinkage(const ESL_MSA *msa, double maxid, 
			     int **opt_c, int **opt_nin, int *opt_nc)

{
  int   status;
  int  *workspace  = NULL;
  int  *assignment = NULL;
  int  *nin        = NULL;
  int   nc;
  int   i;
#ifdef eslAUGMENT_ALPHABET
  struct msa_param_s param;
#endif

  /* Allocations */
  ESL_ALLOC(workspace,  sizeof(int) * msa->nseq * 2);
  ESL_ALLOC(assignment, sizeof(int) * msa->nseq);

  /* call to SLC API: */
  if (! (msa->flags & eslMSA_DIGITAL))
    status = esl_cluster_SingleLinkage((void *) msa->aseq, (size_t) msa->nseq, sizeof(char *),
				       msacluster_clinkage, (void *) &maxid, 
				       workspace, assignment, &nc);
#ifdef eslAUGMENT_ALPHABET
  else {
    param.maxid = maxid;
    param.abc   = msa->abc;
    status = esl_cluster_SingleLinkage((void *) msa->ax, (size_t) msa->nseq, sizeof(ESL_DSQ *),
				       msacluster_xlinkage, (void *) &param, 
				       workspace, assignment, &nc);
  }
#endif

  if (opt_nin != NULL) 
    {
      ESL_ALLOC(nin, sizeof(int) * nc);
      for (i = 0; i < nc; i++) nin[i] = 0;
      for (i = 0; i < msa->nseq; i++)
	nin[assignment[i]]++;
      *opt_nin = nin;
    }

  /* cleanup and return */
  free(workspace);
  if (opt_c  != NULL) *opt_c  = assignment; else free(assignment);
  if (opt_nc != NULL) *opt_nc = nc;
//.........这里部分代码省略.........

/* Function:  esl_hxp_FitCompleteBinned()
 *
 * Purpose:   Given a histogram <g> with binned observations, where each
 *            bin i holds some number of observed samples x with values from 
 *            lower bound l to upper bound u (that is, $l < x \leq u$),
 *            and given a starting guess <h> for hyperexponential parameters;
 *
 *            Find maximum likelihood parameters <h> by conjugate gradient
 *            descent, starting from the initial <h> and leaving the
 *            optimized solution in <h>.
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEMEM> on allocation error, and <h> is left in its
 *            initial state.
 */
int 
esl_hxp_FitCompleteBinned(ESL_HISTOGRAM *g, ESL_HYPEREXP *h)
{
  struct hyperexp_binned_data data;
  int     status;
  double *p   = NULL;
  double *u   = NULL;
  double *wrk = NULL;
  double  fx;
  int     i;
  double  tol = 1e-6;
  int     np;

  np = 0;
  if (! h->fixmix) np = h->K-1;  /* K-1 mix coefficients...      */
  for (i = 0; i < h->K; i++)     /* ...and up to K lambdas free. */
    if (! h->fixlambda[i]) np++;

  ESL_ALLOC(p,   sizeof(double) * np);
  ESL_ALLOC(u,   sizeof(double) * np);
  ESL_ALLOC(wrk, sizeof(double) * np * 4);

  /* Copy shared info into the "data" structure  */
  data.g     = g;
  data.h     = h;

  /* From h, create the parameter vector. */
  hyperexp_pack_paramvector(p, np, h);

  /* Define the step size vector u.
   */
  for (i = 0; i < np; i++) u[i] = 1.0;

  /* Feed it all to the mighty optimizer.
   */
  status = esl_min_ConjugateGradientDescent(p, u, np, 
					    &hyperexp_complete_binned_func, 
					    &hyperexp_complete_binned_gradient,
					    (void *) (&data), tol, wrk, &fx);
  if (status != eslOK) goto ERROR;

  /* Convert the final parameter vector back to a hyperexponential
   */
  hyperexp_unpack_paramvector(p, np, h);
  
  free(p);
  free(u);
  free(wrk);
  esl_hyperexp_SortComponents(h);
  return eslOK;

 ERROR:
  if (p   != NULL) free(p);
  if (u   != NULL) free(u);
  if (wrk != NULL) free(wrk);
  return status;
}

/* mpi_worker()
 * The main control for an MPI worker process.
 */
static void
mpi_worker(ESL_GETOPTS *go, struct cfg_s *cfg)
{
  int             xstatus = eslOK;
  int             status;
  P7_HMM         *hmm     = NULL;
  char           *wbuf    = NULL;
  double         *xv      = NULL; /* result: array of N scores */
  int            *av      = NULL; /* optional result: array of N alignment lengths */
  int             wn      = 0;
  char            errbuf[eslERRBUFSIZE];
  int             pos;
 
  /* Worker initializes */
  if ((status = minimum_mpi_working_buffer(go, cfg->N, &wn)) != eslOK) xstatus = status;
  ESL_ALLOC(wbuf, wn * sizeof(char));
  ESL_ALLOC(xv,   cfg->N * sizeof(double) + 2);	
  if (esl_opt_GetBoolean(go, "-a"))
    ESL_ALLOC(av, cfg->N * sizeof(int));

  /* Main worker loop */
  while (p7_hmm_mpi_Recv(0, 0, MPI_COMM_WORLD, &wbuf, &wn, &(cfg->abc), &hmm) == eslOK) 
    {
      if (esl_opt_GetBoolean(go, "--recal")) {
	if (( status = recalibrate_model(go, cfg, errbuf, hmm))     != eslOK) goto CLEANERROR;
      }
      if ((status = process_workunit(go, cfg, errbuf, hmm, xv, av)) != eslOK) goto CLEANERROR;

      pos = 0;
      MPI_Pack(&status, 1,      MPI_INT,    wbuf, wn, &pos, MPI_COMM_WORLD);
      MPI_Pack(xv,      cfg->N, MPI_DOUBLE, wbuf, wn, &pos, MPI_COMM_WORLD);
      if (esl_opt_GetBoolean(go, "-a"))
	MPI_Pack(av,    cfg->N, MPI_INT,    wbuf, wn, &pos, MPI_COMM_WORLD);
      MPI_Send(wbuf, pos, MPI_PACKED, 0, 0, MPI_COMM_WORLD);
      p7_hmm_Destroy(hmm);
    }

  free(wbuf);
  free(xv);
  if (av != NULL) free(av);
  return;

 CLEANERROR:
  pos = 0;
  MPI_Pack(&status, 1,                MPI_INT,  wbuf, wn, &pos, MPI_COMM_WORLD);
  MPI_Pack(errbuf,  eslERRBUFSIZE,    MPI_CHAR, wbuf, wn, &pos, MPI_COMM_WORLD);
  MPI_Send(wbuf, pos, MPI_PACKED, 0, 0, MPI_COMM_WORLD);
  if (wbuf != NULL) free(wbuf);
  if (hmm  != NULL) p7_hmm_Destroy(hmm);
  if (xv   != NULL) free(xv);
  if (av   != NULL) free(av);
  return;

 ERROR:
  p7_Fail("Allocation error in mpi_worker");
}

/* Function: CreateCP9Matrix()
 * based on  CreatePlan7Matrix() <-- this function's comments below  
 * Purpose:  Create a dynamic programming matrix for standard Forward,
 *           Backward, or Viterbi, with scores kept as scaled log-odds
 *           integers. Keeps 2D arrays compact in RAM in an attempt 
 *           to maximize cache hits. 
 *           
 *           The mx structure can be dynamically grown, if a new
 *           HMM or seq exceeds the currently allocated size. Dynamic
 *           growing is more efficient than an alloc/free of a whole
 *           matrix for every new target. The ResizePlan7Matrix()
 *           call does this reallocation, if needed. Here, in the
 *           creation step, we set up some pads - to inform the resizing
 *           call how much to overallocate when it realloc's. 
 *           
 * Args:     N     - N+1 rows are allocated, usually N == 1 for 
 *                   scanning in memory efficient mode, or N == L, length of sequence.  
 *           M     - size of model in nodes
 *                 
 * Return:   mx
 *           mx is allocated here. Caller frees with FreeCP9Matrix(mx).
 */
CP9_MX *
CreateCP9Matrix(int N, int M)
{
  int status;
  CP9_MX *mx;
  int i;

  ESL_ALLOC(mx,      sizeof(CP9_MX));
  ESL_ALLOC(mx->mmx, sizeof(int *) * (N+1));
  ESL_ALLOC(mx->imx, sizeof(int *) * (N+1));
  ESL_ALLOC(mx->dmx, sizeof(int *) * (N+1));
  ESL_ALLOC(mx->elmx,sizeof(int *) * (N+1)); 
  /* slightly wasteful, some nodes can't go to EL (for ex: right half of MATPs) */
  ESL_ALLOC(mx->erow,    sizeof(int) * (N+1));
  ESL_ALLOC(mx->mmx_mem, sizeof(int) * ((N+1)*(M+1)));
  ESL_ALLOC(mx->imx_mem, sizeof(int) * ((N+1)*(M+1)));
  ESL_ALLOC(mx->dmx_mem, sizeof(int) * ((N+1)*(M+1)));
  ESL_ALLOC(mx->elmx_mem,sizeof(int) * ((N+1)*(M+1)));

  /* The indirect assignment below looks wasteful; it's actually
   * used for aligning data on 16-byte boundaries as a cache 
   * optimization in the fast altivec implementation
   */
  mx->mmx[0] = (int *) mx->mmx_mem;
  mx->imx[0] = (int *) mx->imx_mem;
  mx->dmx[0] = (int *) mx->dmx_mem;
  mx->elmx[0]= (int *) mx->elmx_mem;
  for (i = 1; i <= N; i++)
    {
      mx->mmx[i] = mx->mmx[0] + (i*(M+1));
      mx->imx[i] = mx->imx[0] + (i*(M+1));
      mx->dmx[i] = mx->dmx[0] + (i*(M+1));
      mx->elmx[i]= mx->elmx[0]+ (i*(M+1));
    }

  mx->M = M;
  mx->rows = N;
  mx->kmin = NULL;
  mx->kmax = NULL;
  mx->ncells_allocated = (M+1) * (N+1);
  mx->ncells_valid     = (M+1) * (N+1);
  mx->size_Mb =  (float) sizeof(CP9_MX);
  mx->size_Mb += (float) (sizeof(int *) * (mx->rows+1) * 4); /* mx->*mx ptrs */
  mx->size_Mb += (float) (sizeof(int)   * (mx->rows+1) * (M+1) * 4); /* mx->*mx_mem */
  mx->size_Mb += (float) (sizeof(int)   * (mx->rows+1));             /* mx->erow */
  mx->size_Mb /= 1000000.;

  return mx;

 ERROR:
  cm_Fail("Memory allocation error.");
  return NULL; /* never reached */
}

/* map_sub_msas
 *                   
 * msa1 and msa2 contain the same named sequences, msa1 contains a superset 
 * of the columns in msa2. Determine which of the msa1 columns the msa2
 * columns correspond to.
 */
static int
map_sub_msas(const ESL_GETOPTS *go, char *errbuf, ESL_MSA *msa1, ESL_MSA *msa2, char **ret_msa1_to_msa2_mask)
{
  int status;
  int  apos1, apos2;          /* counters over alignment position in msa1, msa2 respectively */
  int i;
  int *msa1_to_msa2_map;    /* [0..apos1..msa1->alen] msa2 alignment position that apos1 corresponds to */
  char *mask;

  /* contract check */
  if(! (msa1->flags & eslMSA_DIGITAL)) ESL_FAIL(eslEINVAL, errbuf, "in map_sub_msas() msa1 (%s) not digitized.\n", esl_opt_GetArg(go, 1));
  if(! (msa2->flags & eslMSA_DIGITAL)) ESL_FAIL(eslEINVAL, errbuf, "in map_sub_msas() msa2 (%s) not digitized.\n", esl_opt_GetString(go, "--submap"));
  if(msa1->alen <= msa2->alen) ESL_FAIL(eslEINVAL, errbuf, "in map_sub_msas() alignment length for msa1 (%" PRId64 "d) <= length for msa2 (%" PRId64 ")\n", msa1->alen, msa2->alen);
  
  ESL_ALLOC(mask, sizeof(char) * (msa1->alen+1));
  for(apos1 = 0; apos1 < msa1->alen; apos1++) mask[apos1] = '0';
  mask[msa1->alen] = '\0';

  ESL_ALLOC(msa1_to_msa2_map, sizeof(int) * (msa1->alen+1));
  esl_vec_ISet(msa1_to_msa2_map, (msa1->alen+1), -1);

  /* both alignments must have same 'named' sequences in same order */
  if(msa1->nseq != msa2->nseq) ESL_FAIL(eslEINVAL, errbuf, "in map_sub_msas() msa1 has %d sequences, msa2 has %d sequences\n", msa1->nseq, msa2->nseq);
  for(i = 0; i < msa1->nseq; i++) { 
    if(strcmp(msa1->sqname[i], msa2->sqname[i]) != 0) ESL_FAIL(eslEINVAL, errbuf, "in map_sub_msas() msa1 seq %d is named %s, msa2 seq %d is named %s\n", i, msa1->sqname[i], i, msa2->sqname[i]);
  }

  apos1 = 1;
  apos2 = 1;
  while((apos2 <= msa2->alen) || (apos1 <= msa1->alen)) { /* determine which apos1 (alignment column in msa1), apos2 (alignment column in msa2) corresponds to */
    for(i = 0; i < msa1->nseq; i++) { 
      if(msa1->ax[i][apos1] != msa2->ax[i][apos2]) { 
	apos1++; 
	break; /* try next apos1 */ 
      }
    }	
    if(i == msa1->nseq) { /* found a match */
      msa1_to_msa2_map[apos1] = apos2;
      mask[(apos1-1)] = '1';
      apos1++;
      apos2++;
    }
  }
  if((apos1 != (msa1->alen+1)) || (apos2 != (msa2->alen+1))) ESL_FAIL(eslEINVAL, errbuf, "in map_sub_msas(), failure mapping alignments, end of loop apos1-1 = %d (msa1->alen: %" PRId64 ") and apos2-1 = %d (msa2->alen: %" PRId64 ")\n", apos1-1, msa1->alen, apos2-1, msa2->alen);

  free(msa1_to_msa2_map);
  *ret_msa1_to_msa2_mask = mask;
  return eslOK;
  
 ERROR: 
  return status;
}