Loading...
Searching...
No Matches
pairwise_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 <concepts>
#include <ranges>
#include <seqan3/alignment/pairwise/detail/type_traits.hpp>
#include <seqan3/core/detail/empty_type.hpp>
#include <seqan3/utility/container/aligned_allocator.hpp>
#include <seqan3/utility/detail/type_name_as_string.hpp>
#include <seqan3/utility/simd/views/to_simd.hpp>
#include <seqan3/utility/views/elements.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 : protected policies_t...
{
protected:
using traits_type = alignment_configuration_traits<alignment_configuration_t>;
using score_type = typename traits_type::score_type;
using alignment_result_type = typename traits_type::alignment_result_type;
static_assert(!std::same_as<alignment_result_type, empty_type>, "Alignment result type was not configured.");
public:
pairwise_alignment_algorithm() = default;
pairwise_alignment_algorithm(pairwise_alignment_algorithm const &) = default;
pairwise_alignment_algorithm(pairwise_alignment_algorithm &&) = default;
pairwise_alignment_algorithm & operator=(pairwise_alignment_algorithm const &) = default;
pairwise_alignment_algorithm & operator=(pairwise_alignment_algorithm &&) = default;
~pairwise_alignment_algorithm() = default;
pairwise_alignment_algorithm(alignment_configuration_t const & config) : policies_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 const 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);
compute_matrix(get<0>(sequence_pair), 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{};
convert_batch_of_sequences_to_simd_vector(simd_seq1_collection,
seq1_collection,
this->scoring_scheme.padding_symbol);
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);
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 <typename simd_sequence_t, std::ranges::forward_range sequence_collection_t, arithmetic padding_symbol_t>
requires std::ranges::output_range<simd_sequence_t, score_type>
void convert_batch_of_sequences_to_simd_vector(simd_sequence_t & simd_sequence,
sequence_collection_t & sequences,
padding_symbol_t const & padding_symbol)
{
assert(static_cast<size_t>(std::ranges::distance(sequences)) <= traits_type::alignments_per_vector);
simd_sequence.clear();
for (auto && simd_vector_chunk : sequences | views::to_simd<score_type>(padding_symbol))
std::ranges::move(simd_vector_chunk, std::back_inserter(simd_sequence));
}
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();
initialise_column(*alignment_matrix_it, *indexed_matrix_it, sequence2);
// ---------------------------------------------------------------------
// Iteration phase: compute column-wise the alignment matrix.
// ---------------------------------------------------------------------
for (auto alphabet1 : sequence1)
compute_column(*++alignment_matrix_it,
*++indexed_matrix_it,
this->scoring_scheme_profile_column(alphabet1),
sequence2);
// ---------------------------------------------------------------------
// Final phase: track score of last column
// ---------------------------------------------------------------------
auto && alignment_column = *alignment_matrix_it;
auto && cell_index_column = *indexed_matrix_it;
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 ([[maybe_unused]] auto && unused : sequence2)
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::input_range alignment_column_t,
std::ranges::input_range cell_index_column_t,
std::ranges::input_range sequence2_t>
void initialise_column(alignment_column_t && alignment_column,
cell_index_column_t && cell_index_column,
sequence2_t && sequence2)
{
// ---------------------------------------------------------------------
// Initial phase: prepare column and initialise first cell
// ---------------------------------------------------------------------
auto first_column_it = alignment_column.begin();
auto cell_index_column_it = cell_index_column.begin();
*first_column_it = this->track_cell(this->initialise_origin_cell(), *cell_index_column_it);
// ---------------------------------------------------------------------
// Iteration phase: iterate over column and compute each cell
// ---------------------------------------------------------------------
for ([[maybe_unused]] auto const & unused : sequence2)
{
++first_column_it;
*first_column_it =
this->track_cell(this->initialise_first_column_cell(*first_column_it), *++cell_index_column_it);
}
// ---------------------------------------------------------------------
// Final phase: track last cell of initial column
// ---------------------------------------------------------------------
this->track_last_row_cell(*first_column_it, *cell_index_column_it);
}
template <std::ranges::input_range alignment_column_t,
std::ranges::input_range cell_index_column_t,
typename alphabet1_t,
std::ranges::input_range sequence2_t>
requires semialphabet<alphabet1_t> || simd_concept<alphabet1_t>
void compute_column(alignment_column_t && alignment_column,
cell_index_column_t && cell_index_column,
alphabet1_t const & alphabet1,
sequence2_t && sequence2)
{
using score_type = typename traits_type::score_type;
// ---------------------------------------------------------------------
// Initial phase: prepare column and initialise first cell
// ---------------------------------------------------------------------
auto alignment_column_it = alignment_column.begin();
auto cell_index_column_it = cell_index_column.begin();
auto cell = *alignment_column_it;
score_type diagonal = cell.best_score();
*alignment_column_it = this->track_cell(this->initialise_first_row_cell(cell), *cell_index_column_it);
// ---------------------------------------------------------------------
// Iteration phase: iterate over column and compute each cell
// ---------------------------------------------------------------------
for (auto const & alphabet2 : sequence2)
{
auto cell = *++alignment_column_it;
score_type next_diagonal = cell.best_score();
*alignment_column_it = this->track_cell(
this->compute_inner_cell(diagonal, cell, this->scoring_scheme.score(alphabet1, alphabet2)),
*++cell_index_column_it);
diagonal = next_diagonal;
}
// ---------------------------------------------------------------------
// Final phase: track last cell
// ---------------------------------------------------------------------
this->track_last_row_cell(*alignment_column_it, *cell_index_column_it);
}
};
} // namespace seqan3::detail
aligned_allocator.hpp
Provides seqan3::aligned_allocator.
type_traits.hpp
Provides helper type traits for the configuration and execution of the alignment algorithm.
std::back_inserter
T back_inserter(T... args)
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.
elements.hpp
Provides seqan3::views::elements.
empty_type.hpp
Provides seqan3::detail::empty_type.
std::forward
T forward(T... args)
seqan3::views::elements
constexpr auto elements
A view calling get on each element in a range.
Definition elements.hpp:77
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
type_name_as_string.hpp
Provides traits to inspect some information of a type, for example its name.
Hide me