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: Anne-Katrin Emde <emde@fu-berlin.de>
33 // ==========================================================================
34
35
36 // TODO(holtgrew): Bring back random center computation using new random module?
37
38 // (weese:) This header needs improvement
39 //
40 // documentation lacks:
41 // - is createIntervalTree clearing the tree first (what happens after a second call of createIntervalTree?)
42 // - where is clear(itree)?
43 // - some flags are missing in the dddoc entries, e.g. ..cat:, class:, or ..type: of arguments
44 //
45 // interface flaws:
46 // - there are 2 interfaces for createIntervalTree
47 // - based on IntervalTree
48 // - based on Graph and PropertyMap (should be removed)
49 // - addInterval
50
51
52 #ifndef SEQAN_HEADER_MISC_INTERVAL_TREE_H
53 #define SEQAN_HEADER_MISC_INTERVAL_TREE_H
54
55 #include <seqan/graph_types.h>
56
57
58 namespace SEQAN_NAMESPACE_MAIN {
59 //////////////////////////////////////////////////////////////////////////////
60 // Graph - Interval Tree Types
61 //////////////////////////////////////////////////////////////////////////////
62
63 ///---------------------------------------------------------------///
64
65 //////////////////// Interval and ID type ///////////////////
66
67 // TODO(holtgrew): Are these actual the first/last or begin/end entries?
68
69 /*!
70 * @class IntervalAndCargo
71 * @headerfile <seqan/misc/interval_tree.h>
72 * @brief A simple record type that stores an interval and a cargo value.
73 *
74 * @signature template <[typename TValue[, typename TCargo]]>
75 * class IntervalAndCargo;
76 *
77 * @tparam TValue The value type. Default: <tt>int</tt>.
78 * @tparam TCargo The cargo type. Default: <tt>int</tt>.
79 *
80 * @fn IntervalAndCargo::IntervalAndCargo
81 * @brief Constructor
82 *
83 * @signature IntervalAndCargo::IntervalAndCargo();
84 * @signature IntervalAndCargo::IntervalAndCargo(i1, i2, cargo);
85 *
86 * @param[in] i1 The first element in the interval.
87 * @param[in] i2 The last element in the interval.
88 * @param[in] cargo The cargo to store together with the interval.
89 *
90 * @var TValue IntervalAndCargo::i1;
91 * @brief The first element in the interval.
92 *
93 * @var TValue IntervalAndCargo::i2;
94 * @brief The last element in the interval.
95 *
96 * @var TCargo IntervalAndCargo::cargo;
97 * @brief The stored cargo.
98 */
99
100 template <typename TValue = int, typename TCargo = int>
101 class IntervalAndCargo
102 {
103 public:
104 TValue i1;
105
106 TValue i2;
107
108 TCargo cargo;
109
IntervalAndCargo()110 IntervalAndCargo() :
111 i1(), i2(), cargo()
112 {}
113
IntervalAndCargo(TValue i1,TValue i2,TCargo cargo)114 IntervalAndCargo(TValue i1, TValue i2, TCargo cargo) :
115 i1(i1), i2(i2), cargo(cargo)
116 {
117 SEQAN_CHECKPOINT
118 }
119
120 };
121
122
123
124 /////////////////////// Point and ID type ////////////////
125
126 /*!
127 * @class PointAndCargo
128 * @headerfile <seqan/misc/interval_tree.h>
129 * @brief Simple record class storing a point (one-value) interval and cargo.
130 *
131 * @signature template <[typename TValue[, typename TCargo]]>
132 * class PointAndCargo;
133 *
134 * @tparam TValue The value type.
135 * @tparam TCargo The cargo type.
136 *
137 * @fn PointAndCargo::PointAndCargo
138 * @brief Constructor
139 *
140 * @signature PointAndCargo::PointAndCargo();
141 * @signature PointAndCargo::PointAndCargo(point, cargo);
142 *
143 * @param[in] point The point to store.
144 * @param[in] cargo The cargo to store.
145 *
146 * @var TValue PointAndCargo::point
147 * @brief The point to store.
148 *
149 * @var TCargo PointAndCargo::cargo
150 * @brief The cargo to store.
151 */
152
153 template <typename TValue = int, typename TCargo = int>
154 class PointAndCargo
155 {
156 public:
157 TValue point;
158
159 TCargo cargo;
160
PointAndCargo()161 PointAndCargo() :
162 point(), cargo()
163 {}
164
PointAndCargo(TValue point,TCargo cargo)165 PointAndCargo(TValue point, TCargo cargo) :
166 point(point), cargo(cargo)
167 {}
168 };
169
170 ///////////////////////////////////////////////////////////////////////////
171 /////////////////////////// IntervalTreeNode ///////////////////////////
172
173 /*!
174 * @defgroup IntervalTreeNodeTypeTags IntervalTree Node Types Tags
175 * @brief Tags to select the node type for @link IntervalTree @endlink.
176 *
177 * @tag IntervalTreeNodeTypeTags#StorePointsOnly
178 * @headerfile <seqan/misc/interval_tree.h>
179 * @signature struct StorePointsOnly {};
180 * @brief The tree nodes store points.
181 *
182 * @tag IntervalTreeNodeTypeTags#StoreIntervals
183 * @headerfile <seqan/misc/interval_tree.h>
184 * @signature struct StoreIntervals {};
185 * @brief The tree nodes store intervals.
186 */
187
188 struct StorePointsOnly {};
189
190
191 struct StoreIntervals {};
192
193 /*!
194 * @class IntervalTreeNode
195 * @headerfile <seqan/misc/interval_tree.h>
196 * @brief Element of @link IntervalTree @endlink.
197 *
198 * @signature template <typename TInterval[, typename TSpec]>
199 * class IntervalTreeNode;
200 *
201 * @tparam TInterval The type to use for intervals.
202 * @tparam TSpec The specializing tag. Default: @link IntervalTreeNodeTypeTags#StorePointsOnly @endlink.
203 */
204
205
206 template <typename TInterval, typename TSpec = StorePointsOnly>
207 class IntervalTreeNode;
208
209 /*!
210 * @class StoreIntervalsIntervalTreeNode
211 * @extends IntervalTreeNode
212 * @headerfile <seqan/misc/interval_tree.h>
213 * @brief An IntervalTreeNode that stores intervals explicitely in each node.
214 *
215 * @signature template <typename TInterval>
216 * class IntervalTreeNode<TInterval, StoreIntervals>;
217 *
218 * @tparam TInterval The Interval type to use.
219 */
220
221 /*!
222 * @var TValue StoreIntervalsIntervalTreeNode::center;
223 * @brief Center of the interval tree node.
224 *
225 * @var TString StoreIntervalsIntervalTreeNode::list1;
226 * @brief @link AllocString @endlink of intervals sorted by begin point.
227 *
228 * @var TString StoreIntervalsIntervalTreeNode::list;
229 * @brief @link AllocString @endlink of intervals sorted by end point.
230 */
231
232 template <typename TInterval>
233 class IntervalTreeNode<TInterval, StoreIntervals>
234 {
235 public:
236 typedef typename Value<TInterval>::Type TValue;
237
238 TValue center;
239 String<TInterval> list1;
240 String<TInterval> list2;
241
IntervalTreeNode()242 IntervalTreeNode() :
243 center()
244 {}
245 };
246
247 /*!
248 * @class StorePointsOnlyIntervalTreeNode
249 * @extends IntervalTreeNode
250 * @headerfile <seqan/misc/interval_tree.h>
251 * @brief An IntervalTreeNode that stores only the relevant points in each node.
252 *
253 * Only the end points of the intervals in the list sorted by endpoints (list2) and only the begin point of the interval
254 * list sorted by begin points (list1) are stored.
255 *
256 * @signature template <typename TInterval>
257 * class IntervalTreeNode<TInterval, StoreIntervals>;
258 *
259 * @tparam TInterval The Interval type to use.
260 */
261
262 /*!
263 * @var TValue StorePointsOnlyIntervalTreeNode::center;
264 * @brief Center of the interval.
265 *
266 * @var TString StorePointsOnlyIntervalTreeNode::list1;
267 * @brief Points with cargo sorted by the begin points.
268 *
269 * @var TString StorePointsOnlyIntervalTreeNode::list2;
270 * @brief Points with cargo sorted by the end points.
271 */
272
273 template <typename TInterval>
274 class IntervalTreeNode<TInterval, StorePointsOnly>
275 {
276 public:
277 typedef typename Cargo<TInterval>::Type TCargo;
278 typedef typename Value<TInterval>::Type TValue;
279
280 TValue center;
281 String<PointAndCargo<TValue, TCargo> > list1;
282 String<PointAndCargo<TValue, TCargo> > list2;
283
IntervalTreeNode()284 IntervalTreeNode() :
285 center()
286 {
287 SEQAN_CHECKPOINT
288 }
289
IntervalTreeNode(IntervalTreeNode const & other)290 IntervalTreeNode(IntervalTreeNode const & other) :
291 center(other.center),
292 list1(other.list1),
293 list2(other.list2)
294 {
295 SEQAN_CHECKPOINT
296 }
297
298 };
299
300
301
302
303
304 //////////////////////////////////////////////////////////////////////////////
305 // Graph - Interval Tree
306 //////////////////////////////////////////////////////////////////////////////
307
308 /*!
309 * @class IntervalTree
310 * @headerfile <seqan/misc/interval_tree.h>
311 * @brief Data structure for efficient interval storage.
312 *
313 * @signature template <[typename TValue[, typename TCargo]]>
314 * class IntervalTree;
315 *
316 * @tparam TValue The type to use for coordinates. Default: <tt>int</tt>.
317 * @tparam TCargo The type to use for cargo. Default: <tt>unsigned</tt>.
318 *
319 * @section Remarks
320 *
321 * If the intervals are not associated with cargos/IDs, they will be numbered consecutively.
322 *
323 * @section Example
324 *
325 * The following example creates an integer interval tree with string keys. This tree is quired for keys of intervals
326 * that overlap the interval <tt>[550, 990)</tt>.
327 *
328 * @include demos/misc/interval_tree_example.cpp
329 *
330 * The resulting keys are:
331 *
332 * @code{.console}
333 * gene
334 * exon2
335 * coding2
336 * @endcode
337 *
338 *
339 * @fn IntervalTree::IntervalTree
340 * @brief Constructor
341 *
342 * @signature IntervalTree::IntervalTree();
343 * @signature IntervalTree::IntervalTree(intervals);
344 * @signature IntervalTree::IntervalTree(intervals[, center]);
345 * @signature IntervalTree::IntervalTree(intervals[, tag]);
346 * @signature IntervalTree::IntervalTree(intervalBegins, intervalEnds, [intervalCargos,] len);
347 *
348 * @param[in] intervals Container of intervals. A strin gof <tt>IntervalAndCargo<Value, TCargo></tt>
349 * objects, see @link IntervalAndCargo @endlink.
350 * @param[in] intervalBegins
351 * Iterator pointing to begin position of first interval.
352 * @param[in] intervalEnds
353 * Iterator pointing to end position of first interval.
354 * @param[in] intervalCargos
355 * Iterator pointing to cargo/ids for intervals.
356 * @param[in] len Number of intervals to store in tree.
357 * @param[in] tag Tag for tree construction method.
358 */
359
360 template <typename TValue = int, typename TCargo = unsigned int>
361 class IntervalTree
362 {
363 public:
364 typedef Graph<Directed<void, WithoutEdgeId> > TGraph;
365 typedef IntervalAndCargo<TValue, TCargo> TInterval;
366 typedef IntervalTreeNode<TInterval> TNode;
367 typedef String<TNode> TPropertyMap;
368
369 TGraph g;
370 TPropertyMap pm;
371 size_t interval_counter;
372
IntervalTree()373 IntervalTree()
374 {
375 SEQAN_CHECKPOINT
376 interval_counter = 0;
377 }
378
379 template <typename TIterator, typename TCargoIterator>
IntervalTree(TIterator interval_begins,TIterator interval_ends,TCargoIterator interval_cargos,size_t len)380 IntervalTree(TIterator interval_begins,
381 TIterator interval_ends,
382 TCargoIterator interval_cargos,
383 size_t len)
384 {
385 SEQAN_CHECKPOINT
386 String<TInterval> intervals;
387 resize(intervals, len);
388 size_t i = 0;
389 while (i < len)
390 {
391 intervals[i].i1 = value(interval_begins);
392 ++interval_begins;
393 intervals[i].i2 = value(interval_ends);
394 ++interval_ends;
395 intervals[i].cargo = value(interval_cargos);
396 ++interval_cargos;
397 ++i;
398 }
399 createIntervalTree(*this, intervals);
400 }
401
402 template <typename TIterator>
IntervalTree(TIterator interval_begins,TIterator interval_ends,size_t len)403 IntervalTree(TIterator interval_begins,
404 TIterator interval_ends,
405 size_t len)
406 {
407 SEQAN_CHECKPOINT
408 String<TInterval> intervals;
409 resize(intervals, len);
410 size_t i = 0;
411 while (i < len)
412 {
413 intervals[i].i1 = value(interval_begins);
414 ++interval_begins;
415 intervals[i].i2 = value(interval_ends);
416 ++interval_ends;
417 intervals[i].cargo = i;
418 ++i;
419 }
420 createIntervalTree(*this, intervals);
421 }
422
IntervalTree(String<TInterval> intervals)423 IntervalTree(String<TInterval> intervals)
424 {
425 SEQAN_CHECKPOINT
426 createIntervalTree(*this, intervals);
427 }
428
429 template <typename TTagSpec>
IntervalTree(String<TInterval> intervals,Tag<TTagSpec> const tag)430 IntervalTree(String<TInterval> intervals, Tag<TTagSpec> const tag)
431 {
432 SEQAN_CHECKPOINT
433 interval_counter = length(intervals);
434 createIntervalTree(g, pm, intervals, tag);
435 }
436
IntervalTree(String<TInterval> intervals,TValue center)437 IntervalTree(String<TInterval> intervals, TValue center)
438 {
439 SEQAN_CHECKPOINT
440 interval_counter = length(intervals);
441 createIntervalTree(g, pm, intervals, center);
442 }
443
444 };
445
446
447
448 ///////Specs for the way interval centers are determined
449
450 //template <typename TSpec = SpecPointAndCargo>
451 struct TagComputeCenter_;
452 typedef Tag<TagComputeCenter_> const ComputeCenter;
453
454
455 ///////////////////////////////////////////////////////////////////////////
456 ///////////////////// IntervalAndCargo functions //////////////////////////
457 ///////////////////////////////////////////////////////////////////////////
458
459 /*!
460 * @fn IntervalAndCargo#leftBoundary
461 * @brief Access to left boundary.
462 *
463 * @signature TBoundary leftBoundary(interval);
464 *
465 * @param[in] interval The IntervalAndCargo to query for its left boundary.
466 *
467 * @return TBoundary Reference to the left boundary value.
468 */
469
470 template <typename TValue, typename TCargo>
471 TValue &
leftBoundary(IntervalAndCargo<TValue,TCargo> & interval)472 leftBoundary(IntervalAndCargo<TValue, TCargo> & interval)
473 {
474 return interval.i1;
475 }
476
477 template <typename TValue, typename TCargo>
478 TValue const &
leftBoundary(IntervalAndCargo<TValue,TCargo> const & interval)479 leftBoundary(IntervalAndCargo<TValue, TCargo> const & interval)
480 {
481 return interval.i1;
482 }
483
484 /*!
485 * @fn IntervalAndCargo#rightBoundary
486 * @brief Access to right boundary.
487 *
488 * @signature TBoundary rightBoundary(interval);
489 *
490 * @param[in] interval The IntervalAndCargo to query for its right boundary.
491 *
492 * @return TBoundary Reference to the right boundary value.
493 */
494
495 template <typename TValue, typename TCargo>
496 TValue &
rightBoundary(IntervalAndCargo<TValue,TCargo> & interval)497 rightBoundary(IntervalAndCargo<TValue, TCargo> & interval)
498 {
499 return interval.i2;
500 }
501
502 template <typename TValue, typename TCargo>
503 TValue const &
rightBoundary(IntervalAndCargo<TValue,TCargo> const & interval)504 rightBoundary(IntervalAndCargo<TValue, TCargo> const & interval)
505 {
506 return interval.i2;
507 }
508
509 /*!
510 * @fn IntervalAndCargo#getLeftBoundary
511 * @brief Access to getLeft boundary.
512 *
513 * @signature TBoundary getLeftBoundary(interval);
514 *
515 * @param[in] interval The IntervalAndCargo to query for its left boundary.
516 *
517 * @return TBoundary Copy of the left boundary value.
518 */
519
520 template <typename TValue, typename TCargo>
521 TValue
getLeftBoundary(IntervalAndCargo<TValue,TCargo> const & interval)522 getLeftBoundary(IntervalAndCargo<TValue, TCargo> const & interval)
523 {
524 return interval.i1;
525 }
526
527 /*!
528 * @fn IntervalAndCargo#getRightBoundary
529 * @brief Access to getRight boundary.
530 *
531 * @signature TBoundary getRightBoundary(interval);
532 *
533 * @param[in] interval The IntervalAndCargo to query for its right boundary.
534 *
535 * @return TBoundary Copy of the right boundary value.
536 */
537
538 template <typename TValue, typename TCargo>
539 TValue
getRightBoundary(IntervalAndCargo<TValue,TCargo> const & interval)540 getRightBoundary(IntervalAndCargo<TValue, TCargo> const & interval)
541 {
542 return interval.i2;
543 }
544
545 /*!
546 * @fn IntervalAndCargo#cargo
547 * @brief Access to the cargo.
548 *
549 * @signature TCargo cargo(interval);
550 *
551 * @param[in] interval The IntervalAndCargo to query for its cargo.
552 *
553 * @return TCargo Reference to the cargo member.
554 */
555
556 template <typename TValue, typename TCargo>
557 TCargo const &
cargo(IntervalAndCargo<TValue,TCargo> const & interval)558 cargo(IntervalAndCargo<TValue, TCargo> const & interval)
559 {
560 return interval.cargo;
561 }
562
563 template <typename TValue, typename TCargo>
564 TCargo &
cargo(IntervalAndCargo<TValue,TCargo> & interval)565 cargo(IntervalAndCargo<TValue, TCargo> & interval)
566 {
567 return interval.cargo;
568 }
569
570 /*!
571 * @fn IntervalAndCargo#getCargo
572 * @brief Access to the cargo.
573 *
574 * @signature TCargo getCargo(interval);
575 *
576 * @param[in] interval The IntervalAndCargo to query for its cargo.
577 *
578 * @return TCargo Copy of the cargo member.
579 */
580
581 template <typename TValue, typename TCargo>
582 TCargo
getCargo(IntervalAndCargo<TValue,TCargo> const & interval)583 getCargo(IntervalAndCargo<TValue, TCargo> const & interval)
584 {
585 return interval.cargo;
586 }
587
588 /////////////////// Metafunctions //////////////////////
589
590 /*!
591 * @mfn IntervalAndCargo#Value
592 * @brief Return the value type.
593 *
594 * @signature Value<TIntervalAndCargo>::Type;
595 */
596
597 template <typename TValue, typename TCargo>
598 struct Value<IntervalAndCargo<TValue, TCargo> >
599 {
600 typedef TValue Type;
601 };
602
603 /*!
604 * @mfn IntervalAndCargo#Cargo
605 * @brief Return the cargo type.
606 *
607 * @signature Cargo<TIntervalAndCargo>::Type;
608 */
609
610 template <typename TValue, typename TCargo>
611 struct Cargo<IntervalAndCargo<TValue, TCargo> >
612 {
613 typedef TCargo Type;
614 };
615
616
617 ///////////////////////////////////////////////////////////////////////////
618 ///////////////////// PointAndCargo functions /////////////////////////////
619 ///////////////////////////////////////////////////////////////////////////
620
621 /*!
622 * @fn PointAndCargo#leftBoundary
623 * @brief Access to left boundary.
624 *
625 * @signature TBoundary leftBoundary(point);
626 *
627 * @param[in] point The PointAndCargo to query for its left boundary.
628 *
629 * @return TBoundary Reference to the left boundary value.
630 */
631
632 template <typename TValue, typename TCargo>
633 TValue const &
634 leftBoundary(PointAndCargo<TValue, TCargo> const & point)
635 {
636 return point.point;
637 }
638
639 template <typename TValue, typename TCargo>
640 TValue &
641 leftBoundary(PointAndCargo<TValue, TCargo> & point)
642 {
643 return point.point;
644 }
645
646 /*!
647 * @fn PointAndCargo#rightBoundary
648 * @brief Access to right boundary.
649 *
650 * @signature TBoundary rightBoundary(point);
651 *
652 * @param[in] point The PointAndCargo to query for its right boundary.
653 *
654 * @return TBoundary Reference to the right boundary value.
655 */
656
657 template <typename TValue, typename TCargo>
658 TValue const &
659 rightBoundary(PointAndCargo<TValue, TCargo> const & point)
660 {
661 return point.point;
662 }
663
664 template <typename TValue, typename TCargo>
665 TValue &
666 rightBoundary(PointAndCargo<TValue, TCargo> & point)
667 {
668 return point.point;
669 }
670
671 template <typename TValue, typename TCargo>
672 TValue
673 getLeftBoundary(PointAndCargo<TValue, TCargo> const & point)
674 {
675 return point.point;
676 }
677
678 /*!
679 * @fn PointAndCargo#getLeftBoundary
680 * @brief Access to getLeft boundary.
681 *
682 * @signature TBoundary getLeftBoundary(point);
683 *
684 * @param[in] point The PointAndCargo to query for its left boundary.
685 *
686 * @return TBoundary Copy of the left boundary value.
687 */
688
689 template <typename TValue, typename TCargo>
690 TValue
691 getRightBoundary(PointAndCargo<TValue, TCargo> const & point)
692 {
693 return point.point;
694 }
695
696 /*!
697 * @fn PointAndCargo#cargo
698 * @brief Access to the cargo.
699 *
700 * @signature TCargo cargo(point);
701 *
702 * @param[in] point The PointAndCargo to query for its cargo.
703 *
704 * @return TCargo Reference to the cargo member.
705 */
706
707
708 template <typename TValue, typename TCargo>
709 TCargo const &
710 cargo(PointAndCargo<TValue, TCargo> const & point)
711 {
712 return point.cargo;
713 }
714
715 template <typename TValue, typename TCargo>
716 TCargo &
717 cargo(PointAndCargo<TValue, TCargo> & point)
718 {
719 return point.cargo;
720 }
721
722 /*!
723 * @fn PointAndCargo#getCargo
724 * @brief Access to the cargo.
725 *
726 * @signature TCargo getCargo(point);
727 *
728 * @param[in] point The PointAndCargo to query for its cargo.
729 *
730 * @return TCargo Copy of the cargo member.
731 */
732
733 template <typename TValue, typename TCargo>
734 TCargo
735 getCargo(PointAndCargo<TValue, TCargo> const & point)
736 {
737 return point.cargo;
738 }
739
740 ////////////////// Metafunctions //////////////////
741
742 /*!
743 * @mfn PointAndCargo#Value
744 * @brief Return the value type.
745 *
746 * @signature Value<TPointAndCargo>::Type;
747 */
748
749 template <typename TValue, typename TCargo>
750 struct Value<PointAndCargo<TValue, TCargo> >
751 {
752 typedef TValue Type;
753 };
754
755 /*!
756 * @mfn PointAndCargo#Cargo
757 * @brief Return the cargo type.
758 *
759 * @signature Cargo<TPointAndCargo>::Type;
760 */
761
762 template <typename TValue, typename TCargo>
763 struct Cargo<PointAndCargo<TValue, TCargo> >
764 {
765 typedef TCargo Type;
766 };
767
768
769 //// Comparators
770 template <typename TPair>
771 bool _less_compI1_ITree(TPair const & p1, TPair const & p2)
772 {
773 return leftBoundary(p1) < leftBoundary(p2);
774 }
775
776 template <typename TPair>
777 bool _greater_compI2_ITree(TPair const & p1, TPair const & p2)
778 {
779 return rightBoundary(p1) > rightBoundary(p2);
780 }
781
782 ///////////////////////////////////////////////////////////////////////////
783 ///////////////////// IntervalTreeNode functions //////////////////////////
784 ///////////////////////////////////////////////////////////////////////////
785
786
787
788
789
790 // internal set node functions
791 template <typename TValue, typename TInterval>
792 void
793 _setIntervalTreeNode(IntervalTreeNode<TInterval, StoreIntervals> & knot, TValue center, TInterval const & interval)
794 {
795 knot.center = center;
796 appendValue(knot.list1, interval);
797 appendValue(knot.list2, interval);
798 }
799
800 // append intervals to lists in node knot
801 template <typename TInterval>
802 void
803 _appendIntervalTreeNodeLists(IntervalTreeNode<TInterval, StoreIntervals> & knot, TInterval const & interval)
804 {
805 appendValue(knot.list1, interval);
806 appendValue(knot.list2, interval);
807 }
808
809 //internal set node functions
810 template <typename TValue, typename TInterval>
811 void
812 _setIntervalTreeNode(IntervalTreeNode<TInterval, StorePointsOnly> & knot, TValue center, TInterval const & interval)
813 {
814 knot.center = center;
815 appendValue(knot.list1, PointAndCargo<TValue, typename Cargo<TInterval>::Type>(leftBoundary(interval), cargo(interval)));
816 appendValue(knot.list2, PointAndCargo<TValue, typename Cargo<TInterval>::Type>(rightBoundary(interval), cargo(interval)));
817 }
818
819 template <typename TInterval>
820 void
821 _appendIntervalTreeNodeLists(IntervalTreeNode<TInterval, StorePointsOnly> & knot, TInterval const & interval)
822 {
823 appendValue(knot.list1, PointAndCargo<typename Value<TInterval>::Type, typename Cargo<TInterval>::Type>(leftBoundary(interval), cargo(interval)));
824 appendValue(knot.list2, PointAndCargo<typename Value<TInterval>::Type, typename Cargo<TInterval>::Type>(rightBoundary(interval), cargo(interval)));
825
826
827 }
828
829 /////////////////// Metafunctions ///////////////////////
830
831 /*!
832 * @mfn IntervalTreeNode#Value
833 * @brief Return value type.
834 *
835 * @signature Value<TNode>::Type;
836 */
837
838 template <typename TInterval, typename TSpec>
839 struct Value<IntervalTreeNode<TInterval, TSpec> >
840 {
841 typedef typename Value<TInterval>::Type Type;
842 };
843
844 /*!
845 * @mfn IntervalTreeNode#Cargo
846 * @brief Return cargo type.
847 *
848 * @signature Cargo<TNode>::Type;
849 */
850
851 template <typename TInterval, typename TSpec>
852 struct Cargo<IntervalTreeNode<TInterval, TSpec> >
853 {
854 typedef typename Cargo<TInterval>::Type Type;
855 };
856
857
858 /*!
859 * @mfn IntervalTreeNode#ListType
860 * @brief Type of the lists in tree nodes.
861 *
862 * @signature ListType<T>::Type;
863 */
864
865 template <typename T>
866 struct ListType;
867
868
869 template <typename TInterval>
870 struct ListType<IntervalTreeNode<TInterval, StorePointsOnly> >
871 {
872 typedef String<PointAndCargo<typename Value<TInterval>::Type, typename Cargo<TInterval>::Type> > Type;
873
874 };
875
876
877 template <typename TInterval>
878 struct ListType<IntervalTreeNode<TInterval, StoreIntervals> >
879 {
880 typedef String<IntervalAndCargo<typename Value<TInterval>::Type, typename Cargo<TInterval>::Type> > Type;
881
882 };
883
884
885
886 ///////////////////////////////////////////////////////////////////////////
887 /////////////////////// IntervalTree functions ////////////////////////////
888 ///////////////////////////////////////////////////////////////////////////
889
890 /*!
891 * @fn IntervalTree#createIntervalTree
892 * @brief Create an interval tree.
893 *
894 * @signature void createIntervalTree(intervalTree, intervals[, tag]);
895 * @signature void createIntervalTree(g, pm, intervals[, tag]);
896 * @signature void createIntervalTree(g, pm, intervals, center[, tag]]);
897 *
898 * @param[in,out] intervalTree An interval tree Types: IntervalTree
899 * @param[in,out] g DirectedGraph to create interval tree in. Types: @link Graph @endlink.
900 * @param[in,out] pm Property map to use for the created interval tree.
901 * @param[in] tag Tag for tree construction method;
902 * @param[in] intervals Container of intervals. A string of <tt>IntervalAndCargo<TValue, TCargo></tt>
903 * objects, see @link IntervalAndCargo @endlink. Types: @link AllocString @endlink.
904 */
905
906 template <typename TGraph, typename TPropertyMap, typename TIntervals, typename TSpec>
907 inline void
908 createIntervalTree(TGraph & g,
909 TPropertyMap & pm,
910 TIntervals & intervals,
911 Tag<TSpec> const tag)
912 {
913 typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
914 typedef typename Value<TIntervals>::Type TInterval;
915 typedef typename Value<TInterval>::Type TValue;
916
917 reserve(g.data_vertex, length(intervals));
918 reserve(pm, length(intervals));
919
920 TVertexDescriptor root = addVertex(g);
921 resizeVertexMap(pm, g);
922
923 if (length(intervals) > 0u)
924 {
925 TValue center = _calcIntervalTreeRootCenter(intervals);
926
927 std::sort(begin(intervals, Standard()), end(intervals, Standard()), _less_compI1_ITree<TInterval>);
928
929 String<TInterval *> interval_pointers;
930 // interval tree stores pointers to intervals, not original intervals
931 _makePointerInterval(intervals, interval_pointers);
932
933 _createIntervalTree(g, pm, interval_pointers, root, (TValue)0.0, center, length(intervals), tag);
934 reserve(pm, length(pm), Exact());
935 reserve(g.data_vertex, length(g.data_vertex), Exact());
936 }
937 }
938
939 template <typename TGraph, typename TPropertyMap, typename TIntervals>
940 inline void
941 createIntervalTree(TGraph & g,
942 TPropertyMap & pm,
943 TIntervals & intervals)
944 {
945 createIntervalTree(g, pm, intervals, ComputeCenter());
946 }
947
948 template <typename TGraph, typename TPropertyMap, typename TIntervals, typename TSpec>
949 inline void
950 createIntervalTree(TGraph & g,
951 TPropertyMap & pm,
952 TIntervals & intervals,
953 typename Value<typename Value<TIntervals>::Type>::Type center,
954 Tag<TSpec> const tag)
955 {
956 typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
957 typedef typename Value<TIntervals>::Type TInterval;
958 typedef typename Value<typename Value<TIntervals>::Type>::Type TValue;
959
960 reserve(g.data_vertex, length(intervals));
961 reserve(pm, length(intervals));
962
963 TVertexDescriptor root = addVertex(g);
964 resizeVertexMap(pm, g);
965
966 if (length(intervals) > 0u)
967 {
968 TInterval a;
969 typename Iterator<TIntervals, Standard>::Type begin_ = begin(intervals, Standard());
970 typename Iterator<TIntervals, Standard>::Type end_ = end(intervals, Standard());
971 std::sort(begin_, end_, _less_compI1_ITree<TInterval>);
972
973 String<TInterval const *> interval_pointers;
974 _makePointerInterval(intervals, interval_pointers);
975
976 if (length(intervals) == 1) // if there is just one interval -> center = center of this interval
977 center = (rightBoundary(intervals[0]) - leftBoundary(intervals[0])) / (TValue)2.0;
978
979 _createIntervalTree(g, pm, interval_pointers, root, (TValue)0.0, center, length(intervals), tag);
980
981 reserve(pm, length(pm), Exact());
982 reserve(g.data_vertex, length(g.data_vertex), Exact());
983 }
984 }
985
986 // ComputeCenter tag as default construction method
987 template <typename TGraph, typename TPropertyMap, typename TIntervals>
988 inline void
989 createIntervalTree(TGraph & g, TPropertyMap & pm,
990 TIntervals & intervals,
991 typename Value<typename Value<TIntervals>::Type>::Type center)
992 {
993 createIntervalTree(g, pm, intervals, center, ComputeCenter());
994 }
995
996 template <typename TValue, typename TCargo, typename TIntervals, typename TSpec>
997 inline void
998 createIntervalTree(IntervalTree<TValue, TCargo> & it,
999 TIntervals & intervals,
1000 Tag<TSpec> const tag)
1001 {
1002 it.interval_counter = length(intervals);
1003 createIntervalTree(it.g, it.pm, intervals, tag);
1004 }
1005
1006 template <typename TValue, typename TCargo, typename TIntervals>
1007 inline void
1008 createIntervalTree(IntervalTree<TValue, TCargo> & it,
1009 TIntervals & intervals)
1010 {
1011 createIntervalTree(it, intervals, ComputeCenter());
1012 }
1013
1014 //////////////////////////////////////////////////////////////////////////////
1015 //remembers minimum and maximum of point values in intervals and sets the center
1016 //of each node to min+(max-min)/2
1017 template <typename TGraph, typename TPropertyMap, typename TIntervalPointer, typename TValue>
1018 inline void
1019 _createIntervalTree(TGraph & g, TPropertyMap & pm,
1020 String<TIntervalPointer *> & intervals,
1021 typename VertexDescriptor<TGraph>::Type & knot,
1022 TValue,
1023 TValue center,
1024 typename VertexDescriptor<TGraph>::Type len,
1025 Tag<TagComputeCenter_> const tag)
1026 {
1027 // Rekursionsanker
1028 if (len == 1)
1029 {
1030 _setIntervalTreeNode(value(pm, knot), center, *intervals[0]);
1031 return;
1032 }
1033
1034 typedef typename Value<TPropertyMap>::Type TNode;
1035 typedef typename ListType<TNode>::Type TList;
1036 typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
1037 typedef String<TIntervalPointer *> TIntervalPointers;
1038
1039 // one list of interval pointers for the intervals to the left of center
1040 TIntervalPointers S_left;
1041 // one list of interval pointers for the intervals to the right of center
1042 TIntervalPointers S_right;
1043
1044 TValue min1 = maxValue<TValue>();
1045 TValue min2 = maxValue<TValue>();
1046 TValue max1 = minValue<TValue>();
1047 TValue max2 = minValue<TValue>();
1048
1049 value(pm, knot).center = center;
1050
1051
1052 typedef typename Iterator<TIntervalPointers, Standard>::Type TIntervalIterator;
1053 TIntervalIterator it = begin(intervals, Standard());
1054 TIntervalIterator it_end = end(intervals, Standard());
1055
1056 // walk through intervals
1057 while (it != it_end)
1058 {
1059 // interval belongs to the left list
1060 if ((**it).i2 <= center)
1061 {
1062 appendValue(S_left, *it, Generous());
1063 //remember right most and left most point in left list
1064 if ((**it).i2 > max1)
1065 max1 = (**it).i2;
1066 if ((**it).i1 < min1)
1067 min1 = (**it).i1;
1068 }
1069 else
1070 {
1071 // interval belongs to the right list
1072 if ((**it).i1 > center)
1073 {
1074 appendValue(S_right, (*it), Generous());
1075 //remember right most and left most point in right list
1076 if ((**it).i2 > max2)
1077 max2 = (**it).i2;
1078 if ((**it).i1 < min2)
1079 min2 = (**it).i1;
1080 }
1081 else // interval belongs to this node
1082 {
1083 _appendIntervalTreeNodeLists(value(pm, knot), **it);
1084 }
1085 }
1086 ++it;
1087 }
1088
1089 // std::sort(begin(value(pm,knot).list1),end(value(pm,knot).list1),_less_compI1_ITree<typename Value<TList>::Type>);
1090 std::sort(begin(value(pm, knot).list2), end(value(pm, knot).list2), _greater_compI2_ITree<typename Value<TList>::Type>);
1091
1092 // build subtree to the left
1093 if (!empty(S_left))
1094 {
1095 TVertexDescriptor vd = addVertex(g);
1096 resize(pm, vd + 1);
1097 addEdge(g, knot, vd);
1098 _createIntervalTree(g, pm, S_left, vd, center, min1 + (max1 - min1) / 2, length(S_left), tag);
1099 }
1100 // build subtree to the right
1101 if (!empty(S_right))
1102 {
1103 TVertexDescriptor vd = addVertex(g);
1104 resize(pm, vd + 1);
1105 addEdge(g, knot, vd);
1106 _createIntervalTree(g, pm, S_right, vd, center, min2 + (max2 - min2) / 2, length(S_right), tag);
1107 }
1108 }
1109
1110 //////////////////////////////////////////////////////////////////////////////
1111 //createIntervalTree for all specs except CompCenter, the center value of each
1112 //node is determined by functions _calcIntervalTreeNodeCenterLeft and
1113 //_calcIntervalTreeNodeCenterRight
1114 template <typename TGraph, typename TPropertyMap, typename TSpec, typename TInterval, typename TValue>
1115 inline void
1116 _createIntervalTree(TGraph & g, TPropertyMap & pm,
1117 String<TInterval *> & intervals,
1118 typename VertexDescriptor<TGraph>::Type & knot,
1119 TValue last_center, TValue center,
1120 typename VertexDescriptor<TGraph>::Type len,
1121 Tag<TSpec> const tag)
1122 {
1123 SEQAN_CHECKPOINT
1124 // Rekursionsanker
1125 if (len == 1)
1126 {
1127 _setIntervalTreeNode(value(pm, knot), center, *value(intervals, 0));
1128 return;
1129 }
1130
1131 typedef typename Value<TPropertyMap>::Type TNode;
1132 typedef typename ListType<TNode>::Type TList;
1133 typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
1134 typedef String<TInterval *> TIntervalPointers;
1135
1136 // one list of interval pointers for the intervals to the left of center
1137 TIntervalPointers S_left;
1138 // one list of interval pointers for the intervals to the right of center
1139 TIntervalPointers S_right;
1140
1141 value(pm, knot).center = center;
1142
1143 typedef typename Iterator<TIntervalPointers, Standard>::Type TIntervalIterator;
1144 TIntervalIterator it = begin(intervals, Standard());
1145 TIntervalIterator it_end = end(intervals, Standard());
1146
1147 // walk through intervals
1148 while (it != it_end)
1149 {
1150 // interval belongs to the left list
1151 if ((**it).i2 <= center)
1152 {
1153 appendValue(S_left, *it, Generous());
1154 }
1155 else
1156 { // interval belongs to the right list
1157 if ((**it).i1 > center)
1158 appendValue(S_right, (*it), Generous());
1159 else
1160 // interval belongs to the current node
1161 _appendIntervalTreeNodeLists(value(pm, knot), **it);
1162 }
1163 ++it;
1164 }
1165
1166 // std::sort(begin(value(pm,knot).list1),end(value(pm,knot).list1),_less_compI1_ITree<typename Value<TList>::Type>);
1167 std::sort(begin(value(pm, knot).list2), end(value(pm, knot).list2), _greater_compI2_ITree<typename Value<TList>::Type>);
1168
1169 // build subtree to the left
1170 if (!empty(S_left))
1171 {
1172 TVertexDescriptor vd = addVertex(g);
1173 resizeVertexMap(pm, g);
1174 addEdge(g, knot, vd);
1175 TValue next_center = _calcIntervalTreeNodeCenterLeft(S_left, last_center, center, tag);
1176 _createIntervalTree(g, pm, S_left, vd, center, next_center, length(S_left), tag);
1177 }
1178 // build subtree to the right
1179 if (!empty(S_right))
1180 {
1181 TVertexDescriptor vd = addVertex(g);
1182 resizeVertexMap(pm, g);
1183 addEdge(g, knot, vd);
1184 TValue next_center = _calcIntervalTreeNodeCenterRight(S_right, last_center, center, tag);
1185 _createIntervalTree(g, pm, S_right, vd, center, next_center, length(S_right), tag);
1186 }
1187 }
1188
1189 // fill the container interval_pointers with pointers to the corresponding objects in intervals.
1190 // this is done to avoid copying and passing the whole IntervalAndCargo objects during interval tree construction
1191 template <typename TIntervals, typename TIntervalPointers>
1192 void
1193 _makePointerInterval(TIntervals & intervals, TIntervalPointers & interval_pointers)
1194 {
1195 typedef typename Iterator<TIntervals, Standard>::Type TIntervalIterator;
1196 typedef typename Iterator<TIntervalPointers, Standard>::Type TIntervalPointerIterator;
1197
1198 resize(interval_pointers, length(intervals));
1199 if (empty(intervals))
1200 return;
1201
1202 TIntervalIterator it = begin(intervals, Standard());
1203 TIntervalIterator itEnd = end(intervals, Standard());
1204 TIntervalPointerIterator iit = begin(interval_pointers, Standard());
1205
1206 for (; it != itEnd; ++it, ++iit)
1207 *iit = it;
1208 }
1209
1210 // if the center of the root is not given, it is placed in the "ComputeCenter way": in the middle of minValue and maxValue
1211 // where minValue is the minimum left boundary and maxValue is the maximum right boundary of all intervals
1212 template <typename TIntervals>
1213 typename Value<typename Value<TIntervals>::Type>::Type
1214 _calcIntervalTreeRootCenter(TIntervals & intervals)
1215 {
1216 SEQAN_CHECKPOINT;
1217
1218 SEQAN_ASSERT_GT(length(intervals), 0u);
1219
1220 typedef typename Value<typename Value<TIntervals>::Type>::Type TValue;
1221 typedef typename Iterator<TIntervals, Standard>::Type TIntervalIterator;
1222
1223 TIntervalIterator it = begin(intervals);
1224 TIntervalIterator it_end = end(intervals);
1225
1226 TValue min = maxValue<TValue>();
1227 TValue max = minValue<TValue>();
1228
1229 // get min and max
1230 while (it != it_end)
1231 {
1232 if (leftBoundary(*it) < min)
1233 min = leftBoundary(*it);
1234 if (rightBoundary(*it) > max)
1235 max = rightBoundary(*it);
1236 SEQAN_ASSERT_LEQ(min, max);
1237 ++it;
1238 }
1239
1240 SEQAN_ASSERT_LEQ(min, max);
1241
1242 // return middle between max and min
1243 return min + (max - min) / (TValue)2.0;
1244
1245 }
1246
1247 /*!
1248 * @fn IntervalTree#addInterval
1249 *
1250 * @headerfile <seqan/misc/interval_tree.h>
1251 *
1252 * @brief Adds an interval to an interval tree.
1253 *
1254 * @signature void addInterval(intervalTree, interval);
1255 * @signature void addInterval(intervalTree, begin, end[, cargo]);
1256 * @signature void addInterval(graph, propertyMap, interval);
1257 *
1258 * @param[in,out] intervalTree The interval tree to add the interval to. Types: @link IntervalTree @endlink.
1259 * @param[in] interval The interval to be added to the interval tree.
1260 * @param[in] begin Begin position of interval of type TValue.
1261 * @param[in] end End position of interval of type TValue.
1262 * @param[in] cargo Cargo to attach to the interval. Types: @link IntervalAndCargo @endlink.
1263 * @param[in,out] graph The directed graph that contains the topography of the interval tree.
1264 * @param[in,out] propertyMap The property map containing the node properties of the interval tree.
1265 */
1266
1267
1268 template <typename TGraph, typename TPropertyMap, typename TInterval>
1269 void
1270 addInterval(TGraph & g, TPropertyMap & pm, TInterval interval)
1271 {
1272 SEQAN_CHECKPOINT
1273
1274 typedef typename Iterator<TGraph, OutEdgeIterator>::Type TOutEdgeIterator;
1275 typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
1276 typedef typename Value<TPropertyMap>::Type TProperty;
1277 typedef typename Value<TInterval>::Type TValue;
1278 typedef typename ListType<TProperty>::Type TList;
1279
1280
1281 if (empty(pm))
1282 {
1283 TVertexDescriptor vd = addVertex(g);
1284 resizeVertexMap(pm, g);
1285 _setIntervalTreeNode(property(pm, vd), (rightBoundary(interval) + leftBoundary(interval)) / 2, interval);
1286 return;
1287
1288 }
1289 // start at root
1290 TVertexDescriptor act_knot = 0;
1291 TProperty act_prop = property(pm, act_knot);
1292 TProperty next_prop;
1293
1294 // look for the right node to add interval to
1295 while (true)
1296 {
1297 TOutEdgeIterator it(g, act_knot);
1298 act_prop = property(pm, act_knot);
1299 if (act_prop.center < leftBoundary(interval)) // interval to the left of current node?
1300 {
1301 if (atEnd(it))
1302 {
1303 TVertexDescriptor vd = addVertex(g);
1304 resizeVertexMap(pm, g);
1305 addEdge(g, act_knot, vd);
1306 _setIntervalTreeNode(property(pm, vd), (rightBoundary(interval) + leftBoundary(interval)) / (TValue)2.0, interval);
1307 break;
1308 }
1309 else
1310 {
1311 next_prop = property(pm, targetVertex(it));
1312 if (next_prop.center <= act_prop.center)
1313 {
1314 goNext(it);
1315 if (atEnd(it))
1316 {
1317 TVertexDescriptor vd = addVertex(g);
1318 resizeVertexMap(pm, g);
1319 addEdge(g, act_knot, vd);
1320 _setIntervalTreeNode(property(pm, vd), (rightBoundary(interval) + leftBoundary(interval)) / (TValue)2.0, interval);
1321 break;
1322 }
1323 }
1324 }
1325 act_knot = targetVertex(it);
1326 }
1327 else
1328 {
1329 if (rightBoundary(interval) <= act_prop.center) // interval to the right of current node?
1330 {
1331 if (atEnd(it))
1332 {
1333 TVertexDescriptor vd = addVertex(g);
1334 resizeVertexMap(pm, g);
1335 addEdge(g, act_knot, vd);
1336 _setIntervalTreeNode(property(pm, vd), (rightBoundary(interval) + leftBoundary(interval)) / 2, interval);
1337 break;
1338 }
1339 else
1340 {
1341 next_prop = property(pm, targetVertex(it));
1342 if (next_prop.center >= act_prop.center)
1343 {
1344 goNext(it);
1345 if (atEnd(it))
1346 {
1347 TVertexDescriptor vd = addVertex(g);
1348 resizeVertexMap(pm, g);
1349 addEdge(g, act_knot, vd);
1350 _setIntervalTreeNode(property(pm, vd), (rightBoundary(interval) + leftBoundary(interval)) / 2, interval);
1351 break;
1352 }
1353 }
1354 }
1355 act_knot = targetVertex(it);
1356 }
1357 else // need to create new node for interval
1358 {
1359 _appendIntervalTreeNodeLists(property(pm, act_knot), interval);
1360 std::sort(begin(property(pm, act_knot).list1), end(property(pm, act_knot).list1), _less_compI1_ITree<typename Value<TList>::Type>);
1361 std::sort(begin(property(pm, act_knot).list2), end(property(pm, act_knot).list2), _greater_compI2_ITree<typename Value<TList>::Type>);
1362 break;
1363 }
1364 }
1365 }
1366
1367 }
1368
1369 template <typename TValue, typename TCargo, typename TInterval>
1370 void
1371 addInterval(IntervalTree<TValue, TCargo> & itree, TInterval interval)
1372 {
1373 SEQAN_CHECKPOINT
1374
1375 ++ itree.interval_counter;
1376 addInterval(itree.g, itree.pm, interval);
1377
1378 }
1379
1380 template <typename TValue, typename TCargo>
1381 void
1382 addInterval(IntervalTree<TValue, TCargo> & itree, TValue begin, TValue end, TCargo cargo)
1383 {
1384 SEQAN_CHECKPOINT
1385
1386 IntervalAndCargo<TValue, TCargo> interval;
1387 interval.i1 = begin;
1388 interval.i2 = end;
1389 interval.cargo = cargo;
1390 ++itree.interval_counter;
1391 addInterval(itree.g, itree.pm, interval);
1392
1393 }
1394
1395 template <typename TValue, typename TCargo>
1396 void
1397 addInterval(IntervalTree<TValue, TCargo> & itree, TValue begin, TValue end)
1398 {
1399 SEQAN_CHECKPOINT
1400
1401 IntervalAndCargo<TValue, TCargo> interval;
1402 interval.i1 = begin;
1403 interval.i2 = end;
1404 interval.cargo = itree.interval_counter;
1405 ++itree.interval_counter;
1406 addInterval(itree.g, itree.pm, interval);
1407
1408 }
1409
1410 /*!
1411 * @fn IntervalTree#findIntervals
1412 * @brief Find all intervals that contain the query point or overlap with the query interval.
1413 *
1414 * @signature void findIntervals(result, intervalTree, query);
1415 * @signature void findIntervals(result, intervalTree, queryBegin, queryEnd);
1416 * @signature void findIntervals(result, graph, propertyMap, query);
1417 *
1418 * @param[out] result A reference to the result string of <tt>TCargo</tt> objects. Types: @link String @endlink.
1419 * @param[in] intervalTree An IntervalTree.
1420 * @param[in] graph The directed @link Graph graph @endlink that contains the topography of the interval tree.
1421 * @param[in] propertyMap The property map containing the node properties of the interval tree.
1422 * @param[in] query A query point.
1423 * @param[in] queryBegin The begin position of the query interval.
1424 * @param[in] queryEnd The end position of the query interval.
1425 */
1426
1427 template <typename TSpec, typename TPropertyMap, typename TValue, typename TCargo>
1428 inline void
1429 findIntervals(
1430 String<TCargo> & result,
1431 Graph<TSpec> const & g,
1432 TPropertyMap const & pm,
1433 TValue query)
1434 {
1435 typedef Graph<TSpec> const TGraph;
1436 typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
1437 typedef typename Value<TPropertyMap>::Type TProperty;
1438 typedef typename Value<TProperty>::Type TPropertyValue;
1439 typedef typename Iterator<TGraph, OutEdgeIterator>::Type TOutEdgeIterator;
1440
1441 resize(result, 0);
1442 if (empty(g))
1443 return;
1444
1445 // start at root
1446 TVertexDescriptor act_knot = 0;
1447 TProperty act_prop = property(pm, act_knot);
1448 TProperty next_prop;
1449
1450 while (true)
1451 {
1452 typename Iterator<Graph<TSpec>, OutEdgeIterator>::Type it7;
1453 Iter<Graph<TSpec>, GraphIterator<InternalOutEdgeIterator<OutEdgeIterator> > > it5(g, act_knot);
1454 TOutEdgeIterator it4;
1455 TOutEdgeIterator it(g, act_knot);
1456 act_prop = property(pm, act_knot);
1457 if (act_prop.center < (TPropertyValue)query) // look in current node and right subtree
1458 {
1459 unsigned int i = 0;
1460 while (i<length(act_prop.list2) && rightBoundary(value(act_prop.list2, i))>(TPropertyValue) query)
1461 {
1462 appendValue(result, cargo(value(act_prop.list2, i)), Generous());
1463 ++i;
1464 }
1465 if (atEnd(it))
1466 break;
1467
1468 next_prop = property(pm, targetVertex(it));
1469 if (next_prop.center <= act_prop.center)
1470 {
1471 goNext(it);
1472 if (atEnd(it))
1473 break;
1474 }
1475 act_knot = targetVertex(it);
1476 }
1477 else
1478 {
1479 if ((TPropertyValue)query < act_prop.center) // look in current node and left subtree
1480 {
1481 unsigned int i = 0;
1482 while (i < length(act_prop.list1) && leftBoundary(value(act_prop.list1, i)) <= (TPropertyValue)query)
1483 {
1484 appendValue(result, cargo(value(act_prop.list1, i)), Generous());
1485 ++i;
1486 }
1487 if (atEnd(it))
1488 break;
1489
1490 next_prop = property(pm, targetVertex(it));
1491 if (next_prop.center >= act_prop.center)
1492 {
1493 goNext(it);
1494 if (atEnd(it))
1495 break;
1496 }
1497 act_knot = targetVertex(it);
1498 }
1499 else // look in current node only, as query is center
1500 {
1501 for (unsigned int i = 0; i < length(act_prop.list1); ++i)
1502 appendValue(result, cargo(value(act_prop.list1, i)), Generous());
1503 break;
1504 }
1505 }
1506 }
1507
1508 }
1509
1510 template <typename TValue, typename TCargo, typename TValue2>
1511 inline void
1512 findIntervals(
1513 String<TCargo> & result,
1514 IntervalTree<TValue, TCargo> const & it,
1515 TValue2 query)
1516 {
1517 findIntervals(result, it.g, it.pm, query);
1518 }
1519
1520 template <typename TValue, typename TCargo, typename TValue2>
1521 inline void
1522 findIntervals(
1523 String<TCargo> & result,
1524 IntervalTree<TValue, TCargo> const & tree,
1525 TValue2 query_begin,
1526 TValue2 query_end)
1527 {
1528 findIntervals(result, tree.g, tree.pm, query_begin, query_end);
1529 }
1530
1531 template <typename TSpec, typename TPropertyMap, typename TValue, typename TCargo>
1532 inline void
1533 findIntervals(
1534 String<TCargo> & result,
1535 Graph<TSpec> const & g,
1536 TPropertyMap const & pm,
1537 TValue query_begin,
1538 TValue query_end)
1539 {
1540 typedef Graph<TSpec> const TGraph;
1541 typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
1542
1543 resize(result, 0);
1544
1545 // start at root
1546 TVertexDescriptor act_knot = 0;
1547 findIntervals(result, g, pm, act_knot, query_begin, query_end);
1548 }
1549
1550 template <
1551 typename TSpec,
1552 typename TPropertyMap,
1553 typename TVertexDescriptor,
1554 typename TValue,
1555 typename TCargo>
1556 inline void
1557 findIntervals(
1558 String<TCargo> & result,
1559 Graph<TSpec> const & g,
1560 TPropertyMap const & pm,
1561 TVertexDescriptor & act_knot,
1562 TValue query_begin,
1563 TValue query_end)
1564 {
1565 typedef typename Value<TPropertyMap>::Type TProperty;
1566 typedef typename Iterator<Graph<TSpec>, OutEdgeIterator>::Type TOutEdgeIterator;
1567
1568 if (empty(g))
1569 return;
1570
1571 TProperty act_prop = property(pm, act_knot);
1572 TProperty next_prop;
1573
1574 while (true)
1575 {
1576 TOutEdgeIterator it(g, act_knot);
1577 act_prop = property(pm, act_knot);
1578 //
1579 if (act_prop.center < query_begin) // query interval is to the right of node center
1580 {
1581 unsigned int i = 0;
1582 while (i<length(act_prop.list2) && rightBoundary(value(act_prop.list2, i))> query_begin)
1583 {
1584 appendValue(result, cargo(value(act_prop.list2, i)), Generous());
1585 ++i;
1586 }
1587 if (atEnd(it))
1588 break;
1589
1590 next_prop = property(pm, targetVertex(it));
1591 if (next_prop.center <= act_prop.center)
1592 {
1593 goNext(it);
1594 if (atEnd(it))
1595 break;
1596 }
1597 act_knot = targetVertex(it);
1598 }
1599 else
1600 {
1601 if (query_end <= act_prop.center) // query interval is to the left of node center
1602 {
1603 unsigned int i = 0;
1604 while (i < length(act_prop.list1) && leftBoundary(value(act_prop.list1, i)) < query_end)
1605 {
1606 appendValue(result, cargo(value(act_prop.list1, i)), Generous());
1607 ++i;
1608 }
1609
1610 if (atEnd(it))
1611 break;
1612
1613 next_prop = property(pm, targetVertex(it));
1614 if (next_prop.center >= act_prop.center)
1615 {
1616 goNext(it);
1617 if (atEnd(it))
1618 break;
1619 }
1620 act_knot = targetVertex(it);
1621 }
1622 else
1623 { //node center is contained in query interval
1624 for (unsigned int i = 0; i < length(act_prop.list1); ++i)
1625 appendValue(result, cargo(value(act_prop.list1, i)), Generous());
1626
1627 while (!atEnd(it))
1628 {
1629 TVertexDescriptor next_knot = targetVertex(it);
1630 findIntervals(result, g, pm, next_knot, query_begin, query_end);
1631 goNext(it);
1632 }
1633 break;
1634
1635 //break; //dont break! continue in both subtrees!!
1636 }
1637 }
1638 }
1639 }
1640
1641 /*!
1642 * @fn IntervalTree#findIntervalsExcludeTouching
1643 * @brief Find all intervals that contain the query point, exclude intervals that touch the query, i.e. where the query
1644 * point equals the start or end point.
1645 *
1646 * @signature void findIntervalsExcludeTouching(result, intervalTree, query);
1647 * @signature void findIntervalsExcludeTouching(result, graph, propertyMap, query,);
1648 *
1649 * @param[out] result The resulting string of cargos/ids of the intervals that contain the query point. Should
1650 * be a string of TCargo. Types: String
1651 * @param[in] intervalTree An interval tree Types: IntervalTree
1652 * @param[in] graph The directed graph that contains the topography of the interval tree.
1653 * @param[in] query The TValue to query here.
1654 * @param[in] propertyMap The property map containing the node properties of the interval tree
1655 */
1656
1657 template <typename TSpec, typename TPropertyMap, typename TValue, typename TCargo>
1658 inline void
1659 findIntervalsExcludeTouching(
1660 String<TCargo> & result,
1661 Graph<TSpec> const & g,
1662 TPropertyMap const & pm,
1663 TValue query)
1664 {
1665 typedef Graph<TSpec> const TGraph;
1666 typedef typename Iterator<TGraph, OutEdgeIterator>::Type TOutEdgeIterator;
1667 typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
1668 typedef typename Value<TPropertyMap>::Type TProperty;
1669
1670 resize(result, 0);
1671 if (empty(g))
1672 return;
1673
1674 // start at root
1675 TVertexDescriptor act_knot = 0;
1676 TProperty act_prop = property(pm, act_knot);
1677 TProperty next_prop;
1678
1679 while (true)
1680 {
1681 TOutEdgeIterator it(g, act_knot);
1682 act_prop = property(pm, act_knot);
1683 if ((TValue) act_prop.center < query) // look in current node and right subtree
1684 {
1685 int i = 0;
1686 while (i < (int)length(act_prop.list2) && (TValue)rightBoundary(value(act_prop.list2, i)) > query)
1687 {
1688 appendValue(result, cargo(value(act_prop.list2, i)), Generous());
1689 ++i;
1690 }
1691 if (atEnd(it))
1692 break;
1693
1694 next_prop = property(pm, targetVertex(it));
1695 if (next_prop.center <= act_prop.center)
1696 {
1697 goNext(it);
1698 if (atEnd(it))
1699 break;
1700 }
1701 act_knot = targetVertex(it);
1702 }
1703 else
1704 {
1705 if (query < (TValue) act_prop.center) // look in current node and left subtree
1706 {
1707 int i = 0;
1708 while (i < (int)length(act_prop.list1) && (TValue)leftBoundary(value(act_prop.list1, i)) < query)
1709 {
1710 appendValue(result, cargo(value(act_prop.list1, i)), Generous());
1711 ++i;
1712 }
1713 if (atEnd(it))
1714 break;
1715
1716 next_prop = property(pm, targetVertex(it));
1717 if (next_prop.center >= act_prop.center)
1718 {
1719 goNext(it);
1720 if (atEnd(it))
1721 break;
1722 }
1723 act_knot = targetVertex(it);
1724 }
1725 else // look in current node only
1726 {
1727 int i = 0;
1728 while (i < (int)length(act_prop.list1) && (TValue)leftBoundary(value(act_prop.list1, i)) < query)
1729 {
1730 appendValue(result, cargo(value(act_prop.list1, i)), Generous());
1731 ++i;
1732 }
1733 break;
1734 }
1735 }
1736 }
1737
1738 }
1739
1740 template <typename TValue, typename TCargo>
1741 inline void
1742 findIntervalsExcludeTouching(
1743 String<TCargo> & result,
1744 IntervalTree<TValue, TCargo> const & tree,
1745 TValue query)
1746 {
1747 findIntervalsExcludeTouching(result, tree.g, tree.pm, query);
1748 }
1749
1750 /*!
1751 * @fn IntervalTree#removeInterval
1752 * @brief Removes an interval from the interval tree.
1753 *
1754 * @signature bool removeInterval(intervalTree, iBegin, iEnd, iId);
1755 *
1756 * @param[in,out] intervalTree An interval tree Types: IntervalTree
1757 * @param[in] iBegin The begin position of the interval to be removed.
1758 * @param[in] iEnd The end position of the interval to be removed.
1759 * @param[in] iId The ID of the interval to be removed.
1760 *
1761 * @return bool <tt>true</tt> on success, <tt>false</tt> on failure.
1762 */
1763
1764 template <
1765 typename TSpec,
1766 typename TPropertyMap,
1767 typename TVertexDescriptor,
1768 typename TValue,
1769 typename TCargo>
1770 inline bool
1771 removeInterval(
1772 Graph<TSpec> & g,
1773 TPropertyMap & pm,
1774 TVertexDescriptor & act_knot,
1775 TValue i_begin,
1776 TValue i_end,
1777 TCargo i_id)
1778 {
1779 SEQAN_CHECKPOINT
1780
1781 typedef typename Value<TPropertyMap>::Type TProperty;
1782 typedef typename ListType<TProperty>::Type TList;
1783 typedef typename Iterator<TList, Standard>::Type TListIterator;
1784 typedef typename Iterator<Graph<TSpec>, OutEdgeIterator>::Type TOutEdgeIterator;
1785
1786 if (empty(g))
1787 return false;
1788
1789 TProperty act_prop = property(pm, act_knot);
1790 TProperty next_prop;
1791
1792 while (true)
1793 {
1794 TOutEdgeIterator it(g, act_knot);
1795 act_prop = property(pm, act_knot);
1796 //
1797 if (act_prop.center < i_begin) // interval is to the right of node center
1798 {
1799 if (atEnd(it))
1800 break;
1801
1802 next_prop = property(pm, targetVertex(it));
1803 if (next_prop.center <= act_prop.center)
1804 {
1805 goNext(it);
1806 if (atEnd(it))
1807 break;
1808 }
1809 act_knot = targetVertex(it);
1810 }
1811 else
1812 {
1813 if (i_end < act_prop.center) // interval is to the left of node center
1814 {
1815 if (atEnd(it))
1816 break;
1817
1818 next_prop = property(pm, targetVertex(it));
1819 if (next_prop.center >= act_prop.center)
1820 {
1821 goNext(it);
1822 if (atEnd(it))
1823 break;
1824 }
1825 act_knot = targetVertex(it);
1826 }
1827 else //node center is contained in interval, this is where we should find the interval to be removed
1828 { // remove from list1
1829 TProperty & change_prop = property(pm, act_knot);
1830 bool foundInLeft = false;
1831 TListIterator list_it_keep = begin(change_prop.list1, Standard());
1832 TListIterator list_it = begin(change_prop.list1, Standard());
1833 while (list_it != end(change_prop.list1, Standard()))
1834 {
1835 //std::cout << "Element: " << getLeftBoundary(*list_it) << ".." << getRightBoundary(*list_it) << " " << cargo(*list_it) << std::endl;
1836 if (getLeftBoundary(*list_it) == i_begin && cargo(*list_it) == i_id)
1837 {
1838 foundInLeft = true;
1839 ++list_it;
1840 continue;
1841 }
1842 *list_it_keep = *list_it;
1843 //std::cout << "Element keep: " << getLeftBoundary(*list_it_keep) << ".." << getRightBoundary(*list_it_keep) << " " << cargo(*list_it_keep) << std::endl;
1844 ++list_it; ++list_it_keep;
1845 }
1846
1847 bool foundInRight = false;
1848 list_it_keep = begin(change_prop.list2, Standard());
1849 list_it = begin(change_prop.list2, Standard());
1850 while (list_it != end(change_prop.list2, Standard()))
1851 {
1852 //std::cout << "Element: " << getLeftBoundary(*list_it) << ".." << getRightBoundary(*list_it) << " " << cargo(*list_it) << std::endl;
1853 if (getRightBoundary(*list_it) == i_end && cargo(*list_it) == i_id)
1854 {
1855 foundInRight = true;
1856 ++list_it;
1857 continue;
1858 }
1859 *list_it_keep = *list_it;
1860 //std::cout << "Element keep: " << getLeftBoundary(*list_it_keep) << ".." << getRightBoundary(*list_it_keep) << " " << cargo(*list_it_keep) << std::endl;
1861 ++list_it; ++list_it_keep;
1862
1863 }
1864
1865
1866 // TODO: if node is empty and does not have child nodes --> remove node.
1867 // keeping these empty leaf nodes just takes space unnecessarily
1868 if (foundInRight && foundInLeft)
1869 {
1870 resize(change_prop.list2, length(change_prop.list2) - 1);
1871 resize(change_prop.list1, length(change_prop.list1) - 1);
1872 return true;
1873 }
1874
1875 }
1876 }
1877 }
1878 return false;
1879 }
1880
1881 template <
1882 typename TSpec,
1883 typename TPropertyMap,
1884 typename TValue,
1885 typename TCargo>
1886 inline bool
1887 removeInterval(
1888 Graph<TSpec> & g,
1889 TPropertyMap & pm,
1890 TValue i_begin,
1891 TValue i_end,
1892 TCargo i_id)
1893 {
1894 typedef typename VertexDescriptor<Graph<TSpec> >::Type TVertexDescriptor;
1895
1896 // start looking at root
1897 TVertexDescriptor act_knot = 0;
1898 return removeInterval(g, pm, act_knot, i_begin, i_end, i_id);
1899 }
1900
1901 template <typename TValue, typename TCargo>
1902 inline bool
1903 removeInterval(
1904 IntervalTree<TValue, TCargo> & tree,
1905 TValue i_begin,
1906 TValue i_end,
1907 TCargo i_id)
1908 {
1909 return removeInterval(tree.g, tree.pm, i_begin, i_end, i_id);
1910 // we do not decrease the interval_counter of tree, as it would mix up interval IDs
1911 }
1912
1913 /////////////////// Metafunctions ///////////////////////
1914
1915 template <typename TValue, typename TCargo>
1916 struct Value<IntervalTree<TValue, TCargo> >
1917 {
1918 typedef TValue Type;
1919 };
1920
1921
1922 template <typename TValue, typename TCargo>
1923 struct Cargo<IntervalTree<TValue, TCargo> >
1924 {
1925 typedef TCargo Type;
1926 };
1927
1928 } // namespace SEQAN_NAMESPACE_MAIN
1929
1930 #endif //#ifndef SEQAN_MISC_INTERVAL_TREE_H
1931