« back to all changes in this revision
Viewing changes to FMIndex/FMIndex.h
- browse files at revision 6
- compare with another revision
- download diff
- download tarball
- view history from revision 6
- Committer: Package Import Robot
- Author(s): Shaun Jackman
- Date: 2011-10-24 11:32:59 UTC
- mfrom: (1.1.2)
- Revision ID: package-import@ubuntu.com-20111024113259-qtid38pgm9p889ah
Tags: 1.3.1-1
New upstream release.
- files added:
- Common/BitUtil.h
- FilterGraph
- files removed:
- Graph/adjtodot.cc
- files modified:
- .pc/abyss-pe/bin/abyss-pe
added
removed
FMIndex/FMIndex.h
10
10
#include <cassert>
11
11
#include <cstdlib> // for exit
12
12
#include <iostream>
13
#include <iterator>
13
14
#include <limits> // for numeric_limits
14
15
#include <stdint.h>
15
16
#include <string>
16
17
#include <vector>
17
18
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
19
/** An FM index. */
52
20
class FMIndex
53
21
{
59
27
60
28
public:
61
29
/** An index. */
62
typedef boost::uint_t<FMBITS>::exact size_type;
30
typedef boost::uint_t<FMBITS>::least size_type;
63
31
64
32
/** An index for SAIS, which must be signed. */
65
typedef boost::int_t<FMBITS>::exact sais_size_type;
33
typedef boost::int_t<FMBITS>::least sais_size_type;
66
34
67
35
/** The type of a symbol. */
68
36
typedef T value_type;
69
37
70
void read(const char *path, std::vector<T> &s);
71
72
FMIndex() : m_sampleSA(1) { }
38
/** A suffix array interval. */
39
struct SAInterval
40
{
41
size_type l, u;
42
SAInterval(size_type l, size_type u) : l(l), u(u) { }
43
SAInterval(const FMIndex& fm) : l(1), u(fm.m_occ.size()) { }
44
bool empty() const { return l >= u; }
45
bool operator==(const SAInterval& x) const
46
{
47
return x.l == l && x.u == u;
48
}
49
50
size_type size() const
51
{
52
assert(l <= u);
53
return u - l;
54
}
55
};
56
57
/** A match of a substring of a query sequence to an FM index. */
58
struct Match : public SAInterval
59
{
60
unsigned qstart, qend;
61
62
Match(size_type l, size_type u, unsigned qstart, unsigned qend)
63
: SAInterval(l, u), qstart(qstart), qend(qend) { }
64
65
unsigned qspan() const
66
{
67
assert(qstart <= qend);
68
return qend - qstart;
69
}
70
};
71
72
FMIndex() : m_sampleSA(1) { }
73
73
74
74
/** Return the size of the string not counting the sentinel. */
75
75
size_t size() const { return m_occ.size() - 1; }
76
76
77
/** The size of the alphabet. */
78
unsigned alphabetSize() const { return m_alphabet.size(); }
79
77
80
/** Build an FM-index of the specified file. */
78
81
template<typename It>
79
82
void assign(It first, It last)
137
140
assert(!m_sa.empty());
138
141
}
139
142
140
/** Return the position of the specified suffix in the original
141
* string.
142
*/
143
size_t locate(size_t i) const
143
/** Return the specified element of the suffix array. */
144
size_t at(size_t i) const
144
145
{
146
assert(i < m_occ.size());
145
147
size_t n = 0;
146
148
while (i % m_sampleSA != 0) {
147
149
T c = m_occ.at(i);
148
i = c == SENTINEL() ? 0 : m_cf[c] + m_occ.rank(c, i + 1) - 1;
150
i = c == SENTINEL() ? 0 : m_cf[c] + m_occ.rank(c, i);
149
151
n++;
150
152
}
151
153
assert(i / m_sampleSA < m_sa.size());
153
155
return pos < m_occ.size() ? pos : pos - m_occ.size();
154
156
}
155
157
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
}
158
/** Return the specified element of the suffix array. */
159
size_t operator[](size_t i) const
160
{
161
return at(i);
162
}
163
164
/** Return the symbol at the specifed position of the BWT. */
165
T bwtAt(size_t i) const
166
{
167
assert(i < m_occ.size());
168
T c = m_occ.at(i);
169
assert(c != SENTINEL());
170
assert(c < m_alphabet.size());
171
return m_alphabet[c];
172
}
173
174
/** Return the first symbol of the specified suffix. */
175
T symbolAt(size_t i) const
176
{
177
assert(i < m_occ.size());
178
std::vector<size_type>::const_iterator it
179
= std::upper_bound(m_cf.begin(), m_cf.end(), i);
180
assert(it != m_cf.begin());
181
T c = it - m_cf.begin() - 1;
182
assert(c < m_alphabet.size());
183
return m_alphabet[c];
184
}
185
186
/** Decompress the index. */
187
template <typename It>
188
void decompress(It out)
189
{
190
// Construct the original string.
191
std::vector<T> s;
192
for (size_t i = 0;;) {
193
assert(i < m_occ.size());
194
T c = m_occ.at(i);
195
if (c == SENTINEL())
196
break;
197
s.push_back(c);
198
i = m_cf[c] + m_occ.rank(c, i);
199
assert(i > 0);
200
}
201
202
// Translate the character set and output the result.
203
for (std::vector<T>::reverse_iterator it = s.rbegin();
204
it != s.rend(); ++it) {
205
T c = *it;
206
assert(c < m_alphabet.size());
207
*out++ = m_alphabet[c];
208
}
209
}
210
211
/** Extend a suffix array coordinate by one character to the left. */
212
size_type update(size_type i, T c) const
213
{
214
assert(c < m_cf.size());
215
return m_cf[c] + m_occ.rank(c, i);
216
}
217
218
/** Extend a suffix array interval by one character to the left. */
219
SAInterval update(SAInterval sai, T c) const
220
{
221
return SAInterval(update(sai.l, c), update(sai.u, c));
222
}
223
224
/** Search for an exact match. */
225
template <typename It>
226
SAInterval findExact(It first, It last, SAInterval sai) const
227
{
228
assert(first < last);
229
for (It it = last - 1; it >= first && !sai.empty(); --it) {
230
T c = *it;
231
if (c == SENTINEL())
232
return SAInterval(0, 0);
233
sai = update(sai, c);
234
}
235
return sai;
236
}
237
238
/** Search for an exact match. */
239
template <typename String>
240
SAInterval findExact(const String& q) const
241
{
242
String s(q.size());
243
std::transform(q.begin(), q.end(), s.begin(), Translate(*this));
244
return findExact(s.begin(), s.end());
245
}
246
247
/** Search for a suffix of the query that matches a prefix of the
248
* target.
249
*/
250
template <typename It, typename OutIt>
251
OutIt findOverlapSuffix(It first, It last, OutIt out,
252
unsigned minOverlap) const
253
{
254
assert(first < last);
255
256
SAInterval sai(*this);
257
for (It it = last - 1; it >= first && !sai.empty(); --it) {
258
T c = *it;
259
if (c == SENTINEL())
260
break;
261
sai = update(sai, c);
262
if (sai.empty())
263
break;
264
265
if (unsigned(last - it) < minOverlap)
266
continue;
267
268
SAInterval sai1 = update(sai, 0);
269
if (!sai1.empty())
270
*out++ = Match(sai1.l, sai1.u,
271
it - first, last - first);
272
}
273
return out;
274
}
275
276
/** Search for a suffix of the query that matches a prefix of the
277
* target.
278
*/
279
template <typename OutIt>
280
OutIt findOverlapSuffix(const std::string& q, OutIt out,
281
unsigned minOverlap) const
282
{
283
std::string s = q;
284
std::transform(s.begin(), s.end(), s.begin(), Translate(*this));
285
return findOverlapSuffix(s.begin(), s.end(), out, minOverlap);
286
}
287
288
/** Search for a prefix of the query that matches a suffix of the
289
* target.
290
*/
291
template <typename OutIt>
292
OutIt findOverlapPrefix(const std::string& q, OutIt out,
293
unsigned minOverlap) const
294
{
295
std::string s = q;
296
std::transform(s.begin(), s.end(), s.begin(), Translate(*this));
297
typedef std::string::const_iterator It;
298
It first = s.begin();
299
for (It it = first + minOverlap; it <= s.end(); ++it) {
300
SAInterval sai = findExact(first, it,
301
update(SAInterval(*this), 0));
302
if (!sai.empty())
303
*out++ = Match(sai.l, sai.u, 0, it - first);
304
}
305
return out;
306
}
307
308
/** Search for a matching suffix of the query. */
309
template <typename It, typename MemoIt>
310
Match findSuffix(It first, It last, MemoIt memoIt) const
311
{
312
assert(first < last);
313
314
SAInterval sai(*this);
315
It it;
316
for (it = last - 1; it >= first && !sai.empty(); --it) {
317
T c = *it;
318
if (c == SENTINEL())
319
break;
320
SAInterval sai1 = update(sai, c);
321
if (sai1.empty())
322
break;
323
sai = sai1;
324
325
if (*memoIt == sai) {
326
// This vertex of the prefix DAWG has been visited.
327
break;
328
}
329
*memoIt++ = sai;
330
}
331
return Match(sai.l, sai.u, it - first + 1, last - first);
332
}
333
334
/** Search for a matching substring of the query at least k long.
335
* @return the longest match
336
*/
337
template <typename It>
338
Match findSubstring(It first, It last, unsigned k) const
339
{
340
assert(first < last);
341
Match best(0, 0, 0, k > 0 ? k - 1 : 0);
342
343
// Record which vertices of the prefix DAWG have been visited.
344
std::vector<SAInterval> memo(last - first, SAInterval(0, 0));
345
SAInterval* memoIt = &memo[0];
346
for (It it = last; it > first; --it) {
347
if (unsigned(it - first) < best.qspan())
348
return best;
349
Match interval = findSuffix(first, it, memoIt++);
350
if (interval.qspan() > best.qspan())
351
best = interval;
352
}
353
return best;
354
}
355
356
/** Translate from ASCII to the indexed alphabet. */
357
struct Translate {
358
Translate(const FMIndex& fmIndex) : m_fmIndex(fmIndex) { }
359
T operator()(unsigned char c) const
360
{
361
return c < m_fmIndex.m_mapping.size()
362
? m_fmIndex.m_mapping[c] : SENTINEL();
363
}
364
private:
365
const FMIndex& m_fmIndex;
366
};
367
368
/** Search for a matching substring of the query at least k long.
369
* @return the longest match and the number of matches
370
*/
371
Match find(const std::string& q, unsigned k) const
372
{
373
std::string s = q;
374
std::transform(s.begin(), s.end(), s.begin(), Translate(*this));
375
376
Match m = findSubstring(s.begin(), s.end(), k);
377
return m;
378
}
379
380
/** Set the alphabet to [first, last). */
381
template <typename It>
382
void setAlphabet(It first, It last)
383
{
384
std::vector<bool> mask(std::numeric_limits<T>::max());
385
for (It it = first; it < last; ++it) {
386
assert((size_t)*it < mask.size());
387
mask[*it] = true;
388
}
389
390
m_alphabet.clear();
391
m_mapping.clear();
392
for (unsigned c = 0; c < mask.size(); ++c) {
393
if (!mask[c])
394
continue;
395
m_mapping.resize(c + 1, std::numeric_limits<T>::max());
396
m_mapping[c] = m_alphabet.size();
397
m_alphabet.push_back(c);
398
}
399
assert(!m_alphabet.empty());
400
}
401
402
/** Set the alphabet to the characters of s. */
403
void setAlphabet(const std::string& s)
404
{
405
setAlphabet(s.begin(), s.end());
406
}
304
407
305
408
#define STRINGIFY(X) #X
306
409
#define FM_VERSION_BITS(BITS) "FM " STRINGIFY(BITS) " 1"
313
416
<< o.m_sampleSA << '\n';
314
417
315
418
out << o.m_alphabet.size() << '\n';
316
out.write((char*)o.m_alphabet.data(),
419
out.write(reinterpret_cast<const char*>(&o.m_alphabet[0]),
317
420
o.m_alphabet.size() * sizeof o.m_alphabet[0]);
318
421
319
422
out << o.m_sa.size() << '\n';
320
out.write((char*)o.m_sa.data(),
423
out.write(reinterpret_cast<const char*>(&o.m_sa[0]),
321
424
o.m_sa.size() * sizeof o.m_sa[0]);
322
425
323
426
return out << o.m_occ;
348
451
assert(c == '\n');
349
452
assert(n < std::numeric_limits<size_type>::max());
350
453
o.m_alphabet.resize(n);
351
in.read((char*)o.m_alphabet.data(), n * sizeof o.m_alphabet[0]);
454
in.read(reinterpret_cast<char*>(&o.m_alphabet[0]),
455
n * sizeof o.m_alphabet[0]);
352
456
o.setAlphabet(o.m_alphabet.begin(), o.m_alphabet.end());
353
457
354
458
in >> n;
357
461
assert(c == '\n');
358
462
assert(n < std::numeric_limits<size_type>::max());
359
463
o.m_sa.resize(n);
360
in.read((char*)o.m_sa.data(), n * sizeof o.m_sa[0]);
464
in.read(reinterpret_cast<char*>(&o.m_sa[0]),
465
n * sizeof o.m_sa[0]);
361
466
362
467
in >> o.m_occ;
363
468
assert(in);
Loggerhead is a web-based interface for Breezy
Version: 2.0.1