simd_find_optimum_policy.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2020, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2020, 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 <type_traits>
 
#include <seqan3/alignment/matrix/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/core/simd/concept.hpp>
#include <seqan3/range/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,
          typename is_global_alignment_t,
          typename traits_type = default_find_optimum_trait>
class simd_find_optimum_policy : public std::conditional_t<is_global_alignment_t::value,
                                                           simd_global_alignment_state<simd_t>,
                                                           empty_type>
{
private:
    friend alignment_algorithm_t;
 
    static constexpr bool is_global_alignment = is_global_alignment_t::value;
    static constexpr bool search_in_every_cell = traits_type::find_in_every_cell_type::value;
    static constexpr bool search_in_last_row = traits_type::find_in_last_row_type::value || is_global_alignment;
    static constexpr bool search_in_last_column = traits_type::find_in_last_column_type::value || is_global_alignment;
 
    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; 
 
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 constexpr (search_in_every_cell)
            check_and_update(current_cell, state);
    }
 
    template <typename other_alignment_algorithm_t, typename score_t, typename is_local_t>
    friend class simd_affine_gap_policy;
 
    template <typename other_alignment_algorithm_t, typename other_traits_type>
    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 constexpr (!search_in_every_cell && search_in_last_row)
        {
            if constexpr (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 constexpr (!search_in_every_cell && search_in_last_column)
        {
            for (auto && cell : last_column)
            {
                if constexpr (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 constexpr (!search_in_every_cell && !search_in_last_row && !search_in_last_column)
            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 constexpr (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
    {
        static_assert(!is_global_alignment, "This function should not be called for the global alignment. "
                                            "Use the other check_and_update function instead.");
        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
    {
        static_assert(is_global_alignment, "This function should only be called for the global alignment. "
                                           "Use the other check_and_update function instead.");
 
        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