1 // ==========================================================================
2 //                 SeqAn - The Library for Sequence Analysis
3 // ==========================================================================
4 // Copyright (c) 2006-2015, 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: Manuel Holtgrewe <manuel.holtgrewe@fu-berlin.de>
33 // ==========================================================================
34 
35 #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_GAPS_ITERATOR_ANCHOR_H_
36 #define SEQAN_INCLUDE_SEQAN_ALIGN_GAPS_ITERATOR_ANCHOR_H_
37 
38 namespace seqan {
39 
40 // ============================================================================
41 // Forwards
42 // ============================================================================
43 
44 // ============================================================================
45 // Tags, Classes, Enums
46 // ============================================================================
47 
48 // Iterator for AnchorGaps objects.
49 
50 template <typename TGaps_, typename TGapAnchors_> //Gaps<TSource_, AnchorGaps<TGapAnchors_> >
51 class Iter<TGaps_, GapsIterator<AnchorGaps<TGapAnchors_> > >
52 {
53 public:
54     typedef TGaps_                                            TGaps;
55     typedef typename Source<TGaps>::Type                    TSource;
56     typedef TGapAnchors_                                    TGapAnchors;
57 
58     // TODO(holtgrew): Why is the following commented out?
59 //    typedef typename Value<TGapAnchors>::Type                TGapAnchor;
60     typedef typename Size<typename Value<TGapAnchors>::Type>::Type          TGapAnchorSize_;
61     typedef GapAnchor<typename MakeSigned_<TGapAnchorSize_>::Type>          TGapAnchor;
62     typedef typename MakeSigned<typename Position<TGapAnchor>::Type>::Type  TGapPos;
63     typedef typename Iterator<TGapAnchors, Standard>::Type                  TAnchorIter;
64 
65     TGaps *                    data_container;                            //the gaps object
66     TGapPos                 seqLength;
67     mutable TGapAnchor        current;
68     mutable TGapAnchor        prevAnchor;
69     mutable TGapAnchor        nextAnchor;
70     mutable TGapAnchor        viewBegin;
71     mutable TGapAnchor        viewEnd;
72     mutable int                anchorIdx;
73 
74 public:
Iter()75     Iter()
76     {
77 SEQAN_CHECKPOINT
78         data_container = NULL;
79         seqLength = 0;
80     }
81 /*    Iter(Iter const & other_):
82         data_container(other_.data_container),
83         seqLength(other_.seqLength),
84         current(other_.current),
85         prevAnchor(other_.prevAnchor),
86         nextAnchor(other_.nextAnchor),
87         anchorIdx(other_.anchorIdx)
88     {
89 SEQAN_CHECKPOINT
90     }
Iter(TGaps & container_)91 */    Iter(TGaps & container_):
92         data_container(&container_)
93     {
94 SEQAN_CHECKPOINT
95         _assignSourceLength(seqLength, container_);
96         _goToGapAnchorIterator(*this, data_container->data_viewCutBegin + data_container->data_cutBegin);
97         viewBegin = current;
98         viewEnd.gapPos   = _unclippedLength(*data_container) + data_container->data_cutBegin - data_container->data_viewCutEnd;
99         viewEnd.seqPos = positionGapToSeq(*data_container, viewEnd.gapPos);
100     }
Iter(TGaps & container_,TGapPos clippedViewPosition)101     Iter(TGaps & container_, TGapPos clippedViewPosition):
102         data_container(&container_)
103     {
104 SEQAN_CHECKPOINT
105         _assignSourceLength(seqLength, container_);
106         _goToGapAnchorIterator(*this, clippedViewPosition + data_container->data_viewCutBegin + data_container->data_cutBegin);
107         viewBegin.gapPos = data_container->data_viewCutBegin + data_container->data_cutBegin;
108         viewEnd.gapPos   = _unclippedLength(*data_container) + data_container->data_cutBegin - data_container->data_viewCutEnd;
109         viewBegin.seqPos = positionGapToSeq(*data_container, viewBegin.gapPos);
110         viewEnd.seqPos   = positionGapToSeq(*data_container, viewEnd.gapPos);
111     }
~Iter()112     ~Iter()
113     {
114 SEQAN_CHECKPOINT
115     }
116 
117     Iter const & operator = (Iter const & other_)
118     {
119 SEQAN_CHECKPOINT
120         data_container = other_.data_container;
121         seqLength = other_.seqLength;
122         current = other_.current;
123         prevAnchor = other_.prevAnchor;
124         nextAnchor = other_.nextAnchor;
125         anchorIdx = other_.anchorIdx;
126         viewBegin = other_.viewBegin;
127         viewEnd = other_.viewEnd;
128         return *this;
129     }
130 };
131 
132 // ============================================================================
133 // Metafunctions
134 // ============================================================================
135 
136 // ============================================================================
137 // Functions
138 // ============================================================================
139 
140 // ----------------------------------------------------------------------------
141 // Function container()
142 // ----------------------------------------------------------------------------
143 
144 // TODO(holtgrew): Can go if data_container were _container, dupe in _base.h
145 
146 template <typename TGaps, typename TGapsArray>
147 inline TGaps &
container(Iter<TGaps,GapsIterator<AnchorGaps<TGapsArray>>> & me)148 container(Iter<TGaps, GapsIterator<AnchorGaps<TGapsArray> > > & me)
149 {
150     return *me.data_container;
151 }
152 
153 template <typename TGaps, typename TGapsArray>
154 inline TGaps &
container(Iter<TGaps,GapsIterator<AnchorGaps<TGapsArray>>> const & me)155 container(Iter<TGaps, GapsIterator<AnchorGaps<TGapsArray> > > const & me)
156 {
157     return *me.data_container;
158 }
159 
160 // ----------------------------------------------------------------------------
161 // Function source()
162 // ----------------------------------------------------------------------------
163 
164 template <typename TGaps, typename TGapAnchors>
165 inline typename Source<Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const>::Type
source(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> & me)166 source(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > & me)
167 {
168     return begin(source(*me.data_container), Rooted()) + me.current.seqPos;
169 }
170 
171 template <typename TGaps, typename TGapAnchors>
172 inline typename Source<Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > >::Type
source(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me)173 source(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me)
174 {
175     return begin(source(*me.data_container), Rooted()) + me.current.seqPos;
176 }
177 
178 // ----------------------------------------------------------------------------
179 // Function getValue()
180 // ----------------------------------------------------------------------------
181 
182 template <typename TGaps, typename TGapAnchors>
183 inline typename GetValue< Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > >::Type
getValue(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> & me)184 getValue(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > & me)
185 {
186     typedef typename Value<Iter<TGaps, GapsIterator<ArrayGaps> > >::Type TValue;
187     if (isGap(me)) return gapValue<TValue>();
188     else if (isUnknown(me)) return unknownValue<TValue>();
189     else return getValue(source(me));
190 }
191 
192 template <typename TGaps, typename TGapAnchors>
193 inline typename GetValue< Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const>::Type
getValue(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me)194 getValue(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me)
195 {
196     typedef typename Value<Iter<TGaps, GapsIterator<ArrayGaps> > const>::Type TValue;
197     if (isGap(me)) return gapValue<TValue>();
198     else if (isUnknown(me)) return unknownValue<TValue>();
199     else return getValue(source(me));
200 }
201 
202 // ----------------------------------------------------------------------------
203 // Function value()
204 // ----------------------------------------------------------------------------
205 
206 template <typename TGaps, typename TGapAnchors>
207 inline typename Reference< Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > >::Type
value(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> & it)208 value(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > & it)
209 {
210     typedef typename Reference<Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > >::Type TProxy;
211     return TProxy(it);
212 }
213 
214 template <typename TGaps, typename TGapAnchors>
215 inline typename Reference< Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const>::Type
value(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & it)216 value(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & it)
217 {
218     typedef typename Reference<Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const>::Type TProxy;
219     return TProxy(it);
220 }
221 
222 // ----------------------------------------------------------------------------
223 // Function isGap()
224 // ----------------------------------------------------------------------------
225 
226 template <typename TGaps, typename TGapAnchors>
227 inline bool
isGap(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me)228 isGap(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me)
229 {
230     return me.current.seqPos == me.nextAnchor.seqPos;
231 }
232 
233 // ----------------------------------------------------------------------------
234 // Function isUnknown()
235 // ----------------------------------------------------------------------------
236 
237 template <typename TGaps, typename TGapAnchors>
238 inline bool
isUnknown(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me)239 isUnknown(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me)
240 {
241     int len;
242     _assignSourceLength(len, *me.data_container);
243     return me.current.seqPos < 0 || me.current.seqPos >= len;
244 }
245 
246 // ----------------------------------------------------------------------------
247 // Function isClipped()
248 // ----------------------------------------------------------------------------
249 
250 template <typename TGaps, typename TGapAnchors>
251 inline bool
isClipped(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me)252 isClipped(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me)
253 {
254     return me.current.gapPos == me.nextAnchor.gapPos;
255 }
256 
257 // ----------------------------------------------------------------------------
258 // Function countGaps()
259 // ----------------------------------------------------------------------------
260 
261 template <typename TGaps, typename TGapAnchors>
262 inline typename Size<TGaps>::Type
countGaps(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me)263 countGaps(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me)
264 {
265     if (!isGap(me))
266         return 0;
267     if (me.nextAnchor.gapPos > me.viewEnd.gapPos)
268         return me.viewEnd.gapPos - me.current.gapPos;
269     return me.nextAnchor.gapPos - me.current.gapPos;
270 }
271 
272 // ----------------------------------------------------------------------------
273 // Function countCharacters()
274 // ----------------------------------------------------------------------------
275 
276 template <typename TGaps, typename TGapAnchors>
277 inline typename Size<TGaps>::Type
countCharacters(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me)278 countCharacters(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me)
279 {
280     if (isGap(me))
281         return 0;
282     if (me.nextAnchor.seqPos > me.viewEnd.seqPos)
283         return me.viewEnd.seqPos - me.current.seqPos;
284     return me.nextAnchor.seqPos - me.current.seqPos;
285 }
286 
287 // ----------------------------------------------------------------------------
288 // Function blockLength()
289 // ----------------------------------------------------------------------------
290 
291 template <typename TGaps, typename TGapAnchors>
292 inline typename Size<TGaps>::Type
blockLength(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> & me)293 blockLength(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > & me)
294 {
295     if (isGap(me))
296         return countGaps(me);
297     else
298         return countCharacters(me);
299 }
300 
301 // ----------------------------------------------------------------------------
302 // Function atBegin()
303 // ----------------------------------------------------------------------------
304 
305 template <typename TGaps, typename TGapAnchors>
306 inline bool
atBegin(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> & me)307 atBegin(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > & me)
308 {
309 //    return me.current.seqPos == 0 && me.current.gapPos == 0;
310     return me.current <= me.viewBegin;
311 }
312 
313 template <typename TGaps, typename TGapAnchors>
314 inline bool
atBegin(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me)315 atBegin(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me)
316 {
317 //    return me.current.seqPos == 0 && me.current.gapPos == 0;
318     return me.current <= me.viewBegin;
319 }
320 
321 // ----------------------------------------------------------------------------
322 // Function atEnd()
323 // ----------------------------------------------------------------------------
324 
325 template <typename TGaps, typename TGapAnchors>
326 inline bool
atEnd(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> & me)327 atEnd(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > & me)
328 {
329 //    return me.current == me.nextAnchor;
330     return me.current >= me.viewEnd;
331 }
332 
333 template <typename TGaps, typename TGapAnchors>
334 inline bool
atEnd(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me)335 atEnd(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me)
336 {
337 //    return me.current == me.nextAnchor;
338     return me.current >= me.viewEnd;
339 }
340 
341 // ----------------------------------------------------------------------------
342 // Function operator==()
343 // ----------------------------------------------------------------------------
344 
345 template <typename TGaps, typename TGapAnchors>
346 inline bool
347 operator == (
348     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & left,
349     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & right)
350 {
351     return left.current == right.current;
352 }
353 
354 // ----------------------------------------------------------------------------
355 // Function operator!=()
356 // ----------------------------------------------------------------------------
357 
358 template <typename TGaps, typename TGapAnchors>
359 inline bool
360 operator != (
361     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & left,
362     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & right)
363 {
364     return left.current != right.current;
365 }
366 
367 // ----------------------------------------------------------------------------
368 // Function operator<()
369 // ----------------------------------------------------------------------------
370 
371 template <typename TGaps, typename TGapAnchors>
372 inline bool
373 operator < (
374     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & left,
375     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & right)
376 {
377     return left.current < right.current;
378 }
379 
380 // ----------------------------------------------------------------------------
381 // Function operator<=()
382 // ----------------------------------------------------------------------------
383 
384 template <typename TGaps, typename TGapAnchors>
385 inline bool
386 operator<=(
387     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & left,
388     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & right)
389 {
390     return !(left.current > right.current);
391 }
392 
393 // ----------------------------------------------------------------------------
394 // Function operator>()
395 // ----------------------------------------------------------------------------
396 
397 template <typename TGaps, typename TGapAnchors>
398 inline bool
399 operator > (
400     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & left,
401     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & right)
402 {
403     return left.current > right.current;
404 }
405 
406 // ----------------------------------------------------------------------------
407 // Function operator>=()
408 // ----------------------------------------------------------------------------
409 
410 template <typename TGaps, typename TGapAnchors>
411 inline bool
412 operator>=(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & lhs,
413            Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & rhs)
414 {
415     return !(lhs < rhs);
416 }
417 
418 // ----------------------------------------------------------------------------
419 // Function insertGaps()
420 // ----------------------------------------------------------------------------
421 
422 template <typename TGaps, typename TGapAnchors, typename TCount>
423 inline void
insertGaps(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & me,TCount size)424 insertGaps(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & me,
425            TCount size)
426 {
427     TGapAnchors & anchors = _dataAnchors(*me.data_container);
428     typedef typename Iterator<TGapAnchors, Standard>::Type TIter;
429 
430     if (size <= 0) return;
431 
432     // insert a new anchor
433     if (!isGap(me))
434     {
435         if (me.prevAnchor.gapPos == me.current.gapPos)
436         {
437             me.nextAnchor = me.prevAnchor;
438             _getAnchor(me.prevAnchor, *me.data_container, --me.anchorIdx);
439         }
440         else
441         {
442             me.nextAnchor = me.current;
443             insertValue(anchors, me.anchorIdx, me.nextAnchor, Generous());
444         }
445     }
446     else
447     {
448         if (me.anchorIdx >= (int)length(anchors))
449         {
450             // add gap after the sequence and in (or at the right boundary of) the view
451             if (me.current.gapPos <= me.viewEnd.gapPos)
452             {
453                 container(me).data_cutEnd -= size;
454                 me.viewEnd.gapPos += size;
455             }
456             return;
457         }
458         if (empty(anchors))
459             appendValue(anchors, me.nextAnchor, Generous());
460     }
461     if (me.anchorIdx < (int)length(anchors))
462     {
463         if (me.anchorIdx >= 0)
464         {
465             me.nextAnchor.gapPos += size;
466             TIter it = begin(anchors, Standard());
467             TIter itEnd = end(anchors, Standard());
468             if (me.anchorIdx >= 0)
469                 it += me.anchorIdx;
470             for (; it != itEnd; ++it)
471                 (*it).gapPos += size;
472         }
473         else
474             // add gap before the sequence and in (or at the left boundary of) the view
475             if (me.current.gapPos >= me.viewBegin.gapPos)
476             {
477                 container(me).data_cutBegin -= size;
478                 me.viewBegin.gapPos -= size;
479                 me.current.gapPos -= size;
480                 return;
481             }
482     }
483     if (me.current.gapPos <= me.viewEnd.gapPos)
484         me.viewEnd.gapPos += size;
485 
486 /*
487     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > it2 = begin(*me.data_container) + me.current.gapPos;
488     if (me.current != it2.current || me.prevAnchor != it2.prevAnchor || me.nextAnchor != it2.nextAnchor || me.anchorIdx != it2.anchorIdx)
489         std::cout<<"*";
490 */
491 }
492 
493 // ----------------------------------------------------------------------------
494 // Function removeGaps()
495 // ----------------------------------------------------------------------------
496 
497 template <typename TGaps, typename TGapAnchors, typename TCount>
498 inline typename Size<TGaps>::Type
removeGaps(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & it,TCount size_)499 removeGaps(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & it,
500            TCount size_)
501 {
502     TGapAnchors & anchors = _dataAnchors(*it.data_container);
503     typedef typename Iterator<TGapAnchors, Standard>::Type TAnchorsIter;
504 
505     typedef Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > TIter;
506     typedef typename TIter::TGapAnchor TGapAnchor;
507     // typedef typename Value<TGapAnchors>::Type   TGapAnchor;
508     typedef typename Position<TGapAnchor>::Type TPos;
509 
510     if (size_ <= 0 || !isGap(it))
511         return 0;
512     TPos size = size_;
513 
514     // static_cast<TGapAnchor>(Nothing());
515     // static_cast<TPos>(Nothing());
516     if (it.current.gapPos + size > it.nextAnchor.gapPos)
517         size = it.nextAnchor.gapPos - it.current.gapPos;
518 
519     if (it.prevAnchor.gapPos + it.current.seqPos == it.current.gapPos + it.prevAnchor.seqPos &&
520         it.current.gapPos + size == it.nextAnchor.gapPos)
521     {
522         // remove the gap
523         if (it.anchorIdx < (int)length(anchors))
524             erase(anchors, it.anchorIdx);
525         _getAnchor(it.nextAnchor, *it.data_container, it.anchorIdx);
526     }
527 
528     // shift anchors
529     if (it.anchorIdx < (int)length(anchors))
530     {
531         if (it.anchorIdx >= 0)
532         {
533             it.nextAnchor.gapPos -= size;
534             TAnchorsIter itA = begin(anchors, Standard());
535             TAnchorsIter itAEnd = end(anchors, Standard());
536             if (it.anchorIdx >= 0)
537                 itA += it.anchorIdx;
538             for (; itA != itAEnd; ++itA)
539                 (*itA).gapPos -= size;
540         }
541         else
542             // remove gap before the sequence and in (or at the left boundary of) the view
543             if (it.current.gapPos >= it.viewBegin.gapPos)
544             {
545                 // assure that we don't remove more gaps than available
546                 if (size > it.nextAnchor.gapPos - it.current.gapPos)
547                     size = it.nextAnchor.gapPos - it.current.gapPos;
548                 container(it).data_cutBegin += size;
549                 it.viewBegin.gapPos += size;
550                 it.current.gapPos += size;
551                 return size;
552             }
553     }
554     else
555     {
556         if (it.current.gapPos <= it.viewEnd.gapPos)
557             container(it).data_cutEnd += size;
558     }
559     if (it.current.gapPos <= it.viewEnd.gapPos)
560         it.viewEnd.gapPos -= size;
561 
562     return size;
563 /*
564     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > it2 = begin(*me.data_container) + me.current.gapPos;
565     if (me.current != it2.current || me.prevAnchor != it2.prevAnchor || me.nextAnchor != it2.nextAnchor || me.anchorIdx != it2.anchorIdx)
566         std::cout<<"*";
567 */
568 }
569 
570 // ----------------------------------------------------------------------------
571 // Helper Function _goNextGapAnchorIterator()
572 // ----------------------------------------------------------------------------
573 
574 template <typename T>
575 inline void
_goNextGapAnchorIterator(T & me)576 _goNextGapAnchorIterator(T & me)
577 {
578     if (me.current.gapPos < me.nextAnchor.gapPos)
579     {
580         ++me.current.gapPos;
581         if (me.current.seqPos < me.nextAnchor.seqPos)
582             ++me.current.seqPos;
583     }
584     while (me.current.gapPos == me.nextAnchor.gapPos)
585     {
586         me.current = me.prevAnchor = me.nextAnchor;
587         _getAnchor(me.nextAnchor, *me.data_container, ++me.anchorIdx + 1);
588     }
589 }
590 
591 // ----------------------------------------------------------------------------
592 // Helper Function _goPreviousGapAnchorIterator()
593 // ----------------------------------------------------------------------------
594 
595 template <typename T>
596 inline void
_goPreviousGapAnchorIterator(T & me)597 _goPreviousGapAnchorIterator(T & me)
598 {
599     while (me.current.gapPos == me.prevAnchor.gapPos)
600     {
601         me.current = me.nextAnchor = me.prevAnchor;
602         _getAnchor(me.prevAnchor, *me.data_container, --me.anchorIdx);
603     }
604     --me.current.gapPos;
605     if (me.nextAnchor.seqPos - me.prevAnchor.seqPos > me.current.gapPos - me.prevAnchor.gapPos)
606         me.current.seqPos = me.prevAnchor.seqPos + (me.current.gapPos - me.prevAnchor.gapPos);
607     else
608         me.current.seqPos = me.nextAnchor.seqPos;
609 }
610 
611 // ----------------------------------------------------------------------------
612 // Helper Function _goToGapAnchorIterator()
613 // ----------------------------------------------------------------------------
614 
615 template <typename T, typename TPos>
616 inline void
_goToGapAnchorIterator(T & me,TPos pos)617 _goToGapAnchorIterator(T & me, TPos pos)
618 {
619     typedef typename T::TGapAnchors                 TGapAnchors;
620     typedef typename Value<TGapAnchors>::Type       TGapAnchor;
621     typedef typename Position<TGapAnchor>::Type     TAnchorPos;
622     typedef typename MakeSigned<TAnchorPos>::Type   TAnchorSPos;
623 
624     if (isNegative(pos))
625         me.anchorIdx = -1;
626     else
627     {
628         TGapAnchors const & anchors = _dataAnchors(*me.data_container);
629         if (!empty(anchors))
630         {
631             me.anchorIdx = upperBoundGapAnchor(anchors, pos, SortGapPos()) - begin(anchors, Standard());
632             if (me.anchorIdx < (int)length(anchors))
633                 if (anchors[me.anchorIdx].gapPos == (TAnchorPos)pos && anchors[me.anchorIdx].seqPos != (TAnchorPos)me.seqLength)
634                     ++me.anchorIdx;
635         }
636         else
637         {
638             me.anchorIdx = ((TAnchorSPos)pos < me.seqLength)? 0: 1;
639         }
640     }
641     _getAnchor(me.prevAnchor, *me.data_container, me.anchorIdx);
642     _getAnchor(me.nextAnchor, *me.data_container, me.anchorIdx + 1);
643 
644     me.current.gapPos = pos;
645     if (me.nextAnchor.seqPos - me.prevAnchor.seqPos > (int)pos - me.prevAnchor.gapPos)
646         me.current.seqPos = me.prevAnchor.seqPos + ((int)pos - me.prevAnchor.gapPos);
647     else
648         me.current.seqPos = me.nextAnchor.seqPos;
649 }
650 
651 // ----------------------------------------------------------------------------
652 // Function goNext()
653 // ----------------------------------------------------------------------------
654 
655 template <typename TGaps, typename TGapAnchors>
656 inline void
goNext(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> & me)657 goNext(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > & me)
658 {
659     _goNextGapAnchorIterator(me);
660 }
661 
662 // ----------------------------------------------------------------------------
663 // Function goPrevious()
664 // ----------------------------------------------------------------------------
665 
666 template <typename TGaps, typename TGapAnchors>
667 inline void
goPrevious(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> & me)668 goPrevious(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > & me)
669 {
670     _goPreviousGapAnchorIterator(me);
671 }
672 
673 // ----------------------------------------------------------------------------
674 // Function goFurther()
675 // ----------------------------------------------------------------------------
676 
677 template <typename TGaps, typename TGapAnchors, typename TSize>
678 inline void
goFurther(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> & me,TSize steps)679 goFurther(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > & me, TSize steps)
680 {
681     _goToGapAnchorIterator(me, me.current.gapPos + steps);
682 }
683 
684 // ----------------------------------------------------------------------------
685 // Function position()
686 // ----------------------------------------------------------------------------
687 
688 // Returns clipped view position of gaps iterator.
689 
690 template <typename TGaps, typename TGapAnchors>
691 inline typename Position<Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > >::Type
position(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & it)692 position(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & it)
693 {
694     return it.current.gapPos - it.viewBegin.gapPos;
695 }
696 
697 // ----------------------------------------------------------------------------
698 // Function difference()
699 // ----------------------------------------------------------------------------
700 
701 template <typename TGaps, typename TGapAnchors>
702 inline typename Difference<Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > >::Type
difference(Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & lhs,Iter<TGaps,GapsIterator<AnchorGaps<TGapAnchors>>> const & rhs)703 difference(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & lhs,
704            Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & rhs)
705 {
706     return lhs.current.gapPos - rhs.current.gapPos;
707 }
708 
709 // ----------------------------------------------------------------------------
710 // Function operator-()                                            [difference]
711 // ----------------------------------------------------------------------------
712 
713 template <typename TGaps, typename TGapAnchors>
714 inline typename Difference<Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > >::Type
715 operator-(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & lhs,
716           Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & rhs)
717 {
718     return difference(lhs, rhs);
719 }
720 
721 // ----------------------------------------------------------------------------
722 // Function operator-()                                         [copy movement]
723 // ----------------------------------------------------------------------------
724 
725 template <typename TGaps, typename TGapAnchors, typename TDifference>
726 inline Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > >
727 operator-(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & lhs, TDifference d)
728 {
729     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > result = lhs;
730     result -= d;
731     return result;
732 }
733 
734 // ----------------------------------------------------------------------------
735 // Function operator+()                                         [copy movement]
736 // ----------------------------------------------------------------------------
737 
738 template <typename TGaps, typename TGapAnchors, typename TDifference>
739 inline Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > >
740 operator+(Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > const & lhs, TDifference d)
741 {
742     Iter<TGaps, GapsIterator<AnchorGaps<TGapAnchors> > > result = lhs;
743     result += d;
744     return result;
745 }
746 
747 
748 }  // namespace seqan
749 
750 #endif  // #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_GAPS_ITERATOR_ANCHOR_H_
751