/* copyright (c) 1998, 1999 William R. Pearson and the U. of Virginia */

/*  $Name: fa_34_26_5 $ - $Id: dropfs2.c,v 1.40 2007/02/26 21:56:59 wrp Exp $ */

/* changed to return 2.0, rather than -1.0, for failure */

/* Feb 4, 2005 - modifications to allow searches with ktup=2 for very
   long queries.  This is a temporary solution to savemax(), spam()
   which do not preserve exact matches

   do_fasts() has been modified to allow higher maxsav for do_walign
   than for do_work (2*nsegs, 6*nsegs)
 */

/* this code implements the "fasts" algorithm, which compares a set of
   protein fragments to a protein sequence.  Comma's are used to separate
   the sequence fragments, which need not be the same length.

   The expected input is:

   >mgstm1
   MGDAPDFD,
   MILGYW,
   MLLEYTDS

   The fragments do not need to be in the correct order (which is
   presumably unknown from the peptide sequencing.
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>

#include "defs.h"
#include "param.h"
#include "tatstats.h"

#define EOSEQ 0
#define ESS 49
#define NMAP_X 23 /* for 'X' */
#define NMAP_Z 24 /* for '*' */
#define MAXHASH 32
#define NMAP MAXHASH+1

static char *verstr="4.32 Feb 2007";

#define DROP_INTERN
#include "drop_func.h"

int shscore(const unsigned char *aa0, const int n0, int **pam2, int nsq);
static void update_code(char *al_str, int al_str_max, int op, int op_cnt, int fnum);
extern void aancpy(char *to, char *from, int count, struct pstruct pst);

#ifdef TFAST
extern int aatran(const unsigned char *ntseq, unsigned char *aaseq, const int maxs, const int frame);
#endif

void savemax(struct dstruct *, struct f_struct *, int maxsav, int exact,int t_end);

int spam(const unsigned char *, const unsigned char *, int, struct savestr *, int **, struct f_struct *);
int sconn(struct savestr **v,
	  int nsave,
	  struct f_struct *,
	  struct rstruct *,
	  struct pstruct *,
	  const unsigned char *aa0, int n0,
	  const unsigned char *aa1, int n1,
	  int opt_prob);

void kpsort(struct savestr **, int);
void kssort(struct savestr **, int);	/* sort by score */
int sconn_a(unsigned char *, int, 
	    const unsigned char *, int,
	    struct f_struct *, 
	    struct a_res_str *,
	    struct pstruct *);
void kpsort(struct savestr **, int);

/* initialize for fasta */

void
init_work (unsigned char *aa0, const int n0, 
	   struct pstruct *ppst,
	   struct f_struct **f_arg
	   )
{
   int mhv, phv;
   int hmax, nsegs;
   int i0, ib, hv, old_hv;
   int pamfact;
   struct f_struct *f_str;
   /* these used to be globals, but do not need to be */
   int ktup, fact, kt1;

   int maxn0;
   int stmp;	/* temporary score */
   int i, j, q;
   int tat_size;
   int *res;

   unsigned char *query;
   int k, l, m, n, N, length, index;

   double *tatprobptr;

   f_str = (struct f_struct *)calloc(1,sizeof(struct f_struct));

   ppst->param_u.fa.pgap = ppst->gdelval + ppst->ggapval;
   ktup = ppst->param_u.fa.ktup;
   if ( ktup  > 2 ) {
     ktup = ppst->param_u.fa.ktup = 2;
   }
   fact = ppst->param_u.fa.scfact;

   /* fasts3 cannot work with lowercase symbols as low complexity;
      thus, NMAP must be disabled; this depends on aascii['X']  */
   if (ppst->hsq[NMAP_X] == NMAP ) {ppst->hsq[NMAP_X]=1;}
   if (ppst->hsq[NMAP_Z] == NMAP ) {ppst->hsq[NMAP_Z]=1;}
   /* this does not work in a threaded environment */
   /*    else {fprintf(stderr," cannot find 'X'==NMAP\n");} */

   for (i0 = 1, mhv = -1; i0 <= ppst->nsq; i0++)
      if (ppst->hsq[i0] < NMAP && ppst->hsq[i0] > mhv)  mhv = ppst->hsq[i0];

   if (mhv <= 0) {
      fprintf (stderr, " maximum hsq <=0 %d\n", mhv);
      exit (1);
   }

   for (f_str->kshft = 0; mhv > 0; mhv /= 2) f_str->kshft++;

/*      kshft = 2;	*/
   kt1 = ktup-1;
   hv = 1;
   for (i0 = 0; i0 < ktup; i0++) hv = hv << f_str->kshft;
   hmax = hv;
   f_str->hmask = (hmax >> f_str->kshft) - 1;

   if ((f_str->aa0t = (unsigned char *) calloc(n0+1, sizeof(char))) == NULL) {
     fprintf (stderr, " cannot allocate f_str0->aa0t array; %d\n",n0+1);
     exit (1);
   }

   if ((f_str->aa0ti = (int *) calloc(n0+1, sizeof(int))) == NULL) {
     fprintf (stderr, " cannot allocate f_str0->aa0ti array; %d\n",n0+1);
     exit (1);
   }

   if ((f_str->aa0b = (int *) calloc(n0+1, sizeof(int))) == NULL) {
     fprintf (stderr, " cannot allocate f_str0->aa0b array; %d\n",n0+1);
     exit (1);
   }

   if ((f_str->aa0e = (int *) calloc(n0+1, sizeof(int))) == NULL) {
     fprintf (stderr, " cannot allocate f_str0->aa0e array; %d\n",n0+1);
     exit (1);
   }

   if ((f_str->aa0i = (int *) calloc(n0+1, sizeof(int))) == NULL) {
     fprintf (stderr, " cannot allocate f_str0->aa0i array; %d\n",n0+1);
     exit (1);
   }

   if ((f_str->aa0s = (int *) calloc(n0+1, sizeof(int))) == NULL) {
     fprintf (stderr, " cannot allocate f_str0->aa0s array; %d\n",n0+1);
     exit (1);
   }

   if ((f_str->aa0l = (int *) calloc(n0+1, sizeof(int))) == NULL) {
     fprintf (stderr, " cannot allocate f_str0->aa0l array; %d\n",n0+1);
     exit (1);
   }

   if ((f_str->harr = (int *) calloc (hmax, sizeof (int))) == NULL) {
     fprintf (stderr, " cannot allocate hash array: hmax: %d hmask: %d\n",
	      hmax, f_str->hmask);
     exit (1);
   }
   if ((f_str->pamh1 = (int *) calloc (ppst->nsq+1, sizeof (int))) == NULL) {
     fprintf (stderr, " cannot allocate pamh1 array\n");
     exit (1);
   }
   if ((f_str->pamh2 = (int *) calloc (hmax, sizeof (int))) == NULL) {
     fprintf (stderr, " cannot allocate pamh2 array\n");
     exit (1);
   }

   if ((f_str->link = (int *) calloc (n0, sizeof (int))) == NULL) {
     fprintf (stderr, " cannot allocate hash link array");
     exit (1);
   }

   /* for FASTS/FASTM, we want to know when we get to the end of a peptide,
      so we can ensure that we set the end and restart */

   if ((f_str->l_end = (int *) calloc (n0, sizeof (int))) == NULL) {
     fprintf (stderr, " cannot allocate link end array");
     exit (1);
   }

   for (i0 = 0; i0 < hmax; i0++) f_str->harr[i0] = -1;
   for (i0 = 0; i0 < n0; i0++) f_str->link[i0] = -1;
   for (i0 = 0; i0 < n0; i0++) f_str->l_end[i0] = 0;

   /* count the number of peptides */
   nsegs = 1;
   for (i0 = 0; i0 < n0; i0++) {
     if (aa0[i0] == ESS || aa0[i0] == 0) nsegs++;
   }

   /* allocate space for peptides offsets, nm_u */
   if ((f_str->nmoff = (int *)calloc(nsegs+1, sizeof(int)))==NULL) {
     fprintf(stderr, " cannot allocat nmoff array: %d\n", nsegs);
     exit(1);
   }

   if ((f_str->nm_u = (int *)calloc(nsegs+1, sizeof(int)))==NULL) {
     fprintf(stderr, " cannot allocat nm_u array: %d\n", nsegs);
     exit(1);
   }

   phv = hv = 0;
   f_str->nmoff[0] = 0;
   f_str->nm0 = 1;

   /* encode the aa0 array */
   if (kt1 > 0) {
     hv = ppst->hsq[aa0[0]];
     phv = ppst->pam2[0][aa0[0]][aa0[0]];
   }

   for (i0=kt1 ; i0 < n0; i0++) {
     if (aa0[i0] == ESS || aa0[i0] == 0) {
       /*       fprintf(stderr," converted %d to 0\n",aa0[i0]); */
       aa0[i0] = EOSEQ;	/* set ESS to 0 */
       f_str->nmoff[f_str->nm0++] = i0+1; 
       f_str->l_end[i0-1] = 1;
       phv = hv = 0;
       if (kt1 > 0) {
	 i0++;
	 hv = ppst->hsq[aa0[i0]];
	 phv = ppst->pam2[0][aa0[i0]][aa0[i0]];
       }
       continue;
     }

     hv = ((hv & f_str->hmask) << f_str->kshft) + ppst->hsq[aa0[i0]];
     f_str->link[i0] = f_str->harr[hv];
     f_str->harr[hv] = i0;
     f_str->pamh2[hv] = (phv += ppst->pam2[0][aa0[i0]][aa0[i0]]);
     phv -= ppst->pam2[0][aa0[i0 - kt1]][aa0[i0 - kt1]];
   }
   f_str->l_end[n0-1] = 1;

   f_str->nmoff[f_str->nm0] = n0+1;

   /*
#ifdef DEBUG
   fprintf(stderr, ">>%s\n",qtitle);
   for (j=0; j<f_str->nm0; j++) {
     for (i=f_str->nmoff[j]; i < f_str->nmoff[j+1]-1; i++) {
       fprintf(stderr,"%c",ppst->sq[aa0[i]]);
     }
     fprintf(stderr," %d\n",aa0[i]);
   }

   for (j=1; j<=ppst->nsq; j++) {
     fprintf(stderr, "%c %d\n", ppst->sq[j], f_str->harr[j]);
   }

   for (j=0; j<=n0; j++) {
     fprintf(stderr, "%c %d\n", ppst->sq[aa0[j]], f_str->link[j]);
   }

#endif
   */

   /* build an integer array of the max score that can be achieved
      from that position - use in savemax to mark some segments as
      fixed */

   /* setup aa0b[], aa0e[], which specify the begining and end of each
      segment */

   stmp = 0;
   q = -1;
   for (ib = i0 = 0; i0 < n0; i0++) {
     f_str->aa0l[i0] = i0 - q;
     if (aa0[i0]==EOSEQ) {
       f_str->aa0b[i0] = -1;
       f_str->aa0e[i0] = -1;
       f_str->aa0i[i0] = -1;
       f_str->aa0l[i0] = -1;
       q = i0;
       if (i0 > 0)f_str->aa0s[i0-1] = stmp;
       stmp = 0;
       ib++;
     }
     else {
       stmp += ppst->pam2[0][aa0[i0]][aa0[i0]];
     }

     f_str->aa0b[i0] =  f_str->nmoff[ib];
     f_str->aa0e[i0] =  f_str->nmoff[ib+1]-2;
     f_str->aa0i[i0] =  ib;

     /*
     fprintf(stderr,"%2d %c: %2d %2d %2d\n",i0,ppst->sq[aa0[i0]],
	     f_str->aa0b[i0],f_str->aa0e[i0],f_str->aa0i[i0]);
     */
   }
   f_str->aa0s[n0-1]=stmp;	/* save last best possible score */

   /* maxsav - maximum number of peptide alignments saved in search */
   /* maxsav_w - maximum number of peptide alignments saved in
      alignment */

   f_str->maxsav = max(MAXSAV,2*f_str->nm0);
   f_str->maxsav_w = max(MAXSAV,6*f_str->nm0);

   if ((f_str->vmax = (struct savestr *)
	calloc(f_str->maxsav_w,sizeof(struct savestr)))==NULL) {
     fprintf(stderr, "Couldn't allocate vmax[%d].\n",f_str->maxsav_w);
     exit(1);
   }

   if ((f_str->vptr = (struct savestr **)
	calloc(f_str->maxsav_w,sizeof(struct savestr *)))==NULL) {
     fprintf(stderr, "Couldn't allocate vptr[%d].\n",f_str->maxsav_w);
     exit(1);
   }

   if ((f_str->sarr = (struct slink *)
	calloc(f_str->maxsav_w,sizeof(struct slink)))==NULL) {
     fprintf(stderr, "Couldn't allocate sarr[%d].\n",f_str->maxsav_w);
     exit(1);
   }

   /* Tatusov Statistics Setup */

   /* initialize priors array. */
   if((f_str->priors = (double *)calloc(ppst->nsq+1, sizeof(double))) == NULL) {
     fprintf(stderr, "Couldn't allocate priors array.\n");
     exit(1);
   }

   calc_priors(f_str->priors, ppst, f_str, NULL, 0, ppst->pseudocts);

   /* pre-calculate the Tatusov probability array for each full segment */

   if(ppst->zsflag >= 1 && ppst->zsflag <= 3 && f_str->nm0 <= 10) {

     tat_size = (1<<f_str->nm0) -1;
     f_str->dotat = 1;
     f_str->tatprobs = (struct tat_str **) malloc((size_t)tat_size*sizeof(struct tat_str *));
     if (f_str->tatprobs == NULL) {
       fprintf (stderr, " cannot allocate tatprobs array: %ld\n",
		tat_size * sizeof(struct tat_str *));
       exit (1);
     }

     f_str->intprobs = (double **) malloc((size_t)tat_size * sizeof(double *));
     if(f_str->intprobs == NULL) {
       fprintf(stderr, "Couldn't allocate intprobs array.\n");
       exit(1);
     }

     for(k = 0, l = f_str->nm0 ; k < l ; k++) {
       query = &(aa0[f_str->nmoff[k]]);
       length = f_str->nmoff[k+1] - f_str->nmoff[k] - 1;

       /* this segment alone */
       index = (1 << k) - 1;
       generate_tatprobs(query, 0, length - 1, f_str->priors, ppst->pam2[0], ppst->nsq, &(f_str->tatprobs[index]), NULL);

       /* integrate the probabilities */
       N = f_str->tatprobs[index]->highscore - f_str->tatprobs[index]->lowscore;
       tatprobptr = (double *) calloc(N+1, sizeof(double));
       if(tatprobptr == NULL) {
	 fprintf(stderr, "Couldn't calloc tatprobptr.\n");
	 exit(1);
       }
       f_str->intprobs[index] = tatprobptr;

       for (i = 0; i <= N ; i++ ) {
	 tatprobptr[i] = f_str->tatprobs[index]->probs[i];
	 for (j = i + 1 ; j <= N ; j++ ) {
	   tatprobptr[i] += f_str->tatprobs[index]->probs[j];
	 }
       }

       /* this segment built on top of all other subcombinations */
       for(i = 0, j = (1 << k) - 1 ; i < j ; i++) {
	 index = (1 << k) + i;
	 generate_tatprobs(query, 0, length - 1, f_str->priors, ppst->pam2[0], ppst->nsq, &(f_str->tatprobs[index]), f_str->tatprobs[i]);

	 /* integrate the probabilities */
	 N = f_str->tatprobs[index]->highscore - f_str->tatprobs[index]->lowscore;
	 tatprobptr = (double *) calloc(N+1, sizeof(double));
	 if(tatprobptr == NULL) {
	   fprintf(stderr, "Couldn't calloc tatprobptr.\n");
	   exit(1);
	 }
	 f_str->intprobs[index] = tatprobptr;
       
	 for (m = 0; m <= N ; m++ ) {
	   tatprobptr[m] = f_str->tatprobs[index]->probs[m];
	   for (n = m + 1 ; n <= N ; n++ ) {
	     tatprobptr[m] += f_str->tatprobs[index]->probs[n];
	   }
	 }
       }
     }
   } else {
     f_str->dotat = 0;
     f_str->shuff_cnt = ppst->shuff_node;
   }

   /* End of Tatusov Statistics Setup */

   /*
   for (i0=1; i0<=ppst->nsq; i0++) {
     fprintf(stderr," %c: %2d ",ppst->sq[i0],f_str->harr[i0]);
     hv = f_str->harr[i0];
     while (hv >= 0) {
       fprintf(stderr," %2d",f_str->link[hv]);
       hv = f_str->link[hv];
     }
     fprintf(stderr,"\n");
   }
   */

/* this has been modified from 0..<ppst->nsq to 1..<=ppst->nsq because the
   pam2[0][0] is now undefined for consistency with blast
*/
   for (i0 = 1; i0 <= ppst->nsq; i0++)
     f_str->pamh1[i0] = ppst->pam2[0][i0][i0];

   f_str->ndo = 0;
   f_str->noff = n0-1;
   if (f_str->diag==NULL) 
     f_str->diag = (struct dstruct *) calloc ((size_t)MAXDIAG,
					      sizeof (struct dstruct));
   if (f_str->diag == NULL) {
      fprintf (stderr, " cannot allocate diagonal arrays: %ld\n",
	      (long) MAXDIAG * (long) (sizeof (struct dstruct)));
      exit (1);
   }

#ifdef TFAST
   if ((f_str->aa1x =(unsigned char *)calloc((size_t)ppst->maxlen+2,
					     sizeof(unsigned char)))
       == NULL) {
     fprintf (stderr, "cannot allocate aa1x array %d\n", ppst->maxlen+2);
     exit (1);
   }
   f_str->aa1x++;
#endif

   maxn0 = max(3*n0/2,MIN_RES);
   if ((res = (int *)calloc((size_t)maxn0,sizeof(int)))==NULL) {
     fprintf(stderr,"cannot allocate alignment results array %d\n",maxn0);
     exit(1);
   }
   f_str->res = res;
   f_str->max_res = maxn0;

   *f_arg = f_str;
}


/* pstring1 is a message to the manager, currently 512 */
/* pstring2 is the same information, but in a markx==10 format */
void
get_param (struct pstruct *pstr, char *pstring1, char *pstring2)
{
#ifdef FASTS
#ifndef TFAST
  char *pg_str="FASTS";
#else
  char *pg_str="TFASTS";
#endif
#endif

#ifdef FASTM
#ifndef TFAST
  char *pg_str="FASTM";
#else
  char *pg_str="TFASTM";
#endif
#endif

  sprintf (pstring1, "%s (%s) function [%s matrix (%d:%d)] ktup=%d",pg_str,verstr,
	       pstr->pamfile, pstr->pam_h,pstr->pam_l, pstr->param_u.fa.ktup);
  if (pstr->param_u.fa.iniflag) strcat(pstring1," init1");
  /*
  if (pstr->zsflag==0) strcat(pstring1," not-scaled");
  else if (pstr->zsflag==1) strcat(pstring1," reg.-scaled");
  */
  if (pstring2 != NULL) {
    sprintf (pstring2, "; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)\n\
; pg_gap-pen: %d %d\n; pg_ktup: %d\n",
	     pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l, pstr->gdelval,
	     pstr->ggapval,pstr->param_u.fa.ktup);
   }
}

void
close_work (const unsigned char *aa0, const int n0,
	    struct pstruct *ppst,
	    struct f_struct **f_arg)
{
  struct f_struct *f_str;
  int i, j;

  f_str = *f_arg;

  if (f_str != NULL) {

    free(f_str->res);
#ifdef TFAST
    free(f_str->aa1x - 1); /* because f_str->aa1x got ++'ed when allocated! */
#endif
    free(f_str->diag);
    free(f_str->l_end);
    free(f_str->link);
    free(f_str->pamh2); 
    free(f_str->pamh1);
    free(f_str->harr);
    free(f_str->vmax);
    free(f_str->vptr);
    free(f_str->sarr);
    free(f_str->aa0i);
    free(f_str->aa0e);
    free(f_str->aa0b);
    free(f_str->aa0ti);
    free(f_str->aa0t);
    free(f_str->nmoff);
    free(f_str->nm_u);

    if(f_str->dotat) {
      for(i = 0, j = (1 << f_str->nm0) - 1 ; i < j ; i++) {
	free(f_str->tatprobs[i]->probs);
	free(f_str->tatprobs[i]);
	free(f_str->intprobs[i]);
      }
      free(f_str->tatprobs);
      free(f_str->intprobs);
    }

    free(f_str->priors);
    free(f_str);
    *f_arg = NULL;
  }
}

void do_fasts (const unsigned char *aa0, const int n0,
	       const unsigned char *aa1, const int n1,
	       struct pstruct *ppst, struct f_struct *f_str,
	       struct rstruct *rst, int *hoff, int opt_prob,
	       int maxsav)
{
   int     nd;		/* diagonal array size */
   int     lhval;
   int     kfact;
   register struct dstruct *dptr;
   register int tscor;
   register struct dstruct *diagp;
   struct dstruct *dpmax;
   register int lpos;
   int     tpos;
   struct savestr *vmptr, *vmaxmax;
   int     scor, tmp;
   int     im, ib, nsave;
   int     cmps ();		/* comparison routine for ksort */
   int ktup;
   int doffset;


   vmaxmax = &f_str->vmax[maxsav];

   ktup = ppst->param_u.fa.ktup;

   if (n1 < ktup) {
     rst->score[0] = rst->score[1] = rst->score[2] = 0;
     rst->escore = 1.0;
     rst->segnum = 0;
     rst->seglen = 0;
     return;
   }

   if (n0+n1+1 >= MAXDIAG) {
     fprintf(stderr,"n0,n1 too large: %d, %d\n",n0,n1);
     rst->score[0] = rst->score[1] = rst->score[2] = -1;
     rst->escore = 2.0;
     rst->segnum = 0;
     rst->seglen = 0;
     return;
   }

   nd = n0 + n1;

   dpmax = &f_str->diag[nd];
   for (dptr = &f_str->diag[f_str->ndo]; dptr < dpmax;)
   {
      dptr->stop = -1;
      dptr->dmax = NULL;
      dptr++->score = 0;
   }

   for (vmptr = f_str->vmax; vmptr < vmaxmax; vmptr++) {
      vmptr->score = 0;
      vmptr->exact = 0;
   }
   f_str->lowmax = f_str->vmax;
   f_str->lowscor = 0;

   /* start hashing */
   diagp = &f_str->diag[f_str->noff];
   for (lhval=lpos=0; lpos < n1; lpos++, diagp++) {
     if (ppst->hsq[aa1[lpos]]>=NMAP) {	/* skip residue */
       lpos++ ; diagp++;
       while (lpos < n1 && ppst->hsq[aa1[lpos]]>=NMAP) {lpos++; diagp++;}
       if (lpos >= n1) break;
       lhval = 0;
     }

     lhval = ((lhval & f_str->hmask) << f_str->kshft) + ppst->hsq[aa1[lpos]];

     for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {

       dptr = &diagp[-tpos];

       if (f_str->l_end[tpos]) {
	 if (dptr->score + f_str->pamh1[aa0[tpos]] == f_str->aa0s[tpos]) {
	   dptr->stop = lpos;
	   dptr->score = f_str->aa0s[tpos];
	   savemax(dptr, f_str, maxsav, 1, tpos);
	   dptr->dmax = NULL;
	 }

	 else if (dptr->score + f_str->pamh1[aa0[tpos]] > f_str->aa0s[tpos]) {
	   /*
	   fprintf(stderr,"exact match score too high: %d:%d %d < %d + %d - %d:%d - %d > %d\n",
		   tpos, lpos, f_str->aa0s[tpos],dptr->score, f_str->pamh1[aa0[tpos]],
		   dptr->start, dptr->stop,
		   dptr->stop - dptr->start, f_str->aa0l[tpos]);
	   */
	   dptr->stop = lpos;
	   dptr->start = lpos - f_str->aa0l[tpos];
	   dptr->score = f_str->aa0s[tpos];
	   savemax(dptr, f_str, maxsav, 1, tpos);
	   dptr->dmax = NULL;
	 }
       }
       else if ((tscor = dptr->stop) >= 0) {
	 tscor++;	/* tscor is stop of current, increment it */
	 if ((tscor -= lpos) <= 0) {  /* tscor, the end of the current
					 match, is before lpos, so there
					 is a mismatch - this is also the
					 mismatch cost */
	   tscor *= 2;
	   scor = dptr->score;	/* save the run score on the diag */
	   if ((tscor += (kfact = f_str->pamh2[lhval])) < 0 
	       && f_str->lowscor < scor) {
	     /* if what we will get (tscor + kfact) is < 0 and the
		score is better than the worst savemax() score, save
		it */
	     savemax (dptr, f_str, maxsav,0,-1);
	   }

	   /* if extending is better than starting over, extend */
	   if ((tscor += scor) >= kfact) {
	     dptr->score = tscor;
	     dptr->stop = lpos;
	     if (f_str->l_end[tpos]) {
	       if (dptr->score == f_str->aa0s[tpos]) {
		 savemax(dptr, f_str, maxsav,1,tpos);
		 dptr->dmax = NULL;
	       }
	       else if (dptr->score > f_str->lowscor)
		 savemax(dptr, f_str, maxsav,0,tpos);
	     }
	   }
	   else {     /* otherwise, start new */
	     dptr->score = kfact;
	     dptr->start = dptr->stop = lpos;
	   }
	 } 
	 else { /* tscor is after lpos, so extend one residue */
	   dptr->score += f_str->pamh1[aa0[tpos]];
	   dptr->stop = lpos;
	   if (f_str->l_end[tpos]) {
	     if (dptr->score == f_str->aa0s[tpos]) {
	       savemax(dptr, f_str, maxsav,1,tpos);
	       dptr->dmax = NULL;
	     }
	     else if (dptr->score > f_str->lowscor)
	       savemax(dptr, f_str, maxsav,0,tpos);
	   }
	 }
       }
       else {	/* start new */
	 dptr->score = f_str->pamh2[lhval];
	 dptr->start = dptr->stop = lpos;
       }
     }				/* end tpos */
   }				/* end lpos */

   for (dptr = f_str->diag; dptr < dpmax;) {
     if (dptr->score > f_str->lowscor) savemax (dptr, f_str, maxsav,0,-1);
     dptr->stop = -1;
     dptr->dmax = NULL;
     dptr++->score = 0;
   }
   f_str->ndo = nd;

/*
        at this point all of the elements of aa1[lpos]
        have been searched for elements of aa0[tpos]
        with the results in diag[dpos]
*/

   for (nsave=0, vmptr=f_str->vmax; vmptr< vmaxmax; vmptr++) {
      if (vmptr->score > 0) {
	/*

	fprintf(stderr,"%c 0: %4d-%4d  1: %4d-%4d  dp: %d score: %d",
		(vmptr->exact ? 'x' : ' '),
		f_str->noff+vmptr->start-vmptr->dp,
		f_str->noff+vmptr->stop-vmptr->dp,
		vmptr->start,vmptr->stop,
		vmptr->dp,vmptr->score);
	*/
	vmptr->score = spam (aa0, aa1, n1, vmptr, ppst->pam2[0], f_str);
	/*
	fprintf(stderr,"  sscore: %d %d-%d\n",vmptr->score,vmptr->start,vmptr->stop);
	*/
	if (vmptr->score > 0) f_str->vptr[nsave++] = vmptr;
      }
   }

   if (nsave <= 0) {
     rst->score[0] = rst->score[1] = rst->score[2] = 0;
     rst->escore = 1.0;
     rst->segnum = 0;
     rst->seglen = 0;
     f_str->nsave = 0;
     return;
   }

   /*
   fprintf(stderr,"n0: %d; n1: %d; noff: %d\n",n0,n1,f_str->noff);
   for (ib=0; ib<nsave; ib++) {
     fprintf(stderr,"%c 0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
	     f_str->vptr[ib]->exact ? 'x' : ' ',
	     f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
	     f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
	     f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
	     f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
   }

   fprintf(stderr,"---\n");
   */
   kssort(f_str->vptr,nsave);

   /* make certain each seg is used only once */

   for (ib=0; ib<f_str->nm0; ib++) f_str->nm_u[ib]=0;
   for (ib=0; ib < nsave; ib++) {
     doffset = f_str->vptr[ib]->dp - f_str->noff;
     tpos=f_str->aa0i[f_str->vptr[ib]->start - doffset];
     if (f_str->nm_u[tpos] == 0) {
       f_str->nm_u[tpos]=1;
     } else {
       f_str->vptr[ib]->score = -1;
     }
   }

   kssort(f_str->vptr,nsave);
   for (ib = nsave-1; ib >= 0; ib--)
     if (f_str->vptr[ib]->score > -1) break;
   nsave = ib+1;

   scor = sconn (f_str->vptr, nsave, 
		 f_str, rst, ppst, aa0, n0, aa1, n1,
		 opt_prob);

   if (rst->escore < 0.0) rst->escore = 2.0;
   kssort(f_str->vptr,nsave);

   /* here we should use an nsave that is consistent with sconn and nm0 */

   f_str->nsave = nsave;
   if (nsave > f_str->nm0) f_str->nsave = f_str->nm0;

   rst->score[1] = f_str->vptr[0]->score;
   rst->score[0] = rst->score[2] = max(scor, f_str->vptr[0]->score);

}

void do_work (const unsigned char *aa0, const int n0,
	      const unsigned char *aa1, const int n1,
	      int frame,
	      struct pstruct *ppst, struct f_struct *f_str,
	      int qr_flg, struct rstruct *rst)
{
  int opt_prob;
  int hoff, n10, i;

  if (qr_flg==1 && f_str->shuff_cnt <= 0) {
    rst->escore = 2.0;
    rst->score[0]=rst->score[1]=rst->score[2]= -1;
    return;
  }

  if (f_str->dotat || ppst->zsflag == 4 || ppst->zsflag == 14 ) opt_prob=1;
  else opt_prob = 0;
  if (ppst->zsflag == 2 || ppst->zsflag == 12) opt_prob = 0;
  if (qr_flg) {
    opt_prob=1;
    /*    if (frame==1) */
      f_str->shuff_cnt--;
  }

  if (n1 < ppst->param_u.fa.ktup) {
    rst->score[0] = rst->score[1] = rst->score[2] = -1;
    rst->escore = 2.0;
    return;
  }
#ifdef TFAST
  n10=aatran(aa1,f_str->aa1x,n1,frame);
  if (ppst->debug_lib)
    for (i=0; i<n10; i++)
      if (f_str->aa1x[i]>ppst->nsq) {
	fprintf(stderr,
		"residue[%d/%d] %d range (%d)\n",i,n1,
		f_str->aa1x[i],ppst->nsq);
	f_str->aa1x[i]=0;
	n10=i-1;
      }

  do_fasts (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff, opt_prob, f_str->maxsav);
#else	/* FASTA */
  do_fasts (aa0, n0, aa1, n1, ppst, f_str, rst, &hoff, opt_prob, f_str->maxsav);
#endif

  rst->comp = rst->H = -1.0;
}

void do_opt (const unsigned char *aa0, const int n0,
	     const unsigned char *aa1, const int n1,
	     int frame,
	     struct pstruct *ppst, struct f_struct *f_str,
	     struct rstruct *rst)
{
  int lag, tscore, hoff, n10;

#ifdef TFAST
  n10=aatran(aa1,f_str->aa1x,n1,frame);
  do_fasts (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff, 1, f_str->maxsav);
#else	/* FASTA */
  do_fasts(aa0,n0,aa1,n1,ppst,f_str,rst, &hoff, 1, f_str->maxsav);
#endif
}


/* modify savemax() so that full length 100% matches are marked
   so that they cannot be removed - if we have a 100% match, mark "exact"

   modify savemax() to split alignments that include a comma
*/

/* savemax(dptr, f_str, maxsav) takes a current diagonal run (saved in dptr),
   and places it in the set of runs to be saved (in  f_str->vmax[])
*/

void 
savemax (struct dstruct *dptr, struct f_struct *f_str, int maxsav,
	 int exact, int tpos)
{
  register int dpos;	/* position along the diagonal, -n0 .. n1 */
  int i, j, lowj;
  register struct savestr *vmptr;
  struct savestr *vmaxmax;

  vmaxmax = &f_str->vmax[maxsav];

  dpos = (int) (dptr - f_str->diag);	/* current diagonal */

/* check to see if this is the continuation of a run that is already saved */
/* if we are at the end of the query, save it regardless */

/*  if (t_end > 0 && t_end < dptr->stop - dptr->start) {return;} */

  if ((vmptr = dptr->dmax) != NULL	/* have an active run */
      && vmptr->dp == dpos &&		/* on the correct diagonal */
      vmptr->start == dptr->start) {	/* and it starts at the same place */
    vmptr->stop = dptr->stop;	/* update the end of the match in vmax[] */

    if (exact == 1) {
    /*
      fprintf(stderr,"have cont exact match: %d - %d:%d %d:%d = %d\n",
	      dptr->score, dptr->start, dptr->stop,
	      vmptr->start, vmptr->stop, dptr->stop - dptr->start+1);
    */
      exact = 1;
    }


/* if the score is worse, don't update, return - if the score gets bad
   enough, it will restart in the diagonal scan */
    if ((i = dptr->score) <= vmptr->score) { return;} 

/* score is better, update */
    vmptr->score = i;

    vmptr->exact = exact;
/* if the score is not the worst, return */
    if (vmptr != f_str->lowmax) { return;}
  }
  else {	/* not a continuation */
    /* save in the lowest place */
    /*
    fprintf(stderr," Replacing: %d - %d:%d => %d - %d:%d",
	    f_str->lowmax->score, f_str->lowmax->start, f_str->lowmax->stop,
	    dptr->score, dptr->start, dptr->stop);
    */

    vmptr = f_str->lowmax;

    /*
    if (exact == 1) {
      fprintf(stderr,"have new exact match: %d - %d:%d = %d\n",
	      dptr->score, dptr->start, dptr->stop, dptr->stop - dptr->start+1);
    }
    */
    vmptr->exact = exact;

    i = vmptr->score = dptr->score;   /* 'i' is used as a bound */
    vmptr->dp = dpos;
    vmptr->start = dptr->start;
    vmptr->stop = dptr->stop;
    dptr->dmax = vmptr;
  }

  /* rescan the list for the worst score */
  for (vmptr = f_str->vmax;  vmptr < &f_str->vmax[maxsav] ; vmptr++) {
    if (vmptr->score < i && !vmptr->exact) {
      i = vmptr->score;
      f_str->lowmax = vmptr;
    }
  }

  f_str->lowscor = i;
}

/* this version of spam scans the diagonal to find the best local score,
   then resets the boundaries for a global alignment and re-scans */

/* NOOVERHANG allows one to score any overhanging alignment as zero.
   Useful for SAGE alignments.  Normally, one allows overhangs because
   of the possibility of partial sequences.
*/

#undef NOOVERHANG

/* 
   May, 2005 - spam() has an intesting bug that occurs when two
   peptides match in order, separated by one position (the comma).  In
   this case, spam() splits the match, and only returns the better of
   the two matches.  So, if spam splits an alignment at a comma, it
   needs the ability to insert the missing match.

*/

int spam (const unsigned char *aa0, const unsigned char *aa1,int n1,
	  struct savestr *dmax, int **pam2,
	  struct f_struct *f_str)
{
   int     lpos, doffset;
   int     tot, mtot;
   struct {
     int  start, stop, score;
   } curv, maxv;
   register const unsigned char *aa0p, *aa1p;

   doffset = dmax->dp - f_str->noff;
   curv.start = dmax->start;
   aa1p = &aa1[dmax->start];
   aa0p = &aa0[dmax->start - doffset];

   tot = curv.score = maxv.score = 0;
   for (lpos = dmax->start; lpos <= dmax->stop; lpos++) {
     tot += pam2[*aa0p++][*aa1p++];
     if (tot > curv.score) {
       curv.stop = lpos;	/* here, curv.stop is actually curv.max */
       curv.score = tot;
      }
      else if (tot < 0) {
	if (curv.score > maxv.score) {
	  maxv.start = curv.start;
	  maxv.stop = curv.stop;
	  maxv.score = curv.score;
	}
	tot = curv.score = 0;
	curv.start = lpos+1;
      }
   }

   if (curv.score > maxv.score) {
     maxv.start = curv.start;
     maxv.stop = curv.stop;
     maxv.score = curv.score;
   }

   if (maxv.score <= 0) return 0;

   /* now, reset the boundaries of the alignment using aa0b[]
      and aa0e[], which specify the residues that start and end
      the segment */
      
   maxv.start = f_str->aa0b[maxv.stop-doffset] + doffset;
   if (maxv.start < 0) {
     maxv.start = 0;
#ifdef NOOVERHANG
     return 0;
#endif
   }

   maxv.stop = f_str->aa0e[maxv.stop-doffset] + doffset;
   if (maxv.stop > n1) {
     maxv.stop = n1-1;
#ifdef NOOVERHANG
     return 0;
#endif
   }
   aa1p = &aa1[lpos = maxv.start];
   aa0p = &aa0[lpos - doffset];

   for (tot=0; lpos <= maxv.stop; lpos++) {
     tot += pam2[*aa0p++][*aa1p++];
   }

   maxv.score = tot;

/*	if (maxv.start != dmax->start || maxv.stop != dmax->stop)
		printf(" new region: %3d %3d %3d %3d\n",maxv.start,
			dmax->start,maxv.stop,dmax->stop);
*/
   dmax->start = maxv.start;
   dmax->stop = maxv.stop;

   return maxv.score;
}

int sconn (struct savestr **v, int n, 
	   struct f_struct *f_str,
	   struct rstruct *rst, struct pstruct *ppst,
	   const unsigned char *aa0, int n0,
	   const unsigned char *aa1, int n1, int opt_prob)
{
   int     i, si, cmpp ();
   struct slink *start, *sl, *sj, *so, *sarr;
   int     lstart, ltmp, tstart, plstop, ptstop, ptstart, tstop;
   double  tatprob;
   int     dotat;

   sarr = f_str->sarr;

   /*  sort the score left to right in lib pos */
   kpsort (v, n);

   start = NULL;
   rst->score[0] = 0;
   rst->escore = 2.0;

/*  for the remaining runs, see if they fit */
/*  lstart/lstop -> start/stop in library sequence
    tstart/tstop -> start/stop in query sequence
    plstart/plstop ->
*/

   for (i = 0, si = 0; i < n; i++) {

     /* the segment is worth adding; find out where? */
     lstart = v[i]->start;
     ltmp = v[i]->stop;
     tstart = lstart - v[i]->dp + f_str->noff;
     tstop = ltmp - v[i]->dp + f_str->noff;

     /*	put the run in the group */
     sarr[si].vp = v[i];
     sarr[si].score = v[i]->score;
     sarr[si].next = NULL;
     sarr[si].prev = NULL;
     sarr[si].tat = NULL;

/*
  opt_prob for FASTS only has to do with using aa1 for priors,
  i.e. we always calculate tatprobs for segments in FASTS (unlike
  FASTF)
*/
     if(opt_prob) {
       sarr[si].tatprob = 
	 calc_tatusov(NULL, &sarr[si], aa0, n0, aa1, n1, 
		      ppst->pam2[0], ppst->nsq, f_str, 
		      ppst->pseudocts, opt_prob, ppst->zsflag);
       if (sarr[si].tatprob < 0.0) {
	 fprintf(stderr," negative tatprob: %lg\n",sarr[si].tatprob);
	 sarr[si].tatprob = 1.0;
       }
       sarr[si].tat = sarr[si].newtat;
     }

/*  if it fits, then increase the score

    start points to the highest scoring run
    -> next is the second highest, etc.
    put the segment into the highest scoring run that it fits into
*/
     for (sl = start; sl != NULL; sl = sl->next) {
       ltmp = sl->vp->start;
 /* plstop -> previous lstop */
       plstop = sl->vp->stop;
 /* ptstart -> previous t(query) start */
       ptstart = ltmp - sl->vp->dp + f_str->noff;
 /* ptstop -> previous t(query) stop */
       ptstop = plstop - sl->vp->dp + f_str->noff;
#ifndef FASTM
 /* if the previous library stop is before the current library start */
       if (plstop < lstart && ( ptstop < tstart || ptstart > tstop))
#else
 /* if the previous library stop is before the current library start */
       if (plstop < lstart && ptstop < tstart)
#endif
       {
	 if(!opt_prob) {
	    sarr[si].score = sl->score + v[i]->score;
	    sarr[si].prev = sl;
	    break;
	  } else {
	    tatprob = calc_tatusov(sl, &sarr[si], aa0, n0, aa1, n1, 
				   ppst->pam2[0], ppst->nsq, f_str, 
				   ppst->pseudocts, opt_prob, ppst->zsflag);
	    /* if our tatprob gets worse when we add this, forget it */
	    if(tatprob > sarr[si].tatprob) {
	      free(sarr[si].newtat->probs); /* get rid of new tat struct */
	      free(sarr[si].newtat);
	      continue; /* reuse this sarr[si] */
	    } else {
	      sarr[si].tatprob = tatprob;
	      free(sarr[si].tat->probs); /* get rid of old tat struct */
	      free(sarr[si].tat);
	      sarr[si].tat = sarr[si].newtat;
	      sarr[si].prev = sl;
	      sarr[si].score = sl->score + v[i]->score;
	      /*
		fprintf(stderr,"sconn %d added %d/%d getting %d; si: %d, tat: %g\n",
		i,v[i]->start, v[i]->score,sarr[si].score,si, tatprob);
	      */
	      break;
	    }
	  }
	}
      }
      
      /* now recalculate where the score fits */
      if (start == NULL) start = &sarr[si];
      else {
	if(!opt_prob) {
	  for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
	    if (sarr[si].score > sj->score) {
	      sarr[si].next = sj;
	      if (so != NULL)
		so->next = &sarr[si];
	      else
		start = &sarr[si];
	      break;
	    }
	    so = sj;
	  }
	} else {
	  for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
	    if ( sarr[si].tatprob < sj->tatprob ||
		 ((sarr[si].tatprob == sj->tatprob) && sarr[si].score > sj->score) ) {
	      sarr[si].next = sj;
	      if (so != NULL)
		so->next = &sarr[si];
	      else
		start = &sarr[si];
	      break;
	    }
	    so = sj;
	  }
	}
      }

      si++;
   }
      
   if(opt_prob) {
     for (i = 0 ; i < si ; i++) {
       free(sarr[i].tat->probs);
       free(sarr[i].tat);
     }
   }

   if (start != NULL) {
     if(opt_prob) {
       rst->escore = start->tatprob;
     } else {
       rst->escore = 2.0;
     }

     rst->segnum = rst->seglen = 0;
     for(sj = start ; sj != NULL; sj = sj->prev) {
       rst->segnum++;
       rst->seglen += sj->vp->stop - sj->vp->start + 1;
     }
     return (start->score);
   } else {
     rst->escore = 1.0;
   }

   rst->segnum = rst->seglen = 0;
   return (0);
}

void
kssort (v, n)
struct savestr *v[];
int     n;
{
   int     gap, i, j;
   struct savestr *tmp;

   for (gap = n / 2; gap > 0; gap /= 2)
      for (i = gap; i < n; i++)
	 for (j = i - gap; j >= 0; j -= gap)
	 {
	    if (v[j]->score >= v[j + gap]->score)
	       break;
	    tmp = v[j];
	    v[j] = v[j + gap];
	    v[j + gap] = tmp;
	 }
}

void
kpsort (v, n)
struct savestr *v[];
int     n;
{
   int     gap, i, j;
   struct savestr *tmp;

   for (gap = n / 2; gap > 0; gap /= 2)
      for (i = gap; i < n; i++)
	 for (j = i - gap; j >= 0; j -= gap)
	 {
	    if (v[j]->start <= v[j + gap]->start)
	       break;
	    tmp = v[j];
	    v[j] = v[j + gap];
	    v[j + gap] = tmp;
	 }
}

/* calculate the 100% identical score */
int
shscore(const unsigned char *aa0, const int n0, int **pam2, int nsq)
{
  int i, sum;
  for (i=0,sum=0; i<n0; i++)
    if (aa0[i] != EOSEQ && aa0[i]<=nsq) sum += pam2[aa0[i]][aa0[i]];
  return sum;
}

/* sorts alignments from right to left (back to front) based on stop */

void
krsort (v, n)
struct savestr *v[];
int     n;
{
   int     gap, i, j;
   struct savestr *tmp;

   for (gap = n / 2; gap > 0; gap /= 2)
      for (i = gap; i < n; i++)
	 for (j = i - gap; j >= 0; j -= gap)
	 {
	    if (v[j]->stop > v[j + gap]->stop)
	       break;
	    tmp = v[j];
	    v[j] = v[j + gap];
	    v[j + gap] = tmp;
	 }
}

int  do_walign (const unsigned char *aa0, int n0,
		const unsigned char *aa1, int n1,
		int frame,
		struct pstruct *ppst, 
		struct f_struct *f_str, 
		struct a_res_str *a_res,
		int *have_ares)
{
  int hoff, n10;
  struct rstruct rst;
  int ib, i;
  unsigned char *aa0t;
  const unsigned char *aa1p;
  struct savestr *vmptr;

#ifdef TFAST
  f_str->n10 = n10 = aatran(aa1,f_str->aa1x,n1,frame);
  aa1p = f_str->aa1x;
#else
  n10 = n1;
  aa1p = aa1;
#endif

  do_fasts(aa0, n0, aa1p, n10, ppst, f_str, &rst, &hoff, 1, f_str->maxsav_w);

  /* the alignment portion takes advantage of the information left
     over in f_str after do_fasts is done.  in particular, it is
     easy to run a modified sconn() to produce the alignments.

     unfortunately, the alignment display routine wants to have
     things encoded as with bd_align and sw_align, so we need to do that.
     */

  /* unnecessary; do_fasts just did this */
  /*  kssort(f_str->vptr,f_str->nsave);  */

   /* at some point, we want one best score for each of the segments */

   for ( ; f_str->nsave > 0; f_str->nsave--) 
     if (f_str->vptr[f_str->nsave-1]->score >0) break;

   if ((aa0t = (unsigned char *)calloc(n0+1,sizeof(unsigned char)))==NULL) {
     fprintf(stderr," cannot allocate aa0t %d\n",n0+1);
     exit(1);
   }

   /* copy aa0[] into f_str->aa0t[] */
   for (i=0; i<n0; i++) f_str->aa0t[i] = aa0t[i] = aa0[i];
   f_str->aa0t[i] = aa0t[i] = '\0';

   a_res->nres = sconn_a (aa0t,n0,aa1p,n10,f_str, a_res, ppst);

   free(aa0t);

   a_res->res = f_str->res;
   *have_ares = 0;
   return rst.score[0];
}

/* this version of sconn is modified to provide alignment information */
/* in addition, it needs to know whether a segment has been used before */

/* sconn_a fills in the res[nres] array, but this is passed implicitly
   through f_str->res[f_str->nres] */

int sconn_a (unsigned char *aa0, int n0,
	     const unsigned char *aa1, int n1,
	     struct f_struct *f_str,
	     struct a_res_str *a_res,
	     struct pstruct *ppst)
{
   int     i, si, cmpp (), n;
   unsigned char *aa0p;
   int sx, dx, doff, *aa0tip;

   struct savestr **v;
   struct slink *start, *sl, *sj, *so, *sarr;
   int     lstart, lstop, ltmp, plstart, tstart, plstop, ptstop, ptstart, tstop;

   int *res, nres, tres;

   double tatprob;

/*	sort the score left to right in lib pos */

   v = f_str->vptr;
   n = f_str->nsave;
   sarr = f_str->sarr;

   /* set things up in case nothing fits */
   if (n <=0 || v[0]->score <= 0) return 0;

   if (v[0]->score < 0) {
     sarr[0].vp = v[0];
     sarr[0].score = v[0]->score;
     sarr[0].next = NULL;
     sarr[0].prev = NULL;
     start = &sarr[0];
   }
   else {

     krsort (v, n);	/* sort from left to right in library */

     start = NULL;

     /*	for each alignment, see if it fits */


     for (i = 0, si = 0; i < n; i++) {
       /*	if the score is less than the join threshold, skip it */

       if (v[i]->score < 0) continue;

       lstart = v[i]->start;
       lstop = v[i]->stop;
       tstart = lstart - v[i]->dp + f_str->noff;
       tstop = lstop - v[i]->dp + f_str->noff;

       /*	put the alignment in the group */

       sarr[si].vp = v[i];
       sarr[si].score = v[i]->score;
       sarr[si].next = NULL;
       sarr[si].prev = NULL;
       sarr[si].tat = NULL;

       sarr[si].tatprob = 
	 calc_tatusov(NULL, &sarr[si], aa0, n0, aa1, n1, 
		      ppst->pam2[0], ppst->nsq, f_str, 
		      ppst->pseudocts, 1, ppst->zsflag);
       sarr[si].tat = sarr[si].newtat;


       /* 	if it fits, then increase the score */
       /* start points to a sorted (by total score) list of candidate
	  overlaps */

       for (sl = start; sl != NULL; sl = sl->next) { 
	 plstart = sl->vp->start;
	 plstop = sl->vp->stop;
	 ptstart = plstart - sl->vp->dp + f_str->noff;
	 ptstop = plstop - sl->vp->dp + f_str->noff;
#ifndef FASTM
	 if (plstart > lstop && (ptstop < tstart || ptstart > tstop)) {
#else
         if (plstop > lstart && ptstart > tstop) {
#endif
	   /* alignment always uses probabilistic scoring ... */
	   /*   sarr[si].score = sl->score + v[i]->score;
		sarr[si].prev = sl;
		break; */		/* quit as soon as the alignment has been added */

	   tatprob = calc_tatusov(sl, &sarr[si], aa0, n0, aa1, n1, 
				  ppst->pam2[0], ppst->nsq, f_str, 
				  ppst->pseudocts, 1, ppst->zsflag);
	   /* if our tatprob gets worse when we add this, forget it */
	   if(tatprob > sarr[si].tatprob) {
	     free(sarr[si].newtat->probs); /* get rid of new tat struct */
	     free(sarr[si].newtat);
	     continue; /* reuse this sarr[si] */
	   } else {
	     sarr[si].tatprob = tatprob;
	     free(sarr[si].tat->probs); /* get rid of old tat struct */
	     free(sarr[si].tat);
	     sarr[si].tat = sarr[si].newtat;
	     sarr[si].prev = sl;
	     sarr[si].score = sl->score + v[i]->score;
	     /*
	       fprintf(stderr,"sconn %d added %d/%d getting %d; si: %d, tat: %g\n",
	       i,v[i]->start, v[i]->score,sarr[si].score,si, tatprob);
	     */
	     break;
	   }
	 }
       }

       /* now recalculate the list of best scores */
       if (start == NULL)
	 start = &sarr[si];	/* put the first one in the list */
       else
	 for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
	   /* if (sarr[si].score > sj->score) { */ /* new score better than old */
	   if ( sarr[si].tatprob < sj->tatprob ||
		((sarr[si].tatprob == sj->tatprob) && sarr[si].score > sj->score) ) {
	     sarr[si].next = sj;		/* next best after new score */
	     if (so != NULL)
	       so->next = &sarr[si];	/* prev_best->next points to best */
	     else  start = &sarr[si];	/* start points to best */
	     break;			/* stop looking */
	   }
	   so = sj;		/* previous candidate best */
	 }
       si++;				/* increment to next alignment */
     }
   }

   for (i = 0 ; i < si ; i++) {
     free(sarr[i].tat->probs);
     free(sarr[i].tat);
   }

   res = f_str->res;
   tres = nres = 0;
   aa0p = aa0;
   aa0tip = f_str->aa0ti;	/* point to temporary index */
   a_res->min1 = start->vp->start;
   a_res->min0 = 0;

   for (sj = start; sj != NULL; sj = sj->prev ) {
     doff = (int)(aa0p-aa0) - (sj->vp->start-sj->vp->dp+f_str->noff);
     
     /* fprintf(stderr,"doff: %3d\n",doff); */
     
     for (dx=sj->vp->start,sx=sj->vp->start-sj->vp->dp+f_str->noff;
	  dx <= sj->vp->stop; dx++) {
       *aa0tip++ = f_str->aa0i[sx];	/* save index */
       *aa0p++ = f_str->aa0t[sx++];	/* save sequence at index */
       tres++;
       res[nres++] = 0;
     }
     sj->vp->dp -= doff;
     if (sj->prev != NULL) {
       if (sj->prev->vp->start - sj->vp->stop - 1 > 0 )
	 tres += res[nres++] = (sj->prev->vp->start - sj->vp->stop - 1);
     }

     /*
     fprintf(stderr,"t0: %3d, tx: %3d, l0: %3d, lx: %3d, dp: %3d noff: %3d, score: %3d\n",
       sj->vp->start - sj->vp->dp + f_str->noff,
       sj->vp->stop - sj->vp->dp + f_str->noff,
       sj->vp->start,sj->vp->stop,sj->vp->dp,
       f_str->noff,sj->vp->score);

       fprintf(stderr,"%3d - %3d: %3d\n",
       sj->vp->start,sj->vp->stop,sj->vp->score);
     */
     a_res->max1 = sj->vp->stop+1;
     a_res->max0 = a_res->max1 - sj->vp->dp + f_str->noff;
   }

   /*
   fprintf(stderr,"(%3d - %3d):(%3d - %3d)\n",
	   a_res->min0,a_res->max0,a_res->min1,a_res->max1);
   */
   
   /* now replace f_str->aa0t with aa0
      (f_str->aa0t is permanent, aa0 is not)*/
   for (i=0; i<n0; i++) f_str->aa0t[i] = aa0[i];

   return tres;
}

/* for fasts (and fastf), pre_cons needs to set up f_str as well as do
   necessary translations - for right now, simply do do_walign */

void
pre_cons(const unsigned char *aa1, int n1, int frame, struct f_struct *f_str) {

#ifdef TFAST
  f_str->n10=aatran(aa1,f_str->aa1x,n1,frame);
#endif

}

/* aln_func_vals - set up aln.qlfact, qlrev, llfact, llmult, frame, llrev */
/* call from calcons, calc_id, calc_code */
void 
aln_func_vals(int frame, struct a_struct *aln) {

#ifdef TFAST
  aln->qlrev = 0;
  aln->qlfact= 1;
  aln->llfact = aln->llmult = 3;
  if (frame > 3) aln->llrev = 1;
  else aln->llrev = 0;
  aln->frame = 0;
#else	/* FASTS */
  aln->llfact = aln->llmult = aln->qlfact = 1;
  aln->llrev = aln->qlrev = 0;
  aln->frame = 0;
#endif
}

#include "a_mark.h"

int calcons(const unsigned char *aa0, int n0,
	    const unsigned char *aa1, int n1,
	    int *nc,
	    struct a_struct *aln,
	    struct a_res_str a_res,
	    struct pstruct pst,
	    char *seqc0, char *seqc1, char *seqca,
	    struct f_struct *f_str)
{
  int i0, i1, nn1, n0t;
  int op, lenc, len_gap, nd, ns, itmp;
  const unsigned char *aa1p;
  char *sp0, *sp1, *spa;
  int *rp;
  int mins, smins;
  
#ifndef TFAST
  aa1p = aa1;
  nn1 = n1;
#else
  aa1p = f_str->aa1x;
  nn1 = f_str->n10;
#endif

  aln->amin0 = a_res.min0;
  aln->amin1 = a_res.min1;
  aln->amax0 = a_res.max0;
  aln->amax1 = a_res.max1;

  /* first fill in the ends */
  n0 -= (f_str->nm0-1);

  if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1)
    			/* will we show all the start ?*/
    if (a_res.min0>=a_res.min1) {                        /* aa0 extends more to left */
      smins=0;
      if (aln->showall==1) mins=a_res.min0;
      else mins = min(a_res.min0,aln->llen/2);
      aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
      aln->smin0 = a_res.min0-mins;
      if ((mins-a_res.min1)>0) {
	memset(seqc1,' ',mins-a_res.min1);
	aancpy(seqc1+mins-a_res.min1,(char *)aa1p,a_res.min1,pst);
	aln->smin1 = 0;
      }
      else {
	aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
	aln->smin1 = a_res.min1-mins;
      }
    }
    else {
      smins=0;
      if (aln->showall == 1) mins=a_res.min1;
      else mins = min(a_res.min1,aln->llen/2);
      aancpy(seqc1,(char *)(aa1p+a_res.min1-mins),mins,pst);
      aln->smin1 = a_res.min1-mins;
      if ((mins-a_res.min0)>0) {
	memset(seqc0,' ',mins-a_res.min0);
	aancpy(seqc0+mins-a_res.min0,(char *)f_str->aa0t,a_res.min0,pst);
	aln->smin0 = 0;
      }
      else {
	aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
	aln->smin0 = a_res.min0-mins;
      }
    }
  else {
    mins= min(aln->llen/2,min(a_res.min0,a_res.min1));
    smins=mins;
    aln->smin0=a_res.min0;
    aln->smin1=a_res.min1;
    aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
    aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
  }

  memset(seqca,M_BLANK,mins);

/* now get the middle */

  spa = seqca+mins;
  sp0 = seqc0+mins;
  sp1 = seqc1+mins;
  rp = a_res.res;
  n0t = lenc = len_gap = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = op = 0;
  i0 = a_res.min0;
  i1 = a_res.min1;
  
  /* op is the previous "match/insert" operator; *rp is the current
     operator or repeat count */

  while (i0 < a_res.max0 || i1 < a_res.max1) {
    if (op == 0 && *rp == 0) {	/* previous was match (or start), current is match */
      op = *rp++;		/* get the next match/insert operator */

      /* get the alignment symbol */
      if ((itmp=pst.pam2[0][f_str->aa0t[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
      else if (itmp == 0) { *spa = M_ZERO;}
      else {*spa = M_POS;}
      if (*spa == M_ZERO || *spa == M_POS) { aln->nsim++;}

      *sp0 = pst.sq[f_str->aa0t[i0++]];	/* get the residues for the consensus */
      *sp1 = pst.sq[aa1p[i1++]];
      n0t++;
      lenc++;
      if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
      sp0++; sp1++; spa++;
    }
    else {	/* either op != 0 (previous was insert) or *rp != 0
		   (current is insert) */
      if (op==0) { op = *rp++;}	/* previous was match, start insert */
      				/* previous was insert - count through gap */
      *sp0++ = '-';
      *sp1++ = pst.sq[aa1p[i1++]];
      *spa++ = M_DEL;
      op--;
      len_gap++;
      lenc++;
    }
  }

  *spa = '\0';
  *nc = lenc-len_gap;
/*	now we have the middle, get the right end */

  ns = mins + lenc + aln->llen;
  ns -= (itmp = ns %aln->llen);
  if (itmp>aln->llen/2) ns += aln->llen;
  nd = ns - (mins+lenc);
  if (nd > max(n0t-a_res.max0,nn1-a_res.max1)) nd = max(n0t-a_res.max0,nn1-a_res.max1);
  
  if (aln->showall==1) {
    nd = max(n0t-a_res.max0,nn1-a_res.max1);	/* reset for showall=1 */
    /* get right end */
    aancpy(seqc0+mins+lenc,(char *)f_str->aa0t+a_res.max0,n0t-a_res.max0,pst);
    aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
    /* fill with blanks - this is required to use one 'nc' */
    memset(seqc0+mins+lenc+n0t-a_res.max0,' ',nd-(n0t-a_res.max0));
    memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
  }
  else {
    if ((nd-(n0t-a_res.max0))>0) {
      aancpy(seqc0+mins+lenc,(char *)f_str->aa0t+a_res.max0,
	     n0t-a_res.max0,pst);
      memset(seqc0+mins+lenc+n0t-a_res.max0,' ',nd-(n0t-a_res.max0));
    }
    else aancpy(seqc0+mins+lenc,(char *)f_str->aa0t+a_res.max0,nd,pst);
    if ((nd-(nn1-a_res.max1))>0) {
      aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
      memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
    }
    else aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nd,pst);
  }
  
  return mins+lenc+nd;
}

int
calcons_a(const unsigned char *aa0, unsigned char *aa0a, int n0,
	  const unsigned char *aa1, int n1,
	  int *nc,
	  struct a_struct *aln,
	  struct a_res_str a_res,
	  struct pstruct pst,
	  char *seqc0, char *seqc0a, char *seqc1, char *seqca,
	  char *ann_arr, struct f_struct *f_str)
{
  int i0, i1, nn1, n0t;
  int op, lenc, len_gap, nd, ns, itmp, p_ac, fnum, o_fnum;
  const unsigned char *aa1p;
  unsigned char *aa0ap;
  char *sp0, *sp0a, *sp1, *spa;
  int *rp;
  int mins, smins;
  
#ifndef TFAST
  aa1p = aa1;
  nn1 = n1;
#else
  aa1p = f_str->aa1x;
  nn1 = f_str->n10;
#endif

  aln->amin0 = a_res.min0;
  aln->amin1 = a_res.min1;
  aln->amax0 = a_res.max0;
  aln->amax1 = a_res.max1;

  /* first fill in the ends */
  n0 -= (f_str->nm0-1);

  if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1)
    			/* will we show all the start ?*/
    if (a_res.min0>=a_res.min1) {                        /* aa0 extends more to left */
      smins=0;
      if (aln->showall==1) mins=a_res.min0;
      else mins = min(a_res.min0,aln->llen/2);
      aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
      aln->smin0 = a_res.min0-mins;
      if ((mins-a_res.min1)>0) {
	memset(seqc1,' ',mins-a_res.min1);
	aancpy(seqc1+mins-a_res.min1,(char *)aa1p,a_res.min1,pst);
	aln->smin1 = 0;
      }
      else {
	aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
	aln->smin1 = a_res.min1-mins;
      }
    }
    else {
      smins=0;
      if (aln->showall == 1) mins=a_res.min1;
      else mins = min(a_res.min1,aln->llen/2);
      aancpy(seqc1,(char *)(aa1p+a_res.min1-mins),mins,pst);
      aln->smin1 = a_res.min1-mins;
      if ((mins-a_res.min0)>0) {
	memset(seqc0,' ',mins-a_res.min0);
	aancpy(seqc0+mins-a_res.min0,(char *)f_str->aa0t,a_res.min0,pst);
	aln->smin0 = 0;
      }
      else {
	aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
	aln->smin0 = a_res.min0-mins;
      }
    }
  else {
    mins= min(aln->llen/2,min(a_res.min0,a_res.min1));
    smins=mins;
    aln->smin0=a_res.min0;
    aln->smin1=a_res.min1;
    aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
    aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
  }

  memset(seqca,M_BLANK,mins);
  memset(seqc0a,' ', mins);

/* now get the middle */

  spa = seqca+mins;
  sp0 = seqc0+mins;
  sp0a = seqc0a+mins;
  sp1 = seqc1+mins;
  rp = a_res.res;
  n0t=lenc=len_gap=aln->nident=aln->nsim=aln->ngap_q=aln->ngap_l=op=p_ac= 0;
  i0 = a_res.min0;
  i1 = a_res.min1;
  
  /* op is the previous "match/insert" operator; *rp is the current
     operator or repeat count */

  o_fnum = f_str->aa0ti[i0];
  aa0ap = &aa0a[f_str->nmoff[o_fnum]+i0];

  while (i0 < a_res.max0 || i1 < a_res.max1) {
    fnum = f_str->aa0ti[i0];
    if (op == 0 && *rp == 0) {	/* previous was match (or start), current is match */
      if (p_ac == 0) { /* previous code was a match */
	if (fnum != o_fnum) { /* continuing a match, but with a different fragment */
	  aa0ap = &aa0a[f_str->nmoff[fnum]];
	  o_fnum = fnum;
	}
      }
      else {
	p_ac = 0; o_fnum = fnum = f_str->aa0ti[i0];
	aa0ap = &aa0a[f_str->nmoff[fnum]];
      }
      op = *rp++;		/* get the next match/insert operator */

      /* get the alignment symbol */
      if ((itmp=pst.pam2[0][f_str->aa0t[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
      else if (itmp == 0) { *spa = M_ZERO;}
      else {*spa = M_POS;}
      if (*spa == M_ZERO || *spa == M_POS) { aln->nsim++;}

      *sp0 = pst.sq[f_str->aa0t[i0++]];	/* get the residues for the consensus */
      *sp0a++ = ann_arr[*aa0ap++];
      *sp1 = pst.sq[aa1p[i1++]];
      n0t++;
      lenc++;
      if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
      sp0++; sp1++; spa++;
    }
    else {	/* either op != 0 (previous was insert) or *rp != 0
		   (current is insert) */
      if (op==0) { op = *rp++;}	/* previous was match, start insert */
      				/* previous was insert - count through gap */
      if (p_ac != 1) {
	p_ac = 1; fnum = f_str->aa0ti[i0];
      }

      *sp0++ = '-';
      *sp1++ = pst.sq[aa1p[i1++]];
      *spa++ = M_DEL;
      *sp0a++ = ' ';
      op--;
      len_gap++;
      lenc++;
    }
  }

  *sp0a = *spa = '\0';
  *nc = lenc-len_gap;
/*	now we have the middle, get the right end */

  ns = mins + lenc + aln->llen;
  ns -= (itmp = ns %aln->llen);
  if (itmp>aln->llen/2) ns += aln->llen;
  nd = ns - (mins+lenc);
  if (nd > max(n0t-a_res.max0,nn1-a_res.max1)) nd = max(n0t-a_res.max0,nn1-a_res.max1);
  
  if (aln->showall==1) {
    nd = max(n0t-a_res.max0,nn1-a_res.max1);	/* reset for showall=1 */
    /* get right end */
    aancpy(seqc0+mins+lenc,(char *)f_str->aa0t+a_res.max0,n0t-a_res.max0,pst);
    aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
    /* fill with blanks - this is required to use one 'nc' */
    memset(seqc0+mins+lenc+n0t-a_res.max0,' ',nd-(n0t-a_res.max0));
    memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
  }
  else {
    if ((nd-(n0t-a_res.max0))>0) {
      aancpy(seqc0+mins+lenc,(char *)f_str->aa0t+a_res.max0,
	     n0t-a_res.max0,pst);
      memset(seqc0+mins+lenc+n0t-a_res.max0,' ',nd-(n0t-a_res.max0));
    }
    else aancpy(seqc0+mins+lenc,(char *)f_str->aa0t+a_res.max0,nd,pst);
    if ((nd-(nn1-a_res.max1))>0) {
      aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
      memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
    }
    else aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nd,pst);
  }
  return mins+lenc+nd;
}

void aaptrshuffle(unsigned char *res, int n) {

  int i, j;
  unsigned char tmp;

  for( i = n; --i; ) {

    /* j = nrand(i); if (i == j) continue; */ /* shuffle */
    j = (n - 1) - i; if (i <= j ) break; /* reverse */

    tmp = res[i];
    res[i] = res[j];
    res[j] = tmp;
  }
}

void aa0shuffle(unsigned char *aa0, int n0, struct f_struct *f_str) {

  int i;
  int j;

  for(i = 0 ; i < f_str->nm0 ; i++) { /* for each fragment */

    aaptrshuffle(&(aa0[f_str->nmoff[i]]), 
		 f_str->nmoff[i+1] - f_str->nmoff[i] - 1 );

  }

}

/* build an array of match/ins/del - length strings */
int
calc_code(const unsigned char *aa0, const int n0,
	  const unsigned char *aa1, const int n1,
	  struct a_struct *aln,
	  struct a_res_str a_res,
	  struct pstruct pst,
	  char *al_str, int al_str_n, struct f_struct *f_str)
{
  int i0, i1, nn1;
  int op, lenc, len_gap;
  int p_ac, op_cnt;
  const unsigned char *aa1p;
  char tmp_cnt[20];
  char sp0, sp1, *sq;
  int *rp;
  int mins, smins;
  int o_fnum,fnum = 0;

  if (pst.ext_sq_set) {sq = pst.sqx;}
  else {sq = pst.sq;}

#ifndef TFAST
  aa1p = aa1;
  nn1 = n1;
#else
  aa1p = f_str->aa1x;
  nn1 = f_str->n10;
#endif

  aln->amin0 = a_res.min0;
  aln->amin1 = a_res.min1;
  aln->amax0 = a_res.max0;
  aln->amax1 = a_res.max1;

  rp = a_res.res;
  lenc = len_gap =aln->nident=aln->nsim=aln->ngap_q=aln->ngap_l=aln->nfs=op=p_ac = 0;
  op_cnt = 0;

  i0 = a_res.min0;	/* start in aa0 (f_str->aa0t) */
  i1 = a_res.min1;	/* start in aa1 */
  tmp_cnt[0]='\0';
  
  o_fnum = f_str->aa0ti[i0] + 1;	/* fragment number */
  while (i0 < a_res.max0 || i1 < a_res.max1) {
    fnum = f_str->aa0ti[i0]+1;
    if (op == 0 && *rp == 0) {	/* previous was match, this is match */
      if (p_ac == 0) {	/* previous code was a match */
	if (fnum == o_fnum) { op_cnt++;}
	else {		/* continuing a match, but with a different fragment */
	  update_code(al_str,al_str_n-strlen(al_str), p_ac, op_cnt, o_fnum);
	  o_fnum = fnum;
	  op_cnt=1;
	}
      }
      else {
	update_code(al_str,al_str_n-strlen(al_str),p_ac,op_cnt,o_fnum);
	op_cnt = 1; p_ac = 0; o_fnum = fnum = f_str->aa0ti[i0] + 1;
      }
      op = *rp++;
      lenc++;
      if (pst.pam2[0][f_str->aa0t[i0]][aa1p[i1]]>=0) {aln->nsim++;}
      sp0 = pst.sq[f_str->aa0t[i0++]];
      sp1 = pst.sq[aa1p[i1++]];
      if (toupper(sp0) == toupper(sp1)) aln->nident++;
    }
    else {
      if (op==0) op = *rp++;
      if (p_ac == 1) { op_cnt++;}
      else {
	update_code(al_str,al_str_n - strlen(al_str),p_ac,op_cnt,o_fnum);
	op_cnt = 1; p_ac = 1; fnum = f_str->aa0ti[i0] + 1;
      }
      op--; lenc++; i1++; len_gap++;
    }
  }
  update_code(al_str,al_str_n - strlen(al_str),p_ac,op_cnt,o_fnum);

  return lenc - len_gap;
}

/* update_code(): if "op" == 0, this is the end of a match of length
   "op_cnt" involving fragment "fnum"
   otherwise, this is an insertion (op==1) or deletion (op==2)
*/

void
update_code(char *al_str, int al_str_max, int op, int op_cnt, int fnum) {

  char op_char[4]={"=-+"};
  char tmp_cnt[20];

  if (op == 0)
    sprintf(tmp_cnt,"%c%d[%d]",op_char[op],op_cnt,fnum);
  else
    sprintf(tmp_cnt,"%c%d",op_char[op],op_cnt);

  strncat(al_str,tmp_cnt,al_str_max);
}

int
calc_id(const unsigned char *aa0, int n0,
	const unsigned char *aa1, int n1,
	struct a_struct *aln,
	struct a_res_str a_res,
	struct pstruct pst,
	struct f_struct *f_str)
{
  int i0, i1, nn1;
  int op, lenc, len_gap;
  const unsigned char *aa1p;
  int sp0, sp1;
  int *rp;
  int mins, smins;
  
#ifndef TFAST
  aa1p = aa1;
  nn1 = n1;
#else
  aa1p = f_str->aa1x;
  nn1 = f_str->n10;
#endif

  aln->amin0 = a_res.min0;
  aln->amin1 = a_res.min1;
  aln->amax0 = a_res.max0;
  aln->amax1 = a_res.max1;

  /* first fill in the ends */
  n0 -= (f_str->nm0-1);

  /* now get the middle */
  rp = a_res.res;
  lenc=len_gap=aln->nident=aln->nsim=aln->ngap_q = aln->ngap_l = aln->nfs = op = 0;
  i0 = a_res.min0;
  i1 = a_res.min1;
  
  while (i0 < a_res.max0 || i1 < a_res.max1) {
    if (op == 0 && *rp == 0) {
      op = *rp++;

      if (pst.pam2[0][f_str->aa0t[i0]][aa1p[i1]]>=0) {aln->nsim++;}

      sp0 = pst.sq[f_str->aa0t[i0++]];
      sp1 = pst.sq[aa1p[i1++]];
      lenc++;
      if (toupper(sp0) == toupper(sp1)) aln->nident++;
    }
    else {
      if (op==0) { op = *rp++;}
      i1++;
      op--;
      len_gap++;
      lenc++;
    }
  }
  return lenc-len_gap;
}

#ifdef PCOMPLIB

#include "structs.h"
#include "p_mw.h"

void
update_params(struct qmng_str *qm_msg,
	      struct mngmsg *m_msg, struct pstruct *ppst)
{
  m_msg->n0 = ppst->n0 = qm_msg->n0;
  m_msg->nm0 = qm_msg->nm0;
  m_msg->escore_flg = qm_msg->escore_flg;
  m_msg->qshuffle = qm_msg->qshuffle;
}
#endif