SeqAn3 3.4.0-rc.1
The Modern C++ library for sequence analysis.
Loading...
Searching...
No Matches
fm_index.hpp
Go to the documentation of this file.
1// -----------------------------------------------------------------------------------------------------
2// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin
3// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik
4// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
5// shipped with this file and also available at: https://github.com/seqan/seqan3/blob/master/LICENSE.md
6// -----------------------------------------------------------------------------------------------------
7
13#pragma once
14
15#include <algorithm>
16#include <filesystem>
17#include <ranges>
18
19#include <sdsl/suffix_trees.hpp>
20
26
27namespace seqan3::detail
28{
32struct fm_index_validator
33{
48 template <semialphabet alphabet_t, text_layout text_layout_mode_, std::ranges::range text_t>
49 static void validate(text_t && text)
50 {
51 if constexpr (text_layout_mode_ == text_layout::single)
52 {
53 static_assert(std::ranges::bidirectional_range<text_t>, "The text must model bidirectional_range.");
54 static_assert(std::convertible_to<range_innermost_value_t<text_t>, alphabet_t>,
55 "The alphabet of the text collection must be convertible to the alphabet of the index.");
56 static_assert(range_dimension_v<text_t> == 1, "The input cannot be a text collection.");
57
58 if (std::ranges::empty(text))
59 throw std::invalid_argument("The text to index cannot be empty.");
60 }
61 else
62 {
63 static_assert(std::ranges::bidirectional_range<text_t>,
64 "The text collection must model bidirectional_range.");
65 static_assert(std::ranges::bidirectional_range<std::ranges::range_reference_t<text_t>>,
66 "The elements of the text collection must model bidirectional_range.");
67 static_assert(std::convertible_to<range_innermost_value_t<text_t>, alphabet_t>,
68 "The alphabet of the text collection must be convertible to the alphabet of the index.");
69 static_assert(range_dimension_v<text_t> == 2, "The input must be a text collection.");
70
71 if (std::ranges::empty(text))
72 throw std::invalid_argument("The text collection to index cannot be empty.");
73 }
74 static_assert(alphabet_size<range_innermost_value_t<text_t>> <= 256, "The alphabet is too big.");
75 }
76};
77} // namespace seqan3::detail
78
79namespace seqan3
80{
81
83// forward declarations
84template <typename index_t>
85class fm_index_cursor;
86
87template <typename index_t>
88class bi_fm_index_cursor;
89
90namespace detail
91{
92template <semialphabet alphabet_t, text_layout text_layout_mode_, detail::sdsl_index sdsl_index_type_>
93class reverse_fm_index;
94}
96
130 sdsl::csa_wt<sdsl::wt_blcd<sdsl::bit_vector, // Wavelet tree type
131 sdsl::rank_support_v<>,
132 sdsl::select_support_scan<>,
133 sdsl::select_support_scan<0>>,
134 16, // Sampling rate of the suffix array
135 10'000'000, // Sampling rate of the inverse suffix array
136 sdsl::sa_order_sa_sampling<>, // How to sample positions in the suffix array (text VS SA sampling)
137 sdsl::isa_sampling<>, // How to sample positons in the inverse suffix array
138 sdsl::plain_byte_alphabet>; // How to represent the alphabet
139
146
185template <semialphabet alphabet_t,
186 text_layout text_layout_mode_,
187 detail::sdsl_index sdsl_index_type_ = default_sdsl_index_type>
189{
190private:
195 using sdsl_index_type = sdsl_index_type_;
199 using sdsl_char_type = typename sdsl_index_type::alphabet_type::char_type;
201 using sdsl_sigma_type = typename sdsl_index_type::alphabet_type::sigma_type;
203
204 friend class detail::reverse_fm_index<alphabet_t, text_layout_mode_, sdsl_index_type_>;
205
207 sdsl_index_type index;
208
210 sdsl::sd_vector<> text_begin;
212 sdsl::select_support_sd<1> text_begin_ss;
214 sdsl::rank_support_sd<1> text_begin_rs;
215
217 template <typename output_it_t, typename sequence_t>
218 static output_it_t copy_sequence_ranks_shifted_by_one(output_it_t output_it, sequence_t && sequence)
219 {
220 constexpr size_t sigma = alphabet_size<alphabet_t>;
221 constexpr size_t max_sigma = text_layout_mode_ == text_layout::single ? 256u : 255u;
222
223 constexpr auto warn_if_rank_out_of_range = [](uint8_t const rank)
224 {
225 if (rank >= max_sigma - 1) // same as rank + 1 >= max_sigma but without overflow
226 throw std::out_of_range("The input text cannot be indexed, because for full"
227 "character alphabets the last one/two values are reserved"
228 "(single sequence/collection).");
229 };
230
231 return std::ranges::transform(sequence,
232 output_it,
233 [&warn_if_rank_out_of_range](auto const & chr)
234 {
235 uint8_t const rank = seqan3::to_rank(chr);
236 if constexpr (sigma >= max_sigma)
237 warn_if_rank_out_of_range(rank);
238 else
239 (void)warn_if_rank_out_of_range;
240 return rank + 1;
241 })
242 .out;
243 }
244
264 template <std::ranges::range text_t>
265 requires (text_layout_mode_ == text_layout::single)
266 void construct(text_t && text)
267 {
268 detail::fm_index_validator::validate<alphabet_t, text_layout_mode_>(text);
269
270 // TODO:
271 // * check what happens in sdsl when constructed twice!
272 // * choose between in-memory/external and construction algorithms
273 // * sdsl construction currently only works for int_vector, std::string and char *, not ranges in general
274 // uint8_t largest_char = 0;
275 sdsl::int_vector<8> tmp_text(std::ranges::distance(text));
276
277 // copy ranks into tmp_text
278 copy_sequence_ranks_shifted_by_one(std::ranges::begin(tmp_text), text | std::views::reverse);
279
280 sdsl::construct_im(index, tmp_text, 0);
281
282 // TODO: would be nice but doesn't work since it's private and the public member references are const
283 // index.m_C.resize(largest_char);
284 // index.m_C.shrink_to_fit();
285 // index.m_sigma = largest_char;
286 }
287
289 template <std::ranges::range text_t>
290 requires (text_layout_mode_ == text_layout::collection)
291 void construct(text_t && text, bool reverse = false)
292 {
293 detail::fm_index_validator::validate<alphabet_t, text_layout_mode_>(text);
294
295 std::vector<size_t> text_sizes;
296
297 for (auto && t : text)
298 text_sizes.push_back(std::ranges::distance(t));
299
300 size_t const number_of_texts{text_sizes.size()};
301
302 // text size including delimiters
303 size_t const text_size = std::accumulate(text_sizes.begin(), text_sizes.end(), number_of_texts);
304
305 if (number_of_texts == text_size)
306 throw std::invalid_argument("A text collection that only contains empty texts cannot be indexed.");
307
308 constexpr auto sigma = alphabet_size<alphabet_t>;
309
310 // Instead of creating a bitvector of size `text_size`, setting the bits to 1 and then compressing it, we can
311 // use the `sd_vector_builder(text_size, number_of_ones)` because we know the parameters and the 1s we want to
312 // set are in a strictly increasing order. This inplace construction of the compressed vector saves memory.
313 sdsl::sd_vector_builder builder(text_size, number_of_texts);
314 size_t prefix_sum{0};
315
316 for (auto && size : text_sizes)
317 {
318 builder.set(prefix_sum);
319 prefix_sum += size + 1;
320 }
321
322 text_begin = sdsl::sd_vector<>(builder);
323 text_begin_ss = sdsl::select_support_sd<1>(&text_begin);
324 text_begin_rs = sdsl::rank_support_sd<1>(&text_begin);
325
326 // last text in collection needs no delimiter if we have more than one text in the collection
327 sdsl::int_vector<8> tmp_text(text_size - (number_of_texts > 1));
328
329 constexpr uint8_t delimiter = sigma >= 255 ? 255 : sigma + 1;
330
331 auto copy_join_with = [](auto output_it, auto && collection)
332 {
333 // this is basically std::views::join() with a delimiter
334 auto collection_it = std::ranges::begin(collection);
335 auto const collection_sentinel = std::ranges::end(collection);
336 if (collection_it == collection_sentinel)
337 return;
338
339 output_it = copy_sequence_ranks_shifted_by_one(output_it, *collection_it);
340 ++collection_it;
341
342 for (; collection_it != collection_sentinel; ++collection_it)
343 {
344 *output_it = delimiter;
345 ++output_it;
346 output_it = copy_sequence_ranks_shifted_by_one(output_it, *collection_it);
347 }
348 };
349
350 // copy ranks into tmp_text
351 copy_join_with(std::ranges::begin(tmp_text), text);
352
353 if (!reverse)
354 {
355 // we need at least one delimiter
356 if (number_of_texts == 1)
357 tmp_text.back() = delimiter;
358
359 std::ranges::reverse(tmp_text);
360 }
361 else
362 {
363 // If only one text is in the text collection, we still need one delimiter at the end to be able to
364 // conduct rank and select queries when locating hits in the index.
365 // Also, tmp_text looks like [text|0], but after reversing we need [txet|0] to be able to add the delimiter.
366 if (number_of_texts == 1)
367 {
368 std::ranges::reverse(tmp_text.begin(), tmp_text.end() - 1);
369 tmp_text.back() = delimiter;
370 }
371 else
372 {
373 std::ranges::reverse(tmp_text);
374 }
375 }
376
377 sdsl::construct_im(index, tmp_text, 0);
378 }
379
380public:
382 static constexpr text_layout text_layout_mode = text_layout_mode_;
383
388 using alphabet_type = alphabet_t;
390 using size_type = typename sdsl_index_type::size_type;
394
395 template <typename bi_fm_index_t>
396 friend class bi_fm_index_cursor;
397
398 template <typename fm_index_t>
399 friend class fm_index_cursor;
400
401 template <typename fm_index_t>
402 friend struct detail::fm_index_cursor_node;
403
407 fm_index() = default;
408
410 fm_index(fm_index const & rhs) :
411 index{rhs.index},
412 text_begin{rhs.text_begin},
413 text_begin_ss{rhs.text_begin_ss},
414 text_begin_rs{rhs.text_begin_rs}
415 {
416 text_begin_ss.set_vector(&text_begin);
417 text_begin_rs.set_vector(&text_begin);
418 }
419
422 index{std::move(rhs.index)},
423 text_begin{std::move(rhs.text_begin)},
424 text_begin_ss{std::move(rhs.text_begin_ss)},
425 text_begin_rs{std::move(rhs.text_begin_rs)}
426 {
427 text_begin_ss.set_vector(&text_begin);
428 text_begin_rs.set_vector(&text_begin);
429 }
430
433 {
434 index = std::move(rhs.index);
435 text_begin = std::move(rhs.text_begin);
436 text_begin_ss = std::move(rhs.text_begin_ss);
437 text_begin_rs = std::move(rhs.text_begin_rs);
438
439 text_begin_ss.set_vector(&text_begin);
440 text_begin_rs.set_vector(&text_begin);
441
442 return *this;
443 }
444
445 ~fm_index() = default;
446
455 template <std::ranges::bidirectional_range text_t>
456 explicit fm_index(text_t && text)
457 {
458 construct(std::forward<text_t>(text));
459 }
461
473 size_type size() const noexcept
474 {
475 return index.size();
476 }
477
489 bool empty() const noexcept
490 {
491 return size() == 0;
492 }
493
505 bool operator==(fm_index const & rhs) const noexcept
506 {
507 // (void) rhs;
508 return (index == rhs.index) && (text_begin == rhs.text_begin);
509 }
510
522 bool operator!=(fm_index const & rhs) const noexcept
523 {
524 return !(*this == rhs);
525 }
526
541 cursor_type cursor() const noexcept
542 {
543 return {*this};
544 }
545
553 template <cereal_archive archive_t>
554 void CEREAL_SERIALIZE_FUNCTION_NAME(archive_t & archive)
555 {
556 archive(index);
557 archive(text_begin);
558 archive(text_begin_ss);
559 text_begin_ss.set_vector(&text_begin);
560 archive(text_begin_rs);
561 text_begin_rs.set_vector(&text_begin);
562
563 auto sigma = alphabet_size<alphabet_t>;
564 archive(sigma);
565 if (sigma != alphabet_size<alphabet_t>)
566 {
567 throw std::logic_error{"The fm_index was built over an alphabet of size " + std::to_string(sigma)
568 + " but it is being read into an fm_index with an alphabet of size "
569 + std::to_string(alphabet_size<alphabet_t>) + "."};
570 }
571
572 bool tmp = text_layout_mode;
573 archive(tmp);
574 if (tmp != text_layout_mode)
575 {
576 throw std::logic_error{std::string{"The fm_index was built over a "}
577 + (tmp ? "text collection" : "single text")
578 + " but it is being read into an fm_index expecting a "
579 + (text_layout_mode ? "text collection." : "single text.")};
580 }
581 }
583};
584
589template <std::ranges::range text_t>
590fm_index(text_t &&) -> fm_index<range_innermost_value_t<text_t>, text_layout{range_dimension_v<text_t> != 1}>;
592} // namespace seqan3
593
594namespace seqan3::detail
595{
596
610template <semialphabet alphabet_t,
611 text_layout text_layout_mode,
612 detail::sdsl_index sdsl_index_type = default_sdsl_index_type>
613class reverse_fm_index : public fm_index<alphabet_t, text_layout_mode, sdsl_index_type>
614{
615private:
617 template <std::ranges::range text_t>
618 void construct_(text_t && text)
619 {
620 if constexpr (text_layout_mode == text_layout::single)
621 {
622 auto reverse_text = text | std::views::reverse;
623 this->construct(reverse_text);
624 }
625 else
626 {
627 auto reverse_text = text | views::deep{std::views::reverse} | std::views::reverse;
628 this->construct(reverse_text, true);
629 }
630 }
631
632public:
633 using fm_index<alphabet_t, text_layout_mode, sdsl_index_type>::fm_index;
634
636 template <std::ranges::bidirectional_range text_t>
637 explicit reverse_fm_index(text_t && text)
638 {
639 construct_(std::forward<text_t>(text));
640 }
641};
642
643} // namespace seqan3::detail
T accumulate(T... args)
T begin(T... args)
The SeqAn Bidirectional FM Index Cursor.
Definition bi_fm_index_cursor.hpp:55
The SeqAn FM Index Cursor.
Definition fm_index_cursor.hpp:87
The SeqAn FM Index.
Definition fm_index.hpp:189
fm_index & operator=(fm_index rhs)
When copy/move assigning, also update internal data structures.
Definition fm_index.hpp:432
fm_index()=default
Defaulted.
static constexpr text_layout text_layout_mode
Indicates whether index is built over a collection.
Definition fm_index.hpp:382
cursor_type cursor() const noexcept
Returns a seqan3::fm_index_cursor on the index that can be used for searching. .
Definition fm_index.hpp:541
size_type size() const noexcept
Returns the length of the indexed text including sentinel characters.
Definition fm_index.hpp:473
bool operator!=(fm_index const &rhs) const noexcept
Compares two indices.
Definition fm_index.hpp:522
bool empty() const noexcept
Checks whether the index is empty.
Definition fm_index.hpp:489
bool operator==(fm_index const &rhs) const noexcept
Compares two indices.
Definition fm_index.hpp:505
~fm_index()=default
Defaulted.
alphabet_t alphabet_type
The type of the underlying character of the indexed text.
Definition fm_index.hpp:388
fm_index(fm_index &&rhs)
When move constructing, also update internal data structures.
Definition fm_index.hpp:421
typename sdsl_index_type::size_type size_type
Type for representing positions in the indexed text.
Definition fm_index.hpp:390
fm_index(text_t &&text)
Constructor that immediately constructs the index given a range. The range cannot be empty.
Definition fm_index.hpp:456
fm_index(fm_index const &rhs)
When copy constructing, also update internal data structures.
Definition fm_index.hpp:410
Provides various transformation traits used by the range module.
Provides the internal representation of a node of the seqan3::fm_index_cursor.
T end(T... args)
Provides the seqan3::fm_index_cursor for searching in the unidirectional seqan3::fm_index.
constexpr auto to_rank
Return the rank representation of a (semi-)alphabet object.
Definition alphabet/concept.hpp:155
text_layout
The possible text layouts (single, collection) the seqan3::fm_index and seqan3::bi_fm_index can suppo...
Definition search/fm_index/concept.hpp:91
sdsl_wt_index_type default_sdsl_index_type
The default FM Index Configuration.
Definition fm_index.hpp:145
sdsl::csa_wt< sdsl::wt_blcd< sdsl::bit_vector, sdsl::rank_support_v<>, sdsl::select_support_scan<>, sdsl::select_support_scan< 0 > >, 16, 10 '000 '000, sdsl::sa_order_sa_sampling<>, sdsl::isa_sampling<>, sdsl::plain_byte_alphabet > sdsl_wt_index_type
The FM Index Configuration using a Wavelet Tree.
Definition fm_index.hpp:138
@ single
The text is a single range.
Definition search/fm_index/concept.hpp:93
@ collection
The text is a range of ranges.
Definition search/fm_index/concept.hpp:95
The basis for seqan3::alphabet, but requires only rank interface (not char).
The generic concept for a (biological) sequence.
The main SeqAn3 namespace.
Definition aligned_sequence_concept.hpp:29
fm_index(text_t &&) -> fm_index< range_innermost_value_t< text_t >, text_layout
Deduces the alphabet and dimensions of the text.
Definition fm_index.hpp:590
SeqAn specific customisations in the standard namespace.
T push_back(T... args)
Provides the concept for seqan3::detail::sdsl_index.
T size(T... args)
Provides seqan3::views::to_rank.
T to_string(T... args)