开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
greedy_align.c
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：25k
源码类别：

生物技术

开发平台：

C/C++

greedy_align.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: greedy_align.c,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:08:02  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.19
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: greedy_align.c,v 1000.2 2004/06/01 18:08:02 gouriano Exp $
  11.  * ===========================================================================
  12.  *
  13.  *                            PUBLIC DOMAIN NOTICE
  14.  *               National Center for Biotechnology Information
  15.  *
  16.  *  This software/database is a "United States Government Work" under the
  17.  *  terms of the United States Copyright Act.  It was written as part of
  18.  *  the author's official duties as a United States Government employee and
  19.  *  thus cannot be copyrighted.  This software/database is freely available
  20.  *  to the public for use. The National Library of Medicine and the U.S.
  21.  *  Government have not placed any restriction on its use or reproduction.
  22.  *
  23.  *  Although all reasonable efforts have been taken to ensure the accuracy
  24.  *  and reliability of the software and data, the NLM and the U.S.
  25.  *  Government do not and cannot warrant the performance or results that
  26.  *  may be obtained by using this software or data. The NLM and the U.S.
  27.  *  Government disclaim all warranties, express or implied, including
  28.  *  warranties of performance, merchantability or fitness for any particular
  29.  *  purpose.
  30.  *
  31.  *  Please cite the author in any work or product based on this material.
  32.  *
  33.  * ===========================================================================
  34.  *
  35.  * Author: Webb Miller and Co. Adopted for NCBI libraries by Sergey Shavirin
  36.  *
  37.  * Initial Creation Date: 10/27/1999
  38.  *
  39.  */
  41. /** @file greedy_align.c
  42.  * Greedy gapped alignment functions
  43.  */
  45. static char const rcsid[] = 
  46.     "$Id: greedy_align.c,v 1000.2 2004/06/01 18:08:02 gouriano Exp $";
  48. #include <algo/blast/core/greedy_align.h>
  49. #include <algo/blast/core/blast_util.h> /* for NCBI2NA_UNPACK_BASE macros */
  51. enum {
  52.     EDIT_OP_MASK = 0x3,
  53.     EDIT_OP_ERR  = 0x0,
  54.     EDIT_OP_INS  = 0x1,
  55.     EDIT_OP_DEL  = 0x2,
  56.     EDIT_OP_REP  = 0x3
  57. };
  59. enum {         /* half of the (fixed) match score */
  60.     ERROR_FRACTION=2,  /* 1/this */
  61.     MAX_SPACE=1000000,
  62.     sC = 0, sI = 1, sD = 2, LARGE=100000000
  63. };
  66. /* -------- From original file edit.c ------------- */
  68. static Uint4 edit_val_get(edit_op_t op)
  69. {
  70.     return op >> 2;
  71. }
  72. static Uint4 edit_opc_get(edit_op_t op)
  73. {
  74.     return op & EDIT_OP_MASK;
  75. }
  77. static edit_op_t edit_op_cons(Uint4 op, Uint4 val)
  78. {
  79.     return (val << 2) | (op & EDIT_OP_MASK);
  80. }
  82. static edit_op_t edit_op_inc(edit_op_t op, Uint4 n)
  83. {
  84.     return edit_op_cons(edit_opc_get(op), edit_val_get(op) + n);
  85. }
  87. static edit_op_t edit_op_inc_last(MBGapEditScript *es, Uint4 n)
  88. {
  89.     edit_op_t *last;
  90.     ASSERT (es->num > 0);
  91.     last = &(es->op[es->num-1]);
  92.     *last = edit_op_inc(*last, n);
  93.     return *last;
  94. }
  96. static Int4 edit_script_ready(MBGapEditScript *es, Uint4 n)
  97. {
  98.     edit_op_t *p;
  99.     Uint4 m = n + n/2;
  100.     
  101.     if (es->size <= n) {
  102.         p = realloc(es->op, m*sizeof(edit_op_t));
  103.         if (p == 0) {
  104.             return 0;
  105.         } else {
  106.             es->op = p;
  107.             es->size = m;
  108.         }
  109.     }
  110.     return 1;
  111. }
  113. static Int4 edit_script_readyplus(MBGapEditScript *es, Uint4 n)
  114. {
  115.     if (es->size - es->num <= n)
  116.         return edit_script_ready(es, n + es->num);
  117.     return 1;
  118. }
  120. static Int4 edit_script_put(MBGapEditScript *es, Uint4 op, Uint4 n)
  121. {
  122.     if (!edit_script_readyplus(es, 2))
  123.         return 0;
  124.     es->last = op;
  125.     ASSERT(op != 0);
  126.     es->op[es->num] = edit_op_cons(op, n);
  127.     es->num += 1;
  128.     es->op[es->num] = 0; /* sentinal */
  130.     return 1;
  131. }
  133. static MBGapEditScript *edit_script_init(MBGapEditScript *es)
  134. {
  135. es->op = 0;
  136. es->size = es->num = 0;
  137. es->last = 0;
  138. edit_script_ready(es, 8);
  139. return es;
  140. }
  141. static edit_op_t *edit_script_first(MBGapEditScript *es)
  142. {
  143.     return es->num > 0 ? &es->op[0] : 0;
  144. }
  146. static edit_op_t *edit_script_next(MBGapEditScript *es, edit_op_t *op)
  147. {
  148.     /* XXX - assumes flat address space */
  149.     if (&es->op[0] <= op && op < &es->op[es->num-1])
  150.         return op+1;
  151.     else
  152.         return 0;
  153. }
  154. static Int4 edit_script_more(MBGapEditScript *data, Uint4 op, Uint4 k)
  155. {
  156.     if (op == EDIT_OP_ERR) {
  157. #ifdef ERR_POST_EX_DEFINED
  158.         ErrPostEx(SEV_FATAL, 1, 0, 
  159.                   "edit_script_more: bad opcode %d:%d", op, k);
  160. #endif
  161.         return -1;
  162.     }
  163.     
  164.     if (edit_opc_get(data->last) == op)
  165.         edit_op_inc_last(data, k);
  166.     else
  167.         edit_script_put(data, op, k);
  169.     return 0;
  170. }
  172. MBGapEditScript *MBGapEditScriptAppend(MBGapEditScript *es, MBGapEditScript *et)
  173. {
  174.     edit_op_t *op;
  175.     
  176.     for (op = edit_script_first(et); op; op = edit_script_next(et, op))
  177.         edit_script_more(es, edit_opc_get(*op), edit_val_get(*op));
  179.     return es;
  180. }
  182. MBGapEditScript *MBGapEditScriptNew(void)
  183. {
  184.     MBGapEditScript *es = calloc(1, sizeof(*es));
  185.     if (!es)
  186.         return 0;
  188.     return edit_script_init(es);
  189. }
  191. MBGapEditScript *MBGapEditScriptFree(MBGapEditScript *es)
  192. {
  193.     if (es) {
  194.         if (es->op)
  195.             sfree(es->op);
  196.         memset(es, 0, sizeof(*es));
  197.         sfree(es);
  198.     }
  199.     return 0;
  200. }
  202. static Int4 edit_script_del(MBGapEditScript *data, Uint4 k)
  203. {
  204.     return edit_script_more(data, EDIT_OP_DEL, k);
  205. }
  207. static Int4 edit_script_ins(MBGapEditScript *data, Uint4 k)
  208. {
  209.     return edit_script_more(data, EDIT_OP_INS, k);
  210. }
  211. static Int4 edit_script_rep(MBGapEditScript *data, Uint4 k)
  212. {
  213.     return edit_script_more(data, EDIT_OP_REP, k);
  214. }
  216. static MBGapEditScript *edit_script_reverse_inplace(MBGapEditScript *es)
  217. {
  218.     Uint4 i;
  219.     const Uint4 num = es->num;
  220.     const Uint4 mid = num/2;
  221.     const Uint4 end = num-1;
  222.     
  223.     for (i = 0; i < mid; ++i) {
  224.         const edit_op_t t = es->op[i];
  225.         es->op[i] = es->op[end-i];
  226.         es->op[end-i] = t;
  227.     }
  228.     return es;
  229. }
  231. SMBSpace* MBSpaceNew()
  232. {
  233.     SMBSpace* p;
  234.     Int4 amount;
  235.     
  236.     p = (SMBSpace*) malloc(sizeof(SMBSpace));
  237.     amount = MAX_SPACE;
  238.     p->space_array = (SThreeVal*) malloc(sizeof(SThreeVal)*amount);
  239.     if (p->space_array == NULL) {
  240.        sfree(p);
  241.        return NULL;
  242.     }
  243.     p->used = 0; 
  244.     p->size = amount;
  245.     p->next = NULL;
  247.     return p;
  248. }
  250. static void refresh_mb_space(SMBSpace* sp)
  251. {
  252.    while (sp) {
  253.       sp->used = 0;
  254.       sp = sp->next;
  255.    }
  256. }
  258. void MBSpaceFree(SMBSpace* sp)
  259. {
  260.    SMBSpace* next_sp;
  262.    while (sp) {
  263.       next_sp = sp->next;
  264.       sfree(sp->space_array);
  265.       sfree(sp);
  266.       sp = next_sp;
  267.    }
  268. }
  270. static SThreeVal* get_mb_space(SMBSpace* S, Int4 amount)
  271. {
  272.     SThreeVal* s;
  273.     if (amount < 0) 
  274.         return NULL;  
  275.     
  276.     while (S->used+amount > S->size) {
  277.        if (S->next == NULL)
  278.           if ((S->next = MBSpaceNew()) == NULL) {
  279. #ifdef ERR_POST_EX_DEFINED
  280.      ErrPostEx(SEV_WARNING, 0, 0, "Cannot get new space for greedy extension");
  281. #endif
  282.      return NULL;
  283.           }
  284.        S = S->next;
  285.     }
  287.     s = S->space_array+S->used;
  288.     S->used += amount;
  289.     
  290.     return s;
  291. }
  292. /* ----- */
  294. static Int4 gcd(Int4 a, Int4 b)
  295. {
  296.     Int4 c;
  297.     if (a < b) {
  298.         c = a;
  299.         a = b; b = c;
  300.     }
  301.     while ((c = a % b) != 0) {
  302.         a = b; b = c;
  303.     }
  305.     return b;
  306. }
  309. static Int4 gdb3(Int4* a, Int4* b, Int4* c)
  310. {
  311.     Int4 g;
  312.     if (*b == 0) g = gcd(*a, *c);
  313.     else g = gcd(*a, gcd(*b, *c));
  314.     if (g > 1) {
  315.         *a /= g;
  316.         *b /= g;
  317.         *c /= g;
  318.     }
  319.     return g;
  320. }
  322. static Int4 get_lastC(SThreeVal** flast_d, Int4* lower, Int4* upper, 
  323.                       Int4* d, Int4 diag, Int4 Mis_cost, Int4* row1)
  324. {
  325.     Int4 row;
  326.     
  327.     if (diag >= lower[(*d)-Mis_cost] && diag <= upper[(*d)-Mis_cost]) {
  328.         row = flast_d[(*d)-Mis_cost][diag].C;
  329.         if (row >= MAX(flast_d[*d][diag].I, flast_d[*d][diag].D)) {
  330.             *d = *d-Mis_cost;
  331.             *row1 = row;
  332.             return sC;
  333.         }
  334.     }
  335.     if (flast_d[*d][diag].I > flast_d[*d][diag].D) {
  336.         *row1 = flast_d[*d][diag].I;
  337.         return sI;
  338.     } else {
  339.         *row1 = flast_d[*d][diag].D;
  340.         return sD;
  341.     }
  342. }
  344. static Int4 get_last_ID(SThreeVal** flast_d, Int4* lower, Int4* upper, 
  345.                         Int4* d, Int4 diag, Int4 GO_cost, 
  346.                         Int4 GE_cost, Int4 IorD)
  347. {
  348.     Int4 ndiag; 
  349.     Int4 row;
  351.     if (IorD == sI) { ndiag = diag -1;} else ndiag = diag+1;
  352.     if (ndiag >= lower[(*d)-GE_cost] && ndiag <=upper[(*d)-GE_cost]) 
  353.         row = (IorD == sI)? flast_d[(*d)-GE_cost][ndiag].I: flast_d[(*d)-GE_cost][ndiag].D;
  354.     else row = -100;
  355.     if (ndiag >= lower[(*d)-GO_cost-GE_cost] && ndiag <= upper[(*d)-GO_cost-GE_cost] && row < flast_d[(*d)-GO_cost-GE_cost][ndiag].C) {
  356.         *d = (*d)-GO_cost-GE_cost;
  357.         return sC;
  358.     }
  359.     *d = (*d)-GE_cost;
  360.     return IorD;
  361. }
  363. static Int4 get_lastI(SThreeVal** flast_d, Int4* lower, Int4* upper, 
  364.                       Int4* d, Int4 diag, Int4 GO_cost, Int4 GE_cost)
  365. {
  366.     return get_last_ID(flast_d, lower, upper, d, diag, GO_cost, GE_cost, sI);
  367. }
  370. static int get_lastD(SThreeVal** flast_d, Int4* lower, Int4* upper, 
  371.                      Int4* d, Int4 diag, Int4 GO_cost, Int4 GE_cost)
  372. {
  373.     return get_last_ID(flast_d, lower, upper, d, diag, GO_cost, GE_cost, sD);
  374. }
  376. /* --- From file align.c --- */
  377. /* ----- */
  379. static Int4 get_last(Int4 **flast_d, Int4 d, Int4 diag, Int4 *row1)
  380. {
  381.     if (flast_d[d-1][diag-1] > MAX(flast_d[d-1][diag], flast_d[d-1][diag+1])) {
  382.         *row1 = flast_d[d-1][diag-1];
  383.         return diag-1;
  384.     } 
  385.     if (flast_d[d-1][diag] > flast_d[d-1][diag+1]) {
  386.         *row1 = flast_d[d-1][diag];
  387.         return diag;
  388.     }
  389.     *row1 = flast_d[d-1][diag+1];
  390.     return diag+1;
  391. }
  393. /*
  394. Version to search a (possibly) packed nucl. sequence against
  395. an unpacked sequence.  s2 is the packed nucl. sequence.
  396. s1 can be packed or unpacked. If rem == 4, then s1 is unpacked.
  397. len2 corresponds to the unpacked (true) length.
  399.  * Basic O(ND) time, O(N) space, alignment function. 
  400.  * Parameters:
  401.  *   s1, len1        - first sequence and its length
  402.  *   s2, len2        - second sequence and its length
  403.  *   reverse         - direction of alignment
  404.  *   xdrop_threshold - 
  405.  *   mismatch_cost   -
  406.  *   e1, e2          - endpoint of the computed alignment
  407.  *   edit_script     -
  408.  *   rem             - offset shift within a packed sequence
  409.  */
  410. Int4 BLAST_GreedyAlign(const Uint1* s1, Int4 len1,
  411.   const Uint1* s2, Int4 len2,
  412.   Boolean reverse, Int4 xdrop_threshold, 
  413.   Int4 match_cost, Int4 mismatch_cost,
  414.   Int4* e1, Int4* e2, 
  415.   SGreedyAlignMem* gamp, MBGapEditScript *S,
  416.                           Uint1 rem)
  417. {
  418. #define ICEIL(x,y) ((((x)-1)/(y))+1)
  419.     Int4 col, /* column number */
  420.         d, /* current distance */
  421.         k, /* current diagonal */
  422.         flower, fupper,            /* boundaries for searching diagonals */
  423.         row,         /* row number */
  424.         MAX_D,  /* maximum cost */
  425.         ORIGIN,
  426.         return_val = 0;
  427.     Int4** flast_d = gamp->flast_d; /* rows containing the last d */
  428.     Int4* max_row; /* reached for cost d=0, ... len1.  */
  429.     
  430.     Int4 X_pen = xdrop_threshold;
  431.     Int4 M_half = match_cost/2;
  432.     Int4 Op_cost = mismatch_cost + M_half*2;
  433.     Int4 D_diff = ICEIL(X_pen+M_half, Op_cost);
  434.     
  435.     Int4 x, cur_max, b_diag = 0, best_diag = INT4_MAX/2;
  436.     Int4* max_row_free = gamp->max_row_free;
  437.     char nlower = 0, nupper = 0;
  438.     SMBSpace* space = gamp->space;
  439.     Int4 max_len = len2;
  440.  
  441.     MAX_D = (Int4) (len1/ERROR_FRACTION + 1);
  442.     ORIGIN = MAX_D + 2;
  443.     *e1 = *e2 = 0;
  444.     
  445.     if (reverse) {
  446.        if (!(rem & 4)) {
  447.           for (row = 0; row < len2 && row < len1 && 
  448.                   (s2[len2-1-row] ==
  449.                    NCBI2NA_UNPACK_BASE(s1[(len1-1-row)/4], 
  450.                                         3-(len1-1-row)%4)); 
  451.                row++)
  452.              /*empty*/ ;
  453.        } else {
  454.           for (row = 0; row < len2 && row < len1 && (s2[len2-1-row] == s1[len1-1-row]); row++)
  455.              /*empty*/ ;
  456.        }
  457.     } else {
  458.        if (!(rem & 4)) {
  459.           for (row = 0; row < len2 && row < len1 && 
  460.                   (s2[row] == 
  461.                    NCBI2NA_UNPACK_BASE(s1[(row+rem)/4], 
  462.                                         3-(row+rem)%4)); 
  463.                row++)
  464.              /*empty*/ ;
  465.        } else {
  466.           for (row = 0; row < len2 && row < len1 && (s2[row] == s1[row]); row++)
  467.              /*empty*/ ;
  468.        }
  469.     }
  470.     *e1 = row;
  471.     *e2 = row;
  472.     if (row == len1) {
  473.         if (S != NULL)
  474.             edit_script_rep(S, row);
  475. /* hit last row; stop search */
  476. return 0;
  477.     }
  478.     if (S==NULL) {
  479.        space = 0;
  480.     } else if (!space) {
  481.        gamp->space = space = MBSpaceNew();
  482.     } else { 
  483.         refresh_mb_space(space);
  484.     }
  485.     
  486.     max_row = max_row_free + D_diff;
  487.     for (k = 0; k < D_diff; k++)
  488. max_row_free[k] = 0;
  489.     
  490.     flast_d[0][ORIGIN] = row;
  491.     max_row[0] = (row + row)*M_half;
  492.     
  493.     flower = ORIGIN - 1;
  494.     fupper = ORIGIN + 1;
  496.     d = 1;
  497.     while (d <= MAX_D) {
  498. Int4 fl0, fu0;
  499. flast_d[d - 1][flower - 1] = flast_d[d - 1][flower] = -1;
  500. flast_d[d - 1][fupper] = flast_d[d - 1][fupper + 1] = -1;
  501. x = max_row[d - D_diff] + Op_cost * d - X_pen;
  502. x = ICEIL(x, M_half);
  503. cur_max = 0;
  504. fl0 = flower;
  505. fu0 = fupper;
  506. for (k = fl0; k <= fu0; k++) {
  507.     row = MAX(flast_d[d - 1][k + 1], flast_d[d - 1][k]) + 1;
  508.     row = MAX(row, flast_d[d - 1][k - 1]);
  509.     col = row + k - ORIGIN;
  510.     if (row + col >= x)
  511. fupper = k;
  512.     else {
  513. if (k == flower)
  514.     flower++;
  515. else
  516.     flast_d[d][k] = -1;
  517. continue;
  518.     }
  519.             
  520.             if (row > max_len || row < 0) {
  521.                   flower = k+1; nlower = 1;
  522.             } else {
  523.                /* Slide down the diagonal. */
  524.                if (reverse) {
  525.                   if (!(rem & 4)) {
  526.                      while (row < len2 && col < len1 && s2[len2-1-row] == 
  527.                             NCBI2NA_UNPACK_BASE(s1[(len1-1-col)/4],
  528.                                                  3-(len1-1-col)%4)) {
  529.                         ++row;
  530.                         ++col;
  531.                      }
  532.                   } else {
  533.                      while (row < len2 && col < len1 && 
  534.                             s2[len2-1-row] == s1[len1-1-col]) {
  535.                         ++row;
  536.                         ++col;
  537.                      }
  538.                   }
  539.                } else {
  540.                   if (!(rem & 4)) {
  541.                      while (row < len2 && col < len1 && s2[row] == 
  542.                             NCBI2NA_UNPACK_BASE(s1[(col+rem)/4],
  543.                                                  3-(col+rem)%4)) {
  544.                         ++row;
  545.                         ++col;
  546.                      }
  547.                   } else {
  548.                      while (row < len2 && col < len1 && s2[row] == s1[col]) {
  549.                         ++row;
  550.                         ++col;
  551.                      }
  552.                   }
  553.                }
  554.             }
  555.     flast_d[d][k] = row;
  556.     if (row + col > cur_max) {
  557. cur_max = row + col;
  558. b_diag = k;
  559.     }
  560.     if (row == len2) {
  561. flower = k+1; nlower = 1;
  562.     }
  563.     if (col == len1) {
  564. fupper = k-1; nupper = 1;
  565.     }
  566. }
  567. k = cur_max*M_half - d * Op_cost;
  568. if (max_row[d - 1] < k) {
  569.     max_row[d] = k;
  570.     return_val = d;
  571.     best_diag = b_diag;
  572.     *e2 = flast_d[d][b_diag];
  573.     *e1 = (*e2)+b_diag-ORIGIN;
  574. } else {
  575.     max_row[d] = max_row[d - 1];
  576. }
  577. if (flower > fupper)
  578.     break;
  579. d++;
  580. if (!nlower) flower--; 
  581. if (!nupper) fupper++;
  582. if (S==NULL)
  583.    flast_d[d] = flast_d[d - 2];
  584. else {
  585.            /* space array consists of SThreeVal structures which are 
  586.               3 times larger than Int4, so divide requested amount by 3
  587.            */
  588.    flast_d[d] = (Int4*) get_mb_space(space, (fupper-flower+7)/3);
  589.    if (flast_d[d] != NULL)
  590.               flast_d[d] = flast_d[d] - flower + 2;
  591.            else
  592.       return return_val;
  593.         }
  594.     }
  595.     
  596.     if (S!=NULL) { /*trace back*/
  597.         Int4 row1, col1, diag1, diag;
  598.         d = return_val; diag = best_diag;
  599.         row = *e2; col = *e1;
  600.         while (d > 0) {
  601.             diag1 = get_last(flast_d, d, diag, &row1);
  602.             col1 = row1+diag1-ORIGIN;
  603.             if (diag1 == diag) {
  604.                 if (row-row1 > 0) edit_script_rep(S, row-row1);
  605.             } else if (diag1 < diag) {
  606.                 if (row-row1 > 0) edit_script_rep(S, row-row1);
  607.                 edit_script_ins(S,1);
  608.             } else {
  609.                 if (row-row1-1> 0) edit_script_rep(S, row-row1-1);
  610.                 edit_script_del(S, 1);
  611.             }
  612.             d--; diag = diag1; col = col1; row = row1;
  613.         }
  614.         edit_script_rep(S, flast_d[0][ORIGIN]);
  615.         if (!reverse) 
  616.             edit_script_reverse_inplace(S);
  617.     }
  618.     return return_val;
  619. }
  622. Int4 BLAST_AffineGreedyAlign (const Uint1* s1, Int4 len1, 
  623.  const Uint1* s2, Int4 len2,
  624.  Boolean reverse, Int4 xdrop_threshold, 
  625.  Int4 match_score, Int4 mismatch_score, 
  626.  Int4 gap_open, Int4 gap_extend,
  627.  Int4* e1, Int4* e2, 
  628.  SGreedyAlignMem* gamp, MBGapEditScript *S,
  629.                                  Uint1 rem)
  630. {
  631.     Int4 col, /* column number */
  632.         d, /* current distance */
  633.         k, /* current diagonal */
  634.         flower, fupper,         /* boundaries for searching diagonals */
  635.         row,         /* row number */
  636.         MAX_D,  /* maximum cost */
  637.         ORIGIN,
  638.         return_val = 0;
  639.     SThreeVal** flast_d; /* rows containing the last d */
  640.     Int4* max_row_free = gamp->max_row_free;
  641.     Int4* max_row; /* reached for cost d=0, ... len1.  */
  642.     Int4 Mis_cost, GO_cost, GE_cost;
  643.     Int4 D_diff, gd;
  644.     Int4 M_half;
  645.     Int4 max_cost;
  646.     Int4 *lower, *upper;
  647.     
  648.     Int4 x, cur_max, b_diag = 0, best_diag = INT4_MAX/2;
  649.     char nlower = 0, nupper = 0;
  650.     SMBSpace* space = gamp->space;
  651.     Int4 stop_condition;
  652.     Int4 max_d;
  653.     Int4* uplow_free;
  654.     Int4 max_len = len2;
  655.  
  656.     if (match_score % 2 == 1) {
  657.         match_score *= 2;
  658.         mismatch_score *= 2;
  659.         xdrop_threshold *= 2;
  660.         gap_open *= 2;
  661. gap_extend *= 2;
  662.     }
  663.     M_half = match_score/2;
  665.     if (gap_open == 0 && gap_extend == 0) {
  666.        return BLAST_GreedyAlign(s1, len1, s2, len2, reverse, 
  667.                                    xdrop_threshold, match_score, 
  668.                                    mismatch_score, e1, e2, gamp, S, rem);
  669.     }
  670.     
  671.     Mis_cost = mismatch_score + match_score;
  672.     GO_cost = gap_open;
  673.     GE_cost = gap_extend+M_half;
  674.     gd = gdb3(&Mis_cost, &GO_cost, &GE_cost);
  675.     D_diff = ICEIL(xdrop_threshold+M_half, gd);
  676.     
  677.     
  678.     MAX_D = (Int4) (len1/ERROR_FRACTION + 1);
  679.     max_d = MAX_D*GE_cost;
  680.     ORIGIN = MAX_D + 2;
  681.     max_cost = MAX(Mis_cost, GO_cost+GE_cost);
  682.     *e1 = *e2 = 0;
  683.     
  684.     if (reverse) {
  685.        if (!(rem & 4)) {
  686.           for (row = 0; row < len2 && row < len1 && 
  687.                   (s2[len2-1-row] ==
  688.                    NCBI2NA_UNPACK_BASE(s1[(len1-1-row)/4], 
  689.                                         3-(len1-1-row)%4)); 
  690.                row++)
  691.              /*empty*/ ;
  692.        } else {
  693.           for (row = 0; row < len2 && row < len1 && (s2[len2-1-row] == s1[len1-1-row]); row++)
  694.              /*empty*/ ;
  695.        }
  696.     } else {
  697.        if (!(rem & 4)) {
  698.           for (row = 0; row < len2 && row < len1 && 
  699.                   (s2[row] == 
  700.                    NCBI2NA_UNPACK_BASE(s1[(row+rem)/4], 
  701.                                         3-(row+rem)%4)); 
  702.                row++)
  703.              /*empty*/ ;
  704.        } else {
  705.           for (row = 0; row < len2 && row < len1 && (s2[row] == s1[row]); row++)
  706.              /*empty*/ ;
  707.        }
  708.     }
  709.     *e1 = row;
  710.     *e2 = row;
  711.     if (row == len1 || row == len2) {
  712.         if (S != NULL)
  713.             edit_script_rep(S, row);
  714. /* hit last row; stop search */
  715. return row*match_score;
  716.     }
  717.     flast_d = gamp->flast_d_affine;
  718.     if (S==NULL) {
  719.         space = 0;
  720.     } else if (!space) {
  721.        gamp->space = space = MBSpaceNew();
  722.     } else { 
  723.         refresh_mb_space(space);
  724.     }
  725.     max_row = max_row_free + D_diff;
  726.     for (k = 0; k < D_diff; k++)
  727. max_row_free[k] = 0;
  728.     uplow_free = gamp->uplow_free;
  729.     lower = uplow_free;
  730.     upper = uplow_free+max_d+1+max_cost;
  731.     /* next 3 lines set boundary for -1,-2,...,-max_cost+1*/
  732.     for (k = 0; k < max_cost; k++) {lower[k] =LARGE;  upper[k] = -LARGE;}
  733.     lower += max_cost;
  734.     upper += max_cost; 
  735.     
  736.     flast_d[0][ORIGIN].C = row;
  737.     flast_d[0][ORIGIN].I = flast_d[0][ORIGIN].D = -2;
  738.     max_row[0] = (row + row)*M_half;
  739.     lower[0] = upper[0] = ORIGIN;
  740.     
  741.     flower = ORIGIN - 1;
  742.     fupper = ORIGIN + 1;
  743.     
  744.     d = 1;
  745.     stop_condition = 1;
  746.     while (d <= max_d) {
  747. Int4 fl0, fu0;
  748. x = max_row[d - D_diff] + gd * d - xdrop_threshold;
  749. x = ICEIL(x, M_half);
  750. if (x < 0) x=0;
  751. cur_max = 0;
  752. fl0 = flower;
  753. fu0 = fupper;
  754. for (k = fl0; k <= fu0; k++) {
  755.     row = -2;
  756.     if (k+1 <= upper[d-GO_cost-GE_cost] && k+1 >= lower[d-GO_cost-GE_cost]) 
  757.                 row = flast_d[d-GO_cost-GE_cost][k+1].C;
  758.     if (k+1  <= upper[d-GE_cost] && k+1 >= lower[d-GE_cost] &&
  759. row < flast_d[d-GE_cost][k+1].D) 
  760.                 row = flast_d[d-GE_cost][k+1].D;
  761.     row++;
  762.     if (row+row+k-ORIGIN >= x) 
  763.       flast_d[d][k].D = row;
  764.     else flast_d[d][k].D = -2;
  765.     row = -1; 
  766.     if (k-1 <= upper[d-GO_cost-GE_cost] && k-1 >= lower[d-GO_cost-GE_cost]) 
  767.                 row = flast_d[d-GO_cost-GE_cost][k-1].C;
  768.     if (k-1  <= upper[d-GE_cost] && k-1 >= lower[d-GE_cost] &&
  769. row < flast_d[d-GE_cost][k-1].I) 
  770.                 row = flast_d[d-GE_cost][k-1].I;
  771.     if (row+row+k-ORIGIN >= x) 
  772.                 flast_d[d][k].I = row;
  773.     else flast_d[d][k].I = -2;
  774.             
  775.     row = MAX(flast_d[d][k].I, flast_d[d][k].D);
  776.     if (k <= upper[d-Mis_cost] && k >= lower[d-Mis_cost]) 
  777.                 row = MAX(flast_d[d-Mis_cost][k].C+1,row);
  778.             
  779.     col = row + k - ORIGIN;
  780.     if (row + col >= x)
  781. fupper = k;
  782.     else {
  783. if (k == flower)
  784.     flower++;
  785. else
  786.     flast_d[d][k].C = -2;
  787. continue;
  788.     }
  789.             if (row > max_len || row < -2) {
  790.                flower = k; nlower = k+1; 
  791.             } else {
  792.                /* slide down the diagonal */
  793.                if (reverse) {
  794.                   if (!(rem & 4)) {
  795.                      while (row < len2 && col < len1 && s2[len2-1-row] == 
  796.                             NCBI2NA_UNPACK_BASE(s1[(len1-1-col)/4],
  797.                                                  3-(len1-1-col)%4)) {
  798.                         ++row;
  799.                         ++col;
  800.                      }
  801.                   } else {
  802.                      while (row < len2 && col < len1 && s2[len2-1-row] ==
  803.                             s1[len1-1-col]) {
  804.                         ++row;
  805.                         ++col;
  806.                      }
  807.                   }
  808.                } else {
  809.                   if (!(rem & 4)) {
  810.                      while (row < len2 && col < len1 && s2[row] == 
  811.                             NCBI2NA_UNPACK_BASE(s1[(col+rem)/4],
  812.                                                  3-(col+rem)%4)) {
  813.                         ++row;
  814.                         ++col;
  815.                      }
  816.                   } else {
  817.                      while (row < len2 && col < len1 && s2[row] == s1[col]) {
  818.                         ++row;
  819.                         ++col;
  820.                      }
  821.                   }
  822.                }
  823.             }
  824.             flast_d[d][k].C = row;
  825.             if (row + col > cur_max) {
  826.                cur_max = row + col;
  827.                b_diag = k;
  828.             }
  829.             if (row == len2) {
  830.                flower = k; nlower = k+1;
  831.             }
  832.             if (col == len1) {
  833.                fupper = k; nupper = k-1;
  834.             }
  835. }
  836. k = cur_max*M_half - d * gd;
  837. if (max_row[d - 1] < k) {
  838.     max_row[d] = k;
  839.     return_val = d;
  840.     best_diag = b_diag;
  841.     *e2 = flast_d[d][b_diag].C;
  842.     *e1 = (*e2)+b_diag-ORIGIN;
  843. } else {
  844.     max_row[d] = max_row[d - 1];
  845. }
  846. if (flower <= fupper) {
  847.             stop_condition++;
  848.             lower[d] = flower; upper[d] = fupper;
  849. } else {
  850.             lower[d] = LARGE; upper[d] = -LARGE;
  851. }
  852. if (lower[d-max_cost] <= upper[d-max_cost]) stop_condition--;
  853. if (stop_condition == 0) break;
  854. d++;
  855. flower = MIN(lower[d-Mis_cost], MIN(lower[d-GO_cost-GE_cost], lower[d-GE_cost])-1);
  856. if (nlower) flower = MAX(flower, nlower);
  857. fupper = MAX(upper[d-Mis_cost], MAX(upper[d-GO_cost-GE_cost], upper[d-GE_cost])+1);
  858. if (nupper) fupper = MIN(fupper, nupper);
  859. if (d > max_cost) {
  860.    if (S==NULL) {
  861.       /*if (d > max_cost)*/
  862.       flast_d[d] = flast_d[d - max_cost-1];
  863.    } else {
  864.       flast_d[d] = get_mb_space(space, fupper-flower+1)-flower;
  865.       if (flast_d[d] == NULL)
  866.  return return_val;
  867.            }
  868. }
  869.     }
  870.     
  871.     if (S!=NULL) { /*trace back*/
  872.         Int4 row1, diag, state;
  873.         d = return_val; diag = best_diag;
  874.         row = *e2; state = sC;
  875.         while (d > 0) {
  876.             if (state == sC) {
  877.                 /* diag will not be changed*/
  878.                 state = get_lastC(flast_d, lower, upper, &d, diag, Mis_cost, &row1);
  879.                 if (row-row1 > 0) edit_script_rep(S, row-row1);
  880.                 row = row1;
  881.             } else {
  882.                 if (state == sI) {
  883.                     /*row unchanged */
  884.                     state = get_lastI(flast_d, lower, upper, &d, diag, GO_cost, GE_cost);
  885.                     diag--;
  886.                     edit_script_ins(S,1);
  887.                 } else {
  888.                     edit_script_del(S,1);
  889.                     state = get_lastD(flast_d, lower, upper, &d, diag, GO_cost, GE_cost);
  890.                     diag++;
  891.                     row--;
  892.                 }
  893.             }
  894.         }
  895.         edit_script_rep(S, flast_d[0][ORIGIN].C);
  896.         if (!reverse) 
  897.             edit_script_reverse_inplace(S);
  898.     }
  899.     return_val = max_row[return_val];
  900.     return return_val;
  901. }