fm_index.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2019, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2019, Knut Reinert & MPI für molekulare Genetik
// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
// shipped with this file and also available at: https://github.com/seqan/seqan3/blob/master/LICENSE.md
// -----------------------------------------------------------------------------------------------------

#pragma once

#include <sdsl/suffix_trees.hpp>

#include <seqan3/core/type_traits/range.hpp>
#include <seqan3/std/filesystem>
#include <seqan3/range/shortcuts.hpp>
#include <seqan3/range/view/to_rank.hpp>
#include <seqan3/search/fm_index/concept.hpp>
#include <seqan3/search/fm_index/detail/csa_alphabet_strategy.hpp>
#include <seqan3/search/fm_index/detail/fm_index_cursor.hpp>
#include <seqan3/search/fm_index/fm_index_cursor.hpp>
#include <seqan3/std/algorithm>
#include <seqan3/std/ranges>

namespace seqan3
{

// forward declarations
template <typename index_t>
class fm_index_cursor;

template <typename index_t>
class bi_fm_index_cursor;

using sdsl_wt_index_type =
    sdsl::csa_wt<sdsl::wt_blcd<sdsl::bit_vector,
                               sdsl::rank_support_v<>,
                               sdsl::select_support_scan<>,
                               sdsl::select_support_scan<0>>,
                 16,
                 10000000,
                 sdsl::sa_order_sa_sampling<>,
                 sdsl::isa_sampling<>,
                 sdsl::plain_byte_alphabet>;

using default_sdsl_index_type = sdsl_wt_index_type;

template <bool is_collection = false, detail::SdslIndex sdsl_index_type_ = default_sdsl_index_type>
class fm_index
{
protected:
    size_t sigma{0};
    static constexpr bool is_collection_{is_collection};

    using sdsl_index_type = sdsl_index_type_;
    using sdsl_char_type = typename sdsl_index_type::alphabet_type::char_type;
    using sdsl_sigma_type = typename sdsl_index_type::alphabet_type::sigma_type;

    sdsl_index_type index;

    sdsl::sd_vector<> text_begin;
    sdsl::select_support_sd<1> text_begin_ss;
    sdsl::rank_support_sd<1> text_begin_rs;

public:
    using size_type = typename sdsl_index_type::size_type;
    using cursor_type = fm_index_cursor<fm_index<is_collection_, sdsl_index_type>>;

    template <typename bi_fm_index_t>
    friend class bi_fm_index_cursor;

    template <typename fm_index_t>
    friend class fm_index_cursor;

    template <typename fm_index_t>
    friend class detail::fm_index_cursor_node;

    fm_index() = default;                             
    fm_index(fm_index const &) = default;             
    fm_index & operator=(fm_index const &) = default; 
    fm_index(fm_index &&) = default;                  
    fm_index & operator=(fm_index &&) = default;      
    ~fm_index() = default;                            

    template <std::ranges::Range text_t>
    fm_index(text_t && text)
    {
        construct(std::forward<text_t>(text));
    }

    template <std::ranges::Range text_t>
    void construct(text_t && text)
        requires !is_collection_
    {
        static_assert(std::ranges::BidirectionalRange<text_t>, "The text must model BidirectionalRange.");
        static_assert(alphabet_size<innermost_value_type_t<text_t>> <= 256, "The alphabet is too big.");
        static_assert(dimension_v<text_t> == 1, "The input cannot be a text collection.");

        // text must not be empty
        if (std::ranges::begin(text) == std::ranges::end(text))
            throw std::invalid_argument("The text that is indexed cannot be empty.");

        constexpr auto cexpr_sigma = alphabet_size<innermost_value_type_t<text_t>>;
        sigma = cexpr_sigma;
        // TODO:
        // * check what happens in sdsl when constructed twice!
        // * choose between in-memory/external and construction algorithms
        // * sdsl construction currently only works for int_vector, std::string and char *, not ranges in general
        // uint8_t largest_char = 0;
        sdsl::int_vector<8> tmp_text(text.size());

        std::ranges::copy(text
                          | view::to_rank
                          | std::view::transform([] (uint8_t const r)
                          {
                              if constexpr (cexpr_sigma == 256)
                              {
                                  if (r == 255)
                                      throw std::out_of_range("The input text cannot be indexed, because for full"
                                                              "character alphabets the last one/two values are reserved"
                                                              "(single sequence/collection).");
                              }
                              return r + 1;
                          })
                          | std::view::reverse,
                          seqan3::begin(tmp_text)); // reverse and increase rank by one

        sdsl::construct_im(index, tmp_text, 0);

        // TODO: would be nice but doesn't work since it's private and the public member references are const
        // index.m_C.resize(largest_char);
        // index.m_C.shrink_to_fit();
        // index.m_sigma = largest_char;
    }

    template <std::ranges::Range text_t>
    void construct(text_t && text)
        requires is_collection_
    {
        static_assert(std::ranges::BidirectionalRange<text_t>, "The text collection must model BidirectionalRange.");
        static_assert(std::ranges::BidirectionalRange<reference_t<text_t>>,
                      "The elements of the text collection must model BidirectionalRange.");
        static_assert(alphabet_size<innermost_value_type_t<text_t>> <= 256, "The alphabet is too big.");
        static_assert(dimension_v<text_t> == 2, "The input must be a text collection.");

        // text collection must not be empty
        if (std::ranges::begin(text) == std::ranges::end(text))
            throw std::invalid_argument("The text that is indexed cannot be empty.");

        size_t text_size{0}; // text size including delimiters

        // there must be at least one non-empty text in the collection
        bool all_empty = true;

        for (auto && t : text)
        {
            if (std::ranges::begin(t) != std::ranges::end(t))
            {
                all_empty = false;
            }
            text_size += 1 + t.size(); // text size and delimiter (sum will be 1 for empty texts)
        }

        if (all_empty)
            throw std::invalid_argument("A text collection that only contains empty texts cannot be indexed.");

        constexpr auto cexpr_sigma = alphabet_size<innermost_value_type_t<text_t>>;
        sigma = cexpr_sigma;

        // bitvector where 1 marks the begin position of a single text from the collection in the concatenated text
        sdsl::bit_vector pos(text_size, 0);
        size_t prefix_sum{0};

        for (auto && t : text)
        {
            pos[prefix_sum] = 1;
            prefix_sum += t.size() + 1;
        }

        text_begin    = sdsl::sd_vector(pos);
        text_begin_ss = sdsl::select_support_sd<1>(&text_begin);
        text_begin_rs = sdsl::rank_support_sd<1>(&text_begin);

        sdsl::int_vector<8> tmp_text(text_size - 1); // last text in collection needs no delimiter

        constexpr uint8_t delimiter = cexpr_sigma >= 255 ? 255 : cexpr_sigma + 1;

        std::vector<uint8_t> tmp = text
                                   | view::deep{view::to_rank}
                                   | view::deep
                                   {
                                       std::view::transform([] (uint8_t const r)
                                       {
                                           if constexpr (cexpr_sigma >= 255)
                                           {
                                               if (r >= 254)
                                                   throw std::out_of_range("The input text cannot be indexed, because"
                                                                           " for full character alphabets the last one/"
                                                                           "two values are reserved (single sequence/"
                                                                           "collection).");
                                           }
                                           return r + 1;
                                       })
                                   }
                                   | std::view::join(delimiter);

        std::ranges::copy((tmp | std::view::reverse), seqan3::begin(tmp_text));

        // std::ranges::copy(text
        //                   | view::deep{view::to_rank}
        //                   | view::deep{std::view::transform([] (uint8_t const r) { return r + 1; })} // increase rank
        //                   | view::deep{std::view::reverse}
        //                   | std::view::reverse
        //                   | std::view::join(delimiter), // join with delimiter
        //                   seqan3::begin(tmp_text));

        sdsl::construct_im(index, tmp_text, 0);
    }

    size_type size() const noexcept
    {
        return index.size();
    }

    bool empty() const noexcept
    {
        return size() == 0;
    }

    bool operator==(fm_index const & rhs) const noexcept
    {
        // (void) rhs;
        return (index == rhs.index) && (text_begin == rhs.text_begin);
    }

    bool operator!=(fm_index const & rhs) const noexcept
    {
        return !(*this == rhs);
    }

    cursor_type begin() const noexcept
    {
        return {*this};
    }

    template <CerealArchive archive_t>
    void CEREAL_SERIALIZE_FUNCTION_NAME(archive_t & archive)
    {
        archive(index);
        archive(text_begin);
        archive(text_begin_ss);
        text_begin_ss.set_vector(&text_begin);
        archive(text_begin_rs);
        text_begin_rs.set_vector(&text_begin);
        archive(sigma);
        bool tmp = is_collection_;
        archive(tmp);
        assert(tmp == is_collection_);
    }

};

template <std::ranges::Range text_t>
fm_index(text_t &&) -> fm_index<dimension_v<text_t> != 1>;

} // namespace seqan3