simd_affine_gap_policy.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 <limits>
#include <tuple>
 
#include <seqan3/alignment/configuration/align_config_gap_cost_affine.hpp>
#include <seqan3/alignment/matrix/detail/alignment_optimum.hpp>
#include <seqan3/alignment/matrix/detail/trace_directions.hpp>
#include <seqan3/alignment/pairwise/detail/alignment_algorithm_state.hpp>
#include <seqan3/core/configuration/configuration.hpp>
#include <seqan3/utility/concept.hpp>
#include <seqan3/utility/simd/algorithm.hpp>
#include <seqan3/utility/simd/concept.hpp>
#include <seqan3/utility/simd/simd_traits.hpp>
 
namespace seqan3::detail
{
 
// ----------------------------------------------------------------------------
// simd_affine_gap_policy
// ----------------------------------------------------------------------------
 
template <typename alignment_algorithm_t, simd_concept score_t, typename align_local_t = std::false_type>
class simd_affine_gap_policy
{
private:
    friend alignment_algorithm_t;
 
    using alignment_state_t = alignment_algorithm_state<score_t>;
 
    constexpr simd_affine_gap_policy() noexcept = default;                                           
    constexpr simd_affine_gap_policy(simd_affine_gap_policy const &) noexcept = default;             
    constexpr simd_affine_gap_policy(simd_affine_gap_policy &&) noexcept = default;                  
    constexpr simd_affine_gap_policy & operator=(simd_affine_gap_policy const &) noexcept = default; 
    constexpr simd_affine_gap_policy & operator=(simd_affine_gap_policy &&) noexcept = default;      
    ~simd_affine_gap_policy() noexcept = default;                                                    
 
    template <typename configuration_t>
    simd_affine_gap_policy(configuration_t const & /*config*/)
    {}
 
    template <typename cell_t>
    constexpr void
    compute_cell(cell_t && current_cell, alignment_algorithm_state<score_t> & state, score_t const score) const noexcept
    {
        // score_cell = seqan3::detail::alignment_score_matrix_proxy
        // trace_cell = seqan3::detail::alignment_trace_matrix_proxy
        auto & [score_cell, trace_cell] = current_cell;
        constexpr bool with_trace = !decays_to_ignore_v<std::remove_reference_t<decltype(trace_cell.current)>>;
        // Precompute the diagonal score.
        score_t tmp = score_cell.diagonal + score;
 
        if constexpr (with_trace)
        {
            auto mask = tmp < score_cell.up;
            tmp = (mask) ? score_cell.up : tmp;
            trace_cell.current = (mask) ? trace_cell.up : convert_to_simd(trace_directions::diagonal) | trace_cell.up;
 
            mask = tmp < score_cell.r_left;
            tmp = (mask) ? score_cell.r_left : tmp;
            trace_cell.current = (mask) ? trace_cell.r_left : trace_cell.current | trace_cell.r_left;
        }
        else
        {
            tmp = (tmp < score_cell.up) ? score_cell.up : tmp;
            tmp = (tmp < score_cell.r_left) ? score_cell.r_left : tmp;
        }
 
        if constexpr (align_local_t::value)
        {
            tmp = (tmp < simd::fill<score_t>(0))
                    /*then*/
                    ? (trace_cell.current = convert_to_simd(trace_directions::none), simd::fill<score_t>(0))
                    /*else*/
                    : tmp;
        }
 
        // Store the current max score.
        score_cell.current = tmp;
        // Check if this was the optimum. Possibly a noop.
        static_cast<alignment_algorithm_t const &>(*this).check_score_of_cell(current_cell, state);
 
        // Prepare horizontal and vertical score for next column.
        tmp += state.gap_open_score;
        score_cell.up += state.gap_extension_score;
        score_cell.w_left = score_cell.r_left + state.gap_extension_score;
 
        auto mask = score_cell.up < tmp;
        score_cell.up = (mask) ? tmp : score_cell.up;
        trace_cell.up = (mask) ? convert_to_simd(trace_directions::up_open) : convert_to_simd(trace_directions::up);
        mask = score_cell.w_left < tmp;
        score_cell.w_left = (mask) ? tmp : score_cell.w_left;
        trace_cell.w_left =
            (mask) ? convert_to_simd(trace_directions::left_open) : convert_to_simd(trace_directions::left);
    }
 
    template <typename alignment_configuration_t>
    constexpr void initialise_alignment_state(alignment_configuration_t const & config) noexcept
    {
        using scalar_t = typename simd_traits<score_t>::scalar_type;
        auto scheme =
            config.get_or(align_cfg::gap_cost_affine{align_cfg::open_score{-10}, align_cfg::extension_score{-1}});
 
        alignment_state.gap_extension_score = simd::fill<score_t>(static_cast<scalar_t>(scheme.extension_score));
        alignment_state.gap_open_score =
            simd::fill<score_t>(static_cast<scalar_t>(scheme.extension_score + scheme.open_score));
    }
 
    constexpr score_t convert_to_simd(trace_directions const direction) const noexcept
    {
        using scalar_t = typename simd_traits<score_t>::scalar_type;
 
        return simd::fill<score_t>(static_cast<scalar_t>(direction));
    }
 
    alignment_state_t alignment_state{}; 
};
 
} // namespace seqan3::detail