File: FrugalBooMap.hpp

package info (click to toggle)
rapmap 0.15.0%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 6,228 kB
  • sloc: cpp: 48,810; ansic: 4,686; sh: 215; python: 82; makefile: 15
file content (323 lines) | stat: -rw-r--r-- 11,924 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
#ifndef __BOO_MAP_FRUGAL__
#define __BOO_MAP_FRUGAL__

#include "BooPHF.hpp"
#include "RapMapUtils.hpp"

#include "cereal/types/vector.hpp"
#include "cereal/types/utility.hpp"
#include "cereal/archives/binary.hpp"

#include <fstream>
#include <vector>
#include <iterator>
#include <type_traits>

#include <sys/stat.h>


// adapted from :
// http://stackoverflow.com/questions/34875315/implementation-my-own-list-and-iterator-stl-c
template <typename Iter, typename IndexT, typename HashMapT>
class KeyProxyIterator {
public:
    typedef KeyProxyIterator<Iter, IndexT, HashMapT> self_type;
    typedef uint64_t value_type;//std::iterator_traits<Iter>::value_type::first_type value_type;
    typedef value_type& reference;
    typedef value_type* pointer;
    typedef std::forward_iterator_tag iterator_category;
    typedef int64_t difference_type;

    KeyProxyIterator(Iter first, HashMapT* hm) : 
        curr_(first), hm_(hm) {}
    KeyProxyIterator operator++() { KeyProxyIterator i = *this; curr_++; return i; }
    KeyProxyIterator operator++(int) { ++curr_; return *this; }
    reference operator*() { 
        intRep_ = hm_->getKmerFromPos_(*curr_, mer_);
        return intRep_; 
    }
    /*
    pointer operator->() { 
        return &(curr_->first); 
    }
    */
    bool operator==(const self_type& rhs) { return curr_ == rhs.curr_; }
    bool operator!=(const self_type& rhs) { return curr_ != rhs.curr_; }
    bool operator<(const self_type& rhs) { return curr_ < rhs.curr_; }
    bool operator<=(const self_type& rhs) { return curr_ <= rhs.curr_; }
    
private:
    rapmap::utils::my_mer mer_;
    uint64_t intRep_;
    Iter curr_;
    HashMapT* hm_{nullptr}; 
};


template <typename IterT, typename ProxyValueT>
class KVProxy {
public:
    typedef KVProxy<IterT, ProxyValueT> self_type;
    typedef typename std::iterator_traits<IterT>::value_type ValueT;
    typedef std::pair<uint64_t, ValueT> value_type;
    typedef std::pair<uint64_t, ProxyValueT>& reference;
    typedef std::pair<uint64_t, ProxyValueT>* pointer;

    KVProxy(uint64_t mer, IterT it, ValueT len, bool isEnd = false) : curr_(it) {
        if(!isEnd) {ProxyValueT x{*it, len}; pair_ = std::make_pair(mer, x);}
    }
    reference operator*() { return pair_; }
    pointer operator->() { return &pair_; }
    bool operator==(const self_type& rhs) { return curr_ == rhs.curr_; }
    bool operator!=(const self_type& rhs) { return curr_ != rhs.curr_; }
    bool operator<(const self_type& rhs) { return curr_ < rhs.curr_; }
    bool operator<=(const self_type& rhs) { return curr_ <= rhs.curr_; }
private:
    IterT curr_;
    std::pair<uint64_t, ProxyValueT> pair_;
};


// Unlike the standard "generic" BooMap, the frugal variant
// *does not* store the key.  Rather, it assumes we have a
// pointer to the suffix array, and it "spot checks" the index 
// returned by the perfect hash by ensuring that the suffix at
// the corresponding offset starts with the query k-mer.
template <typename KeyT, typename ValueT>
class FrugalBooMap {
public:
    using self_type = FrugalBooMap<KeyT, ValueT>;
    using HasherT = boomphf::SingleHashFunctor<KeyT>;
    using BooPHFT = boomphf::mphf<KeyT, HasherT>;
    typedef typename ValueT::index_type IndexT;
    using IteratorT = KVProxy<typename std::vector<IndexT>::iterator, ValueT>;

    //using IteratorT = typename std::vector<std::pair<KeyT, ValueT>>::iterator;

    FrugalBooMap() : built_(false) {}
    void setSAPtr(std::vector<IndexT>* saPtr) { saPtr_ = saPtr; }
    void setTextPtr(const char* txtPtr, size_t textLen) { txtPtr_ = txtPtr; textLen_ = textLen; }

    void add(KeyT&& k, ValueT&& v) {
        // In the frugal map, we don't even keep the key!
        data_.emplace_back(v.begin());
        IndexT l = v.end() - v.begin();
        if (l >= std::numeric_limits<uint8_t>::max()) {
            overflow_[v.begin()] = l;
            lens_.emplace_back(std::numeric_limits<uint8_t>::max());
        } else {
            lens_.emplace_back(static_cast<uint8_t>(l));
        }
    }


    bool validate_hash(){
        for( auto& e : data_ ) {
            rapmap::utils::my_mer kmer(txtPtr_ + (*saPtr_)[e]);
            auto ind = boophf_->lookup(kmer.word(0));
            if (ind >= data_.size()) { 
                rapmap::utils::my_mer km(txtPtr_ + (*saPtr_)[e]);
                std::cerr << "index for " << km << " was " << ind << ", outside bounds of data_ (" << data_.size() << ")\n";
                return false;
            }
            auto mer = getKmerFromInterval_(e);
            if (mer != getKmerFromInterval_(data_[ind]) ) {
                std::cerr << "lookup of " << mer << " failed!\n";
            }
        }
        return true;
    }


    bool build(int nthreads=1) {
        size_t numElem = data_.size();
        KeyProxyIterator<decltype(data_.begin()), IndexT, self_type> kb(data_.begin(), this);
        KeyProxyIterator<decltype(data_.begin()), IndexT, self_type> ke(data_.end(), this);
        auto keyIt = boomphf::range(kb, ke);
        BooPHFT* ph = new BooPHFT(numElem, keyIt, nthreads);
        boophf_.reset(ph);
        std::cerr << "reordering keys and values to coincide with phf ... ";
        reorder_fn_();
        //std::cerr << "validating hash\n";
        //validate_hash();
        std::cerr << "done\n";
        std::cerr << "size of overflow table is " << overflow_.size() << '\n';
        built_ = true;
        return built_;
    }

    inline IteratorT find(const KeyT& k) {
        auto intervalIndex = boophf_->lookup(k);
        if (intervalIndex >= data_.size()) return end();
        auto ind = data_[intervalIndex];
        auto textInd = (*saPtr_)[ind];
        rapmap::utils::my_mer m(txtPtr_ + textInd);

        // If what we find matches the key, return the iterator
        // otherwise we don't have the key (it must have been here if it
        // existed).
        if (m.word(0) == k) {
            IndexT l = *(lens_.begin() + intervalIndex);
            if (l == std::numeric_limits<uint8_t>::max()) {
                l = overflow_[ind];
            }
            return IteratorT(m.word(0), data_.begin() + intervalIndex, ind + l);
        }
        return end();
    }
    
    /**
     * NOTE: This function *assumes* that the key is in the hash.
     * If it isn't, you'll get back a random element!
     */
    /*
    inline ValueT& operator[](const KeyT& k) {
        auto ind = boophf_->lookup(k);
        return (ind < data_.size() ? data_[ind].second : data_[0].second);
    }
    */
    
    inline IteratorT begin() { return IteratorT(0, data_.begin(), lens_.front()); }
    inline IteratorT end() { return IteratorT(0, data_.end(), 0, true); }
    inline IteratorT cend() const { return IteratorT(0, data_.cend(), 0, true); }
    inline IteratorT cbegin() const { return IteratorT(0, data_.cbegin(), lens_.front()); }
    
    void save(const std::string& ofileBase) {
        if (built_) {
            std::string hashFN = ofileBase + ".bph";
            // save the perfect hash function
            {
                std::ofstream os(hashFN, std::ios::binary);
                if (!os.is_open()) {
                    std::cerr << "BooM: unable to open output file [" << hashFN << "]; exiting!\n";
                    std::exit(1);
                }
                boophf_->save(os);
                os.close();
            }
            // and the values
            std::string dataFN = ofileBase + ".val";
            {
                std::ofstream valStream(dataFN, std::ios::binary);
                if (!valStream.is_open()) {
                    std::cerr << "BooM: unable to open output file [" << dataFN << "]; exiting!\n";
                    std::exit(1);
                }
                {
                    cereal::BinaryOutputArchive outArchive(valStream);
                    outArchive(data_);
                    outArchive(lens_);
                    overflow_.serialize(typename spp_utils::pod_hash_serializer<IndexT, IndexT>(), &valStream);
                }
                valStream.close();
            }
        }
    }
    
    void load(const std::string& ofileBase) {
        std::string hashFN = ofileBase + ".bph";
        std::string dataFN = ofileBase + ".val";

        if ( !FileExists_(hashFN.c_str()) ) {
            std::cerr << "BooM: Looking for perfect hash function file [" << hashFN << "], which doesn't exist! exiting.\n";
            std::exit(1);
        }
        if ( !FileExists_(dataFN.c_str()) ) {
            std::cerr << "BooM: Looking for key-value file [" << dataFN << "], which doesn't exist! exiting.\n";
            std::exit(1);
        }

        // load the perfect hash function
        {
            boophf_.reset(new BooPHFT);
            std::ifstream is(hashFN, std::ios::binary);
            boophf_->load(is);
            is.close();
        }
        // and the values
        {
            std::ifstream dataStream(dataFN, std::ios::binary);
            {
                cereal::BinaryInputArchive inArchive(dataStream);
                inArchive(data_);
                inArchive(lens_);
                overflow_.unserialize(typename spp_utils::pod_hash_serializer<IndexT, IndexT>(), &dataStream);
            }
            dataStream.close();
        }

        built_ = true;
    }

    inline KeyT getKmerFromInterval_(ValueT& ival) {
        rapmap::utils::my_mer m;// copy the global mer to get k-mer object
        m.fromChars(txtPtr_ + (*saPtr_)[ival.begin()]);
        return m.word(0);
    }

    // variant where we provide an existing mer object
    inline KeyT getKmerFromInterval_(ValueT& ival, rapmap::utils::my_mer& m) {
        m.fromChars(txtPtr_ + (*saPtr_)[ival.begin()]);
        return m.word(0);
    }

    // variant where we provide an existing mer object
    inline KeyT getKmerFromPos_(IndexT pos, rapmap::utils::my_mer& m) {
        m.fromChars(static_cast<decltype(txtPtr_)>(txtPtr_ + (*saPtr_)[pos]));
        return m.word(0);
    }

private:
    // Taken from http://stackoverflow.com/questions/12774207/fastest-way-to-check-if-a-file-exist-using-standard-c-c11-c
    bool FileExists_(const char *path) {
        struct stat fileStat;
        if ( stat(path, &fileStat) ) {
            return false;
        }
        if ( !S_ISREG(fileStat.st_mode) ) {
            return false;
        }
        return true;
    }

    void reorder_fn_()  {
        /* Adapted from code at: http://blog.merovius.de/2014/08/12/applying-permutation-in-constant.html */
        // Note, we can actually do this with out the bitvector by using the high-order bit 
        // of the start of the suffix array intervals (since they are signed integers and negative
        // positions are forbidden). 
        rapmap::utils::my_mer mer;
        std::vector<bool> bits(data_.size(), false);
        for ( size_t i = 0; i < data_.size(); ++i ) {
            if (!bits[i]) {
                decltype(data_.front()) v = data_[i];
                decltype(lens_.front()) v2 = lens_[i];
                auto j = boophf_->lookup(getKmerFromPos_(data_[i], mer));
                while (i != j) {
                    auto pj = boophf_->lookup(getKmerFromPos_(data_[j], mer));
                    std::swap(data_[j], v);
                    std::swap(lens_[j], v2);
                    bits[j] = 1;
                    j = pj; 
                }
                data_[i] = v;
                lens_[i] = v2;
            }
        }
    }

    std::vector<IndexT>* saPtr_;
    const char* txtPtr_; 
    size_t textLen_;
    rapmap::utils::my_mer mer_;
    bool built_;
    // Starting offset in the suffix array
    std::vector<IndexT> data_;
    // Length of the interval
    std::vector<uint8_t> lens_;
    // Overflow table if interval is >= std::numeric_limits<uint8_t>::max()
    spp::sparse_hash_map<IndexT, IndexT> overflow_;
    std::unique_ptr<BooPHFT> boophf_{nullptr};
};


#endif // __BOO_MAP_FRUGAL__