← Back to branch summary

~ubuntu-branches/ubuntu/karmic/ncbi-tools6/karmic

« back to all changes in this revision

Viewing changes to algo/blast/core/greedy_align.c

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2005-03-27 12:00:15 UTC
  • mfrom: (2.1.2 hoary)
  • Revision ID: james.westby@ubuntu.com-20050327120015-embhesp32nj73p9r
Tags: 6.1.20041020-3
http://bugs.debian.org/295110
* Fix FTBFS under GCC 4.0 caused by inconsistent use of "static" on
  functions.  (Closes: #295110.)
* Add a watch file, now that we can.  (Upstream's layout needs version=3.)

Show diffs side-by-side

added added

removed removed

Lines of Context:
algo/blast/core/greedy_align.c
 
1
/* $Id: greedy_align.c,v 1.24 2004/10/07 14:05:53 papadopo Exp $
 
2
 * ===========================================================================
 
3
 *
 
4
 *                            PUBLIC DOMAIN NOTICE
 
5
 *               National Center for Biotechnology Information
 
6
 *
 
7
 *  This software/database is a "United States Government Work" under the
 
8
 *  terms of the United States Copyright Act.  It was written as part of
 
9
 *  the author's official duties as a United States Government employee and
 
10
 *  thus cannot be copyrighted.  This software/database is freely available
 
11
 *  to the public for use. The National Library of Medicine and the U.S.
 
12
 *  Government have not placed any restriction on its use or reproduction.
 
13
 *
 
14
 *  Although all reasonable efforts have been taken to ensure the accuracy
 
15
 *  and reliability of the software and data, the NLM and the U.S.
 
16
 *  Government do not and cannot warrant the performance or results that
 
17
 *  may be obtained by using this software or data. The NLM and the U.S.
 
18
 *  Government disclaim all warranties, express or implied, including
 
19
 *  warranties of performance, merchantability or fitness for any particular
 
20
 *  purpose.
 
21
 *
 
22
 *  Please cite the author in any work or product based on this material.
 
23
 *
 
24
 * ===========================================================================
 
25
 *
 
26
 * Author: Webb Miller and Co. Adopted for NCBI libraries by Sergey Shavirin
 
27
 *
 
28
 * Initial Creation Date: 10/27/1999
 
29
 *
 
30
 */
 
31
 
 
32
/** @file greedy_align.c
 
33
 * Greedy gapped alignment functions
 
34
 */
 
35
 
 
36
static char const rcsid[] = 
 
37
    "$Id: greedy_align.c,v 1.24 2004/10/07 14:05:53 papadopo Exp $";
 
38
 
 
39
#include <algo/blast/core/greedy_align.h>
 
40
#include <algo/blast/core/blast_util.h> /* for NCBI2NA_UNPACK_BASE macros */
 
41
 
 
42
/** Constants used during greedy alignment */
 
43
enum {
 
44
    MAX_SPACE = 1000000, /**< The number of SThreeVal structures available
 
45
                              in a single link of the greedy memory pool */
 
46
    LARGE = 100000000,   /**< used to represent infinity */
 
47
};
 
48
 
 
49
/** Ensures that an edit script has enough memory allocated
 
50
    to hold a given number of edit operations
 
51
 
 
52
    @param script The script to examine [in/modified]
 
53
    @param total_ops Number of operations the script must support [in]
 
54
    @return 0 if successful, nonzero otherwise
 
55
*/
 
56
static Int2 
 
57
edit_script_realloc(MBGapEditScript *script, Uint4 total_ops)
 
58
{
 
59
    if (script->num_ops_allocated <= total_ops) {
 
60
        Uint4 new_size = total_ops * 3 / 2;
 
61
        MBEditOp *new_ops;
 
62
    
 
63
        new_ops = realloc(script->edit_ops, new_size * 
 
64
                                sizeof(MBEditOp));
 
65
        if (new_ops == NULL)
 
66
            return -1;
 
67
 
 
68
        script->edit_ops = new_ops;
 
69
        script->num_ops_allocated = new_size;
 
70
    }
 
71
    return 0;
 
72
}
 
73
 
 
74
/** Add an edit operation to an edit script
 
75
 
 
76
  @param script The edit script to update [in/modified]
 
77
  @param op_type The edit operation to add [in]
 
78
  @param num_ops The number of operations of the specified type [in]
 
79
  @return 0 on success, nonzero otherwise
 
80
*/
 
81
static Int2 
 
82
edit_script_add_new(MBGapEditScript *script, 
 
83
                enum EOpType op_type, Uint4 num_ops)
 
84
{
 
85
    if (edit_script_realloc(script, script->num_ops + 2) != 0)
 
86
        return -1;
 
87
 
 
88
    ASSERT(op_type != eEditOpError);
 
89
 
 
90
    script->last_op = op_type;
 
91
    script->edit_ops[script->num_ops].op_type = op_type;
 
92
    script->edit_ops[script->num_ops].num_ops = num_ops;
 
93
    script->num_ops++;
 
94
 
 
95
    return 0;
 
96
}
 
97
 
 
98
/** Initialize an edit script structure
 
99
 
 
100
  @param script Pointer to an uninitialized edit script
 
101
  @return The initialized edit script
 
102
*/
 
103
static MBGapEditScript *
 
104
edit_script_init(MBGapEditScript *script)
 
105
{
 
106
        script->edit_ops = NULL;
 
107
        script->num_ops_allocated = 0;
 
108
        script->num_ops = 0;
 
109
        script->last_op = eEditOpError;
 
110
        edit_script_realloc(script, 8);
 
111
        return script;
 
112
}
 
113
 
 
114
/** Add a new operation to an edit script, possibly combining
 
115
    it with the last operation if the two operations are identical
 
116
 
 
117
    @param script The edit script to update [in/modified]
 
118
    @param op_type The operation type to add [in]
 
119
    @param num_ops The number of the specified type of operation to add [in]
 
120
    @return 0 on success, nonzero otherwise
 
121
*/
 
122
static Int2 
 
123
edit_script_add(MBGapEditScript *script, 
 
124
                 enum EOpType op_type, Uint4 num_ops)
 
125
{
 
126
    if (op_type == eEditOpError) {
 
127
#ifdef ERR_POST_EX_DEFINED
 
128
        ErrPostEx(SEV_FATAL, 1, 0, 
 
129
                  "edit_script_more: bad opcode %d:%d", op_type, num_ops);
 
130
#endif
 
131
        return -1;
 
132
    }
 
133
    
 
134
    if (script->last_op == op_type)
 
135
        script->edit_ops[script->num_ops-1].num_ops += num_ops;
 
136
    else
 
137
        edit_script_add_new(script, op_type, num_ops);
 
138
 
 
139
    return 0;
 
140
}
 
141
 
 
142
/** See greedy_align.h for description */
 
143
MBGapEditScript *
 
144
MBGapEditScriptAppend(MBGapEditScript *dest_script, 
 
145
                      MBGapEditScript *src_script)
 
146
{
 
147
    Int4 i;
 
148
 
 
149
    for (i = 0; i < src_script->num_ops; i++) {
 
150
        MBEditOp *edit_op = src_script->edit_ops + i;
 
151
 
 
152
        edit_script_add(dest_script, edit_op->op_type, edit_op->num_ops);
 
153
    }
 
154
 
 
155
    return dest_script;
 
156
}
 
157
 
 
158
/** See greedy_align.h for description */
 
159
MBGapEditScript *
 
160
MBGapEditScriptNew(void)
 
161
{
 
162
    MBGapEditScript *script = calloc(1, sizeof(MBGapEditScript));
 
163
    if (script != NULL)
 
164
        return edit_script_init(script);
 
165
    return NULL;
 
166
}
 
167
 
 
168
/** See greedy_align.h for description */
 
169
MBGapEditScript *
 
170
MBGapEditScriptFree(MBGapEditScript *script)
 
171
{
 
172
    if (script == NULL)
 
173
        return NULL;
 
174
 
 
175
    if (script->edit_ops)
 
176
        sfree(script->edit_ops);
 
177
 
 
178
    sfree(script);
 
179
    return NULL;
 
180
}
 
181
 
 
182
/** Reverse the order of the operations in an edit script
 
183
 
 
184
    @param script The script to be reversed [in/modified]
 
185
    @return Pointer to the updated edit script
 
186
*/
 
187
static MBGapEditScript *
 
188
edit_script_reverse_inplace(MBGapEditScript *script)
 
189
{
 
190
    Uint4 i;
 
191
    const Uint4 num_ops = script->num_ops;
 
192
    const Uint4 mid = num_ops / 2;
 
193
    const Uint4 end = num_ops - 1;
 
194
    
 
195
    for (i = 0; i < mid; ++i) {
 
196
        const MBEditOp t = script->edit_ops[i];
 
197
        script->edit_ops[i] = script->edit_ops[end - i];
 
198
        script->edit_ops[end - i] = t;
 
199
    }
 
200
    return script;
 
201
}
 
202
 
 
203
/** see greedy_align.h for description */
 
204
SMBSpace* 
 
205
MBSpaceNew()
 
206
{
 
207
    SMBSpace* new_space;
 
208
    
 
209
    /** @todo FIXME: Later code assumes that a request will
 
210
       never be made for more than MAX_SPACE structures at once */
 
211
 
 
212
    new_space = (SMBSpace*)malloc(sizeof(SMBSpace));
 
213
    if (new_space == NULL)
 
214
        return NULL;
 
215
 
 
216
    new_space->space_array = (SThreeVal*)malloc(
 
217
                                   MAX_SPACE * sizeof(SThreeVal));
 
218
    if (new_space->space_array == NULL) {
 
219
        sfree(new_space);
 
220
        return NULL;
 
221
    }
 
222
    new_space->space_used = 0; 
 
223
    new_space->space_allocated = MAX_SPACE;
 
224
    new_space->next = NULL;
 
225
 
 
226
    return new_space;
 
227
}
 
228
 
 
229
/** Mark the input collection of space structures as unused
 
230
   @param space The space to mark
 
231
*/
 
232
static void 
 
233
refresh_mb_space(SMBSpace* space)
 
234
{
 
235
    while (space != NULL) {
 
236
        space->space_used = 0;
 
237
        space = space->next;
 
238
    }
 
239
}
 
240
 
 
241
/** see greedy_align.h for description */
 
242
void MBSpaceFree(SMBSpace* space)
 
243
{
 
244
    SMBSpace* next_space;
 
245
 
 
246
    while (space != NULL) {
 
247
        next_space = space->next;
 
248
        sfree(space->space_array);
 
249
        sfree(space);
 
250
        space = next_space;
 
251
    }
 
252
}
 
253
 
 
254
/** Allocate a specified number of SThreeVal structures from
 
255
    a memory pool
 
256
 
 
257
    @param pool The memory pool [in]
 
258
    @param num_alloc The number of structures to allocate [in]
 
259
    @return Pointer to the allocated memory, or NULL in case of error
 
260
*/
 
261
static SThreeVal* 
 
262
get_mb_space(SMBSpace* pool, Int4 num_alloc)
 
263
{
 
264
    SThreeVal* out_ptr;
 
265
    if (num_alloc < 0) 
 
266
        return NULL;  
 
267
    
 
268
    while (pool->space_used + num_alloc > pool->space_allocated) {
 
269
       if (pool->next == NULL) {
 
270
          pool->next = MBSpaceNew();
 
271
          if (pool->next == NULL) {
 
272
#ifdef ERR_POST_EX_DEFINED
 
273
             ErrPostEx(SEV_WARNING, 0, 0, "Cannot get new space for greedy extension");
 
274
#endif
 
275
             return NULL;
 
276
          }
 
277
       }
 
278
       pool = pool->next;
 
279
    }
 
280
 
 
281
    out_ptr = pool->space_array + pool->space_used;
 
282
    pool->space_used += num_alloc;
 
283
    
 
284
    return out_ptr;
 
285
}
 
286
 
 
287
/** Compute the greatest common divisor of two numbers
 
288
    @param a First number [in]
 
289
    @param b Second number [in]
 
290
    @return The computed gcd
 
291
*/
 
292
static Int4 
 
293
gcd(Int4 a, Int4 b)
 
294
{
 
295
    Int4 c;
 
296
    if (a < b) {
 
297
        c = a;
 
298
        a = b; b = c;
 
299
    }
 
300
    while ((c = a % b) != 0) {
 
301
        a = b; b = c;
 
302
    }
 
303
 
 
304
    return b;
 
305
}
 
306
 
 
307
/** Compute the greatest common divisor of three numbers,
 
308
    divide that quantity out of each number, and return the computed GCD
 
309
    @param a Pointer to first number [in/modified]
 
310
    @param b Pointer to second number [in/modified]
 
311
    @param c Pointer to third number [in/modified]
 
312
    @return The computed gcd
 
313
*/
 
314
static Int4 
 
315
gdb3(Int4* a, Int4* b, Int4* c)
 
316
{
 
317
    Int4 g;
 
318
    if (*b == 0) 
 
319
        g = gcd(*a, *c);
 
320
    else 
 
321
        g = gcd(*a, gcd(*b, *c));
 
322
    if (g > 1) {
 
323
        *a /= g;
 
324
        *b /= g;
 
325
        *c /= g;
 
326
    }
 
327
    return g;
 
328
}
 
329
 
 
330
/** During the traceback for a greedy affine alignment,
 
331
    determine the state of the traceback that results
 
332
    from moving a specified number of diagonals upwards
 
333
    in the traceback array from a match
 
334
 
 
335
    @param flast_d 2-D Array of traceback scores [in]
 
336
    @param lower Array of per-diagonal lower bounds [in]
 
337
    @param upper Array of per-diagonal upper bounds [in]
 
338
    @param d Starting score / ending score(?) [in/modified]
 
339
    @param diag Starting diagonal(?) [in]
 
340
    @param Mis_cost Cost of a mismatch [in]
 
341
    @param row1 The row of the traceback array where
 
342
                the next score occurs(?) [out]
 
343
    @return The state for the next traceback operation
 
344
*/
 
345
static enum EOpType 
 
346
get_lastC(SThreeVal** flast_d, Int4* lower, Int4* upper, 
 
347
          Int4* d, Int4 diag, Int4 Mis_cost, Int4* row1)
 
348
{
 
349
    Int4 row;
 
350
    
 
351
    if (diag >= lower[(*d)-Mis_cost] && 
 
352
        diag <= upper[(*d)-Mis_cost]) {
 
353
 
 
354
        row = flast_d[(*d)-Mis_cost][diag].C;
 
355
        if (row >= MAX(flast_d[*d][diag].I, flast_d[*d][diag].D)) {
 
356
            *d = *d-Mis_cost;
 
357
            *row1 = row;
 
358
            return eEditOpReplace;
 
359
        }
 
360
    }
 
361
    if (flast_d[*d][diag].I > flast_d[*d][diag].D) {
 
362
        *row1 = flast_d[*d][diag].I;
 
363
        return eEditOpInsert;
 
364
    } 
 
365
    else {
 
366
        *row1 = flast_d[*d][diag].D;
 
367
        return eEditOpDelete;
 
368
    }
 
369
}
 
370
 
 
371
/** During the traceback for a greedy affine alignment,
 
372
    determine the state of the traceback that results
 
373
    from moving a specified number of diagonals upwards
 
374
    in the traceback array from an insertion or deletion
 
375
 
 
376
    @param flast_d 2-D Array of scores [in]
 
377
    @param lower Array of per-diagonal lower bounds [in]
 
378
    @param upper Array of per-diagonal upper bounds [in]
 
379
    @param d Starting score / ending score(?) [in/modified]
 
380
    @param diag Starting diagonal(?) [in]
 
381
    @param GO_cost Cost to open a gap [in]
 
382
    @param GE_cost Cost to extend a gap [in]
 
383
    @param IorD The state of the traceback at present [in]
 
384
    @return The state for the next traceback operation
 
385
*/
 
386
static enum EOpType 
 
387
get_last_ID(SThreeVal** flast_d, Int4* lower, Int4* upper, 
 
388
            Int4* d, Int4 diag, Int4 GO_cost, 
 
389
            Int4 GE_cost, enum EOpType IorD)
 
390
{
 
391
    Int4 ndiag; 
 
392
    Int4 row;
 
393
 
 
394
    if (IorD == eEditOpInsert)
 
395
        ndiag = diag - 1;
 
396
    else 
 
397
        ndiag = diag + 1;
 
398
 
 
399
    if (ndiag >= lower[(*d)-GE_cost] && 
 
400
        ndiag <= upper[(*d)-GE_cost]) {
 
401
 
 
402
        if (IorD == eEditOpInsert)
 
403
            row = flast_d[(*d)-GE_cost][ndiag].I;
 
404
        else
 
405
            row = flast_d[(*d)-GE_cost][ndiag].D;
 
406
    }
 
407
    else {
 
408
        row = -100;
 
409
    }
 
410
 
 
411
    if (ndiag >= lower[(*d)-GO_cost-GE_cost] && 
 
412
        ndiag <= upper[(*d)-GO_cost-GE_cost] && 
 
413
        row < flast_d[(*d)-GO_cost-GE_cost][ndiag].C) {
 
414
 
 
415
        *d = (*d) - GO_cost - GE_cost;
 
416
        return eEditOpReplace;
 
417
    }
 
418
    *d = (*d) - GE_cost;
 
419
    return IorD;
 
420
}
 
421
 
 
422
/** During the traceback for a non-affine greedy alignment,
 
423
    compute the diagonal that the next traceback operation
 
424
    will use
 
425
 
 
426
    @param flast_d 2-D Array of scores [in]
 
427
    @param d Starting score(?) [in]
 
428
    @param diag Starting diagonal(?) [in]
 
429
    @param row1 The next traceback row to examine(?) [out]
 
430
    @return The next diagonal in the traceback
 
431
*/
 
432
static Int4 
 
433
get_last(Int4 **flast_d, Int4 d, Int4 diag, Int4 *row1)
 
434
{
 
435
    if (flast_d[d-1][diag-1] > MAX(flast_d[d-1][diag], flast_d[d-1][diag+1])) {
 
436
        *row1 = flast_d[d-1][diag-1];
 
437
        return diag - 1;
 
438
    } 
 
439
    if (flast_d[d-1][diag] > flast_d[d-1][diag+1]) {
 
440
        *row1 = flast_d[d-1][diag];
 
441
        return diag;
 
442
    }
 
443
    *row1 = flast_d[d-1][diag+1];
 
444
    return diag + 1;
 
445
}
 
446
 
 
447
/** see greedy_align.h for description */
 
448
Int4 BLAST_GreedyAlign(const Uint1* seq1, Int4 len1,
 
449
                       const Uint1* seq2, Int4 len2,
 
450
                       Boolean reverse, Int4 xdrop_threshold, 
 
451
                       Int4 match_cost, Int4 mismatch_cost,
 
452
                       Int4* extent1, Int4* extent2, 
 
453
                       SGreedyAlignMem* aux_data, 
 
454
                       MBGapEditScript *script, Uint1 rem)
 
455
{
 
456
    Int4 col;               /* column number */
 
457
    Int4 row;               /* row number */
 
458
    Int4 d;                 /* current distance */
 
459
    Int4 k;                 /* current diagonal */
 
460
    Int4 flower, fupper;    /* boundaries for searching diagonals */
 
461
    Int4 MAX_D;             /* maximum cost */
 
462
    Int4 ORIGIN;
 
463
    Int4 final_dist = 0;
 
464
    Int4** flast_d = aux_data->flast_d; /* rows containing the last d */
 
465
    Int4* max_row;                      /* reached for cost d=0, ... len1.  */
 
466
    
 
467
    Int4 match_cost_half = match_cost / 2;
 
468
    Int4 op_cost = mismatch_cost + match_cost;
 
469
    Int4 d_dropoff = ICEIL(xdrop_threshold + match_cost_half, op_cost);
 
470
    
 
471
    Int4 cur_max; 
 
472
    Int4 b_diag = 0; 
 
473
    Int4 best_diag = INT4_MAX / 2;
 
474
    Int4* max_row_free = aux_data->max_row_free;
 
475
    char nlower, nupper;
 
476
    SMBSpace* mem_pool = aux_data->space;
 
477
    Int4 max_len = len2;
 
478
 
 
479
    MAX_D = (Int4) (len1 / GREEDY_MAX_COST_FRACTION + 1);
 
480
    ORIGIN = MAX_D + 2;
 
481
    *extent1 = 0;
 
482
    *extent2 = 0;
 
483
    
 
484
    /* find the offset of the first mismatch between
 
485
       seq1 and seq2 */
 
486
 
 
487
    row = 0;
 
488
    if (reverse) {
 
489
        if (rem == 4) {
 
490
            while (row < len1 && row < len2) {
 
491
                if (seq2[len2-1-row] != seq1[len1-1-row])
 
492
                    break;
 
493
                row++;
 
494
            }
 
495
        } 
 
496
        else {
 
497
            while (row < len1 && row < len2) {
 
498
                if (seq2[len2-1-row] != 
 
499
                                NCBI2NA_UNPACK_BASE(seq1[(len1-1-row) / 4], 
 
500
                                                    3 - (len1-1-row) % 4)) 
 
501
                    break;
 
502
                row++;
 
503
            }
 
504
        }
 
505
    } 
 
506
    else {
 
507
        if (rem == 4) {
 
508
            while (row < len1 && row < len2) {
 
509
                if (seq2[row] != seq1[row])
 
510
                    break; 
 
511
                row++;
 
512
            }
 
513
        } 
 
514
        else {
 
515
            while (row < len1 && row < len2) {
 
516
                if (seq2[row] != NCBI2NA_UNPACK_BASE(seq1[(row+rem) / 4], 
 
517
                                                     3 - (row+rem) % 4))
 
518
                    break;
 
519
                row++;
 
520
            }
 
521
        }
 
522
    }
 
523
 
 
524
    /* update the extents of the alignment, and bail out
 
525
       early if no further work is needed. */
 
526
 
 
527
    *extent1 = row;
 
528
    *extent2 = row;
 
529
    if (row == len1 || row == len2) {
 
530
        if (script != NULL)
 
531
            edit_script_add(script, eEditOpReplace, row);
 
532
        return final_dist;
 
533
    }
 
534
 
 
535
    /* set up the memory pool */
 
536
 
 
537
    if (script == NULL) {
 
538
       mem_pool = NULL;
 
539
    } 
 
540
    else if (mem_pool == NULL) {
 
541
       aux_data->space = mem_pool = MBSpaceNew();
 
542
    } 
 
543
    else { 
 
544
        refresh_mb_space(mem_pool);
 
545
    }
 
546
    
 
547
    max_row = max_row_free + d_dropoff;
 
548
    for (k = 0; k < d_dropoff; k++)
 
549
        max_row_free[k] = 0;
 
550
    
 
551
    flast_d[0][ORIGIN] = row;
 
552
    max_row[0] = row * match_cost;
 
553
    flower = ORIGIN - 1;
 
554
    fupper = ORIGIN + 1;
 
555
    d = 1;
 
556
    nupper = nlower = 0;
 
557
 
 
558
    while (d <= MAX_D) {
 
559
        Int4 x, flower0, fupper0;
 
560
 
 
561
        flast_d[d - 1][flower-1] = -1;
 
562
        flast_d[d - 1][flower] = -1;
 
563
        flast_d[d - 1][fupper] = -1;
 
564
        flast_d[d - 1][fupper+1] = -1;
 
565
 
 
566
        x = max_row[d - d_dropoff] + op_cost * d - xdrop_threshold;
 
567
        x = ICEIL(x, match_cost_half);        
 
568
        cur_max = 0;
 
569
        flower0 = flower;
 
570
        fupper0 = fupper;
 
571
 
 
572
        for (k = flower0; k <= fupper0; k++) {
 
573
            row = MAX(flast_d[d - 1][k + 1], flast_d[d - 1][k]) + 1;
 
574
            row = MAX(row, flast_d[d - 1][k - 1]);
 
575
            col = row + k - ORIGIN;
 
576
            if (row + col >= x) {
 
577
                fupper = k;
 
578
            }
 
579
            else {
 
580
                if (k == flower)
 
581
                    flower++;
 
582
                else
 
583
                    flast_d[d][k] = -1;
 
584
                continue;
 
585
            }
 
586
            
 
587
            if (row > max_len || row < 0) {
 
588
                flower = k + 1; 
 
589
                nlower = 1;
 
590
            } 
 
591
            else {
 
592
                /* Slide down the diagonal. */
 
593
                if (reverse) {
 
594
                    if (rem == 4) {
 
595
                        while (row < len2 && col < len1 && 
 
596
                               seq2[len2-1-row] == seq1[len1-1-col]) {
 
597
                            ++row; ++col;
 
598
                        }
 
599
                    } 
 
600
                    else {
 
601
                        while (row < len2 && col < len1 && seq2[len2-1-row] == 
 
602
                               NCBI2NA_UNPACK_BASE(seq1[(len1-1-col) / 4],
 
603
                                                    3 - (len1-1-col) % 4)) {
 
604
                            ++row; ++col;
 
605
                        }
 
606
                    }
 
607
                } 
 
608
                else {
 
609
                    if (rem == 4) {
 
610
                        while (row < len2 && col < len1 && 
 
611
                               seq2[row] == seq1[col]) {
 
612
                            ++row; ++col;
 
613
                        }
 
614
                    } 
 
615
                    else {
 
616
                        while (row < len2 && col < len1 && seq2[row] == 
 
617
                               NCBI2NA_UNPACK_BASE(seq1[(col+rem) / 4],
 
618
                                                    3 - (col+rem) % 4)) {
 
619
                            ++row; ++col;
 
620
                        }
 
621
                    }
 
622
                }
 
623
            }
 
624
            flast_d[d][k] = row;
 
625
            if (row + col > cur_max) {
 
626
                cur_max = row + col;
 
627
                b_diag = k;
 
628
            }
 
629
            if (row == len2) {
 
630
                flower = k + 1; 
 
631
                nlower = 1;
 
632
            }
 
633
            if (col == len1) {
 
634
                fupper = k - 1; 
 
635
                nupper = 1;
 
636
            }
 
637
        }
 
638
 
 
639
        k = cur_max * match_cost_half - d * op_cost;
 
640
        if (max_row[d - 1] < k) {
 
641
            max_row[d] = k;
 
642
            final_dist = d;
 
643
            best_diag = b_diag;
 
644
            *extent2 = flast_d[d][b_diag];
 
645
            *extent1 = (*extent2) + b_diag - ORIGIN;
 
646
        } 
 
647
        else {
 
648
            max_row[d] = max_row[d - 1];
 
649
        }
 
650
        if (flower > fupper)
 
651
            break;
 
652
 
 
653
        d++;
 
654
        if (nlower == 0) 
 
655
            flower--; 
 
656
        if (nupper == 0) 
 
657
            fupper++;
 
658
        if (script == NULL) {
 
659
           flast_d[d] = flast_d[d - 2];
 
660
        }
 
661
        else {
 
662
           /* space array consists of SThreeVal structures which are 
 
663
              3 times larger than Int4, so divide requested amount by 3
 
664
           */
 
665
           flast_d[d] = (Int4*) get_mb_space(mem_pool, 
 
666
                                             (fupper - flower + 7) / 3);
 
667
           if (flast_d[d] != NULL)
 
668
              flast_d[d] = flast_d[d] - flower + 2;
 
669
           else
 
670
              return final_dist;
 
671
        }
 
672
    }
 
673
    
 
674
    /* perform traceback if desired */
 
675
 
 
676
    if (script != NULL) {
 
677
        Int4 diag;
 
678
 
 
679
        d = final_dist; 
 
680
        diag = best_diag;
 
681
        row = *extent2; 
 
682
        col = *extent1;
 
683
 
 
684
        while (d > 0) {
 
685
            Int4 row1, col1, diag1;
 
686
            diag1 = get_last(flast_d, d, diag, &row1);
 
687
            col1 = row1 + diag1 - ORIGIN;
 
688
            if (diag1 == diag) {
 
689
                if (row - row1 > 0) {
 
690
                    edit_script_add(script, eEditOpReplace, row - row1);
 
691
                }
 
692
            } 
 
693
            else if (diag1 < diag) {
 
694
                if (row - row1 > 0) {
 
695
                    edit_script_add(script, eEditOpReplace, row - row1);
 
696
                }
 
697
                edit_script_add(script, eEditOpInsert, 1);
 
698
            } 
 
699
            else {
 
700
                if (row - row1 - 1 > 0) {
 
701
                    edit_script_add(script, eEditOpReplace, row - row1 - 1);
 
702
                }
 
703
                edit_script_add(script, eEditOpDelete, 1);
 
704
            }
 
705
            d--; 
 
706
            diag = diag1; 
 
707
            col = col1; 
 
708
            row = row1;
 
709
        }
 
710
        edit_script_add(script, eEditOpReplace, flast_d[0][ORIGIN]);
 
711
        if (!reverse) 
 
712
            edit_script_reverse_inplace(script);
 
713
    }
 
714
    return final_dist;
 
715
}
 
716
 
 
717
/** See greedy_align.h for description */
 
718
Int4 BLAST_AffineGreedyAlign (const Uint1* seq1, Int4 len1, 
 
719
                              const Uint1* seq2, Int4 len2,
 
720
                              Boolean reverse, Int4 xdrop_threshold, 
 
721
                              Int4 match_score, Int4 mismatch_score, 
 
722
                              Int4 gap_open, Int4 gap_extend,
 
723
                              Int4* extent1, Int4* extent2, 
 
724
                              SGreedyAlignMem* aux_data, 
 
725
                              MBGapEditScript *script, Uint1 rem)
 
726
{
 
727
    Int4 col;                        /* column number */
 
728
    Int4 row;                        /* row number */
 
729
    Int4 d;                        /* current distance */
 
730
    Int4 k;                        /* current diagonal */
 
731
    Int4 flower, fupper;         /* boundaries for searching diagonals */
 
732
    Int4 MAX_D;                         /* maximum cost */
 
733
    Int4 ORIGIN;
 
734
    Int4 final_score = 0;
 
735
    SThreeVal** flast_d;        /* rows containing the last d */
 
736
    Int4* max_row_free;
 
737
    Int4* max_row;                /* reached for cost d=0, ... len1.  */
 
738
    Int4 Mis_cost, GO_cost, GE_cost;
 
739
    Int4 D_diff, gd;
 
740
    Int4 match_score_half;
 
741
    Int4 max_cost;
 
742
    Int4 *lower, *upper;
 
743
    
 
744
    Int4 cur_max; 
 
745
    Int4 b_diag = 0; 
 
746
    Int4 best_diag = INT4_MAX/2;
 
747
    char nlower = 0, nupper = 0;
 
748
    SMBSpace* mem_pool = aux_data->space;
 
749
    Int4 stop_condition;
 
750
    Int4 max_d;
 
751
    Int4* uplow_free;
 
752
    Int4 max_len = len2;
 
753
 
 
754
    if (match_score % 2 == 1) {
 
755
        match_score *= 2;
 
756
        mismatch_score *= 2;
 
757
        xdrop_threshold *= 2;
 
758
        gap_open *= 2;
 
759
        gap_extend *= 2;
 
760
    }
 
761
 
 
762
    if (gap_open == 0 && gap_extend == 0) {
 
763
       return BLAST_GreedyAlign(seq1, len1, seq2, len2, reverse, 
 
764
                                   xdrop_threshold, match_score, 
 
765
                                   mismatch_score, extent1, extent2, 
 
766
                                   aux_data, script, rem);
 
767
    }
 
768
    
 
769
    match_score_half = match_score / 2;
 
770
    Mis_cost = mismatch_score + match_score;
 
771
    GO_cost = gap_open;
 
772
    GE_cost = gap_extend + match_score_half;
 
773
    gd = gdb3(&Mis_cost, &GO_cost, &GE_cost);
 
774
    D_diff = ICEIL(xdrop_threshold + match_score_half, gd);
 
775
    
 
776
    MAX_D = (Int4) (len1/GREEDY_MAX_COST_FRACTION + 1);
 
777
    max_d = MAX_D * GE_cost;
 
778
    ORIGIN = MAX_D + 2;
 
779
    max_cost = MAX(Mis_cost, GO_cost + GE_cost);
 
780
    *extent1 = 0;
 
781
    *extent2 = 0;
 
782
    
 
783
    /* find the offset of the first mismatch between
 
784
       seq1 and seq2 */
 
785
 
 
786
    row = 0;
 
787
    if (reverse) {
 
788
        if (rem == 4) {
 
789
            while (row < len1 && row < len2) {
 
790
                if (seq2[len2-1-row] != seq1[len1-1-row])
 
791
                    break;
 
792
                row++;
 
793
            }
 
794
        } 
 
795
        else {
 
796
            while (row < len1 && row < len2) {
 
797
                if (seq2[len2-1-row] != 
 
798
                                NCBI2NA_UNPACK_BASE(seq1[(len1-1-row) / 4], 
 
799
                                                    3 - (len1-1-row) % 4)) 
 
800
                    break;
 
801
                row++;
 
802
            }
 
803
        }
 
804
    } 
 
805
    else {
 
806
        if (rem == 4) {
 
807
            while (row < len1 && row < len2) {
 
808
                if (seq2[row] != seq1[row])
 
809
                    break; 
 
810
                row++;
 
811
            }
 
812
        } 
 
813
        else {
 
814
            while (row < len1 && row < len2) {
 
815
                if (seq2[row] != NCBI2NA_UNPACK_BASE(seq1[(row+rem) / 4], 
 
816
                                                     3 - (row+rem) % 4))
 
817
                    break;
 
818
                row++;
 
819
            }
 
820
        }
 
821
    }
 
822
 
 
823
    /* update the extents of the alignment, and bail out
 
824
       early if no further work is needed */
 
825
 
 
826
    *extent1 = row;
 
827
    *extent2 = row;
 
828
    if (row == len1 || row == len2) {
 
829
        if (script != NULL)
 
830
            edit_script_add(script, eEditOpReplace, row);
 
831
        return row * match_score;
 
832
    }
 
833
 
 
834
    /* set up the memory pool */
 
835
 
 
836
    if (script == NULL) {
 
837
        mem_pool = NULL;
 
838
    } 
 
839
    else if (!mem_pool) {
 
840
       aux_data->space = mem_pool = MBSpaceNew();
 
841
    } 
 
842
    else { 
 
843
        refresh_mb_space(mem_pool);
 
844
    }
 
845
 
 
846
    flast_d = aux_data->flast_d_affine;
 
847
    max_row_free = aux_data->max_row_free;
 
848
    max_row = max_row_free + D_diff;
 
849
    for (k = 0; k < D_diff; k++)
 
850
        max_row_free[k] = 0;
 
851
 
 
852
    uplow_free = aux_data->uplow_free;
 
853
    lower = uplow_free;
 
854
    upper = uplow_free + max_d + 1 + max_cost;
 
855
 
 
856
    /* set boundary for -1,-2,...,-max_cost+1*/
 
857
    for (k = 0; k < max_cost; k++) {
 
858
        lower[k] = LARGE;  
 
859
        upper[k] = -LARGE;
 
860
    }
 
861
    lower += max_cost;
 
862
    upper += max_cost; 
 
863
    
 
864
    flast_d[0][ORIGIN].C = row;
 
865
    flast_d[0][ORIGIN].I = -2;
 
866
    flast_d[0][ORIGIN].D = -2;
 
867
    max_row[0] = row * match_score;
 
868
    lower[0] = ORIGIN;
 
869
    upper[0] = ORIGIN;
 
870
    flower = ORIGIN - 1;
 
871
    fupper = ORIGIN + 1;
 
872
    
 
873
    d = 1;
 
874
    stop_condition = 1;
 
875
    while (d <= max_d) {
 
876
        Int4 x, flower0, fupper0;
 
877
 
 
878
        x = max_row[d - D_diff] + gd * d - xdrop_threshold;
 
879
        x = ICEIL(x, match_score_half);
 
880
        if (x < 0) 
 
881
            x = 0;
 
882
        cur_max = 0;
 
883
        flower0 = flower;
 
884
        fupper0 = fupper;
 
885
 
 
886
        for (k = flower0; k <= fupper0; k++) {
 
887
            row = -2;
 
888
            if (k+1 <= upper[d-GO_cost-GE_cost] && 
 
889
                k+1 >= lower[d-GO_cost-GE_cost]) {
 
890
                row = flast_d[d-GO_cost-GE_cost][k+1].C;
 
891
            }
 
892
            if (k+1  <= upper[d-GE_cost] && k+1 >= lower[d-GE_cost] &&
 
893
                row < flast_d[d-GE_cost][k+1].D) {
 
894
                row = flast_d[d-GE_cost][k+1].D;
 
895
            }
 
896
            row++;
 
897
 
 
898
            if (2 * row + k - ORIGIN >= x) {
 
899
                flast_d[d][k].D = row;
 
900
            }
 
901
            else {
 
902
                flast_d[d][k].D = -2;
 
903
            }
 
904
            row = -1; 
 
905
 
 
906
            if (k-1 <= upper[d-GO_cost-GE_cost] && 
 
907
                k-1 >= lower[d-GO_cost-GE_cost]) {
 
908
                row = flast_d[d-GO_cost-GE_cost][k-1].C;
 
909
            }
 
910
            if (k-1 <= upper[d-GE_cost] && 
 
911
                k-1 >= lower[d-GE_cost] &&
 
912
                row < flast_d[d-GE_cost][k-1].I) {
 
913
                row = flast_d[d-GE_cost][k-1].I;
 
914
            }
 
915
            if (2 * row + k - ORIGIN >= x) {
 
916
                flast_d[d][k].I = row;
 
917
            }
 
918
            else {
 
919
                flast_d[d][k].I = -2;
 
920
            }
 
921
            
 
922
            row = MAX(flast_d[d][k].I, flast_d[d][k].D);
 
923
            if (k <= upper[d-Mis_cost] && 
 
924
                k >= lower[d-Mis_cost]) {
 
925
                row = MAX(flast_d[d-Mis_cost][k].C+1,row);
 
926
            }
 
927
            
 
928
            col = row + k - ORIGIN;
 
929
            if (row + col >= x) {
 
930
                fupper = k;
 
931
            }
 
932
            else {
 
933
                if (k == flower)
 
934
                    flower++;
 
935
                else
 
936
                    flast_d[d][k].C = -2;
 
937
                continue;
 
938
            }
 
939
 
 
940
            if (row > max_len || row < -2) {
 
941
                flower = k; nlower = k+1; 
 
942
            } 
 
943
            else {
 
944
                /* slide down the diagonal */
 
945
                if (reverse) {
 
946
                    if (rem == 4) {
 
947
                        while (row < len2 && col < len1 && 
 
948
                               seq2[len2-1-row] == seq1[len1-1-col]) {
 
949
                            ++row; ++col;
 
950
                        }
 
951
                    } 
 
952
                    else {
 
953
                        while (row < len2 && col < len1 && seq2[len2-1-row] == 
 
954
                               NCBI2NA_UNPACK_BASE(seq1[(len1-1-col) / 4],
 
955
                                                    3 - (len1-1-col) % 4)) {
 
956
                            ++row; ++col;
 
957
                        }
 
958
                    }
 
959
                } 
 
960
                else {
 
961
                    if (rem == 4) {
 
962
                        while (row < len2 && col < len1 && 
 
963
                               seq2[row] == seq1[col]) {
 
964
                            ++row; ++col;
 
965
                        }
 
966
                    } 
 
967
                    else {
 
968
                        while (row < len2 && col < len1 && seq2[row] == 
 
969
                               NCBI2NA_UNPACK_BASE(seq1[(col+rem) / 4],
 
970
                                                    3 - (col+rem) % 4)) {
 
971
                            ++row; ++col;
 
972
                        }
 
973
                    }
 
974
                }
 
975
            }
 
976
            flast_d[d][k].C = row;
 
977
            if (row + col > cur_max) {
 
978
                cur_max = row + col;
 
979
                b_diag = k;
 
980
            }
 
981
            if (row == len2) {
 
982
                flower = k; 
 
983
                nlower = k+1;
 
984
            }
 
985
            if (col == len1) {
 
986
                fupper = k; 
 
987
                nupper = k-1;
 
988
            }
 
989
        }
 
990
 
 
991
        k = cur_max * match_score_half - d * gd;
 
992
        if (max_row[d - 1] < k) {
 
993
            max_row[d] = k;
 
994
            final_score = d;
 
995
            best_diag = b_diag;
 
996
            *extent2 = flast_d[d][b_diag].C;
 
997
            *extent1 = (*extent2) + b_diag - ORIGIN;
 
998
        } 
 
999
        else {
 
1000
            max_row[d] = max_row[d - 1];
 
1001
        }
 
1002
 
 
1003
        if (flower <= fupper) {
 
1004
            stop_condition++;
 
1005
            lower[d] = flower; 
 
1006
            upper[d] = fupper;
 
1007
        } 
 
1008
        else {
 
1009
            lower[d] = LARGE; 
 
1010
            upper[d] = -LARGE;
 
1011
        }
 
1012
 
 
1013
        if (lower[d-max_cost] <= upper[d-max_cost]) 
 
1014
            stop_condition--;
 
1015
        if (stop_condition == 0) 
 
1016
            break;
 
1017
        
 
1018
        d++;
 
1019
        flower = MIN(lower[d-Mis_cost], 
 
1020
                     MIN(lower[d-GO_cost-GE_cost], lower[d-GE_cost])-1);
 
1021
        if (nlower > 0) 
 
1022
            flower = MAX(flower, nlower);
 
1023
 
 
1024
        fupper = MAX(upper[d-Mis_cost], 
 
1025
                     MAX(upper[d-GO_cost-GE_cost], upper[d-GE_cost])+1);
 
1026
        if (nupper > 0) 
 
1027
            fupper = MIN(fupper, nupper);
 
1028
 
 
1029
        if (d > max_cost) {
 
1030
           if (script == NULL) {
 
1031
               flast_d[d] = flast_d[d - max_cost-1];
 
1032
           } 
 
1033
           else {
 
1034
               flast_d[d] = get_mb_space(mem_pool, fupper-flower+1)-flower;
 
1035
               if (flast_d[d] == NULL)
 
1036
                   return final_score;
 
1037
           }
 
1038
        }
 
1039
    }
 
1040
    
 
1041
    /* compute the traceback if desired */
 
1042
 
 
1043
    if (script != NULL) { 
 
1044
        Int4 row1, diag; 
 
1045
        enum EOpType state;
 
1046
 
 
1047
        d = final_score; 
 
1048
        diag = best_diag;
 
1049
        row = *extent2; 
 
1050
        state = eEditOpReplace;
 
1051
        while (d > 0) {
 
1052
            if (state == eEditOpReplace) {
 
1053
                /* diag unchanged */
 
1054
                state = get_lastC(flast_d, lower, upper, 
 
1055
                                  &d, diag, Mis_cost, &row1);
 
1056
                if (row - row1 > 0) 
 
1057
                    edit_script_add(script, eEditOpReplace, row-row1);
 
1058
                row = row1;
 
1059
            } 
 
1060
            else {
 
1061
                if (state == eEditOpInsert) {
 
1062
                    /*row unchanged */
 
1063
                    state = get_last_ID(flast_d, lower, upper, &d, 
 
1064
                                  diag, GO_cost, GE_cost, eEditOpInsert);
 
1065
                    diag--;
 
1066
                    edit_script_add(script, eEditOpInsert, 1);
 
1067
                } 
 
1068
                else {
 
1069
                    edit_script_add(script, eEditOpDelete, 1);
 
1070
                    state = get_last_ID(flast_d, lower, upper, &d, 
 
1071
                                  diag, GO_cost, GE_cost, eEditOpDelete);
 
1072
                    diag++;
 
1073
                    row--;
 
1074
                }
 
1075
            }
 
1076
        }
 
1077
        edit_script_add(script, eEditOpReplace, flast_d[0][ORIGIN].C);
 
1078
        if (!reverse) 
 
1079
            edit_script_reverse_inplace(script);
 
1080
    }
 
1081
    final_score = max_row[final_score];
 
1082
    return final_score;
 
1083
}