Loading...
Searching...
No Matches
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
concept.hpp
Core alphabet concept and free function/type trait wrappers.
std::forward
T forward(T... args)
seqan3::alphabet_size
constexpr auto alphabet_size
A type trait that holds the size of a (semi-)alphabet.
Definition alphabet/concept.hpp:834
seqan3::to_rank
constexpr auto to_rank
Return the rank representation of a (semi-)alphabet object.
Definition alphabet/concept.hpp:152
seqan3::views::kmer_hash
constexpr auto kmer_hash
Computes hash values for each position of a range via a given shape.
Definition kmer_hash.hpp:766
seqan3::pow
base_t pow(base_t base, exp_t exp)
Computes the value of base raised to the power exp.
Definition math.hpp:120
math.hpp
Provides math related functionality.
std::max
T max(T... args)
seqan3::views
The SeqAn namespace for views.
Definition char_strictly_to.hpp:19
std::next
T next(T... args)
std::rel_ops::operator!=
T operator!=(T... args)
hash.hpp
Provides overloads for std::hash.
shape.hpp
Provides seqan3::shape.
std::tie
T tie(T... args)
std::to_string
T to_string(T... args)
Hide me