Loading...
Searching...
No Matches
pairwise_alignment_algorithm_banded.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 <concepts>
#include <ranges>
#include <seqan3/alignment/pairwise/detail/pairwise_alignment_algorithm.hpp>
#include <seqan3/utility/views/slice.hpp>
namespace seqan3::detail
{
template <typename alignment_configuration_t, typename... policies_t>
requires is_type_specialisation_of_v<alignment_configuration_t, configuration>
class pairwise_alignment_algorithm_banded :
protected pairwise_alignment_algorithm<alignment_configuration_t, policies_t...>
{
protected:
using base_algorithm_t = pairwise_alignment_algorithm<alignment_configuration_t, policies_t...>;
// Import types from base class.
using typename base_algorithm_t::alignment_result_type;
using typename base_algorithm_t::score_type;
using typename base_algorithm_t::traits_type;
static_assert(!std::same_as<alignment_result_type, empty_type>, "Alignment result type was not configured.");
static_assert(traits_type::is_banded, "Alignment configuration must have band configured.");
public:
pairwise_alignment_algorithm_banded() = default;
pairwise_alignment_algorithm_banded(pairwise_alignment_algorithm_banded const &) = default;
pairwise_alignment_algorithm_banded(pairwise_alignment_algorithm_banded &&) = default;
pairwise_alignment_algorithm_banded &
operator=(pairwise_alignment_algorithm_banded const &) = default;
pairwise_alignment_algorithm_banded & operator=(pairwise_alignment_algorithm_banded &&) = default;
~pairwise_alignment_algorithm_banded() = default;
pairwise_alignment_algorithm_banded(alignment_configuration_t const & config) : base_algorithm_t(config)
{}
template <indexed_sequence_pair_range indexed_sequence_pairs_t, typename callback_t>
requires std::invocable<callback_t, alignment_result_type>
void operator()(indexed_sequence_pairs_t && indexed_sequence_pairs, callback_t && callback)
{
using std::get;
for (auto && [sequence_pair, idx] : indexed_sequence_pairs)
{
size_t sequence1_size = std::ranges::distance(get<0>(sequence_pair));
size_t const sequence2_size = std::ranges::distance(get<1>(sequence_pair));
auto && [alignment_matrix, index_matrix] =
this->acquire_matrices(sequence1_size, sequence2_size, this->lowest_viable_score());
// Initialise the cell updater with the dimensions of the regular matrix.
this->compare_and_set_optimum.set_target_indices(row_index_type{sequence2_size},
column_index_type{sequence1_size});
// Shrink the first sequence if the band ends before its actual end.
sequence1_size = std::min(sequence1_size, this->upper_diagonal + sequence2_size);
using sequence1_difference_t = std::ranges::range_difference_t<decltype(get<0>(sequence_pair))>;
compute_matrix(std::views::take(get<0>(sequence_pair), static_cast<sequence1_difference_t>(sequence1_size)),
get<1>(sequence_pair),
alignment_matrix,
index_matrix);
this->make_result_and_invoke(std::forward<decltype(sequence_pair)>(sequence_pair),
std::move(idx),
this->optimal_score,
this->optimal_coordinate,
alignment_matrix,
callback);
}
}
template <indexed_sequence_pair_range indexed_sequence_pairs_t, typename callback_t>
requires traits_type::is_vectorised && std::invocable<callback_t, alignment_result_type>
auto operator()(indexed_sequence_pairs_t && indexed_sequence_pairs, callback_t && callback)
{
using simd_collection_t = std::vector<score_type, aligned_allocator<score_type, alignof(score_type)>>;
using original_score_t = typename traits_type::original_score_type;
// Extract the batch of sequences for the first and the second sequence.
auto seq1_collection = indexed_sequence_pairs | views::elements<0> | views::elements<0>;
auto seq2_collection = indexed_sequence_pairs | views::elements<0> | views::elements<1>;
this->initialise_tracker(seq1_collection, seq2_collection);
// Convert batch of sequences to sequence of simd vectors.
thread_local simd_collection_t simd_seq1_collection{};
thread_local simd_collection_t simd_seq2_collection{};
this->convert_batch_of_sequences_to_simd_vector(simd_seq1_collection,
seq1_collection,
this->scoring_scheme.padding_symbol);
this->convert_batch_of_sequences_to_simd_vector(simd_seq2_collection,
seq2_collection,
this->scoring_scheme.padding_symbol);
size_t const sequence1_size = std::ranges::distance(simd_seq1_collection);
size_t const sequence2_size = std::ranges::distance(simd_seq2_collection);
auto && [alignment_matrix, index_matrix] =
this->acquire_matrices(sequence1_size, sequence2_size, this->lowest_viable_score());
compute_matrix(simd_seq1_collection, simd_seq2_collection, alignment_matrix, index_matrix);
size_t index = 0;
for (auto && [sequence_pair, idx] : indexed_sequence_pairs)
{
original_score_t score = this->optimal_score[index]
- (this->padding_offsets[index] * this->scoring_scheme.padding_match_score());
matrix_coordinate coordinate{row_index_type{size_t{this->optimal_coordinate.row[index]}},
column_index_type{size_t{this->optimal_coordinate.col[index]}}};
this->make_result_and_invoke(std::forward<decltype(sequence_pair)>(sequence_pair),
std::move(idx),
std::move(score),
std::move(coordinate),
alignment_matrix,
callback);
++index;
}
}
protected:
template <std::ranges::forward_range sequence1_t,
std::ranges::forward_range sequence2_t,
std::ranges::input_range alignment_matrix_t,
std::ranges::input_range index_matrix_t>
requires std::ranges::forward_range<std::ranges::range_reference_t<alignment_matrix_t>>
&& std::ranges::forward_range<std::ranges::range_reference_t<index_matrix_t>>
void compute_matrix(sequence1_t && sequence1,
sequence2_t && sequence2,
alignment_matrix_t && alignment_matrix,
index_matrix_t && index_matrix)
{
// ---------------------------------------------------------------------
// Initialisation phase: allocate memory and initialise first column.
// ---------------------------------------------------------------------
this->reset_optimum(); // Reset the tracker for the new alignment computation.
auto alignment_matrix_it = alignment_matrix.begin();
auto indexed_matrix_it = index_matrix.begin();
using row_index_t = std::ranges::range_difference_t<sequence2_t>; // row_size = |sequence2| + 1
using column_index_t = std::ranges::range_difference_t<sequence1_t>; // column_size = |sequence1| + 1
row_index_t row_size = std::max<int32_t>(0, -this->lower_diagonal);
column_index_t const column_size = std::max<int32_t>(0, this->upper_diagonal);
this->initialise_column(*alignment_matrix_it, *indexed_matrix_it, std::views::take(sequence2, row_size));
// ---------------------------------------------------------------------
// 1st recursion phase: band intersects with the first row.
// ---------------------------------------------------------------------
for (auto alphabet1 : std::views::take(sequence1, column_size))
{
this->compute_column(*++alignment_matrix_it,
*++indexed_matrix_it,
alphabet1,
std::views::take(sequence2, ++row_size));
}
// ---------------------------------------------------------------------
// 2nd recursion phase: iterate until the end of the matrix.
// ---------------------------------------------------------------------
row_index_t first_row_index = 0u;
for (auto alphabet1 : std::views::drop(sequence1, column_size))
{
compute_band_column(*++alignment_matrix_it,
std::views::drop(*++indexed_matrix_it, first_row_index + 1),
alphabet1,
views::slice(sequence2, first_row_index, ++row_size));
++first_row_index;
}
// ---------------------------------------------------------------------
// Final phase: track score of last column
// ---------------------------------------------------------------------
auto alignment_column = *alignment_matrix_it;
auto cell_index_column = std::views::drop(*indexed_matrix_it, first_row_index);
auto alignment_column_it = alignment_column.begin();
auto cell_index_column_it = cell_index_column.begin();
this->track_last_column_cell(*alignment_column_it, *cell_index_column_it);
for (row_index_t last_row = std::min<row_index_t>(std::ranges::distance(sequence2), row_size);
first_row_index < last_row;
++first_row_index)
this->track_last_column_cell(*++alignment_column_it, *++cell_index_column_it);
this->track_final_cell(*alignment_column_it, *cell_index_column_it);
}
template <std::ranges::forward_range alignment_column_t,
std::ranges::input_range cell_index_column_t,
typename alphabet1_t,
std::ranges::input_range sequence2_t>
void compute_band_column(alignment_column_t && alignment_column,
cell_index_column_t && cell_index_column,
alphabet1_t const & alphabet1,
sequence2_t && sequence2)
{
// ---------------------------------------------------------------------
// Initial phase: prepare column and initialise first cell
// ---------------------------------------------------------------------
auto current_alignment_column_it = alignment_column.begin();
auto cell_index_column_it = cell_index_column.begin();
// Points to the last valid cell in the column.
decltype(current_alignment_column_it) next_alignment_column_it{current_alignment_column_it};
auto cell = *current_alignment_column_it;
cell = this->track_cell(
this->initialise_band_first_cell(cell.best_score(),
*++next_alignment_column_it,
this->scoring_scheme.score(alphabet1, *std::ranges::begin(sequence2))),
*cell_index_column_it);
// ---------------------------------------------------------------------
// Iteration phase: iterate over column and compute each cell
// ---------------------------------------------------------------------
for (auto && alphabet2 : std::views::drop(sequence2, 1))
{
current_alignment_column_it = next_alignment_column_it;
auto cell = *current_alignment_column_it;
cell = this->track_cell(this->compute_inner_cell(cell.best_score(),
*++next_alignment_column_it,
this->scoring_scheme.score(alphabet1, alphabet2)),
*++cell_index_column_it);
}
// ---------------------------------------------------------------------
// Final phase: track last cell
// ---------------------------------------------------------------------
this->track_last_row_cell(*current_alignment_column_it, *cell_index_column_it);
}
};
} // namespace seqan3::detail
seqan3::align_cfg::scoring_scheme::scoring_scheme
scoring_scheme(scheme_t) -> scoring_scheme< std::remove_cvref_t< scheme_t > >
Deduces the scoring scheme type from the constructor argument.
std::forward
T forward(T... args)
seqan3::views::slice
constexpr auto slice
A view adaptor that returns a half-open interval on the underlying range.
Definition slice.hpp:137
seqan3::views::elements
constexpr auto elements
A view calling get on each element in a range.
Definition elements.hpp:77
std::max
T max(T... args)
std::min
T min(T... args)
seqan3::get
constexpr auto const & get(configuration< configs_t... > const &config) noexcept
This is an overloaded member function, provided for convenience. It differs from the above function o...
Definition configuration.hpp:412
pairwise_alignment_algorithm.hpp
Provides seqan3::detail::pairwise_alignment_algorithm.
slice.hpp
Provides seqan3::views::slice.
Hide me