Imported Debian patch 0.12.7-3
[bowtie.git] / range_cache.h
1 /*
2  * range_cache.h
3  *
4  * Classes that encapsulate the caching of
5  */
6
7 #ifndef RANGE_CACHE_H_
8 #define RANGE_CACHE_H_
9
10 #include <stdint.h>
11 #include <utility>
12 #include <iostream>
13 #include <stdexcept>
14 #include <map>
15 #include "ebwt.h"
16 #include "row_chaser.h"
17
18 #define RANGE_NOT_SET 0xffffffff
19 #define RANGE_CACHE_BAD_ALLOC 0xffffffff
20
21 /**
22  * Manages a pool of memory used exclusively for range cache entries.
23  * This manager is allocate-only; it exists mainly so that we can avoid
24  * lots of new[]s and delete[]s.
25  *
26  * A given stretch of words may be one of two types: a cache entry, or
27  * a cache entry wrapper.  A cache entry has a length and a list of
28  * already-resolved reference positions.  A cache entry wrapper has a
29  * pointer to a cache entry for a different range, along with an
30  * integer indicating how many "jumps" to the left that range is from
31  * the one that owns the wrapper.
32  */
33 class RangeCacheMemPool {
34 public:
35         RangeCacheMemPool(uint32_t lim /* max cache size in bytes */) :
36                 lim_(lim >> 2 /* convert to words */), occ_(0), buf_(NULL),
37                 closed_(false)
38         {
39                 if(lim_ > 0) {
40                         try {
41                                 buf_ = new uint32_t[lim_];
42                                 if(buf_ == NULL) throw std::bad_alloc();
43                         } catch(std::bad_alloc& e) {
44                                 cerr << "Allocation error allocating " << lim
45                                          << " words of range-cache memory" << endl;
46                                 throw 1;
47                         }
48                         assert(buf_ != NULL);
49                         // Fill with 1s to signal that these elements are
50                         // uninitialized
51                         memset(buf_, 0xff, lim_ << 2 /* convert back to bytes */);
52                 }
53         }
54
55         ~RangeCacheMemPool() {
56                 // Release all word memory!
57                 if(lim_ > 0) delete[] buf_;
58         }
59
60         /**
61          * Allocate numElts elements from the word pool.
62          */
63         uint32_t alloc(uint32_t numElts) {
64                 assert_gt(numElts, 0);
65                 assert_leq(occ_, lim_);
66                 if(occ_ + numElts > lim_ || numElts >= 0x80000000) {
67                         return RANGE_CACHE_BAD_ALLOC;
68                 }
69                 assert_gt(lim_, 0);
70                 uint32_t ret = occ_;
71                 assert(allocs_.find(ret) == allocs_.end());
72                 ASSERT_ONLY(allocs_.insert(ret));
73                 // Clear the first elt so that we don't think there's already
74                 // something there
75 #ifndef NDEBUG
76                 for(size_t i = 0; i < numElts; i++) {
77                         assert_eq(0xffffffff, buf_[occ_ + i]);
78                 }
79 #endif
80                 buf_[occ_] = 0;
81                 occ_ += numElts;
82                 assert_leq(occ_, lim_);
83                 if(lim_ - occ_ < 10) {
84                         // No more room - don't try anymore
85                         closed_ = true;
86                 }
87                 return ret;
88         }
89
90         /**
91          * Turn a pool-array index into a pointer; check that it doesn't
92          * fall outside the pool first.
93          */
94         inline uint32_t *get(uint32_t off) {
95                 assert_gt(lim_, 0);
96                 assert_lt(off, lim_);
97                 assert(allocs_.find(off) != allocs_.end());
98                 uint32_t *ret = buf_ + off;
99                 assert_neq(0x80000000, ret[0]);
100                 assert_neq(0xffffffff, ret[0]);
101                 return ret;
102         }
103
104         /**
105          * Return true iff there's no more room in the cache.
106          */
107         inline bool closed() {
108                 return closed_;
109         }
110
111 private:
112         uint32_t lim_;  /// limit on number of 32-bit words to dish out in total
113         uint32_t occ_;  /// number of occupied words
114         uint32_t *buf_; /// buffer of 32-bit words
115         bool closed_;   ///
116 #ifndef NDEBUG
117         std::set<uint32_t> allocs_; // elements allocated
118 #endif
119 };
120
121 /**
122  * A view to a range of cached reference positions.
123  */
124 class RangeCacheEntry {
125
126         typedef Ebwt<String<Dna> > TEbwt;
127         typedef std::pair<uint32_t,uint32_t> U32Pair;
128         typedef RowChaser<String<Dna> > TRowChaser;
129
130 public:
131         /**
132          *
133          */
134         RangeCacheEntry(bool sanity = false) :
135                 top_(0xffffffff), jumps_(0), len_(0), ents_(NULL), ebwt_(NULL),
136                 sanity_(sanity)
137         { }
138
139         /**
140          * Create a new RangeCacheEntry from the data in the pool at 'ents'.
141          */
142         RangeCacheEntry(RangeCacheMemPool& pool, uint32_t top,
143                         uint32_t ent, TEbwt* ebwt, bool sanity = false) :
144             sanity_(sanity)
145         {
146                 init(pool, top, ent, ebwt);
147         }
148
149         /**
150          * Initialize a RangeCacheEntry from the data in the pool at 'ents'.
151          */
152         void init(RangeCacheMemPool& pool, uint32_t top, uint32_t ent, TEbwt* ebwt) {
153                 assert(ebwt != NULL);
154                 top_ = top;
155                 ebwt_ = ebwt;
156                 uint32_t *ents = pool.get(ent);
157                 assert_neq(0x80000000, ents[0]);
158                 // Is hi bit set?
159                 if((ents[0] & 0x80000000) != 0) {
160                         // If so, the target is a wrapper and the non-hi bits
161                         // contain the # jumps
162                         jumps_ = (ents[0] & ~0x80000000);
163                         assert_gt(jumps_, 0);
164                         assert_leq(jumps_, ebwt_->_eh._len);
165                         // Get the target entry
166                         uint32_t *dest = pool.get(ents[1]);
167                         // Get the length from the target entry
168                         len_ = dest[0];
169                         assert_leq(top_ + len_, ebwt_->_eh._len);
170                         assert_gt(len_, 0);
171                         assert_leq(len_, ebwt_->_eh._len);
172                         // Get the pointer to the entries themselves
173                         ents_ = dest + 1;
174                 } else {
175                         // Not a wrapper, so there are no jumps
176                         jumps_ = 0;
177                         // Get the length from the target entry
178                         len_  = ents[0];
179                         assert_leq(top_ + len_, ebwt_->_eh._len);
180                         assert_gt(len_, 0);
181                         assert_leq(len_, ebwt_->_eh._len);
182                         // Get the pointer to the entries themselves
183                         ents_ = ents + 1;
184                 }
185                 assert(sanityCheckEnts());
186         }
187
188         /**
189          * Initialize a wrapper with given number of jumps and given target
190          * entry index.
191          */
192         void init(RangeCacheMemPool& pool, uint32_t top, uint32_t jumps,
193                   uint32_t ent, TEbwt* ebwt)
194         {
195                 assert(ebwt != NULL);
196                 ebwt_ = ebwt;
197                 top_ = top;
198                 jumps_ = jumps;
199                 uint32_t *ents = pool.get(ent);
200                 // Must not be a wrapper
201                 assert_eq(0, ents[0] & 0x80000000);
202                 // Get the length from the target entry
203                 len_ = ents[0];
204                 assert_gt(len_, 0);
205                 assert_leq(len_, ebwt_->_eh._len);
206                 // Get the pointer to the entries themselves
207                 ents_ = ents + 1;
208                 assert_leq(top_ + len_, ebwt_->_eh._len);
209                 assert(sanityCheckEnts());
210         }
211
212         uint32_t len() const   {
213                 assert(ents_ != NULL);
214                 assert(ebwt_ != NULL);
215                 return len_;
216         }
217
218         uint32_t jumps() const {
219                 assert(ents_ != NULL);
220                 assert(ebwt_ != NULL);
221                 return jumps_;
222         }
223
224         /**
225          *
226          */
227         void reset() {
228                 ents_ = NULL;
229         }
230
231         /**
232          * Return true iff this object represents a valid cache entry.
233          */
234         bool valid() const {
235                 return ents_ != NULL;
236         }
237
238         TEbwt *ebwt() {
239                 return ebwt_;
240         }
241
242         /**
243          * Install a result obtained by a client of this cache; be sure to
244          * adjust for how many jumps down the tunnel the cache entry is
245          * situated.
246          */
247         void install(uint32_t elt, uint32_t val) {
248                 if(ents_ == NULL) {
249                         // This is not a valid cache entry; do nothing
250                         return;
251                 }
252                 assert(ents_ != NULL);
253                 assert(ebwt_ != NULL);
254                 assert_leq(jumps_, val);
255                 assert_neq(0xffffffff, val);
256                 assert_leq(top_ + len_, ebwt_->_eh._len);
257                 if(elt < len_) {
258                         val -= jumps_;
259                         if(verbose_) cout << "Installed reference offset: " << (top_ + elt) << endl;
260                         ASSERT_ONLY(uint32_t sanity = TRowChaser::toFlatRefOff(ebwt_, 1, top_ + elt));
261                         assert_eq(sanity, val);
262 #ifndef NDEBUG
263                         for(size_t i = 0; i < len_; i++) {
264                                 if(i == elt) continue;
265                                 assert_neq(val, ents_[i]);
266                         }
267 #endif
268                         ents_[elt] = val;
269                 } else {
270                         // ignore install request
271                         if(verbose_) cout << "Fell off end of cache entry for install: " << (top_ + elt) << endl;
272                 }
273         }
274
275         /**
276          * Get an element from the cache, adjusted for tunnel jumps.
277          */
278         inline uint32_t get(uint32_t elt) const {
279                 if(ents_ == NULL) {
280                         // This is not a valid cache entry; do nothing
281                         return RANGE_NOT_SET;
282                 }
283                 assert(ents_ != NULL);
284                 assert(ebwt_ != NULL);
285                 assert_leq(top_ + len_, ebwt_->_eh._len);
286                 if(elt < len_ && ents_[elt] != RANGE_NOT_SET) {
287                         if(verbose_) cout << "Retrieved result from cache: " << (top_ + elt) << endl;
288                         uint32_t ret = ents_[elt] + jumps_;
289                         ASSERT_ONLY(uint32_t sanity = TRowChaser::toFlatRefOff(ebwt_, 1, top_ + elt));
290                         assert_eq(sanity, ret);
291                         return ret;
292                 } else {
293                         if(verbose_) cout << "Cache entry not set: " << (top_ + elt) << endl;
294                         return RANGE_NOT_SET;
295                 }
296         }
297
298         /**
299          * Check that len_ and the ents_ array both make sense.
300          */
301         static bool sanityCheckEnts(uint32_t len, uint32_t *ents, TEbwt* ebwt) {
302                 assert_gt(len, 0);
303                 assert_leq(len, ebwt->_eh._len);
304                 if(len < 10) {
305                         for(size_t i = 0; i < len; i++) {
306                                 if(ents[i] == 0xffffffff) continue;
307                                 assert_leq(ents[i], ebwt->_eh._len);
308                                 for(size_t j = i+1; j < len; j++) {
309                                         if(ents[j] == 0xffffffff) continue;
310                                         assert_neq(ents[i], ents[j]);
311                                 }
312                         }
313                 } else {
314                         std::set<uint32_t> seen;
315                         for(size_t i = 0; i < len; i++) {
316                                 if(ents[i] == 0xffffffff) continue;
317                                 assert(seen.find(ents[i]) == seen.end());
318                                 seen.insert(ents[i]);
319                         }
320                 }
321                 return true;
322         }
323
324         /**
325          * Check that len_ and the ents_ array both make sense.
326          */
327         bool sanityCheckEnts() {
328                 return RangeCacheEntry::sanityCheckEnts(len_, ents_, ebwt_);
329         }
330
331 private:
332
333         uint32_t top_;   /// top pointer for this range
334         uint32_t jumps_; /// how many tunnel-jumps it is away from the requester
335         uint32_t len_;   /// # of entries in cache entry
336         uint32_t *ents_; /// ptr to entries, which are flat offs within joined ref
337         //U32Pair *ents_;  /// pointer to entries, which are tidx,toff pairs
338         TEbwt    *ebwt_; /// index that alignments are in
339         bool     verbose_; /// be talkative?
340         bool     sanity_;  /// do consistency checks?
341 };
342
343 /**
344  *
345  */
346 class RangeCache {
347
348         typedef Ebwt<String<Dna> > TEbwt;
349         typedef std::vector<uint32_t> TU32Vec;
350         typedef std::map<uint32_t, uint32_t> TMap;
351         typedef std::map<uint32_t, uint32_t>::iterator TMapItr;
352
353 public:
354         RangeCache(uint32_t lim, TEbwt* ebwt) :
355                 lim_(lim), map_(), pool_(lim), closed_(false), ebwt_(ebwt), sanity_(true) { }
356
357         /**
358          * Given top and bot offsets, retrieve the canonical cache entry
359          * that best covers that range.  The cache entry may not directly
360          * correspond to the top offset provided, rather, it might be an
361          * entry that lies "at the end of the tunnel" when top and bot are
362          * walked backward.
363          */
364         bool lookup(uint32_t top, uint32_t bot, RangeCacheEntry& ent) {
365                 if(ebwt_ == NULL || lim_ == 0) return false;
366                 assert_gt(bot, top);
367                 ent.reset();
368                 TMapItr itr = map_.find(top);
369                 if(itr == map_.end()) {
370                         // No cache entry for the given 'top' offset
371                         if(closed_) {
372                                 return false; // failed to get cache entry
373                         } else {
374                                 if(pool_.closed()) {
375                                         closed_ = true;
376                                         return false; // failed to get cache entry
377                                 }
378                         }
379                         // Use the tunnel
380                         bool ret = tunnel(top, bot, ent);
381                         return ret;
382                 } else {
383                         // There is a cache entry for the given 'top' offset
384                         uint32_t ret = itr->second;
385                         ent.init(pool_, top, ret, ebwt_);
386                         return true; // success
387                 }
388         }
389
390         /**
391          * Exhaustively check all entries linked to from map_ to ensure
392          * they're well-formed.
393          */
394         bool repOk() {
395 #ifndef NDEBUG
396                 for(TMapItr itr = map_.begin(); itr != map_.end(); itr++) {
397                         uint32_t top = itr->first;
398                         uint32_t idx = itr->second;
399                         uint32_t jumps = 0;
400                         assert_leq(top, ebwt_->_eh._len);
401                         uint32_t *ents = pool_.get(idx);
402                         if((ents[0] & 0x80000000) != 0) {
403                                 jumps = ents[0] & ~0x80000000;
404                                 assert_leq(jumps, ebwt_->_eh._len);
405                                 idx = ents[1];
406                                 ents = pool_.get(idx);
407                         }
408                         uint32_t len = ents[0];
409                         assert_leq(top + len, ebwt_->_eh._len);
410                         RangeCacheEntry::sanityCheckEnts(len, ents + 1, ebwt_);
411                 }
412 #endif
413                 return true;
414         }
415
416 protected:
417
418         /**
419          * Tunnel through to the first range that 1) includes all the same
420          * suffixes (though longer) as the given range, and 2) has a cache
421          * entry for it.
422          */
423         bool tunnel(uint32_t top, uint32_t bot, RangeCacheEntry& ent) {
424                 assert_gt(bot, top);
425                 TU32Vec tops;
426                 const uint32_t spread = bot - top;
427                 SideLocus tloc, bloc;
428                 SideLocus::initFromTopBot(top, bot, ebwt_->_eh, ebwt_->_ebwt, tloc, bloc);
429                 uint32_t newtop = top, newbot = bot;
430                 uint32_t jumps = 0;
431                 // Walk left through the tunnel
432                 while(true) {
433                         if(ebwt_->rowL(tloc) != ebwt_->rowL(bloc)) {
434                                 // Different characters at top and bot positions of
435                                 // BWT; this means that the calls to mapLF below are
436                                 // guaranteed to yield rows in two different character-
437                                 // sections of the BWT.
438                                 break;
439                         }
440                         // Advance top and bot
441                         newtop = ebwt_->mapLF(tloc);
442                         newbot = ebwt_->mapLF(bloc);
443                         assert_geq(newbot, newtop);
444                         assert_leq(newbot - newtop, spread);
445                         // If the new spread is the same as the old spread, we can
446                         // be confident that the new range includes all of the same
447                         // suffixes as the last range (though longer by 1 char)
448                         if((newbot - newtop) == spread) {
449                                 // Check if newtop is already cached
450                                 TMapItr itr = map_.find(newtop);
451                                 jumps++;
452                                 if(itr != map_.end()) {
453                                         // This range, which is further to the left in the
454                                         // same tunnel as the query range, has a cache
455                                         // entry already, so use that
456                                         uint32_t idx = itr->second;
457                                         uint32_t *ents = pool_.get(idx);
458                                         if((ents[0] & 0x80000000) != 0) {
459                                                 // The cache entry we found was a wrapper; make
460                                                 // a new wrapper that points to that wrapper's
461                                                 // target, with the appropriate number of jumps
462                                                 jumps += (ents[0] & ~0x80000000);
463                                                 idx = ents[1];
464                                         }
465                                         // Allocate a new wrapper
466                                         uint32_t newentIdx = pool_.alloc(2);
467                                         if(newentIdx != RANGE_CACHE_BAD_ALLOC) {
468                                                 // We successfully made a new wrapper entry;
469                                                 // now populate it and install it in map_
470                                                 uint32_t *newent = pool_.get(newentIdx); // get ptr to it
471                                                 assert_eq(0, newent[0]);
472                                                 newent[0] = 0x80000000 | jumps; // set jumps
473                                                 newent[1] = idx;                // set target
474                                                 assert(map_.find(top) == map_.end());
475                                                 map_[top] = newentIdx;
476                                                 if(sanity_) assert(repOk());
477                                         }
478                                         // Initialize the entry
479                                         ent.init(pool_, top, jumps, idx, ebwt_);
480                                         return true;
481                                 }
482                                 // Save this range
483                                 tops.push_back(newtop);
484                                 SideLocus::initFromTopBot(newtop, newbot, ebwt_->_eh, ebwt_->_ebwt, tloc, bloc);
485                         } else {
486                                 // Not all the suffixes were preserved, so we can't
487                                 // link the source range's cached result to this
488                                 // range's cached results
489                                 break;
490                         }
491                         assert_eq(jumps, tops.size());
492                 }
493                 assert_eq(jumps, tops.size());
494                 // Try to create a new cache entry for the leftmost range in
495                 // the tunnel (which might be the query range)
496                 uint32_t newentIdx = pool_.alloc(spread + 1);
497                 if(newentIdx != RANGE_CACHE_BAD_ALLOC) {
498                         // Successfully allocated new range cache entry; install it
499                         uint32_t *newent = pool_.get(newentIdx);
500                         assert_eq(0, newent[0]);
501                         // Store cache-range length in first word
502                         newent[0] = spread;
503                         assert_lt(newent[0], 0x80000000);
504                         assert_eq(spread, newent[0]);
505                         uint32_t entTop = top;
506                         uint32_t jumps = 0;
507                         if(tops.size() > 0) {
508                                 entTop = tops.back();
509                                 jumps = tops.size();
510                         }
511                         // Cache the entry for the end of the tunnel
512                         assert(map_.find(entTop) == map_.end());
513                         map_[entTop] = newentIdx;
514                         if(sanity_) assert(repOk());
515                         ent.init(pool_, entTop, jumps, newentIdx, ebwt_);
516                         assert_eq(spread, newent[0]);
517                         if(jumps > 0) {
518                                 assert_neq(entTop, top);
519                                 // Cache a wrapper entry for the query range (if possible)
520                                 uint32_t wrapentIdx = pool_.alloc(2);
521                                 if(wrapentIdx != RANGE_CACHE_BAD_ALLOC) {
522                                         uint32_t *wrapent = pool_.get(wrapentIdx);
523                                         assert_eq(0, wrapent[0]);
524                                         wrapent[0] = 0x80000000 | jumps;
525                                         wrapent[1] = newentIdx;
526                                         assert(map_.find(top) == map_.end());
527                                         map_[top] = wrapentIdx;
528                                         if(sanity_) assert(repOk());
529                                 }
530                         }
531                         return true;
532                 } else {
533                         // Could not allocate new range cache entry
534                         return false;
535                 }
536         }
537
538         uint32_t lim_;           /// Total number of key/val bytes to keep in cache
539         TMap map_;               ///
540         RangeCacheMemPool pool_; /// Memory pool
541         bool closed_;            /// Out of space; no new entries
542         TEbwt* ebwt_;            /// Index that alignments are in
543         bool sanity_;
544 };
545
546 #endif /* RANGE_CACHE_H_ */