1 // ==========================================================================
2 // SeqAn - The Library for Sequence Analysis
3 // ==========================================================================
4 // Copyright (c) 2006-2010, 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: Andreas Gogol-Doering <andreas.doering@mdc-berlin.de>
33 // ==========================================================================
34
35 #ifndef SEQAN_HEADER_GAPS_ARRAY_H
36 #define SEQAN_HEADER_GAPS_ARRAY_H
37
38 namespace SEQAN_NAMESPACE_MAIN
39 {
40
41 //////////////////////////////////////////////////////////////////////////////
42 // Tag
43
44 struct ArrayGaps;
45
46
47 //////////////////////////////////////////////////////////////////////////////
48 // Gaps - ArrayGaps Spec
49 //////////////////////////////////////////////////////////////////////////////
50
51 /**
52 .Spec.ArrayGaps:
53 ..cat:Alignments
54 ..general:Class.Gaps
55 ..summary:Stores gaps sizes in an array.
56 ..signature:Gaps<TSource, ArrayGaps >
57 ..param.TSource:Type of the ungapped sequence.
58 ...metafunction:Metafunction.Source
59 ..include:seqan/align.h
60 */
61 template <typename TSource>
62 class Gaps<TSource, ArrayGaps>
63 {
64 public:
65 typedef typename Size<Gaps>::Type TSize;
66 typedef String<TSize> TArr;
67 typedef typename Position<TSource>::Type TSourcePosition;
68 typedef typename Position<Gaps>::Type TViewPosition;
69
70 public:
71 String<TSize> data_arr; //a list of gap and non-gap region lengths
72 TViewPosition data_end_position;
73 TViewPosition data_unclipped_end_position;
74
75 Holder<TSource> data_source;
76 TSourcePosition clipped_source_begin;
77 TSourcePosition clipped_source_end;
78
79 public:
Gaps()80 Gaps():
81 data_unclipped_end_position(0),
82 clipped_source_begin(0),
83 clipped_source_end(0)
84 {
85 SEQAN_CHECKPOINT
86 }
Gaps(TSize _size)87 Gaps(TSize _size):
88 data_unclipped_end_position(0),
89 clipped_source_begin(0),
90 clipped_source_end(0)
91 {
92 _initToResize(*this, _size);
93 }
Gaps(TSource & source_)94 Gaps(TSource & source_):
95 data_unclipped_end_position(0),
96 data_source(source_),
97 clipped_source_begin(beginPosition(source_)),
98 clipped_source_end(endPosition(source_))
99 {
100 SEQAN_CHECKPOINT
101 _initToResize(*this, length(source_));
102 }
103
104 template <typename TSource2>
Gaps(TSource2 const & source_)105 Gaps(TSource2 const & source_):
106 data_unclipped_end_position(0),
107 clipped_source_begin(0),
108 clipped_source_end(length(source_))
109 //clipped_source_begin(beginPosition(source_)),
110 //clipped_source_end(endPosition(source_))
111 {
112 SEQAN_CHECKPOINT
113 data_source = source_;
114 _initToResize(*this, length(source_));
115 }
116
Gaps(Gaps const & other_)117 Gaps(Gaps const & other_):
118 data_arr(other_.data_arr),
119 data_end_position(other_.data_end_position),
120 data_unclipped_end_position(other_.data_unclipped_end_position),
121 // data_source(value(other_.data_source)), //variant: setValue => Align benutzen gemeinsame Sources
122 data_source(other_.data_source), //variant: assignValue => Align kopieren Sources
123 clipped_source_begin(other_.clipped_source_begin),
124 clipped_source_end(other_.clipped_source_end)
125 {
126 SEQAN_CHECKPOINT
127 }
128 Gaps & operator = (Gaps const & other_)
129 {
130 SEQAN_CHECKPOINT
131 data_arr = other_.data_arr;
132 data_end_position = other_.data_end_position;
133 data_unclipped_end_position = other_.data_unclipped_end_position,
134 setValue(data_source, source(other_));
135 clipped_source_begin = other_.clipped_source_begin;
136 clipped_source_end = other_.clipped_source_end;
137 return *this;
138 }
139
~Gaps()140 ~Gaps()
141 {
142 SEQAN_CHECKPOINT
143 }
144
145 inline typename Reference<Gaps>::Type
146 operator[](TViewPosition view_pos)
147 {
148 SEQAN_CHECKPOINT
149 return value(*this, view_pos);
150 }
151 inline typename Reference<Gaps const>::Type
152 operator[](TViewPosition view_pos) const
153 {
154 SEQAN_CHECKPOINT
155 return value(*this, view_pos);
156 }
157 //____________________________________________________________________________
158 };
159
160 //////////////////////////////////////////////////////////////////////////////
161
162 template <typename TSource>
163 inline String<typename Size<Gaps<TSource, ArrayGaps> >::Type> &
_dataArr(Gaps<TSource,ArrayGaps> & me)164 _dataArr(Gaps<TSource, ArrayGaps> & me)
165 {
166 SEQAN_CHECKPOINT
167 return me.data_arr;
168 }
169 template <typename TSource>
170 inline String<typename Size<Gaps<TSource, ArrayGaps> >::Type> const &
_dataArr(Gaps<TSource,ArrayGaps> const & me)171 _dataArr(Gaps<TSource, ArrayGaps> const & me)
172 {
173 SEQAN_CHECKPOINT
174 return me.data_arr;
175 }
176
177 //////////////////////////////////////////////////////////////////////////////
178
179 template <typename TSource>
180 inline Holder<TSource> &
_dataSource(Gaps<TSource,ArrayGaps> & me)181 _dataSource(Gaps<TSource, ArrayGaps> & me)
182 {
183 SEQAN_CHECKPOINT
184 return me.data_source;
185 }
186 template <typename TSource>
187 inline Holder<TSource> const &
_dataSource(Gaps<TSource,ArrayGaps> const & me)188 _dataSource(Gaps<TSource, ArrayGaps> const & me)
189 {
190 SEQAN_CHECKPOINT
191 return me.data_source;
192 }
193
194 //////////////////////////////////////////////////////////////////////////////
195
196 template <typename TSource>
197 inline typename Position<Gaps<TSource, ArrayGaps> >::Type
endPosition(Gaps<TSource,ArrayGaps> & me)198 endPosition(Gaps<TSource, ArrayGaps> & me)
199 {
200 SEQAN_CHECKPOINT
201 return me.data_end_position;
202 }
203 template <typename TSource>
204 inline typename Position<Gaps<TSource, ArrayGaps> >::Type
endPosition(Gaps<TSource,ArrayGaps> const & me)205 endPosition(Gaps<TSource, ArrayGaps> const & me)
206 {
207 SEQAN_CHECKPOINT
208 return me.data_end_position;
209 }
210
211 //////////////////////////////////////////////////////////////////////////////
212
213 template <typename TSource>
214 inline typename Position<Gaps<TSource, ArrayGaps> >::Type
_getTrailingGaps(Gaps<TSource,ArrayGaps> const & me)215 _getTrailingGaps(Gaps<TSource, ArrayGaps> const & me)
216 {
217 SEQAN_CHECKPOINT
218 return me.data_unclipped_end_position;
219 }
220
221 template <typename TSource, typename TSize>
222 inline void
_setTrailingGaps(Gaps<TSource,ArrayGaps> & me,TSize const & size)223 _setTrailingGaps(Gaps<TSource, ArrayGaps> & me, TSize const & size)
224 {
225 TSize zero = 0;
226 if (size >= zero)
227 {
228 me.data_unclipped_end_position = size;
229 }
230 else
231 {
232 me.data_unclipped_end_position = 0;
233 }
234 }
235
236 //////////////////////////////////////////////////////////////////////////////
237 /**
238 * Returns the length of the gaps sequence including the trailing gaps
239 */
240 template <typename TSource>
241 inline typename Position<Gaps<TSource, ArrayGaps> >::Type
_unclippedLength(Gaps<TSource,ArrayGaps> & me)242 _unclippedLength(Gaps<TSource, ArrayGaps> & me)
243 {
244 SEQAN_CHECKPOINT
245 return me.data_end_position + me.data_unclipped_end_position;
246 }
247
248 //////////////////////////////////////////////////////////////////////////////
249
250 template <typename TSource, typename TPosition>
251 inline void
_setEndPosition(Gaps<TSource,ArrayGaps> & me,TPosition _pos)252 _setEndPosition(Gaps<TSource, ArrayGaps> & me, TPosition _pos)
253 {
254 SEQAN_CHECKPOINT
255 me.data_end_position = _pos;
256 }
257
258 //////////////////////////////////////////////////////////////////////////////
259
260 /**
261 .Function.clippedBeginPosition:
262 ..summary:Begin position of the source segment.
263 ..cat:Alignments
264 ..signature:clippedBeginPosition(object)
265 ..param.object:An object that has a source
266 ...type:Class.Gaps
267 ..returns:The position of the first item in @Function.source.source(object)@ that is used in object.
268 ..see:Function.begin
269 ..see:Function.source
270 ..see:Function.clippedEndPosition
271 ..include:seqan/align.h
272 */
273 template <typename TSource>
274 inline typename Position<TSource>::Type
clippedBeginPosition(Gaps<TSource,ArrayGaps> const & me)275 clippedBeginPosition(Gaps<TSource, ArrayGaps> const & me)
276 {
277 SEQAN_CHECKPOINT
278 return me.clipped_source_begin;
279 }
280
281 //////////////////////////////////////////////////////////////////////////////
282
283 template <typename TSource, typename TSourcePosition>
284 inline void
_setClippedBeginPosition(Gaps<TSource,ArrayGaps> & me,TSourcePosition _pos)285 _setClippedBeginPosition(Gaps<TSource, ArrayGaps> & me, TSourcePosition _pos)
286 {
287 SEQAN_CHECKPOINT
288 me.clipped_source_begin = _pos;
289 }
290
291 //////////////////////////////////////////////////////////////////////////////
292
293 /**
294 .Function.clippedEndPosition:
295 ..summary:Position of the end of the source segment.
296 ..cat:Alignments
297 ..signature:clippedEndPosition(object)
298 ..param.object:An object that has a source
299 ...type:Class.Gaps
300 ..returns:The position behind the last element of the source segment.
301 ..see:Function.end
302 ..see:Function.sourceEnd
303 ..see:Function.clippedBeginPosition
304 ..include:seqan/align.h
305 */
306 template <typename TSource>
307 inline typename Position<TSource>::Type
clippedEndPosition(Gaps<TSource,ArrayGaps> const & me)308 clippedEndPosition(Gaps<TSource, ArrayGaps> const & me)
309 {
310 SEQAN_CHECKPOINT
311 return me.clipped_source_end;
312 }
313 //////////////////////////////////////////////////////////////////////////////
314
315 template <typename TSource, typename TSourcePosition>
316 inline void
_setClippedEndPosition(Gaps<TSource,ArrayGaps> & me,TSourcePosition _pos)317 _setClippedEndPosition(Gaps<TSource, ArrayGaps> & me, TSourcePosition _pos)
318 {
319 SEQAN_CHECKPOINT
320 me.clipped_source_end = _pos;
321 }
322
323 //////////////////////////////////////////////////////////////////////////////
324
325 template <typename TSource, typename TSize>
326 inline void
_initToResize(Gaps<TSource,ArrayGaps> & me,TSize _size)327 _initToResize(Gaps<TSource, ArrayGaps> & me,
328 TSize _size)
329 {
330 SEQAN_CHECKPOINT
331 resize(_dataArr(me), 2);
332 _dataArr(me)[0] = 0;
333 _dataArr(me)[1] = _size;
334 _setEndPosition(me, _size);
335 _setTrailingGaps(me, 0);
336 }
337
338 //////////////////////////////////////////////////////////////////////////////
339
340 template <typename TSource>
341 inline void
clearGaps(Gaps<TSource,ArrayGaps> & me)342 clearGaps(Gaps<TSource, ArrayGaps> & me)
343 {
344 SEQAN_CHECKPOINT
345 _initToResize(me, clippedEndPosition(me) - clippedBeginPosition(me));
346 }
347
348 //////////////////////////////////////////////////////////////////////////////
349
350 template <typename TSource>
351 inline void
clear(Gaps<TSource,ArrayGaps> & me)352 clear(Gaps<TSource, ArrayGaps> & me)
353 {
354 SEQAN_CHECKPOINT
355 _initToResize(me, 0);
356 _setClippedBeginPosition(me, 0);
357 _setClippedEndPosition(me, 0);
358 }
359
360 //////////////////////////////////////////////////////////////////////////////
361
362 template <typename TSource>
363 inline void
clearClipping(Gaps<TSource,ArrayGaps> & me)364 clearClipping(Gaps<TSource, ArrayGaps> & me)
365 {
366 SEQAN_CHECKPOINT
367 _setEndPosition(me, length(source(me)));
368 _setTrailingGaps(me, 0);
369 _setClippedBeginPosition(me, 0);
370 _setClippedEndPosition(me, length(source(me)));
371 }
372
373 //////////////////////////////////////////////////////////////////////////////
374
375 /**
376 .Function.toViewPosition:
377 ..summary:Transforms source to view position.
378 ..cat:Alignments
379 ..signature:toViewPosition(gaps, pos)
380 ..param.gap:A Gaps object, e.g. a row in the alignment.
381 ...type:Class.Gaps
382 ..param.pos:Position in the original sequence to get the view position for.
383 ..returns:The position in the view/gaps position.
384 ..remarks:If $gap$ is a clipped alignment row, gaps in the clipped part will influence the result. The position $pos$ is counted from the unclipped begin position and must be greater or equal the clipped begin position of $gap$.
385 ..see:Function.toSourcePosition
386 ..include:seqan/align.h
387 */
388 template <typename TSource>
389 inline typename Position< Gaps<TSource, ArrayGaps> >::Type
toViewPosition(Gaps<TSource,ArrayGaps> const & gaps,typename Position<TSource>::Type pos)390 toViewPosition(Gaps<TSource, ArrayGaps> const & gaps,
391 typename Position<TSource>::Type pos)
392 {
393 SEQAN_CHECKPOINT
394
395 SEQAN_ASSERT_GEQ(pos, clippedBeginPosition(gaps));
396 pos -= clippedBeginPosition(gaps);
397
398 typedef Gaps<TSource, ArrayGaps> const TGaps;
399 typedef typename Size<TGaps>::Type TSize;
400 typedef String<TSize> const TArr;
401
402 typename Iterator<TArr, Standard>::Type arr_begin = begin(_dataArr(gaps));
403 typename Iterator<TArr, Standard>::Type arr_end = end(_dataArr(gaps));
404 typename Position<TGaps>::Type view_pos = pos;
405 typename Position<TSource>::Type source_pos = pos;
406
407 while (true)
408 {
409 if (arr_begin == arr_end) return view_pos;
410 TSize step = *arr_begin;
411 view_pos += step;
412
413 ++arr_begin;
414
415 step = *arr_begin;
416 if (source_pos < step) return view_pos;
417 source_pos -= step;
418
419 ++arr_begin;
420 }
421 }
422
423 //____________________________________________________________________________
424
425 /**
426 .Function.toSourcePosition:
427 ..summary:Transforms view to source position, if the view position is a gap, the original position of the next non-gap entry is returned.
428 ..cat:Alignments
429 ..signature:toSourcePosition(gaps, pos)
430 ..param.gap:A Gaps object, e.g. a row in the alignment.
431 ...type:Class.Gaps
432 ..param.pos:Position in the view sequence (this includes gaps) to get the original position for.
433 ..returns:The position in the source sequence.
434 ..remarks:If $gap$ is a clipped alignment row, gaps in the clipped part will influence the result. The position $pos$ is counted from the unclipped begin position.
435 ..see:Function.toViewPosition
436 ..include:seqan/align.h
437 */
438 template <typename TSource>
439 inline typename Position<TSource>::Type
toSourcePosition(Gaps<TSource,ArrayGaps> const & gaps,typename Position<Gaps<TSource,ArrayGaps>>::Type pos)440 toSourcePosition(Gaps<TSource, ArrayGaps> const & gaps,
441 typename Position< Gaps<TSource, ArrayGaps> >::Type pos)
442 {
443 SEQAN_CHECKPOINT
444 typedef Gaps<TSource, ArrayGaps> const TGaps;
445 typedef typename Size<TGaps>::Type TSize;
446 typedef String<TSize> const TArr;
447
448 typename Iterator<TArr, Standard>::Type arr_begin = begin(_dataArr(gaps));
449 typename Iterator<TArr, Standard>::Type arr_end = end(_dataArr(gaps));
450 typename Position<TSource>::Type source_pos = clippedBeginPosition(gaps);
451 typename Position<TGaps>::Type view_pos = pos;
452
453 while (true)
454 {
455 if (arr_begin == arr_end) return source_pos;
456
457 TSize step = *arr_begin;
458 if (view_pos <= step) return source_pos;
459 view_pos -= step;
460
461 ++arr_begin;
462 step = *arr_begin;
463 if (view_pos <= step ) return source_pos + view_pos;
464 source_pos += step;
465 view_pos -= step;
466
467 ++arr_begin;
468 }
469 }
470
471 //////////////////////////////////////////////////////////////////////////////
472 // returns an iterator to view-Position
473
474 template <typename TGaps, typename TPosition>
475 inline typename Iterator<TGaps, Standard>::Type
_iteratorGapsArray(TGaps & gaps,TPosition view_position)476 _iteratorGapsArray(TGaps & gaps,
477 TPosition view_position)
478 {
479 typedef typename Size<TGaps>::Type TSize;
480 typedef typename TGaps::TArr const TArr;
481 typedef typename Source<TGaps>::Type TSource;
482 typedef typename Iterator<TArr, Standard>::Type TArrIterator;
483 typedef typename Value<TArrIterator>::Type TArrIteratorType;
484
485 TArrIterator arr_begin = begin(_dataArr(gaps), Standard());
486 TArrIterator arr_end = end(_dataArr(gaps), Standard());
487 typename Position<TGaps>::Type block = 0;
488
489 TArrIteratorType view_pos = view_position;
490
491 if (emptySource(gaps))
492 {
493 for (; arr_begin != arr_end; ++arr_begin)
494 {
495 if (view_pos < *arr_begin) break;
496 ++block;
497 view_pos -= *arr_begin;
498 }
499 return typename Iterator<TGaps, Standard>::Type(
500 gaps,
501 block,
502 view_pos);
503 }
504 else
505 {
506 typename Position<TSource>::Type source_pos = clippedBeginPosition(gaps);
507
508 while (true)
509 {
510 if ((arr_begin == arr_end) || (view_pos < *arr_begin))
511 {
512 break;
513 }
514 ++block;
515 view_pos -= *arr_begin;
516 ++arr_begin;
517
518 if (view_pos < *arr_begin)
519 {
520 source_pos += view_pos;
521 break;
522 }
523 ++block;
524 view_pos -= *arr_begin;
525 source_pos += *arr_begin;
526 ++arr_begin;
527 }
528
529 return typename Iterator<TGaps, Standard>::Type(
530 gaps,
531 iter(source(gaps), source_pos),
532 block,
533 view_pos);
534 }
535
536 }
537
538 template <typename TSource, typename TPosition, typename TTag>
539 inline typename Iterator<Gaps<TSource, ArrayGaps>, Tag<TTag> const>::Type
iter(Gaps<TSource,ArrayGaps> & gaps,TPosition view_pos,Tag<TTag> const)540 iter(Gaps<TSource, ArrayGaps> & gaps,
541 TPosition view_pos,
542 Tag<TTag> const)
543 {
544 SEQAN_CHECKPOINT
545 return _iteratorGapsArray(gaps, view_pos);
546 }
547 template <typename TSource, typename TPosition, typename TTag>
548 inline typename Iterator<Gaps<TSource, ArrayGaps> const, Tag<TTag> const>::Type
iter(Gaps<TSource,ArrayGaps> const & gaps,TPosition view_pos,Tag<TTag> const)549 iter(Gaps<TSource, ArrayGaps> const & gaps,
550 TPosition view_pos,
551 Tag<TTag> const)
552 {
553 SEQAN_CHECKPOINT
554 return _iteratorGapsArray<Gaps<TSource, ArrayGaps> const>(gaps, view_pos);
555 }
556
557 //////////////////////////////////////////////////////////////////////////////
558
559 template <typename TSource>
560 inline typename Position< Gaps<TSource, ArrayGaps> >::Type
beginPosition(Gaps<TSource,ArrayGaps> & gaps)561 beginPosition(Gaps<TSource, ArrayGaps> & gaps)
562 {
563 SEQAN_CHECKPOINT
564 return value(_dataArr(gaps), 0);
565 }
566 template <typename TSource>
567 inline typename Position< Gaps<TSource, ArrayGaps> const>::Type
beginPosition(Gaps<TSource,ArrayGaps> const & gaps)568 beginPosition(Gaps<TSource, ArrayGaps> const & gaps)
569 {
570 SEQAN_CHECKPOINT
571 return value(_dataArr(gaps), 0);
572 }
573
574 //////////////////////////////////////////////////////////////////////////////
575 //////////////////////////////////////////////////////////////////////////////
576
577
578 template <typename TSource, typename TPosition>
579 inline void
setBeginPosition(Gaps<TSource,ArrayGaps> & me,TPosition view_position)580 setBeginPosition(Gaps<TSource, ArrayGaps> & me,
581 TPosition view_position)
582 {
583 SEQAN_CHECKPOINT
584 TPosition old_pos = beginPosition(me);
585 if (length(_dataArr(me)) == 0)
586 {
587 _initToResize(me, length(source(me)));
588 }
589 _dataArr(me)[0] = view_position;
590 _setEndPosition(me, endPosition(me) + view_position - old_pos);
591 }
592
593 //////////////////////////////////////////////////////////////////////////////
594
595 template <typename TSource, typename TPosition>
596 inline void
setClippedBeginPosition(Gaps<TSource,ArrayGaps> & me,TPosition source_position)597 setClippedBeginPosition(Gaps<TSource, ArrayGaps> & me,
598 TPosition source_position)
599 {
600 SEQAN_ASSERT(length(_dataArr(me)))
601
602 typedef Gaps<TSource, ArrayGaps> TGaps;
603 typedef typename Position<TGaps>::Type TViewPosition;
604 typedef typename TGaps::TArr TArr;
605
606 TPosition old_clipped_begin_pos = clippedBeginPosition(me);
607 if (old_clipped_begin_pos == source_position) return;
608 else if (source_position < old_clipped_begin_pos)
609 {
610 SEQAN_CHECKPOINT
611 TViewPosition delta = old_clipped_begin_pos - source_position;
612 TViewPosition view_begin_pos = beginPosition(me);
613 if (view_begin_pos <= delta)
614 {
615 setBeginPosition(me, 0);
616 }
617 else
618 {
619 setBeginPosition(me, view_begin_pos - delta);
620 }
621 _setEndPosition(me, endPosition(me) + delta);
622 _dataArr(me)[1] += delta;
623 _setClippedBeginPosition(me, source_position);
624 }
625 else //(source_position > old_clipped_begin_pos)
626 {
627 SEQAN_CHECKPOINT
628 TViewPosition view_pos = toViewPosition(me, source_position);
629 TViewPosition source_pos_left = source_position - old_clipped_begin_pos;
630 TViewPosition gaps_count = 0;
631
632 typename Iterator<TArr, Standard>::Type it_arr_begin = begin(_dataArr(me));
633 typename Iterator<TArr, Standard>::Type it_arr_end = end(_dataArr(me));
634 typename Iterator<TArr, Standard>::Type it_arr = it_arr_begin;
635
636 while (it_arr != it_arr_end)
637 {
638 gaps_count += *it_arr;
639 ++it_arr;
640 if (*it_arr > source_pos_left)
641 {
642 *it_arr_begin = view_pos;
643 *it_arr -= source_pos_left;
644 if (it_arr - it_arr_begin > 1)
645 {
646 replace(_dataArr(me), 1, it_arr - it_arr_begin, "");
647 }
648 _setClippedBeginPosition(me, source_position);
649 return;
650 }
651 gaps_count += *it_arr;
652 source_pos_left -= *it_arr;
653 ++it_arr;
654 }
655 //alignment is empty
656 clear(me);
657 }
658 }
659
660 //////////////////////////////////////////////////////////////////////////////
661
662 template <typename TSource, typename TPosition>
663 inline void
setClippedEndPosition(Gaps<TSource,ArrayGaps> & me,TPosition source_position)664 setClippedEndPosition(Gaps<TSource, ArrayGaps> & me,
665 TPosition source_position)
666 {
667 typedef Gaps<TSource, ArrayGaps> TGaps;
668 typedef typename Position<TGaps>::Type TViewPosition;
669 typedef typename TGaps::TArr TArr;
670
671 TArr arr = _dataArr(me);
672 SEQAN_ASSERT(length(arr));
673
674 TPosition old_end_begin_pos = clippedEndPosition(me);
675 if (old_end_begin_pos == source_position) return;
676 else if (source_position < old_end_begin_pos)
677 {
678 SEQAN_CHECKPOINT
679 typename Iterator<TArr, Standard>::Type it_arr_begin = begin(_dataArr(me));
680 typename Iterator<TArr, Standard>::Type it_arr_end = end(_dataArr(me));
681 typename Iterator<TArr, Standard>::Type it_arr = it_arr_begin;
682 TViewPosition end_pos = 0;
683 TViewPosition chars_to_scan = source_position - clippedBeginPosition(me);
684
685 while (it_arr != it_arr_end)
686 {
687 end_pos += *it_arr;
688 ++it_arr;
689 if (*it_arr >= chars_to_scan)
690 {
691 resize(_dataArr(me), it_arr - it_arr_begin + 1);
692 *it_arr = chars_to_scan;
693 _setEndPosition(me, end_pos + chars_to_scan);
694 break;
695 }
696 end_pos += *it_arr;
697 chars_to_scan -= *it_arr;
698 ++it_arr;
699 }
700 }
701 else //(source_position > old_end_begin_pos)
702 {
703 SEQAN_CHECKPOINT
704 *(end(_dataArr(me)) - 1) += (source_position - old_end_begin_pos);
705 _setEndPosition(me, endPosition(me) + source_position - old_end_begin_pos);
706 }
707 _setClippedEndPosition(me, source_position);
708 }
709
710 //////////////////////////////////////////////////////////////////////////////
711
712 template <typename TSource>
713 inline typename Size<TSource>::Type
sourceLength(Gaps<TSource,ArrayGaps> & me)714 sourceLength(Gaps<TSource, ArrayGaps> & me)
715 {
716 SEQAN_CHECKPOINT
717 return clippedEndPosition(me) - clippedBeginPosition(me);
718 }
719
720 //////////////////////////////////////////////////////////////////////////////
721 //////////////////////////////////////////////////////////////////////////////
722 //////////////////////////////////////////////////////////////////////////////
723 // Gaps Iterator
724 //////////////////////////////////////////////////////////////////////////////
725
726
727 template <typename TGaps>
728 class Iter<TGaps, GapsIterator<ArrayGaps> >
729 {
730 public:
731 typedef typename Size<TGaps>::Type TGapsSize;
732 typedef String<TGapsSize> TArr;
733 typedef typename Position<TArr>::Type TArrPosition;
734 typedef typename Source<Iter>::Type TSourceIterator;
735
736 TGaps * data_container; //the gaps object
737 mutable TSourceIterator data_source;
738 mutable TArrPosition data_block; //block in array of container
739 mutable TGapsSize data_sub; //position block
740
741 public:
Iter()742 Iter()
743 {
744 SEQAN_CHECKPOINT
745 }
Iter(Iter const & other_)746 Iter(Iter const & other_):
747 data_container(other_.data_container),
748 data_source(other_.data_source),
749 data_block(other_.data_block),
750 data_sub(other_.data_sub)
751 {
752 SEQAN_CHECKPOINT
753 }
754 Iter(TGaps & container_,
755 TSourceIterator source_iterator,
756 TArrPosition block_ = 0,
757 TGapsSize sub_ = 0):
758 data_container(&container_),
759 data_source(source_iterator),
760 data_block(block_),
761 data_sub(sub_)
762 {
763 SEQAN_CHECKPOINT
764 }
765
766 Iter(TGaps & container_,
767 TArrPosition block_,
768 TGapsSize sub_ = 0):
769 data_container(&container_),
770 data_block(block_),
771 data_sub(sub_)
772 {
773 SEQAN_CHECKPOINT
774 }
775
~Iter()776 ~Iter()
777 {
778 SEQAN_CHECKPOINT
779 }
780
781 Iter const & operator = (Iter const & other_)
782 {
783 SEQAN_CHECKPOINT
784 data_container = other_.data_container;
785 data_block = other_.data_block;
786 data_sub = other_.data_sub;
787 data_source = other_.data_source;
788 return *this;
789 }
790 };
791
792 //////////////////////////////////////////////////////////////////////////////
793
794 template <typename TGaps>
795 inline typename GetValue< Iter<TGaps, GapsIterator<ArrayGaps> > >::Type
getValue(Iter<TGaps,GapsIterator<ArrayGaps>> & me)796 getValue(Iter<TGaps, GapsIterator<ArrayGaps> > & me)
797 {
798 SEQAN_CHECKPOINT
799 typedef typename Value<Iter<TGaps, GapsIterator<ArrayGaps> > >::Type TValue;
800 if (isGap(me)) return gapValue<TValue>();
801 else return getValue(source(me));
802 }
803 template <typename TGaps>
804 inline typename GetValue< Iter<TGaps, GapsIterator<ArrayGaps> > const>::Type
getValue(Iter<TGaps,GapsIterator<ArrayGaps>> const & me)805 getValue(Iter<TGaps, GapsIterator<ArrayGaps> > const & me)
806 {
807 SEQAN_CHECKPOINT
808 typedef typename Value<Iter<TGaps, GapsIterator<ArrayGaps> > const>::Type TValue;
809 if (isGap(me)) return gapValue<TValue>();
810 else return getValue(source(me));
811 }
812
813 //____________________________________________________________________________
814
815 template <typename TGaps>
816 inline bool
isGap(Iter<TGaps,GapsIterator<ArrayGaps>> const & me)817 isGap(Iter<TGaps, GapsIterator<ArrayGaps> > const & me)
818 {
819 SEQAN_CHECKPOINT
820 return !(me.data_block & 1);
821 }
822
823 //____________________________________________________________________________
824
825 template <typename TGaps>
826 inline typename Size<TGaps>::Type
countGaps(Iter<TGaps,GapsIterator<ArrayGaps>> const & me)827 countGaps(Iter<TGaps, GapsIterator<ArrayGaps> > const & me)
828 {
829 if (isGap(me))
830 {
831 typedef typename Size<TGaps>::Type TGapsSize;
832 typedef String<TGapsSize> const TArr;
833
834 TArr & arr = _dataArr(container(me));
835 typename Position<TArr>::Type pos = me.data_block;
836 typename Position<TArr>::Type block_pos = me.data_sub;
837
838 if (length(arr) == pos)
839 {//counting trailing gaps
840 if (block_pos <= _getTrailingGaps(container(me)))
841 {//is in range of trailing gaps
842 SEQAN_CHECKPOINT
843 return _getTrailingGaps(container(me)) - block_pos;
844 }
845 else
846 {//no trailing gaps here
847 SEQAN_CHECKPOINT
848 return 0;
849 }
850 }
851 else
852 {
853 SEQAN_CHECKPOINT
854 return arr[pos] - me.data_sub;
855 }
856 }
857 else
858 {
859 SEQAN_CHECKPOINT
860 return 0;
861 }
862 }
863
864 template <typename TGaps>
865 inline typename Size<TGaps>::Type
countCharacters(Iter<TGaps,GapsIterator<ArrayGaps>> const & me)866 countCharacters(Iter<TGaps, GapsIterator<ArrayGaps> > const & me)
867 {
868 if (isGap(me))
869 {// only count chars
870 SEQAN_CHECKPOINT
871 return 0;
872 }
873 else
874 {//count chars
875 SEQAN_CHECKPOINT
876 typedef typename Size<TGaps>::Type TGapsSize;
877 typedef String<TGapsSize> const TArr;
878
879 TArr & arr = _dataArr(container(me));
880 typename Position<TArr>::Type pos = me.data_block;
881 typename Position<TArr>::Type block_pos = me.data_sub;
882
883 return arr[pos] - block_pos;
884 }
885 }
886
887 //____________________________________________________________________________
888
889 template <typename T>
890 inline void
_goNextArrayGapsIterator(T const & me)891 _goNextArrayGapsIterator(T const & me)
892 {
893 if (!isGap(me))
894 {
895 SEQAN_CHECKPOINT
896 goNext(me.data_source);
897 }
898
899 ++me.data_sub;
900 if ((me.data_block < length(_dataArr(container(me)))) && (me.data_sub >= _dataArr(container(me))[me.data_block]))
901 {
902 SEQAN_CHECKPOINT
903 ++me.data_block;
904 me.data_sub = 0;
905 }
906 }
907 template <typename TGaps>
908 inline void
goNext(Iter<TGaps,GapsIterator<ArrayGaps>> & me)909 goNext(Iter<TGaps, GapsIterator<ArrayGaps> > & me)
910 {
911 _goNextArrayGapsIterator(me);
912 }
913 template <typename TGaps>
914 inline void
goNext(Iter<TGaps,GapsIterator<ArrayGaps>> const & me)915 goNext(Iter<TGaps, GapsIterator<ArrayGaps> > const & me)
916 {
917 _goNextArrayGapsIterator(me);
918 }
919
920 //____________________________________________________________________________
921
922 template <typename T>
923 inline void
_goPreviousArrayGapsIterator(T const & me)924 _goPreviousArrayGapsIterator(T const & me)
925 {
926 if (me.data_sub > 0)
927 {
928 SEQAN_CHECKPOINT
929 --me.data_sub;
930 }
931 else
932 {
933 if (me.data_block > 0)
934 {
935 SEQAN_CHECKPOINT
936 --me.data_block;
937 me.data_sub = _dataArr(container(me))[me.data_block] - 1;
938 }
939 }
940 if (!isGap(me))
941 {
942 SEQAN_CHECKPOINT
943 goPrevious(me.data_source);
944 }
945 }
946 template <typename TGaps>
947 inline void
goPrevious(Iter<TGaps,GapsIterator<ArrayGaps>> & me)948 goPrevious(Iter<TGaps, GapsIterator<ArrayGaps> > & me)
949 {
950 _goPreviousArrayGapsIterator(me);
951 }
952 template <typename TGaps>
953 inline void
goPrevious(Iter<TGaps,GapsIterator<ArrayGaps>> const & me)954 goPrevious(Iter<TGaps, GapsIterator<ArrayGaps> > const & me)
955 {
956 _goPreviousArrayGapsIterator(me);
957 }
958 //____________________________________________________________________________
959
960 template <typename TGaps>
961 inline bool
atBegin(Iter<TGaps,GapsIterator<ArrayGaps>> const & me)962 atBegin(Iter<TGaps, GapsIterator<ArrayGaps> > const & me)
963 {
964 SEQAN_CHECKPOINT
965 return ((me.data_block == 1) && (me.data_sub == 0));
966 }
967
968 //____________________________________________________________________________
969
970 template <typename TGaps>
971 inline bool
atEnd(Iter<TGaps,GapsIterator<ArrayGaps>> const & me)972 atEnd(Iter<TGaps, GapsIterator<ArrayGaps> > const & me)
973 {
974 SEQAN_CHECKPOINT
975 return ((me.data_block == length(_dataArr(container(me)))) && (me.data_sub == 0));
976 }
977
978 //____________________________________________________________________________
979
980 template <typename TGaps, typename TCount>
981 inline void
insertGaps(Iter<TGaps,GapsIterator<ArrayGaps>> const & me,TCount size)982 insertGaps(Iter<TGaps, GapsIterator<ArrayGaps> > const & me,
983 TCount size)
984 {
985 typedef typename Size<TGaps>::Type TGapsSize;
986 typedef String<TGapsSize> TArr;
987
988 TArr & arr = _dataArr(container(me));
989 typename Position<TArr>::Type pos = me.data_block;
990 typename Position<TArr>::Type block_pos = me.data_sub;
991 typename Iterator<TArr, Rooted>::Type it = iter(arr, pos);
992
993 if (isGap(me))
994 {
995 if (pos < length(arr))
996 {//here is a gap: expand it
997 SEQAN_CHECKPOINT
998 value(it) += size;
999 }
1000 else if (pos == length(arr))
1001 {//gap at end: add trailing gaps
1002 if (block_pos <= _getTrailingGaps(container(me)))
1003 {
1004 SEQAN_CHECKPOINT
1005 _setTrailingGaps(container(me), _getTrailingGaps(container(me)) + size);
1006 }
1007 return; // do not adjust right position
1008 }
1009 }
1010 else
1011 {
1012 if (me.data_sub)
1013 {//here is no gap: insert one
1014 SEQAN_CHECKPOINT
1015 _clearSpace(arr, 2, pos + 1, pos + 1, Generous());
1016
1017 it = iter(arr, pos); //reload, since iterator could be invalid
1018 it[2] = it[0] - me.data_sub;
1019 it[0] = me.data_sub;
1020 it[1] = size;
1021
1022 ++me.data_block;
1023 me.data_sub = 0;
1024 }
1025 else
1026 {//insert gaps at begin of a char block
1027 SEQAN_CHECKPOINT
1028 me.data_sub = arr[pos - 1];
1029 arr[pos - 1] += size; //note: pos > 0 because this is a char block
1030 --me.data_block;
1031 }
1032 }
1033
1034 //adjust right position
1035 _setEndPosition(container(me), endPosition(container(me)) + size);
1036 }
1037
1038 //____________________________________________________________________________
1039
1040 //delete up to size gaps
1041
1042 template <typename TGaps, typename TCount>
1043 inline void
removeGaps(Iter<TGaps,GapsIterator<ArrayGaps>> const & me,TCount _size)1044 removeGaps(Iter<TGaps, GapsIterator<ArrayGaps> > const & me,
1045 TCount _size)
1046 {
1047 typedef typename Size<TGaps>::Type TGapsSize;
1048 typedef String<TGapsSize> TArr;
1049
1050 TGapsSize size = _size;
1051
1052 if (isGap(me))
1053 {//here is a gap
1054 TArr & arr = _dataArr(container(me));
1055 if (me.data_block < length(arr))
1056 {//not trailing gaps
1057 typename Iterator<TArr, Standard>::Type it = iter(arr, me.data_block);
1058
1059 TGapsSize rest_size = *it - me.data_sub;
1060
1061 if (rest_size <= size)
1062 {//remove rest of gap block
1063 if (me.data_sub || !me.data_block)
1064 {//keep this gap block
1065 SEQAN_CHECKPOINT
1066 *it = me.data_sub;
1067 ++me.data_block;
1068 me.data_sub = 0;
1069 }
1070 else
1071 {//remove complete gap block, merge two char blocks
1072 SEQAN_CHECKPOINT
1073 TGapsSize data_sub = *(it-1); //(note: this well defined, since "it" was not the first block)
1074
1075 *(it-1) += *(it+1);
1076 replace(_dataArr(container(me)), me.data_block, me.data_block + 2, "");
1077
1078 --me.data_block;
1079 me.data_sub = data_sub;
1080 }
1081 _setEndPosition(container(me), endPosition(container(me)) - rest_size);
1082 }
1083 else
1084 {//remove a part of this block
1085 SEQAN_CHECKPOINT
1086 *it -= size;
1087 _setEndPosition(container(me), endPosition(container(me)) - size);
1088 }
1089 }
1090 else if (me.data_block == length(arr))
1091 {
1092 if (me.data_sub <= _getTrailingGaps(container(me)))
1093 {//remove trailing gaps
1094 if (size > (_getTrailingGaps(container(me)) - me.data_sub))
1095 {//remove more than exists
1096 SEQAN_CHECKPOINT
1097 _setTrailingGaps(container(me), me.data_sub);
1098 }
1099 else
1100 {//remove less or equal than exists
1101 SEQAN_CHECKPOINT
1102 _setTrailingGaps(container(me), _getTrailingGaps(container(me)) - size);
1103 }
1104 }
1105 }//else: trailing gaps: infinite gaps here - nothing to do
1106 }
1107 //else: here is no gap - nothing to do
1108 }
1109
1110 //____________________________________________________________________________
1111
1112 template <typename TGaps1, typename TGaps2>
1113 inline bool
1114 operator == (Iter<TGaps1, GapsIterator<ArrayGaps> > & _left,
1115 Iter<TGaps2, GapsIterator<ArrayGaps> > & _right)
1116 {
1117 SEQAN_CHECKPOINT
1118 return (_left.data_block == _right.data_block) && (_left.data_sub == _right.data_sub);
1119 }
1120 template <typename TGaps1, typename TGaps2>
1121 inline bool
1122 operator == (Iter<TGaps1, GapsIterator<ArrayGaps> > const & _left,
1123 Iter<TGaps2, GapsIterator<ArrayGaps> > & _right)
1124 {
1125 SEQAN_CHECKPOINT
1126 return (_left.data_block == _right.data_block) && (_left.data_sub == _right.data_sub);
1127 }
1128 template <typename TGaps1, typename TGaps2>
1129 inline bool
1130 operator == (Iter<TGaps1, GapsIterator<ArrayGaps> > & _left,
1131 Iter<TGaps2, GapsIterator<ArrayGaps> > const & _right)
1132 {
1133 SEQAN_CHECKPOINT
1134 return (_left.data_block == _right.data_block) && (_left.data_sub == _right.data_sub);
1135 }
1136 template <typename TGaps1, typename TGaps2>
1137 inline bool
1138 operator == (Iter<TGaps1, GapsIterator<ArrayGaps> > const & _left,
1139 Iter<TGaps2, GapsIterator<ArrayGaps> > const & _right)
1140 {
1141 SEQAN_CHECKPOINT
1142 return (_left.data_block == _right.data_block) && (_left.data_sub == _right.data_sub);
1143 }
1144 //____________________________________________________________________________
1145
1146 template <typename TGaps1, typename TGaps2>
1147 inline bool
1148 operator != (Iter<TGaps1, GapsIterator<ArrayGaps> > & _left,
1149 Iter<TGaps2, GapsIterator<ArrayGaps> > & _right)
1150 {
1151 SEQAN_CHECKPOINT
1152 return (_left.data_block != _right.data_block) || (_left.data_sub != _right.data_sub);
1153 }
1154 template <typename TGaps1, typename TGaps2>
1155 inline bool
1156 operator != (Iter<TGaps1, GapsIterator<ArrayGaps> > const & _left,
1157 Iter<TGaps2, GapsIterator<ArrayGaps> > & _right)
1158 {
1159 SEQAN_CHECKPOINT
1160 return (_left.data_block != _right.data_block) || (_left.data_sub != _right.data_sub);
1161 }
1162 template <typename TGaps1, typename TGaps2>
1163 inline bool
1164 operator != (Iter<TGaps1, GapsIterator<ArrayGaps> > & _left,
1165 Iter<TGaps2, GapsIterator<ArrayGaps> > const & _right)
1166 {
1167 SEQAN_CHECKPOINT
1168 return (_left.data_block != _right.data_block) || (_left.data_sub != _right.data_sub);
1169 }
1170 template <typename TGaps1, typename TGaps2>
1171 inline bool
1172 operator != (Iter<TGaps1, GapsIterator<ArrayGaps> > const & _left,
1173 Iter<TGaps2, GapsIterator<ArrayGaps> > const & _right)
1174 {
1175 SEQAN_CHECKPOINT
1176 return (_left.data_block != _right.data_block) || (_left.data_sub != _right.data_sub);
1177 }
1178
1179
1180 //////////////////////////////////////////////////////////////////////////////
1181 }// namespace SEQAN_NAMESPACE_MAIN
1182
1183 //////////////////////////////////////////////////////////////////////////////
1184 #endif //#ifndef SEQAN_HEADER_...
1185