alignment_algorithm.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 <iterator>
#include <memory>
#include <optional>
#include <ranges>
#include <type_traits>
 
#include <seqan3/alignment/configuration/align_config_band.hpp>
#include <seqan3/alignment/configuration/align_config_scoring_scheme.hpp>
#include <seqan3/alignment/exception.hpp>
#include <seqan3/alignment/matrix/detail/aligned_sequence_builder.hpp>
#include <seqan3/alignment/matrix/detail/trace_directions.hpp>
#include <seqan3/alignment/pairwise/align_result_selector.hpp>
#include <seqan3/alignment/pairwise/detail/concept.hpp>
#include <seqan3/alignment/pairwise/detail/type_traits.hpp>
#include <seqan3/core/detail/deferred_crtp_base.hpp>
#include <seqan3/core/detail/empty_type.hpp>
#include <seqan3/utility/container/aligned_allocator.hpp>
#include <seqan3/utility/simd/concept.hpp>
#include <seqan3/utility/simd/simd.hpp>
#include <seqan3/utility/simd/simd_traits.hpp>
#include <seqan3/utility/simd/views/to_simd.hpp>
#include <seqan3/utility/type_traits/function_traits.hpp>
#include <seqan3/utility/views/elements.hpp>
 
namespace seqan3::detail
{
 
template <typename config_t, typename... algorithm_policies_t>
class alignment_algorithm :
    public invoke_deferred_crtp_base<algorithm_policies_t, alignment_algorithm<config_t, algorithm_policies_t...>>...
{
private:
    using traits_t = alignment_configuration_traits<config_t>;
 
    template <typename alignment_algorithm_t = alignment_algorithm>
    static auto _alignment_column_t() -> decltype(std::declval<alignment_algorithm_t>().current_alignment_column());
 
    using alignment_column_t = decltype(_alignment_column_t());
    using alignment_column_iterator_t = std::ranges::iterator_t<alignment_column_t>;
    using alignment_result_t = typename traits_t::alignment_result_type;
 
    static_assert(!std::same_as<alignment_result_t, empty_type>, "Alignment result type was not configured.");
 
    using score_debug_matrix_t =
        std::conditional_t<traits_t::is_debug,
                           two_dimensional_matrix<std::optional<typename traits_t::original_score_type>,
                                                  std::allocator<std::optional<typename traits_t::original_score_type>>,
                                                  matrix_major_order::column>,
                           empty_type>;
    using trace_debug_matrix_t =
        std::conditional_t<traits_t::is_debug,
                           two_dimensional_matrix<std::optional<trace_directions>,
                                                  std::allocator<std::optional<trace_directions>>,
                                                  matrix_major_order::column>,
                           empty_type>;
 
public:
    constexpr alignment_algorithm() = default;                                        
    constexpr alignment_algorithm(alignment_algorithm const &) = default;             
    constexpr alignment_algorithm(alignment_algorithm &&) = default;                  
    constexpr alignment_algorithm & operator=(alignment_algorithm const &) = default; 
    constexpr alignment_algorithm & operator=(alignment_algorithm &&) = default;      
    ~alignment_algorithm() = default;                                                 
 
    explicit constexpr alignment_algorithm(config_t const & cfg) :
        invoke_deferred_crtp_base<algorithm_policies_t, alignment_algorithm<config_t, algorithm_policies_t...>>{cfg}...,
        cfg_ptr{std::make_shared<config_t>(cfg)}
    {
        this->scoring_scheme = seqan3::get<align_cfg::scoring_scheme>(*cfg_ptr).scheme;
        this->initialise_alignment_state(*cfg_ptr);
    }
 
    template <indexed_sequence_pair_range indexed_sequence_pairs_t, typename callback_t>
        requires (!traits_t::is_vectorised) && std::invocable<callback_t, alignment_result_t>
    void operator()(indexed_sequence_pairs_t && indexed_sequence_pairs, callback_t && callback)
    {
        using std::get;
 
        for (auto && [sequence_pair, idx] : indexed_sequence_pairs)
            compute_single_pair(idx, get<0>(sequence_pair), get<1>(sequence_pair), callback);
    }
 
    template <indexed_sequence_pair_range indexed_sequence_pairs_t, typename callback_t>
        requires traits_t::is_vectorised && std::invocable<callback_t, alignment_result_t>
    void operator()(indexed_sequence_pairs_t && indexed_sequence_pairs, callback_t && callback)
    {
        assert(cfg_ptr != nullptr);
 
        static_assert(simd_concept<typename traits_t::score_type>, "Expected simd score type.");
        static_assert(simd_concept<typename traits_t::trace_type>, "Expected simd trace type.");
 
        // Extract the batch of sequences for the first and the second sequence.
        auto sequence1_range = indexed_sequence_pairs | views::elements<0> | views::elements<0>;
        auto sequence2_range = indexed_sequence_pairs | views::elements<0> | views::elements<1>;
 
        // Initialise the find_optimum policy in the simd case.
        this->initialise_find_optimum_policy(sequence1_range,
                                             sequence2_range,
                                             this->scoring_scheme.padding_match_score());
 
        // Convert batch of sequences to sequence of simd vectors.
        auto simd_sequences1 = convert_batch_of_sequences_to_simd_vector(sequence1_range);
        auto simd_sequences2 = convert_batch_of_sequences_to_simd_vector(sequence2_range);
 
        max_size_in_collection = std::pair{simd_sequences1.size(), simd_sequences2.size()};
        // Reset the alignment state's optimum between executions of the alignment algorithm.
        this->alignment_state.reset_optimum();
 
        compute_matrix(simd_sequences1, simd_sequences2);
 
        make_alignment_result(indexed_sequence_pairs, callback);
    }
 
private:
    template <typename sequence_range_t>
    constexpr auto convert_batch_of_sequences_to_simd_vector(sequence_range_t & sequences)
    {
        assert(static_cast<size_t>(std::ranges::distance(sequences)) <= traits_t::alignments_per_vector);
 
        using simd_score_t = typename traits_t::score_type;
 
        std::vector<simd_score_t, aligned_allocator<simd_score_t, alignof(simd_score_t)>> simd_sequence{};
 
        for (auto && simd_vector_chunk : sequences | views::to_simd<simd_score_t>(this->scoring_scheme.padding_symbol))
            for (auto && simd_vector : simd_vector_chunk)
                simd_sequence.push_back(std::move(simd_vector));
 
        return simd_sequence;
    }
 
    template <std::ranges::forward_range sequence1_t, std::ranges::forward_range sequence2_t, typename callback_t>
    constexpr void
    compute_single_pair(size_t const idx, sequence1_t && sequence1, sequence2_t && sequence2, callback_t & callback)
    {
        assert(cfg_ptr != nullptr);
 
        if constexpr (traits_t::is_debug)
            initialise_debug_matrices(sequence1, sequence2);
 
        // Reset the alignment state's optimum between executions of the alignment algorithm.
        this->alignment_state.reset_optimum();
 
        if constexpr (traits_t::is_banded)
        {
            using seqan3::get;
            // Get the band and check if band configuration is valid.
            auto const & band = get<align_cfg::band_fixed_size>(*cfg_ptr);
            check_valid_band_parameter(sequence1, sequence2, band);
            auto && [subsequence1, subsequence2] = this->slice_sequences(sequence1, sequence2, band);
            // It would be great to use this interface here instead
            compute_matrix(subsequence1, subsequence2, band);
            make_alignment_result(idx, subsequence1, subsequence2, callback);
        }
        else
        {
            compute_matrix(sequence1, sequence2);
            make_alignment_result(idx, sequence1, sequence2, callback);
        }
    }
 
    template <typename sequence1_t, typename sequence2_t>
    constexpr void check_valid_band_parameter(sequence1_t && sequence1,
                                              sequence2_t && sequence2,
                                              align_cfg::band_fixed_size const & band)
    {
        static_assert(config_t::template exists<align_cfg::band_fixed_size>(),
                      "The band configuration is required for the banded alignment algorithm.");
 
        using diff_type = std::iter_difference_t<std::ranges::iterator_t<sequence1_t>>;
        static_assert(std::is_signed_v<diff_type>, "Only signed types can be used to test the band parameters.");
 
        if (static_cast<diff_type>(band.lower_diagonal) > std::ranges::distance(sequence1))
        {
            throw invalid_alignment_configuration{
                "Invalid band error: The lower diagonal excludes the whole alignment matrix."};
        }
 
        if (static_cast<diff_type>(band.upper_diagonal) < -std::ranges::distance(sequence2))
        {
            throw invalid_alignment_configuration{
                "Invalid band error: The upper diagonal excludes the whole alignment matrix."};
        }
    }
 
    template <typename sequence1_t, typename sequence2_t>
    constexpr void initialise_debug_matrices(sequence1_t & sequence1, sequence2_t & sequence2)
    {
        size_t rows = std::ranges::distance(sequence2) + 1;
        size_t cols = std::ranges::distance(sequence1) + 1;
 
        score_debug_matrix = score_debug_matrix_t{number_rows{rows}, number_cols{cols}};
        trace_debug_matrix = trace_debug_matrix_t{number_rows{rows}, number_cols{cols}};
    }
 
    template <typename sequence1_t, typename sequence2_t>
    void compute_matrix(sequence1_t & sequence1, sequence2_t & sequence2)
        requires (!traits_t::is_banded)
    {
        // ----------------------------------------------------------------------------
        // Initialisation phase: allocate memory and initialise first column.
        // ----------------------------------------------------------------------------
 
        this->allocate_matrix(sequence1, sequence2);
        initialise_first_alignment_column(sequence2);
 
        // ----------------------------------------------------------------------------
        // Recursion phase: compute column-wise the alignment matrix.
        // ----------------------------------------------------------------------------
 
        for (auto const & alphabet1 : sequence1)
        {
            compute_alignment_column<true>(this->scoring_scheme_profile_column(alphabet1), sequence2);
            finalise_last_cell_in_column(true);
        }
 
        // ----------------------------------------------------------------------------
        // Wrap up phase: track score in last column and prepare the alignment result.
        // ----------------------------------------------------------------------------
 
        finalise_alignment();
    }
 
    template <typename sequence1_t, typename sequence2_t>
    void compute_matrix(sequence1_t & sequence1, sequence2_t & sequence2, align_cfg::band_fixed_size const & band)
        requires (traits_t::is_banded)
    {
        // ----------------------------------------------------------------------------
        // Initialisation phase: allocate memory and initialise first column.
        // ----------------------------------------------------------------------------
 
        // Allocate and initialise first column.
        this->allocate_matrix(sequence1, sequence2, band, this->alignment_state);
        using row_index_t = std::ranges::range_difference_t<sequence2_t>;
        row_index_t last_row_index = this->score_matrix.band_row_index;
        initialise_first_alignment_column(std::views::take(sequence2, last_row_index));
 
        // ----------------------------------------------------------------------------
        // 1st recursion phase: iterate as long as the band intersects with the first row.
        // ----------------------------------------------------------------------------
 
        row_index_t sequence2_size = std::ranges::distance(sequence2);
        for (auto const & seq1_value : std::views::take(sequence1, this->score_matrix.band_col_index))
        {
            compute_alignment_column<true>(seq1_value, std::views::take(sequence2, ++last_row_index));
            // Only if band reached last row of matrix the last cell might be tracked.
            finalise_last_cell_in_column(last_row_index >= sequence2_size);
        }
 
        // ----------------------------------------------------------------------------
        // 2nd recursion phase: iterate until the end of the matrix.
        // ----------------------------------------------------------------------------
 
        size_t first_row_index = 0;
        for (auto const & seq1_value : std::views::drop(sequence1, this->score_matrix.band_col_index))
        {
            // In the second phase the band moves in every column one base down on the second sequence.
            compute_alignment_column<false>(seq1_value, sequence2 | views::slice(first_row_index++, ++last_row_index));
            // Only if band reached last row of matrix the last cell might be tracked.
            finalise_last_cell_in_column(last_row_index >= sequence2_size);
        }
 
        // ----------------------------------------------------------------------------
        // Wrap up phase: track score in last column and prepare the alignment result.
        // ----------------------------------------------------------------------------
 
        finalise_alignment();
    }
 
    template <typename sequence2_t>
    auto initialise_first_alignment_column(sequence2_t && sequence2)
    {
        // Get the initial column.
        alignment_column = this->current_alignment_column();
        assert(!alignment_column.empty()); // Must contain at least one element.
 
        // Initialise first cell.
        alignment_column_it = alignment_column.begin();
        this->init_origin_cell(*alignment_column_it, this->alignment_state);
 
        // Initialise the remaining cells of this column.
        for (auto it = std::ranges::begin(sequence2); it != std::ranges::end(sequence2); ++it)
            this->init_column_cell(*++alignment_column_it, this->alignment_state);
 
        // Finalise the last cell of the initial column.
        bool at_last_row = true;
        if constexpr (traits_t::is_banded) // If the band reaches until the last row of the matrix.
            at_last_row = static_cast<size_t>(this->score_matrix.band_row_index) == this->score_matrix.num_rows - 1;
 
        finalise_last_cell_in_column(at_last_row);
    }
 
    template <bool initialise_first_cell, typename sequence1_value_t, typename sequence2_t>
    void compute_alignment_column(sequence1_value_t const & seq1_value, sequence2_t && sequence2)
    {
        this->next_alignment_column(); // move to next column and set alignment column iterator accordingly.
        alignment_column = this->current_alignment_column();
        alignment_column_it = alignment_column.begin();
 
        auto seq2_it = std::ranges::begin(sequence2);
 
        if constexpr (initialise_first_cell) // Initialise first cell if it intersects with the first row of the matrix.
        {
            this->init_row_cell(*alignment_column_it, this->alignment_state);
        }
        else // Compute first cell of banded column if it does not intersect with the first row of the matrix.
        {
            this->compute_first_band_cell(*alignment_column_it,
                                          this->alignment_state,
                                          this->scoring_scheme.score(seq1_value, *seq2_it));
            ++seq2_it;
        }
 
        for (; seq2_it != std::ranges::end(sequence2); ++seq2_it)
            this->compute_cell(*++alignment_column_it,
                               this->alignment_state,
                               this->scoring_scheme.score(seq1_value, *seq2_it));
    }
 
    constexpr void finalise_last_cell_in_column(bool const at_last_row) noexcept
    {
        if (at_last_row)
            this->check_score_of_last_row_cell(*alignment_column_it, this->alignment_state);
 
        if constexpr (traits_t::is_debug)
            dump_alignment_column();
    }
 
    constexpr void finalise_alignment() noexcept
    {
        // ----------------------------------------------------------------------------
        // Check for the optimum in last cell/column.
        // ----------------------------------------------------------------------------
 
        this->check_score_of_cells_in_last_column(alignment_column, this->alignment_state);
        this->check_score_of_last_cell(*alignment_column_it, this->alignment_state);
    }
 
    template <typename index_t, typename sequence1_t, typename sequence2_t, typename callback_t>
        requires (!traits_t::is_vectorised)
    constexpr void make_alignment_result([[maybe_unused]] index_t const idx,
                                         [[maybe_unused]] sequence1_t & sequence1,
                                         [[maybe_unused]] sequence2_t & sequence2,
                                         callback_t & callback)
    {
        using result_value_t = typename alignment_result_value_type_accessor<alignment_result_t>::type;
 
        // ----------------------------------------------------------------------------
        // Build the alignment result
        // ----------------------------------------------------------------------------
 
        static_assert(seqan3::detail::alignment_configuration_traits<config_t>::has_output_configuration,
                      "The configuration must contain at least one align_cfg::output_* element.");
 
        result_value_t res{};
 
        if constexpr (traits_t::output_sequence1_id)
            res.sequence1_id = idx;
 
        if constexpr (traits_t::output_sequence2_id)
            res.sequence2_id = idx;
 
        // Choose what needs to be computed.
        if constexpr (traits_t::compute_score)
            res.score = this->alignment_state.optimum.score;
 
        if constexpr (traits_t::compute_end_positions)
        {
            using alignment_coordinate_t = detail::advanceable_alignment_coordinate<>;
            res.end_positions = alignment_coordinate_t{column_index_type{this->alignment_state.optimum.column_index},
                                                       row_index_type{this->alignment_state.optimum.row_index}};
            // At some point this needs to be refactored so that it is not necessary to adapt the coordinate.
            if constexpr (traits_t::is_banded)
            {
                res.end_positions.second += res.end_positions.first - this->trace_matrix.band_col_index;
                res.end_positions.first = this->to_original_sequence1_position(res.end_positions.first);
                res.end_positions.second = this->to_original_sequence2_position(res.end_positions.second);
            }
        }
 
        if constexpr (traits_t::compute_begin_positions)
        {
            // Get a aligned sequence builder for banded or un-banded case.
            aligned_sequence_builder builder{sequence1, sequence2};
 
            detail::matrix_coordinate const optimum_coordinate{
                detail::row_index_type{this->alignment_state.optimum.row_index},
                detail::column_index_type{this->alignment_state.optimum.column_index}};
            auto trace_res = builder(this->trace_matrix.trace_path(optimum_coordinate));
            res.begin_positions.first =
                this->to_original_sequence1_position(trace_res.first_sequence_slice_positions.first);
            res.begin_positions.second =
                this->to_original_sequence2_position(trace_res.second_sequence_slice_positions.first);
 
            if constexpr (traits_t::compute_sequence_alignment)
                res.alignment = std::move(trace_res.alignment);
        }
 
        // Store the matrices in debug mode.
        if constexpr (traits_t::is_debug)
        {
            res.score_debug_matrix = std::move(score_debug_matrix);
            if constexpr (traits_t::compute_sequence_alignment) // compute alignment
                res.trace_debug_matrix = std::move(trace_debug_matrix);
        }
 
        callback(std::move(res));
    }
 
    template <typename indexed_sequence_pair_range_t, typename callback_t>
        requires traits_t::is_vectorised
    constexpr auto make_alignment_result(indexed_sequence_pair_range_t && index_sequence_pairs, callback_t & callback)
    {
        using result_value_t = typename alignment_result_value_type_accessor<alignment_result_t>::type;
 
        size_t simd_index = 0;
        for (auto && [sequence_pairs, alignment_index] : index_sequence_pairs)
        {
            (void)sequence_pairs;
            result_value_t res{};
 
            if constexpr (traits_t::output_sequence1_id)
                res.sequence1_id = alignment_index;
 
            if constexpr (traits_t::output_sequence2_id)
                res.sequence2_id = alignment_index;
 
            if constexpr (traits_t::compute_score)
                res.score = this->alignment_state.optimum.score[simd_index]; // Just take this
 
            if constexpr (traits_t::compute_end_positions)
            {
                res.end_positions.first =
                    this->to_original_sequence1_position(this->alignment_state.optimum.column_index[simd_index]);
                res.end_positions.second =
                    this->to_original_sequence2_position(this->alignment_state.optimum.row_index[simd_index]);
            }
 
            callback(std::move(res));
            ++simd_index;
        }
    }
 
    void dump_alignment_column()
    {
        using std::get;
 
        auto column = this->current_alignment_column();
 
        auto coord = get<1>(column.front()).coordinate;
        if constexpr (traits_t::is_banded)
            coord.second += coord.first - this->score_matrix.band_col_index;
 
        matrix_offset offset{row_index_type{static_cast<std::ptrdiff_t>(coord.second)},
                             column_index_type{static_cast<std::ptrdiff_t>(coord.first)}};
 
        std::ranges::copy(column
                              | std::views::transform(
                                  [](auto const & tpl)
                                  {
                                      using std::get;
                                      return get<0>(tpl).current;
                                  }),
                          score_debug_matrix.begin() + offset);
 
        // if traceback is enabled.
        if constexpr (traits_t::compute_sequence_alignment)
        {
            std::ranges::copy(
                column
                    | std::views::transform(
                        [](auto const & tpl)
                        {
                            using std::get;
                            auto trace = get<1>(tpl).current;
 
                            if (auto _up = (trace & trace_directions::up_open); _up == trace_directions::carry_up_open)
                                trace = trace ^ trace_directions::carry_up_open; // remove silent up open signal
                            else if (_up == trace_directions::up_open)
                                trace = trace ^ trace_directions::up; // display up open only with single bit.
 
                            if (auto _left = (trace & trace_directions::left_open);
                                _left == trace_directions::carry_left_open)
                                trace = trace ^ trace_directions::carry_left_open; // remove silent left open signal
                            else if (_left == trace_directions::left_open)
                                trace = trace ^ trace_directions::left; // display left open only with single bit.
 
                            return trace;
                        }),
                trace_debug_matrix.begin() + offset);
        }
    }
 
    std::shared_ptr<config_t> cfg_ptr{};
    alignment_column_t alignment_column{};
    alignment_column_iterator_t alignment_column_it{};
    score_debug_matrix_t score_debug_matrix{};
    trace_debug_matrix_t trace_debug_matrix{};
    std::pair<size_t, size_t> max_size_in_collection{};
};
 
} // namespace seqan3::detail