Loading...
Searching...
No Matches
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
align_config_band.hpp
Provides seqan3::detail::align_config_band.
align_config_scoring_scheme.hpp
Provides seqan3::align_cfg::scoring_scheme.
align_result_selector.hpp
Provides seqan3::detail::align_result_selector.
aligned_allocator.hpp
Provides seqan3::aligned_allocator.
aligned_sequence_builder.hpp
Provides seqan3::detail::aligned_sequence_builder.
exception.hpp
Includes customized exception types for the alignment module .
concept.hpp
Provides concepts needed internally for the alignment algorithms.
type_traits.hpp
Provides helper type traits for the configuration and execution of the alignment algorithm.
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::declval
T declval(T... args)
deferred_crtp_base.hpp
Provides seqan3::detail::deferred_crtp_base.
elements.hpp
Provides seqan3::views::elements.
empty_type.hpp
Provides seqan3::detail::empty_type.
function_traits.hpp
Provides various type traits for use on functions.
seqan3::field::offset
@ offset
Sequence (seqan3::field::seq) relative start position (0-based), unsigned value.
Definition record.hpp:66
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::is_signed_v
T is_signed_v
iterator
std::make_shared
T make_shared(T... args)
memory
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
optional
std::vector::push_back
T push_back(T... args)
trace_directions.hpp
Provides the declaration of seqan3::detail::trace_directions.
type_traits
Hide me