1 // ==========================================================================
2 // SeqAn - The Library for Sequence Analysis
3 // ==========================================================================
4 // Copyright (c) 2006-2018, Knut Reinert, FU Berlin
5 // All rights reserved.
6 //
7 // Redistribution and use in source and binary forms, with or without
8 // modification, are permitted provided that the following conditions are met:
9 //
10 // * Redistributions of source code must retain the above copyright
11 // notice, this list of conditions and the following disclaimer.
12 // * Redistributions in binary form must reproduce the above copyright
13 // notice, this list of conditions and the following disclaimer in the
14 // documentation and/or other materials provided with the distribution.
15 // * Neither the name of Knut Reinert or the FU Berlin nor the names of
16 // its contributors may be used to endorse or promote products derived
17 // from this software without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 // ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE
23 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
25 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
26 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27 // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
29 // DAMAGE.
30 //
31 // ==========================================================================
32 // Author: Stephan Aiche <stephan.aiche@fu-berlin.de>
33 // ==========================================================================
34
35 #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_GLOBAL_ALIGNMENT_MYERS_HIRSCHBERG_IMPL_H_
36 #define SEQAN_INCLUDE_SEQAN_ALIGN_GLOBAL_ALIGNMENT_MYERS_HIRSCHBERG_IMPL_H_
37
38 #include <bitset>
39
40 namespace seqan {
41
42 // ============================================================================
43 // Forwards
44 // ============================================================================
45
46 // ============================================================================
47 // Tags, Classes, Enums
48 // ============================================================================
49
50 // ============================================================================
51 // Metafunctions
52 // ============================================================================
53
54 // ============================================================================
55 // Functions
56 // ============================================================================
57
58 // ----------------------------------------------------------------------------
59 // Function _writeDebugMatrix()
60 // ----------------------------------------------------------------------------
61
62 #ifdef MYERS_HIRSCHBERG_VERBOSE
63 template<typename TSource>
_writeDebugMatrix(TSource s1,TSource s2)64 void _writeDebugMatrix(TSource s1,TSource s2)
65 {
66 //IOREV _notio_ not relevant for iorev
67 int l1 = length(s1);
68 int l2 = length(s2);
69
70 int i,j,sg,sd;
71
72 String<String<int> > fMatrix,rMatrix,tMatrix;
73
74 resize(fMatrix,l1 + 1);
75 resize(rMatrix,l1 + 1);
76 resize(tMatrix,l1 + 1);
77
78 for(i = 0;i <= l1;++i)
79 {
80 resize(fMatrix[i],l2 + 1);
81 resize(rMatrix[i],l2 + 1);
82 resize(tMatrix[i],l2 + 1);
83 }
84
85 for(i = 0;i <= l1;++i)
86 fMatrix[i][0] = i * (-1);
87
88 for(i = l1;i >= 0;--i)
89 rMatrix[i][l2] = (l1 - i) * (-1);
90
91 // calculate forward matrix
92 for(j = 1;j <= l2;++j)
93 {
94 fMatrix[0][j] = j*(-1);
95 for(i = 1;i <= l1;++i)
96 {
97 sg = -1 + ((fMatrix[i-1][j] > fMatrix[i][j-1]) ? fMatrix[i-1][j] : fMatrix[i][j-1]);
98 sd = fMatrix[i-1][j-1] + ((s1[i - 1] == s2[j-1]) ? 0 : -1 );
99
100 fMatrix[i][j] = ((sg > sd) ? sg : sd);
101 }
102 }
103
104 // calculate reverse matrix
105 for(j = l2 - 1;j >= 0;--j)
106 {
107 rMatrix[l1][j] = (l2 - j)*(-1);
108 for(i = l1 - 1;i >= 0;--i)
109 {
110 sg = -1 + ((rMatrix[i+1][j] > rMatrix[i][j+1]) ? rMatrix[i+1][j] : rMatrix[i][j+1]);
111 sd = rMatrix[i+1][j+1] + ((s1[i] == s2[j]) ? 0 : -1 );
112
113 rMatrix[i][j] = ((sg > sd) ? sg : sd);
114 }
115 }
116
117 // print fMatrix
118 std::cout << ";-;";
119 for(i = 0;i < l1;++i)
120 std::cout << s1[i] << ";";
121
122 std::cout << std::endl << "-;";
123 for(j = 0;j <= l2;++j)
124 {
125 if(j != 0) std::cout << s2[j-1] << ";";
126 for(i = 0;i <= l1;++i)
127 {
128 std::cout << fMatrix[i][j] << ";";
129 }
130 std::cout << std::endl;
131 }
132 // print rMatrix
133 std::cout << ";";
134 for(i = 0;i < l1;++i)
135 std::cout << s1[i] << ";";
136 std::cout << "-;" << std::endl;
137
138 for(j = 0;j <= l2;++j)
139 {
140 if(j != l2) std::cout << s2[j] << ";";
141 else std::cout << "-;";
142 for(i = 0;i <= l1;++i)
143 {
144 std::cout << rMatrix[i][j] << ";";
145 }
146 std::cout << std::endl;
147 }
148
149 // fill and print target matrix
150 std::cout << ";-;";
151 for(i = 0;i < l1;++i)
152 std::cout << s1[i] << ";";
153
154 std::cout << std::endl << "-;";
155 for(j = 0;j <= l2;++j)
156 {
157 if(j != 0) std::cout << s2[j-1] << ";";
158 for(i = 0;i <= l1;++i)
159 {
160 tMatrix[i][j] = fMatrix[i][j] + rMatrix[i][j];
161 std::cout << tMatrix[i][j] << ";";
162 }
163 std::cout << std::endl;
164 }
165 }
166 #endif
167
168 // ----------------------------------------------------------------------------
169 // Function _printBinary()
170 // ----------------------------------------------------------------------------
171
172 // Debug integer types in binary format.
173 template <typename TStream, typename T>
174 inline void
_printBinary(TStream & stream,T const n)175 _printBinary(TStream & stream, T const n)
176 {
177 std::bitset<BitsPerValue<T>::VALUE> bits(n);
178 stream << bits;
179 stream << '\n';
180 }
181
182 // ----------------------------------------------------------------------------
183 // Function globalAlignment()
184 // ----------------------------------------------------------------------------
185
186 // When using different alphabets, we will internally use the pattern alphabet
187 // for the comparison. This means that the text character is converted to the
188 // pattern alphabet. Here, the "pattern" is the shorter source sequence.
189
190 template <typename TSequenceH, typename TGapsSpecH, typename TSequenceV, typename TGapsSpecV>
191 int
_globalAlignment(Gaps<TSequenceH,TGapsSpecH> & gapsH,Gaps<TSequenceV,TGapsSpecV> & gapsV,MyersHirschberg const & algorithmTag)192 _globalAlignment(Gaps<TSequenceH, TGapsSpecH> & gapsH,
193 Gaps<TSequenceV, TGapsSpecV> & gapsV,
194 MyersHirschberg const & algorithmTag)
195 {
196 // Switch horizontal and vertical gap roles, gapsV should be the shorter one
197 // to fit into less words.
198 if (length(source(gapsH)) < length(source(gapsV)))
199 return _globalAlignment(gapsV, gapsH, algorithmTag);
200
201 clearGaps(gapsH);
202 clearGaps(gapsV);
203 clearClipping(gapsH);
204 clearClipping(gapsV);
205
206 typedef int TScoreValue;
207
208 // use size of unsigned int as blocksize for bit-vectors
209 const unsigned int BLOCK_SIZE = BitsPerValue<unsigned int>::VALUE;
210
211 // saves the score value that will be returned
212 TScoreValue score,total_score = 0;
213
214 typedef typename Value<TSequenceV>::Type TPatternAlphabet;
215 typedef typename Size<TSequenceH>::Type TStringSize;
216
217 typedef typename Iterator<TSequenceH const, Rooted>::Type TSequenceHIterator;
218 typedef typename Iterator<TSequenceV const, Rooted>::Type TSequenceVIterator;
219 typedef Gaps<TSequenceH, TGapsSpecH> TGapsH;
220 typedef Gaps<TSequenceV, TGapsSpecV> TGapsV;
221 typedef typename Iterator<TGapsH, Rooted>::Type TGapsHIterator;
222 typedef typename Iterator<TGapsV, Rooted>::Type TGapsVIterator;
223
224 typedef typename Iterator<Matrix<TScoreValue>, Rooted>::Type TMatrixIterator;
225
226 TGapsHIterator target_0 = begin(gapsH);
227 TGapsVIterator target_1 = begin(gapsV);
228
229 TSequenceH const & x = source(gapsH);
230 TSequenceV const & y = source(gapsV);
231
232 TStringSize len_x = length(x);
233 TStringSize len_y = length(y);
234
235 // string to store the score values for the currently active cell
236 String<TScoreValue> c_score;
237 resize(c_score, len_x + 1, 0);
238
239 // scoring-scheme specific score values
240 TScoreValue score_match = 0;
241 TScoreValue score_mismatch = -1;
242 TScoreValue score_gap = -1;
243
244 // additional vars
245 unsigned i;
246
247 // stack with parts of matrix that have to be processed
248 std::stack<HirschbergSet_> to_process;
249 HirschbergSet_ target;
250
251 // myers specific vars and preprocessing
252 unsigned int patternAlphabetSize = ValueSize<TPatternAlphabet>::VALUE;
253 unsigned int blockCount = (len_y + BLOCK_SIZE - 1) / BLOCK_SIZE; // maximal count of blocks
254
255 String<unsigned> VP;
256 String<unsigned> VN;
257 String<unsigned> forwardBitMask;
258 String<unsigned> reverseBitMask;
259
260 resize(VP, blockCount, std::numeric_limits<unsigned>::max());
261 resize(VN, blockCount, 0);
262
263 // first bitMask will be constructed from the shorter sequence
264 resize(forwardBitMask, patternAlphabetSize * blockCount, 0);
265 resize(reverseBitMask, patternAlphabetSize * blockCount, 0);
266
267 // encoding the letters as bit-vectors
268 for (unsigned int j = 0; j < len_y; j++){
269 forwardBitMask[blockCount * ordValue(getValue(y,j)) + j/BLOCK_SIZE] = forwardBitMask[blockCount * ordValue(getValue(y,j)) + j/BLOCK_SIZE] | 1 << (j%BLOCK_SIZE);
270 reverseBitMask[blockCount * ordValue(getValue(y,len_y - j - 1)) + j/BLOCK_SIZE] = reverseBitMask[blockCount * ordValue(getValue(y,len_y - j - 1)) + j/BLOCK_SIZE] | 1 << (j%BLOCK_SIZE);
271 }
272
273 HirschbergSet_ hs_complete{0, static_cast<unsigned>(len_x), 0, static_cast<unsigned>(len_y), 1};
274 to_process.push(hs_complete);
275
276 while(!to_process.empty())
277 {
278 target = to_process.top();
279 to_process.pop();
280 /* if score is zero, the whole part of the sequence can be simply skipped */
281 if(_score(target) == 0)
282 {
283 /* coukd work faster */
284 for(i = 0;i < (_end1(target) - _begin1(target));++i)
285 {
286 ++target_0;
287 ++target_1;
288 }
289
290 #ifdef MYERS_HIRSCHBERG_VERBOSE
291 printf("skipped %i to %i in first sequence\n",_begin1(target),_end1(target));
292 #endif
293 }
294 else if(_begin1(target) == _end1(target))
295 {
296
297 #ifdef MYERS_HIRSCHBERG_VERBOSE
298 std::cout << "align y " << _begin2(target) << " to " << _end2(target) << std::endl;
299 std::cout << "align " << infix(y,_begin2(target),_end2(target)) << std::endl << std::endl;
300 #endif
301 for(i = 0;i < (_end2(target) - _begin2(target));++i)
302 {
303 insertGap(target_0);
304 ++target_0;
305 ++target_1;
306 }
307 }
308 else if(_begin2(target) + 1 == _end2(target))
309 {
310 /* ALIGN */
311 #ifdef MYERS_HIRSCHBERG_VERBOSE
312 std::cout << "align x " << _begin1(target) << " to " << _end1(target) << " and y " << _begin2(target) << " to " << _end2(target) << std::endl;
313 std::cout << "align " << infix(x,_begin1(target),_end1(target)) << " and " << infix(y,_begin2(target),_end2(target)) << std::endl << std::endl;
314 #endif
315
316 TStringSize len_1 = _end1(target) - _begin1(target);
317 TStringSize len_2 = _end2(target) - _begin2(target);
318
319 Matrix<TScoreValue> matrix_;
320
321 setDimension(matrix_, 2);
322 setLength(matrix_, 0, len_1 + 1);
323 setLength(matrix_, 1, len_2 + 1);
324 resize(matrix_);
325
326 /* init matrix */
327 TSequenceHIterator xs_begin = iter(x,_begin1(target)) - 1;
328 TSequenceHIterator xs_end = iter(x,_end1(target)) - 1;
329 TSequenceVIterator ys_begin = iter(y,_begin2(target)) - 1;
330 TSequenceVIterator ys_end = iter(y,_end2(target)) - 1;
331
332 TSequenceHIterator xs = xs_end;
333 TSequenceVIterator ys;
334
335 TMatrixIterator col_ = end(matrix_) - 1;
336 TMatrixIterator finger1;
337 TMatrixIterator finger2;
338
339
340 TScoreValue h = 0;
341 TScoreValue border_ = score_gap;
342 TScoreValue v = border_;
343
344
345 //-------------------------------------------------------------------------
346 // init
347
348 finger1 = col_;
349 *finger1 = 0;
350 for (xs = xs_end; xs != xs_begin; --xs)
351 {
352 goPrevious(finger1, 0);
353 *finger1 = border_;
354 border_ += score_gap;
355 }
356
357 //-------------------------------------------------------------------------
358 //fill matrix
359
360 border_ = 0;
361 for (ys = ys_end; ys != ys_begin; --ys)
362 {
363 TPatternAlphabet cy = *ys;
364 h = border_;
365 border_ += score_gap;
366 v = border_;
367
368 finger2 = col_;
369 goPrevious(col_, 1);
370 finger1 = col_;
371
372 *finger1 = v;
373
374 for (xs = xs_end; xs != xs_begin; --xs)
375 {
376 goPrevious(finger1, 0);
377 goPrevious(finger2, 0);
378 if (*xs == cy)
379 {
380 v = h + score_match;
381 h = *finger2;
382 }
383 else
384 {
385 TScoreValue s1 = h + score_mismatch;
386 h = *finger2;
387 TScoreValue s2 = score_gap + ((h > v) ? h : v);
388 v = (s1 > s2) ? s1 : s2;
389 }
390 *finger1 = v;
391 }
392 }
393
394 // if computed the whole matrix last value of v = alignment score
395 if(target == hs_complete) total_score = v;
396
397 /* TRACE BACK */
398 finger1 = begin(matrix_);
399 xs = iter(x,_begin1(target));
400 ys = iter(y,_begin2(target));
401 xs_end = iter(x,_end1(target));
402 ys_end = iter(y,_end2(target));
403
404 while ((xs != xs_end) && (ys != ys_end))
405 {
406 bool gv;
407 bool gh;
408
409 if (*xs == *ys)
410 {
411 gv = gh = true;
412 }
413 else
414 {
415 TMatrixIterator it_ = finger1;
416
417 goNext(it_, 0);
418 TScoreValue v = *it_;
419
420 goNext(it_, 1);
421 TScoreValue d = *it_;
422
423 it_ = finger1;
424 goNext(it_, 1);
425 TScoreValue h = *it_;
426
427 gv = (v >= h) | (d >= h);
428 gh = (h >= v) | (d >= v);
429 }
430
431 if (gv)
432 {
433 ++xs;
434 goNext(finger1, 0);
435 }
436 else
437 {
438 insertGap(target_0);
439 }
440
441 if (gh)
442 {
443 ++ys;
444 goNext(finger1, 1);
445 }
446 else
447 {
448 insertGap(target_1);
449 }
450
451 ++target_0;
452 ++target_1;
453 }
454
455 // if x or y did not reached there end position, fill the rest with gaps
456 while(xs != xs_end)
457 {
458 insertGap(target_1);
459 ++target_0;
460 ++target_1;
461 ++xs;
462 }
463
464 while(ys != ys_end)
465 {
466 insertGap(target_0);
467 ++target_0;
468 ++target_1;
469 ++ys;
470 }
471 /* END ALIGN */
472
473
474 #ifdef MYERS_HIRSCHBERG_VERBOSE
475 std::cout << std::endl << align_ << std::endl << std::endl;
476 #endif
477
478 }
479 else
480 {
481 /*
482 ---------------------------------------------------------------
483 Calculate cut position using extended Myers-Bitvector-Algorithm
484 ---------------------------------------------------------------
485 */
486
487 /* declare variables */
488 unsigned int X, D0, HN, HP;
489
490 /* compute cut position */
491 unsigned mid = static_cast<int>(floor( static_cast<double>((_begin2(target) + _end2(target))/2) ));
492
493 /* debug infos */
494 #ifdef MYERS_HIRSCHBERG_VERBOSE
495 std::cout << "calculate cut for x " << _begin1(target) << " to " << _end1(target) << " and y " << _begin2(target) << " to " << _end2(target) << std::endl;
496 std::cout << "calculate cut for " << infix(x,_begin1(target),_end1(target)) << " and " << infix(y,_begin2(target),_end2(target)) << std::endl;
497 std::cout << "cut is in row " << mid << " symbol is " << getValue(x,mid-1) << std::endl << std::endl;
498
499 std::cout << std::endl;
500 _writeDebugMatrix(infix(x,_begin1(target),_end1(target)),infix(y,_begin2(target),_end2(target)));
501 std::cout << std::endl;
502 #endif
503 /* compute blocks and score masks */
504 unsigned fStartBlock = _begin2(target) / BLOCK_SIZE;
505 unsigned fEndBlock = (mid - 1) / BLOCK_SIZE;
506 unsigned fSpannedBlocks = (fEndBlock - fStartBlock) + 1;
507
508 unsigned int fScoreMask = 1 << ((mid - 1) % BLOCK_SIZE);
509
510 unsigned int fOffSet = _begin2(target) % BLOCK_SIZE;
511 unsigned int fSilencer = ~0;
512 fSilencer <<= fOffSet;
513
514 /* reset v-bitvectors */
515 std::fill(begin(VP, Standard()) + fStartBlock, begin(VP, Standard()) + fEndBlock + 1, std::numeric_limits<unsigned>::max());
516 std::fill(begin(VN, Standard()) + fStartBlock, begin(VN, Standard()) + fEndBlock + 1, 0);
517
518 /* determine start-position and start-score */
519 auto pos = _begin1(target);
520 score = (mid - _begin2(target)) * score_gap;
521 c_score[pos] = score;
522
523 /* compute with myers - forward - begin */
524 if(fSpannedBlocks == 1)
525 {
526 while (pos < _end1(target)) {
527 X = (fSilencer & forwardBitMask[(blockCount * ordValue(static_cast<TPatternAlphabet>(getValue(x,pos)))) + fStartBlock]) | VN[fStartBlock];
528
529 D0 = ((VP[fStartBlock] + (X & VP[fStartBlock])) ^ VP[fStartBlock]) | X;
530 HN = VP[fStartBlock] & D0;
531 HP = VN[fStartBlock] | ~(VP[fStartBlock] | D0);
532
533 X = (HP << 1) | (1 << fOffSet);
534 VN[fStartBlock] = X & D0;
535 VP[fStartBlock] = (HN << 1) | ~(X | D0);
536
537 if (HP & fScoreMask)
538 score--;
539 else if (HN & fScoreMask)
540 score++;
541
542 c_score[pos + 1] = score;
543
544 ++pos;
545 }
546 } /* end - short patten */
547 else
548 {
549 unsigned shift, currentBlock;
550 unsigned int temp, carryD0, carryHP, carryHN;
551
552 while (pos < _end1(target))
553 {
554 carryD0 = carryHP = carryHN = 0;
555 shift = blockCount * ordValue(static_cast<TPatternAlphabet>(getValue(x,pos)));
556
557 // computing first the top most block
558 X = (fSilencer & forwardBitMask[shift + fStartBlock]) | VN[fStartBlock];
559
560 temp = VP[fStartBlock] + (X & VP[fStartBlock]);
561 carryD0 = temp < VP[fStartBlock];
562
563 D0 = (temp ^ VP[fStartBlock]) | X;
564 HN = VP[fStartBlock] & D0;
565 HP = VN[fStartBlock] | ~(VP[fStartBlock] | D0);
566
567 X = (HP << 1) | (1 << fOffSet);
568 carryHP = HP >> (BLOCK_SIZE - 1);
569
570 VN[fStartBlock] = X & D0;
571
572 temp = (HN << 1);
573 carryHN = HN >> (BLOCK_SIZE - 1);
574
575 VP[fStartBlock] = temp | ~(X | D0);
576
577 // compute the remaining blocks
578 for (currentBlock = fStartBlock + 1; currentBlock <= fEndBlock; currentBlock++) {
579 X = forwardBitMask[shift + currentBlock] | VN[currentBlock];
580
581 temp = VP[currentBlock] + (X & VP[currentBlock]) + carryD0;
582
583 carryD0 = ((carryD0) ? temp <= VP[currentBlock] : temp < VP[currentBlock]);
584
585 D0 = (temp ^ VP[currentBlock]) | X;
586 HN = VP[currentBlock] & D0;
587 HP = VN[currentBlock] | ~(VP[currentBlock] | D0);
588
589 X = (HP << 1) | carryHP;
590 carryHP = HP >> (BLOCK_SIZE-1);
591
592 VN[currentBlock] = X & D0;
593
594 temp = (HN << 1) | carryHN;
595 carryHN = HN >> (BLOCK_SIZE - 1);
596
597 VP[currentBlock] = temp | ~(X | D0);
598 }
599
600 /* update score */
601 if (HP & fScoreMask)
602 score--;
603 else if (HN & fScoreMask)
604 score++;
605
606 c_score[pos + 1] = score;
607
608 ++pos;
609 }
610
611 } /* end - long patten */
612 /* compute with myers - forward - end */
613
614 /* compute blocks and score masks */
615 unsigned rStartBlock = (len_y - _end2(target)) / BLOCK_SIZE;
616 unsigned rEndBlock = (len_y - mid - 1) / BLOCK_SIZE;
617 unsigned rSpannedBlocks = (rEndBlock - rStartBlock) + 1;
618
619 unsigned int rScoreMask = 1 << ((len_y - mid - 1) % BLOCK_SIZE);
620 unsigned int rOffSet = (len_y - _end2(target)) % BLOCK_SIZE;
621 unsigned int rSilencer = ~0;
622 rSilencer <<= rOffSet;
623
624 /* reset v-bitvectors */
625 std::fill(begin(VP, Standard()) + rStartBlock, begin(VP, Standard()) + rEndBlock + 1, std::numeric_limits<unsigned>::max());
626 std::fill(begin(VN, Standard()) + rStartBlock, begin(VN, Standard()) + rEndBlock + 1, 0);
627
628 /* determine start-position and start-score */
629 pos = _end1(target);
630 score = (_end2(target) - mid) * score_gap;
631
632 /* set start score */
633 c_score[_end1(target)] += score;
634
635 /* determine optimal cut position -- score extension */
636 TScoreValue max = c_score[_end1(target)];
637 TScoreValue rmax = score;
638 unsigned int pos_max = _end1(target);
639
640 /* compute with myers - reverse - begin */
641 if(rSpannedBlocks == 1)
642 {
643 while (pos-- > _begin1(target)) {
644 X = (rSilencer & reverseBitMask[(blockCount * ordValue(static_cast<TPatternAlphabet>(getValue(x,pos)))) + rStartBlock]) | VN[rStartBlock];
645
646 D0 = ((VP[rStartBlock] + (X & VP[rStartBlock])) ^ VP[rStartBlock]) | X;
647 HN = VP[rStartBlock] & D0;
648 HP = VN[rStartBlock] | ~(VP[rStartBlock] | D0);
649
650 X = (HP << 1) | (1 << rOffSet);
651 VN[rStartBlock] = X & D0;
652 VP[rStartBlock] = (HN << 1) | ~(X | D0);
653
654 if (HP & rScoreMask)
655 --score;
656 else if (HN & rScoreMask)
657 ++score;
658
659 c_score[pos] += score;
660
661 /* check for optimality -- score extension */
662 if(c_score[pos]> max)
663 {
664 pos_max = pos;
665 max = c_score[pos];
666 rmax = score;
667 }
668 }
669 } /* end - short pattern */
670 else
671 {
672 unsigned shift, currentBlock;
673 unsigned int temp, carryD0, carryHP, carryHN;
674
675 while (pos-- > _begin1(target))
676 {
677 carryD0 = carryHP = carryHN = 0;
678 shift = blockCount * ordValue(static_cast<TPatternAlphabet>(getValue(x,pos)));
679
680 // compute first the top most block
681 X = (rSilencer & reverseBitMask[shift + rStartBlock]) | VN[rStartBlock];
682
683 temp = VP[rStartBlock] + (X & VP[rStartBlock]);
684 carryD0 = temp < VP[rStartBlock];
685
686 D0 = (temp ^ VP[rStartBlock]) | X;
687 HN = VP[rStartBlock] & D0;
688 HP = VN[rStartBlock] | ~(VP[rStartBlock] | D0);
689
690 X = (HP << 1) | (1 << rOffSet);
691 carryHP = HP >> (BLOCK_SIZE - 1);
692
693 VN[rStartBlock] = X & D0;
694
695 temp = (HN << 1);
696 carryHN = HN >> (BLOCK_SIZE - 1);
697
698 VP[rStartBlock] = temp | ~(X | D0);
699
700 // compute the remaining blocks
701 for (currentBlock = rStartBlock + 1; currentBlock <= rEndBlock; currentBlock++) {
702 X = reverseBitMask[shift + currentBlock] | VN[currentBlock];
703
704 temp = VP[currentBlock] + (X & VP[currentBlock]) + carryD0;
705
706 carryD0 = ((carryD0) ? temp <= VP[currentBlock] : temp < VP[currentBlock]);
707
708 D0 = (temp ^ VP[currentBlock]) | X;
709 HN = VP[currentBlock] & D0;
710 HP = VN[currentBlock] | ~(VP[currentBlock] | D0);
711
712 X = (HP << 1) | carryHP;
713 carryHP = HP >> (BLOCK_SIZE-1);
714
715 VN[currentBlock] = X & D0;
716
717 temp = (HN << 1) | carryHN;
718 carryHN = HN >> (BLOCK_SIZE - 1);
719
720 VP[currentBlock] = temp | ~(X | D0);
721 }
722
723 if (HP & rScoreMask)
724 --score;
725 else if (HN & rScoreMask)
726 ++score;
727
728 c_score[pos] += score;
729
730 /* check for optimality -- score extension*/
731 if(c_score[pos] > max)
732 {
733 pos_max = pos;
734 max = c_score[pos];
735 rmax = score;
736 }
737 }
738
739 } /* end - long pattern */
740 /* compute with myers - reverse - end */
741
742 // if computed the whole matrix max = alignment score
743 if(target == hs_complete)
744 total_score = max;
745
746 #ifdef MYERS_HIRSCHBERG_VERBOSE
747 printf("Optimal cut is at %i and %i with forward score %i and reverse score %i\n\n",mid,pos_max,(max - rmax),rmax);
748 #endif
749 /* push the two computed parts of the dp-matrix on process stack */
750 to_process.push(HirschbergSet_{pos_max, _end1(target), mid, _end2(target), rmax});
751 to_process.push(HirschbergSet_{_begin1(target), pos_max, _begin2(target), mid, max - rmax});
752
753 }
754 /* END CUT */
755 }
756
757 return total_score;
758 }
759
760 } // namespace seqan
761
762 #endif // #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_GLOBAL_ALIGNMENT_MYERS_HIRSCHBERG_IMPL_H_
763