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

/*  - dropffa.c,v 1.1.1.1 1999/10/22 20:55:59 wrp Exp */

/* this code implements the "fastf" algorithm, which is designed to
   deconvolve mixtures of protein sequences derived from mixed-peptide
   Edman sequencing.  The expected input is:

   >test | 40001 90043 | mgstm1
   MGCEN,
   MIDYP,
   MLLAY,
   MLLGY

   Where the ','s indicate the length/end of the sequencing cycle
   data.  Thus, in this example, the sequence is from a mixture of 4
   peptides, M was found in the first position, G,I, and L(2) at the second,
   C,D, L(2) at the third, etc.

   Because the sequences are derived from mixtures, there need not be
   any partial sequence "MGCEN", the actual deconvolved sequence might be
   "MLDGN".
*/

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

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

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

#ifndef MAXSAV
#define MAXSAV 10
#endif

#define DROP_INTERN
#include "drop_func.h"

static char *verstr="4.21 May 2006 (ajm/wrp)";

int shscore(unsigned char *aa0, const int n0, int **pam2, int nsq);
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

struct hlstr { int next, pos;};

void savemax(struct dstruct *, struct f_struct *);

static int m0_spam(unsigned char *, const unsigned char *, int, struct savestr *,
	    int **, struct f_struct *);
static int m1_spam(unsigned char *, int,
		   const unsigned char *, int,
		   struct savestr *, int **, int, struct f_struct *);

int sconn(struct savestr **v, int nsave, int cgap,
	  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);
void kpsort(struct savestr **, int);

int
sconn_a(unsigned char *, int, int, struct f_struct *,
	struct a_res_str *);

/* initialize for fasta */

void
init_work (unsigned char *aa0, int n0, 
	   struct pstruct *ppst,
	   struct f_struct **f_arg)
{
   int mhv, phv;
   int hmax;
   int i0, ii0, hv;
   struct f_struct *f_str;

   int maxn0;
   int i, j, q;
   struct savestr *vmptr;
   int *res;

   f_str = (struct f_struct *) calloc(1, sizeof(struct f_struct));
   if(f_str == NULL) {
     fprintf(stderr, "Couldn't calloc f_str\n");
     exit(1);
   }

   ppst->sw_flag = 0;

   /* fastf3 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 for share ppst structs, as in threads */
   /*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;	*/
   hmax = hv = (1 << f_str->kshft);
   f_str->hmask = (hmax >> f_str->kshft) - 1;

   if ((f_str->aa0 = (unsigned char *) calloc(n0+1, sizeof(char))) == NULL) {
     fprintf (stderr, " cannot allocate f_str->aa0 array; %d\n",n0+1);
     exit (1);
   }
   for (i=0; i<n0; i++) f_str->aa0[i] = aa0[i];
   aa0 = f_str->aa0;

   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);
   }
   f_str->aa0ix = 0;

   if ((f_str->harr = (struct hlstr *) calloc (hmax, sizeof (struct hlstr))) == 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 = (struct hlstr *) calloc (n0, sizeof (struct hlstr))) == NULL) {
     fprintf (stderr, " cannot allocate hash link array");
     exit (1);
   }

   for (i0 = 0; i0 < hmax; i0++) {
      f_str->harr[i0].next = -1;
      f_str->harr[i0].pos = -1;
   }

   for (i0 = 0; i0 < n0; i0++) {
      f_str->link[i0].next = -1;
      f_str->link[i0].pos = -1;
   }

   /* encode the aa0 array */
   /*
     this code has been modified to allow for mixed peptide sequences
      aa0[] = 5 8 9 3 4 NULL 5 12 3 7 2 NULL
      the 'NULL' character resets the hash position counter, to indicate that
      any of several residues can be in the same position.
      We also need to keep track of the number of times this has happened, so that
      we can redivide the sequence later

      i0 counts through the sequence
      ii0 counts through the hashed sequence

      */

   f_str->nm0 = 1;
   f_str->nmoff = -1;
   phv = hv = 0;
   for (i0= ii0 = 0; i0 < n0; i0++, ii0++) {
     /* reset the counter and start hashing again */
     if (aa0[i0] == ESS || aa0[i0] == 0) {
       aa0[i0] = 0;	/* set ESS to 0 */
       /*       fprintf(stderr," converted ',' to 0\n");*/
       i0++;	/* skip over the blank */
       f_str->nm0++;
       if (f_str->nmoff < 0) f_str->nmoff = i0;
       phv = hv = 0;
       ii0 = 0;
     }
     hv = ppst->hsq[aa0[i0]];
     f_str->link[i0].next = f_str->harr[hv].next;
     f_str->link[i0].pos = f_str->harr[hv].pos;
     f_str->harr[hv].next = i0;
     f_str->harr[hv].pos = ii0;
     f_str->pamh2[hv] = ppst->pam2[0][aa0[i0]][aa0[i0]];
   }
   if (f_str-> nmoff < 0) f_str->nmoff = n0;


#ifdef DEBUG
   /*
   fprintf(stderr," nmoff: %d/%d nm0: %d\n", f_str->nmoff, n0,f_str->nm0);
   */
#endif

/*
#ifdef DEBUG
   fprintf(stderr," hmax: %d\n",hmax);
   for ( hv=0; hv<hmax; hv++)
       fprintf(stderr,"%2d %c %3d %3d\n",hv,
	       (hv > 0 && hv < ppst->nsq ) ? ppst->sq[ppst->hsq[hv]] : ' ',
	       f_str->harr[hv].pos,f_str->harr[hv].next);
   fprintf(stderr,"----\n");
   for ( hv=0; hv<n0; hv++)
       fprintf(stderr,"%2d: %3d %3d\n",hv,
	       f_str->link[hv].pos,f_str->link[hv].next);
#endif
*/

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

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

   for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++) {
     vmptr->used = (int *) calloc(n0, sizeof(int));
     if(vmptr->used == NULL) {
       fprintf(stderr, "Couldn't alloc vmptr->used\n");
       exit(1);
     }
   }

/* 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];

   ppst->param_u.fa.cgap = shscore(aa0,f_str->nmoff-1,ppst->pam2[0],ppst->nsq)/3;
   if (ppst->param_u.fa.cgap > ppst->param_u.fa.bestmax/4)
     ppst->param_u.fa.cgap = ppst->param_u.fa.bestmax/4;

   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

   /* allocate space for the scoring arrays */
   maxn0 = n0 + 4;

   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;

   /* 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);

   f_str->dotat = 0;
   f_str->shuff_cnt = ppst->shuff_node;

   /* End of Tatusov Statistics Setup */

   *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)
{
#ifndef TFAST
  char *pg_str="FASTF";
#else
  char *pg_str="TFASTF";
#endif

  sprintf (pstring1, "%s (%s) function [%s matrix (%d:%d)] join: %d",pg_str,verstr,
	   pstr->pamfile, pstr->pam_h,pstr->pam_l,pstr->param_u.fa.cgap);

  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_join: %d\n",
	     pg_str,verstr, pstr->pamfile, pstr->pam_h,pstr->pam_l,
	     pstr->param_u.fa.cgap);
   }
}

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

  f_str = *f_arg;

  if (f_str != NULL) {

    for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
      free(vmptr->used);

    free(f_str->res);
#ifdef TFAST
    free(f_str->aa1x - 1); /* allocated, then aa1x++'ed */
#endif
    free(f_str->diag);
    free(f_str->link);
    free(f_str->pamh2); 
    free(f_str->pamh1);
    free(f_str->harr);
    free(f_str->aa0t);
    free(f_str->aa0);
    free(f_str->priors);
    free(f_str);
    *f_arg = NULL;
  }
}

int do_fastf (unsigned char *aa0, int n0,
	      const unsigned char *aa1, int n1,
	      struct pstruct *ppst, struct f_struct *f_str,
	      struct rstruct *rst, int *hoff, int opt_prob)
{
   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, npos;
   struct savestr *vmptr;
   int     scor, tmp;
   int     im, ib, nsave;
   int     cmps ();		/* comparison routine for ksort */
   int *hsq;

   hsq = ppst->hsq;

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

   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 -1;
   }

   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;
   }

   /* initialize the saved segment structures */
   for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++) {
      vmptr->score = 0;
      memset(vmptr->used, 0, n0 * sizeof(int));
   }

   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 (hsq[aa1[lpos]]>=NMAP) {
       lpos++ ; diagp++;
       while (lpos < n1 && hsq[aa1[lpos]]>=NMAP) {lpos++; diagp++;}
       if (lpos >= n1) break;
       lhval = 0;
     }
     lhval = hsq[aa1[lpos]];
     for (tpos = f_str->harr[lhval].pos, npos = f_str->harr[lhval].next;
	  tpos >= 0; tpos = f_str->link[npos].pos, npos = f_str->link[npos].next) {
       /* tscor gets position of end of current lpos diag run */
       if ((tscor = (dptr = &diagp[-tpos])->stop) >= 0) {
	 tscor++;		 /* move forward one */
	 if ((tscor -= lpos) <= 0) { /* check for size of gap to this hit - */
	                             /* includes implicit -1 mismatch penalty */
	   scor = dptr->score;	     /* current score of this run */
	   if ((tscor += (kfact = f_str->pamh2[lhval])) < 0 &&
	       f_str->lowscor < scor)	/* if updating tscor makes run worse, */
	     savemax (dptr, f_str);     /* save it */

	   if ((tscor += scor) >= kfact) {  /* add to current run if continuing */
	     				    /* is better than restart (kfact) */
	       dptr->score = tscor;
	       dptr->stop = lpos;
	     }
	     else {
	       dptr->score = kfact;	/* starting over is better */
	       dptr->start = (dptr->stop = lpos);
	     }
	 }
	 else {				/* continue current run */
	   dptr->score += f_str->pamh1[aa0[tpos]]; 
	   dptr->stop = lpos;
	 }
       }
       else {				/* no diagonal run yet */
	 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);
     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]
*/

   /* set up pointers for sorting */

   for (nsave = 0, vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++) {
     if (vmptr->score > 0) {
       vmptr->score = m0_spam (aa0, aa1, n1, vmptr, ppst->pam2[0], f_str);
       f_str->vptr[nsave++] = vmptr;
     }
   }

   /* sort them */
   kssort (f_str->vptr, nsave);

   
#ifdef DEBUG
   /*
   for (ib=0; ib<nsave; ib++) {
     fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
  	     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);
       for (im=f_str->vptr[ib]->start; im<=f_str->vptr[ib]->stop; im++)
         fprintf(stderr," %c:%c",ppst->sq[aa0[f_str->noff+im-f_str->vptr[ib]->dp]],
		 ppst->sq[aa1[im]]);
       fputc('\n',stderr);
   }
   fprintf(stderr,"---\n");
   */
   /* now use m_spam to re-evaluate */
   /*
   for (tpos = 0; tpos < n0; tpos++) {
     fprintf(stderr,"%c:%2d ",ppst->sq[aa0[tpos]],aa0[tpos]);
     if (tpos %10 == 9) fputc('\n',stderr);
   }
   fputc('\n',stderr);
   */
#endif   
   
   f_str->aa0ix = 0;
   for (ib=0; ib < nsave; ib++) {
     if ((vmptr=f_str->vptr[ib])->score > 0) {
       vmptr->score = m1_spam (aa0, n0, aa1, n1, vmptr,
			       ppst->pam2[0], ppst->pam_l, f_str);
     }
   }
   /* reset aa0 - modified by m1_spam */
   for (tpos = 0; tpos < n0; tpos++) {
     if (aa0[tpos] >= 32) aa0[tpos] -= 32;
   }

   kssort(f_str->vptr,nsave);

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

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

     return 1;
   }
   else f_str->nsave = nsave;

   
#ifdef DEBUG
   /*
   fprintf(stderr,"n0: %d; n1: %d; noff: %d\n",n0,n1,f_str->noff);
   for (ib=0; ib<nsave; ib++) {
     fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
  	     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);
     for (im=f_str->vptr[ib]->start; im<=f_str->vptr[ib]->stop; im++)
       fprintf(stderr," %c:%c",ppst->sq[aa0[f_str->noff+im-f_str->vptr[ib]->dp]],
  	       ppst->sq[aa1[im]]);
     fputc('\n',stderr);
   }

   fprintf(stderr,"---\n");
   */
#endif   

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

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

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

   return 1;
}

void do_work (const unsigned char *aa0, int n0,
	      const unsigned char *aa1, 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 < 1) {
    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_fastf (f_str->aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff, opt_prob);
#else	/* FASTF */
  do_fastf (f_str->aa0, n0, aa1, n1, ppst, f_str, rst, &hoff, opt_prob);
#endif

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

}

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

  optflag = ppst->param_u.fa.optflag;
  ppst->param_u.fa.optflag = 1;

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

void
savemax (dptr, f_str)
  register struct dstruct *dptr;
  struct f_struct *f_str;
{
   register int dpos;
   register struct savestr *vmptr;
   register int i;

   dpos = (int) (dptr - f_str->diag);

/* check to see if this is the continuation of a run that is already saved */

   if ((vmptr = dptr->dmax) != NULL && vmptr->dp == dpos &&
	 vmptr->start == dptr->start)
   {
      vmptr->stop = dptr->stop;
      if ((i = dptr->score) <= vmptr->score)
	 return;
      vmptr->score = i;
      if (vmptr != f_str->lowmax)
	 return;
   }
   else
   {
      i = f_str->lowmax->score = dptr->score;
      f_str->lowmax->dp = dpos;
      f_str->lowmax->start = dptr->start;
      f_str->lowmax->stop = dptr->stop;
      dptr->dmax = f_str->lowmax;
   }

   for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
      if (vmptr->score < i)
      {
	 i = vmptr->score;
	 f_str->lowmax = vmptr;
      }
   f_str->lowscor = i;
}

/* this version of spam() is designed to work with a collection of
   subfragments, selecting the best amino acid at each position so
   that, from each subfragment, each position is only used once.

   As a result, m_spam needs to know the number of fragments.

   In addition, it now requires a global alignment to the fragment
   and resets the start and stop positions

   */

static int
m1_spam (unsigned char *aa0, int n0,
	 const unsigned char *aa1, int n1,
	 struct savestr *dmax, int **pam2, int pam_l,
	 struct f_struct *f_str)
{
  int     tpos, lpos, im, ii, nm, ci;
  int     tot, ctot, pv;

  struct {
    int     start, stop, score;
  } curv, maxv;
  unsigned char *aa0p;
  const unsigned char *aa1p;

  lpos = dmax->start;                   /* position in library sequence */
   tpos = lpos - dmax->dp + f_str->noff; /* position in query sequence */
   /* force global alignment, reset start*/
   if (tpos < lpos) {
     lpos = dmax->start -= tpos;
     tpos = 0;
   }
   else {
     tpos -= lpos;
     lpos = dmax->start = 0;
   }

   dmax->stop = dmax->start + (f_str->nmoff -2 - tpos);
   if (dmax->stop > n1) dmax->stop = n1;

   /*
   if (dmax->start < 0) {
     tpos = -dmax->start;
     lpos = dmax->start=0;
   }
   else tpos = 0;
   */

   aa1p = &aa1[lpos];
   aa0p = &aa0[tpos];

   nm = f_str->nm0;

   tot = curv.score = maxv.score = 0;
   for (; lpos <= dmax->stop; lpos++,aa0p++,aa1p++) {
     ctot = pam_l;
     ci = -1;
     for (im = 0, ii=0; im < nm; im++,ii+=f_str->nmoff) {
       if (aa0p[ii] < 32 && (pv = pam2[aa0p[ii]][*aa1p]) > ctot) {
	 ctot = pv;
	 ci = ii;
/*  	 fprintf(stderr, "lpos: %d im: %d ii: %d ci: %d ctot: %d pi: %d pv: %d\n", lpos, im, ii, ci, ctot, aa0p[ii], pam2[aa0p[ii]][*aa1p]); */
       }
     }
     tot += ctot;
     if (ci >= 0 && aa0p[ci] < 32) {
#ifdef DEBUG
/*         fprintf(stderr, "used: lpos: %d ci: %d : %c\n", lpos, ci, sq[aa0p[ci]]); */
#endif
       aa0p[ci] +=  32;
       dmax->used[&aa0p[ci] - aa0] = 1;
     }
   }
   return tot;
}

int ma_spam (unsigned char *aa0, int n0, const unsigned char *aa1,
	     struct savestr *dmax, struct pstruct *ppst,
	     struct f_struct *f_str)
{
  int **pam2;
  int     tpos, lpos, im, ii, nm, ci, lp0;
  int     tot, ctot, pv;
  struct {
    int     start, stop, score;
  } curv, maxv;
   const unsigned char *aa1p;
   unsigned char *aa0p, *aa0pt;
   int aa0t_flg;

   pam2 = ppst->pam2[0];
   aa0t_flg = 0;

   lpos = dmax->start;			/* position in library sequence */
   tpos = lpos - dmax->dp + f_str->noff; /* position in query sequence */
   lp0 = lpos = dmax->start;
   aa1p = &aa1[lpos];
   aa0p = &aa0[tpos];			/* real aa0 sequence */

			/* the destination aa0 sequence (without nulls) */
   aa0pt = &f_str->aa0t[f_str->aa0ix];

   curv.start = lpos;
   nm = f_str->nm0;

   /* sometimes, tpos may be > 0, with lpos = 0 - fill with 'X' */
   if (lpos == 0 && tpos > 0)
     for (ii = 0; ii < tpos; ii++) *aa0pt++ = 31;  /* filler character */

   tot = curv.score = maxv.score = 0;
   for (; lpos <= dmax->stop; lpos++) {
     ctot = ppst->pam_l;
     ci = -1;
     for (im = 0, ii=0; im < nm; im++,ii+=f_str->nmoff) {
       if (aa0p[ii] < 32 && (pv = pam2[aa0p[ii]][*aa1p]) > ctot) {
	 ctot = pv;
	 ci = ii;
       }
     }
     tot += ctot;
     if (ci >= 0) {
       if (ci >= n0) {fprintf(stderr," warning - ci off end %d/%d\n",ci,n0);}
       else {
	 *aa0pt++ = aa0p[ci];
	 aa0p[ci] +=  32;
	 aa0t_flg=1;
       }
     }
     aa0p++; aa1p++;
   }

   if (aa0t_flg) {
     dmax->dp -= f_str->aa0ix;		/* shift ->dp for aa0t */
     if ((ci=(int)(aa0pt-f_str->aa0t)) > n0) {
       fprintf(stderr," warning - aapt off %d/%d end\n",ci,n0);
     }
     else 
       *aa0pt++ = 0;			/* skip over NULL */

     aa0pt = &f_str->aa0t[f_str->aa0ix];
     aa1p = &aa1[lp0];

     /*
     for (im = 0; im < f_str->nmoff; im++)
       fprintf(stderr,"%c:%c,",ppst->sq[aa0pt[im]],ppst->sq[aa1p[im]]);
     fprintf(stderr,"- %3d (%3d:%3d)\n",dmax->score,f_str->aa0ix,lp0);
     */

     f_str->aa0ix += f_str->nmoff;	/* update offset into aa0t */
   }
   /*
      fprintf(stderr," ma_spam returning: %d\n",tot);
   */
   return tot;
}

static int
m0_spam (unsigned char *aa0, const unsigned char *aa1, int n1,
	 struct savestr *dmax, int **pam2,
	 struct f_struct *f_str)
{
   int tpos, lpos, lend, im, ii, nm;
   int     tot, ctot, pv;
   struct {
     int     start, stop, score;
   } curv, maxv;
   const unsigned char *aa0p, *aa1p;

   lpos = dmax->start;			/* position in library sequence */
   tpos = lpos - dmax->dp + f_str->noff; /* position in query sequence */
   if (tpos > 0) {
     if (lpos-tpos >= 0) {
       lpos = dmax->start -= tpos;    /* force global alignment, reset start*/
       tpos = 0;
     }
     else {
       tpos -= lpos;
       lpos = dmax->start = 0;
     }
   }

   nm = f_str->nm0;
   lend = dmax->stop;
   if (n1 - (lpos + f_str->nmoff-2) < 0 ) {
     lend = dmax->stop = (lpos - tpos) + f_str->nmoff-2;
     if (lend >= n1) lend = n1-1;
   }

   aa1p = &aa1[lpos];
   aa0p = &aa0[tpos];

   curv.start = lpos;

   tot = curv.score = maxv.score = 0;
   for (; lpos <= lend; lpos++) {
     ctot = -10000;
     for (im = 0, ii=0; im < nm; im++,ii+=f_str->nmoff) {
       if ((pv = pam2[aa0p[ii]][*aa1p]) > ctot) {
	 ctot = pv;
       }
     }
     tot += ctot;
     aa0p++; aa1p++;
   }

   /* reset dmax if necessary */

   return tot;
}

/* sconn links up non-overlapping alignments and calculates the score */

int sconn (struct savestr **v, int n, int cgap, 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[MAXSAV];
   int     lstart, plstop;
   double tatprob;

   /* sarr[] saves each alignment score/position, and provides a link
      back to the previous alignment that maximizes the score */

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

   start = NULL;

   /* for the remaining runs, see if they fit */
   for (i = 0, si = 0; i < n; i++) {

     /* if the score is less than the gap penalty, it never helps */
     if (!opt_prob && (v[i]->score < cgap) ){ continue; }

     lstart = v[i]->start;

     /* 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;
     
     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);
       sarr[si].tat = sarr[si].newtat;
     }

     /* if it fits, then increase the score */
     for (sl = start; sl != NULL; sl = sl->next) {
       plstop = sl->vp->stop;
       /* if end < start or start > end, add score */
       if (plstop < lstart ) {
	 if(!opt_prob) {
	   sarr[si].score = sl->score + v[i]->score;
	   sarr[si].prev = sl;
	   /*
	     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, 2.0);
	   */
	   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;
	   } 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 TAT %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 - resort the scores */
     if (start == NULL) {
       start = &sarr[si];
     } else {
       if(!opt_prob) { /* sort by scores */
	 for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
	   if (sarr[si].score > sj->score) { /* if new score > best score */
	     sarr[si].next = sj;	     /* previous best linked to best */
	     if (so != NULL)		
	       so->next = &sarr[si];	     /* old best points to new best */
	     else
	       start = &sarr[si];
	     break;
	   }
	   so = sj;			     /* old-best saved in so */
	 }
       } else { /* sort by tatprobs */
	 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 {

     if(opt_prob)
       rst->escore = 1.0;
     else
       rst->escore = 2.0;

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

void
kssort (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;
	 }
}

/* 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;
  unsigned char *aa0t;
  const unsigned char *aa1p;

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

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

  /* the alignment portion takes advantage of the information left
     over in f_str after do_fastf 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.
     */

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

   kssort (f_str->vptr, f_str->nsave);
   f_str->aa0ix = 0;
   if (f_str->nsave > f_str->nm0) f_str->nsave = f_str->nm0;
   for (ib=0; ib < f_str->nm0; ib++) {
     if (f_str->vptr[ib]->score > 0) {
       f_str->vptr[ib]->score = 
	 ma_spam (f_str->aa0, n0, aa1p, f_str->vptr[ib], ppst, f_str);
     }
   }

   /* after ma_spam is over, we need to reset aa0 */
   for (ib = 0; ib < n0; ib++) {
     if (f_str->aa0[ib] >= 32) f_str->aa0[ib] -= 32;
   }

   kssort(f_str->vptr,f_str->nsave);

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

  a_res->nres = sconn_a (aa0t,n0, ppst->param_u.fa.cgap, f_str,a_res);
  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 */

int sconn_a (unsigned char *aa0, int n0, int cgap, 
	     struct f_struct *f_str,
	     struct a_res_str *a_res)
{
   int     i, si, cmpp (), n;
   unsigned char *aa0p;
   int sx, dx, doff;

   struct savestr **v;
   struct slink {
     int     score;
     struct savestr *vp;
     struct slink *snext;
     struct slink *aprev;
   } *start, *sl, *sj, *so, sarr[MAXSAV];
   int     lstop, plstart;
   int *res, nres, tres;

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

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

   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 < cgap) continue;

     lstop = v[i]->stop;		/* have right-most lstart */

/*	put the alignment in the group */

     sarr[si].vp = v[i];
     sarr[si].score = v[i]->score;
     sarr[si].snext = NULL;
     sarr[si].aprev = NULL;

/* 	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->snext) { 
       plstart = sl->vp->start;
       if (plstart > lstop ) {
	 sarr[si].score = sl->score + v[i]->score;
	 sarr[si].aprev = sl;
	 break;		/* quit as soon as the alignment has been added */
       }
     }

/* 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->snext) {
	 if (sarr[si].score > sj->score) { /* new score better than old */
	   sarr[si].snext = sj;		/* snext best after new score */
	   if (so != NULL)
	     so->snext = &sarr[si];	/* prev_best->snext points to best */
	   else  start = &sarr[si];	/* start points to best */
	   break;			/* stop looking */
	 }
	 so = sj;		/* previous candidate best */
       }
     si++;				/* increment to snext alignment */
   }

   /* we have the best set of alignments, write them to *res */
   if (start != NULL) {
     res = f_str->res;	/* set a destination for the alignment ops */
     tres = nres = 0;	/* alignment op length = 0 */
     aa0p = aa0;	/* point into query (needed for calcons later) */
     a_res->min1 = start->vp->start;	/* start in library */
     a_res->min0 = 0;			/* start in query */
     for (sj = start; sj != NULL; sj = sj->aprev ) {
       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++) {
	 *aa0p++ = f_str->aa0t[sx++];	/* copy residue into aa0 */
	 tres++;			/* bump alignment counter */
	 res[nres++] = 0;		/* put 0-op in res */
       }
       sj->vp->dp -= doff;
       if (sj->aprev != NULL) {
	 if (sj->aprev->vp->start - sj->vp->stop - 1 > 0 )
	 /* put an insert op into res to get to next aligned block */
	   tres += res[nres++] = (sj->aprev->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;
       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 */
     for (i=0; i<n0; i++) f_str->aa0t[i] = aa0[i];

     return tres;
   }
   else return (0);
}

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

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;
  aln->frame = 0;
  if (frame > 3) aln->llrev = 1;
#else	/* FASTF */
  aln->llfact = aln->qlfact = aln->llmult = 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, *sq, *spa;
  int *rp;
  int mins, smins;

  /* do not allow low complexity */
  sq = pst.sq;
  
#ifndef TFAST
  aa1p = aa1;
  nn1 = n1;
#else
  aa1p = f_str->aa1x;
  nn1 = f_str->n10;
#endif

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

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

  if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1) {
    /* will we show all the start ?*/ 
    smins=0;
    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;
  
  while (i0 < a_res.max0 || i1 < a_res.max1) {
    if (op == 0 && *rp == 0) {
      op = *rp++;

      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 = sq[f_str->aa0t[i0++]];
      *sp1 = sq[aa1p[i1++]];
      n0t++;
      lenc++;
      if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
      sp0++; sp1++; spa++;
    }
    else {
      if (op==0) { op = *rp++;}
      if (op>0) {
	*sp0++ = '-';
	*sp1++ = sq[aa1p[i1++]];
	*spa++ = M_DEL;
	op--;
	len_gap++;
	lenc++;
      }
      else {
	*sp0++ = sq[f_str->aa0t[i0++]];
	*sp1++ = '-';
	*spa++ = M_DEL;
	op++;
	n0t++;
	len_gap++;
	lenc++;
      }
    }
  }

  *spa = '\0';
  *nc = lenc-len_gap;

  /* now we have the middle, get the right end */
  /* ns is amount to be shown */
  /* nd is amount remaining to be shown */
  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 */
    /* there isn't any aa0 to get */
    memset(seqc0+mins+lenc,' ',n0t-a_res.max0);
    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 {
    memset(seqc0+mins+lenc,' ',nd);
    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;
  const unsigned char *aa1p;
  char *sp0, *sp0a, *sp1, *sq, *spa;
  int *rp;
  int mins, smins;

  /* do not allow low complexity */
  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;

  /* 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 ?*/ 
    smins=0;
    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 = 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 ((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++;}

      *sp0a++ = ' ';
      *sp0 = sq[f_str->aa0t[i0++]];
      *sp1 = sq[aa1p[i1++]];
      n0t++;
      lenc++;
      if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
      sp0++; sp1++; spa++;
    }
    else {
      if (op==0) { op = *rp++;}
      if (op>0) {
	*sp0++ = '-';
	*sp0a++ = ' ';
	*sp1++ = sq[aa1p[i1++]];
	*spa++ = M_DEL;
	op--;
	len_gap++;
	lenc++;
      }
      else {
	*sp0++ = sq[f_str->aa0t[i0++]];
	*sp0a++ = ' ';
	*sp1++ = '-';
	*spa++ = M_DEL;
	op++;
	n0t++;
	len_gap++;
	lenc++;
      }
    }
  }

  *sp0a = *spa = '\0';
  *nc = lenc-len_gap;

  /* now we have the middle, get the right end */
  /* ns is amount to be shown */
  /* nd is amount remaining to be shown */
  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 */
    /* there isn't any aa0 to get */
    memset(seqc0+mins+lenc,' ',n0t-a_res.max0);
    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 {
    memset(seqc0+mins+lenc,' ',nd);
    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 aa0shuffle(unsigned char *aa0, int n0, struct f_struct *f_str) {

  int i, j, k;
  unsigned char tmp;

  for (i = f_str->nmoff-1 ; --i ; ) {

    /* j = nrand(i); if (i == j) continue;*/       /* shuffle columns */ 
    j = (f_str->nmoff - 2) - i; if (i <= j) break; /* reverse columns */

    /* swap all i'th column residues for all j'th column residues */
    for(k = 0 ; k < f_str->nm0 ; k++) {
      tmp = aa0[(k * (f_str->nmoff)) + i];
      aa0[(k * (f_str->nmoff)) + i] = aa0[(k * (f_str->nmoff)) + j];
      aa0[(k * (f_str->nmoff)) + j] = tmp;
    }
  }
}

/* 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_op, op_cnt;
  const unsigned char *aa1p;
  char tmp_cnt[20];
  char sp0, sp1, *sq;
  int *rp;
  int mins, smins;
  int 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

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

  i0 = a_res.min0;
  i1 = a_res.min1;
  tmp_cnt[0]='\0';
  
  fnum = f_str->aa0ti[i0] + 1;
  while (i0 < a_res.max0 || i1 < a_res.max1) {
    if (op == 0 && *rp == 0) {
      if (p_op == 0) { 	op_cnt++;}
      else {
	update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt,fnum);
	op_cnt = 1; p_op = 0;
	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 (op>0) {
	if (p_op == 1) { op_cnt++;}
	else {
	  update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt,fnum);
	  op_cnt = 1; p_op = 1; fnum = f_str->aa0ti[i0] + 1;
	}
	op--; lenc++; i1++; len_gap++;
      }
      else {
	if (p_op == 2) { op_cnt++;}
	else {
	  update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt,fnum);
	  op_cnt = 1; p_op = 2; fnum = f_str->aa0ti[i0] + 1;
	}
	op++; lenc++; i0++; len_gap++;
      }
    }
  }
  update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt,fnum);

  return lenc - len_gap;
}

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, n0t;
  int op, lenc, len_gap;
  const unsigned char *aa1p;
  int sp0, sp1;
  int *rp;
  int mins, smins;
  char *sq;

  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

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

  /* now get the middle */
  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;
  
  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 = sq[f_str->aa0t[i0++]];
      sp1 = sq[aa1p[i1++]];
      n0t++;
      lenc++;
      if (toupper(sp0) == toupper(sp1)) aln->nident++;
    }
    else {
      if (op==0) { op = *rp++;}
      if (op>0) {
        i1++;
        op--;
        len_gap++;
        lenc++;
      }
      else {
        i0++;
        op++;
        n0t++;
        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