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

/* $Name: fa_34_26_5 $ - $Id: dropnfa.c,v 1.81 2007/04/26 18:37:19 wrp Exp $ */

/* 18-Sep-2006 - removed global variables for alignment from nw_align
   and bg_align */

/* 18-Oct-2005 - converted to use a_res and aln for alignment coordinates */

/* 14-May-2003 - modified to return alignment start at 0, rather than
   1, for begin:end alignments
*/

/*
  implements the fasta algorithm, see:

  W. R. Pearson, D. J. Lipman (1988) "Improved tools for biological
  sequence comparison" Proc. Natl. Acad. Sci. USA 85:2444-2448

  This version uses Smith-Waterman for final protein alignments

  W. R. Pearson (1996) "Effective protein sequence comparison"
  Methods Enzymol. 266:227-258


  26-April-2001 - -DGAP_OPEN redefines -f, as gap open penalty

  4-Nov-2001 - modify spam() while(1).
*/

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

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

/* this must be consistent with upam.h */
#define MAXHASH 32
#define NMAP MAXHASH+1

/* globals for fasta */
#define MAXWINDOW 64

#ifndef MAXSAV
#define MAXSAV 10
#endif

#ifndef ALLOCN0
static char *verstr="3.5 Sept 2006";
#else
static char *verstr="3.5an0 Sept 2006";
#endif

extern void w_abort(char *, char *);
int shscore(const unsigned char *aa0, int n0, int **pam2);
extern void init_karlin(const unsigned char *aa0, int n0, struct pstruct *ppst,
			double *aa0_f, double **kp);
extern void init_karlin_a(struct pstruct *, double *, double **);
extern int do_karlin(const unsigned char *, int n1, int **,
		     struct pstruct *, double *, double *,
		     double *, double *);
extern void aancpy(char *to, char *from, int count, struct pstruct pst);
char *ckalloc(size_t);

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

#include "dropnfa.h"

#define DROP_INTERN
#include "drop_func.h"

struct swstr { int H, E;};

static int
dmatch (const unsigned char *aa0, int n0,
	const unsigned char *aa1, int n1,
	int hoff, int window, 
	int **pam2, int gdelval, int ggapval,
	struct f_struct *f_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, hv;
   int pamfact;
   int btemp;
   struct f_struct *f_str;
   /* these used to be globals, but do not need to be */
   int ktup;		/* word size examined */
   int fact;		/* factor used to scale ktup match value */
   int kt1;		/* ktup-1 */
   int lkt;		/* last ktup - initiall kt1, but can be increased
			   for hsq >= NMAP */

   int maxn0;		/* used in band alignment */
   int *pwaa;		/* pam[aa0[]] profile */
   int i, j;
   struct swstr *ss;
   int *waa;
   int nsq, ip, *hsq;

   if (ppst->ext_sq_set) {
     nsq = ppst->nsqx; ip = 1;
     hsq = ppst->hsqx;
   }
   else {
     nsq = ppst->nsq; ip = 0;
     hsq = ppst->hsq;
   }

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

#ifndef TFASTA
  if((ppst->zsflag%10) == 6) {
    f_str->kar_p = NULL;
    init_karlin(aa0, n0, ppst, &f_str->aa0_f[0], &f_str->kar_p);
  }
#endif

  btemp = 2 * ppst->param_u.fa.bestoff / 3 +
    n0 / ppst->param_u.fa.bestscale +
    ppst->param_u.fa.bkfact *
    (ppst->param_u.fa.bktup - ppst->param_u.fa.ktup);

  if (ppst->nt_align)
     btemp = (btemp*ppst->pam_h)/5;  /* normalize to standard +5/-4 */

   btemp = min (btemp, ppst->param_u.fa.bestmax);
   if (btemp > 3 * n0) btemp = 3 * shscore(aa0,n0,ppst->pam2[0]) / 5;

   ppst->param_u.fa.cgap = btemp + ppst->param_u.fa.bestoff / 3;

   if (ppst->param_u.fa.optcut_set != 1)
#ifndef TFASTA
      ppst->param_u.fa.optcut = btemp;
#else
      ppst->param_u.fa.optcut = (btemp*3)/2;
#endif

#ifndef OLD_FASTA_GAP
   ppst->param_u.fa.pgap = ppst->gdelval + 2*ppst->ggapval;
#else
   ppst->param_u.fa.pgap = ppst->gdelval + ppst->ggapval;
#endif
   pamfact = ppst->param_u.fa.pamfact;
   ktup = ppst->param_u.fa.ktup;
   fact = ppst->param_u.fa.scfact * ktup;

   if (pamfact == -1) pamfact = 0;
   else if (pamfact == -2) pamfact = 1;

   for (i0 = 1, mhv = -1; i0 <= ppst->nsq; i0++)
      if (hsq[i0] < NMAP && hsq[i0] > mhv) mhv = 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->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 (nsq+1, sizeof (int))) == NULL) {
     fprintf (stderr, " cannot allocate pamh1 array nsq=%d\n",nsq);
     exit (1);
   }

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

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

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

   /* encode the aa0 array */
   phv = hv = 0;
   lkt = kt1;
   /* restart hv, phv calculation */
   for (i0 = 0; i0 < min(lkt,n0); i0++) {
     if (hsq[aa0[i0]] >= NMAP) {hv=phv=0; lkt = i0+ ktup; continue;}
     hv = (hv << f_str->kshft) + hsq[aa0[i0]];
     phv += ppst->pam2[ip][aa0[i0]][aa0[i0]]*ktup;
   }

   for (; i0 < n0; i0++) {
     if (hsq[aa0[i0]] >= NMAP) {
       hv=phv=0;
       /* restart hv, phv calculation */
       for (lkt = i0+kt1; (i0 < lkt || hsq[aa0[i0]]>=NMAP) && i0<n0; i0++) {
	 if (hsq[aa0[i0]] >= NMAP) {
	   hv=phv=0; 
	   lkt = i0+ktup;
	   continue;
	 }
	 hv = (hv << f_str->kshft) + hsq[aa0[i0]];
	 phv += ppst->pam2[ip][aa0[i0]][aa0[i0]]*ktup;
       }
     }
     if (i0 >= n0) break;
     hv = ((hv & f_str->hmask) << f_str->kshft) + hsq[aa0[i0]];
     f_str->link[i0] = f_str->harr[hv];
     f_str->harr[hv] = i0;
     if (pamfact) {
       f_str->pamh2[hv] = (phv += ppst->pam2[ip][aa0[i0]][aa0[i0]] * ktup);
       /* this check should always be true, but just in case */
       if (hsq[aa0[i0-kt1]]<NMAP)
	 phv -= ppst->pam2[ip][aa0[i0 - kt1]][aa0[i0 - kt1]] * ktup;
     }
     else f_str->pamh2[hv] = fact * ktup;
   }

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

   if (pamfact)
      for (i0 = 1; i0 <= nsq; i0++)
	 f_str->pamh1[i0] = ppst->pam2[ip][i0][i0] * ktup;
   else
      for (i0 = 1; i0 <= nsq; i0++)
	 f_str->pamh1[i0] = fact;

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


#ifdef TFASTA
   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

   f_str->bss = (struct bdstr *) calloc((size_t)ppst->param_u.fa.optwid*2+4,
					sizeof(struct bdstr));
   f_str->bss++;
  
   /* allocate space for the scoring arrays */
   maxn0 = n0 + 4;
   if ((ss = (struct swstr *) calloc (maxn0, sizeof (struct swstr)))
	 == NULL) {
     fprintf (stderr, "cannot allocate ss array %3d\n", n0);
     exit (1);
   }
   ss++;
   f_str->ss = ss;

   /* initialize the "variable" pam array */

   if ((waa= (int *)calloc ((size_t)(nsq+1)*n0,sizeof(int))) == NULL) {
     fprintf(stderr,"cannot allocate waa struct %3d\n",nsq*n0);
     exit(1);
   }

   pwaa = waa;
   for (i=0; i<=nsq; i++) {
     for (j=0;j<n0; j++) {
       *pwaa = ppst->pam2[ip][aa0[j]][i];
       pwaa++;
     }
   }
   f_str->waa0 = waa;

   /* initialize the "conventional" pam array used for alignments */

   if ((waa= (int *)calloc ((size_t)(nsq+1)*n0,sizeof(int))) == NULL) {
     fprintf(stderr,"cannot allocate waa struct %3d\n",nsq*n0);
     exit(1);
   }

   pwaa = waa;
   for (i=0; i<=nsq; i++) {
     for (j=0;j<n0; j++) {
       *pwaa = ppst->pam2[0][aa0[j]][i];
       pwaa++;
     }
   }
   f_str->waa1 = waa;

   f_str->max_res = max(3*n0/2,MIN_RES);

   /* now we need alignment storage - get it */
   if ((f_str->res = (int *)calloc((size_t)f_str->max_res,sizeof(int)))==NULL) {
     fprintf(stderr,"cannot allocate alignment results array %d\n",f_str->max_res);
     exit(1);
   }

   *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 TFASTA
  char *pg_str="FASTA";
#else
  char *pg_str="TFASTA";
#endif

   if (!pstr->param_u.fa.optflag)
#ifdef OLD_FASTA_GAP
     sprintf (pstring1, "%s (%s) function [%s matrix, (%d:%d)%s] ktup: %d\n join: %d, gap-pen: %d/%d, width: %3d",
#else
     sprintf (pstring1, "%s (%s) function [%s matrix, (%d:%d)%s] ktup: %d\n join: %d, open/ext: %d/%d, width: %3d",
#endif
	      pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l,
	      (pstr->ext_sq_set) ? "xS":"\0",
	      pstr->param_u.fa.ktup, pstr->param_u.fa.cgap,
	       pstr->gdelval, pstr->ggapval, pstr->param_u.fa.optwid);
   else
#ifdef OLD_FASTA_GAP
      sprintf (pstring1, "%s (%s) function [optimized, %s matrix (%d:%d)%s] ktup: %d\n join: %d, opt: %d, gap-pen: %d/%d, width: %3d",
#else
      sprintf (pstring1, "%s (%s) function [optimized, %s matrix (%d:%d)%s] ktup: %d\n join: %d, opt: %d, open/ext: %d/%d, width: %3d",
#endif
	       pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l,
	       (pstr->ext_sq_set) ? "xS":"\0",
	       pstr->param_u.fa.ktup, pstr->param_u.fa.cgap,
	       pstr->param_u.fa.optcut, pstr->gdelval, pstr->ggapval,
	       pstr->param_u.fa.optwid);
   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) {
#ifdef OLD_FASTA_GAP
     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_optcut: %d\n; pg_cgap: %d\n",
#else
     sprintf (pstring2, "; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)\n\
; pg_open-ext: %d %d\n; pg_ktup: %d\n; pg_optcut: %d\n; pg_cgap: %d\n",
#endif
	      pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l, pstr->gdelval,
              pstr->ggapval,pstr->param_u.fa.ktup,pstr->param_u.fa.optcut,
	      pstr->param_u.fa.cgap);
   }
}

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


  f_str = *f_arg;


  if (f_str != NULL) {
    if (f_str->kar_p!=NULL) free(f_str->kar_p);
    f_str->ss--;
    f_str->bss--;

    free(f_str->res);
    free(f_str->waa1);
    free(f_str->waa0);
    free(f_str->ss);
    free(f_str->bss);
    free(f_str->diag);
    free(f_str->link);
    free(f_str->pamh2); 
    free(f_str->pamh1);
    free(f_str->harr);

    free(f_str);
    *f_arg = NULL;
  }
}

#ifdef ALLOCN0
void savemax (struct dstruct *, int, struct f_struct *);
#else
void savemax (struct dstruct *, struct f_struct *);
#endif

int spam (const unsigned char *, const unsigned char *, struct savestr *,
	 int **, int, int, int);
int sconn(struct savestr **, int nsave, int cgap, int pgap, int noff);
void kpsort(struct savestr **, int);

static int
ALIGN(const unsigned char *, const unsigned char *, int, int, 
      int **, int, int, int *, int *, struct f_struct *);

static int
LOCAL_ALIGN(const unsigned char *, const unsigned char *,
	    int, int, int, int,
	    int **, int, int, int *, int *, int *, int *, int,
	    struct f_struct *);

static int
B_ALIGN(const unsigned char *A, const unsigned char *B, int M,
	int N, int low, int up, int **W, int G, int H, int *S,
	int *nS, int MW, int MX, struct bdstr *bss);

static void
do_fasta (const 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     nd;		/* diagonal array size */
   int     lhval;
   int     kfact;
   register struct dstruct *dptr;
   register int tscor;

#ifndef ALLOCN0
   register struct dstruct *diagp;
#else
   register int dpos;
   int     lposn0;
#endif
   int noff;
   struct dstruct *dpmax;
   register int lpos;
   int     tpos;
   struct savestr *vmptr;
   int     scor, ib, nsave;
   int xdrop, do_extend;
   int ktup, kt1, lkt, *hsq, ip;

  if (ppst->ext_sq_set) {
    ip = 1;
    hsq = ppst->hsqx;
  }
  else {
    ip = 0;
    hsq = ppst->hsq;
  }

  xdrop = -ppst->pam_l;
  /* do extended alignment in spam iff protein or short sequences */
  do_extend = !ppst->nt_align || (n0 < 50) || (n1 < 50);

  ktup = ppst->param_u.fa.ktup;
  kt1 = ktup-1;

  if (n1 < ktup) {
     rst->score[0] = rst->score[1] = rst->score[2] = 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;
     return;
   }

#ifdef ALLOCN0
   nd = n0;
#else
   nd = n0 + n1;
#endif

   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 < &f_str->vmax[MAXSAV]; vmptr++)
      vmptr->score = 0;
   f_str->lowmax = f_str->vmax;
   f_str->lowscor = 0;

   /* start hashing */
   lhval = 0;
   lkt = kt1;
   for (lpos = 0; (lpos < lkt || hsq[aa1[lpos]]>=NMAP) && lpos <n1; lpos++) {
     /* restart lhval calculation */
     if (hsq[aa1[lpos]]>=NMAP) {
       lhval = 0; lkt = lpos + ktup;
       continue;
     }
     lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
   }

   noff = f_str->noff;
#ifndef ALLOCN0
   diagp = &f_str->diag[noff + lkt];
   for (; 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 = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
     for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
       if ((tscor = (dptr = &diagp[-tpos])->stop) >= 0) {
#else
   lposn0 = noff + lpos;
   for (; lpos < n1; lpos++, lposn0++) {
     if (hsq[aa1[lpos]]>=NMAP) {lhval = 0; goto loopl;}
     /*
     if (hsq[aa1[lpos]]>=NMAP) {
       lpos++; lposn0++;
       while (lpos < n1 && hsq[aa1[lpos]]>=NMAP) {lpos++; lposn0++;}
     }
     */
     lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
     for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
       dpos = lposn0 - tpos;
       if ((tscor = (dptr = &f_str->diag[dpos % nd])->stop) >= 0) {
#endif
	 tscor += ktup;
	 if ((tscor -= lpos) <= 0) {
	   scor = dptr->score;
	   if ((tscor += (kfact = f_str->pamh2[lhval])) < 0 && f_str->lowscor < scor)
#ifdef ALLOCN0
	     savemax (dptr, dpos, f_str);
#else
	     savemax (dptr, f_str);
#endif
	     if ((tscor += scor) >= kfact) {
	       dptr->score = tscor;
	       dptr->stop = lpos;
	     }
	     else {
	       dptr->score = kfact;
	       dptr->start = (dptr->stop = lpos) - kt1;
	     }
	 }
	 else {
	   dptr->score += f_str->pamh1[aa0[tpos]];
	   dptr->stop = lpos;
	 }
       }
       else {
	 dptr->score = f_str->pamh2[lhval];
	 dptr->start = (dptr->stop = lpos) - kt1;
       }
     }				/* end tpos */

#ifdef ALLOCN0
      /* reinitialize diag structure */
   loopl:
     if ((dptr = &f_str->diag[lpos % nd])->score > f_str->lowscor)
       savemax (dptr, lpos, f_str);
     dptr->stop = -1;
     dptr->dmax = NULL;
     dptr->score = 0;
#endif
   }				/* end lpos */

#ifdef ALLOCN0
   for (tpos = 0, dpos = noff + n1 - 1; tpos < n0; tpos++, dpos--) {
     if ((dptr = &f_str->diag[dpos % nd])->score > f_str->lowscor)
       savemax (dptr, dpos, f_str);
   }
#else
   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;
#endif

/*
        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 < &f_str->vmax[MAXSAV]; vmptr++)    {
     /*
     fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
	     noff+vmptr->start-vmptr->dp,
	     noff+vmptr->stop-vmptr->dp,
	     vmptr->start,vmptr->stop,
	     vmptr->dp,vmptr->score);

     */
      if (vmptr->score > 0) {
	vmptr->score = spam (aa0, aa1, vmptr, ppst->pam2[ip], xdrop,
			     noff,do_extend);
	f_str->vptr[nsave++] = vmptr;
      }
   }

   if (nsave <= 0) {
     rst->score[0] = rst->score[1] = rst->score[2] = 0;
     return;
   }
       
   /*
   fprintf(stderr,"n0: %d; n1: %d; noff: %d\n",n0,n1,noff);
   for (ib=0; ib<nsave; ib++) {
     fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
	     noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
	     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");
   */

   scor = sconn (f_str->vptr, nsave, ppst->param_u.fa.cgap, 
		 ppst->param_u.fa.pgap, noff);

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

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

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

   if (ppst->param_u.fa.optflag) {
     if (rst->score[0] > ppst->param_u.fa.optcut)
       rst->score[2] = dmatch (aa0, n0, aa1, n1, *hoff=noff - vmptr->dp,
			     ppst->param_u.fa.optwid, ppst->pam2[ip],
			     ppst->gdelval,ppst->ggapval,f_str);
   }
}

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 hoff, n10;

  double lambda, H;
  
  rst->score[0] = rst->score[1] = rst->score[2] = 0;
  rst->escore = 1.0;
  rst->segnum = rst->seglen = 1;

  if (n1 < ppst->param_u.fa.ktup) return;

#ifdef TFASTA  
  n10=aatran(aa1,f_str->aa1x,n1,frame);
  do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff);
#else	/* FASTA */
  do_fasta (aa0, n0, aa1, n1, ppst, f_str, rst, &hoff);
#endif

#ifndef TFASTA
  if((ppst->zsflag%10) == 6 && 
     (do_karlin(aa1, n1, ppst->pam2[0], ppst,f_str->aa0_f, 
		f_str->kar_p, &lambda, &H)>0)) {
    rst->comp = 1.0/lambda;
    rst->H = H;
  }
  else {rst->comp = rst->H = -1.0;}
#else
  rst->comp = rst->H = -1.0;
#endif
}

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 TFASTA  
  n10=aatran(aa1,f_str->aa1x,n1,frame);
  do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff);
#else	/* FASTA */
  do_fasta(aa0,n0,aa1,n1,ppst,f_str,rst, &hoff);
#endif
  ppst->param_u.fa.optflag = optflag;
}

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

#else
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);

#endif

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

int spam (const unsigned char *aa0, const unsigned char *aa1,
	  struct savestr *dmax, int **pam2, int xdrop,
	  int noff, int do_extend)
{
   register int lpos, tot;
   register const unsigned char *aa0p, *aa1p;

   int drop_thresh;

   struct {
     int     start, stop, score;
   } curv, maxv;

   aa1p = &aa1[lpos= dmax->start];	/* get the start of lib seq */
   aa0p = &aa0[lpos - dmax->dp + noff];	/* start of query */
#ifdef DEBUG
   /* also add check in calling routine */
   if (aa0p < aa0) { return -99; }
#endif
   curv.start = lpos;			/* start index in lib seq */

   tot = curv.score = maxv.score = 0;

   for (; lpos <= dmax->stop; lpos++) {
     tot += pam2[*aa0p++][*aa1p++];
     if (tot > curv.score) {		/* update current score */
       curv.stop = lpos;
       curv.score = tot;
      }
      else if (tot < 0) {
	if (curv.score > maxv.score) {	/* save score, start, stop */
	  maxv.start = curv.start;
	  maxv.stop = curv.stop;
	  maxv.score = curv.score;
	}
	tot = curv.score = 0;		/* reset running score */
	curv.start = lpos+1;		/* reset start */
	if(lpos >= dmax->stop) break;	/* if the zero is beyond stop, quit */
      }
   }

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

#ifndef NOSPAM_EXT

   /* now check to see if the score gets better by extending */
   if (do_extend && maxv.score > xdrop) {

     if (maxv.stop == dmax->stop) {
       tot = maxv.score;
       drop_thresh = maxv.score - xdrop;
       aa1p = &aa1[lpos= dmax->stop];
       aa0p = &aa0[lpos - dmax->dp + noff];
       while (tot > drop_thresh ) {
	 ++lpos;
	 tot += pam2[*(++aa0p)][*(++aa1p)];
	 if (tot > maxv.score) {
	   maxv.start = lpos;
	   maxv.score = tot;
	   drop_thresh = tot - xdrop;
	 }
       }
     }

     /* scan backwards now */

     if (maxv.start == dmax->start) {
       tot = maxv.score;
       drop_thresh = maxv.score - xdrop;
       aa1p = &aa1[lpos= dmax->start];
       aa0p = &aa0[lpos - dmax->dp + noff];
       while (tot > drop_thresh) {
	 --lpos;
	 tot += pam2[*(--aa0p)][*(--aa1p)];
	 if (tot > maxv.score) {
	   maxv.start = lpos;
	   maxv.score = tot;
	   drop_thresh = tot - xdrop;
	 }
       }
     }
   }
#endif

/*	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, int cgap, int pgap, int noff)
{
   int     i, si;
   struct slink
   {
      int     score;
      struct savestr *vp;
      struct slink *next;
   }      *start, *sl, *sj, *so, sarr[MAXSAV];
   int     lstart, tstart, plstop, ptstop;

/*	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 (v[i]->score < cgap)
	 continue;
      lstart = v[i]->start;
      tstart = lstart - v[i]->dp + noff;

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

/* 	if it fits, then increase the score */
      for (sl = start; sl != NULL; sl = sl->next)
      {
	 plstop = sl->vp->stop;
	 ptstop = plstop - sl->vp->dp + noff;
	 if (plstop < lstart && ptstop < tstart)
	 {
	    sarr[si].score = sl->score + v[i]->score + pgap;
	    break;
	 }
      }

/*	now recalculate where the score fits */
      if (start == NULL)
	 start = &sarr[si];
      else
	 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;
	 }
      si++;
   }

   if (start != NULL)
      return (start->score);
   else
      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;
	 }
}

static int dmatch (const unsigned char *aa0, int n0,
		   const unsigned char *aa1, int n1,
		   int hoff, int window, 
		   int **pam2, int gdelval, int ggapval,
		   struct f_struct *f_str)
{
   int low, up;

   window = min (n1, window);
   /* hoff is the offset found from aa1 to seq 2 by hmatch */

   low = -window/2-hoff;
   up = low+window;

   return FLOCAL_ALIGN(aa0-1,aa1-1,n0,n1, low, up,
		      pam2,
#ifdef OLD_FASTA_GAP
		       -(gdelval-ggapval),
#else
		       -gdelval,
#endif
		       -ggapval,window,f_str);
 }


/* A PACKAGE FOR LOCALLY ALIGNING TWO SEQUENCES WITHIN A BAND:

   To invoke, call LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,MW).
   The parameters are explained as follows:
	A, B : two sequences to be aligned
	M : the length of sequence A
	N : the length of sequence B
	L : lower bound of the band
	U : upper bound of the band
	W : scoring table for matches and mismatches
	G : gap-opening penalty
	H : gap-extension penalty
	MW  : maximum window size
*/

#include <stdio.h>

#define MININT -9999999

static int
FLOCAL_ALIGN(const unsigned char *A, const unsigned char *B,
	     int M, int N, int low, int up,
	     int **W, int G,int H, int MW,
	     struct f_struct *f_str)
{
  int band;
  register struct bdstr *bssp;
  int i, j, si, ei;
  int c, d, e, m;
  int leftd, rightd;
  int best_score;
  int *wa, curd;
  int ib;
  
  bssp = f_str->bss;

  m = G+H;
  low = max(-M, low);
  up = min(N, up);
  
  if (N <= 0) return 0;

  if (M <= 0) return 0;

  band = up-low+1;
  if (band < 1) {
    fprintf(stderr,"low > up is unacceptable!: M: %d N: %d l/u: %d/%d\n",
	    M, N, low, up);
    return 0;
  }

  if (low > 0) leftd = 1;
  else if (up < 0) leftd = band;
  else leftd = 1-low;
  rightd = band;
  si = max(0,-up);	/* start index -1 */
  ei = min(M,N-low);	/* end index */
  bssp[leftd].CC = 0;
  for (j = leftd+1; j <= rightd; j++) {
    bssp[j].CC = 0;
    bssp[j].DD = -G;
  }

  bssp[rightd+1].CC = MININT;
  bssp[rightd+1].DD = MININT;

  best_score = 0;
  bssp[leftd-1].CC = MININT;
  bssp[leftd].DD = -G;

  for (i = si+1; i <= ei; i++) {
    if (i > N-up) rightd--;
    if (leftd > 1) leftd--;
    wa = W[A[i]];
    if ((c = bssp[leftd+1].CC-m) > (d = bssp[leftd+1].DD-H)) d = c;
    if ((ib = leftd+low-1+i ) > 0) c = bssp[leftd].CC+wa[B[ib]];

    if (d > c) c = d;
    if (c < 0) c = 0;
    e = c-G;
    bssp[leftd].DD = d;
    bssp[leftd].CC = c;
    if (c > best_score) best_score = c;

    for (curd=leftd+1; curd <= rightd; curd++) {
      if ((c = c-m) > (e = e-H)) e = c;
      if ((c = bssp[curd+1].CC-m) > (d = bssp[curd+1].DD-H)) d = c;
      c = bssp[curd].CC + wa[B[curd+low-1+i]];
      if (e > c) c = e;
      if (d > c) c = d;
      if (c < 0) c = 0;
      bssp[curd].CC = c;
      bssp[curd].DD = d;
      if (c > best_score) best_score = c;
    }
  }

  return best_score;
}

/* ckalloc - allocate space; check for success */
char *ckalloc(size_t amount)
{
  char *p;

  if ((p = malloc( (unsigned) amount)) == NULL)
    w_abort("Ran out of memory.","");
  return(p);
}

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

int  sw_walign (const unsigned char *aa0, int n0,
		const unsigned char *aa1, int n1,
		struct pstruct *ppst, 
		struct f_struct *f_str,
		struct a_res_str *a_res
		)
{
   register const unsigned char *aa0p, *aa1p;
   register int *pwaa;
   register int i, j;
   register struct swstr *ssj;
   struct swstr *ss;
   int *waa;
   int e, f, h, p;
   int     q, r, m;
   int     score;
   int cost, I, J, K, L;

   ss = f_str->ss;
   waa = f_str->waa1;

#ifdef OLD_FASTA_GAP
   q = -(ppst->gdelval - ppst->ggapval);
#else
   q = -ppst->gdelval;
#endif
   r = -ppst->ggapval;
   m = q + r;

   /* initialize 0th row */
   for (ssj=ss; ssj<ss+n0; ssj++) {
     ssj->H = 0;
     ssj->E = -q;
   }

   score = I = J = 0;
   aa1p = aa1;
   i = 0;
   while (*aa1p) {
     h = p = 0;
     f = -q;
     pwaa = waa + (*aa1p++ * n0);
     for (ssj = ss, aa0p = aa0; ssj < ss+n0; ssj++) {
       if ((h =   h     - m) > (f =   f     - r)) f = h;
       if ((h = ssj->H - m) > (e = ssj->E - r)) e = h;
       h = p + *pwaa++;
       if (h < 0 ) h = 0;
       if (h < f ) h = f;
       if (h < e ) h = e;
       p = ssj->H;
       ssj->H = h;
       ssj->E = e;
       if (h > score) {
	 score = h;
	 I = i;
	 J = (int)(ssj-ss);
       }
     }
     i++;
   }				/* done with forward pass */
   if (score <= 0) return 0;  

  /* to get the start point, go backwards */
  
  cost = K = L = 0;
  for (ssj=ss+J; ssj>=ss; ssj--) ssj->H= ssj->E= -1;
  
  for (i=I; i>=0; i--) {
    h = f = -1;
    p = (i == I) ? 0 : -1;
    ssj = ss+J;  /* bug in compiler */
    for (aa0p = &aa0[J]; ssj>=ss; ssj--,aa0p--) {
      f = max (f,h-q)-r;
      ssj->E=max(ssj->E,ssj->H-q)-r;
      h = max(max(ssj->E,f),p+ppst->pam2[0][*aa0p][aa1[i]]);
      p = ssj->H;
      ssj->H=h;
      if (h > cost) {
	cost = h;
	K = i;
	L = (int)(ssj-ss);
	if (cost >= score) goto found;
      }
    }
  }
  
found:	

  /*  printf(" %d: L: %3d-%3d/%3d; K: %3d-%3d/%3d\n",score,L,J,n0,K,I,n1); */

/* in the f_str version, the *res array is already allocated at 4*n0/3 */

  a_res->max0 = J+1; a_res->min0 = L; a_res->max1 = I+1; a_res->min1 = K;
  
  /* the seq array arguments in this call have been reversed to allow
     assymetric scoring matrices - this affects the score decoding, 
     and allocation of the score row matrix */
  ALIGN(&aa0[L-1],&aa1[K-1],J-L+1,I-K+1,ppst->pam2[0],q,r,a_res->res,&a_res->nres,f_str);

  /*  DISPLAY(&aa1[K-1],&aa0[L-1],I-K+1,J-L+1,res,L,K,ppst->sq); */

   return score;
}

static int CHECK_SCORE(const unsigned char *A, const unsigned char *B, 
		       int M, int N,
		       int *S, int **W, int G, int H, int *nres);

#define gap(k)  ((k) <= 0 ? 0 : g+h*(k))	/* k-symbol indel cost */

/* static int *sapp; */				/* Current script append ptr */
/* static int  last; */				/* Last script op appended */

						/* Append "Delete k" op */
#define DEL(k)				\
{ if (*last < 0)				\
    *last = (*sapp)[-1] -= (k);		\
  else	{				\
    *last = (*sapp)[0] = -(k);		\
    (*sapp)++;				\
  }					\
}
						/* Append "Insert k" op */
#define INS(k)				\
{ if (*last > 0)			\
     *last = (*sapp)[-1] += (k); 	\
  else	{				\
    *last = (*sapp)[0] = (k);		\
    (*sapp)++;				\
  }					\
}

#define REP { *last = (*sapp)[0] = 0; (*sapp)++;} /* Append "Replace" op */

/* align(A,B,M,N,tb,te) returns the cost of an optimum conversion between
   A[1..M] and B[1..N] that begins(ends) with a delete if tb(te) is zero
   and appends such a conversion to the current script.                   */

static int
nw_align(const unsigned char *A, const unsigned char *B,
	 int M, int N,
	 int tb, int te, int **w, int g, int h, 
	 struct f_struct *f_str,
	 int **sapp, int *last)
{
  int   midi, midj, type;	/* Midpoint, type, and cost */
  int midc;

  register int   i, j;
  register int c, e, d, s;
  int m, t, *wa;
  struct swstr *f_ss, *r_ss;

  m = g + h;

  f_ss = f_str->f_ss;
  r_ss = f_str->r_ss;

/* Boundary cases: M <= 1 or N == 0 */

  if (N <= 0) {
    if (M > 0) {DEL(M)}
     return -gap(M);
  }

  if (M <= 1) {
    if (M <= 0) {
      INS(N);
      return -gap(N);
    }
    if (tb < te) tb = te;
    midc = (tb-h) - gap(N);
    midj = 0;
    wa = w[A[1]];	/* in the original version of this code, A[]
			   is the second sequence */
    for (j = 1; j <= N; j++) {
      c = -gap(j-1) + wa[B[j]] - gap(N-j);
      if (c > midc) {
	midc = c;
	midj = j;
      }
    }
    if (midj == 0) { DEL(1) INS(N) }
    else {
      if (midj > 1) { INS(midj-1) }
      REP
      if (midj < N) { INS(N-midj) }
    }
    return midc;
  }

/* Divide: Find optimum midpoint (midi,midj) of cost midc */

  midi = M/2;		/* Forward phase:                          */
  f_ss[0].H = 0;		/*   Compute H(M/2,k) & E(M/2,k) for all k */
  t = -g;
  for (j = 1; j <= N; j++) {
    f_ss[j].H = t = t-h;
    f_ss[j].E = t-g;
  }
  t = tb;
  for (i = 1; i <= midi; i++) {
    s = f_ss[0].H;
    f_ss[0].H = c = t = t-h;
    e = t-g;
    wa = w[A[i]];
    for (j = 1; j <= N; j++) {
      if ((c =   c   - m) > (e =   e   - h)) e = c;
      if ((c = f_ss[j].H - m) > (d = f_ss[j].E - h)) d = c;
      c = s + wa[B[j]];
      if (e > c) c = e;
      if (d > c) c = d;
      s = f_ss[j].H;
      f_ss[j].H = c;
      f_ss[j].E = d;
    }
  }
  f_ss[0].E = f_ss[0].H;

  r_ss[N].H = 0;		/* Reverse phase:                          */
  t = -g;			/*   Compute R(M/2,k) & S(M/2,k) for all k */
  for (j = N-1; j >= 0; j--)
    { r_ss[j].H = t = t-h;
      r_ss[j].E = t-g;
    }
  t = te;
  for (i = M-1; i >= midi; i--)
    { s = r_ss[N].H;
      r_ss[N].H = c = t = t-h;
      e = t-g;
      wa = w[A[i+1]];
      for (j = N-1; j >= 0; j--)
        { if ((c =   c   - m) > (e =   e   - h)) e = c;
          if ((c = r_ss[j].H - m) > (d = r_ss[j].E - h)) d = c;
          c = s + wa[B[j+1]];
          if (e > c) c = e;
          if (d > c) c = d;
          s = r_ss[j].H;
          r_ss[j].H = c;
          r_ss[j].E = d;
        }
    }
  r_ss[N].E = r_ss[N].H;

  midc = f_ss[0].H+r_ss[0].H;		/* Find optimal midpoint */
  midj = 0;
  type = 1;
  for (j = 0; j <= N; j++)
    if ((c = f_ss[j].H + r_ss[j].H) >= midc)
      if (c > midc || (f_ss[j].H != f_ss[j].E && r_ss[j].H == r_ss[j].E))
        { midc = c;
          midj = j;
        }
  for (j = N; j >= 0; j--)
    if ((c = f_ss[j].E + r_ss[j].E + g) > midc)
      { midc = c;
        midj = j;
        type = 2;
      }


/* Conquer: recursively around midpoint */

  if (type == 1) {
    nw_align(A,B,midi,midj,tb,-g,w,g,h,f_str, sapp, last);
    nw_align(A+midi,B+midj,M-midi,N-midj,-g,te,w,g,h,f_str,sapp, last);
  }
  else {
    nw_align(A,B,midi-1,midj,tb,0,w,g,h,f_str, sapp, last);
    DEL(2);
    nw_align(A+midi+1,B+midj,M-midi-1,N-midj,0,te,w,g,h,f_str, sapp, last);
  }
  return midc;
}

/* Interface and top level of comparator */

static int
ALIGN(const unsigned char *A, const unsigned char *B,
      int M, int N,
      int **W, int G, int H, int *S, int *nS,
      struct f_struct *f_str)
{ 
  int c, ck;
  struct swstr *f_ss, *r_ss;
  int *sapp, last;

  sapp = S;
  last = 0;

   if ((f_ss = (struct swstr *) calloc (N+2, sizeof (struct swstr)))
       == NULL) {
     fprintf (stderr, "cannot allocate f_ss array %3d\n", N+2);
     exit (1);
   }
   f_ss++;
   f_str->f_ss = f_ss;

   if ((r_ss = (struct swstr *) calloc (N+2, sizeof (struct swstr)))
       == NULL) {
     fprintf (stderr, "cannot allocate r_ss array %3d\n", N+2);
     exit (1);
   }
   r_ss++;
   f_str->r_ss = r_ss;

  c = nw_align(A,B,M,N,-G,-G,W,G,H,f_str,&sapp, &last);	/* OK, do it */

  ck = CHECK_SCORE(A,B,M,N,S,W,G,H,nS);
  if (c != ck) fprintf(stderr,"Check_score error %d != %d\n",c,ck);

  f_ss--; r_ss--;
  free(r_ss); free(f_ss);

  return c;
}

/* Alignment display routine */

static char ALINE[51], BLINE[51], CLINE[51];

void DISPLAY(unsigned char *A, unsigned char *B, int M, int N,
	    int *S, int AP, int BP, char *sq)
{ register char *a, *b, *c;
  register int   i,  j, op;
           int   lines, ap, bp;

  i = j = op = lines = 0;
  ap = AP;
  bp = BP;
  a = ALINE;
  b = BLINE;
  c = CLINE;
  while (i < M || j < N)
    { if (op == 0 && *S == 0)
        { op = *S++;
          *a = sq[A[++i]];
          *b = sq[B[++j]];
          *c++ = (*a++ == *b++) ? '|' : ' ';
        }
      else
        { if (op == 0)
            op = *S++;
          if (op > 0)
            { *a++ = ' ';
              *b++ = sq[B[++j]];
              op--;
            }
          else
            { *a++ = sq[A[++i]];
              *b++ = ' ';
              op++;
            }
          *c++ = '-';
        }
      if (a >= ALINE+50 || (i >= M && j >= N))
        { *a = *b = *c = '\0';
          printf("\n%5d ",50*lines++);
          for (b = ALINE+10; b <= a; b += 10)
            printf("    .    :");
          if (b <= a+5)
            printf("    .");
          printf("\n%5d %s\n      %s\n%5d %s\n",ap,ALINE,CLINE,bp,BLINE);
	  ap = AP + i;
	  bp = BP + j;
          a = ALINE;
          b = BLINE;
          c = CLINE;
        }
    }
}

/* CHECK_SCORE - return the score of the alignment stored in S */

static int CHECK_SCORE(const unsigned char *A, const unsigned char *B,
		       int M, int N, int *S, int **w, int g, int h,
		       int *nres)
{ 
  register int   i,  j, op, nc;
  int score;

  score = i = j = op = nc = 0;
  while (i < M || j < N) {
	op = *S++;
	if (op == 0) {
	  score = w[A[++i]][B[++j]] + score;
	  nc++;
  /*  	  fprintf(stderr,"=%4d %4d %4d %4d\n",i,j,w[A[i]][B[i]],score); */
	}
	else if (op > 0) {
	  score = score - (g+op*h);
  /*  	  fprintf(stderr,">%4d %4d %4d %4d\n",i,j,-(g+op*h),score); */
	  j += op;
	  nc += op;
	} else {
	  score = score - (g-op*h);
  /*  	  fprintf(stderr,"<%4d %4d %4d %4d\n",i,j,-(g-op*h),score); */
	  i -= op;
	  nc -= op;
	}
  }
  *nres = nc;
  return score;
}


static int
BCHECK_SCORE(const unsigned char *A, const unsigned char *B,
		       int M, int N, int *S, int **w, int g, int h,
		       int *nres)
{ 
  register int   i,  j, op, nc;
  int *Ssave;
  int score;

  score = i = j = op = nc = 0;
  Ssave = S;
  while (i < M || j < N) {
	op = *S++;
	if (op == 0) {
	  score = w[A[++i]][B[++j]] + score;
	  nc++;
/*  	  fprintf(stderr,"op0 %4d %4d %4d %4d\n",i,j,w[A[i]][B[i]],score); */
	}
	else if (op > 0) {
	  score = score - (g+op*h);
/*  	  fprintf(stderr,"op> %4d %4d %4d %4d %4d\n",i,j,op,-(g+op*h),score); */
	  j += op;
	  nc += op;
	} else {
	  score = score - (g-op*h);
/*  	  fprintf(stderr,"op< %4d %4d %4d %4d %4d\n",i,j,op,-(g-op*h),score); */
	  i -= op;
	  nc -= op;
	}
  }
  *nres = nc;
  return score;
}


/* A PACKAGE FOR LOCALLY ALIGNING TWO SEQUENCES WITHIN A BAND:

   To invoke, call LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,S,dflag,&SI,&SJ,&EI,&EJ,MW).
   The parameters are explained as follows:
	A, B : two sequences to be aligned
	M : the length of sequence A
	N : the length of sequence B
	L : lower bound of the band
	U : upper bound of the band
	W : scoring table for matches and mismatches
	G : gap-opening penalty
	H : gap-extension penalty
	dflag : 0 - no display or backward pass
	*SI : starting position of sequence A in the optimal local alignment
	*SJ : starting position of sequence B in the optimal local alignment
	*EI : ending position of sequence A in the optimal local alignment
	*EJ : ending position of sequence B in the optimal local alignment
	MW  : maximum window size
*/

int bd_walign (const unsigned char *aa0, int n0,
		 const unsigned char *aa1, int n1,
		 struct pstruct *ppst, 
		 struct f_struct *f_str, int hoff,
		 struct a_res_str *a_res)
{
   int low, up, score;
   int min0, min1, max0, max1;
   int window;

   window = min (n1, ppst->param_u.fa.optwid);
   /* hoff is the offset found from aa1 to seq 2 by hmatch */

   low = -window/2-hoff;
   up = low+window;

   score=LOCAL_ALIGN(aa0-1,aa1-1,n0,n1, low, up,
		    ppst->pam2[0],
#ifdef OLD_FASTA_GAP
		     -(ppst->gdelval-ppst->ggapval),
#else
		     -ppst->gdelval,
#endif
		     -ppst->ggapval,
		    &min0,&min1,&max0,&max1,ppst->param_u.fa.optwid,f_str);
  
   if (score <=0) {
     fprintf(stderr,"n0/n1: %d/%d hoff: %d window: %d\n",
	     n0, n1, hoff, window);
     return 0;
   }
  
/*
  fprintf(stderr," ALIGN: start0: %d start1: %d stop0: %d stop1: %d, bot: %d top: %d, win: %d MX %d\n",
  min0-1,min1-1,max0-min0+1,max1-min1+1,low-(min1-min0),up-(min1-min0),
  ppst->param_u.fa.optwid,n0);
*/

   a_res->min0 = min0-1; a_res->min1 = min1-1;
   a_res->max0 = max0; a_res->max1 = max1; 

   B_ALIGN(aa0-1+min0-1,aa1-1+min1-1,max0-min0+1,max1-min1+1,
	   low-(min1-min0),up-(min1-min0),
	   ppst->pam2[0],
#ifdef OLD_FASTA_GAP
	   -(ppst->gdelval-ppst->ggapval),
#else
	   -ppst->gdelval,
#endif
	   -ppst->ggapval,
	   a_res->res,&a_res->nres,ppst->param_u.fa.optwid,n0,f_str->bss);

   return score;
}

static int
LOCAL_ALIGN(const unsigned char *A, const unsigned char *B,
	    int M, int N,
	    int low, int up, int **W, int G,int H,
	    int *psi, int *psj, int *pei, int *pej, int MW,
	    struct f_struct *f_str)
{ 
  int band;
  register struct bdstr *bssp;
  int i, j, si, ei;
  int c, d, e, t, m;
  int leftd, rightd;
  int best_score, starti, startj, endi, endj;
  int *wa, curd;
  int ib;
  char flag;
  
  bssp = f_str->bss;

  m = G+H;
  low = max(-M, low);
  up = min(N, up);
  
  if (N <= 0) { 
    *psi = *psj = *pei = *pej;
    return 0;
  }
  if (M <= 0) {
    *psi = *psj = *pei = *pej;
    return 0;
  }
  band = up-low+1;
  if (band < 1) {
    fprintf(stderr,"low > up is unacceptable!: M: %d N: %d l/u: %d/%d\n",
	    M, N, low, up);
    return -1;
  }

  j = (MW + 2 + 2) * sizeof(struct bdstr);

  /* already done by init_work(); 
  if (f_str->bss==NULL) f_str->bss = (struct bdstr *) ckalloc(j);
  */

  if (low > 0) leftd = 1;
  else if (up < 0) leftd = band;
  else leftd = 1-low;
  rightd = band;
  si = max(0,-up);
  ei = min(M,N-low);
  bssp[leftd].CC = 0;
  for (j = leftd+1; j <= rightd; j++) {
    bssp[j].CC = 0;
    bssp[j].DD = -G;
  }
  bssp[rightd+1].CC = MININT;
  bssp[rightd+1].DD = MININT;
  best_score = 0;
  endi = si;
  endj = si+low;
  bssp[leftd-1].CC = MININT;
  bssp[leftd].DD = -G;
  for (i = si+1; i <= ei; i++) {
    if (i > N-up) rightd--;
    if (leftd > 1) leftd--;
    wa = W[A[i]];
    if ((c = bssp[leftd+1].CC-m) > (d = bssp[leftd+1].DD-H)) d = c;
    if ((ib = leftd+low-1+i ) > 0) c = bssp[leftd].CC+wa[B[ib]];
/*
    if (ib > N) fprintf(stderr,"B[%d] out of range %d\n",ib,N);
*/
    if (d > c) c = d;
    if (c < 0) c = 0;
    e = c-G;
    bssp[leftd].DD = d;
    bssp[leftd].CC = c;
    if (c > best_score) {
      best_score = c;
      endi = i;
      endj = ib;
    }
    for (curd=leftd+1; curd <= rightd; curd++) {
      if ((c = c-m) > (e = e-H)) e = c;
      if ((c = bssp[curd+1].CC-m) > (d = bssp[curd+1].DD-H)) d = c;
/*
      if ((ib=curd+low-1+i) <= 0 || ib > N)
	fprintf(stderr,"B[%d]:%d\n",ib,B[ib]);
*/
      c = bssp[curd].CC + wa[B[curd+low-1+i]];
      if (e > c) c = e;
      if (d > c) c = d;
      if (c < 0) c = 0;
      bssp[curd].CC = c;
      bssp[curd].DD = d;
      if (c > best_score) {
	best_score = c;
	endi = i;
	endj = curd+low-1+i;
      }
    }
  }
  
  leftd = max(1,-endi-low+1);
  rightd = band-(up-(endj-endi));
  bssp[rightd].CC = 0;
  t = -G;
  for (j = rightd-1; j >= leftd; j--) {
    bssp[j].CC = t = t-H;
    bssp[j].DD = t-G;
  }
  for (j = rightd+1; j <= band; ++j) bssp[j].CC = MININT;
  bssp[leftd-1].CC = bssp[leftd-1].DD = MININT;
  bssp[rightd].DD = -G;
  flag = 0;
  for (i = endi; i >= 1; i--) {
    if (i+low <= 0) leftd++;
    if (rightd < band) rightd++;
    wa = W[A[i]];
    if ((c = bssp[rightd-1].CC-m) > (d = bssp[rightd-1].DD-H)) d = c;
    if ((ib = rightd+low-1+i) <= N) c = bssp[rightd].CC+wa[B[ib]];

/*
    if (ib <= 0) fprintf(stderr,"rB[%d] <1\n",ib);
*/
    if (d > c) c = d;
    e = c-G;
    bssp[rightd].DD = d;
    bssp[rightd].CC = c;
    if (c == best_score) {
      starti = i;
      startj = ib;
      flag = 1;
      break;
    }
    for (curd=rightd-1; curd >= leftd; curd--) {
      if ((c = c-m) > (e = e-H)) e = c;
      if ((c = bssp[curd-1].CC-m) > (d = bssp[curd-1].DD-H)) d = c;

/*
      if ((ib=curd+low-1+i) <= 0 || ib > N)
	fprintf(stderr,"i: %d, B[%d]:%d\n",i,ib,B[ib]);
*/
      c = bssp[curd].CC + wa[B[curd+low-1+i]];
      if (e > c) c = e;
      if (d > c) c = d;
      bssp[curd].CC = c;
      bssp[curd].DD = d;
      if (c == best_score) {
	starti = i;
	startj = curd+low-1+i;
	flag = 1;
	break;
      }
    }
    if (flag == 1) break;
  }
  
  if (starti < 0 || starti > M || startj < 0 || startj > N) {
    printf("starti=%d, startj=%d\n",starti,startj);
    *psi = *psj = *pei = *pej;
    exit(1);
  }
  *psi = starti;
  *psj = startj;
  *pei = endi;
  *pej = endj;
  return best_score;
}

/* A PACKAGE FOR GLOBALLY ALIGNING TWO SEQUENCES WITHIN A BAND:

   To invoke, call B_ALIGN(A,B,M,N,L,U,W,G,H,S,MW,MX).
   The parameters are explained as follows:
	A, B : two sequences to be aligned
	M : the length of sequence A
	N : the length of sequence B
	L : lower bound of the band
	U : upper bound of the band
	W : scoring table for matches and mismatches
	G : gap-opening penalty
	H : gap-extension penalty
	S : script for DISPLAY routine
	MW : maximum window size
	MX : maximum length sequence M to be aligned
*/

static int IP;
static int *MP[3];		/* save crossing points */
static int *FP;			/* forward dividing points */
static char *MT[3];		/* 0: rep, 1: del, 2: ins */
static char *FT;

/* bg_align(A,B,M,N,up,low,tb,te) returns the cost of an optimum conversion between
  A[1..M] and B[1..N] and appends such a conversion to the current script.
  tb(te)= 1  no gap-open penalty if the conversion begins(ends) with a delete.
  tb(te)= 2  no gap-open penalty if the conversion begins(ends) with an insert.
*/
static int
bg_align(const unsigned char *A, const unsigned char *B, 
	 int M, int N,
	 int low, int up, int tb, int te,
	 int **w, int g, int h,
	 struct bdstr *bss, int **sapp, int *last)
{
  int rmid, k, l, r, v, kt;
  int t1, t2, t3;

  {
  int band, midd;
  int leftd, rightd;	/* for CC, DD, CP and DP */
  register int curd;	/* current index for CC, DD CP and DP */
  register int i, j;
  register int c, d, e;
  int t, fr, *wa, ib, m;

  /* Boundary cases: M <= 0 , N <= 0, or up-low <= 0 */
  if (N <= 0) { 
    if (M > 0) { DEL(M) }
    return 0;
  }
  if (M <= 0) {
    INS(N)
    return 0;
  }
  if ((band = up-low+1) <= 1) {
    for (i = 1; i <= M; i++) { REP }
    return 0;
  }

  /* Divide: Find all crossing points */

  /* Initialization */
  m = g + h;

  midd = band/2 + 1;
  rmid = low + midd - 1;
  leftd = 1-low;
  rightd = up-low+1;
  if (leftd < midd) {
    fr = -1;
    for (j = 0; j < midd; j++) 
      bss[j].CP = bss[j].DP = -1;
    for (j = midd; j <= rightd; j++) {
      bss[j].CP = bss[j].DP = 0;
    }
    MP[0][0] = -1;
    MP[1][0] = -1;
    MP[2][0] = -1;
    MT[0][0] = MT[1][0] = MT[2][0] = 0;
  } else if (leftd > midd) {
    fr = leftd-midd;
    for (j = 0; j <= midd; j++) {
      bss[j].CP = bss[j].DP = fr;
    }
    for (j = midd+1; j <= rightd; j++) 
      bss[j].CP = bss[j].DP = -1;
    MP[0][fr] = -1;
    MP[1][fr] = -1;
    MP[2][fr] = -1;
    MT[0][fr] = MT[1][fr] = MT[2][fr] = 0;
  } else {
    fr = 0;
    for (j = 0; j < midd; j++) {
      bss[j].CP = bss[j].DP = 0;
    }
    for (j = midd; j <= rightd; j++) {
      bss[j].CP = bss[j].DP = 0;
    }
    MP[0][0] = -1;
    MP[1][0] = -1;
    MP[2][0] = -1;
    MT[0][0] = MT[1][0] = MT[2][0] = 0;
  }

  bss[leftd].CC = 0;
  if (tb == 2) t = 0;
  else t = -g;
  for (j = leftd+1; j <= rightd; j++) {
    bss[j].CC = t = t-h;
    bss[j].DD = t-g;
  }
  bss[rightd+1].CC = MININT;
  bss[rightd+1].DD = MININT;
  if (tb == 1) bss[leftd].DD = 0;
  else bss[leftd].DD = -g;
  bss[leftd-1].CC = MININT;
  for (i = 1; i <= M; i++) {
    if (i > N-up) rightd--;
    if (leftd > 1) leftd--;
    wa = w[A[i]];
    if ((c = bss[leftd+1].CC-m) > (d = bss[leftd+1].DD-h)) {
      d = c;
      bss[leftd].DP = bss[leftd+1].CP;
    } else bss[leftd].DP = bss[leftd+1].DP;
    if ((ib = leftd+low-1+i) > 0) c = bss[leftd].CC+wa[B[ib]];
    if (d > c || ib <= 0) {
      c = d;
      bss[leftd].CP = bss[leftd].DP;
    }
    e = c-g;
    bss[leftd].DD = d;
    bss[leftd].CC = c;
    IP = bss[leftd].CP;
    if (leftd == midd) bss[leftd].CP = bss[leftd].DP = IP = i;
    for (curd=leftd+1; curd <= rightd; curd++) {
      if (curd != midd) {
	if ((c = c-m) > (e = e-h)) {
	  e = c;
	  IP = bss[curd-1].CP;
	}  /* otherwise, IP is unchanged */
	if ((c = bss[curd+1].CC-m) > (d = bss[curd+1].DD-h)) {
	  d = c;
	  bss[curd].DP = bss[curd+1].CP;
	} else {
	  bss[curd].DP = bss[curd+1].DP;
	}
	c = bss[curd].CC + wa[B[curd+low-1+i]];
	if (c < d || c < e) {
	  if (e > d) {
	    c = e;
	    bss[curd].CP = IP;
	  } else {
	    c = d;
	    bss[curd].CP = bss[curd].DP;
	  }
	} /* otherwise, CP is unchanged */
	bss[curd].CC = c;
	bss[curd].DD = d;
      } else {
	if ((c = c-m) > (e = e-h)) {
	  e = c;
	  MP[1][i] = bss[curd-1].CP;
	  MT[1][i] = 2;
	} else {
	  MP[1][i] = IP;
	  MT[1][i] = 2;
	}
	if ((c = bss[curd+1].CC-m) > (d = bss[curd+1].DD-h)) {
	  d = c;
	  MP[2][i] = bss[curd+1].CP;
	  MT[2][i] = 1;
	} else {
	  MP[2][i] = bss[curd+1].DP;
	  MT[2][i] = 1;
	}
	c = bss[curd].CC + wa[B[curd+low-1+i]];
	if (c < d || c < e) {
	  if (e > d) {
	    c = e;
	    MP[0][i] = MP[1][i];
	    MT[0][i] = 2;
	  } else {
	    c = d;
	    MP[0][i] = MP[2][i];
	    MT[0][i] = 1;
	  }
	} else {
	  MP[0][i] = i-1;
	  MT[0][i] = 0;
	}
	if (c-g > e) {
	  MP[1][i] = MP[0][i];
	  MT[1][i] = MT[0][i];
	}
	if (c-g > d) {
	  MP[2][i] = MP[0][i];
	  MT[2][i] = MT[0][i];
	}
	bss[curd].CP = bss[curd].DP = IP = i;
	bss[curd].CC = c;
	bss[curd].DD = d;
      }
    }
  }

  /* decide which path to be traced back */
  if (te == 1 && d+g > c) {
    k = bss[rightd].DP;
    l = 2;
  } else if (te == 2 && e+g > c) {
    k = IP;
    l = 1;
  } else {
    k = bss[rightd].CP;
    l = 0;
  }
  if (rmid > N-M) l = 2;
  else if (rmid < N-M) l = 1;
  v = c;
  }
  /* Conquer: Solve subproblems recursively */

  /* trace back */
  r = -1;	
  for (; k > -1; r=k, k=MP[l][r], l=MT[l][r]){
    FP[k] = r;
    FT[k] = l;	/* l=0,1,2 */
  }
  /* forward dividing */
  if (r == -1) { /* optimal alignment did not cross the middle diagonal */
    if (rmid < 0) {
      bg_align(A,B,M,N,rmid+1,up,tb,te,w,g,h,bss, sapp, last);
    }
    else {
      bg_align(A,B,M,N,low,rmid-1,tb,te,w,g,h,bss, sapp, last);
    }
  } else {
    k = r;
    l = FP[k];
    kt = FT[k];

    /* first block */
    if (rmid < 0) {
      bg_align(A,B,r-1,r+rmid,rmid+1,min(up,r+rmid),tb,1,w,g,h,bss,sapp,last);
      DEL(1)
    } else if (rmid > 0) {
      bg_align(A,B,r,r+rmid-1,max(-r,low),rmid-1,tb,2,w,g,h,bss,sapp,last);
      INS(1)
    }

    /* intermediate blocks */
    t2 = up-rmid-1;
    t3 = low-rmid+1;
    for (; l > -1; k = l, l = FP[k], kt = FT[k]) {
      if (kt == 0) { REP }
      else if (kt == 1) { /* right-hand side triangle */
	INS(1)
	t1 = l-k-1;
	bg_align(A+k,B+k+rmid+1,t1,t1,0,min(t1,t2),2,1,w,g,h,bss,sapp,last);
	DEL(1)
      }
      else { /* kt == 2, left-hand side triangle */
	DEL(1)
	t1 = l-k-1;
	bg_align(A+k+1,B+k+rmid,t1,t1,max(-t1,t3),0,1,2,w,g,h,bss,sapp,last);
	INS(1)
      }
    }

    /* last block */
    if (N-M > rmid) {
      INS(1)
      t1 = k+rmid+1;
      bg_align(A+k,B+t1,M-k,N-t1,0,min(N-t1,t2),2,te,w,g,h,bss,sapp,last);
    } else if (N-M < rmid) {
      DEL(1)
      t1 = M-(k+1);
      bg_align(A+k+1,B+k+rmid,t1,N-(k+rmid),max(-t1,t3),0,1,te,w,g,h,
	       bss,sapp,last);
    }
  }
  return(v);
}

int B_ALIGN(const unsigned char *A, const unsigned char *B,
	    int M, int N,
	    int low, int up, int **W, int G, int H, int *S, int *nS,
	    int MW, int MX, struct bdstr *bss)
{ 
  int c, i, j;
  int g, h;
  size_t mj;
  int check_score;
  int **sapp, *sapp_v, *last, last_v;

  g = G;
  h = H;
  sapp_v = S;
  sapp = &sapp_v;

  last_v = 0;
  last = &last_v;

  low = min(max(-M, low),min(N-M,0));
  up = max(min(N, up),max(N-M,0));

  if (N <= 0) { 
    if (M > 0) { DEL(M); }
    return -gap(M);
  }
  if (M <= 0) {
    INS(N);
    return -gap(N);
  }
  if (up-low+1 <= 1) {
    c = 0;
    for (i = 1; i <= M; i++) {
      REP;
      c += W[A[i]][B[i]];
    }
    return c;
  }

  if (MT[0]==NULL) {
    mj = MX+1;
    MT[0] = (char *) ckalloc(mj);
    MT[1] = (char *) ckalloc(mj);
    MT[2] = (char *) ckalloc(mj);
    FT = (char *) ckalloc(mj);

    mj *= sizeof(int);
    MP[0] = (int *) ckalloc(mj);
    MP[1] = (int *) ckalloc(mj);
    MP[2] = (int *) ckalloc(mj);
    FP = (int *) ckalloc(mj);
  }

  c = bg_align(A,B,M,N,low,up,0,0,W,G,H,bss, sapp, last);

  check_score = BCHECK_SCORE(A,B,M,N,S,W,G,H,nS);

  free(FP); free(MP[2]); free(MP[1]); free(MP[0]);
  free(FT); free(MT[2]); free(MT[1]); free(MT[0]);
  MT[0]=NULL;

  if (check_score != c)
    printf("\nBCheck_score=%d != %d\n", check_score,c);
  return c;
}

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, optflag_s, optcut_s, optwid_s, n10, score;
  const unsigned char *aa1p;
  struct rstruct rst;

#ifdef TFASTA  
  f_str->n10 = n10=aatran(aa1,f_str->aa1x,n1,frame);
  do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, &rst, &hoff);
  aa1p = f_str->aa1x;

#else
  n10 = n1;
  aa1p = aa1;
#endif

  a_res->res = f_str->res;
  *have_ares = 1;

  if (ppst->sw_flag)
    return sw_walign(aa0, n0, aa1p, n10, ppst, f_str, a_res);
  else {
    optflag_s = ppst->param_u.fa.optflag;
    optcut_s = ppst->param_u.fa.optcut;
    optwid_s = ppst->param_u.fa.optwid;
    ppst->param_u.fa.optflag = 1;
    ppst->param_u.fa.optcut = 0;
    ppst->param_u.fa.optwid *= 2;

    do_fasta(aa0, n0, aa1p, n10, ppst, f_str, &rst, &hoff);

    if (rst.score[0]>0) {
      score=bd_walign(aa0, n0, aa1p, n10, ppst, f_str, hoff, a_res);
    }
    else {
      a_res->nres = 0;
      score=0;
    }

    ppst->param_u.fa.optflag = optflag_s;
    ppst->param_u.fa.optcut = optcut_s;
    ppst->param_u.fa.optwid = optwid_s;
    return score;
  }
}

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

#ifdef TFASTA
  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 TFASTA
  aln->qlfact = 1;
  aln->llfact = 3;
  aln->llmult = 3;
  aln->qlrev = 0;
  aln->frame = frame;
  if (frame > 2) {
    aln->llrev = 1;
    aln->frame = 3 - frame;
  }
  else aln->llrev = 0;
#else	/* FASTA */
  aln->llfact = aln->qlfact = aln->llmult = 1;
  aln->llrev = 0;
  if (frame > 0) aln->qlrev = 1;
  else 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;
  int op, lenc, nd, ns, itmp;
  const unsigned char *aa1p;
  char *sp0, *sp1, *spa, *sq;
  int *rp;
  int smins, mins;

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

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

  aln->amin0 = a_res.min0;
  aln->amax0 = a_res.max0;
  aln->amin1 = a_res.min1;
  aln->amax1 = a_res.max1;
  					/* will we show all the start ?*/
  if (min(a_res.min0,a_res.min1) < aln->llen || aln->showall==1)
    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->llcntx);
      aancpy(seqc0,(char *)aa0+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->llcntx);
      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 *)aa0,a_res.min0,pst);
	aln->smin0 = 0;
      }
      else {
	aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
	aln->smin0 = a_res.min0-mins;
      }
    }
  else {
    mins= min(aln->llcntx,min(a_res.min0,a_res.min1));
    smins=mins;
    aln->smin0=a_res.min0 - mins;
    aln->smin1=a_res.min1 - mins;
    aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
    aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
  }
  /* set the alignment code to zero for context */
  memset(seqca,0,mins);

  /* TFASTA
  smins = mins = 0;
  aln->smin0=a_res.min0;
  aln->smin1=a_res.min1;
  */

/* now get the middle */

  spa = seqca+mins;
  sp0 = seqc0+mins;
  sp1 = seqc1+mins;
  rp = a_res.res;
  lenc = 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++;
      lenc++;
      if ((itmp=pst.pam2[0][aa0[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
      else if (itmp == 0) { *spa = M_ZERO;}
      else {*spa = M_POS;}
      if (*spa == M_POS || *spa==M_ZERO) aln->nsim++;


      *sp0 = sq[aa0[i0++]];
      *sp1 = sq[aa1p[i1++]];
      if (toupper(*sp0) == toupper(*sp1)) {
	aln->nident++;
	*spa = M_IDENT;
      }
      else if (pst.nt_align) {
	if ((toupper(*sp0) == 'T' && toupper(*sp1) == 'U') ||
	    (toupper(*sp0)=='U' && toupper(*sp1)=='T')) {
	  aln->nident++;
	  *spa = M_IDENT;
	}
	else if (toupper(*sp0) == 'N') aln->ngap_q++;
	else if (toupper(*sp1) == 'N') aln->ngap_l++;
      }
      sp0++; sp1++; spa++;
    }
    else {
      if (op==0) op = *rp++;
      if (op > 0) {
	*sp0++ = '-';
	*sp1++ = sq[aa1p[i1++]];
	*spa++ = M_DEL;
	op--;
	lenc++;
	aln->ngap_q++;
      }
      else {
	*sp0++ = sq[aa0[i0++]];
	*sp1++ = '-';
	*spa++ = M_DEL;
	op++;
	lenc++;
	aln->ngap_l++;
      }
    }
  }

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

/*	now we have the middle, get the right end */
  if (!aln->llcntx_flg) {
    ns = mins + lenc + aln->llen;	/* show an extra line? */
    ns -= (itmp = ns %aln->llen);	/* itmp = left over on last line */
    if (itmp>aln->llen/2) ns += aln->llen;  /* more than 1/2 , use another*/
    nd = ns - (mins+lenc);		/* this much extra */
  }
  else nd = aln->llcntx;

  if (nd > max(n0-a_res.max0,nn1-a_res.max1)) 
    nd = max(n0-a_res.max0,nn1-a_res.max1);
  
  if (aln->showall==1) {
    nd = max(n0-a_res.max0,nn1-a_res.max1);	/* reset for showall=1 */
    /* get right end */
    aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-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+n0-a_res.max0,' ',nd-(n0-a_res.max0));
    memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
  }
  else {
    if ((nd-(n0-a_res.max0))>0) {
      aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,(n0-a_res.max0),pst);
      memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
    }
    else aancpy(seqc0+mins+lenc,(char *)aa0+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);
  }

  /*  fprintf(stderr,"%d\n",mins+lenc+nd); */

  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;
  int op, lenc, nd, ns, itmp;
  const unsigned char *aa1p;
  char *sp0, *sp0a, *sp1, *spa, *sq;
  int *rp;
  int smins, mins;

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

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

  aln->amin0 = a_res.min0;
  aln->amax0 = a_res.max0;
  aln->amin1 = a_res.min1;
  aln->amax1 = a_res.max1;
  					/* will we show all the start ?*/
  /* will we show all the start ?*/
  if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1)
    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->llcntx);
      aancpy(seqc0,(char *)aa0+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->llcntx);
      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 *)aa0,a_res.min0,pst);
	aln->smin0 = 0;
      }
      else {
	aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
	aln->smin0 = a_res.min0-mins;
      }
    }
  else {
    mins= min(aln->llcntx,min(a_res.min0,a_res.min1));
    smins=mins;
    aln->smin0=a_res.min0 - smins;
    aln->smin1=a_res.min1 - smins;
    aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
    aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
  }
  /* set the alignment code to zero for context */
  memset(seqca,0,mins);
  memset(seqc0a,' ',mins);

  /* TFASTA
  smins = mins = 0;
  aln->smin0=a_res.min0;
  aln->smin1=a_res.min1;
  */

/* now get the middle */

  spa = seqca+mins;
  sp0 = seqc0+mins;
  sp0a = seqc0a+mins;
  sp1 = seqc1+mins;
  rp = a_res.res;
  lenc = 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++;
      lenc++;
      if ((itmp=pst.pam2[0][aa0[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
      else if (itmp == 0) { *spa = M_ZERO;}
      else {*spa = M_POS;}
      if (*spa == M_POS || *spa==M_ZERO) aln->nsim++;

      *sp0a++ = ann_arr[aa0a[i0]];

      *sp0 = sq[aa0[i0++]];
      *sp1 = sq[aa1p[i1++]];

      if (toupper(*sp0) == toupper(*sp1)) {
	aln->nident++;
	*spa = M_IDENT;
      }
      else if (pst.nt_align) {
	if ((toupper(*sp0) == 'T' && toupper(*sp1) == 'U') ||
	    (toupper(*sp0)=='U' && toupper(*sp1)=='T')) {
	  aln->nident++;
	  *spa = M_IDENT;
	}
	else if (toupper(*sp0) == 'N') aln->ngap_q++;
	else if (toupper(*sp1) == 'N') aln->ngap_l++;
      }
      sp0++; sp1++; spa++;
    }
    else {
      if (op==0) op = *rp++;
      if (op>0) {
	*sp0++ = '-';
	*sp1++ = sq[aa1p[i1++]];
	*spa++ = M_DEL;
	*sp0a++ = ' ';
	op--;
	lenc++;
	aln->ngap_q++;
      }
      else {
	*sp0a++ = ann_arr[aa0a[i0]];
	*sp0++ = sq[aa0[i0++]];
	*sp1++ = '-';
	*spa++ = M_DEL;
	op++;
	lenc++;
	aln->ngap_l++;
      }
    }
  }

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

/*	now we have the middle, get the right end */
  if (!aln->llcntx_flg) {
    ns = mins + lenc + aln->llen;	/* show an extra line? */
    ns -= (itmp = ns %aln->llen);	/* itmp = left over on last line */
    if (itmp>aln->llen/2) ns += aln->llen;  /* more than 1/2 , use another*/
    nd = ns - (mins+lenc);		/* this much extra */
  }
  else nd = aln->llcntx;

  if (nd > max(n0-a_res.max0,nn1-a_res.max1)) 
    nd = max(n0-a_res.max0,nn1-a_res.max1);
  
  if (aln->showall==1) {
    nd = max(n0-a_res.max0,nn1-a_res.max1);	/* reset for showall=1 */
    /* get right end */
    aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-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+n0-a_res.max0,' ',nd-(n0-a_res.max0));
    memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
  }
  else {
    if ((nd-(n0-a_res.max0))>0) {
      aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,(n0-a_res.max0),pst);
      memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
    }
    else aancpy(seqc0+mins+lenc,(char *)aa0+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);
  }

  /*  fprintf(stderr,"%d\n",mins+lenc+nd); */

  return mins+lenc+nd;
}

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

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

  sprintf(tmp_cnt,"%c%d",op_char[op],op_cnt);
  strncat(al_str,tmp_cnt,al_str_max);
}


/* build an array of match/ins/del - length strings */
int calc_code(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,
	      char *al_str, int al_str_n, struct f_struct *f_str)
{
  int i0, i1, nn1;
  int op, lenc;
  int p_op, op_cnt;
  const unsigned char *aa1p;
  char sp0, sp1, *sq;
  int *rp;

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

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

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

  rp = a_res.res;
  lenc = 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;
  
  while (i0 < a_res.max0 || i1 < a_res.max1) {
    if (op == 0 && *rp == 0) {

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

      sp0 = sq[aa0[i0++]];
      sp1 = sq[aa1p[i1++]];

      if (p_op == 0 || p_op==3) {
	if (sp0 != '*' && sp1 != '*') {
	  if (p_op == 3) {
	    update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
	    op_cnt = 1; p_op = 0;
	  }
	  else {op_cnt++;}
	}
	else {
	  update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
	  op_cnt = 1; p_op = 3;
	}
      }
      else {
	update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
	op_cnt = 1; p_op = 0;
      }

      op = *rp++;
      lenc++;

      if (toupper(sp0) == toupper(sp1)) aln->nident++;
      else if (pst.nt_align) {
	if ((toupper(sp0) == 'T' && toupper(sp1) == 'U') ||
	    (toupper(sp0)=='U' && toupper(sp1)=='T')) aln->nident++;
	else if (toupper(sp0) == 'N') aln->ngap_q++;
	else if (toupper(sp1) == 'N') aln->ngap_l++;
      }
    }
    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);
	  op_cnt = 1; p_op = 1;
	}
	op--; lenc++; i1++; aln->ngap_q++;
      }
      else {
	if (p_op == 2) { op_cnt++;}
	else {
	  update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
	  op_cnt = 1; p_op = 2;
	}
	op++; lenc++; i0++; aln->ngap_l++;
      }
    }
  }
  update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);

  return lenc;
}

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;
  int sp0, sp1;
  const unsigned char *aa1p;
  int *rp;
  char *sq;
  
  if (pst.ext_sq_set) { sq = pst.sqx; }
  else { sq = pst.sq; }

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

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

  rp = a_res.res;
  lenc = 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++;
      lenc++;

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

      sp0 = sq[aa0[i0++]];
      sp1 = sq[aa1p[i1++]];
      if (toupper(sp0) == toupper(sp1)) {aln->nident++;}
      else if (pst.nt_align) {
	if ((toupper(sp0)=='T' && toupper(sp1)== 'U')||
	  (toupper(sp0)=='U' && toupper(sp1)=='T')) {aln->nident++;}
	else if (toupper(sp0) == 'N') aln->ngap_q++;
	else if (toupper(sp1) == 'N') aln->ngap_l++;
      }
    }
    else {
      if (op==0) op = *rp++;
      if (op>0) {op--; lenc++; i1++; aln->ngap_q++;}
      else {op++; lenc++; i0++; aln->ngap_l++; }
    }
  }
  return lenc;
}

#ifdef PCOMPLIB

#include "w_mw.h"

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