← Back to branch summary

~ubuntu-branches/ubuntu/saucy/abyss/saucy

« back to all changes in this revision

Viewing changes to FMIndex/FMIndex.h

  • Committer: Package Import Robot
  • Author(s): Shaun Jackman
  • Date: 2011-09-11 10:00:13 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: package-import@ubuntu.com-20110911100013-oa4m5fi036mjuwc8
Tags: 1.3.0-1
* New upstream release.
* Add libboost-graph-dev to Build-Depends.

Show diffs side-by-side

added added

removed removed

Lines of Context:
FMIndex/FMIndex.h
 
1
#ifndef FMINDEX_H
 
2
#define FMINDEX_H 1
 
3
 
 
4
#include "config.h"
 
5
#include "BitArrays.h"
 
6
#include "IOUtil.h"
 
7
#include "sais.hxx"
 
8
#include <boost/integer.hpp>
 
9
#include <algorithm>
 
10
#include <cassert>
 
11
#include <cstdlib> // for exit
 
12
#include <iostream>
 
13
#include <limits> // for numeric_limits
 
14
#include <stdint.h>
 
15
#include <string>
 
16
#include <vector>
 
17
 
 
18
/** A match of a substring of a query sequence to an FM index. */
 
19
struct FMInterval
 
20
{
 
21
        FMInterval(size_t l, size_t u, unsigned qstart, unsigned qend)
 
22
                : l(l), u(u), qstart(qstart), qend(qend) { }
 
23
        size_t l, u;
 
24
        unsigned qstart, qend;
 
25
 
 
26
        unsigned qspan() const
 
27
        {
 
28
                assert(qstart <= qend);
 
29
                return qend - qstart;
 
30
        }
 
31
};
 
32
 
 
33
/** A match of a substring of a query sequence to a target sequence.
 
34
 */
 
35
struct Match
 
36
{
 
37
        unsigned qstart, qend;
 
38
        size_t tstart;
 
39
        unsigned count;
 
40
        Match(unsigned qstart, unsigned qend,
 
41
                        size_t tstart, unsigned count)
 
42
                : qstart(qstart), qend(qend), tstart(tstart), count(count) { }
 
43
 
 
44
        unsigned qspan() const
 
45
        {
 
46
                assert(qstart <= qend);
 
47
                return qend - qstart;
 
48
        }
 
49
};
 
50
 
 
51
/** An FM index. */
 
52
class FMIndex
 
53
{
 
54
        /** A symbol. */
 
55
        typedef uint8_t T;
 
56
 
 
57
        /** The sentinel symbol. */
 
58
        static T SENTINEL() { return std::numeric_limits<T>::max(); }
 
59
 
 
60
  public:
 
61
        /** An index. */
 
62
        typedef boost::uint_t<FMBITS>::exact size_type;
 
63
 
 
64
        /** An index for SAIS, which must be signed. */
 
65
        typedef boost::int_t<FMBITS>::exact sais_size_type;
 
66
 
 
67
        /** The type of a symbol. */
 
68
        typedef T value_type;
 
69
 
 
70
        void read(const char *path, std::vector<T> &s);
 
71
 
 
72
        FMIndex() : m_sampleSA(1) { }
 
73
 
 
74
/** Return the size of the string not counting the sentinel. */
 
75
size_t size() const { return m_occ.size() - 1; }
 
76
 
 
77
/** Build an FM-index of the specified file. */
 
78
template<typename It>
 
79
void assign(It first, It last)
 
80
{
 
81
        assert(first < last);
 
82
        assert(size_t(last - first)
 
83
                        < std::numeric_limits<size_type>::max());
 
84
 
 
85
        m_sampleSA = 1;
 
86
 
 
87
        // Translate the alphabet.
 
88
        if (m_alphabet.empty())
 
89
                setAlphabet(first, last);
 
90
        std::transform(first, last, first, Translate(*this));
 
91
        std::replace(first, last, SENTINEL(), T(0));
 
92
 
 
93
        std::cerr << "The alphabet has "
 
94
                << m_alphabet.size() << " symbols.\n";
 
95
 
 
96
        // Construct the suffix array.
 
97
        std::cerr << "Building the suffix array...\n";
 
98
        size_t n = last - first;
 
99
        m_sa.resize(n + 1);
 
100
        m_sa[0] = n;
 
101
 
 
102
        assert(sizeof (size_type) == sizeof (sais_size_type));
 
103
        int status = saisxx(first,
 
104
                        reinterpret_cast<sais_size_type*>(&m_sa[1]),
 
105
                        (sais_size_type)n,
 
106
                        (sais_size_type)m_alphabet.size());
 
107
        assert(status == 0);
 
108
        if (status != 0)
 
109
                abort();
 
110
 
 
111
        // Construct the Burrows-Wheeler transform.
 
112
        std::vector<T> bwt;
 
113
        std::cerr << "Building the Burrows-Wheeler transform...\n";
 
114
        bwt.resize(m_sa.size());
 
115
        for (size_t i = 0; i < m_sa.size(); i++)
 
116
                bwt[i] = m_sa[i] == 0 ? SENTINEL() : first[m_sa[i] - 1];
 
117
 
 
118
        std::cerr << "Building the character occurence table...\n";
 
119
        m_occ.assign(bwt);
 
120
        countOccurences();
 
121
}
 
122
 
 
123
/** Sample the suffix array. */
 
124
void sampleSA(unsigned period)
 
125
{
 
126
        assert(period > 0);
 
127
        if (period == m_sampleSA)
 
128
                return;
 
129
        assert(m_sampleSA == 1);
 
130
        m_sampleSA = period;
 
131
        if (m_sampleSA == 1)
 
132
                return;
 
133
        std::vector<size_type>::iterator out = m_sa.begin();
 
134
        for (size_t i = 0; i < m_sa.size(); i += m_sampleSA)
 
135
                *out++ = m_sa[i];
 
136
        m_sa.erase(out, m_sa.end());
 
137
        assert(!m_sa.empty());
 
138
}
 
139
 
 
140
/** Return the position of the specified suffix in the original
 
141
 * string.
 
142
 */
 
143
size_t locate(size_t i) const
 
144
{
 
145
        size_t n = 0;
 
146
        while (i % m_sampleSA != 0) {
 
147
                T c = m_occ.at(i);
 
148
                i = c == SENTINEL() ? 0 : m_cf[c] + m_occ.rank(c, i + 1) - 1;
 
149
                n++;
 
150
        }
 
151
        assert(i / m_sampleSA < m_sa.size());
 
152
        size_t pos = m_sa[i / m_sampleSA] + n;
 
153
        return pos < m_occ.size() ? pos : pos - m_occ.size();
 
154
}
 
155
 
 
156
        /** Decompress the index. */
 
157
        template <typename It>
 
158
        void decompress(It out)
 
159
        {
 
160
                // Construct the original string.
 
161
                std::vector<T> s;
 
162
                for (size_t i = 0;;) {
 
163
                        assert(i < m_occ.size());
 
164
                        T c = m_occ.at(i);
 
165
                        if (c == SENTINEL())
 
166
                                break;
 
167
                        s.push_back(c);
 
168
                        i = m_cf[c] + m_occ.rank(c, i);
 
169
                        assert(i > 0);
 
170
                }
 
171
 
 
172
                // Translate the character set and output the result.
 
173
                for (std::vector<T>::const_reverse_iterator it = s.rbegin();
 
174
                                it != s.rend(); ++it) {
 
175
                        T c = *it;
 
176
                        assert(c < m_alphabet.size());
 
177
                        *out++ = m_alphabet[c];
 
178
                }
 
179
        }
 
180
 
 
181
        /** Search for an exact match. */
 
182
        template <typename It>
 
183
        std::pair<size_t, size_t> findExact(It first, It last) const
 
184
        {
 
185
                assert(first < last);
 
186
                size_t l = 1, u = m_occ.size();
 
187
                It it;
 
188
                for (it = last - 1; it >= first && l < u; --it) {
 
189
                        T c = *it;
 
190
                        if (c == SENTINEL())
 
191
                                return std::make_pair(0, 0);
 
192
                        l = m_cf[c] + m_occ.rank(c, l);
 
193
                        u = m_cf[c] + m_occ.rank(c, u);
 
194
                }
 
195
                return std::make_pair(l, u);
 
196
        }
 
197
 
 
198
        /** Search for an exact match. */
 
199
        template <typename String>
 
200
        std::pair<size_t, size_t> findExact(const String& q) const
 
201
        {
 
202
                String s(q.size());
 
203
                std::transform(q.begin(), q.end(), s.begin(),
 
204
                                Translate(*this));
 
205
                return findExact(s.begin(), s.end());
 
206
        }
 
207
 
 
208
        /** Search for a matching suffix of the query. */
 
209
        template <typename It>
 
210
        FMInterval findSuffix(It first, It last) const
 
211
        {
 
212
                assert(first < last);
 
213
                size_t l = 1, u = m_occ.size();
 
214
                It it;
 
215
                for (it = last - 1; it >= first && l < u; --it) {
 
216
                        T c = *it;
 
217
                        if (c == SENTINEL())
 
218
                                break;
 
219
                        size_t l1 = m_cf[c] + m_occ.rank(c, l);
 
220
                        size_t u1 = m_cf[c] + m_occ.rank(c, u);
 
221
                        if (l1 >= u1)
 
222
                                break;
 
223
                        l = l1;
 
224
                        u = u1;
 
225
                }
 
226
                return FMInterval(l, u, it - first + 1, last - first);
 
227
        }
 
228
 
 
229
        /** Search for a matching substring of the query at least k long.
 
230
         * @return the longest match
 
231
         */
 
232
        template <typename It>
 
233
        FMInterval findSubstring(It first, It last, unsigned k) const
 
234
        {
 
235
                assert(first < last);
 
236
                FMInterval best(0, 0, 0, k > 0 ? k - 1 : 0);
 
237
                for (It it = last; it > first; --it) {
 
238
                        if (unsigned(it - first) < best.qspan())
 
239
                                return best;
 
240
                        FMInterval interval = findSuffix(first, it);
 
241
                        if (interval.qspan() > best.qspan())
 
242
                                best = interval;
 
243
                }
 
244
                return best;
 
245
        }
 
246
 
 
247
        /** Translate from ASCII to the indexed alphabet. */
 
248
        struct Translate {
 
249
                Translate(const FMIndex& fmIndex) : m_fmIndex(fmIndex) { }
 
250
                T operator()(unsigned char c) const
 
251
                {
 
252
                        return c < m_fmIndex.m_mapping.size()
 
253
                                ? m_fmIndex.m_mapping[c] : SENTINEL();
 
254
                }
 
255
          private:
 
256
                const FMIndex& m_fmIndex;
 
257
        };
 
258
 
 
259
        /** Search for a matching substring of the query at least k long.
 
260
         * @return the longest match and the number of matches
 
261
         */
 
262
        Match find(const std::string& q, unsigned k) const
 
263
        {
 
264
                std::string s = q;
 
265
                std::transform(s.begin(), s.end(), s.begin(),
 
266
                                Translate(*this));
 
267
 
 
268
                FMInterval interval = findSubstring(s.begin(), s.end(), k);
 
269
                assert(interval.l <= interval.u);
 
270
                size_t count = interval.u - interval.l;
 
271
                if (count == 0)
 
272
                        return Match(0, 0, 0, 0);
 
273
                return Match(interval.qstart, interval.qend,
 
274
                                locate(interval.l), count);
 
275
        }
 
276
 
 
277
        /** Set the alphabet to [first, last). */
 
278
        template <typename It>
 
279
        void setAlphabet(It first, It last)
 
280
        {
 
281
                std::vector<bool> mask(std::numeric_limits<T>::max());
 
282
                for (It it = first; it < last; ++it) {
 
283
                        assert((size_t)*it < mask.size());
 
284
                        mask[*it] = true;
 
285
                }
 
286
 
 
287
                m_alphabet.clear();
 
288
                m_mapping.clear();
 
289
                for (unsigned c = 0; c < mask.size(); ++c) {
 
290
                        if (!mask[c])
 
291
                                continue;
 
292
                        m_mapping.resize(c + 1, std::numeric_limits<T>::max());
 
293
                        m_mapping[c] = m_alphabet.size();
 
294
                        m_alphabet.push_back(c);
 
295
                }
 
296
                assert(!m_alphabet.empty());
 
297
        }
 
298
 
 
299
        /** Set the alphabet to the characters of s. */
 
300
        void setAlphabet(const std::string& s)
 
301
        {
 
302
                setAlphabet(s.begin(), s.end());
 
303
        }
 
304
 
 
305
#define STRINGIFY(X) #X
 
306
#define FM_VERSION_BITS(BITS) "FM " STRINGIFY(BITS) " 1"
 
307
#define FM_VERSION FM_VERSION_BITS(FMBITS)
 
308
 
 
309
/** Store an index. */
 
310
friend std::ostream& operator<<(std::ostream& out, const FMIndex& o)
 
311
{
 
312
        out << FM_VERSION << '\n'
 
313
                << o.m_sampleSA << '\n';
 
314
 
 
315
        out << o.m_alphabet.size() << '\n';
 
316
        out.write((char*)o.m_alphabet.data(),
 
317
                        o.m_alphabet.size() * sizeof o.m_alphabet[0]);
 
318
 
 
319
        out << o.m_sa.size() << '\n';
 
320
        out.write((char*)o.m_sa.data(),
 
321
                o.m_sa.size() * sizeof o.m_sa[0]);
 
322
 
 
323
        return out << o.m_occ;
 
324
}
 
325
 
 
326
/** Load an index. */
 
327
friend std::istream& operator>>(std::istream& in, FMIndex& o)
 
328
{
 
329
        std::string version;
 
330
        std::getline(in, version);
 
331
        assert(in);
 
332
        if (version != FM_VERSION) {
 
333
                std::cerr << "error: the version of this FM-index, `"
 
334
                        << version << "', does not match the version required "
 
335
                        "by this program, `" FM_VERSION "'.\n";
 
336
                exit(EXIT_FAILURE);
 
337
        }
 
338
 
 
339
        in >> o.m_sampleSA;
 
340
        assert(in);
 
341
 
 
342
        size_t n;
 
343
        char c;
 
344
 
 
345
        in >> n;
 
346
        assert(in);
 
347
        c = in.get();
 
348
        assert(c == '\n');
 
349
        assert(n < std::numeric_limits<size_type>::max());
 
350
        o.m_alphabet.resize(n);
 
351
        in.read((char*)o.m_alphabet.data(), n * sizeof o.m_alphabet[0]);
 
352
        o.setAlphabet(o.m_alphabet.begin(), o.m_alphabet.end());
 
353
 
 
354
        in >> n;
 
355
        assert(in);
 
356
        c = in.get();
 
357
        assert(c == '\n');
 
358
        assert(n < std::numeric_limits<size_type>::max());
 
359
        o.m_sa.resize(n);
 
360
        in.read((char*)o.m_sa.data(), n * sizeof o.m_sa[0]);
 
361
 
 
362
        in >> o.m_occ;
 
363
        assert(in);
 
364
        o.countOccurences();
 
365
 
 
366
        return in;
 
367
}
 
368
 
 
369
private:
 
370
 
 
371
/** Build the cumulative frequency table m_cf from m_occ. */
 
372
void countOccurences()
 
373
{
 
374
        assert(!m_alphabet.empty());
 
375
        m_cf.resize(m_alphabet.size());
 
376
        // The sentinel character occurs once.
 
377
        m_cf[0] = 1;
 
378
        for (unsigned i = 0; i < m_cf.size() - 1; ++i)
 
379
                m_cf[i + 1] = m_cf[i] + m_occ.count(i);
 
380
}
 
381
 
 
382
        unsigned m_sampleSA;
 
383
        std::vector<T> m_alphabet;
 
384
        std::vector<T> m_mapping;
 
385
        std::vector<size_type> m_cf;
 
386
        std::vector<size_type> m_sa;
 
387
        BitArrays m_occ;
 
388
};
 
389
 
 
390
#endif