29#include <seqan3/utility/simd/concept.hpp>
30#include <seqan3/utility/simd/simd.hpp>
31#include <seqan3/utility/simd/simd_traits.hpp>
32#include <seqan3/utility/simd/views/to_simd.hpp>
36namespace seqan3::detail
72template <
typename config_t,
typename... algorithm_policies_t>
73class alignment_algorithm :
74 public invoke_deferred_crtp_base<algorithm_policies_t, alignment_algorithm<config_t, algorithm_policies_t...>>...
78 using traits_t = alignment_configuration_traits<config_t>;
90 template <
typename alignment_algorithm_t = alignment_algorithm>
91 static auto _alignment_column_t() ->
decltype(std::declval<alignment_algorithm_t>().current_alignment_column());
94 using alignment_column_t =
decltype(_alignment_column_t());
96 using alignment_column_iterator_t = std::ranges::iterator_t<alignment_column_t>;
98 using alignment_result_t =
typename traits_t::alignment_result_type;
100 static_assert(!std::same_as<alignment_result_t, empty_type>,
"Alignment result type was not configured.");
103 using score_debug_matrix_t =
105 two_dimensional_matrix<std::optional<typename traits_t::original_score_type>,
107 matrix_major_order::column>,
110 using trace_debug_matrix_t =
112 two_dimensional_matrix<std::optional<trace_directions>,
114 matrix_major_order::column>,
121 constexpr alignment_algorithm() =
default;
122 constexpr alignment_algorithm(alignment_algorithm
const &) =
default;
123 constexpr alignment_algorithm(alignment_algorithm &&) =
default;
124 constexpr alignment_algorithm & operator=(alignment_algorithm
const &) =
default;
125 constexpr alignment_algorithm & operator=(alignment_algorithm &&) =
default;
126 ~alignment_algorithm() =
default;
136 explicit constexpr alignment_algorithm(config_t
const & cfg) :
137 invoke_deferred_crtp_base<algorithm_policies_t, alignment_algorithm<config_t, algorithm_policies_t...>>{cfg}...,
138 cfg_ptr{std::make_shared<config_t>(cfg)}
140 this->scoring_scheme = seqan3::get<align_cfg::scoring_scheme>(*cfg_ptr).scheme;
141 this->initialise_alignment_state(*cfg_ptr);
190 template <indexed_sequence_pair_range indexed_sequence_pairs_t,
typename callback_t>
191 requires (!traits_t::is_vectorised) && std::invocable<callback_t, alignment_result_t>
192 void operator()(indexed_sequence_pairs_t && indexed_sequence_pairs, callback_t && callback)
196 for (
auto && [sequence_pair, idx] : indexed_sequence_pairs)
197 compute_single_pair(idx,
get<0>(sequence_pair),
get<1>(sequence_pair), callback);
201 template <indexed_sequence_pair_range indexed_sequence_pairs_t,
typename callback_t>
202 requires traits_t::is_vectorised && std::invocable<callback_t, alignment_result_t>
203 void operator()(indexed_sequence_pairs_t && indexed_sequence_pairs, callback_t && callback)
205 assert(cfg_ptr !=
nullptr);
207 static_assert(simd_concept<typename traits_t::score_type>,
"Expected simd score type.");
208 static_assert(simd_concept<typename traits_t::trace_type>,
"Expected simd trace type.");
211 auto sequence1_range = indexed_sequence_pairs | views::elements<0> | views::elements<0>;
212 auto sequence2_range = indexed_sequence_pairs | views::elements<0> | views::elements<1>;
215 this->initialise_find_optimum_policy(sequence1_range,
217 this->scoring_scheme.padding_match_score());
220 auto simd_sequences1 = convert_batch_of_sequences_to_simd_vector(sequence1_range);
221 auto simd_sequences2 = convert_batch_of_sequences_to_simd_vector(sequence2_range);
223 max_size_in_collection =
std::pair{simd_sequences1.size(), simd_sequences2.size()};
225 this->alignment_state.reset_optimum();
227 compute_matrix(simd_sequences1, simd_sequences2);
229 make_alignment_result(indexed_sequence_pairs, callback);
247 template <
typename sequence_range_t>
248 constexpr auto convert_batch_of_sequences_to_simd_vector(sequence_range_t & sequences)
250 assert(
static_cast<size_t>(std::ranges::distance(sequences)) <= traits_t::alignments_per_vector);
252 using simd_score_t =
typename traits_t::score_type;
254 std::vector<simd_score_t, aligned_allocator<simd_score_t,
alignof(simd_score_t)>> simd_sequence{};
256 for (
auto && simd_vector_chunk : sequences | views::to_simd<simd_score_t>(this->scoring_scheme.padding_symbol))
257 for (auto && simd_vector : simd_vector_chunk)
258 simd_sequence.push_back(
std::
move(simd_vector));
260 return simd_sequence;
280 template <std::ranges::forward_range sequence1_t, std::ranges::forward_range sequence2_t,
typename callback_t>
282 compute_single_pair(
size_t const idx, sequence1_t && sequence1, sequence2_t && sequence2, callback_t & callback)
284 assert(cfg_ptr !=
nullptr);
286 if constexpr (traits_t::is_debug)
287 initialise_debug_matrices(sequence1, sequence2);
290 this->alignment_state.reset_optimum();
292 if constexpr (traits_t::is_banded)
296 auto const & band = get<align_cfg::band_fixed_size>(*cfg_ptr);
297 check_valid_band_parameter(sequence1, sequence2, band);
298 auto && [subsequence1, subsequence2] = this->slice_sequences(sequence1, sequence2, band);
300 compute_matrix(subsequence1, subsequence2, band);
301 make_alignment_result(idx, subsequence1, subsequence2, callback);
305 compute_matrix(sequence1, sequence2);
306 make_alignment_result(idx, sequence1, sequence2, callback);
326 template <
typename sequence1_t,
typename sequence2_t>
327 constexpr void check_valid_band_parameter(sequence1_t && sequence1,
328 sequence2_t && sequence2,
329 align_cfg::band_fixed_size
const & band)
331 static_assert(config_t::template exists<align_cfg::band_fixed_size>(),
332 "The band configuration is required for the banded alignment algorithm.");
335 static_assert(std::is_signed_v<diff_type>,
"Only signed types can be used to test the band parameters.");
337 if (
static_cast<diff_type
>(band.lower_diagonal) > std::ranges::distance(sequence1))
339 throw invalid_alignment_configuration{
340 "Invalid band error: The lower diagonal excludes the whole alignment matrix."};
343 if (
static_cast<diff_type
>(band.upper_diagonal) < -std::ranges::distance(sequence2))
345 throw invalid_alignment_configuration{
346 "Invalid band error: The upper diagonal excludes the whole alignment matrix."};
362 template <
typename sequence1_t,
typename sequence2_t>
363 constexpr void initialise_debug_matrices(sequence1_t & sequence1, sequence2_t & sequence2)
365 size_t rows = std::ranges::distance(sequence2) + 1;
366 size_t cols = std::ranges::distance(sequence1) + 1;
368 score_debug_matrix = score_debug_matrix_t{number_rows{rows}, number_cols{cols}};
369 trace_debug_matrix = trace_debug_matrix_t{number_rows{rows}, number_cols{cols}};
379 template <
typename sequence1_t,
typename sequence2_t>
380 void compute_matrix(sequence1_t & sequence1, sequence2_t & sequence2)
381 requires (!traits_t::is_banded)
387 this->allocate_matrix(sequence1, sequence2);
388 initialise_first_alignment_column(sequence2);
394 for (
auto const & alphabet1 : sequence1)
396 compute_alignment_column<true>(this->scoring_scheme_profile_column(alphabet1), sequence2);
397 finalise_last_cell_in_column(
true);
404 finalise_alignment();
408 template <
typename sequence1_t,
typename sequence2_t>
409 void compute_matrix(sequence1_t & sequence1, sequence2_t & sequence2, align_cfg::band_fixed_size
const & band)
410 requires (traits_t::is_banded)
417 this->allocate_matrix(sequence1, sequence2, band, this->alignment_state);
418 using row_index_t = std::ranges::range_difference_t<sequence2_t>;
419 row_index_t last_row_index = this->score_matrix.band_row_index;
420 initialise_first_alignment_column(std::views::take(sequence2, last_row_index));
426 row_index_t sequence2_size = std::ranges::distance(sequence2);
427 for (
auto const & seq1_value :
std::views::
take(sequence1, this->score_matrix.band_col_index))
429 compute_alignment_column<true>(seq1_value, std::views::take(sequence2, ++last_row_index));
431 finalise_last_cell_in_column(last_row_index >= sequence2_size);
438 size_t first_row_index = 0;
439 for (
auto const & seq1_value :
std::views::
drop(sequence1, this->score_matrix.band_col_index))
442 compute_alignment_column<false>(seq1_value, sequence2 |
views::slice(first_row_index++, ++last_row_index));
444 finalise_last_cell_in_column(last_row_index >= sequence2_size);
451 finalise_alignment();
466 template <
typename sequence2_t>
467 auto initialise_first_alignment_column(sequence2_t && sequence2)
470 alignment_column = this->current_alignment_column();
471 assert(!alignment_column.empty());
474 alignment_column_it = alignment_column.begin();
475 this->init_origin_cell(*alignment_column_it, this->alignment_state);
478 for (
auto it =
std::ranges::begin(sequence2); it != std::ranges::end(sequence2); ++it)
479 this->init_column_cell(*++alignment_column_it, this->alignment_state);
482 bool at_last_row =
true;
483 if constexpr (traits_t::is_banded)
484 at_last_row =
static_cast<size_t>(this->score_matrix.band_row_index) == this->score_matrix.num_rows - 1;
486 finalise_last_cell_in_column(at_last_row);
504 template <
bool initialise_first_cell,
typename sequence1_value_t,
typename sequence2_t>
505 void compute_alignment_column(sequence1_value_t
const & seq1_value, sequence2_t && sequence2)
507 this->next_alignment_column();
508 alignment_column = this->current_alignment_column();
509 alignment_column_it = alignment_column.begin();
513 if constexpr (initialise_first_cell)
515 this->init_row_cell(*alignment_column_it, this->alignment_state);
519 this->compute_first_band_cell(*alignment_column_it,
520 this->alignment_state,
521 this->scoring_scheme.score(seq1_value, *seq2_it));
525 for (; seq2_it != std::ranges::end(sequence2); ++seq2_it)
526 this->compute_cell(*++alignment_column_it,
527 this->alignment_state,
528 this->scoring_scheme.score(seq1_value, *seq2_it));
541 constexpr void finalise_last_cell_in_column(
bool const at_last_row)
noexcept
544 this->check_score_of_last_row_cell(*alignment_column_it, this->alignment_state);
546 if constexpr (traits_t::is_debug)
547 dump_alignment_column();
551 constexpr void finalise_alignment() noexcept
557 this->check_score_of_cells_in_last_column(alignment_column, this->alignment_state);
558 this->check_score_of_last_cell(*alignment_column_it, this->alignment_state);
587 template <
typename index_t,
typename sequence1_t,
typename sequence2_t,
typename callback_t>
588 requires (!traits_t::is_vectorised)
589 constexpr void make_alignment_result([[maybe_unused]] index_t
const idx,
590 [[maybe_unused]] sequence1_t & sequence1,
591 [[maybe_unused]] sequence2_t & sequence2,
592 callback_t & callback)
594 using result_value_t =
typename alignment_result_value_type_accessor<alignment_result_t>::type;
600 static_assert(seqan3::detail::alignment_configuration_traits<config_t>::has_output_configuration,
601 "The configuration must contain at least one align_cfg::output_* element.");
603 result_value_t res{};
605 if constexpr (traits_t::output_sequence1_id)
606 res.sequence1_id = idx;
608 if constexpr (traits_t::output_sequence2_id)
609 res.sequence2_id = idx;
612 if constexpr (traits_t::compute_score)
613 res.score = this->alignment_state.optimum.score;
615 if constexpr (traits_t::compute_end_positions)
617 using alignment_coordinate_t = detail::advanceable_alignment_coordinate<>;
618 res.end_positions = alignment_coordinate_t{column_index_type{this->alignment_state.optimum.column_index},
619 row_index_type{this->alignment_state.optimum.row_index}};
621 if constexpr (traits_t::is_banded)
622 res.end_positions.second += res.end_positions.first - this->trace_matrix.band_col_index;
625 if constexpr (traits_t::compute_begin_positions)
628 aligned_sequence_builder builder{sequence1, sequence2};
630 detail::matrix_coordinate
const optimum_coordinate{
631 detail::row_index_type{this->alignment_state.optimum.row_index},
632 detail::column_index_type{this->alignment_state.optimum.column_index}};
633 auto trace_res = builder(this->trace_matrix.trace_path(optimum_coordinate));
634 res.begin_positions.first = trace_res.first_sequence_slice_positions.first;
635 res.begin_positions.second = trace_res.second_sequence_slice_positions.first;
637 if constexpr (traits_t::compute_sequence_alignment)
638 res.alignment = std::move(trace_res.alignment);
642 if constexpr (traits_t::is_debug)
644 res.score_debug_matrix = std::move(score_debug_matrix);
645 if constexpr (traits_t::compute_sequence_alignment)
646 res.trace_debug_matrix = std::move(trace_debug_matrix);
649 callback(std::move(res));
677 template <
typename indexed_sequence_pair_range_t,
typename callback_t>
678 requires traits_t::is_vectorised
679 constexpr auto make_alignment_result(indexed_sequence_pair_range_t && index_sequence_pairs, callback_t & callback)
681 using result_value_t =
typename alignment_result_value_type_accessor<alignment_result_t>::type;
683 size_t simd_index = 0;
684 for (
auto && [sequence_pairs, alignment_index] : index_sequence_pairs)
686 (void)sequence_pairs;
687 result_value_t res{};
689 if constexpr (traits_t::output_sequence1_id)
690 res.sequence1_id = alignment_index;
692 if constexpr (traits_t::output_sequence2_id)
693 res.sequence2_id = alignment_index;
695 if constexpr (traits_t::compute_score)
696 res.score = this->alignment_state.optimum.score[simd_index];
698 if constexpr (traits_t::compute_end_positions)
700 res.end_positions.first = this->alignment_state.optimum.column_index[simd_index];
701 res.end_positions.second = this->alignment_state.optimum.row_index[simd_index];
704 callback(std::move(res));
717 void dump_alignment_column()
721 auto column = this->current_alignment_column();
723 auto coord = get<1>(column.front()).coordinate;
724 if constexpr (traits_t::is_banded)
725 coord.second += coord.first - this->score_matrix.band_col_index;
731 | std::views::transform(
735 return get<0>(tpl).current;
737 score_debug_matrix.begin() +
offset);
740 if constexpr (traits_t::compute_sequence_alignment)
744 | std::views::transform(
748 auto trace = get<1>(tpl).current;
750 if (
auto _up = (trace & trace_directions::up_open); _up == trace_directions::carry_up_open)
751 trace = trace ^ trace_directions::carry_up_open;
752 else if (_up == trace_directions::up_open)
753 trace = trace ^ trace_directions::up;
755 if (
auto _left = (trace & trace_directions::left_open);
756 _left == trace_directions::carry_left_open)
757 trace = trace ^ trace_directions::carry_left_open;
758 else if (_left == trace_directions::left_open)
759 trace = trace ^ trace_directions::left;
763 trace_debug_matrix.begin() +
offset);
770 alignment_column_t alignment_column{};
772 alignment_column_iterator_t alignment_column_it{};
774 score_debug_matrix_t score_debug_matrix{};
776 trace_debug_matrix_t trace_debug_matrix{};
Provides seqan3::detail::align_config_band.
Provides seqan3::align_cfg::scoring_scheme.
Provides seqan3::detail::align_result_selector.
Provides seqan3::aligned_allocator.
Provides seqan3::detail::aligned_sequence_builder.
Includes customized exception types for the alignment module .
Provides concepts needed internally for the alignment algorithms.
Provides helper type traits for the configuration and execution of the alignment algorithm.
Provides seqan3::detail::deferred_crtp_base.
Provides seqan3::views::elements.
Provides seqan3::detail::empty_type.
Provides various type traits for use on functions.
@ offset
Sequence (seqan3::field::seq) relative start position (0-based), unsigned value.
typename decltype(detail::split_after< i >(list_t{}))::first_type take
Return a seqan3::type_list of the first n types in the input type list.
Definition type_list/traits.hpp:374
typename decltype(detail::split_after< i >(list_t{}))::second_type drop
Return a seqan3::type_list of the types in the input type list, except the first n.
Definition type_list/traits.hpp:392
constexpr auto slice
A view adaptor that returns a half-open interval on the underlying range.
Definition slice.hpp:175
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
SeqAn specific customisations in the standard namespace.
Provides the declaration of seqan3::detail::trace_directions.