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