binaries/src/fasta34/dropnfa.c

   1
   2 /* copyright (c) 1998, 1999 William R. Pearson and the U. of Virginia */
   3
   4 /* $Name: fa_34_26_5 $ - $Id: dropnfa.c,v 1.81 2007/04/26 18:37:19 wrp Exp $ */
   5
   6 /* 18-Sep-2006 - removed global variables for alignment from nw_align
   7    and bg_align */
   8
   9 /* 18-Oct-2005 - converted to use a_res and aln for alignment coordinates */
  10
  11 /* 14-May-2003 - modified to return alignment start at 0, rather than
  12    1, for begin:end alignments
  13 */
  14
  15 /*
  16   implements the fasta algorithm, see:
  17
  18   W. R. Pearson, D. J. Lipman (1988) "Improved tools for biological
  19   sequence comparison" Proc. Natl. Acad. Sci. USA 85:2444-2448
  20
  21   This version uses Smith-Waterman for final protein alignments
  22
  23   W. R. Pearson (1996) "Effective protein sequence comparison"
  24   Methods Enzymol. 266:227-258
  25
  26
  27   26-April-2001 - -DGAP_OPEN redefines -f, as gap open penalty
  28
  29   4-Nov-2001 - modify spam() while(1).
  30 */
  31
  32 #include <stdio.h>
  33 #include <stdlib.h>
  34 #include <string.h>
  35 #include <ctype.h>
  36 #include <math.h>
  37
  38 #include "defs.h"
  39 #include "param.h"
  40
  41 /* this must be consistent with upam.h */
  42 #define MAXHASH 32
  43 #define NMAP MAXHASH+1
  44
  45 /* globals for fasta */
  46 #define MAXWINDOW 64
  47
  48 #ifndef MAXSAV
  49 #define MAXSAV 10
  50 #endif
  51
  52 #ifndef ALLOCN0
  53 static char *verstr="3.5 Sept 2006";
  54 #else
  55 static char *verstr="3.5an0 Sept 2006";
  56 #endif
  57
  58 extern void w_abort(char *, char *);
  59 int shscore(const unsigned char *aa0, int n0, int **pam2);
  60 extern void init_karlin(const unsigned char *aa0, int n0, struct pstruct *ppst,
  61                         double *aa0_f, double **kp);
  62 extern void init_karlin_a(struct pstruct *, double *, double **);
  63 extern int do_karlin(const unsigned char *, int n1, int **,
  64                      struct pstruct *, double *, double *,
  65                      double *, double *);
  66 extern void aancpy(char *to, char *from, int count, struct pstruct pst);
  67 char *ckalloc(size_t);
  68
  69 #ifdef TFASTA
  70 extern int aatran(const unsigned char *ntseq, unsigned char *aaseq, int maxs, int frame);
  71 #endif
  72
  73 #include "dropnfa.h"
  74
  75 #define DROP_INTERN
  76 #include "drop_func.h"
  77
  78 struct swstr { int H, E;};
  79
  80 static int
  81 dmatch (const unsigned char *aa0, int n0,
  82         const unsigned char *aa1, int n1,
  83         int hoff, int window,
  84         int **pam2, int gdelval, int ggapval,
  85         struct f_struct *f_str);
  86
  87 /* initialize for fasta */
  88
  89 void
  90 init_work (unsigned char *aa0, int n0,
  91            struct pstruct *ppst,
  92            struct f_struct **f_arg)
  93 {
  94    int mhv, phv;
  95    int hmax;
  96    int i0, hv;
  97    int pamfact;
  98    int btemp;
  99    struct f_struct *f_str;
 100    /* these used to be globals, but do not need to be */
 101    int ktup;            /* word size examined */
 102    int fact;            /* factor used to scale ktup match value */
 103    int kt1;             /* ktup-1 */
 104    int lkt;             /* last ktup - initiall kt1, but can be increased
 105                            for hsq >= NMAP */
 106
 107    int maxn0;           /* used in band alignment */
 108    int *pwaa;           /* pam[aa0[]] profile */
 109    int i, j;
 110    struct swstr *ss;
 111    int *waa;
 112    int nsq, ip, *hsq;
 113
 114    if (ppst->ext_sq_set) {
 115      nsq = ppst->nsqx; ip = 1;
 116      hsq = ppst->hsqx;
 117    }
 118    else {
 119      nsq = ppst->nsq; ip = 0;
 120      hsq = ppst->hsq;
 121    }
 122
 123   f_str = (struct f_struct *)calloc(1,sizeof(struct f_struct));
 124
 125 #ifndef TFASTA
 126   if((ppst->zsflag%10) == 6) {
 127     f_str->kar_p = NULL;
 128     init_karlin(aa0, n0, ppst, &f_str->aa0_f[0], &f_str->kar_p);
 129   }
 130 #endif
 131
 132   btemp = 2 * ppst->param_u.fa.bestoff / 3 +
 133     n0 / ppst->param_u.fa.bestscale +
 134     ppst->param_u.fa.bkfact *
 135     (ppst->param_u.fa.bktup - ppst->param_u.fa.ktup);
 136
 137   if (ppst->nt_align)
 138      btemp = (btemp*ppst->pam_h)/5;  /* normalize to standard +5/-4 */
 139
 140    btemp = min (btemp, ppst->param_u.fa.bestmax);
 141    if (btemp > 3 * n0) btemp = 3 * shscore(aa0,n0,ppst->pam2[0]) / 5;
 142
 143    ppst->param_u.fa.cgap = btemp + ppst->param_u.fa.bestoff / 3;
 144
 145    if (ppst->param_u.fa.optcut_set != 1)
 146 #ifndef TFASTA
 147       ppst->param_u.fa.optcut = btemp;
 148 #else
 149       ppst->param_u.fa.optcut = (btemp*3)/2;
 150 #endif
 151
 152 #ifndef OLD_FASTA_GAP
 153    ppst->param_u.fa.pgap = ppst->gdelval + 2*ppst->ggapval;
 154 #else
 155    ppst->param_u.fa.pgap = ppst->gdelval + ppst->ggapval;
 156 #endif
 157    pamfact = ppst->param_u.fa.pamfact;
 158    ktup = ppst->param_u.fa.ktup;
 159    fact = ppst->param_u.fa.scfact * ktup;
 160
 161    if (pamfact == -1) pamfact = 0;
 162    else if (pamfact == -2) pamfact = 1;
 163
 164    for (i0 = 1, mhv = -1; i0 <= ppst->nsq; i0++)
 165       if (hsq[i0] < NMAP && hsq[i0] > mhv) mhv = hsq[i0];
 166
 167    if (mhv <= 0) {
 168       fprintf (stderr, " maximum hsq <=0 %d\n", mhv);
 169       exit (1);
 170    }
 171
 172    for (f_str->kshft = 0; mhv > 0; mhv /= 2) f_str->kshft++;
 173
 174 /*      kshft = 2;      */
 175    kt1 = ktup - 1;
 176    hv = 1;
 177    for (i0 = 0; i0 < ktup; i0++) hv = hv << f_str->kshft;
 178    hmax = hv;
 179    f_str->hmask = (hmax >> f_str->kshft) - 1;
 180
 181    if ((f_str->harr = (int *) calloc (hmax, sizeof (int))) == NULL) {
 182      fprintf (stderr, " cannot allocate hash array: hmax: %d hmask: %d\n",
 183               hmax, f_str->hmask);
 184      exit (1);
 185    }
 186
 187    if ((f_str->pamh1 = (int *) calloc (nsq+1, sizeof (int))) == NULL) {
 188      fprintf (stderr, " cannot allocate pamh1 array nsq=%d\n",nsq);
 189      exit (1);
 190    }
 191
 192    if ((f_str->pamh2 = (int *) calloc (hmax, sizeof (int))) == NULL) {
 193      fprintf (stderr, " cannot allocate pamh2 array hmax=%d\n",hmax);
 194      exit (1);
 195    }
 196
 197    if ((f_str->link = (int *) calloc (n0, sizeof (int))) == NULL) {
 198      fprintf (stderr, " cannot allocate hash link array n0=%d",n0);
 199      exit (1);
 200    }
 201
 202    for (i0 = 0; i0 < hmax; i0++) f_str->harr[i0] = -1;
 203    for (i0 = 0; i0 < n0; i0++) f_str->link[i0] = -1;
 204
 205    /* encode the aa0 array */
 206    phv = hv = 0;
 207    lkt = kt1;
 208    /* restart hv, phv calculation */
 209    for (i0 = 0; i0 < min(lkt,n0); i0++) {
 210      if (hsq[aa0[i0]] >= NMAP) {hv=phv=0; lkt = i0+ ktup; continue;}
 211      hv = (hv << f_str->kshft) + hsq[aa0[i0]];
 212      phv += ppst->pam2[ip][aa0[i0]][aa0[i0]]*ktup;
 213    }
 214
 215    for (; i0 < n0; i0++) {
 216      if (hsq[aa0[i0]] >= NMAP) {
 217        hv=phv=0;
 218        /* restart hv, phv calculation */
 219        for (lkt = i0+kt1; (i0 < lkt || hsq[aa0[i0]]>=NMAP) && i0<n0; i0++) {
 220          if (hsq[aa0[i0]] >= NMAP) {
 221            hv=phv=0;
 222            lkt = i0+ktup;
 223            continue;
 224          }
 225          hv = (hv << f_str->kshft) + hsq[aa0[i0]];
 226          phv += ppst->pam2[ip][aa0[i0]][aa0[i0]]*ktup;
 227        }
 228      }
 229      if (i0 >= n0) break;
 230      hv = ((hv & f_str->hmask) << f_str->kshft) + hsq[aa0[i0]];
 231      f_str->link[i0] = f_str->harr[hv];
 232      f_str->harr[hv] = i0;
 233      if (pamfact) {
 234        f_str->pamh2[hv] = (phv += ppst->pam2[ip][aa0[i0]][aa0[i0]] * ktup);
 235        /* this check should always be true, but just in case */
 236        if (hsq[aa0[i0-kt1]]<NMAP)
 237          phv -= ppst->pam2[ip][aa0[i0 - kt1]][aa0[i0 - kt1]] * ktup;
 238      }
 239      else f_str->pamh2[hv] = fact * ktup;
 240    }
 241
 242 /* this has been modified from 0..<ppst->nsq to 1..<=ppst->nsq because the
 243    pam2[0][0] is now undefined for consistency with blast
 244 */
 245
 246    if (pamfact)
 247       for (i0 = 1; i0 <= nsq; i0++)
 248          f_str->pamh1[i0] = ppst->pam2[ip][i0][i0] * ktup;
 249    else
 250       for (i0 = 1; i0 <= nsq; i0++)
 251          f_str->pamh1[i0] = fact;
 252
 253    f_str->ndo = 0;
 254    f_str->noff = n0-1;
 255 #ifndef ALLOCN0
 256    if ((f_str->diag = (struct dstruct *) calloc ((size_t)MAXDIAG,
 257                                                  sizeof (struct dstruct)))==NULL) {
 258       fprintf (stderr," cannot allocate diagonal arrays: %lu\n",
 259               MAXDIAG *sizeof (struct dstruct));
 260       exit (1);
 261      };
 262 #else
 263    if ((f_str->diag = (struct dstruct *) calloc ((size_t)n0,
 264                                               sizeof (struct dstruct)))==NULL) {
 265       fprintf (stderr," cannot allocate diagonal arrays: %ld\n",
 266               (long)n0*sizeof (struct dstruct));
 267       exit (1);
 268      };
 269 #endif
 270
 271
 272 #ifdef TFASTA
 273    if ((f_str->aa1x =(unsigned char *)calloc((size_t)ppst->maxlen+2,
 274                                              sizeof(unsigned char)))
 275        == NULL) {
 276      fprintf (stderr, "cannot allocate aa1x array %d\n", ppst->maxlen+2);
 277      exit (1);
 278    }
 279    f_str->aa1x++;
 280 #endif
 281
 282    f_str->bss = (struct bdstr *) calloc((size_t)ppst->param_u.fa.optwid*2+4,
 283                                         sizeof(struct bdstr));
 284    f_str->bss++;
 285
 286    /* allocate space for the scoring arrays */
 287    maxn0 = n0 + 4;
 288    if ((ss = (struct swstr *) calloc (maxn0, sizeof (struct swstr)))
 289          == NULL) {
 290      fprintf (stderr, "cannot allocate ss array %3d\n", n0);
 291      exit (1);
 292    }
 293    ss++;
 294    f_str->ss = ss;
 295
 296    /* initialize the "variable" pam array */
 297
 298    if ((waa= (int *)calloc ((size_t)(nsq+1)*n0,sizeof(int))) == NULL) {
 299      fprintf(stderr,"cannot allocate waa struct %3d\n",nsq*n0);
 300      exit(1);
 301    }
 302
 303    pwaa = waa;
 304    for (i=0; i<=nsq; i++) {
 305      for (j=0;j<n0; j++) {
 306        *pwaa = ppst->pam2[ip][aa0[j]][i];
 307        pwaa++;
 308      }
 309    }
 310    f_str->waa0 = waa;
 311
 312    /* initialize the "conventional" pam array used for alignments */
 313
 314    if ((waa= (int *)calloc ((size_t)(nsq+1)*n0,sizeof(int))) == NULL) {
 315      fprintf(stderr,"cannot allocate waa struct %3d\n",nsq*n0);
 316      exit(1);
 317    }
 318
 319    pwaa = waa;
 320    for (i=0; i<=nsq; i++) {
 321      for (j=0;j<n0; j++) {
 322        *pwaa = ppst->pam2[0][aa0[j]][i];
 323        pwaa++;
 324      }
 325    }
 326    f_str->waa1 = waa;
 327
 328    f_str->max_res = max(3*n0/2,MIN_RES);
 329
 330    /* now we need alignment storage - get it */
 331    if ((f_str->res = (int *)calloc((size_t)f_str->max_res,sizeof(int)))==NULL) {
 332      fprintf(stderr,"cannot allocate alignment results array %d\n",f_str->max_res);
 333      exit(1);
 334    }
 335
 336    *f_arg = f_str;
 337 }
 338
 339
 340 /* pstring1 is a message to the manager, currently 512 */
 341 /* pstring2 is the same information, but in a markx==10 format */
 342 void
 343 get_param (struct pstruct *pstr, char *pstring1, char *pstring2)
 344 {
 345 #ifndef TFASTA
 346   char *pg_str="FASTA";
 347 #else
 348   char *pg_str="TFASTA";
 349 #endif
 350
 351    if (!pstr->param_u.fa.optflag)
 352 #ifdef OLD_FASTA_GAP
 353      sprintf (pstring1, "%s (%s) function [%s matrix, (%d:%d)%s] ktup: %d\n join: %d, gap-pen: %d/%d, width: %3d",
 354 #else
 355      sprintf (pstring1, "%s (%s) function [%s matrix, (%d:%d)%s] ktup: %d\n join: %d, open/ext: %d/%d, width: %3d",
 356 #endif
 357               pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l,
 358               (pstr->ext_sq_set) ? "xS":"\0",
 359               pstr->param_u.fa.ktup, pstr->param_u.fa.cgap,
 360                pstr->gdelval, pstr->ggapval, pstr->param_u.fa.optwid);
 361    else
 362 #ifdef OLD_FASTA_GAP
 363       sprintf (pstring1, "%s (%s) function [optimized, %s matrix (%d:%d)%s] ktup: %d\n join: %d, opt: %d, gap-pen: %d/%d, width: %3d",
 364 #else
 365       sprintf (pstring1, "%s (%s) function [optimized, %s matrix (%d:%d)%s] ktup: %d\n join: %d, opt: %d, open/ext: %d/%d, width: %3d",
 366 #endif
 367                pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l,
 368                (pstr->ext_sq_set) ? "xS":"\0",
 369                pstr->param_u.fa.ktup, pstr->param_u.fa.cgap,
 370                pstr->param_u.fa.optcut, pstr->gdelval, pstr->ggapval,
 371                pstr->param_u.fa.optwid);
 372    if (pstr->param_u.fa.iniflag) strcat(pstring1," init1");
 373    /*
 374    if (pstr->zsflag==0) strcat(pstring1," not-scaled");
 375    else if (pstr->zsflag==1) strcat(pstring1," reg.-scaled");
 376    */
 377
 378    if (pstring2 != NULL) {
 379 #ifdef OLD_FASTA_GAP
 380      sprintf (pstring2, "; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)\n\
 381 ; pg_gap-pen: %d %d\n; pg_ktup: %d\n; pg_optcut: %d\n; pg_cgap: %d\n",
 382 #else
 383      sprintf (pstring2, "; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)\n\
 384 ; pg_open-ext: %d %d\n; pg_ktup: %d\n; pg_optcut: %d\n; pg_cgap: %d\n",
 385 #endif
 386               pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l, pstr->gdelval,
 387               pstr->ggapval,pstr->param_u.fa.ktup,pstr->param_u.fa.optcut,
 388               pstr->param_u.fa.cgap);
 389    }
 390 }
 391
 392 void
 393 close_work (const unsigned char *aa0, int n0,
 394             struct pstruct *ppst,
 395             struct f_struct **f_arg)
 396 {
 397   struct f_struct *f_str;
 398
 399
 400   f_str = *f_arg;
 401
 402
 403   if (f_str != NULL) {
 404     if (f_str->kar_p!=NULL) free(f_str->kar_p);
 405     f_str->ss--;
 406     f_str->bss--;
 407
 408     free(f_str->res);
 409     free(f_str->waa1);
 410     free(f_str->waa0);
 411     free(f_str->ss);
 412     free(f_str->bss);
 413     free(f_str->diag);
 414     free(f_str->link);
 415     free(f_str->pamh2);
 416     free(f_str->pamh1);
 417     free(f_str->harr);
 418
 419     free(f_str);
 420     *f_arg = NULL;
 421   }
 422 }
 423
 424 #ifdef ALLOCN0
 425 void savemax (struct dstruct *, int, struct f_struct *);
 426 #else
 427 void savemax (struct dstruct *, struct f_struct *);
 428 #endif
 429
 430 int spam (const unsigned char *, const unsigned char *, struct savestr *,
 431          int **, int, int, int);
 432 int sconn(struct savestr **, int nsave, int cgap, int pgap, int noff);
 433 void kpsort(struct savestr **, int);
 434
 435 static int
 436 ALIGN(const unsigned char *, const unsigned char *, int, int,
 437       int **, int, int, int *, int *, struct f_struct *);
 438
 439 static int
 440 LOCAL_ALIGN(const unsigned char *, const unsigned char *,
 441             int, int, int, int,
 442             int **, int, int, int *, int *, int *, int *, int,
 443             struct f_struct *);
 444
 445 static int
 446 B_ALIGN(const unsigned char *A, const unsigned char *B, int M,
 447         int N, int low, int up, int **W, int G, int H, int *S,
 448         int *nS, int MW, int MX, struct bdstr *bss);
 449
 450 static void
 451 do_fasta (const unsigned char *aa0, int n0,
 452           const unsigned char *aa1, int n1,
 453           struct pstruct *ppst, struct f_struct *f_str,
 454           struct rstruct *rst, int *hoff)
 455 {
 456    int     nd;          /* diagonal array size */
 457    int     lhval;
 458    int     kfact;
 459    register struct dstruct *dptr;
 460    register int tscor;
 461
 462 #ifndef ALLOCN0
 463    register struct dstruct *diagp;
 464 #else
 465    register int dpos;
 466    int     lposn0;
 467 #endif
 468    int noff;
 469    struct dstruct *dpmax;
 470    register int lpos;
 471    int     tpos;
 472    struct savestr *vmptr;
 473    int     scor, ib, nsave;
 474    int xdrop, do_extend;
 475    int ktup, kt1, lkt, *hsq, ip;
 476
 477   if (ppst->ext_sq_set) {
 478     ip = 1;
 479     hsq = ppst->hsqx;
 480   }
 481   else {
 482     ip = 0;
 483     hsq = ppst->hsq;
 484   }
 485
 486   xdrop = -ppst->pam_l;
 487   /* do extended alignment in spam iff protein or short sequences */
 488   do_extend = !ppst->nt_align || (n0 < 50) || (n1 < 50);
 489
 490   ktup = ppst->param_u.fa.ktup;
 491   kt1 = ktup-1;
 492
 493   if (n1 < ktup) {
 494      rst->score[0] = rst->score[1] = rst->score[2] = 0;
 495      return;
 496    }
 497
 498    if (n0+n1+1 >= MAXDIAG) {
 499      fprintf(stderr,"n0,n1 too large: %d, %d\n",n0,n1);
 500      rst->score[0] = rst->score[1] = rst->score[2] = -1;
 501      return;
 502    }
 503
 504 #ifdef ALLOCN0
 505    nd = n0;
 506 #else
 507    nd = n0 + n1;
 508 #endif
 509
 510    dpmax = &f_str->diag[nd];
 511    for (dptr = &f_str->diag[f_str->ndo]; dptr < dpmax;)
 512    {
 513       dptr->stop = -1;
 514       dptr->dmax = NULL;
 515       dptr++->score = 0;
 516    }
 517
 518    for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
 519       vmptr->score = 0;
 520    f_str->lowmax = f_str->vmax;
 521    f_str->lowscor = 0;
 522
 523    /* start hashing */
 524    lhval = 0;
 525    lkt = kt1;
 526    for (lpos = 0; (lpos < lkt || hsq[aa1[lpos]]>=NMAP) && lpos <n1; lpos++) {
 527      /* restart lhval calculation */
 528      if (hsq[aa1[lpos]]>=NMAP) {
 529        lhval = 0; lkt = lpos + ktup;
 530        continue;
 531      }
 532      lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
 533    }
 534
 535    noff = f_str->noff;
 536 #ifndef ALLOCN0
 537    diagp = &f_str->diag[noff + lkt];
 538    for (; lpos < n1; lpos++, diagp++) {
 539      if (hsq[aa1[lpos]]>=NMAP) {
 540        lpos++ ; diagp++;
 541        while (lpos < n1 && hsq[aa1[lpos]]>=NMAP) {lpos++; diagp++;}
 542        if (lpos >= n1) break;
 543        lhval = 0;
 544      }
 545      lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
 546      for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
 547        if ((tscor = (dptr = &diagp[-tpos])->stop) >= 0) {
 548 #else
 549    lposn0 = noff + lpos;
 550    for (; lpos < n1; lpos++, lposn0++) {
 551      if (hsq[aa1[lpos]]>=NMAP) {lhval = 0; goto loopl;}
 552      /*
 553      if (hsq[aa1[lpos]]>=NMAP) {
 554        lpos++; lposn0++;
 555        while (lpos < n1 && hsq[aa1[lpos]]>=NMAP) {lpos++; lposn0++;}
 556      }
 557      */
 558      lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
 559      for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
 560        dpos = lposn0 - tpos;
 561        if ((tscor = (dptr = &f_str->diag[dpos % nd])->stop) >= 0) {
 562 #endif
 563          tscor += ktup;
 564          if ((tscor -= lpos) <= 0) {
 565            scor = dptr->score;
 566            if ((tscor += (kfact = f_str->pamh2[lhval])) < 0 && f_str->lowscor < scor)
 567 #ifdef ALLOCN0
 568              savemax (dptr, dpos, f_str);
 569 #else
 570              savemax (dptr, f_str);
 571 #endif
 572              if ((tscor += scor) >= kfact) {
 573                dptr->score = tscor;
 574                dptr->stop = lpos;
 575              }
 576              else {
 577                dptr->score = kfact;
 578                dptr->start = (dptr->stop = lpos) - kt1;
 579              }
 580          }
 581          else {
 582            dptr->score += f_str->pamh1[aa0[tpos]];
 583            dptr->stop = lpos;
 584          }
 585        }
 586        else {
 587          dptr->score = f_str->pamh2[lhval];
 588          dptr->start = (dptr->stop = lpos) - kt1;
 589        }
 590      }                          /* end tpos */
 591
 592 #ifdef ALLOCN0
 593       /* reinitialize diag structure */
 594    loopl:
 595      if ((dptr = &f_str->diag[lpos % nd])->score > f_str->lowscor)
 596        savemax (dptr, lpos, f_str);
 597      dptr->stop = -1;
 598      dptr->dmax = NULL;
 599      dptr->score = 0;
 600 #endif
 601    }                            /* end lpos */
 602
 603 #ifdef ALLOCN0
 604    for (tpos = 0, dpos = noff + n1 - 1; tpos < n0; tpos++, dpos--) {
 605      if ((dptr = &f_str->diag[dpos % nd])->score > f_str->lowscor)
 606        savemax (dptr, dpos, f_str);
 607    }
 608 #else
 609    for (dptr = f_str->diag; dptr < dpmax;) {
 610      if (dptr->score > f_str->lowscor) savemax (dptr, f_str);
 611      dptr->stop = -1;
 612      dptr->dmax = NULL;
 613      dptr++->score = 0;
 614    }
 615    f_str->ndo = nd;
 616 #endif
 617
 618 /*
 619         at this point all of the elements of aa1[lpos]
 620         have been searched for elements of aa0[tpos]
 621         with the results in diag[dpos]
 622 */
 623    for (nsave = 0, vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)    {
 624      /*
 625      fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
 626              noff+vmptr->start-vmptr->dp,
 627              noff+vmptr->stop-vmptr->dp,
 628              vmptr->start,vmptr->stop,
 629              vmptr->dp,vmptr->score);
 630
 631      */
 632       if (vmptr->score > 0) {
 633         vmptr->score = spam (aa0, aa1, vmptr, ppst->pam2[ip], xdrop,
 634                              noff,do_extend);
 635         f_str->vptr[nsave++] = vmptr;
 636       }
 637    }
 638
 639    if (nsave <= 0) {
 640      rst->score[0] = rst->score[1] = rst->score[2] = 0;
 641      return;
 642    }
 643
 644    /*
 645    fprintf(stderr,"n0: %d; n1: %d; noff: %d\n",n0,n1,noff);
 646    for (ib=0; ib<nsave; ib++) {
 647      fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
 648              noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
 649              noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
 650              f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
 651              f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
 652    }
 653    fprintf(stderr,"---\n");
 654    */
 655
 656    scor = sconn (f_str->vptr, nsave, ppst->param_u.fa.cgap,
 657                  ppst->param_u.fa.pgap, noff);
 658
 659    for (vmptr=f_str->vptr[0],ib=1; ib<nsave; ib++)
 660      if (f_str->vptr[ib]->score > vmptr->score) vmptr=f_str->vptr[ib];
 661
 662 /*  kssort (f_str->vptr, nsave); */
 663
 664    rst->score[1] = vmptr->score;
 665    rst->score[0] = max (scor, vmptr->score);
 666    rst->score[2] = rst->score[0];               /* initn */
 667
 668    if (ppst->param_u.fa.optflag) {
 669      if (rst->score[0] > ppst->param_u.fa.optcut)
 670        rst->score[2] = dmatch (aa0, n0, aa1, n1, *hoff=noff - vmptr->dp,
 671                              ppst->param_u.fa.optwid, ppst->pam2[ip],
 672                              ppst->gdelval,ppst->ggapval,f_str);
 673    }
 674 }
 675
 676 void do_work (const unsigned char *aa0, int n0,
 677               const unsigned char *aa1, int n1,
 678               int frame,
 679               struct pstruct *ppst, struct f_struct *f_str,
 680               int qr_flg, struct rstruct *rst)
 681 {
 682   int hoff, n10;
 683
 684   double lambda, H;
 685
 686   rst->score[0] = rst->score[1] = rst->score[2] = 0;
 687   rst->escore = 1.0;
 688   rst->segnum = rst->seglen = 1;
 689
 690   if (n1 < ppst->param_u.fa.ktup) return;
 691
 692 #ifdef TFASTA
 693   n10=aatran(aa1,f_str->aa1x,n1,frame);
 694   do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff);
 695 #else   /* FASTA */
 696   do_fasta (aa0, n0, aa1, n1, ppst, f_str, rst, &hoff);
 697 #endif
 698
 699 #ifndef TFASTA
 700   if((ppst->zsflag%10) == 6 &&
 701      (do_karlin(aa1, n1, ppst->pam2[0], ppst,f_str->aa0_f,
 702                 f_str->kar_p, &lambda, &H)>0)) {
 703     rst->comp = 1.0/lambda;
 704     rst->H = H;
 705   }
 706   else {rst->comp = rst->H = -1.0;}
 707 #else
 708   rst->comp = rst->H = -1.0;
 709 #endif
 710 }
 711
 712 void do_opt (const unsigned char *aa0, int n0,
 713              const unsigned char *aa1, int n1,
 714              int frame,
 715              struct pstruct *ppst,
 716              struct f_struct *f_str,
 717              struct rstruct *rst)
 718 {
 719   int optflag, tscore, hoff, n10;
 720
 721   optflag = ppst->param_u.fa.optflag;
 722   ppst->param_u.fa.optflag = 1;
 723
 724 #ifdef TFASTA
 725   n10=aatran(aa1,f_str->aa1x,n1,frame);
 726   do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff);
 727 #else   /* FASTA */
 728   do_fasta(aa0,n0,aa1,n1,ppst,f_str,rst, &hoff);
 729 #endif
 730   ppst->param_u.fa.optflag = optflag;
 731 }
 732
 733 #ifdef ALLOCN0
 734 void
 735 savemax (dptr, dpos, f_str)
 736   register struct dstruct *dptr;
 737   int  dpos;
 738   struct f_struct *f_str;
 739 {
 740    register struct savestr *vmptr;
 741    register int i;
 742
 743 #else
 744 void
 745 savemax (dptr, f_str)
 746   register struct dstruct *dptr;
 747   struct f_struct *f_str;
 748 {
 749    register int dpos;
 750    register struct savestr *vmptr;
 751    register int i;
 752
 753    dpos = (int) (dptr - f_str->diag);
 754
 755 #endif
 756
 757 /* check to see if this is the continuation of a run that is already saved */
 758
 759    if ((vmptr = dptr->dmax) != NULL && vmptr->dp == dpos &&
 760          vmptr->start == dptr->start)
 761    {
 762       vmptr->stop = dptr->stop;
 763       if ((i = dptr->score) <= vmptr->score)
 764          return;
 765       vmptr->score = i;
 766       if (vmptr != f_str->lowmax)
 767          return;
 768    }
 769    else
 770    {
 771       i = f_str->lowmax->score = dptr->score;
 772       f_str->lowmax->dp = dpos;
 773       f_str->lowmax->start = dptr->start;
 774       f_str->lowmax->stop = dptr->stop;
 775       dptr->dmax = f_str->lowmax;
 776    }
 777
 778    for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
 779       if (vmptr->score < i)
 780       {
 781          i = vmptr->score;
 782          f_str->lowmax = vmptr;
 783       }
 784    f_str->lowscor = i;
 785 }
 786
 787 int spam (const unsigned char *aa0, const unsigned char *aa1,
 788           struct savestr *dmax, int **pam2, int xdrop,
 789           int noff, int do_extend)
 790 {
 791    register int lpos, tot;
 792    register const unsigned char *aa0p, *aa1p;
 793
 794    int drop_thresh;
 795
 796    struct {
 797      int     start, stop, score;
 798    } curv, maxv;
 799
 800    aa1p = &aa1[lpos= dmax->start];      /* get the start of lib seq */
 801    aa0p = &aa0[lpos - dmax->dp + noff]; /* start of query */
 802 #ifdef DEBUG
 803    /* also add check in calling routine */
 804    if (aa0p < aa0) { return -99; }
 805 #endif
 806    curv.start = lpos;                   /* start index in lib seq */
 807
 808    tot = curv.score = maxv.score = 0;
 809
 810    for (; lpos <= dmax->stop; lpos++) {
 811      tot += pam2[*aa0p++][*aa1p++];
 812      if (tot > curv.score) {            /* update current score */
 813        curv.stop = lpos;
 814        curv.score = tot;
 815       }
 816       else if (tot < 0) {
 817         if (curv.score > maxv.score) {  /* save score, start, stop */
 818           maxv.start = curv.start;
 819           maxv.stop = curv.stop;
 820           maxv.score = curv.score;
 821         }
 822         tot = curv.score = 0;           /* reset running score */
 823         curv.start = lpos+1;            /* reset start */
 824         if(lpos >= dmax->stop) break;   /* if the zero is beyond stop, quit */
 825       }
 826    }
 827
 828    if (curv.score > maxv.score) {
 829      maxv.start = curv.start;
 830      maxv.stop = curv.stop;
 831      maxv.score = curv.score;
 832    }
 833
 834 #ifndef NOSPAM_EXT
 835
 836    /* now check to see if the score gets better by extending */
 837    if (do_extend && maxv.score > xdrop) {
 838
 839      if (maxv.stop == dmax->stop) {
 840        tot = maxv.score;
 841        drop_thresh = maxv.score - xdrop;
 842        aa1p = &aa1[lpos= dmax->stop];
 843        aa0p = &aa0[lpos - dmax->dp + noff];
 844        while (tot > drop_thresh ) {
 845          ++lpos;
 846          tot += pam2[*(++aa0p)][*(++aa1p)];
 847          if (tot > maxv.score) {
 848            maxv.start = lpos;
 849            maxv.score = tot;
 850            drop_thresh = tot - xdrop;
 851          }
 852        }
 853      }
 854
 855      /* scan backwards now */
 856
 857      if (maxv.start == dmax->start) {
 858        tot = maxv.score;
 859        drop_thresh = maxv.score - xdrop;
 860        aa1p = &aa1[lpos= dmax->start];
 861        aa0p = &aa0[lpos - dmax->dp + noff];
 862        while (tot > drop_thresh) {
 863          --lpos;
 864          tot += pam2[*(--aa0p)][*(--aa1p)];
 865          if (tot > maxv.score) {
 866            maxv.start = lpos;
 867            maxv.score = tot;
 868            drop_thresh = tot - xdrop;
 869          }
 870        }
 871      }
 872    }
 873 #endif
 874
 875 /*      if (maxv.start != dmax->start || maxv.stop != dmax->stop)
 876                 printf(" new region: %3d %3d %3d %3d\n",maxv.start,
 877                         dmax->start,maxv.stop,dmax->stop);
 878 */
 879    dmax->start = maxv.start;
 880    dmax->stop = maxv.stop;
 881
 882    return maxv.score;
 883 }
 884
 885 int sconn (struct savestr **v, int n, int cgap, int pgap, int noff)
 886 {
 887    int     i, si;
 888    struct slink
 889    {
 890       int     score;
 891       struct savestr *vp;
 892       struct slink *next;
 893    }      *start, *sl, *sj, *so, sarr[MAXSAV];
 894    int     lstart, tstart, plstop, ptstop;
 895
 896 /*      sort the score left to right in lib pos */
 897
 898    kpsort (v, n);
 899
 900    start = NULL;
 901
 902 /*      for the remaining runs, see if they fit */
 903
 904    for (i = 0, si = 0; i < n; i++)
 905    {
 906
 907 /*      if the score is less than the gap penalty, it never helps */
 908       if (v[i]->score < cgap)
 909          continue;
 910       lstart = v[i]->start;
 911       tstart = lstart - v[i]->dp + noff;
 912
 913 /*      put the run in the group */
 914       sarr[si].vp = v[i];
 915       sarr[si].score = v[i]->score;
 916       sarr[si].next = NULL;
 917
 918 /*      if it fits, then increase the score */
 919       for (sl = start; sl != NULL; sl = sl->next)
 920       {
 921          plstop = sl->vp->stop;
 922          ptstop = plstop - sl->vp->dp + noff;
 923          if (plstop < lstart && ptstop < tstart)
 924          {
 925             sarr[si].score = sl->score + v[i]->score + pgap;
 926             break;
 927          }
 928       }
 929
 930 /*      now recalculate where the score fits */
 931       if (start == NULL)
 932          start = &sarr[si];
 933       else
 934          for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
 935             if (sarr[si].score > sj->score) {
 936                sarr[si].next = sj;
 937                if (so != NULL) so->next = &sarr[si];
 938                else  start = &sarr[si];
 939                break;
 940             }
 941             so = sj;
 942          }
 943       si++;
 944    }
 945
 946    if (start != NULL)
 947       return (start->score);
 948    else
 949       return (0);
 950 }
 951
 952 void
 953 kssort (v, n)
 954 struct savestr *v[];
 955 int     n;
 956 {
 957    int     gap, i, j;
 958    struct savestr *tmp;
 959
 960    for (gap = n / 2; gap > 0; gap /= 2)
 961       for (i = gap; i < n; i++)
 962          for (j = i - gap; j >= 0; j -= gap)
 963          {
 964             if (v[j]->score >= v[j + gap]->score)
 965                break;
 966             tmp = v[j];
 967             v[j] = v[j + gap];
 968             v[j + gap] = tmp;
 969          }
 970 }
 971
 972 void
 973 kpsort (v, n)
 974 struct savestr *v[];
 975 int     n;
 976 {
 977    int     gap, i, j;
 978    struct savestr *tmp;
 979
 980    for (gap = n / 2; gap > 0; gap /= 2)
 981       for (i = gap; i < n; i++)
 982          for (j = i - gap; j >= 0; j -= gap)
 983          {
 984             if (v[j]->start <= v[j + gap]->start)
 985                break;
 986             tmp = v[j];
 987             v[j] = v[j + gap];
 988             v[j + gap] = tmp;
 989          }
 990 }
 991
 992 static int dmatch (const unsigned char *aa0, int n0,
 993                    const unsigned char *aa1, int n1,
 994                    int hoff, int window,
 995                    int **pam2, int gdelval, int ggapval,
 996                    struct f_struct *f_str)
 997 {
 998    int low, up;
 999
1000    window = min (n1, window);
1001    /* hoff is the offset found from aa1 to seq 2 by hmatch */
1002
1003    low = -window/2-hoff;
1004    up = low+window;
1005
1006    return FLOCAL_ALIGN(aa0-1,aa1-1,n0,n1, low, up,
1007                       pam2,
1008 #ifdef OLD_FASTA_GAP
1009                        -(gdelval-ggapval),
1010 #else
1011                        -gdelval,
1012 #endif
1013                        -ggapval,window,f_str);
1014  }
1015
1016
1017 /* A PACKAGE FOR LOCALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
1018
1019    To invoke, call LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,MW).
1020    The parameters are explained as follows:
1021         A, B : two sequences to be aligned
1022         M : the length of sequence A
1023         N : the length of sequence B
1024         L : lower bound of the band
1025         U : upper bound of the band
1026         W : scoring table for matches and mismatches
1027         G : gap-opening penalty
1028         H : gap-extension penalty
1029         MW  : maximum window size
1030 */
1031
1032 #include <stdio.h>
1033
1034 #define MININT -9999999
1035
1036 static int
1037 FLOCAL_ALIGN(const unsigned char *A, const unsigned char *B,
1038              int M, int N, int low, int up,
1039              int **W, int G,int H, int MW,
1040              struct f_struct *f_str)
1041 {
1042   int band;
1043   register struct bdstr *bssp;
1044   int i, j, si, ei;
1045   int c, d, e, m;
1046   int leftd, rightd;
1047   int best_score;
1048   int *wa, curd;
1049   int ib;
1050
1051   bssp = f_str->bss;
1052
1053   m = G+H;
1054   low = max(-M, low);
1055   up = min(N, up);
1056
1057   if (N <= 0) return 0;
1058
1059   if (M <= 0) return 0;
1060
1061   band = up-low+1;
1062   if (band < 1) {
1063     fprintf(stderr,"low > up is unacceptable!: M: %d N: %d l/u: %d/%d\n",
1064             M, N, low, up);
1065     return 0;
1066   }
1067
1068   if (low > 0) leftd = 1;
1069   else if (up < 0) leftd = band;
1070   else leftd = 1-low;
1071   rightd = band;
1072   si = max(0,-up);      /* start index -1 */
1073   ei = min(M,N-low);    /* end index */
1074   bssp[leftd].CC = 0;
1075   for (j = leftd+1; j <= rightd; j++) {
1076     bssp[j].CC = 0;
1077     bssp[j].DD = -G;
1078   }
1079
1080   bssp[rightd+1].CC = MININT;
1081   bssp[rightd+1].DD = MININT;
1082
1083   best_score = 0;
1084   bssp[leftd-1].CC = MININT;
1085   bssp[leftd].DD = -G;
1086
1087   for (i = si+1; i <= ei; i++) {
1088     if (i > N-up) rightd--;
1089     if (leftd > 1) leftd--;
1090     wa = W[A[i]];
1091     if ((c = bssp[leftd+1].CC-m) > (d = bssp[leftd+1].DD-H)) d = c;
1092     if ((ib = leftd+low-1+i ) > 0) c = bssp[leftd].CC+wa[B[ib]];
1093
1094     if (d > c) c = d;
1095     if (c < 0) c = 0;
1096     e = c-G;
1097     bssp[leftd].DD = d;
1098     bssp[leftd].CC = c;
1099     if (c > best_score) best_score = c;
1100
1101     for (curd=leftd+1; curd <= rightd; curd++) {
1102       if ((c = c-m) > (e = e-H)) e = c;
1103       if ((c = bssp[curd+1].CC-m) > (d = bssp[curd+1].DD-H)) d = c;
1104       c = bssp[curd].CC + wa[B[curd+low-1+i]];
1105       if (e > c) c = e;
1106       if (d > c) c = d;
1107       if (c < 0) c = 0;
1108       bssp[curd].CC = c;
1109       bssp[curd].DD = d;
1110       if (c > best_score) best_score = c;
1111     }
1112   }
1113
1114   return best_score;
1115 }
1116
1117 /* ckalloc - allocate space; check for success */
1118 char *ckalloc(size_t amount)
1119 {
1120   char *p;
1121
1122   if ((p = malloc( (unsigned) amount)) == NULL)
1123     w_abort("Ran out of memory.","");
1124   return(p);
1125 }
1126
1127 /* calculate the 100% identical score */
1128 int
1129 shscore(const unsigned char *aa0, int n0, int **pam2)
1130 {
1131   int i, sum;
1132   for (i=0,sum=0; i<n0; i++)
1133     sum += pam2[aa0[i]][aa0[i]];
1134   return sum;
1135 }
1136
1137 int  sw_walign (const unsigned char *aa0, int n0,
1138                 const unsigned char *aa1, int n1,
1139                 struct pstruct *ppst,
1140                 struct f_struct *f_str,
1141                 struct a_res_str *a_res
1142                 )
1143 {
1144    register const unsigned char *aa0p, *aa1p;
1145    register int *pwaa;
1146    register int i, j;
1147    register struct swstr *ssj;
1148    struct swstr *ss;
1149    int *waa;
1150    int e, f, h, p;
1151    int     q, r, m;
1152    int     score;
1153    int cost, I, J, K, L;
1154
1155    ss = f_str->ss;
1156    waa = f_str->waa1;
1157
1158 #ifdef OLD_FASTA_GAP
1159    q = -(ppst->gdelval - ppst->ggapval);
1160 #else
1161    q = -ppst->gdelval;
1162 #endif
1163    r = -ppst->ggapval;
1164    m = q + r;
1165
1166    /* initialize 0th row */
1167    for (ssj=ss; ssj<ss+n0; ssj++) {
1168      ssj->H = 0;
1169      ssj->E = -q;
1170    }
1171
1172    score = I = J = 0;
1173    aa1p = aa1;
1174    i = 0;
1175    while (*aa1p) {
1176      h = p = 0;
1177      f = -q;
1178      pwaa = waa + (*aa1p++ * n0);
1179      for (ssj = ss, aa0p = aa0; ssj < ss+n0; ssj++) {
1180        if ((h =   h     - m) > (f =   f     - r)) f = h;
1181        if ((h = ssj->H - m) > (e = ssj->E - r)) e = h;
1182        h = p + *pwaa++;
1183        if (h < 0 ) h = 0;
1184        if (h < f ) h = f;
1185        if (h < e ) h = e;
1186        p = ssj->H;
1187        ssj->H = h;
1188        ssj->E = e;
1189        if (h > score) {
1190          score = h;
1191          I = i;
1192          J = (int)(ssj-ss);
1193        }
1194      }
1195      i++;
1196    }                            /* done with forward pass */
1197    if (score <= 0) return 0;
1198
1199   /* to get the start point, go backwards */
1200
1201   cost = K = L = 0;
1202   for (ssj=ss+J; ssj>=ss; ssj--) ssj->H= ssj->E= -1;
1203
1204   for (i=I; i>=0; i--) {
1205     h = f = -1;
1206     p = (i == I) ? 0 : -1;
1207     ssj = ss+J;  /* bug in compiler */
1208     for (aa0p = &aa0[J]; ssj>=ss; ssj--,aa0p--) {
1209       f = max (f,h-q)-r;
1210       ssj->E=max(ssj->E,ssj->H-q)-r;
1211       h = max(max(ssj->E,f),p+ppst->pam2[0][*aa0p][aa1[i]]);
1212       p = ssj->H;
1213       ssj->H=h;
1214       if (h > cost) {
1215         cost = h;
1216         K = i;
1217         L = (int)(ssj-ss);
1218         if (cost >= score) goto found;
1219       }
1220     }
1221   }
1222
1223 found:
1224
1225   /*  printf(" %d: L: %3d-%3d/%3d; K: %3d-%3d/%3d\n",score,L,J,n0,K,I,n1); */
1226
1227 /* in the f_str version, the *res array is already allocated at 4*n0/3 */
1228
1229   a_res->max0 = J+1; a_res->min0 = L; a_res->max1 = I+1; a_res->min1 = K;
1230
1231   /* the seq array arguments in this call have been reversed to allow
1232      assymetric scoring matrices - this affects the score decoding,
1233      and allocation of the score row matrix */
1234   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);
1235
1236   /*  DISPLAY(&aa1[K-1],&aa0[L-1],I-K+1,J-L+1,res,L,K,ppst->sq); */
1237
1238    return score;
1239 }
1240
1241 static int CHECK_SCORE(const unsigned char *A, const unsigned char *B,
1242                        int M, int N,
1243                        int *S, int **W, int G, int H, int *nres);
1244
1245 #define gap(k)  ((k) <= 0 ? 0 : g+h*(k))        /* k-symbol indel cost */
1246
1247 /* static int *sapp; */                         /* Current script append ptr */
1248 /* static int  last; */                         /* Last script op appended */
1249
1250                                                 /* Append "Delete k" op */
1251 #define DEL(k)                          \
1252 { if (*last < 0)                                \
1253     *last = (*sapp)[-1] -= (k);         \
1254   else  {                               \
1255     *last = (*sapp)[0] = -(k);          \
1256     (*sapp)++;                          \
1257   }                                     \
1258 }
1259                                                 /* Append "Insert k" op */
1260 #define INS(k)                          \
1261 { if (*last > 0)                        \
1262      *last = (*sapp)[-1] += (k);        \
1263   else  {                               \
1264     *last = (*sapp)[0] = (k);           \
1265     (*sapp)++;                          \
1266   }                                     \
1267 }
1268
1269 #define REP { *last = (*sapp)[0] = 0; (*sapp)++;} /* Append "Replace" op */
1270
1271 /* align(A,B,M,N,tb,te) returns the cost of an optimum conversion between
1272    A[1..M] and B[1..N] that begins(ends) with a delete if tb(te) is zero
1273    and appends such a conversion to the current script.                   */
1274
1275 static int
1276 nw_align(const unsigned char *A, const unsigned char *B,
1277          int M, int N,
1278          int tb, int te, int **w, int g, int h,
1279          struct f_struct *f_str,
1280          int **sapp, int *last)
1281 {
1282   int   midi, midj, type;       /* Midpoint, type, and cost */
1283   int midc;
1284
1285   register int   i, j;
1286   register int c, e, d, s;
1287   int m, t, *wa;
1288   struct swstr *f_ss, *r_ss;
1289
1290   m = g + h;
1291
1292   f_ss = f_str->f_ss;
1293   r_ss = f_str->r_ss;
1294
1295 /* Boundary cases: M <= 1 or N == 0 */
1296
1297   if (N <= 0) {
1298     if (M > 0) {DEL(M)}
1299      return -gap(M);
1300   }
1301
1302   if (M <= 1) {
1303     if (M <= 0) {
1304       INS(N);
1305       return -gap(N);
1306     }
1307     if (tb < te) tb = te;
1308     midc = (tb-h) - gap(N);
1309     midj = 0;
1310     wa = w[A[1]];       /* in the original version of this code, A[]
1311                            is the second sequence */
1312     for (j = 1; j <= N; j++) {
1313       c = -gap(j-1) + wa[B[j]] - gap(N-j);
1314       if (c > midc) {
1315         midc = c;
1316         midj = j;
1317       }
1318     }
1319     if (midj == 0) { DEL(1) INS(N) }
1320     else {
1321       if (midj > 1) { INS(midj-1) }
1322       REP
1323       if (midj < N) { INS(N-midj) }
1324     }
1325     return midc;
1326   }
1327
1328 /* Divide: Find optimum midpoint (midi,midj) of cost midc */
1329
1330   midi = M/2;           /* Forward phase:                          */
1331   f_ss[0].H = 0;                /*   Compute H(M/2,k) & E(M/2,k) for all k */
1332   t = -g;
1333   for (j = 1; j <= N; j++) {
1334     f_ss[j].H = t = t-h;
1335     f_ss[j].E = t-g;
1336   }
1337   t = tb;
1338   for (i = 1; i <= midi; i++) {
1339     s = f_ss[0].H;
1340     f_ss[0].H = c = t = t-h;
1341     e = t-g;
1342     wa = w[A[i]];
1343     for (j = 1; j <= N; j++) {
1344       if ((c =   c   - m) > (e =   e   - h)) e = c;
1345       if ((c = f_ss[j].H - m) > (d = f_ss[j].E - h)) d = c;
1346       c = s + wa[B[j]];
1347       if (e > c) c = e;
1348       if (d > c) c = d;
1349       s = f_ss[j].H;
1350       f_ss[j].H = c;
1351       f_ss[j].E = d;
1352     }
1353   }
1354   f_ss[0].E = f_ss[0].H;
1355
1356   r_ss[N].H = 0;                /* Reverse phase:                          */
1357   t = -g;                       /*   Compute R(M/2,k) & S(M/2,k) for all k */
1358   for (j = N-1; j >= 0; j--)
1359     { r_ss[j].H = t = t-h;
1360       r_ss[j].E = t-g;
1361     }
1362   t = te;
1363   for (i = M-1; i >= midi; i--)
1364     { s = r_ss[N].H;
1365       r_ss[N].H = c = t = t-h;
1366       e = t-g;
1367       wa = w[A[i+1]];
1368       for (j = N-1; j >= 0; j--)
1369         { if ((c =   c   - m) > (e =   e   - h)) e = c;
1370           if ((c = r_ss[j].H - m) > (d = r_ss[j].E - h)) d = c;
1371           c = s + wa[B[j+1]];
1372           if (e > c) c = e;
1373           if (d > c) c = d;
1374           s = r_ss[j].H;
1375           r_ss[j].H = c;
1376           r_ss[j].E = d;
1377         }
1378     }
1379   r_ss[N].E = r_ss[N].H;
1380
1381   midc = f_ss[0].H+r_ss[0].H;           /* Find optimal midpoint */
1382   midj = 0;
1383   type = 1;
1384   for (j = 0; j <= N; j++)
1385     if ((c = f_ss[j].H + r_ss[j].H) >= midc)
1386       if (c > midc || (f_ss[j].H != f_ss[j].E && r_ss[j].H == r_ss[j].E))
1387         { midc = c;
1388           midj = j;
1389         }
1390   for (j = N; j >= 0; j--)
1391     if ((c = f_ss[j].E + r_ss[j].E + g) > midc)
1392       { midc = c;
1393         midj = j;
1394         type = 2;
1395       }
1396
1397
1398 /* Conquer: recursively around midpoint */
1399
1400   if (type == 1) {
1401     nw_align(A,B,midi,midj,tb,-g,w,g,h,f_str, sapp, last);
1402     nw_align(A+midi,B+midj,M-midi,N-midj,-g,te,w,g,h,f_str,sapp, last);
1403   }
1404   else {
1405     nw_align(A,B,midi-1,midj,tb,0,w,g,h,f_str, sapp, last);
1406     DEL(2);
1407     nw_align(A+midi+1,B+midj,M-midi-1,N-midj,0,te,w,g,h,f_str, sapp, last);
1408   }
1409   return midc;
1410 }
1411
1412 /* Interface and top level of comparator */
1413
1414 static int
1415 ALIGN(const unsigned char *A, const unsigned char *B,
1416       int M, int N,
1417       int **W, int G, int H, int *S, int *nS,
1418       struct f_struct *f_str)
1419 {
1420   int c, ck;
1421   struct swstr *f_ss, *r_ss;
1422   int *sapp, last;
1423
1424   sapp = S;
1425   last = 0;
1426
1427    if ((f_ss = (struct swstr *) calloc (N+2, sizeof (struct swstr)))
1428        == NULL) {
1429      fprintf (stderr, "cannot allocate f_ss array %3d\n", N+2);
1430      exit (1);
1431    }
1432    f_ss++;
1433    f_str->f_ss = f_ss;
1434
1435    if ((r_ss = (struct swstr *) calloc (N+2, sizeof (struct swstr)))
1436        == NULL) {
1437      fprintf (stderr, "cannot allocate r_ss array %3d\n", N+2);
1438      exit (1);
1439    }
1440    r_ss++;
1441    f_str->r_ss = r_ss;
1442
1443   c = nw_align(A,B,M,N,-G,-G,W,G,H,f_str,&sapp, &last); /* OK, do it */
1444
1445   ck = CHECK_SCORE(A,B,M,N,S,W,G,H,nS);
1446   if (c != ck) fprintf(stderr,"Check_score error %d != %d\n",c,ck);
1447
1448   f_ss--; r_ss--;
1449   free(r_ss); free(f_ss);
1450
1451   return c;
1452 }
1453
1454 /* Alignment display routine */
1455
1456 static char ALINE[51], BLINE[51], CLINE[51];
1457
1458 void DISPLAY(unsigned char *A, unsigned char *B, int M, int N,
1459             int *S, int AP, int BP, char *sq)
1460 { register char *a, *b, *c;
1461   register int   i,  j, op;
1462            int   lines, ap, bp;
1463
1464   i = j = op = lines = 0;
1465   ap = AP;
1466   bp = BP;
1467   a = ALINE;
1468   b = BLINE;
1469   c = CLINE;
1470   while (i < M || j < N)
1471     { if (op == 0 && *S == 0)
1472         { op = *S++;
1473           *a = sq[A[++i]];
1474           *b = sq[B[++j]];
1475           *c++ = (*a++ == *b++) ? '|' : ' ';
1476         }
1477       else
1478         { if (op == 0)
1479             op = *S++;
1480           if (op > 0)
1481             { *a++ = ' ';
1482               *b++ = sq[B[++j]];
1483               op--;
1484             }
1485           else
1486             { *a++ = sq[A[++i]];
1487               *b++ = ' ';
1488               op++;
1489             }
1490           *c++ = '-';
1491         }
1492       if (a >= ALINE+50 || (i >= M && j >= N))
1493         { *a = *b = *c = '\0';
1494           printf("\n%5d ",50*lines++);
1495           for (b = ALINE+10; b <= a; b += 10)
1496             printf("    .    :");
1497           if (b <= a+5)
1498             printf("    .");
1499           printf("\n%5d %s\n      %s\n%5d %s\n",ap,ALINE,CLINE,bp,BLINE);
1500           ap = AP + i;
1501           bp = BP + j;
1502           a = ALINE;
1503           b = BLINE;
1504           c = CLINE;
1505         }
1506     }
1507 }
1508
1509 /* CHECK_SCORE - return the score of the alignment stored in S */
1510
1511 static int CHECK_SCORE(const unsigned char *A, const unsigned char *B,
1512                        int M, int N, int *S, int **w, int g, int h,
1513                        int *nres)
1514 {
1515   register int   i,  j, op, nc;
1516   int score;
1517
1518   score = i = j = op = nc = 0;
1519   while (i < M || j < N) {
1520         op = *S++;
1521         if (op == 0) {
1522           score = w[A[++i]][B[++j]] + score;
1523           nc++;
1524   /*      fprintf(stderr,"=%4d %4d %4d %4d\n",i,j,w[A[i]][B[i]],score); */
1525         }
1526         else if (op > 0) {
1527           score = score - (g+op*h);
1528   /*      fprintf(stderr,">%4d %4d %4d %4d\n",i,j,-(g+op*h),score); */
1529           j += op;
1530           nc += op;
1531         } else {
1532           score = score - (g-op*h);
1533   /*      fprintf(stderr,"<%4d %4d %4d %4d\n",i,j,-(g-op*h),score); */
1534           i -= op;
1535           nc -= op;
1536         }
1537   }
1538   *nres = nc;
1539   return score;
1540 }
1541
1542
1543 static int
1544 BCHECK_SCORE(const unsigned char *A, const unsigned char *B,
1545                        int M, int N, int *S, int **w, int g, int h,
1546                        int *nres)
1547 {
1548   register int   i,  j, op, nc;
1549   int *Ssave;
1550   int score;
1551
1552   score = i = j = op = nc = 0;
1553   Ssave = S;
1554   while (i < M || j < N) {
1555         op = *S++;
1556         if (op == 0) {
1557           score = w[A[++i]][B[++j]] + score;
1558           nc++;
1559 /*        fprintf(stderr,"op0 %4d %4d %4d %4d\n",i,j,w[A[i]][B[i]],score); */
1560         }
1561         else if (op > 0) {
1562           score = score - (g+op*h);
1563 /*        fprintf(stderr,"op> %4d %4d %4d %4d %4d\n",i,j,op,-(g+op*h),score); */
1564           j += op;
1565           nc += op;
1566         } else {
1567           score = score - (g-op*h);
1568 /*        fprintf(stderr,"op< %4d %4d %4d %4d %4d\n",i,j,op,-(g-op*h),score); */
1569           i -= op;
1570           nc -= op;
1571         }
1572   }
1573   *nres = nc;
1574   return score;
1575 }
1576
1577
1578 /* A PACKAGE FOR LOCALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
1579
1580    To invoke, call LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,S,dflag,&SI,&SJ,&EI,&EJ,MW).
1581    The parameters are explained as follows:
1582         A, B : two sequences to be aligned
1583         M : the length of sequence A
1584         N : the length of sequence B
1585         L : lower bound of the band
1586         U : upper bound of the band
1587         W : scoring table for matches and mismatches
1588         G : gap-opening penalty
1589         H : gap-extension penalty
1590         dflag : 0 - no display or backward pass
1591         *SI : starting position of sequence A in the optimal local alignment
1592         *SJ : starting position of sequence B in the optimal local alignment
1593         *EI : ending position of sequence A in the optimal local alignment
1594         *EJ : ending position of sequence B in the optimal local alignment
1595         MW  : maximum window size
1596 */
1597
1598 int bd_walign (const unsigned char *aa0, int n0,
1599                  const unsigned char *aa1, int n1,
1600                  struct pstruct *ppst,
1601                  struct f_struct *f_str, int hoff,
1602                  struct a_res_str *a_res)
1603 {
1604    int low, up, score;
1605    int min0, min1, max0, max1;
1606    int window;
1607
1608    window = min (n1, ppst->param_u.fa.optwid);
1609    /* hoff is the offset found from aa1 to seq 2 by hmatch */
1610
1611    low = -window/2-hoff;
1612    up = low+window;
1613
1614    score=LOCAL_ALIGN(aa0-1,aa1-1,n0,n1, low, up,
1615                     ppst->pam2[0],
1616 #ifdef OLD_FASTA_GAP
1617                      -(ppst->gdelval-ppst->ggapval),
1618 #else
1619                      -ppst->gdelval,
1620 #endif
1621                      -ppst->ggapval,
1622                     &min0,&min1,&max0,&max1,ppst->param_u.fa.optwid,f_str);
1623
1624    if (score <=0) {
1625      fprintf(stderr,"n0/n1: %d/%d hoff: %d window: %d\n",
1626              n0, n1, hoff, window);
1627      return 0;
1628    }
1629
1630 /*
1631   fprintf(stderr," ALIGN: start0: %d start1: %d stop0: %d stop1: %d, bot: %d top: %d, win: %d MX %d\n",
1632   min0-1,min1-1,max0-min0+1,max1-min1+1,low-(min1-min0),up-(min1-min0),
1633   ppst->param_u.fa.optwid,n0);
1634 */
1635
1636    a_res->min0 = min0-1; a_res->min1 = min1-1;
1637    a_res->max0 = max0; a_res->max1 = max1;
1638
1639    B_ALIGN(aa0-1+min0-1,aa1-1+min1-1,max0-min0+1,max1-min1+1,
1640            low-(min1-min0),up-(min1-min0),
1641            ppst->pam2[0],
1642 #ifdef OLD_FASTA_GAP
1643            -(ppst->gdelval-ppst->ggapval),
1644 #else
1645            -ppst->gdelval,
1646 #endif
1647            -ppst->ggapval,
1648            a_res->res,&a_res->nres,ppst->param_u.fa.optwid,n0,f_str->bss);
1649
1650    return score;
1651 }
1652
1653 static int
1654 LOCAL_ALIGN(const unsigned char *A, const unsigned char *B,
1655             int M, int N,
1656             int low, int up, int **W, int G,int H,
1657             int *psi, int *psj, int *pei, int *pej, int MW,
1658             struct f_struct *f_str)
1659 {
1660   int band;
1661   register struct bdstr *bssp;
1662   int i, j, si, ei;
1663   int c, d, e, t, m;
1664   int leftd, rightd;
1665   int best_score, starti, startj, endi, endj;
1666   int *wa, curd;
1667   int ib;
1668   char flag;
1669
1670   bssp = f_str->bss;
1671
1672   m = G+H;
1673   low = max(-M, low);
1674   up = min(N, up);
1675
1676   if (N <= 0) {
1677     *psi = *psj = *pei = *pej;
1678     return 0;
1679   }
1680   if (M <= 0) {
1681     *psi = *psj = *pei = *pej;
1682     return 0;
1683   }
1684   band = up-low+1;
1685   if (band < 1) {
1686     fprintf(stderr,"low > up is unacceptable!: M: %d N: %d l/u: %d/%d\n",
1687             M, N, low, up);
1688     return -1;
1689   }
1690
1691   j = (MW + 2 + 2) * sizeof(struct bdstr);
1692
1693   /* already done by init_work();
1694   if (f_str->bss==NULL) f_str->bss = (struct bdstr *) ckalloc(j);
1695   */
1696
1697   if (low > 0) leftd = 1;
1698   else if (up < 0) leftd = band;
1699   else leftd = 1-low;
1700   rightd = band;
1701   si = max(0,-up);
1702   ei = min(M,N-low);
1703   bssp[leftd].CC = 0;
1704   for (j = leftd+1; j <= rightd; j++) {
1705     bssp[j].CC = 0;
1706     bssp[j].DD = -G;
1707   }
1708   bssp[rightd+1].CC = MININT;
1709   bssp[rightd+1].DD = MININT;
1710   best_score = 0;
1711   endi = si;
1712   endj = si+low;
1713   bssp[leftd-1].CC = MININT;
1714   bssp[leftd].DD = -G;
1715   for (i = si+1; i <= ei; i++) {
1716     if (i > N-up) rightd--;
1717     if (leftd > 1) leftd--;
1718     wa = W[A[i]];
1719     if ((c = bssp[leftd+1].CC-m) > (d = bssp[leftd+1].DD-H)) d = c;
1720     if ((ib = leftd+low-1+i ) > 0) c = bssp[leftd].CC+wa[B[ib]];
1721 /*
1722     if (ib > N) fprintf(stderr,"B[%d] out of range %d\n",ib,N);
1723 */
1724     if (d > c) c = d;
1725     if (c < 0) c = 0;
1726     e = c-G;
1727     bssp[leftd].DD = d;
1728     bssp[leftd].CC = c;
1729     if (c > best_score) {
1730       best_score = c;
1731       endi = i;
1732       endj = ib;
1733     }
1734     for (curd=leftd+1; curd <= rightd; curd++) {
1735       if ((c = c-m) > (e = e-H)) e = c;
1736       if ((c = bssp[curd+1].CC-m) > (d = bssp[curd+1].DD-H)) d = c;
1737 /*
1738       if ((ib=curd+low-1+i) <= 0 || ib > N)
1739         fprintf(stderr,"B[%d]:%d\n",ib,B[ib]);
1740 */
1741       c = bssp[curd].CC + wa[B[curd+low-1+i]];
1742       if (e > c) c = e;
1743       if (d > c) c = d;
1744       if (c < 0) c = 0;
1745       bssp[curd].CC = c;
1746       bssp[curd].DD = d;
1747       if (c > best_score) {
1748         best_score = c;
1749         endi = i;
1750         endj = curd+low-1+i;
1751       }
1752     }
1753   }
1754
1755   leftd = max(1,-endi-low+1);
1756   rightd = band-(up-(endj-endi));
1757   bssp[rightd].CC = 0;
1758   t = -G;
1759   for (j = rightd-1; j >= leftd; j--) {
1760     bssp[j].CC = t = t-H;
1761     bssp[j].DD = t-G;
1762   }
1763   for (j = rightd+1; j <= band; ++j) bssp[j].CC = MININT;
1764   bssp[leftd-1].CC = bssp[leftd-1].DD = MININT;
1765   bssp[rightd].DD = -G;
1766   flag = 0;
1767   for (i = endi; i >= 1; i--) {
1768     if (i+low <= 0) leftd++;
1769     if (rightd < band) rightd++;
1770     wa = W[A[i]];
1771     if ((c = bssp[rightd-1].CC-m) > (d = bssp[rightd-1].DD-H)) d = c;
1772     if ((ib = rightd+low-1+i) <= N) c = bssp[rightd].CC+wa[B[ib]];
1773
1774 /*
1775     if (ib <= 0) fprintf(stderr,"rB[%d] <1\n",ib);
1776 */
1777     if (d > c) c = d;
1778     e = c-G;
1779     bssp[rightd].DD = d;
1780     bssp[rightd].CC = c;
1781     if (c == best_score) {
1782       starti = i;
1783       startj = ib;
1784       flag = 1;
1785       break;
1786     }
1787     for (curd=rightd-1; curd >= leftd; curd--) {
1788       if ((c = c-m) > (e = e-H)) e = c;
1789       if ((c = bssp[curd-1].CC-m) > (d = bssp[curd-1].DD-H)) d = c;
1790
1791 /*
1792       if ((ib=curd+low-1+i) <= 0 || ib > N)
1793         fprintf(stderr,"i: %d, B[%d]:%d\n",i,ib,B[ib]);
1794 */
1795       c = bssp[curd].CC + wa[B[curd+low-1+i]];
1796       if (e > c) c = e;
1797       if (d > c) c = d;
1798       bssp[curd].CC = c;
1799       bssp[curd].DD = d;
1800       if (c == best_score) {
1801         starti = i;
1802         startj = curd+low-1+i;
1803         flag = 1;
1804         break;
1805       }
1806     }
1807     if (flag == 1) break;
1808   }
1809
1810   if (starti < 0 || starti > M || startj < 0 || startj > N) {
1811     printf("starti=%d, startj=%d\n",starti,startj);
1812     *psi = *psj = *pei = *pej;
1813     exit(1);
1814   }
1815   *psi = starti;
1816   *psj = startj;
1817   *pei = endi;
1818   *pej = endj;
1819   return best_score;
1820 }
1821
1822 /* A PACKAGE FOR GLOBALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
1823
1824    To invoke, call B_ALIGN(A,B,M,N,L,U,W,G,H,S,MW,MX).
1825    The parameters are explained as follows:
1826         A, B : two sequences to be aligned
1827         M : the length of sequence A
1828         N : the length of sequence B
1829         L : lower bound of the band
1830         U : upper bound of the band
1831         W : scoring table for matches and mismatches
1832         G : gap-opening penalty
1833         H : gap-extension penalty
1834         S : script for DISPLAY routine
1835         MW : maximum window size
1836         MX : maximum length sequence M to be aligned
1837 */
1838
1839 static int IP;
1840 static int *MP[3];              /* save crossing points */
1841 static int *FP;                 /* forward dividing points */
1842 static char *MT[3];             /* 0: rep, 1: del, 2: ins */
1843 static char *FT;
1844
1845 /* bg_align(A,B,M,N,up,low,tb,te) returns the cost of an optimum conversion between
1846   A[1..M] and B[1..N] and appends such a conversion to the current script.
1847   tb(te)= 1  no gap-open penalty if the conversion begins(ends) with a delete.
1848   tb(te)= 2  no gap-open penalty if the conversion begins(ends) with an insert.
1849 */
1850 static int
1851 bg_align(const unsigned char *A, const unsigned char *B,
1852          int M, int N,
1853          int low, int up, int tb, int te,
1854          int **w, int g, int h,
1855          struct bdstr *bss, int **sapp, int *last)
1856 {
1857   int rmid, k, l, r, v, kt;
1858   int t1, t2, t3;
1859
1860   {
1861   int band, midd;
1862   int leftd, rightd;    /* for CC, DD, CP and DP */
1863   register int curd;    /* current index for CC, DD CP and DP */
1864   register int i, j;
1865   register int c, d, e;
1866   int t, fr, *wa, ib, m;
1867
1868   /* Boundary cases: M <= 0 , N <= 0, or up-low <= 0 */
1869   if (N <= 0) {
1870     if (M > 0) { DEL(M) }
1871     return 0;
1872   }
1873   if (M <= 0) {
1874     INS(N)
1875     return 0;
1876   }
1877   if ((band = up-low+1) <= 1) {
1878     for (i = 1; i <= M; i++) { REP }
1879     return 0;
1880   }
1881
1882   /* Divide: Find all crossing points */
1883
1884   /* Initialization */
1885   m = g + h;
1886
1887   midd = band/2 + 1;
1888   rmid = low + midd - 1;
1889   leftd = 1-low;
1890   rightd = up-low+1;
1891   if (leftd < midd) {
1892     fr = -1;
1893     for (j = 0; j < midd; j++)
1894       bss[j].CP = bss[j].DP = -1;
1895     for (j = midd; j <= rightd; j++) {
1896       bss[j].CP = bss[j].DP = 0;
1897     }
1898     MP[0][0] = -1;
1899     MP[1][0] = -1;
1900     MP[2][0] = -1;
1901     MT[0][0] = MT[1][0] = MT[2][0] = 0;
1902   } else if (leftd > midd) {
1903     fr = leftd-midd;
1904     for (j = 0; j <= midd; j++) {
1905       bss[j].CP = bss[j].DP = fr;
1906     }
1907     for (j = midd+1; j <= rightd; j++)
1908       bss[j].CP = bss[j].DP = -1;
1909     MP[0][fr] = -1;
1910     MP[1][fr] = -1;
1911     MP[2][fr] = -1;
1912     MT[0][fr] = MT[1][fr] = MT[2][fr] = 0;
1913   } else {
1914     fr = 0;
1915     for (j = 0; j < midd; j++) {
1916       bss[j].CP = bss[j].DP = 0;
1917     }
1918     for (j = midd; j <= rightd; j++) {
1919       bss[j].CP = bss[j].DP = 0;
1920     }
1921     MP[0][0] = -1;
1922     MP[1][0] = -1;
1923     MP[2][0] = -1;
1924     MT[0][0] = MT[1][0] = MT[2][0] = 0;
1925   }
1926
1927   bss[leftd].CC = 0;
1928   if (tb == 2) t = 0;
1929   else t = -g;
1930   for (j = leftd+1; j <= rightd; j++) {
1931     bss[j].CC = t = t-h;
1932     bss[j].DD = t-g;
1933   }
1934   bss[rightd+1].CC = MININT;
1935   bss[rightd+1].DD = MININT;
1936   if (tb == 1) bss[leftd].DD = 0;
1937   else bss[leftd].DD = -g;
1938   bss[leftd-1].CC = MININT;
1939   for (i = 1; i <= M; i++) {
1940     if (i > N-up) rightd--;
1941     if (leftd > 1) leftd--;
1942     wa = w[A[i]];
1943     if ((c = bss[leftd+1].CC-m) > (d = bss[leftd+1].DD-h)) {
1944       d = c;
1945       bss[leftd].DP = bss[leftd+1].CP;
1946     } else bss[leftd].DP = bss[leftd+1].DP;
1947     if ((ib = leftd+low-1+i) > 0) c = bss[leftd].CC+wa[B[ib]];
1948     if (d > c || ib <= 0) {
1949       c = d;
1950       bss[leftd].CP = bss[leftd].DP;
1951     }
1952     e = c-g;
1953     bss[leftd].DD = d;
1954     bss[leftd].CC = c;
1955     IP = bss[leftd].CP;
1956     if (leftd == midd) bss[leftd].CP = bss[leftd].DP = IP = i;
1957     for (curd=leftd+1; curd <= rightd; curd++) {
1958       if (curd != midd) {
1959         if ((c = c-m) > (e = e-h)) {
1960           e = c;
1961           IP = bss[curd-1].CP;
1962         }  /* otherwise, IP is unchanged */
1963         if ((c = bss[curd+1].CC-m) > (d = bss[curd+1].DD-h)) {
1964           d = c;
1965           bss[curd].DP = bss[curd+1].CP;
1966         } else {
1967           bss[curd].DP = bss[curd+1].DP;
1968         }
1969         c = bss[curd].CC + wa[B[curd+low-1+i]];
1970         if (c < d || c < e) {
1971           if (e > d) {
1972             c = e;
1973             bss[curd].CP = IP;
1974           } else {
1975             c = d;
1976             bss[curd].CP = bss[curd].DP;
1977           }
1978         } /* otherwise, CP is unchanged */
1979         bss[curd].CC = c;
1980         bss[curd].DD = d;
1981       } else {
1982         if ((c = c-m) > (e = e-h)) {
1983           e = c;
1984           MP[1][i] = bss[curd-1].CP;
1985           MT[1][i] = 2;
1986         } else {
1987           MP[1][i] = IP;
1988           MT[1][i] = 2;
1989         }
1990         if ((c = bss[curd+1].CC-m) > (d = bss[curd+1].DD-h)) {
1991           d = c;
1992           MP[2][i] = bss[curd+1].CP;
1993           MT[2][i] = 1;
1994         } else {
1995           MP[2][i] = bss[curd+1].DP;
1996           MT[2][i] = 1;
1997         }
1998         c = bss[curd].CC + wa[B[curd+low-1+i]];
1999         if (c < d || c < e) {
2000           if (e > d) {
2001             c = e;
2002             MP[0][i] = MP[1][i];
2003             MT[0][i] = 2;
2004           } else {
2005             c = d;
2006             MP[0][i] = MP[2][i];
2007             MT[0][i] = 1;
2008           }
2009         } else {
2010           MP[0][i] = i-1;
2011           MT[0][i] = 0;
2012         }
2013         if (c-g > e) {
2014           MP[1][i] = MP[0][i];
2015           MT[1][i] = MT[0][i];
2016         }
2017         if (c-g > d) {
2018           MP[2][i] = MP[0][i];
2019           MT[2][i] = MT[0][i];
2020         }
2021         bss[curd].CP = bss[curd].DP = IP = i;
2022         bss[curd].CC = c;
2023         bss[curd].DD = d;
2024       }
2025     }
2026   }
2027
2028   /* decide which path to be traced back */
2029   if (te == 1 && d+g > c) {
2030     k = bss[rightd].DP;
2031     l = 2;
2032   } else if (te == 2 && e+g > c) {
2033     k = IP;
2034     l = 1;
2035   } else {
2036     k = bss[rightd].CP;
2037     l = 0;
2038   }
2039   if (rmid > N-M) l = 2;
2040   else if (rmid < N-M) l = 1;
2041   v = c;
2042   }
2043   /* Conquer: Solve subproblems recursively */
2044
2045   /* trace back */
2046   r = -1;
2047   for (; k > -1; r=k, k=MP[l][r], l=MT[l][r]){
2048     FP[k] = r;
2049     FT[k] = l;  /* l=0,1,2 */
2050   }
2051   /* forward dividing */
2052   if (r == -1) { /* optimal alignment did not cross the middle diagonal */
2053     if (rmid < 0) {
2054       bg_align(A,B,M,N,rmid+1,up,tb,te,w,g,h,bss, sapp, last);
2055     }
2056     else {
2057       bg_align(A,B,M,N,low,rmid-1,tb,te,w,g,h,bss, sapp, last);
2058     }
2059   } else {
2060     k = r;
2061     l = FP[k];
2062     kt = FT[k];
2063
2064     /* first block */
2065     if (rmid < 0) {
2066       bg_align(A,B,r-1,r+rmid,rmid+1,min(up,r+rmid),tb,1,w,g,h,bss,sapp,last);
2067       DEL(1)
2068     } else if (rmid > 0) {
2069       bg_align(A,B,r,r+rmid-1,max(-r,low),rmid-1,tb,2,w,g,h,bss,sapp,last);
2070       INS(1)
2071     }
2072
2073     /* intermediate blocks */
2074     t2 = up-rmid-1;
2075     t3 = low-rmid+1;
2076     for (; l > -1; k = l, l = FP[k], kt = FT[k]) {
2077       if (kt == 0) { REP }
2078       else if (kt == 1) { /* right-hand side triangle */
2079         INS(1)
2080         t1 = l-k-1;
2081         bg_align(A+k,B+k+rmid+1,t1,t1,0,min(t1,t2),2,1,w,g,h,bss,sapp,last);
2082         DEL(1)
2083       }
2084       else { /* kt == 2, left-hand side triangle */
2085         DEL(1)
2086         t1 = l-k-1;
2087         bg_align(A+k+1,B+k+rmid,t1,t1,max(-t1,t3),0,1,2,w,g,h,bss,sapp,last);
2088         INS(1)
2089       }
2090     }
2091
2092     /* last block */
2093     if (N-M > rmid) {
2094       INS(1)
2095       t1 = k+rmid+1;
2096       bg_align(A+k,B+t1,M-k,N-t1,0,min(N-t1,t2),2,te,w,g,h,bss,sapp,last);
2097     } else if (N-M < rmid) {
2098       DEL(1)
2099       t1 = M-(k+1);
2100       bg_align(A+k+1,B+k+rmid,t1,N-(k+rmid),max(-t1,t3),0,1,te,w,g,h,
2101                bss,sapp,last);
2102     }
2103   }
2104   return(v);
2105 }
2106
2107 int B_ALIGN(const unsigned char *A, const unsigned char *B,
2108             int M, int N,
2109             int low, int up, int **W, int G, int H, int *S, int *nS,
2110             int MW, int MX, struct bdstr *bss)
2111 {
2112   int c, i, j;
2113   int g, h;
2114   size_t mj;
2115   int check_score;
2116   int **sapp, *sapp_v, *last, last_v;
2117
2118   g = G;
2119   h = H;
2120   sapp_v = S;
2121   sapp = &sapp_v;
2122
2123   last_v = 0;
2124   last = &last_v;
2125
2126   low = min(max(-M, low),min(N-M,0));
2127   up = max(min(N, up),max(N-M,0));
2128
2129   if (N <= 0) {
2130     if (M > 0) { DEL(M); }
2131     return -gap(M);
2132   }
2133   if (M <= 0) {
2134     INS(N);
2135     return -gap(N);
2136   }
2137   if (up-low+1 <= 1) {
2138     c = 0;
2139     for (i = 1; i <= M; i++) {
2140       REP;
2141       c += W[A[i]][B[i]];
2142     }
2143     return c;
2144   }
2145
2146   if (MT[0]==NULL) {
2147     mj = MX+1;
2148     MT[0] = (char *) ckalloc(mj);
2149     MT[1] = (char *) ckalloc(mj);
2150     MT[2] = (char *) ckalloc(mj);
2151     FT = (char *) ckalloc(mj);
2152
2153     mj *= sizeof(int);
2154     MP[0] = (int *) ckalloc(mj);
2155     MP[1] = (int *) ckalloc(mj);
2156     MP[2] = (int *) ckalloc(mj);
2157     FP = (int *) ckalloc(mj);
2158   }
2159
2160   c = bg_align(A,B,M,N,low,up,0,0,W,G,H,bss, sapp, last);
2161
2162   check_score = BCHECK_SCORE(A,B,M,N,S,W,G,H,nS);
2163
2164   free(FP); free(MP[2]); free(MP[1]); free(MP[0]);
2165   free(FT); free(MT[2]); free(MT[1]); free(MT[0]);
2166   MT[0]=NULL;
2167
2168   if (check_score != c)
2169     printf("\nBCheck_score=%d != %d\n", check_score,c);
2170   return c;
2171 }
2172
2173 int  do_walign (const unsigned char *aa0, int n0,
2174                 const unsigned char *aa1, int n1,
2175                 int frame,
2176                 struct pstruct *ppst,
2177                 struct f_struct *f_str,
2178                 struct a_res_str *a_res,
2179                 int *have_ares)
2180 {
2181   int hoff, optflag_s, optcut_s, optwid_s, n10, score;
2182   const unsigned char *aa1p;
2183   struct rstruct rst;
2184
2185 #ifdef TFASTA
2186   f_str->n10 = n10=aatran(aa1,f_str->aa1x,n1,frame);
2187   do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, &rst, &hoff);
2188   aa1p = f_str->aa1x;
2189
2190 #else
2191   n10 = n1;
2192   aa1p = aa1;
2193 #endif
2194
2195   a_res->res = f_str->res;
2196   *have_ares = 1;
2197
2198   if (ppst->sw_flag)
2199     return sw_walign(aa0, n0, aa1p, n10, ppst, f_str, a_res);
2200   else {
2201     optflag_s = ppst->param_u.fa.optflag;
2202     optcut_s = ppst->param_u.fa.optcut;
2203     optwid_s = ppst->param_u.fa.optwid;
2204     ppst->param_u.fa.optflag = 1;
2205     ppst->param_u.fa.optcut = 0;
2206     ppst->param_u.fa.optwid *= 2;
2207
2208     do_fasta(aa0, n0, aa1p, n10, ppst, f_str, &rst, &hoff);
2209
2210     if (rst.score[0]>0) {
2211       score=bd_walign(aa0, n0, aa1p, n10, ppst, f_str, hoff, a_res);
2212     }
2213     else {
2214       a_res->nres = 0;
2215       score=0;
2216     }
2217
2218     ppst->param_u.fa.optflag = optflag_s;
2219     ppst->param_u.fa.optcut = optcut_s;
2220     ppst->param_u.fa.optwid = optwid_s;
2221     return score;
2222   }
2223 }
2224
2225 void
2226 pre_cons(const unsigned char *aa1, int n1, int frame, struct f_struct *f_str) {
2227
2228 #ifdef TFASTA
2229   f_str->n10 = aatran(aa1,f_str->aa1x,n1,frame);
2230 #endif
2231 }
2232
2233 /* aln_func_vals - set up aln.qlfact, qlrev, llfact, llmult, frame, llrev */
2234 /* call from calcons, calc_id, calc_code */
2235 void
2236 aln_func_vals(int frame, struct a_struct *aln) {
2237
2238 #ifdef TFASTA
2239   aln->qlfact = 1;
2240   aln->llfact = 3;
2241   aln->llmult = 3;
2242   aln->qlrev = 0;
2243   aln->frame = frame;
2244   if (frame > 2) {
2245     aln->llrev = 1;
2246     aln->frame = 3 - frame;
2247   }
2248   else aln->llrev = 0;
2249 #else   /* FASTA */
2250   aln->llfact = aln->qlfact = aln->llmult = 1;
2251   aln->llrev = 0;
2252   if (frame > 0) aln->qlrev = 1;
2253   else aln->qlrev = 0;
2254   aln->frame = 0;
2255 #endif
2256 }
2257
2258 #include "a_mark.h"
2259
2260 int calcons(const unsigned char *aa0, int n0,
2261             const unsigned char *aa1, int n1,
2262             int *nc,
2263             struct a_struct *aln, struct a_res_str a_res,
2264             struct pstruct pst,
2265             char *seqc0, char *seqc1, char *seqca,
2266             struct f_struct *f_str)
2267 {
2268   int i0, i1, nn1;
2269   int op, lenc, nd, ns, itmp;
2270   const unsigned char *aa1p;
2271   char *sp0, *sp1, *spa, *sq;
2272   int *rp;
2273   int smins, mins;
2274
2275   if (pst.ext_sq_set) { sq = pst.sqx; }
2276   else { sq = pst.sq; }
2277
2278 #ifndef TFASTA
2279   aa1p = aa1;
2280   nn1 = n1;
2281 #else
2282   aa1p = f_str->aa1x;
2283   nn1 = f_str->n10;
2284 #endif
2285
2286   aln->amin0 = a_res.min0;
2287   aln->amax0 = a_res.max0;
2288   aln->amin1 = a_res.min1;
2289   aln->amax1 = a_res.max1;
2290                                         /* will we show all the start ?*/
2291   if (min(a_res.min0,a_res.min1) < aln->llen || aln->showall==1)
2292     if (a_res.min0 >= a_res.min1) {     /* aa0 extends more to left */
2293       smins=0;
2294       if (aln->showall==1) mins = a_res.min0;
2295       else mins = min(a_res.min0,aln->llcntx);
2296       aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2297       aln->smin0 = a_res.min0-mins;
2298       if ((mins-a_res.min1)>0) {
2299         memset(seqc1,' ',mins-a_res.min1);
2300         aancpy(seqc1+mins-a_res.min1,(char *)aa1p,a_res.min1,pst);
2301         aln->smin1 = 0;
2302       }
2303       else {
2304         aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
2305         aln->smin1 = a_res.min1-mins;
2306       }
2307     }
2308     else {
2309       smins=0;
2310       if (aln->showall == 1) mins=a_res.min1;
2311       else mins = min(a_res.min1,aln->llcntx);
2312       aancpy(seqc1,(char *)(aa1p+a_res.min1-mins),mins,pst);
2313       aln->smin1 = a_res.min1-mins;
2314       if ((mins-a_res.min0)>0) {
2315         memset(seqc0,' ',mins-a_res.min0);
2316         aancpy(seqc0+mins-a_res.min0,(char *)aa0,a_res.min0,pst);
2317         aln->smin0 = 0;
2318       }
2319       else {
2320         aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2321         aln->smin0 = a_res.min0-mins;
2322       }
2323     }
2324   else {
2325     mins= min(aln->llcntx,min(a_res.min0,a_res.min1));
2326     smins=mins;
2327     aln->smin0=a_res.min0 - mins;
2328     aln->smin1=a_res.min1 - mins;
2329     aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2330     aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
2331   }
2332   /* set the alignment code to zero for context */
2333   memset(seqca,0,mins);
2334
2335   /* TFASTA
2336   smins = mins = 0;
2337   aln->smin0=a_res.min0;
2338   aln->smin1=a_res.min1;
2339   */
2340
2341 /* now get the middle */
2342
2343   spa = seqca+mins;
2344   sp0 = seqc0+mins;
2345   sp1 = seqc1+mins;
2346   rp = a_res.res;
2347   lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = 0;
2348   i0 = a_res.min0;
2349   i1 = a_res.min1;
2350
2351   while (i0 < a_res.max0 || i1 < a_res.max1) {
2352     if (op == 0 && *rp == 0) {
2353       op = *rp++;
2354       lenc++;
2355       if ((itmp=pst.pam2[0][aa0[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
2356       else if (itmp == 0) { *spa = M_ZERO;}
2357       else {*spa = M_POS;}
2358       if (*spa == M_POS || *spa==M_ZERO) aln->nsim++;
2359
2360
2361       *sp0 = sq[aa0[i0++]];
2362       *sp1 = sq[aa1p[i1++]];
2363       if (toupper(*sp0) == toupper(*sp1)) {
2364         aln->nident++;
2365         *spa = M_IDENT;
2366       }
2367       else if (pst.nt_align) {
2368         if ((toupper(*sp0) == 'T' && toupper(*sp1) == 'U') ||
2369             (toupper(*sp0)=='U' && toupper(*sp1)=='T')) {
2370           aln->nident++;
2371           *spa = M_IDENT;
2372         }
2373         else if (toupper(*sp0) == 'N') aln->ngap_q++;
2374         else if (toupper(*sp1) == 'N') aln->ngap_l++;
2375       }
2376       sp0++; sp1++; spa++;
2377     }
2378     else {
2379       if (op==0) op = *rp++;
2380       if (op > 0) {
2381         *sp0++ = '-';
2382         *sp1++ = sq[aa1p[i1++]];
2383         *spa++ = M_DEL;
2384         op--;
2385         lenc++;
2386         aln->ngap_q++;
2387       }
2388       else {
2389         *sp0++ = sq[aa0[i0++]];
2390         *sp1++ = '-';
2391         *spa++ = M_DEL;
2392         op++;
2393         lenc++;
2394         aln->ngap_l++;
2395       }
2396     }
2397   }
2398
2399   *nc = lenc;
2400   *spa = '\0';
2401
2402 /*      now we have the middle, get the right end */
2403   if (!aln->llcntx_flg) {
2404     ns = mins + lenc + aln->llen;       /* show an extra line? */
2405     ns -= (itmp = ns %aln->llen);       /* itmp = left over on last line */
2406     if (itmp>aln->llen/2) ns += aln->llen;  /* more than 1/2 , use another*/
2407     nd = ns - (mins+lenc);              /* this much extra */
2408   }
2409   else nd = aln->llcntx;
2410
2411   if (nd > max(n0-a_res.max0,nn1-a_res.max1))
2412     nd = max(n0-a_res.max0,nn1-a_res.max1);
2413
2414   if (aln->showall==1) {
2415     nd = max(n0-a_res.max0,nn1-a_res.max1);     /* reset for showall=1 */
2416     /* get right end */
2417     aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-a_res.max0,pst);
2418     aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
2419     /* fill with blanks - this is required to use one 'nc' */
2420     memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
2421     memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
2422   }
2423   else {
2424     if ((nd-(n0-a_res.max0))>0) {
2425       aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,(n0-a_res.max0),pst);
2426       memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
2427     }
2428     else aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,nd,pst);
2429
2430     if ((nd-(nn1-a_res.max1))>0) {
2431       aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
2432       memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
2433     }
2434     else aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nd,pst);
2435   }
2436
2437   /*  fprintf(stderr,"%d\n",mins+lenc+nd); */
2438
2439   return mins+lenc+nd;
2440 }
2441
2442 int calcons_a(const unsigned char *aa0, unsigned char *aa0a, int n0,
2443               const unsigned char *aa1, int n1,
2444               int *nc,
2445               struct a_struct *aln,
2446               struct a_res_str a_res,
2447               struct pstruct pst,
2448               char *seqc0, char *seqc0a, char *seqc1, char *seqca,
2449               char *ann_arr, struct f_struct *f_str)
2450 {
2451   int i0, i1, nn1;
2452   int op, lenc, nd, ns, itmp;
2453   const unsigned char *aa1p;
2454   char *sp0, *sp0a, *sp1, *spa, *sq;
2455   int *rp;
2456   int smins, mins;
2457
2458   if (pst.ext_sq_set) {
2459     sq = pst.sqx;
2460   }
2461   else {
2462     sq = pst.sq;
2463   }
2464
2465 #ifndef TFASTA
2466   aa1p = aa1;
2467   nn1 = n1;
2468 #else
2469   aa1p = f_str->aa1x;
2470   nn1 = f_str->n10;
2471 #endif
2472
2473   aln->amin0 = a_res.min0;
2474   aln->amax0 = a_res.max0;
2475   aln->amin1 = a_res.min1;
2476   aln->amax1 = a_res.max1;
2477                                         /* will we show all the start ?*/
2478   /* will we show all the start ?*/
2479   if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1)
2480     if (a_res.min0>=a_res.min1) {              /* aa0 extends more to left */
2481       smins=0;
2482       if (aln->showall==1) mins = a_res.min0;
2483       else mins = min(a_res.min0,aln->llcntx);
2484       aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2485       aln->smin0 = a_res.min0-mins;
2486       if ((mins-a_res.min1)>0) {
2487         memset(seqc1,' ',mins-a_res.min1);
2488         aancpy(seqc1+mins-a_res.min1,(char *)aa1p,a_res.min1,pst);
2489         aln->smin1 = 0;
2490       }
2491       else {
2492         aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
2493         aln->smin1 = a_res.min1-mins;
2494       }
2495     }
2496     else {
2497       smins=0;
2498       if (aln->showall == 1) mins=a_res.min1;
2499       else mins = min(a_res.min1,aln->llcntx);
2500       aancpy(seqc1,(char *)(aa1p+a_res.min1-mins),mins,pst);
2501       aln->smin1 = a_res.min1-mins;
2502       if ((mins-a_res.min0)>0) {
2503         memset(seqc0,' ',mins-a_res.min0);
2504         aancpy(seqc0+mins-a_res.min0,(char *)aa0,a_res.min0,pst);
2505         aln->smin0 = 0;
2506       }
2507       else {
2508         aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2509         aln->smin0 = a_res.min0-mins;
2510       }
2511     }
2512   else {
2513     mins= min(aln->llcntx,min(a_res.min0,a_res.min1));
2514     smins=mins;
2515     aln->smin0=a_res.min0 - smins;
2516     aln->smin1=a_res.min1 - smins;
2517     aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2518     aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
2519   }
2520   /* set the alignment code to zero for context */
2521   memset(seqca,0,mins);
2522   memset(seqc0a,' ',mins);
2523
2524   /* TFASTA
2525   smins = mins = 0;
2526   aln->smin0=a_res.min0;
2527   aln->smin1=a_res.min1;
2528   */
2529
2530 /* now get the middle */
2531
2532   spa = seqca+mins;
2533   sp0 = seqc0+mins;
2534   sp0a = seqc0a+mins;
2535   sp1 = seqc1+mins;
2536   rp = a_res.res;
2537   lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = 0;
2538   i0 = a_res.min0;
2539   i1 = a_res.min1;
2540
2541   while (i0 < a_res.max0 || i1 < a_res.max1) {
2542     if (op == 0 && *rp == 0) {
2543       op = *rp++;
2544       lenc++;
2545       if ((itmp=pst.pam2[0][aa0[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
2546       else if (itmp == 0) { *spa = M_ZERO;}
2547       else {*spa = M_POS;}
2548       if (*spa == M_POS || *spa==M_ZERO) aln->nsim++;
2549
2550       *sp0a++ = ann_arr[aa0a[i0]];
2551
2552       *sp0 = sq[aa0[i0++]];
2553       *sp1 = sq[aa1p[i1++]];
2554
2555       if (toupper(*sp0) == toupper(*sp1)) {
2556         aln->nident++;
2557         *spa = M_IDENT;
2558       }
2559       else if (pst.nt_align) {
2560         if ((toupper(*sp0) == 'T' && toupper(*sp1) == 'U') ||
2561             (toupper(*sp0)=='U' && toupper(*sp1)=='T')) {
2562           aln->nident++;
2563           *spa = M_IDENT;
2564         }
2565         else if (toupper(*sp0) == 'N') aln->ngap_q++;
2566         else if (toupper(*sp1) == 'N') aln->ngap_l++;
2567       }
2568       sp0++; sp1++; spa++;
2569     }
2570     else {
2571       if (op==0) op = *rp++;
2572       if (op>0) {
2573         *sp0++ = '-';
2574         *sp1++ = sq[aa1p[i1++]];
2575         *spa++ = M_DEL;
2576         *sp0a++ = ' ';
2577         op--;
2578         lenc++;
2579         aln->ngap_q++;
2580       }
2581       else {
2582         *sp0a++ = ann_arr[aa0a[i0]];
2583         *sp0++ = sq[aa0[i0++]];
2584         *sp1++ = '-';
2585         *spa++ = M_DEL;
2586         op++;
2587         lenc++;
2588         aln->ngap_l++;
2589       }
2590     }
2591   }
2592
2593   *nc = lenc;
2594   *sp0a = *spa = '\0';
2595
2596 /*      now we have the middle, get the right end */
2597   if (!aln->llcntx_flg) {
2598     ns = mins + lenc + aln->llen;       /* show an extra line? */
2599     ns -= (itmp = ns %aln->llen);       /* itmp = left over on last line */
2600     if (itmp>aln->llen/2) ns += aln->llen;  /* more than 1/2 , use another*/
2601     nd = ns - (mins+lenc);              /* this much extra */
2602   }
2603   else nd = aln->llcntx;
2604
2605   if (nd > max(n0-a_res.max0,nn1-a_res.max1))
2606     nd = max(n0-a_res.max0,nn1-a_res.max1);
2607
2608   if (aln->showall==1) {
2609     nd = max(n0-a_res.max0,nn1-a_res.max1);     /* reset for showall=1 */
2610     /* get right end */
2611     aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-a_res.max0,pst);
2612     aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
2613     /* fill with blanks - this is required to use one 'nc' */
2614     memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
2615     memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
2616   }
2617   else {
2618     if ((nd-(n0-a_res.max0))>0) {
2619       aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,(n0-a_res.max0),pst);
2620       memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
2621     }
2622     else aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,nd,pst);
2623
2624     if ((nd-(nn1-a_res.max1))>0) {
2625       aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
2626       memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
2627     }
2628     else aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nd,pst);
2629   }
2630
2631   /*  fprintf(stderr,"%d\n",mins+lenc+nd); */
2632
2633   return mins+lenc+nd;
2634 }
2635
2636 static void
2637 update_code(char *al_str, int al_str_max, int op, int op_cnt) {
2638
2639   char op_char[5]={"=-+*"};
2640   char tmp_cnt[20];
2641
2642   sprintf(tmp_cnt,"%c%d",op_char[op],op_cnt);
2643   strncat(al_str,tmp_cnt,al_str_max);
2644 }
2645
2646
2647 /* build an array of match/ins/del - length strings */
2648 int calc_code(const unsigned char *aa0, int n0,
2649               const unsigned char *aa1, int n1,
2650               struct a_struct *aln, struct a_res_str a_res,
2651               struct pstruct pst,
2652               char *al_str, int al_str_n, struct f_struct *f_str)
2653 {
2654   int i0, i1, nn1;
2655   int op, lenc;
2656   int p_op, op_cnt;
2657   const unsigned char *aa1p;
2658   char sp0, sp1, *sq;
2659   int *rp;
2660
2661   if (pst.ext_sq_set) {
2662     sq = pst.sqx;
2663   }
2664   else {
2665     sq = pst.sq;
2666   }
2667
2668 #ifndef TFASTA
2669   aa1p = aa1;
2670   nn1 = n1;
2671 #else
2672   aa1p = f_str->aa1x;
2673   nn1 = f_str->n10;
2674 #endif
2675
2676   aln->amin0 = a_res.min0;
2677   aln->amax0 = a_res.max0;
2678   aln->amin1 = a_res.min1;
2679   aln->amax1 = a_res.max1;
2680
2681   rp = a_res.res;
2682   lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = p_op = 0;
2683   op_cnt = 0;
2684
2685   i0 = a_res.min0;
2686   i1 = a_res.min1;
2687
2688   while (i0 < a_res.max0 || i1 < a_res.max1) {
2689     if (op == 0 && *rp == 0) {
2690
2691       if (pst.pam2[0][aa0[i0]][aa1p[i1]]>=0) { aln->nsim++;}
2692
2693       sp0 = sq[aa0[i0++]];
2694       sp1 = sq[aa1p[i1++]];
2695
2696       if (p_op == 0 || p_op==3) {
2697         if (sp0 != '*' && sp1 != '*') {
2698           if (p_op == 3) {
2699             update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
2700             op_cnt = 1; p_op = 0;
2701           }
2702           else {op_cnt++;}
2703         }
2704         else {
2705           update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
2706           op_cnt = 1; p_op = 3;
2707         }
2708       }
2709       else {
2710         update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
2711         op_cnt = 1; p_op = 0;
2712       }
2713
2714       op = *rp++;
2715       lenc++;
2716
2717       if (toupper(sp0) == toupper(sp1)) aln->nident++;
2718       else if (pst.nt_align) {
2719         if ((toupper(sp0) == 'T' && toupper(sp1) == 'U') ||
2720             (toupper(sp0)=='U' && toupper(sp1)=='T')) aln->nident++;
2721         else if (toupper(sp0) == 'N') aln->ngap_q++;
2722         else if (toupper(sp1) == 'N') aln->ngap_l++;
2723       }
2724     }
2725     else {
2726       if (op==0) op = *rp++;
2727       if (op>0) {
2728         if (p_op == 1) { op_cnt++;}
2729         else {
2730           update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
2731           op_cnt = 1; p_op = 1;
2732         }
2733         op--; lenc++; i1++; aln->ngap_q++;
2734       }
2735       else {
2736         if (p_op == 2) { op_cnt++;}
2737         else {
2738           update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
2739           op_cnt = 1; p_op = 2;
2740         }
2741         op++; lenc++; i0++; aln->ngap_l++;
2742       }
2743     }
2744   }
2745   update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
2746
2747   return lenc;
2748 }
2749
2750 int calc_id(const unsigned char *aa0, int n0,
2751             const unsigned char *aa1, int n1,
2752             struct a_struct *aln,
2753             struct a_res_str a_res,
2754             struct pstruct pst,
2755             struct f_struct *f_str)
2756 {
2757   int i0, i1, nn1;
2758   int op, lenc;
2759   int sp0, sp1;
2760   const unsigned char *aa1p;
2761   int *rp;
2762   char *sq;
2763
2764   if (pst.ext_sq_set) { sq = pst.sqx; }
2765   else { sq = pst.sq; }
2766
2767 #ifndef TFASTA
2768   aa1p = aa1;
2769   nn1 = n1;
2770 #else
2771   aa1p = f_str->aa1x;
2772   nn1 = f_str->n10;
2773 #endif
2774
2775   aln->amin0 = a_res.min0;
2776   aln->amax0 = a_res.max0;
2777   aln->amin1 = a_res.min1;
2778   aln->amax1 = a_res.max1;
2779
2780   rp = a_res.res;
2781   lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = 0;
2782   i0 = a_res.min0;
2783   i1 = a_res.min1;
2784
2785   while (i0 < a_res.max0 || i1 < a_res.max1) {
2786     if (op == 0 && *rp == 0) {
2787       op = *rp++;
2788       lenc++;
2789
2790       if (pst.pam2[0][aa0[i0]][aa1p[i1]]>=0) { aln->nsim++;}
2791
2792       sp0 = sq[aa0[i0++]];
2793       sp1 = sq[aa1p[i1++]];
2794       if (toupper(sp0) == toupper(sp1)) {aln->nident++;}
2795       else if (pst.nt_align) {
2796         if ((toupper(sp0)=='T' && toupper(sp1)== 'U')||
2797           (toupper(sp0)=='U' && toupper(sp1)=='T')) {aln->nident++;}
2798         else if (toupper(sp0) == 'N') aln->ngap_q++;
2799         else if (toupper(sp1) == 'N') aln->ngap_l++;
2800       }
2801     }
2802     else {
2803       if (op==0) op = *rp++;
2804       if (op>0) {op--; lenc++; i1++; aln->ngap_q++;}
2805       else {op++; lenc++; i0++; aln->ngap_l++; }
2806     }
2807   }
2808   return lenc;
2809 }
2810
2811 #ifdef PCOMPLIB
2812
2813 #include "w_mw.h"
2814
2815 void
2816 update_params(struct qmng_str *qm_msg, struct pstruct *ppst)
2817 {
2818   ppst->n0 = qm_msg->n0;
2819 }
2820 #endif