simd_match_mismatch_scoring_scheme.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2021, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2021, Knut Reinert & MPI für molekulare Genetik
// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
// shipped with this file and also available at: https://github.com/seqan/seqan3/blob/master/LICENSE.md
// -----------------------------------------------------------------------------------------------------
 
#pragma once
 
#include <seqan3/std/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/detail/bits_of.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 (seqan3::alphabet_size<alphabet_t> > 1) &&
             (std::same_as<alignment_t, align_cfg::method_local> ||
              std::same_as<alignment_t, align_cfg::method_global>)
class simd_match_mismatch_scoring_scheme
{
private:
    using scalar_type = typename simd_traits<simd_score_t>::scalar_type;
    using alphabet_ranks_type = simd_score_t;
 
public:
    static constexpr scalar_type padding_symbol = static_cast<scalar_type>(1u << (bits_of<scalar_type> - 1));
 
    constexpr simd_match_mismatch_scoring_scheme() = default;
    constexpr simd_match_mismatch_scoring_scheme(simd_match_mismatch_scoring_scheme const &) = default;
    constexpr simd_match_mismatch_scoring_scheme(simd_match_mismatch_scoring_scheme &&) = default;
    constexpr simd_match_mismatch_scoring_scheme & operator=(simd_match_mismatch_scoring_scheme const &) = default;
    constexpr simd_match_mismatch_scoring_scheme & operator=(simd_match_mismatch_scoring_scheme &&) = default;
    ~simd_match_mismatch_scoring_scheme() = default;
 
    template <typename scoring_scheme_t>
        requires scoring_scheme_for<scoring_scheme_t, alphabet_t>
    constexpr explicit simd_match_mismatch_scoring_scheme(scoring_scheme_t const & scoring_scheme)
    {
        initialise_from_scalar_scoring_scheme(scoring_scheme);
    }
 
    template <typename scoring_scheme_t>
        requires scoring_scheme_for<scoring_scheme_t, alphabet_t>
    constexpr simd_match_mismatch_scoring_scheme & operator=(scoring_scheme_t const & scoring_scheme)
    {
        initialise_from_scalar_scoring_scheme(scoring_scheme);
        return *this;
    }
 
    constexpr simd_score_t score(alphabet_ranks_type const & ranks1, alphabet_ranks_type const & ranks2)
        const noexcept
    {
        typename simd_traits<simd_score_t>::mask_type mask;
        // For global and local alignment there are slightly different formulas because
        // in global alignment padded characters always match
        if constexpr (std::same_as<alignment_t, align_cfg::method_global>)
            mask = (ranks1 ^ ranks2) <= simd::fill<simd_score_t>(0);
        else // and in local alignment type padded characters always mismatch.
            mask = (ranks1 ^ ranks2) == simd::fill<simd_score_t>(0);
 
        return mask ? match_score : mismatch_score;
    }
 
    constexpr auto padding_match_score() noexcept
    {
        return match_score[0];
    }
 
    template <typename alphabet_ranks_t>
        requires std::same_as<std::remove_cvref_t<alphabet_ranks_t>, alphabet_ranks_type>
    constexpr alphabet_ranks_t make_score_profile(alphabet_ranks_t && ranks) const noexcept
    {
        return std::forward<alphabet_ranks_t>(ranks);
    }
 
private:
    template <typename scoring_scheme_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{}));
        using simd_scalar_t = typename simd_traits<simd_score_t>::scalar_type;
 
        score_t scalar_match_score = scoring_scheme.score(seqan3::assign_rank_to(0, alphabet_t{}),
                                                          seqan3::assign_rank_to(0, alphabet_t{}));
        score_t scalar_mismatch_score = scoring_scheme.score(seqan3::assign_rank_to(0, alphabet_t{}),
                                                             seqan3::assign_rank_to(1, alphabet_t{}));
 
        // Check if the scoring scheme match and mismatch scores do not overflow with the respective scalar type.
        if constexpr (sizeof(simd_scalar_t) < sizeof(score_t))
        {
            if (scalar_match_score > static_cast<score_t>(std::numeric_limits<simd_scalar_t>::max()) ||
                scalar_mismatch_score < static_cast<score_t>(std::numeric_limits<simd_scalar_t>::lowest()))
            {
                throw std::invalid_argument{"The selected scoring scheme score overflows "
                                            "for the selected scalar type of the simd type."};
            }
        }
 
        match_score = simd::fill<simd_score_t>(static_cast<simd_scalar_t>(scalar_match_score));
        mismatch_score = simd::fill<simd_score_t>(static_cast<simd_scalar_t>(scalar_mismatch_score));
    }
 
    simd_score_t match_score; 
    simd_score_t mismatch_score; 
};
 
} // namespace seqan3::detail