kmer_hash.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 <seqan3/alphabet/concept.hpp>
#include <seqan3/alphabet/range/hash.hpp>
#include <seqan3/search/kmer_index/shape.hpp>
#include <seqan3/utility/math.hpp>
 
namespace seqan3::detail
{
// ---------------------------------------------------------------------------------------------------------------------
// kmer_hash_view class
// ---------------------------------------------------------------------------------------------------------------------
 
template <std::ranges::view urng_t>
class kmer_hash_view : public std::ranges::view_interface<kmer_hash_view<urng_t>>
{
private:
    static_assert(std::ranges::forward_range<urng_t>, "The kmer_hash_view only works on forward_ranges");
    static_assert(semialphabet<std::ranges::range_reference_t<urng_t>>,
                  "The reference type of the underlying range must model seqan3::semialphabet.");
 
    urng_t urange;
 
    shape shape_;
 
    template <bool const_range>
    class basic_iterator;
 
    static inline int max_shape_count = 64 / std::log2(alphabet_size<std::ranges::range_reference_t<urng_t>>);
 
    [[noreturn]] void shape_too_long_error() const
    {
        std::string message{"The shape is too long for the given alphabet.\n"};
        message += "Alphabet: ";
        // Note: Since we want the alphabet type name, std::ranges::range_value_t is the better choice.
        // For seqan3::bitpacked_sequence<seqan3::dna4>:
        // reference_t: seqan3::bitpacked_sequence<seqan3::dna4>::reference_proxy_type
        // value_t: seqan3::dna4
        message += detail::type_name_as_string<std::remove_cvref_t<std::ranges::range_value_t<urng_t>>>;
        message += "\nMaximum shape count: ";
        message += std::to_string(max_shape_count);
        message += "\nGiven shape count: ";
        message += std::to_string(shape_.count());
        throw std::invalid_argument{message};
    }
 
    inline void validate_shape() const
    {
        if (shape_.count() > max_shape_count)
            shape_too_long_error();
    }
 
public:
    kmer_hash_view()
        requires std::default_initializable<urng_t>
    = default;                                                        
    kmer_hash_view(kmer_hash_view const & rhs) = default;             
    kmer_hash_view(kmer_hash_view && rhs) = default;                  
    kmer_hash_view & operator=(kmer_hash_view const & rhs) = default; 
    kmer_hash_view & operator=(kmer_hash_view && rhs) = default;      
    ~kmer_hash_view() = default;                                      
 
    explicit kmer_hash_view(urng_t urange_, shape const & s_) : urange{std::move(urange_)}, shape_{s_}
    {
        validate_shape();
    }
 
    template <typename rng_t>
        requires (!std::same_as<std::remove_cvref_t<rng_t>, kmer_hash_view>) && std::ranges::viewable_range<rng_t>
                  && std::constructible_from<urng_t, std::ranges::ref_view<std::remove_reference_t<rng_t>>>
    explicit kmer_hash_view(rng_t && urange_, shape const & s_) :
        urange{std::views::all(std::forward<rng_t>(urange_))},
        shape_{s_}
    {
        validate_shape();
    }
 
    auto begin() noexcept
    {
        return basic_iterator<false>{std::ranges::begin(urange), std::ranges::end(urange), shape_};
    }
 
    auto begin() const noexcept
        requires const_iterable_range<urng_t>
    {
        return basic_iterator<true>{std::ranges::begin(urange), std::ranges::end(urange), shape_};
    }
 
    auto end() noexcept
    {
        // Assigning the end iterator to the text_right iterator of the basic_iterator only works for common ranges.
        if constexpr (std::ranges::common_range<urng_t>)
            return basic_iterator<false>{std::ranges::begin(urange), std::ranges::end(urange), shape_, true};
        else
            return std::ranges::end(urange);
    }
 
    auto end() const noexcept
        requires const_iterable_range<urng_t>
    {
        // Assigning the end iterator to the text_right iterator of the basic_iterator only works for common ranges.
        if constexpr (std::ranges::common_range<urng_t const>)
            return basic_iterator<true>{std::ranges::begin(urange), std::ranges::end(urange), shape_, true};
        else
            return std::ranges::cend(urange);
    }
 
    auto size()
        requires std::ranges::sized_range<urng_t>
    {
        using size_type = std::ranges::range_size_t<urng_t>;
        return std::max<size_type>(std::ranges::size(urange) + 1, shape_.size()) - shape_.size();
    }
 
    auto size() const
        requires std::ranges::sized_range<urng_t const>
    {
        using size_type = std::ranges::range_size_t<urng_t const>;
        return std::max<size_type>(std::ranges::size(urange) + 1, shape_.size()) - shape_.size();
    }
};
 
template <std::ranges::view urng_t>
template <bool const_range>
class kmer_hash_view<urng_t>::basic_iterator :
    public maybe_iterator_category<maybe_const_iterator_t<const_range, urng_t>>
{
private:
    using it_t = maybe_const_iterator_t<const_range, urng_t>;
    using sentinel_t = maybe_const_sentinel_t<const_range, urng_t>;
 
    template <bool other_const_range>
    friend class basic_iterator;
 
public:
    using difference_type = typename std::iter_difference_t<it_t>;
    using value_type = size_t;
    using pointer = void;
    using reference = value_type;
    using iterator_concept = std::conditional_t<std::contiguous_iterator<it_t>,
                                                typename std::random_access_iterator_tag,
                                                detail::iterator_concept_tag_t<it_t>>;
 
    constexpr basic_iterator() = default;                                   
    constexpr basic_iterator(basic_iterator const &) = default;             
    constexpr basic_iterator(basic_iterator &&) = default;                  
    constexpr basic_iterator & operator=(basic_iterator const &) = default; 
    constexpr basic_iterator & operator=(basic_iterator &&) = default;      
    ~basic_iterator() = default;                                            
 
    constexpr basic_iterator(basic_iterator<!const_range> const & it) noexcept
        requires const_range
        :
        hash_value{std::move(it.hash_value)},
        roll_factor{std::move(it.roll_factor)},
        shape_{std::move(it.shape_)},
        text_left{std::move(it.text_left)},
        text_right{std::move(it.text_right)}
    {}
 
    basic_iterator(it_t it_start, sentinel_t it_end, shape s_) :
        shape_{s_},
        text_left{it_start},
        text_right{std::ranges::next(text_left, shape_.size() - 1, it_end)}
    {
        assert(std::ranges::size(shape_) > 0);
 
        // shape size = 3
        // Text:      1 2 3 4 5 6 7 8 9
        // text_left: ^
        // text_right:    ^
        // distance(text_left, text_right) = 2
        if (shape_.size() <= std::ranges::distance(text_left, text_right) + 1)
        {
            roll_factor = pow(sigma, static_cast<size_t>(std::ranges::size(shape_) - 1));
            hash_full();
        }
    }
 
    basic_iterator(it_t it_start, sentinel_t it_end, shape s_, bool SEQAN3_DOXYGEN_ONLY(is_end)) : shape_{s_}
    {
        assert(std::ranges::size(shape_) > 0);
 
        auto urange_size = std::ranges::distance(it_start, it_end);
        auto step = (shape_.size() > urange_size + 1) ? 0 : urange_size - shape_.size() + 1;
        text_left = std::ranges::next(it_start, step, it_end);
 
        // shape size = 3
        // Text:      1 2 3 4 5 6 7 8 9
        // text_left: ^
        // text_right:    ^
        // distance(text_left, text_right) = 2
        if (shape_.size() <= std::ranges::distance(text_left, it_end) + 1)
        {
            roll_factor = pow(sigma, static_cast<size_t>(std::ranges::size(shape_) - 1));
            hash_full();
        }
 
        text_right = it_end;
    }
 
 
    friend bool operator==(basic_iterator const & lhs, sentinel_t const & rhs) noexcept
    {
        return lhs.text_right == rhs;
    }
 
    friend bool operator==(sentinel_t const & lhs, basic_iterator const & rhs) noexcept
    {
        return lhs == rhs.text_right;
    }
 
    friend bool operator==(basic_iterator const & lhs, basic_iterator const & rhs) noexcept
    {
        return std::tie(lhs.text_right, lhs.shape_) == std::tie(rhs.text_right, rhs.shape_);
    }
 
    friend bool operator!=(basic_iterator const & lhs, sentinel_t const & rhs) noexcept
    {
        return !(lhs == rhs);
    }
 
    friend bool operator!=(sentinel_t const & lhs, basic_iterator const & rhs) noexcept
    {
        return !(lhs == rhs);
    }
 
    friend bool operator!=(basic_iterator const & lhs, basic_iterator const & rhs) noexcept
    {
        return !(lhs == rhs);
    }
 
    friend bool operator<(basic_iterator const & lhs, basic_iterator const & rhs) noexcept
    {
        return (lhs.shape_ <= rhs.shape_) && (lhs.text_right < rhs.text_right);
    }
 
    friend bool operator>(basic_iterator const & lhs, basic_iterator const & rhs) noexcept
    {
        return (lhs.shape_ >= rhs.shape_) && (lhs.text_right > rhs.text_right);
    }
 
    friend bool operator<=(basic_iterator const & lhs, basic_iterator const & rhs) noexcept
    {
        return (lhs.shape_ <= rhs.shape_) && (lhs.text_right <= rhs.text_right);
    }
 
    friend bool operator>=(basic_iterator const & lhs, basic_iterator const & rhs) noexcept
    {
        return (lhs.shape_ >= rhs.shape_) && (lhs.text_right >= rhs.text_right);
    }
 
 
    basic_iterator & operator++() noexcept
    {
        hash_forward();
        return *this;
    }
 
    basic_iterator operator++(int) noexcept
    {
        basic_iterator tmp{*this};
        hash_forward();
        return tmp;
    }
 
    basic_iterator & operator--() noexcept
        requires std::bidirectional_iterator<it_t>
    {
        hash_backward();
        return *this;
    }
 
    basic_iterator operator--(int) noexcept
        requires std::bidirectional_iterator<it_t>
    {
        basic_iterator tmp{*this};
        hash_backward();
        return tmp;
    }
 
    basic_iterator & operator+=(difference_type const skip) noexcept
        requires std::random_access_iterator<it_t>
    {
        hash_forward(skip);
        return *this;
    }
 
    basic_iterator operator+(difference_type const skip) const noexcept
        requires std::random_access_iterator<it_t>
    {
        basic_iterator tmp{*this};
        return tmp += skip;
    }
 
    friend basic_iterator operator+(difference_type const skip, basic_iterator const & it) noexcept
        requires std::random_access_iterator<it_t>
    {
        return it + skip;
    }
 
    basic_iterator & operator-=(difference_type const skip) noexcept
        requires std::random_access_iterator<it_t>
    {
        hash_backward(skip);
        return *this;
    }
 
    basic_iterator operator-(difference_type const skip) const noexcept
        requires std::random_access_iterator<it_t>
    {
        basic_iterator tmp{*this};
        return tmp -= skip;
    }
 
    friend basic_iterator operator-(difference_type const skip, basic_iterator const & it) noexcept
        requires std::random_access_iterator<it_t>
    {
        return it - skip;
    }
 
    friend difference_type operator-(basic_iterator const & lhs, basic_iterator const & rhs) noexcept
        requires std::sized_sentinel_for<it_t, it_t>
    {
        return static_cast<difference_type>(lhs.text_right - rhs.text_right);
    }
 
    friend difference_type operator-(sentinel_t const & lhs, basic_iterator const & rhs) noexcept
        requires std::sized_sentinel_for<sentinel_t, it_t>
    {
        return static_cast<difference_type>(lhs - rhs.text_right);
    }
 
    friend difference_type operator-(basic_iterator const & lhs, sentinel_t const & rhs) noexcept
        requires std::sized_sentinel_for<it_t, sentinel_t>
    {
        return static_cast<difference_type>(lhs.text_right - rhs);
    }
 
    reference operator[](difference_type const n) const
        requires std::random_access_iterator<it_t>
    {
        return *(*this + n);
    }
 
    value_type operator*() const noexcept
    {
        return hash_value + to_rank(*text_right);
    }
 
private:
    using alphabet_t = std::iter_value_t<it_t>;
 
    static constexpr auto const sigma{alphabet_size<alphabet_t>};
 
    size_t hash_value{0};
 
    size_t roll_factor{0};
 
    shape shape_;
 
    it_t text_left;
 
    it_t text_right;
 
    void hash_forward()
    {
        if (shape_.all())
        {
            hash_roll_forward();
        }
        else
        {
            std::ranges::advance(text_left, 1);
            hash_full();
        }
    }
 
    void hash_forward(difference_type const skip)
        requires std::random_access_iterator<it_t>
    {
        std::ranges::advance(text_left, skip);
        hash_full();
    }
 
    void hash_backward()
        requires std::bidirectional_iterator<it_t>
    {
        if (shape_.all())
        {
            hash_roll_backward();
        }
        else
        {
            std::ranges::advance(text_left, -1);
            hash_full();
        }
    }
 
    void hash_backward(difference_type const skip)
    {
        std::ranges::advance(text_left, -skip);
        hash_full();
    }
 
    void hash_full()
    {
        text_right = text_left;
        hash_value = 0;
 
        for (size_t i{0}; i < shape_.size() - 1u; ++i)
        {
            hash_value += shape_[i] * to_rank(*text_right);
            hash_value *= shape_[i] ? sigma : 1;
            std::ranges::advance(text_right, 1);
        }
    }
 
    void hash_roll_forward()
    {
        hash_value -= to_rank(*(text_left)) * roll_factor;
        hash_value += to_rank(*(text_right));
        hash_value *= sigma;
 
        std::ranges::advance(text_left, 1);
        std::ranges::advance(text_right, 1);
    }
 
    void hash_roll_backward()
        requires std::bidirectional_iterator<it_t>
    {
        std::ranges::advance(text_left, -1);
        std::ranges::advance(text_right, -1);
 
        hash_value /= sigma;
        hash_value -= to_rank(*(text_right));
        hash_value += to_rank(*(text_left)) * roll_factor;
    }
};
 
template <std::ranges::viewable_range rng_t>
kmer_hash_view(rng_t &&, shape const & shape_) -> kmer_hash_view<std::views::all_t<rng_t>>;
 
// ---------------------------------------------------------------------------------------------------------------------
// kmer_hash_fn (adaptor definition)
// ---------------------------------------------------------------------------------------------------------------------
 
struct kmer_hash_fn
{
    constexpr auto operator()(shape const & shape_) const
    {
        return adaptor_from_functor{*this, shape_};
    }
 
    template <std::ranges::range urng_t>
    constexpr auto operator()(urng_t && urange, shape const & shape_) const
    {
        static_assert(std::ranges::viewable_range<urng_t>,
                      "The range parameter to views::kmer_hash cannot be a temporary of a non-view range.");
        static_assert(std::ranges::forward_range<urng_t>,
                      "The range parameter to views::kmer_hash must model std::ranges::forward_range.");
        static_assert(semialphabet<std::ranges::range_reference_t<urng_t>>,
                      "The range parameter to views::kmer_hash must be over elements of seqan3::semialphabet.");
 
        return kmer_hash_view{std::forward<urng_t>(urange), shape_};
    }
};
 
} // namespace seqan3::detail
 
namespace seqan3::views
{
inline constexpr auto kmer_hash = detail::kmer_hash_fn{};
 
} // namespace seqan3::views