Loading...
Searching...
No Matches
simd_matrix_scoring_scheme.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 <concepts>
#include <seqan3/alignment/configuration/align_config_method.hpp>
#include <seqan3/alignment/scoring/scoring_scheme_concept.hpp>
#include <seqan3/alphabet/concept.hpp>
#include <seqan3/core/concept/cereal.hpp>
#include <seqan3/utility/simd/algorithm.hpp>
#include <seqan3/utility/simd/concept.hpp>
namespace seqan3::detail
{
template <simd_concept simd_score_t, semialphabet alphabet_t, typename alignment_t>
requires (std::same_as<alignment_t, align_cfg::method_local> || std::same_as<alignment_t, align_cfg::method_global>)
class simd_matrix_scoring_scheme
{
private:
using scalar_type = typename simd_traits<simd_score_t>::scalar_type;
using simd_score_profile_type = simd_score_t;
using simd_alphabet_ranks_type = simd_score_t;
static_assert(seqan3::alphabet_size<alphabet_t> <= std::numeric_limits<scalar_type>::max(),
"The selected simd scalar type is not large enough to represent the given alphabet including an "
"additional padding symbol!");
static_assert(seqan3::alphabet_size<alphabet_t> < std::numeric_limits<scalar_type>::max(),
"The selected simd scalar type is not large enough to represent the given alphabet including an "
"additional padding symbol!");
static constexpr bool is_global = std::same_as<alignment_t, align_cfg::method_global>;
static constexpr size_t index_offset = seqan3::alphabet_size<alphabet_t> + 1; // scheme is extended by one.
static constexpr scalar_type score_for_padding_symbol = (is_global) ? 1 : -1;
std::vector<scalar_type> scoring_scheme_data{};
public:
static constexpr scalar_type padding_symbol = static_cast<scalar_type>(seqan3::alphabet_size<alphabet_t>);
constexpr simd_matrix_scoring_scheme() = default;
constexpr simd_matrix_scoring_scheme(simd_matrix_scoring_scheme const &) = default;
constexpr simd_matrix_scoring_scheme(simd_matrix_scoring_scheme &&) = default;
constexpr simd_matrix_scoring_scheme & operator=(simd_matrix_scoring_scheme const &) = default;
constexpr simd_matrix_scoring_scheme & operator=(simd_matrix_scoring_scheme &&) = default;
~simd_matrix_scoring_scheme() = default;
template <typename scoring_scheme_t>
constexpr explicit simd_matrix_scoring_scheme(scoring_scheme_t const & scoring_scheme)
{
initialise_from_scalar_scoring_scheme(scoring_scheme);
}
template <typename scoring_scheme_t>
constexpr simd_matrix_scoring_scheme & operator=(scoring_scheme_t const & scoring_scheme)
{
initialise_from_scalar_scoring_scheme(scoring_scheme);
return *this;
}
constexpr simd_score_t score(simd_score_profile_type const & score_profile,
simd_alphabet_ranks_type const & ranks) const noexcept
{
simd_score_t const matrix_index = score_profile + ranks; // Compute the matrix indices for the lookup.
simd_score_t result{};
for (size_t idx = 0; idx < simd_traits<simd_score_t>::length; ++idx)
result[idx] = scoring_scheme_data.data()[matrix_index[idx]];
return result;
}
constexpr scalar_type padding_match_score() const noexcept
{
return score_for_padding_symbol;
}
constexpr simd_score_profile_type make_score_profile(simd_alphabet_ranks_type const & ranks) const noexcept
{
return ranks * simd::fill<simd_score_t>(index_offset);
}
private:
template <typename scoring_scheme_t>
requires scoring_scheme_for<scoring_scheme_t, alphabet_t>
constexpr void initialise_from_scalar_scoring_scheme(scoring_scheme_t const & scoring_scheme)
{
using score_t = decltype(std::declval<scoring_scheme_t const &>().score(alphabet_t{}, alphabet_t{}));
// Helper function to check if the matrix score is in the value range representable by the selected simd vector.
[[maybe_unused]] auto check_score_range = [&]([[maybe_unused]] score_t score)
{
// Note only if the size of the scalar type of the simd vector is smaller than the size of the score type
// of the original scoring scheme, the score might exceed the valid value range of the scalar type. In this
// case an exception will be thrown.
if constexpr (sizeof(scalar_type) < sizeof(score_t))
{
constexpr score_t max_score_value = static_cast<score_t>(std::numeric_limits<scalar_type>::max());
constexpr score_t min_score_value = static_cast<score_t>(std::numeric_limits<scalar_type>::lowest());
if (score > max_score_value || score < min_score_value)
throw std::invalid_argument{"The selected scoring scheme score overflows "
"for the selected scalar type of the simd type."};
}
};
// For the global alignment we extend the alphabet by one symbol to handle sequences with different size.
scoring_scheme_data.resize(index_offset * index_offset, score_for_padding_symbol);
// Convert the scoring matrix into a linear vector to allow gather operations later on.
using alphabet_size_t = std::remove_const_t<decltype(seqan3::alphabet_size<alphabet_t>)>;
auto data_it = scoring_scheme_data.begin();
for (alphabet_size_t lhs_rank = 0; lhs_rank < seqan3::alphabet_size<alphabet_t>; ++lhs_rank)
{
for (alphabet_size_t rhs_rank = 0; rhs_rank < seqan3::alphabet_size<alphabet_t>; ++rhs_rank, ++data_it)
{
score_t tmp_score = scoring_scheme.score(seqan3::assign_rank_to(lhs_rank, alphabet_t{}),
seqan3::assign_rank_to(rhs_rank, alphabet_t{}));
check_score_range(tmp_score);
*data_it = tmp_score;
}
++data_it; // skip one for the padded symbol.
}
}
};
} // namespace seqan3::detail
align_config_method.hpp
Provides global and local alignment configurations.
concept.hpp
Core alphabet concept and free function/type trait wrappers.
std::vector::begin
T begin(T... args)
cereal.hpp
Adaptions of concepts from the Cereal library.
seqan3::align_cfg::scoring_scheme::scoring_scheme
scoring_scheme(scheme_t) -> scoring_scheme< std::remove_cvref_t< scheme_t > >
Deduces the scoring scheme type from the constructor argument.
std::vector::data
T data(T... args)
std::declval
T declval(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::assign_rank_to
constexpr auto assign_rank_to
Assign a rank to an alphabet object.
Definition alphabet/concept.hpp:288
std::numeric_limits::lowest
T lowest(T... args)
std::numeric_limits::max
T max(T... args)
std::vector::resize
T resize(T... args)
scoring_scheme_concept.hpp
Provides seqan3::scoring_scheme_for.
Hide me