simd_find_optimum_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 <type_traits>
 
#include <seqan3/alignment/matrix/detail/alignment_optimum.hpp>
#include <seqan3/alignment/pairwise/detail/alignment_algorithm_state.hpp>
#include <seqan3/alignment/pairwise/policy/find_optimum_policy.hpp>
#include <seqan3/core/detail/empty_type.hpp>
#include <seqan3/utility/simd/concept.hpp>
#include <seqan3/utility/views/zip.hpp>
 
namespace seqan3::detail
{
 
template <simd::simd_concept simd_t>
struct simd_global_alignment_state
{
    simd_t score_offset{};
    simd_t coordinate_offset{};
    simd_t last_column_mask{};
    simd_t last_row_mask{};
};
 
template <typename alignment_algorithm_t, simd::simd_concept simd_t>
class simd_find_optimum_policy : public simd_global_alignment_state<simd_t>
{
private:
    friend alignment_algorithm_t;
 
    bool is_global_alignment{false};
    bool test_every_cell{false};
    bool test_last_row_cell{false};
    bool test_last_column_cell{false};
 
    constexpr simd_find_optimum_policy() = default;                                             
    constexpr simd_find_optimum_policy(simd_find_optimum_policy const &) = default;             
    constexpr simd_find_optimum_policy(simd_find_optimum_policy &&) = default;                  
    constexpr simd_find_optimum_policy & operator=(simd_find_optimum_policy const &) = default; 
    constexpr simd_find_optimum_policy & operator=(simd_find_optimum_policy &&) = default;      
    ~simd_find_optimum_policy() = default;                                                      
 
    template <typename configuration_t>
    simd_find_optimum_policy(configuration_t const & config)
    {
        if constexpr (configuration_t::template exists<align_cfg::method_local>())
            test_every_cell = true;
 
        is_global_alignment = configuration_t::template exists<align_cfg::method_global>();
 
        auto method_global_config = config.get_or(align_cfg::method_global{});
 
        test_last_row_cell = method_global_config.free_end_gaps_sequence1_trailing || is_global_alignment;
        test_last_column_cell = method_global_config.free_end_gaps_sequence2_trailing || is_global_alignment;
    }
 
protected:
    template <typename cell_t, typename score_t>
    constexpr void check_score_of_cell([[maybe_unused]] cell_t const & current_cell,
                                       [[maybe_unused]] alignment_algorithm_state<score_t> & state) const noexcept
    {
        if (test_every_cell)
            check_and_update(current_cell, state);
    }
 
    template <typename other_alignment_algorithm_t, simd_concept score_t, typename is_local_t>
    friend class simd_affine_gap_policy;
 
    template <typename other_alignment_algorithm_t>
    friend class affine_gap_init_policy;
 
    template <typename cell_t, typename score_t>
    constexpr void
    check_score_of_last_row_cell([[maybe_unused]] cell_t const & last_row_cell,
                                 [[maybe_unused]] alignment_algorithm_state<score_t> & state) const noexcept
    {
        // Only search in last row if requested and not done already.}
        if (!test_every_cell && test_last_row_cell)
        {
            if (is_global_alignment)
                check_and_update<true>(last_row_cell, state);
            else
                check_and_update(last_row_cell, state);
        }
    }
 
    template <typename alignment_column_t, typename score_t>
    constexpr void
    check_score_of_cells_in_last_column([[maybe_unused]] alignment_column_t && last_column,
                                        [[maybe_unused]] alignment_algorithm_state<score_t> & state) const noexcept
    {
        // Only check last cell if not done before.
        if (!test_every_cell && test_last_column_cell)
        {
            for (auto && cell : last_column)
            {
                if (is_global_alignment)
                    check_and_update<false>(cell, state);
                else
                    check_and_update(cell, state);
            }
        }
    }
 
    template <typename cell_t, typename score_t>
    constexpr void check_score_of_last_cell([[maybe_unused]] cell_t const & last_cell,
                                            [[maybe_unused]] alignment_algorithm_state<score_t> & state) const noexcept
    {
        // Only check last cell if not done before.
        if (!(test_every_cell || test_last_row_cell || test_last_column_cell))
            check_and_update(last_cell, state);
    }
 
    template <std::ranges::forward_range sequence1_collection_t,
              std::ranges::forward_range sequence2_collection_t,
              arithmetic score_t>
    void initialise_find_optimum_policy([[maybe_unused]] sequence1_collection_t && sequence1_collection,
                                        [[maybe_unused]] sequence2_collection_t && sequence2_collection,
                                        [[maybe_unused]] score_t const padding_score)
    {
        if (is_global_alignment)
        {
            assert(std::ranges::distance(sequence1_collection) == std::ranges::distance(sequence2_collection));
 
            constexpr size_t simd_size = simd_traits<simd_t>::length;
            // First get global size.
            std::array<size_t, simd_size> sequence1_sizes{};
            std::array<size_t, simd_size> sequence2_sizes{};
 
            std::ptrdiff_t max_sequence1_size{};
            std::ptrdiff_t max_sequence2_size{};
 
            size_t array_index{};
            for (auto && [sequence1, sequence2] : views::zip(sequence1_collection, sequence2_collection))
            {
                sequence1_sizes[array_index] = std::ranges::distance(sequence1);
                sequence2_sizes[array_index] = std::ranges::distance(sequence2);
                max_sequence1_size = std::max<std::ptrdiff_t>(sequence1_sizes[array_index], max_sequence1_size);
                max_sequence2_size = std::max<std::ptrdiff_t>(sequence2_sizes[array_index], max_sequence2_size);
                ++array_index;
            }
 
            // The global diagonal ending in the sink of the outer alignment matrix.
            std::ptrdiff_t global_diagonal = max_sequence1_size - max_sequence2_size;
 
            for (size_t simd_index = 0; simd_index < simd_size; ++simd_index)
            {
                if (std::ptrdiff_t local_diagonal = sequence1_sizes[simd_index] - sequence2_sizes[simd_index];
                    local_diagonal < global_diagonal)
                { // optimum is stored in last row.
                    this->last_row_mask[simd_index] = max_sequence1_size - (global_diagonal - local_diagonal);
                    this->last_column_mask[simd_index] = max_sequence2_size + 1;
                    this->coordinate_offset[simd_index] = max_sequence2_size - sequence2_sizes[simd_index];
                    this->score_offset[simd_index] = padding_score * this->coordinate_offset[simd_index];
                }
                else // optimum is stored in last column
                {
                    this->last_column_mask[simd_index] = max_sequence2_size - (local_diagonal - global_diagonal);
                    this->last_row_mask[simd_index] = max_sequence1_size + 1;
                    this->coordinate_offset[simd_index] = max_sequence1_size - sequence1_sizes[simd_index];
                    this->score_offset[simd_index] = padding_score * this->coordinate_offset[simd_index];
                }
            }
        }
        // else no-op
    }
 
private:
    template <typename cell_t, typename score_t>
    constexpr void check_and_update(cell_t const & cell, alignment_algorithm_state<score_t> & state) const noexcept
    {
        assert(!is_global_alignment); // This function should not be called for the global alignment.
 
        auto const & [score_cell, trace_cell] = cell;
        state.optimum.update_if_new_optimal_score(score_cell.current,
                                                  column_index_type{trace_cell.coordinate.first},
                                                  row_index_type{trace_cell.coordinate.second});
    }
 
    template <bool in_last_row, typename cell_t, typename score_t>
    constexpr void check_and_update(cell_t const & cell, alignment_algorithm_state<score_t> & state) const noexcept
    {
        assert(is_global_alignment); // This function should only be called for the global alignment.
 
        using simd_mask_t = typename simd_traits<simd_t>::mask_type;
        auto const & [score_cell, trace_cell] = cell;
        simd_t column_positions = simd::fill<simd_t>(trace_cell.coordinate.first);
        simd_t row_positions = simd::fill<simd_t>(trace_cell.coordinate.second);
 
        simd_mask_t mask{};
 
        if constexpr (in_last_row) // Check if column was masked
            mask = (column_positions == this->last_row_mask);
        else // Check if row was masked
            mask = (row_positions == this->last_column_mask);
 
        // In global alignment we are only interested in the position not the max of the scores.
        // In addition, the scores need to be corrected in order to track the right score.
        state.optimum.score = mask ? score_cell.current - this->score_offset : state.optimum.score;
        state.optimum.column_index = mask ? column_positions - this->coordinate_offset : state.optimum.column_index;
        state.optimum.row_index = mask ? row_positions - this->coordinate_offset : state.optimum.row_index;
    }
};
 
} // namespace seqan3::detail