website/archive/binaries/mac/src/fasta34/dropnsw.c

   1 /* copyright (c) 1994, 1995, 1996 William R. Pearson */
   2
   3 /* $Name: fa_34_26_5 $ - $Id: dropnsw.c,v 1.35 2006/10/19 14:49:14 wrp Exp $ */
   4
   5 /*
   6   this is a slower version of dropgsw.c that implements the Smith-Waterman
   7   algorithm.  It lacks the shortcuts in dropgsw.c that prevent scores less
   8   than the penalty for the first residue in a gap from being generated.
   9
  10   Thus, dropnsw.c should be used for tests with very large gap penalties,
  11   and is more appropriate for programs like prss3, which are interested
  12   in accurate low scores.
  13 */
  14
  15 /* the do_walign() code in this file is not thread_safe */
  16 /* init_work(), do_work(), are thread safe */
  17
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20 #include <string.h>
  21 #include <ctype.h>
  22 #include <math.h>
  23
  24 #include "defs.h"
  25 #include "param.h"
  26
  27 static char *verstr="3.5 Sept 2006";
  28
  29 struct swstr { int H, E;};
  30
  31 struct f_struct {
  32   struct swstr *ss;
  33   struct swstr *f_ss;
  34   struct swstr *r_ss;
  35   int *waa_s, *waa_a;
  36   int **pam2p[2];
  37   int *res;
  38   double aa0_f[MAXSQ];
  39   double *kar_p;
  40 };
  41
  42 #define DROP_INTERN
  43 #include "drop_func.h"
  44
  45 extern int do_karlin(const unsigned char *aa1, int n1,
  46                      int **pam2, struct pstruct *ppst,
  47                      double *aa0_f, double *kar_p, double *lambda, double *H);
  48 extern void aancpy(char *to, char *from, int count, struct pstruct pst);
  49 int ALIGN(const unsigned char *A, const unsigned char *B, int M, int N,
  50           int **W, int IW, int G, int H, int *S, int *NC,
  51           struct f_struct *f_str);
  52
  53 /* initialize for Smith-Waterman optimal score */
  54
  55 void init_work (unsigned char *aa0, int n0,
  56                 struct pstruct *ppst,
  57                 struct f_struct **f_arg)
  58 {
  59    int maxn0;
  60    int *pwaa_s, *pwaa_a;
  61    int e, f, i, j, q;
  62    int *res;
  63    struct f_struct *f_str;
  64    int **pam2p;
  65    struct swstr *ss, *f_ss, *r_ss;
  66    int nsq, ip;
  67
  68    if (ppst->ext_sq_set) {
  69      nsq = ppst->nsqx; ip = 1;
  70    }
  71    else {
  72      nsq = ppst->nsq; ip = 0;
  73    }
  74
  75    f_str = (struct f_struct *)calloc(1,sizeof(struct f_struct));
  76
  77    /* allocate space for the scoring arrays */
  78    maxn0 = n0 + 2;
  79    if ((ss = (struct swstr *) calloc (maxn0, sizeof (struct swstr)))
  80          == NULL) {
  81      fprintf (stderr, "cannot allocate ss array %3d\n", n0);
  82      exit (1);
  83    }
  84    ss++;
  85    f_str->ss = ss;
  86
  87    if ((f_ss = (struct swstr *) calloc (maxn0, sizeof (struct swstr)))
  88        == NULL) {
  89      fprintf (stderr, "cannot allocate f_ss array %3d\n", n0);
  90      exit (1);
  91    }
  92    f_ss++;
  93    f_str->f_ss = f_ss;
  94
  95    if ((r_ss = (struct swstr *) calloc (n0+2, sizeof (struct swstr)))
  96        == NULL) {
  97      fprintf (stderr, "cannot allocate r_ss array %3d\n", n0);
  98      exit (1);
  99    }
 100    r_ss++;
 101    f_str->r_ss = r_ss;
 102
 103    /* initialize variable (-S) pam matrix */
 104    if ((f_str->waa_s= (int *)calloc((nsq+1)*(n0+1),sizeof(int))) == NULL) {
 105      fprintf(stderr,"cannot allocate waa_s array %3d\n",nsq*n0);
 106      exit(1);
 107    }
 108
 109    if ((f_str->pam2p[1]= (int **)calloc((n0+1),sizeof(int *))) == NULL) {
 110      fprintf(stderr,"cannot allocate pam2p[1] array %3d\n",n0);
 111      exit(1);
 112    }
 113
 114    pam2p = f_str->pam2p[1];
 115    if ((pam2p[0]=(int *)calloc((nsq+1)*(n0+1),sizeof(int))) == NULL) {
 116      fprintf(stderr,"cannot allocate pam2p[1][] array %3d\n",nsq*n0);
 117      exit(1);
 118    }
 119
 120    for (i=1; i<n0; i++) {
 121      pam2p[i]= pam2p[0] + (i*(nsq+1));
 122    }
 123
 124    /* initialize universal (alignment) matrix */
 125    if ((f_str->waa_a= (int *)calloc((nsq+1)*(n0+1),sizeof(int))) == NULL) {
 126      fprintf(stderr,"cannot allocate waa_a struct %3d\n",nsq*n0);
 127      exit(1);
 128    }
 129
 130    if ((f_str->pam2p[0]= (int **)calloc((n0+1),sizeof(int *))) == NULL) {
 131      fprintf(stderr,"cannot allocate pam2p[1] array %3d\n",n0);
 132      exit(1);
 133    }
 134
 135    pam2p = f_str->pam2p[0];
 136    if ((pam2p[0]=(int *)calloc((nsq+1)*(n0+1),sizeof(int))) == NULL) {
 137      fprintf(stderr,"cannot allocate pam2p[1][] array %3d\n",nsq*n0);
 138      exit(1);
 139    }
 140
 141    for (i=1; i<n0; i++) {
 142      pam2p[i]= pam2p[0] + (i*(nsq+1));
 143    }
 144
 145    /*
 146       pwaa effectively has a sequence profile --
 147        pwaa[0..n0-1] has pam score for residue 0 (-BIGNUM)
 148        pwaa[n0..2n0-1] has pam scores for residue 1 (A)
 149        pwaa[2n0..3n-1] has pam scores for residue 2 (R), ...
 150
 151        thus: pwaa = f_str->waa_s + (*aa1p++)*n0; sets up pwaa so that
 152        *pwaa++ rapidly moves though the scores of the aa1p[] position
 153        without further indexing
 154
 155        For a real sequence profile, pwaa[0..n0-1] vs ['A'] could have
 156        a different score in each position.
 157    */
 158
 159    if (ppst->pam_pssm) {
 160      pwaa_s = f_str->waa_s;
 161      pwaa_a = f_str->waa_a;
 162      for (e = 0; e <=nsq; e++)  {       /* for each residue in the alphabet */
 163        for (f = 0; f < n0; f++) {       /* for each position in aa0 */
 164          *pwaa_s++ = f_str->pam2p[ip][f][e] = ppst->pam2p[ip][f][e];
 165          *pwaa_a++ = f_str->pam2p[0][f][e]  = ppst->pam2p[0][f][e];
 166        }
 167      }
 168    }
 169    else {       /* initialize scanning matrix */
 170      pwaa_s = f_str->waa_s;
 171      pwaa_a = f_str->waa_a;
 172      for (e = 0; e <=nsq; e++)  /* for each residue in the alphabet */
 173        for (f = 0; f < n0; f++) {       /* for each position in aa0 */
 174          *pwaa_s++ = f_str->pam2p[ip][f][e]= ppst->pam2[ip][e][aa0[f]];
 175          *pwaa_a++ = f_str->pam2p[0][f][e] = ppst->pam2[0][e][aa0[f]];
 176        }
 177    }
 178
 179    maxn0 = max(3*n0/2,MIN_RES);
 180    if ((res = (int *)calloc((size_t)maxn0,sizeof(int)))==NULL) {
 181      fprintf(stderr,"cannot allocate alignment results array %d\n",maxn0);
 182      exit(1);
 183    }
 184    f_str->res = res;
 185
 186    *f_arg = f_str;
 187 }
 188
 189 void close_work (const unsigned char *aa0, int n0,
 190                  struct pstruct *ppst, struct f_struct **f_arg)
 191 {
 192   struct f_struct *f_str;
 193
 194   f_str = *f_arg;
 195
 196   if (f_str != NULL) {
 197     if (f_str->kar_p !=NULL) free(f_str->kar_p);
 198     f_str->ss--;
 199     free(f_str->ss);
 200     free(f_str->res);
 201     free(f_str->waa_a);
 202     free(f_str->pam2p[0][0]);
 203     free(f_str->pam2p[0]);
 204     free(f_str->waa_s);
 205     free(f_str->pam2p[1][0]);
 206     free(f_str->pam2p[1]);
 207
 208     free(f_str);
 209     *f_arg = NULL;
 210   }
 211 }
 212
 213
 214 /* pstring1 is a message to the manager, currently 512 */
 215 /*void get_param(struct pstruct *pstr,char *pstring1)*/
 216 void    get_param (struct pstruct *pstr, char *pstring1, char *pstring2)
 217 {
 218   char psi_str[120];
 219
 220   char *pg_str="Smith-Waterman";
 221
 222   if (pstr->pam_pssm) { strncpy(psi_str,"-PSI",sizeof(psi_str));}
 223   else { psi_str[0]='\0';}
 224
 225 #ifdef OLD_FASTA_GAP
 226    sprintf (pstring1, " %s (%s) function [%s matrix%s (%d:%d)%s], gap-penalty: %d/%d",
 227 #else
 228    sprintf (pstring1, " %s (%s) function [%s matrix%s (%d:%d)%s], open/ext: %d/%d",
 229 #endif
 230             pg_str, verstr, pstr->pamfile, psi_str, pstr->pam_h,pstr->pam_l,
 231             (pstr->ext_sq_set)?"xS":"\0", pstr->gdelval, pstr->ggapval);
 232
 233    if (pstring2 != NULL) {
 234 #ifdef OLD_FASTA_GAP
 235      sprintf(pstring2,"; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)%s\n; pg_gap-pen: %d %d\n",
 236 #else
 237      sprintf(pstring2,"; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)%s\n; pg_open-ext: %d %d\n",
 238 #endif
 239              pg_str,verstr,psi_str,pstr->pam_h,pstr->pam_l,
 240              (pstr->ext_sq_set)?"xS":"\0",pstr->gdelval,pstr->ggapval);
 241    }
 242 }
 243
 244
 245 void do_work (const unsigned char *aa0, int n0,
 246               const unsigned char *aa1, int n1,
 247               int frame,
 248               struct pstruct *ppst, struct f_struct *f_str,
 249               int qr_flg,
 250               struct rstruct *rst)
 251 {
 252    const unsigned char *aa0p, *aa1p;
 253    register struct swstr *ssj;
 254    struct swstr *ss, *f_ss, *r_ss;
 255    register int *pwaa;
 256    int *waa;
 257    register int i, j;
 258    int     e, f, h, p;
 259    int     q, r, m;
 260    int     score;
 261
 262    double lambda, H, K;
 263
 264    rst->escore = 1.0;
 265    rst->segnum = rst->seglen = 1;
 266
 267    waa = f_str->waa_s;
 268    ss = f_str->ss;
 269    f_ss = f_str->f_ss;
 270    r_ss = f_str->r_ss;
 271
 272 #ifdef OLD_FASTA_GAP
 273    q = -(ppst->gdelval - ppst->ggapval);
 274 #else
 275    q = -ppst->gdelval;
 276 #endif
 277    r = -ppst->ggapval;
 278    m = q + r;
 279
 280    /* initialize 0th row */
 281    for (ssj=ss; ssj<&ss[n0]; ssj++) {
 282      ssj->H = 0;
 283      ssj->E = -q;
 284    }
 285
 286    score = 0;
 287    aa1p = aa1;
 288    while (*aa1p) {
 289      h = p = 0;
 290      f = -q;
 291      pwaa = waa + (*aa1p++ * n0);
 292      for (ssj = ss, aa0p = aa0; ssj < ss+n0; ssj++) {
 293        if ((h =   h     - m) > (f =   f     - r)) f = h;
 294        if ((h = ssj->H - m) > (e = ssj->E - r)) e = h;
 295        h = p + *pwaa++;
 296        if (h < 0 ) h = 0;
 297        if (h < f ) h = f;
 298        if (h < e ) h = e;
 299        p = ssj->H;
 300        ssj->H = h;
 301        ssj->E = e;
 302        if (h > score) score = h;
 303      }
 304    }                            /* done with forward pass */
 305
 306    rst->score[0] = score;
 307
 308    if(ppst->zsflag == 6 || ppst->zsflag == 16 &&
 309      (do_karlin(aa1, n1, ppst->pam2[0], ppst,f_str->aa0_f,
 310                 f_str->kar_p, &lambda, &H)>0)) {
 311      rst->comp = 1.0/lambda;
 312      rst->H = H;
 313    }
 314   else {rst->comp = rst->H = -1.0;}
 315 }                               /* here we should be all done */
 316
 317 void    do_opt (const unsigned char *aa0, int n0,
 318                 const unsigned char *aa1, int n1,
 319                 int frame,
 320                 struct pstruct *pst, struct f_struct *f_str,
 321                 struct rstruct *rstr)
 322 {
 323 }
 324
 325 int do_walign (const unsigned char *aa0, const int n0,
 326                const unsigned char *aa1, const int n1,
 327                int frame,
 328                struct pstruct *ppst,
 329                struct f_struct *f_str,
 330                struct a_res_str *a_res,
 331                int *have_ares )
 332 {
 333    const unsigned char *aa0p, *aa1p;
 334    register int *pwaa;
 335    register int i, j;
 336    register struct swstr *ssj;
 337    struct swstr *f_ss, *r_ss, *ss;
 338    int *res, *waa;
 339    int e, f, h, p;
 340    int     q, r, m;
 341    int     score;
 342    int cost, I, J, K, L;
 343
 344    ss = f_str->ss;
 345
 346    res = f_str->res;
 347    waa = f_str->waa_a;  /* this time use universal pam2[0] */
 348
 349 #ifdef OLD_FASTA_GAP
 350    q = -(ppst->gdelval - ppst->ggapval);
 351 #else
 352    q = -ppst->gdelval;
 353 #endif
 354
 355    r = -ppst->ggapval;
 356    m = q + r;
 357
 358    /* initialize 0th row */
 359    for (ssj=ss; ssj<ss+n0; ssj++) {
 360      ssj->H = 0;
 361      ssj->E = -q;
 362    }
 363
 364    score = 0;
 365    aa1p = aa1;
 366    i = 0;
 367    while (*aa1p) {
 368      h = p = 0;
 369      f = -q;
 370      pwaa = waa + (*aa1p++ * n0);
 371      for (ssj = ss, aa0p = aa0; ssj < ss+n0; ssj++) {
 372        if ((h =   h     - m) > /* gap open from left best */
 373            /* gap extend from left gapped */
 374            (f =   f     - r)) f = h;    /* if better, use new gap opened */
 375        if ((h = ssj->H - m) >   /* gap open from up best */
 376            /* gap extend from up gap */
 377            (e = ssj->E - r)) e = h;     /* if better, use new gap opened */
 378        h = p + *pwaa++;         /* diagonal match */
 379        if (h < 0 ) h = 0;       /* ?  < 0, reset to 0 */
 380        if (h < f ) h = f;       /* left gap better, reset */
 381        if (h < e ) h = e;       /* up gap better, reset */
 382        p = ssj->H;              /* save previous best score */
 383        ssj->H = h;              /* save (new) up diag-matched */
 384        ssj->E = e;              /* save upper gap opened */
 385        if (h > score) {         /* ? new best score */
 386          score = h;             /* save best */
 387          I = i;                 /* row */
 388          J = (int)(ssj-ss);     /* column */
 389        }
 390      }
 391      i++;
 392    }                            /* done with forward pass */
 393    if (score <= 0) return 0;
 394
 395   /* to get the start point, go backwards */
 396
 397    /* 18-June-2003 fix bug in backtracking code to identify start of
 398       alignment.  Code used pam2[0][aa0[j]][aa1[i]] instead of
 399       pam2p[0][j][aa1[i]].  Ideally, it would use waa_a.
 400    */
 401
 402   cost = K = L = 0;
 403   for (ssj=ss+J; ssj>=ss; ssj--) ssj->H= ssj->E= -1;
 404
 405   for (i=I; i>=0; i--) {
 406     h = f = -1;
 407     p = (i == I) ? 0 : -1;
 408     for (ssj=ss+J, j= J; ssj>=ss; ssj--,j--) {
 409       f = max (f,h-q)-r;
 410       ssj->E=max(ssj->E,ssj->H-q)-r;
 411       h = max(max(ssj->E,f),p+f_str->pam2p[0][j][aa1[i]]);
 412       p = ssj->H;
 413       ssj->H=h;
 414       if (h > cost) {
 415         cost = h;
 416         K = i;
 417         L = (int)(ssj-ss);
 418         if (cost >= score) goto found;
 419       }
 420     }
 421   }
 422
 423 found:
 424
 425   /*  printf(" %d: L: %3d-%3d/%3d; K: %3d-%3d/%3d\n",score,L,J,n0,K,I,n1); */
 426
 427 /* in the f_str version, the *res array is already allocated at 4*n0/3 */
 428
 429   a_res->res = f_str->res;
 430   *have_ares = 1;
 431   a_res->max0 = J+1; a_res->min0 = L; a_res->max1 = I+1; a_res->min1 = K;
 432
 433 /*  ALIGN(&aa1[K-1],&aa0[L-1],I-K+1,J-L+1,ppst->pam2[0],q,r,res,nres,f_str); */
 434
 435 /* this code no longer refers to aa0[], it used pam2p[0][L] instead */
 436   ALIGN(&aa0[L-1],&aa1[K-1],J-L+1,I-K+1,f_str->pam2p[0],L,q,r,
 437         a_res->res,&a_res->nres,f_str);
 438
 439 /*   DISPLAY(&aa0[L-1],&aa1[K-1],J-L+1,I-K+1,res,L,K,ppst->sq); */
 440
 441   return score;
 442 }
 443
 444 static int CHECK_SCORE(const unsigned char *A, const unsigned char *B, int M, int N,
 445                        int *S, int **W, int IW, int G, int H, int *nres);
 446
 447 #define gap(k)  ((k) <= 0 ? 0 : g+h*(k))        /* k-symbol indel cost */
 448
 449 /* Append "Delete k" op */
 450 #define DEL(k)                          \
 451 { if (*last < 0)                        \
 452     *last = (*sapp)[-1] -= (k);         \
 453   else {                                \
 454     *last = (*sapp)[0] = -(k);          \
 455     (*sapp)++;                          \
 456   }                                     \
 457 }
 458
 459 /* Append "Insert k" op */
 460 #define INS(k)                          \
 461 { if (*last > 0)                        \
 462     *last = (*sapp)[-1] += (k);         \
 463   else {                                \
 464     *last = (*sapp)[0] = (k);           \
 465     (*sapp)++;                          \
 466   }                                     \
 467 }
 468
 469 #define REP { *last = (*sapp)[0] = 0; (*sapp)++; } /* Append "Replace" op */
 470
 471 /*
 472 #define XTERNAL
 473 #include "upam.h"
 474
 475 void
 476 print_seq_prof(unsigned char *A, int M,
 477                unsigned char *B, int N,
 478                int **w, int iw, int dir) {
 479   char c_max;
 480   int i_max, j_max, i,j;
 481
 482   char *c_dir="LRlr";
 483
 484   for (i=1; i<=min(60,M); i++) {
 485     fprintf(stderr,"%c",aa[A[i]]);
 486   }
 487   fprintf(stderr, - %d\n,M);
 488
 489   for (i=0; i<min(60,M); i++) {
 490     i_max = -1;
 491     for (j=1; j<21; j++) {
 492       if (w[iw+i][j]> i_max) {
 493         i_max = w[iw+i][j];
 494         j_max = j;
 495       }
 496     }
 497     fprintf(stderr,"%c",aa[j_max]);
 498   }
 499   fputc(':',stderr);
 500   for (i=1; i<=min(60,N); i++) {
 501     fprintf(stderr,"%c",aa[B[i]]);
 502   }
 503   fprintf(stderr," -%c: %d,%d\n",c_dir[dir],M,N);
 504 }
 505 */
 506
 507 /* align(A,B,M,N,tb,te) returns the cost of an optimum conversion between
 508    A[1..M] and B[1..N] that begins(ends) with a delete if tb(te) is zero
 509    and appends such a conversion to the current script.                   */
 510
 511 static int
 512 align(const unsigned char *A, const unsigned char *B, int M, int N,
 513       int tb, int te, int **w, int iw, int g, int h,
 514       struct f_struct *f_str, int dir,
 515       int **sapp, int *last)
 516 {
 517   int   midi, midj, type;       /* Midpoint, type, and cost */
 518   int midc;
 519   int c1, c2;
 520
 521 { register int   i, j;
 522   register int c, e, d, s;
 523   int m, t, *wa;
 524   struct swstr *f_ss, *r_ss;
 525
 526 /*   print_seq_prof(A,M,B,N,w,iw,dir); */
 527
 528   m = g + h;
 529
 530   f_ss = f_str->f_ss;
 531   r_ss = f_str->r_ss;
 532
 533 /* Boundary cases: M <= 1 or N == 0 */
 534
 535   if (N <= 0) {
 536     if (M > 0) {
 537       DEL(M)
 538     }
 539     return -gap(M);
 540   }
 541
 542   if (M <= 1) {
 543     if (M <= 0){
 544       INS(N)
 545       return -gap(N); }
 546     if (tb < te) tb = te;
 547     midc = (tb-h) - gap(N);
 548     midj = 0;
 549 /*  wa = w[A[1]]; */
 550     wa = w[iw];
 551     for (j = 1; j <= N; j++) {
 552       c = -gap(j-1) + wa[B[j]] - gap(N-j);
 553       if (c > midc) { midc = c; midj = j;}
 554     }
 555     if (midj == 0) {
 556       DEL(1)
 557       INS(N)
 558     }
 559     else  {
 560       if (midj > 1) { INS(midj-1)}
 561       REP
 562       if (midj < N) { INS(N-midj)}
 563     }
 564     return midc;
 565   }
 566
 567 /* Divide: Find optimum midpoint (midi,midj) of cost midc */
 568
 569   midi = M/2;           /* Forward phase:                          */
 570   f_ss[0].H = 0;        /*   Compute H(M/2,k) & E(M/2,k) for all k */
 571   t = -g;
 572   for (j = 1; j <= N; j++)
 573     { f_ss[j].H = t = t-h;
 574       f_ss[j].E = t-g;
 575     }
 576   t = tb;
 577   for (i = 1; i <= midi; i++)
 578     { s = f_ss[0].H;
 579       f_ss[0].H = c = t = t-h;
 580       e = t-g;
 581 /*    wa = w[A[i]]; */
 582       wa = w[iw+i-1];
 583       for (j = 1; j <= N; j++)
 584         { if ((c =   c   - m) > (e =   e   - h)) e = c;
 585           if ((c = f_ss[j].H - m) > (d = f_ss[j].E - h)) d = c;
 586           c = s + wa[B[j]];
 587           if (e > c) c = e;
 588           if (d > c) c = d;
 589           s = f_ss[j].H;
 590           f_ss[j].H = c;
 591           f_ss[j].E = d;
 592         }
 593     }
 594   f_ss[0].E = f_ss[0].H;
 595
 596   r_ss[N].H = 0;                /* Reverse phase:                  */
 597   t = -g;                       /*   Compute R(M/2,k) & S(M/2,k) for all k */
 598   for (j = N-1; j >= 0; j--)
 599     { r_ss[j].H = t = t-h;
 600       r_ss[j].E = t-g;
 601     }
 602   t = te;
 603   for (i = M-1; i >= midi; i--)
 604     { s = r_ss[N].H;
 605       r_ss[N].H = c = t = t-h;
 606       e = t-g;
 607 /*    wa = w[A[i+1]]; */
 608       wa = w[iw+i];
 609       for (j = N-1; j >= 0; j--)
 610         { if ((c =   c   - m) > (e =   e   - h)) e = c;
 611           if ((c = r_ss[j].H - m) > (d = r_ss[j].E - h)) d = c;
 612           c = s + wa[B[j+1]];
 613           if (e > c) c = e;
 614           if (d > c) c = d;
 615           s = r_ss[j].H;
 616           r_ss[j].H = c;
 617           r_ss[j].E = d;
 618         }
 619     }
 620   r_ss[N].E = r_ss[N].H;
 621
 622   midc = f_ss[0].H+r_ss[0].H;           /* Find optimal midpoint */
 623   midj = 0;
 624   type = 1;
 625   for (j = 0; j <= N; j++)
 626     if ((c = f_ss[j].H + r_ss[j].H) >= midc)
 627       if (c > midc || f_ss[j].H != f_ss[j].E && r_ss[j].H == r_ss[j].E)
 628         { midc = c;
 629           midj = j;
 630         }
 631   for (j = N; j >= 0; j--)
 632     if ((c = f_ss[j].E + r_ss[j].E + g) > midc)
 633       { midc = c;
 634         midj = j;
 635         type = 2;
 636       }
 637   }
 638
 639 /* Conquer: recursively around midpoint */
 640
 641   if (type == 1)
 642     { c1 = align(A,B,midi,midj,tb,-g,w,iw,g,h,f_str,0, sapp, last);
 643       c2 = align(A+midi,B+midj,M-midi,N-midj,-g,te,w,iw+midi,g,h,f_str,1,sapp,last);
 644     }
 645   else
 646     { align(A,B,midi-1,midj,tb,0,w,iw,g,h,f_str,2,sapp, last);
 647       DEL(2);
 648       align(A+midi+1,B+midj,M-midi-1,N-midj,0,te,w,iw+midi+1,g,h,f_str,3,sapp,last);
 649     }
 650   return midc;
 651 }
 652
 653 /* Interface and top level of comparator */
 654
 655 int ALIGN(const unsigned char *A, const unsigned char *B, int M, int N,
 656           int **W, int IW, int G, int H, int *S, int *NC,
 657           struct f_struct *f_str)
 658 {
 659   struct swstr *f_ss, *r_ss;
 660   int *sapp, last;
 661   int c, ck;
 662
 663   sapp = S;
 664   last = 0;
 665
 666    if ((f_ss = (struct swstr *) calloc (N+2, sizeof (struct swstr)))
 667        == NULL) {
 668      fprintf (stderr, "cannot allocate f_ss array %3d\n", N+2);
 669      exit (1);
 670    }
 671    f_ss++;
 672    f_str->f_ss = f_ss;
 673
 674    if ((r_ss = (struct swstr *) calloc (N+2, sizeof (struct swstr)))
 675        == NULL) {
 676      fprintf (stderr, "cannot allocate r_ss array %3d\n", N+2);
 677      exit (1);
 678    }
 679    r_ss++;
 680    f_str->r_ss = r_ss;
 681
 682   /*   print_seq_prof(A,M,W,IW); */
 683   c = align(A,B,M,N,-G,-G,W,IW,G,H,f_str,0,&sapp, &last);       /* OK, do it */
 684
 685   ck = CHECK_SCORE(A,B,M,N,S,W,IW,G,H,NC);
 686   if (c != ck) printf("Check_score error. %d != %d\n",c,ck);
 687
 688   f_ss--; r_ss--;
 689   free(r_ss); free(f_ss);
 690
 691   return c;
 692 }
 693
 694 /* Alignment display routine */
 695
 696 static char ALINE[51], BLINE[51], CLINE[51];
 697
 698 void DISPLAY(unsigned char *A, unsigned char *B, int M, int N,
 699             int *S, int AP, int BP, char *sq)
 700 { register char *a, *b, *c;
 701   register int   i,  j, op;
 702            int   lines, ap, bp;
 703
 704   i = j = op = lines = 0;
 705   ap = AP;
 706   bp = BP;
 707   a = ALINE;
 708   b = BLINE;
 709   c = CLINE;
 710   while (i < M || j < N)
 711     { if (op == 0 && *S == 0)
 712         { op = *S++;
 713           *a = sq[A[++i]];
 714           *b = sq[B[++j]];
 715           *c++ = (*a++ == *b++) ? '|' : ' ';
 716         }
 717       else
 718         { if (op == 0)
 719             op = *S++;
 720           if (op > 0)
 721             { *a++ = ' ';
 722               *b++ = sq[B[++j]];
 723               op--;
 724             }
 725           else
 726             { *a++ = sq[A[++i]];
 727               *b++ = ' ';
 728               op++;
 729             }
 730           *c++ = '-';
 731         }
 732       if (a >= ALINE+50 || i >= M && j >= N)
 733         { *a = *b = *c = '\0';
 734           printf("\n%5d ",50*lines++);
 735           for (b = ALINE+10; b <= a; b += 10)
 736             printf("    .    :");
 737           if (b <= a+5)
 738             printf("    .");
 739           printf("\n%5d %s\n      %s\n%5d %s\n",ap,ALINE,CLINE,bp,BLINE);
 740           ap = AP + i;
 741           bp = BP + j;
 742           a = ALINE;
 743           b = BLINE;
 744           c = CLINE;
 745         }
 746     }
 747 }
 748
 749 /* CHECK_SCORE - return the score of the alignment stored in S */
 750
 751 static int CHECK_SCORE(const unsigned char *A, const unsigned char *B,
 752                        int M, int N, int *S, int **w, int iw,
 753                        int g, int h, int *NC)
 754 {
 755   register int   i,  j, op, nc;
 756   int score;
 757
 758   /*  print_seq_prof(A,M,w,iw); */
 759
 760   score = i = j = op = nc = 0;
 761   while (i < M || j < N) {
 762     op = *S++;
 763     if (op == 0) {
 764       score = w[iw+i][B[++j]] + score;
 765       i++;
 766       nc++;
 767     }
 768     else if (op > 0) {
 769       score = score - (g+op*h);
 770       j += op;
 771       nc += op;
 772     } else {
 773       score = score - (g-op*h);
 774       i -= op;
 775       nc -= op;
 776     }
 777   }
 778   *NC = nc;
 779   return score;
 780 }
 781
 782 void
 783 pre_cons(const unsigned char *aa1, int n1, int frame, struct f_struct *f_str) {}
 784
 785 /* aln_func_vals - set up aln.qlfact, qlrev, llfact, llmult, frame, llrev */
 786 /* call from calcons, calc_id, calc_code */
 787 void
 788 aln_func_vals(int frame, struct a_struct *aln) {
 789
 790   aln->llfact = aln->llmult = aln->qlfact = 1;
 791   aln->qlrev = aln->llrev = 0;
 792   aln->frame = 0;
 793 }
 794
 795 /* 29-June-2003 this version has been modified to use pst.pam2p
 796    instead of pam2 to indicate similarity */
 797
 798 #include "a_mark.h"
 799
 800 int calcons(const unsigned char *aa0, int n0,
 801             const unsigned char *aa1, int n1,
 802             int *nc,
 803             struct a_struct *aln,
 804             struct a_res_str a_res,
 805             struct pstruct pst,
 806             char *seqc0, char *seqc1, char *seqca,
 807             struct f_struct *f_str)
 808 {
 809   int i0, i1;
 810   int op, lenc, nd, ns, itmp;
 811   char *sp0, *sp1, *spa, *sq;
 812   int *rp;
 813   int mins, smins;
 814
 815   if (pst.ext_sq_set) { sq = pst.sqx;}
 816   else {sq = pst.sq;}
 817
 818   aln->amin0 = a_res.min0;
 819   aln->amax0 = a_res.max0;
 820   aln->amin1 = a_res.min1;
 821   aln->amax1 = a_res.max1;
 822
 823   /* #define LFASTA */
 824 #ifndef LFASTA
 825   if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1)     /* will we show all the start ?*/
 826     if (a_res.min0>=a_res.min1) {              /* aa0 extends more to left */
 827       smins=0;
 828       if (aln->showall==1) mins=a_res.min0;
 829       else mins = min(a_res.min0,aln->llcntx);
 830       aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
 831       aln->smin0 = a_res.min0-mins;
 832       if ((mins-a_res.min1)>0) {
 833         memset(seqc1,' ',mins-a_res.min1);
 834         aancpy(seqc1+mins-a_res.min1,(char *)aa1,a_res.min1,pst);
 835         aln->smin1 = 0;
 836       }
 837       else {
 838         aancpy(seqc1,(char *)aa1+a_res.min1-mins,mins,pst);
 839         aln->smin1 = a_res.min1-mins;
 840       }
 841     }
 842     else {
 843       smins=0;
 844       if (aln->showall == 1) mins=a_res.min1;
 845       else mins = min(a_res.min1,aln->llcntx);
 846       aancpy(seqc1,(char *)(aa1+a_res.min1-mins),mins,pst);
 847       aln->smin1 = a_res.min1-mins;
 848       if ((mins-a_res.min0)>0) {
 849         memset(seqc0,' ',mins-a_res.min0);
 850         aancpy(seqc0+mins-a_res.min0,(char *)aa0,a_res.min0,pst);
 851         aln->smin0 = 0;
 852       }
 853       else {
 854         aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
 855         aln->smin0 = a_res.min0-mins;
 856       }
 857     }
 858   else {
 859     mins= min(aln->llcntx,min(a_res.min0,a_res.min1));
 860     smins=mins;
 861     aln->smin0=a_res.min0;
 862     aln->smin1=a_res.min1;
 863     aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
 864     aancpy(seqc1,(char *)aa1+a_res.min1-mins,mins,pst);
 865   }
 866 #else
 867   aln->smin0 = a_res.min0;
 868   aln->smin1 = a_res.min1;
 869   smins = mins = 0;
 870 #endif
 871
 872 /* now get the middle */
 873
 874   memset(seqca,M_BLANK,mins);
 875
 876   spa = seqca+mins;
 877   sp0 = seqc0+mins;
 878   sp1 = seqc1+mins;
 879   rp = a_res.res;
 880   lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs =op = 0;
 881   i0 = a_res.min0;
 882   i1 = a_res.min1;
 883
 884   while (i0 < a_res.max0 || i1 < a_res.max1) {
 885     if (op == 0 && *rp == 0) {
 886       op = *rp++;
 887       lenc++;
 888       if ((itmp=f_str->pam2p[0][i0][aa1[i1]])<0) { *spa = M_NEG; }
 889       else if (itmp == 0) { *spa = M_ZERO;}
 890       else {*spa = M_POS;}
 891       if (*spa == M_POS || *spa==M_ZERO) aln->nsim++;
 892
 893       *sp0 = sq[aa0[i0++]];
 894       *sp1 = sq[aa1[i1++]];
 895
 896       if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
 897       else if (pst.dnaseq==1 && ((*sp0 == 'T' && *sp1 == 'U') ||
 898                                  (*sp0=='U' && *sp1=='T'))) {
 899         aln->nident++; *spa=M_IDENT;
 900       }
 901
 902       sp0++; sp1++; spa++;
 903     }
 904     else {
 905       if (op==0) op = *rp++;
 906       if (op>0) {
 907         *sp0++ = '-';
 908         *sp1++ = sq[aa1[i1++]];
 909         *spa++ = M_DEL;
 910         op--;
 911         lenc++;
 912         aln->ngap_q++;
 913       }
 914       else {
 915         *sp0++ = sq[aa0[i0++]];
 916         *sp1++ = '-';
 917         *spa++ = M_DEL;
 918         op++;
 919         lenc++;
 920         aln->ngap_l++;
 921       }
 922     }
 923   }
 924
 925   *nc = lenc;
 926   *spa = '\0';
 927 /*      now we have the middle, get the right end */
 928
 929 #ifndef LFASTA
 930   /* how much extra to show at end ? */
 931   if (!aln->llcntx_flg) {
 932     ns = mins + lenc + aln->llen;       /* show an extra line? */
 933     ns -= (itmp = ns %aln->llen);       /* itmp = left over on last line */
 934     if (itmp>aln->llen/2) ns += aln->llen;  /* more than 1/2 , use another*/
 935     nd = ns - (mins+lenc);              /* this much extra */
 936   }
 937   else nd = aln->llcntx;
 938
 939   if (nd > max(n0-a_res.max0,n1-a_res.max1))
 940     nd = max(n0-a_res.max0,n1-a_res.max1);
 941
 942   if (aln->showall==1) {
 943     nd = max(n0-a_res.max0,n1-a_res.max1);      /* reset for showall=1 */
 944     /* get right end */
 945     aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-a_res.max0,pst);
 946     aancpy(seqc1+mins+lenc,(char *)aa1+a_res.max1,n1-a_res.max1,pst);
 947     /* fill with blanks - this is required to use one 'nc' */
 948     memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
 949     memset(seqc1+mins+lenc+n1-a_res.max1,' ',nd-(n1-a_res.max1));
 950   }
 951   else {
 952      if ((nd-(n0-a_res.max0))>0) {
 953        aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-a_res.max0,pst);
 954        memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
 955      }
 956      else aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,nd,pst);
 957
 958      if ((nd-(n1-a_res.max1))>0) {
 959        aancpy(seqc1+mins+lenc,(char *)aa1+a_res.max1,n1-a_res.max1,pst);
 960        memset(seqc1+mins+lenc+n1-a_res.max1,' ',nd-(n1-a_res.max1));
 961      }
 962      else aancpy(seqc1+mins+lenc,(char *)aa1+a_res.max1,nd,pst);
 963  }
 964
 965 #else   /* LFASTA */
 966   nd = 0;
 967 #endif
 968   /* #undef LFASTA */
 969   return mins+lenc+nd;
 970 }
 971
 972 int calcons_a(const unsigned char *aa0, unsigned char *aa0a, int n0,
 973               const unsigned char *aa1, int n1,
 974               int *nc,
 975               struct a_struct *aln,
 976               struct a_res_str a_res,
 977               struct pstruct pst,
 978               char *seqc0, char *seqc0a, char *seqc1, char *seqca,
 979               char *ann_arr, struct f_struct *f_str)
 980 {
 981   int i0, i1;
 982   int op, lenc, nd, ns, itmp;
 983   char *sp0, *sp0a, *sp1, *spa, *sq;
 984   int *rp;
 985   int mins, smins;
 986
 987   if (pst.ext_sq_set) {sq = pst.sqx;}
 988   else {sq = pst.sq;}
 989
 990   aln->amin0 = a_res.min0;
 991   aln->amax0 = a_res.max0;
 992   aln->amin1 = a_res.min1;
 993   aln->amax1 = a_res.max1;
 994
 995   /* first fill in the ends */
 996
 997   /* #define LFASTA */
 998 #ifndef LFASTA
 999   if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1)     /* will we show all the start ?*/
1000     if (a_res.min0>=a_res.min1) {              /* aa0 extends more to left */
1001       smins=0;
1002       if (aln->showall==1) mins=a_res.min0;
1003       else mins = min(a_res.min0,aln->llcntx);
1004       aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
1005       aln->smin0 = a_res.min0-mins;
1006       if ((mins-a_res.min1)>0) {
1007         memset(seqc1,' ',mins-a_res.min1);
1008         aancpy(seqc1+mins-a_res.min1,(char *)aa1,a_res.min1,pst);
1009         aln->smin1 = 0;
1010       }
1011       else {
1012         aancpy(seqc1,(char *)aa1+a_res.min1-mins,mins,pst);
1013         aln->smin1 = a_res.min1-mins;
1014       }
1015     }
1016     else {
1017       smins=0;
1018       if (aln->showall == 1) mins=a_res.min1;
1019       else mins = min(a_res.min1,aln->llcntx);
1020       aancpy(seqc1,(char *)(aa1+a_res.min1-mins),mins,pst);
1021       aln->smin1 = a_res.min1-mins;
1022       if ((mins-a_res.min0)>0) {
1023         memset(seqc0,' ',mins-a_res.min0);
1024         aancpy(seqc0+mins-a_res.min0,(char *)aa0,a_res.min0,pst);
1025         aln->smin0 = 0;
1026       }
1027       else {
1028         aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
1029         aln->smin0 = a_res.min0-mins;
1030       }
1031     }
1032   else {
1033     mins= min(aln->llcntx,min(a_res.min0,a_res.min1));
1034     smins=mins;
1035     aln->smin0=a_res.min0;
1036     aln->smin1=a_res.min1;
1037     aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
1038     aancpy(seqc1,(char *)aa1+a_res.min1-mins,mins,pst);
1039   }
1040 #else
1041   aln->smin0 = a_res.min0;
1042   aln->smin1 = a_res.min1;
1043   smins = mins = 0;
1044 #endif
1045
1046 /* now get the middle */
1047
1048   memset(seqca,M_BLANK,mins);
1049   memset(seqc0a,' ',mins);
1050
1051   spa = seqca+mins;
1052   sp0 = seqc0+mins;
1053   sp0a = seqc0a+mins;
1054   sp1 = seqc1+mins;
1055   rp = a_res.res;
1056   lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs =op = 0;
1057   i0 = a_res.min0;
1058   i1 = a_res.min1;
1059
1060   while (i0 < a_res.max0 || i1 < a_res.max1) {
1061     if (op == 0 && *rp == 0) {
1062       op = *rp++;
1063       lenc++;
1064       if ((itmp=f_str->pam2p[0][i0][aa1[i1]])<0) { *spa = M_NEG; }
1065       else if (itmp == 0) { *spa = M_ZERO;}
1066       else {*spa = M_POS;}
1067       if (*spa == M_POS || *spa==M_ZERO) aln->nsim++;
1068
1069       *sp0a++ = ann_arr[aa0a[i0]];
1070       *sp0 = sq[aa0[i0++]];
1071       *sp1 = sq[aa1[i1++]];
1072
1073       if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
1074       else if (pst.dnaseq==1 && ((*sp0 == 'T' && *sp1 == 'U') ||
1075                                  (*sp0=='U' && *sp1=='T'))) {
1076         aln->nident++; *spa=M_IDENT;
1077       }
1078
1079       sp0++; sp1++; spa++;
1080     }
1081     else {
1082       if (op==0) op = *rp++;
1083       if (op>0) {
1084         *sp0++ = '-';
1085         *sp0a++ = ' ';
1086         *sp1++ = sq[aa1[i1++]];
1087         *spa++ = M_DEL;
1088         op--;
1089         lenc++;
1090         aln->ngap_q++;
1091       }
1092       else {
1093         *sp0a++ = ann_arr[aa0a[i0]];
1094         *sp0++ = sq[aa0[i0++]];
1095         *sp1++ = '-';
1096         *spa++ = M_DEL;
1097         op++;
1098         lenc++;
1099         aln->ngap_l++;
1100       }
1101     }
1102   }
1103
1104   *nc = lenc;
1105   *spa = '\0';
1106 /*      now we have the middle, get the right end */
1107
1108 #ifndef LFASTA
1109   /* how much extra to show at end ? */
1110   if (!aln->llcntx_flg) {
1111     ns = mins + lenc + aln->llen;       /* show an extra line? */
1112     ns -= (itmp = ns %aln->llen);       /* itmp = left over on last line */
1113     if (itmp>aln->llen/2) ns += aln->llen;  /* more than 1/2 , use another*/
1114     nd = ns - (mins+lenc);              /* this much extra */
1115   }
1116   else nd = aln->llcntx;
1117
1118   if (nd > max(n0-a_res.max0,n1-a_res.max1))
1119     nd = max(n0-a_res.max0,n1-a_res.max1);
1120
1121   if (aln->showall==1) {
1122     nd = max(n0-a_res.max0,n1-a_res.max1);      /* reset for showall=1 */
1123     /* get right end */
1124     aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-a_res.max0,pst);
1125     aancpy(seqc1+mins+lenc,(char *)aa1+a_res.max1,n1-a_res.max1,pst);
1126     /* fill with blanks - this is required to use one 'nc' */
1127     memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
1128     memset(seqc1+mins+lenc+n1-a_res.max1,' ',nd-(n1-a_res.max1));
1129   }
1130   else {
1131      if ((nd-(n0-a_res.max0))>0) {
1132        aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-a_res.max0,pst);
1133        memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
1134      }
1135      else aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,nd,pst);
1136
1137      if ((nd-(n1-a_res.max1))>0) {
1138        aancpy(seqc1+mins+lenc,(char *)aa1+a_res.max1,n1-a_res.max1,pst);
1139        memset(seqc1+mins+lenc+n1-a_res.max1,' ',nd-(n1-a_res.max1));
1140      }
1141      else aancpy(seqc1+mins+lenc,(char *)aa1+a_res.max1,nd,pst);
1142  }
1143
1144 #else   /* LFASTA */
1145   nd = 0;
1146 #endif
1147   /* #undef LFASTA */
1148   return mins+lenc+nd;
1149 }
1150
1151 static void
1152 update_code(char *al_str, int al_str_max, int op, int op_cnt);
1153
1154 /* build an array of match/ins/del - length strings */
1155 int calc_code(const unsigned char *aa0, const int n0,
1156               const unsigned char *aa1, const int n1,
1157               struct a_struct *aln,
1158               struct a_res_str a_res,
1159               struct pstruct pst,
1160               char *al_str, int al_str_n, struct f_struct *f_str)
1161 {
1162   int i0, i1, nn1;
1163   int op, lenc, nd, ns, itmp;
1164   int p_op, op_cnt;
1165   const unsigned char *aa1p;
1166   char tmp_cnt[20];
1167   char sp0, sp1, *sq;
1168   int *rp;
1169   int mins, smins;
1170
1171   if (pst.ext_sq_set) {
1172     sq = pst.sqx;
1173   }
1174   else {
1175     sq = pst.sq;
1176   }
1177
1178 #ifndef TFASTA
1179   aa1p = aa1;
1180   nn1 = n1;
1181 #else
1182   aa1p = f_str->aa1x;
1183   nn1 = f_str->n10;
1184 #endif
1185
1186   aln->amin0 = a_res.min0;
1187   aln->amax0 = a_res.max0;
1188   aln->amin1 = a_res.min1;
1189   aln->amax1 = a_res.max1;
1190
1191   rp = a_res.res;
1192   lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = p_op = 0;
1193   op_cnt = 0;
1194
1195   i0 = a_res.min0;
1196   i1 = a_res.min1;
1197   tmp_cnt[0]='\0';
1198
1199   while (i0 < a_res.max0 || i1 < a_res.max1) {
1200     if (op == 0 && *rp == 0) {
1201
1202       if (pst.pam2[0][aa0[i0]][aa1p[i1]]>=0) { aln->nsim++;}
1203
1204       sp0 = sq[aa0[i0++]];
1205       sp1 = sq[aa1p[i1++]];
1206
1207       if (p_op == 0 || p_op==3) {
1208         if (sp0 != '*' && sp1 != '*') {
1209           if (p_op == 3) {
1210             update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
1211             op_cnt = 1; p_op = 0;
1212           }
1213           else {op_cnt++;}
1214         }
1215         else {
1216           update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
1217           op_cnt = 1; p_op = 3;
1218         }
1219       }
1220       else {
1221         update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
1222         op_cnt = 1; p_op = 0;
1223       }
1224
1225       op = *rp++;
1226       lenc++;
1227
1228       if (toupper(sp0) == toupper(sp1)) aln->nident++;
1229       else if (pst.dnaseq==1) {
1230         if ((toupper(sp0) == 'T' && toupper(sp1) == 'U') ||
1231             (toupper(sp0)=='U' && toupper(sp1)=='T')) aln->nident++;
1232         else if (toupper(sp0) == 'N') aln->ngap_q++;
1233         else if (toupper(sp1) == 'N') aln->ngap_l++;
1234       }
1235     }
1236     else {
1237       if (op==0) op = *rp++;
1238       if (op>0) {
1239         if (p_op == 1) { op_cnt++;}
1240         else {
1241           update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
1242           op_cnt = 1; p_op = 1;
1243         }
1244         op--; lenc++; i1++; aln->ngap_q++;
1245       }
1246       else {
1247         if (p_op == 2) { op_cnt++;}
1248         else {
1249           update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
1250           op_cnt = 1; p_op = 2;
1251         }
1252         op++; lenc++; i0++; aln->ngap_l++;
1253       }
1254     }
1255   }
1256   update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
1257
1258   return lenc;
1259 }
1260
1261 static void
1262 update_code(char *al_str, int al_str_max, int op, int op_cnt) {
1263
1264   char op_char[5]={"=-+*"};
1265   char tmp_cnt[20];
1266
1267   sprintf(tmp_cnt,"%c%d",op_char[op],op_cnt);
1268   strncat(al_str,tmp_cnt,al_str_max);
1269 }
1270
1271 int calc_id(const unsigned char *aa0, const int n0,
1272             const unsigned char *aa1, const int n1,
1273             struct a_struct *aln,
1274             struct a_res_str a_res,
1275             struct pstruct pst,
1276             struct f_struct *f_str)
1277 {
1278   int i0, i1, nn1, n_id;
1279   int op, lenc, nd, ns, itmp;
1280   int sp0, sp1;
1281   const unsigned char *aa1p;
1282   int *rp;
1283   char *sq;
1284
1285   if (pst.ext_sq_set) {
1286     sq = pst.sqx;
1287   }
1288   else {
1289     sq = pst.sq;
1290   }
1291
1292 #ifndef TFASTA
1293   aa1p = aa1;
1294   nn1 = n1;
1295 #else
1296   aa1p = f_str->aa1x;
1297   nn1 = f_str->n10;
1298 #endif
1299
1300   rp = a_res.res;
1301   lenc = n_id = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = 0;
1302   i0 = a_res.min0;
1303   i1 = a_res.min1;
1304
1305   while (i0 < a_res.max0 || i1 < a_res.max1) {
1306     if (op == 0 && *rp == 0) {
1307       op = *rp++;
1308       lenc++;
1309       if (pst.pam2[0][aa0[i0]][aa1p[i1]]>=0) { aln->nsim++;}
1310
1311       sp0 = sq[aa0[i0++]];
1312       sp1 = sq[aa1p[i1++]];
1313       if (toupper(sp0) == toupper(sp1)) n_id++;
1314       else if (pst.dnaseq==1 &&
1315                ((sp0=='T' && sp1== 'U')||(sp0=='U' && sp1=='T'))) n_id++;
1316     }
1317     else {
1318       if (op==0) op = *rp++;
1319       if (op>0) {op--; lenc++; i1++; aln->ngap_q++; }
1320       else {op++; lenc++; i0++; aln->ngap_l++; }
1321     }
1322   }
1323   aln->nident = n_id;
1324   return lenc;
1325 }
1326
1327 #ifdef PCOMPLIB
1328 #include "p_mw.h"
1329 void
1330 update_params(struct qmng_str *qm_msg, struct pstruct *ppst)
1331 {
1332   ppst->n0 = qm_msg->n0;
1333 }
1334 #endif