fm_index_cursor.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <array>
#include <ranges>
#include <type_traits>
 
#include <seqan3/alphabet/adaptation/char.hpp>
#include <seqan3/alphabet/concept.hpp>
#include <seqan3/contrib/sdsl-lite.hpp>
#include <seqan3/core/range/type_traits.hpp>
#include <seqan3/search/fm_index/concept.hpp>
#include <seqan3/search/fm_index/detail/fm_index_cursor.hpp>
#include <seqan3/utility/views/slice.hpp>
 
namespace seqan3
{
struct suffix_array_interval
{
    size_t begin_position{};
    size_t end_position{};
 
    friend bool operator==(suffix_array_interval const & lhs, suffix_array_interval const & rhs) noexcept
    {
        return lhs.begin_position == rhs.begin_position && lhs.end_position == rhs.end_position;
    }
 
    friend bool operator!=(suffix_array_interval const & lhs, suffix_array_interval const & rhs) noexcept
    {
        return !(lhs == rhs);
    }
};
 
template <typename index_t>
class fm_index_cursor
{
public:
    using index_type = index_t;
    using size_type = typename index_type::size_type;
 
private:
    using node_type = detail::fm_index_cursor_node<index_t>;
    using sdsl_char_type = typename index_type::sdsl_char_type;
    using sdsl_index_type = typename index_t::sdsl_index_type;
    using sdsl_sigma_type = typename index_type::sdsl_sigma_type;
    using index_alphabet_type = typename index_t::alphabet_type;
    using locate_result_value_type = std::pair<size_type, size_type>;
    using locate_result_type = std::vector<locate_result_value_type>;
 
    index_type const * index{nullptr};
    size_type parent_lb{};
    size_type parent_rb{};
    node_type node{};
    sdsl_sigma_type sigma{};
 
    template <typename _index_t>
    friend class bi_fm_index_cursor;
 
    size_type offset() const noexcept
    {
        assert(index->index.size() > query_length());
        return index->index.size() - query_length() - 1; // since the string is reversed during construction
    }
 
    bool
    backward_search(sdsl_index_type const & csa, sdsl_char_type const c, size_type & l, size_type & r) const noexcept
    {
        assert(l <= r && r < csa.size());
 
        size_type _l, _r;
 
        size_type cc = c;
        if constexpr (!std::same_as<index_alphabet_type, seqan3::contrib::sdsl::plain_byte_alphabet>)
        {
            cc = csa.char2comp[c];
            if (cc == 0 && c > 0) // [[unlikely]]
                return false;
        }
 
        size_type const c_begin = csa.C[cc];
        if (l == 0 && r + 1 == csa.size()) // [[unlikely]]
        {
            _l = c_begin;
            _r = csa.C[cc + 1] - 1;
            // if we use not the plain_byte_alphabet, we could return always return true here
        }
        else
        {
            _l = c_begin + csa.bwt.rank(l, c);         // count c in bwt[0..l-1]
            _r = c_begin + csa.bwt.rank(r + 1, c) - 1; // count c in bwt[0..r]
        }
 
        if (_r >= _l)
        {
            r = _r;
            l = _l;
            assert(r + 1 - l >= 0);
            return true;
        }
        return false;
    }
 
public:
    //        Default construction is necessary to make this class semi-regular and e.g., to allow construction of
    //        std::array of iterators.
    fm_index_cursor() noexcept = default;                                    
    fm_index_cursor(fm_index_cursor const &) noexcept = default;             
    fm_index_cursor & operator=(fm_index_cursor const &) noexcept = default; 
    fm_index_cursor(fm_index_cursor &&) noexcept = default;                  
    fm_index_cursor & operator=(fm_index_cursor &&) noexcept = default;      
    ~fm_index_cursor() = default;                                            
 
    fm_index_cursor(index_t const & _index) noexcept :
        index(&_index),
        node({0, _index.index.size() - 1, 0, 0}),
        sigma(_index.index.sigma - index_t::text_layout_mode)
    {
        assert(_index.index.size() != 0);
    }
    //\}
 
    bool operator==(fm_index_cursor const & rhs) const noexcept
    {
        assert(index != nullptr);
        assert(node != rhs.node || (query_length() == 0 || (parent_lb == rhs.parent_lb && parent_rb == rhs.parent_rb)));
 
        // position in the implicit suffix tree is defined by the SA interval and depth.
        // No need to compare parent intervals
        return node == rhs.node;
    }
 
    bool operator!=(fm_index_cursor const & rhs) const noexcept
    {
        assert(index != nullptr);
 
        return !(*this == rhs);
    }
 
    bool extend_right() noexcept
    {
        // TODO: specialize extend_right() and cycle_back() for EPR-dictionaries
        // store all cursors at once in a private std::array of cursors
        assert(index != nullptr);
 
        sdsl_sigma_type c = 1; // NOTE: start with 0 or 1 depending on implicit_sentintel
        size_type _lb = node.lb, _rb = node.rb;
        while (c < sigma && !backward_search(index->index, index->index.comp2char[c], _lb, _rb))
        {
            ++c;
        }
 
        if (c != sigma)
        {
            parent_lb = node.lb;
            parent_rb = node.rb;
            assert(c <= std::numeric_limits<sdsl_char_type>::max());
            node = {_lb, _rb, node.depth + 1, static_cast<sdsl_char_type>(c)};
            return true;
        }
        return false;
    }
 
    template <typename char_t>
        requires std::convertible_to<char_t, index_alphabet_type>
    bool extend_right(char_t const c) noexcept
    {
        assert(index != nullptr);
        // The rank cannot exceed 255 for single text and 254 for text collections as they are reserved as sentinels
        // for the indexed text.
        assert(seqan3::to_rank(static_cast<index_alphabet_type>(c))
               < ((index_type::text_layout_mode == text_layout::single) ? 255 : 254));
 
        size_type _lb = node.lb, _rb = node.rb;
 
        sdsl_char_type c_char = seqan3::to_rank(static_cast<index_alphabet_type>(c)) + 1;
 
        if (backward_search(index->index, c_char, _lb, _rb))
        {
            parent_lb = node.lb;
            parent_rb = node.rb;
            node = {_lb, _rb, node.depth + 1, c_char};
            return true;
        }
        return false;
    }
 
    template <typename char_type>
        requires detail::is_char_adaptation_v<char_type>
    bool extend_right(char_type const * cstring) noexcept
    {
        return extend_right(std::basic_string_view<char_type>{cstring});
    }
 
    template <std::ranges::range seq_t>
    bool extend_right(seq_t && seq) noexcept
    {
        static_assert(std::ranges::forward_range<seq_t>, "The query must model forward_range.");
        static_assert(std::convertible_to<range_innermost_value_t<seq_t>, index_alphabet_type>,
                      "The alphabet of the sequence must be convertible to the alphabet of the index.");
 
        assert(index != nullptr); // range must not be empty!
 
        size_type _lb = node.lb, _rb = node.rb;
        size_type new_parent_lb = parent_lb, new_parent_rb = parent_rb;
 
        sdsl_char_type c{};
        size_t len{0};
 
        for (auto it = std::ranges::begin(seq); it != std::ranges::end(seq); ++len, ++it)
        {
            // The rank cannot exceed 255 for single text and 254 for text collections as they are reserved as sentinels
            // for the indexed text.
            assert(seqan3::to_rank(static_cast<index_alphabet_type>(*it))
                   < ((index_type::text_layout_mode == text_layout::single) ? 255 : 254));
 
            c = seqan3::to_rank(static_cast<index_alphabet_type>(*it)) + 1;
 
            new_parent_lb = _lb;
            new_parent_rb = _rb;
            if (!backward_search(index->index, c, _lb, _rb))
                return false;
        }
 
        parent_lb = new_parent_lb;
        parent_rb = new_parent_rb;
        node = {_lb, _rb, len + node.depth, c};
        return true;
    }
 
    bool cycle_back() noexcept
    {
        assert(index != nullptr && query_length() > 0);
        // parent_lb > parent_rb --> invalid interval
        assert(parent_lb <= parent_rb);
 
        sdsl_sigma_type c = node.last_char + 1;
        size_type _lb = parent_lb, _rb = parent_rb;
 
        while (c < sigma && !backward_search(index->index, index->index.comp2char[c], _lb, _rb))
        {
            ++c;
        }
 
        if (c != sigma) // Collection has additional sentinel as delimiter
        {
            assert(c <= std::numeric_limits<sdsl_char_type>::max());
            node = {_lb, _rb, node.depth, static_cast<sdsl_char_type>(c)};
            return true;
        }
        return false;
    }
 
    size_type last_rank() const noexcept
    {
        // parent_lb > parent_rb --> invalid interval
        assert(index != nullptr && query_length() > 0 && parent_lb <= parent_rb);
 
        return index->index.comp2char[node.last_char] - 1; // text is not allowed to contain ranks of 0
    }
 
    seqan3::suffix_array_interval suffix_array_interval() const noexcept
    {
        assert(index != nullptr);
 
        return {node.lb, node.rb + 1};
    }
 
    size_type query_length() const noexcept
    {
        assert(index != nullptr);
        assert(node.depth != 0 || (node.lb == 0 && node.rb == index->size() - 1)); // depth == 0 -> root node
 
        return node.depth;
    }
 
    template <std::ranges::range text_t>
    auto path_label(text_t && text) const noexcept
        requires (index_t::text_layout_mode == text_layout::single)
    {
        static_assert(std::ranges::input_range<text_t>, "The text must model input_range.");
        static_assert(range_dimension_v<text_t> == 1, "The input cannot be a text collection.");
        static_assert(std::same_as<range_innermost_value_t<text_t>, index_alphabet_type>,
                      "The alphabet types of the given text and index differ.");
        assert(index != nullptr);
 
        size_type const query_begin = offset() - index->index[node.lb];
        return text | views::slice(query_begin, query_begin + query_length());
    }
 
    template <std::ranges::range text_t>
    auto path_label(text_t && text) const noexcept
        requires (index_t::text_layout_mode == text_layout::collection)
    {
        static_assert(std::ranges::input_range<text_t>, "The text collection must model input_range.");
        static_assert(range_dimension_v<text_t> == 2, "The input must be a text collection.");
        static_assert(std::same_as<range_innermost_value_t<text_t>, index_alphabet_type>,
                      "The alphabet types of the given text and index differ.");
        assert(index != nullptr);
 
        // Position of query in concatenated text.
        size_type const location = offset() - index->index[node.lb];
 
        // The rank represents the number of start positions of the individual sequences/texts in the collection
        // before position `location + 1` and thereby also the number of delimiters.
        size_type const rank = index->text_begin_rs.rank(location + 1);
        assert(rank > 0);
        size_type const text_id = rank - 1;
 
        // The start location of the `text_id`-th text in the sequence (position of the `rank`-th 1 in the bitvector).
        size_type const start_location = index->text_begin_ss.select(rank);
        // Substract lengths of previous sequences.
        size_type const query_begin = location - start_location;
 
        // Take subtext, slice query out of it
        return text[text_id] | views::slice(query_begin, query_begin + query_length());
    }
 
    size_type count() const noexcept
    {
        assert(index != nullptr);
 
        return 1 + node.rb - node.lb;
    }
 
    locate_result_type locate() const
        requires (index_t::text_layout_mode == text_layout::single)
    {
        assert(index != nullptr);
 
        locate_result_type occ{};
        occ.reserve(count());
        for (size_type i = 0; i < count(); ++i)
        {
            occ.emplace_back(0, offset() - index->index[node.lb + i]);
        }
 
        return occ;
    }
 
    locate_result_type locate() const
        requires (index_t::text_layout_mode == text_layout::collection)
    {
        assert(index != nullptr);
 
        locate_result_type occ;
        occ.reserve(count());
        for (size_type i = 0; i < count(); ++i)
        {
            size_type loc = offset() - index->index[node.lb + i];
            size_type sequence_rank = index->text_begin_rs.rank(loc + 1);
            size_type sequence_position = loc - index->text_begin_ss.select(sequence_rank);
            occ.emplace_back(sequence_rank - 1, sequence_position);
        }
        return occ;
    }
 
    auto lazy_locate() const
        requires (index_t::text_layout_mode == text_layout::single)
    {
        assert(index != nullptr);
 
        return std::views::iota(node.lb, node.lb + count())
             | std::views::transform(
                   [*this, _offset = offset()](auto sa_pos)
                   {
                       return locate_result_value_type{0u, _offset - index->index[sa_pos]};
                   });
    }
 
    auto lazy_locate() const
        requires (index_t::text_layout_mode == text_layout::collection)
    {
        assert(index != nullptr);
 
        return std::views::iota(node.lb, node.lb + count())
             | std::views::transform(
                   [*this, _offset = offset()](auto sa_pos)
                   {
                       auto loc = _offset - index->index[sa_pos];
                       size_type sequence_rank = index->text_begin_rs.rank(loc + 1);
                       size_type sequence_position = loc - index->text_begin_ss.select(sequence_rank);
                       return locate_result_value_type{sequence_rank - 1, sequence_position};
                   });
    }
 
    template <cereal_archive archive_t>
    void CEREAL_SERIALIZE_FUNCTION_NAME(archive_t & archive)
    {
        archive(parent_lb);
        archive(parent_rb);
        archive(node);
        archive(sigma);
    }
};
 
} // namespace seqan3