Imported Debian patch 0.12.7-3
[bowtie.git] / row_chaser.h
1 /*
2  * row_chaser.h
3  */
4
5 #ifndef ROW_CHASER_H_
6 #define ROW_CHASER_H_
7
8 #include <iostream>
9 #include "seqan/sequence.h"
10 #include "ebwt.h"
11 #include "aligner_metrics.h"
12
13 /**
14  * A class that statefully converts a row index to a reference
15  * location.  There is a large memory-latency penalty usually
16  * associated with calling the Ebwt object's mapLF method, which this
17  * object does repeatedly in order to resolve the reference offset.
18  * The "statefulness" in how the computation is organized here allows
19  * some or all of that penalty to be hidden using prefetching.
20  */
21 template<typename TStr>
22 class RowChaser {
23
24         typedef std::pair<uint32_t,uint32_t> U32Pair;
25         typedef Ebwt<TStr> TEbwt;
26
27 public:
28         RowChaser(AlignerMetrics *metrics = NULL) :
29                 done(false),
30                 prepped_(false),
31                 ebwt_(NULL),
32                 qlen_(0),
33                 eh_(NULL),
34                 row_(0xffffffff),
35                 jumps_(0),
36                 sideloc_(),
37                 off_(0xffffffff),
38                 tlen_(0),
39                 metrics_(metrics)
40         { }
41
42         /**
43          * Convert a row to a joined reference offset.  This has to be
44          * converted to understand where it is w/r/t the reference hit and
45          * offset within it.
46          */
47         static uint32_t toFlatRefOff(const TEbwt* ebwt, uint32_t qlen, uint32_t row) {
48                 RowChaser rc;
49                 rc.setRow(row, qlen, ebwt);
50                 while(!rc.done) {
51                         rc.advance();
52                 }
53                 return rc.flatOff();
54         }
55
56         /**
57          * Convert a row to a reference offset.
58          */
59         static U32Pair toRefOff(const TEbwt* ebwt, uint32_t qlen, uint32_t row) {
60                 RowChaser rc;
61                 rc.setRow(row, qlen, ebwt);
62                 while(!rc.done) {
63                         rc.advance();
64                 }
65                 return rc.off();
66         }
67
68         /**
69          * Set the next row for us to "chase" (i.e. map to a reference
70          * location using the BWT step-left operation).
71          */
72         void setRow(uint32_t row, uint32_t qlen, const TEbwt* ebwt) {
73                 assert_neq(0xffffffff, row);
74                 assert_gt(qlen, 0);
75                 assert(ebwt != NULL);
76                 ebwt_ = ebwt;
77                 eh_ = &ebwt->_eh;
78                 row_ = row;
79                 qlen_ = qlen;
80                 ASSERT_ONLY(sideloc_.invalidate());
81                 if(row_ == ebwt_->_zOff) {
82                         // We arrived at the extreme left-hand end of the reference
83                         off_ = 0;
84                         done = true;
85                         return;
86                 } else if((row_ & eh_->_offMask) == row_) {
87                         // We arrived at a marked row
88                         off_ = ebwt_->_offs[row_ >> eh_->_offRate];
89                         done = true;
90                         return;
91                 }
92                 done = false;
93                 jumps_ = 0;
94                 off_ = 0xffffffff;
95                 prepped_ = false;
96                 prep();
97         }
98
99         /**
100          * Advance the step-left process by one step.  Check if we're done.
101          */
102         void advance() {
103                 // Advance by 1
104                 assert(!done);
105                 while(!done) {
106                         assert(prepped_);
107                         prepped_ = false;
108                         if(metrics_ != NULL) metrics_->curBwtOps_++;
109                         uint32_t newrow = ebwt_->mapLF(sideloc_);
110                         ASSERT_ONLY(sideloc_.invalidate());
111                         jumps_++;
112                         assert_neq(newrow, row_);
113                         // Update row_ field
114                         row_ = newrow;
115                         if(row_ == ebwt_->_zOff) {
116                                 // We arrived at the extreme left-hand end of the reference
117                                 off_ = jumps_;
118                                 done = true;
119                         } else if((row_ & eh_->_offMask) == row_) {
120                                 // We arrived at a marked row
121                                 off_ = ebwt_->_offs[row_ >> eh_->_offRate] + jumps_;
122                                 done = true;
123                         }
124                         prep();
125                 }
126         }
127
128         /**
129          * Prepare for the next call to advance() by prefetching the
130          * appropriate portions of the index.  The caller should make sure
131          * that the
132          */
133         void prep() {
134                 if(!done) {
135                         assert(!prepped_);
136                         assert(!sideloc_.valid());
137                         assert_leq(row_, eh_->_len);
138                         sideloc_.initFromRow(row_, *eh_, (const uint8_t*)ebwt_->_ebwt);
139                         assert(sideloc_.valid());
140                 }
141                 prepped_ = true;
142         }
143
144         /**
145          * Get the calculated offset.  This has to be converted with a call
146          * to Ebwt::joinedToTextOff() to understand where it is w/r/t the
147          * reference hit and offset within it.
148          */
149         uint32_t flatOff() const {
150                 return off_;
151         }
152
153         /**
154          * Get the calculated offset.
155          */
156         U32Pair off() {
157                 uint32_t off = flatOff();
158                 assert_neq(0xffffffff, off);
159                 uint32_t tidx;
160                 uint32_t textoff = 0xffffffff;
161                 ebwt_->joinedToTextOff(qlen_, off, tidx, textoff, tlen_);
162                 // Note: tidx may be 0xffffffff, if alignment overlaps a
163                 // reference boundary
164                 return make_pair(tidx, textoff);
165         }
166
167         uint32_t tlen() const {
168                 return tlen_;
169         }
170
171         bool done;               /// true = chase is done & answer is in off_
172         bool prepped_; /// true = prefetch is issued and it's OK to call advance()
173
174 protected:
175
176         const TEbwt* ebwt_;      /// index to resolve row in
177         uint32_t qlen_;          /// length of read; needed to convert to ref. coordinates
178         const EbwtParams* eh_;   /// eh field from index
179         uint32_t row_;           /// current row
180         uint32_t jumps_;         /// # steps so far
181         SideLocus sideloc_;      /// current side locus
182         uint32_t off_;           /// calculated offset (0xffffffff if not done)
183         uint32_t tlen_;          /// hit text length
184         AlignerMetrics *metrics_;
185 };
186
187 #endif /* ROW_CHASER_H_ */