policy_optimum_tracker_simd.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 <ranges>
 
#include <seqan3/alignment/pairwise/detail/policy_optimum_tracker.hpp>
#include <seqan3/utility/simd/algorithm.hpp>
#include <seqan3/utility/simd/simd.hpp>
#include <seqan3/utility/simd/simd_traits.hpp>
#include <seqan3/utility/type_traits/lazy_conditional.hpp>
#include <seqan3/utility/views/zip.hpp>
 
namespace seqan3::detail
{
 
struct max_score_updater_simd_global
{
 
    template <typename score_t, typename coordinate_t>
        requires (std::assignable_from<score_t &, score_t const &> &&
                  requires (coordinate_t coordinate)
                  {
                      requires simd_concept<decltype(coordinate.col)>;
                      requires simd_concept<decltype(coordinate.row)>;
                  })
    void operator()(score_t & optimal_score,
                    coordinate_t const & optimal_coordinate,
                    score_t current_score,
                    coordinate_t const & current_coordinate) const noexcept
    {
        auto mask =
            (optimal_coordinate.col == current_coordinate.col) && (optimal_coordinate.row == current_coordinate.row);
        optimal_score = (mask) ? std::move(current_score) : optimal_score;
    }
};
 
template <typename alignment_configuration_t, std::semiregular optimum_updater_t>
    requires is_type_specialisation_of_v<alignment_configuration_t, configuration>
          && std::invocable<
                 optimum_updater_t,
                 typename alignment_configuration_traits<alignment_configuration_t>::score_type &,
                 typename alignment_configuration_traits<alignment_configuration_t>::matrix_coordinate_type &,
                 typename alignment_configuration_traits<alignment_configuration_t>::score_type,
                 typename alignment_configuration_traits<alignment_configuration_t>::matrix_coordinate_type>
class policy_optimum_tracker_simd : protected policy_optimum_tracker<alignment_configuration_t, optimum_updater_t>
{
protected:
    using base_policy_t = policy_optimum_tracker<alignment_configuration_t, optimum_updater_t>;
 
    // Import the configured score type.
    using typename base_policy_t::score_type;
    using typename base_policy_t::traits_type;
 
    using scalar_type = typename simd::simd_traits<score_type>::scalar_type;
    using original_score_type = typename traits_type::original_score_type;
 
    static_assert(simd_concept<score_type>, "Must be a simd type!");
 
    // Import base variables into class scope.
    using base_policy_t::compare_and_set_optimum;
    using base_policy_t::optimal_coordinate;
    using base_policy_t::optimal_score;
    std::array<original_score_type, simd_traits<score_type>::length> padding_offsets{};
 
    policy_optimum_tracker_simd() = default;                                                
    policy_optimum_tracker_simd(policy_optimum_tracker_simd const &) = default;             
    policy_optimum_tracker_simd(policy_optimum_tracker_simd &&) = default;                  
    policy_optimum_tracker_simd & operator=(policy_optimum_tracker_simd const &) = default; 
    policy_optimum_tracker_simd & operator=(policy_optimum_tracker_simd &&) = default;      
    ~policy_optimum_tracker_simd() = default;                                               
 
    policy_optimum_tracker_simd(alignment_configuration_t const & config) : base_policy_t{config}
    {
        base_policy_t::test_last_row_cell = true;
        base_policy_t::test_last_column_cell = true;
    }
 
    void reset_optimum()
    {
        optimal_score = simd::fill<score_type>(std::numeric_limits<scalar_type>::lowest());
    }
 
    template <std::ranges::input_range sequence1_collection_t, std::ranges::input_range sequence2_collection_t>
    void initialise_tracker(sequence1_collection_t & sequence1_collection,
                            sequence2_collection_t & sequence2_collection)
    {
        using index_t = typename traits_type::matrix_index_type;
        using scalar_index_t = typename simd_traits<index_t>::scalar_type;
 
        scalar_index_t largest_sequence1_size{};
        scalar_index_t largest_sequence2_size{};
        alignas(alignof(index_t)) std::array<scalar_index_t, traits_type::alignments_per_vector> sequence1_sizes{};
        alignas(alignof(index_t)) std::array<scalar_index_t, traits_type::alignments_per_vector> sequence2_sizes{};
 
        // First, get all dimensions from the sequences and keep track of the maximal size in either dimension.
        size_t sequence_count{};
        for (auto && [sequence1, sequence2] : views::zip(sequence1_collection, sequence2_collection))
        {
            sequence1_sizes[sequence_count] = std::ranges::distance(sequence1);
            sequence2_sizes[sequence_count] = std::ranges::distance(sequence2);
            largest_sequence1_size = std::max(largest_sequence1_size, sequence1_sizes[sequence_count]);
            largest_sequence2_size = std::max(largest_sequence2_size, sequence2_sizes[sequence_count]);
            ++sequence_count;
        }
 
        // Second, determine the offset for each individual end-coordinate which is used to project the cell to the
        // last row or column of the global alignment matrix. Choose the smallest distance as the correct offset
        // to the projected cell.
        for (size_t index = 0; index != sequence_count; ++index)
        {
            assert(sequence1_sizes[index] <= largest_sequence1_size);
            assert(sequence2_sizes[index] <= largest_sequence2_size);
 
            padding_offsets[index] = std::min(largest_sequence1_size - sequence1_sizes[index],
                                              largest_sequence2_size - sequence2_sizes[index]);
            sequence1_sizes[index] += padding_offsets[index];
            sequence2_sizes[index] += padding_offsets[index];
        }
 
        // Load the target coordinate indices from the respective arrays.
        optimal_coordinate.col = simd::load<index_t>(sequence1_sizes.data());
        optimal_coordinate.row = simd::load<index_t>(sequence2_sizes.data());
    }
};
} // namespace seqan3::detail