SeqAn3  3.1.0
The Modern C++ library for sequence analysis.
pairwise_alignment_algorithm_banded.hpp
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------------------------------
2 // Copyright (c) 2006-2021, Knut Reinert & Freie Universität Berlin
3 // Copyright (c) 2016-2021, Knut Reinert & MPI für molekulare Genetik
4 // This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
5 // shipped with this file and also available at: https://github.com/seqan/seqan3/blob/master/LICENSE.md
6 // -----------------------------------------------------------------------------------------------------
7 
13 #pragma once
14 
15 #include <seqan3/std/concepts>
16 #include <seqan3/std/ranges>
17 
20 
21 namespace seqan3::detail
22 {
23 
29 template <typename alignment_configuration_t, typename ...policies_t>
31  requires is_type_specialisation_of_v<alignment_configuration_t, configuration>
33 class pairwise_alignment_algorithm_banded :
34  protected pairwise_alignment_algorithm<alignment_configuration_t, policies_t...>
35 {
36 protected:
38  using base_algorithm_t = pairwise_alignment_algorithm<alignment_configuration_t, policies_t...>;
39 
40  // Import types from base class.
41  using typename base_algorithm_t::traits_type;
42  using typename base_algorithm_t::alignment_result_type;
43  using typename base_algorithm_t::score_type;
44 
45  static_assert(!std::same_as<alignment_result_type, empty_type>, "Alignment result type was not configured.");
46  static_assert(traits_type::is_banded, "Alignment configuration must have band configured.");
47 
48 public:
52  pairwise_alignment_algorithm_banded() = default;
53  pairwise_alignment_algorithm_banded(pairwise_alignment_algorithm_banded const &) = default;
54  pairwise_alignment_algorithm_banded(pairwise_alignment_algorithm_banded &&) = default;
55  pairwise_alignment_algorithm_banded & operator=(pairwise_alignment_algorithm_banded const &) = default;
56  pairwise_alignment_algorithm_banded & operator=(pairwise_alignment_algorithm_banded &&) = default;
57  ~pairwise_alignment_algorithm_banded() = default;
58 
68  pairwise_alignment_algorithm_banded(alignment_configuration_t const & config) : base_algorithm_t(config)
69  {}
71 
76  template <indexed_sequence_pair_range indexed_sequence_pairs_t, typename callback_t>
78  requires std::invocable<callback_t, alignment_result_type>
80  void operator()(indexed_sequence_pairs_t && indexed_sequence_pairs, callback_t && callback)
81  {
82  using std::get;
83 
84  for (auto && [sequence_pair, idx] : indexed_sequence_pairs)
85  {
86  size_t sequence1_size = std::ranges::distance(get<0>(sequence_pair));
87  size_t const sequence2_size = std::ranges::distance(get<1>(sequence_pair));
88 
89  auto && [alignment_matrix, index_matrix] = this->acquire_matrices(sequence1_size,
90  sequence2_size,
91  this->lowest_viable_score());
92 
93  // Initialise the cell updater with the dimensions of the regular matrix.
94  this->compare_and_set_optimum.set_target_indices(row_index_type{sequence2_size},
95  column_index_type{sequence1_size});
96 
97  // Shrink the first sequence if the band ends before its actual end.
98  sequence1_size = std::min(sequence1_size, this->upper_diagonal + sequence2_size);
99 
100  using sequence1_difference_t = std::ranges::range_difference_t<decltype(get<0>(sequence_pair))>;
101 
102  compute_matrix(std::views::take(get<0>(sequence_pair), static_cast<sequence1_difference_t>(sequence1_size)),
103  get<1>(sequence_pair),
104  alignment_matrix,
105  index_matrix);
106  this->make_result_and_invoke(std::forward<decltype(sequence_pair)>(sequence_pair),
107  std::move(idx),
108  this->optimal_score,
109  this->optimal_coordinate,
110  alignment_matrix,
111  callback);
112  }
113  }
114 
116  template <indexed_sequence_pair_range indexed_sequence_pairs_t, typename callback_t>
118  requires traits_type::is_vectorised && std::invocable<callback_t, alignment_result_type>
120  auto operator()(indexed_sequence_pairs_t && indexed_sequence_pairs, callback_t && callback)
121  {
122  using simd_collection_t = std::vector<score_type, aligned_allocator<score_type, alignof(score_type)>>;
123  using original_score_t = typename traits_type::original_score_type;
124 
125  // Extract the batch of sequences for the first and the second sequence.
126  auto seq1_collection = indexed_sequence_pairs | views::elements<0> | views::elements<0>;
127  auto seq2_collection = indexed_sequence_pairs | views::elements<0> | views::elements<1>;
128 
129  this->initialise_tracker(seq1_collection, seq2_collection);
130 
131  // Convert batch of sequences to sequence of simd vectors.
132  thread_local simd_collection_t simd_seq1_collection{};
133  thread_local simd_collection_t simd_seq2_collection{};
134 
135  this->convert_batch_of_sequences_to_simd_vector(simd_seq1_collection,
136  seq1_collection,
137  this->scoring_scheme.padding_symbol);
138  this->convert_batch_of_sequences_to_simd_vector(simd_seq2_collection,
139  seq2_collection,
140  this->scoring_scheme.padding_symbol);
141 
142  size_t const sequence1_size = std::ranges::distance(simd_seq1_collection);
143  size_t const sequence2_size = std::ranges::distance(simd_seq2_collection);
144 
145  auto && [alignment_matrix, index_matrix] = this->acquire_matrices(sequence1_size,
146  sequence2_size,
147  this->lowest_viable_score());
148 
149  compute_matrix(simd_seq1_collection, simd_seq2_collection, alignment_matrix, index_matrix);
150 
151  size_t index = 0;
152  for (auto && [sequence_pair, idx] : indexed_sequence_pairs)
153  {
154  original_score_t score = this->optimal_score[index] -
155  (this->padding_offsets[index] * this->scoring_scheme.padding_match_score());
156  matrix_coordinate coordinate{row_index_type{size_t{this->optimal_coordinate.row[index]}},
157  column_index_type{size_t{this->optimal_coordinate.col[index]}}};
158  this->make_result_and_invoke(std::forward<decltype(sequence_pair)>(sequence_pair),
159  std::move(idx),
160  std::move(score),
161  std::move(coordinate),
162  alignment_matrix,
163  callback);
164  ++index;
165  }
166  }
168 
169 protected:
222  template <std::ranges::forward_range sequence1_t,
223  std::ranges::forward_range sequence2_t,
224  std::ranges::input_range alignment_matrix_t,
225  std::ranges::input_range index_matrix_t>
227  requires std::ranges::forward_range<std::ranges::range_reference_t<alignment_matrix_t>> &&
228  std::ranges::forward_range<std::ranges::range_reference_t<index_matrix_t>>
230  void compute_matrix(sequence1_t && sequence1,
231  sequence2_t && sequence2,
232  alignment_matrix_t && alignment_matrix,
233  index_matrix_t && index_matrix)
234  {
235  // ---------------------------------------------------------------------
236  // Initialisation phase: allocate memory and initialise first column.
237  // ---------------------------------------------------------------------
238 
239  this->reset_optimum(); // Reset the tracker for the new alignment computation.
240 
241  auto alignment_matrix_it = alignment_matrix.begin();
242  auto indexed_matrix_it = index_matrix.begin();
243 
244  using row_index_t = std::ranges::range_difference_t<sequence2_t>; // row_size = |sequence2| + 1
245  using column_index_t = std::ranges::range_difference_t<sequence1_t>; // column_size = |sequence1| + 1
246 
247  row_index_t row_size = std::max<int32_t>(0, -this->lower_diagonal);
248  column_index_t const column_size = std::max<int32_t>(0, this->upper_diagonal);
249  this->initialise_column(*alignment_matrix_it, *indexed_matrix_it, std::views::take(sequence2, row_size));
250 
251  // ---------------------------------------------------------------------
252  // 1st recursion phase: band intersects with the first row.
253  // ---------------------------------------------------------------------
254 
255  for (auto alphabet1 : std::views::take(sequence1, column_size))
256  {
257  this->compute_column(*++alignment_matrix_it,
258  *++indexed_matrix_it,
259  alphabet1,
260  std::views::take(sequence2, ++row_size));
261  }
262 
263  // ---------------------------------------------------------------------
264  // 2nd recursion phase: iterate until the end of the matrix.
265  // ---------------------------------------------------------------------
266 
267  row_index_t first_row_index = 0u;
268  for (auto alphabet1 : std::views::drop(sequence1, column_size))
269  {
270  compute_band_column(*++alignment_matrix_it,
271  std::views::drop(*++indexed_matrix_it, first_row_index + 1),
272  alphabet1,
273  views::slice(sequence2, first_row_index, ++row_size));
274  ++first_row_index;
275  }
276 
277  // ---------------------------------------------------------------------
278  // Final phase: track score of last column
279  // ---------------------------------------------------------------------
280 
281  auto alignment_column = *alignment_matrix_it;
282  auto cell_index_column = std::views::drop(*indexed_matrix_it, first_row_index);
283 
284  auto alignment_column_it = alignment_column.begin();
285  auto cell_index_column_it = cell_index_column.begin();
286 
287  this->track_last_column_cell(*alignment_column_it, *cell_index_column_it);
288 
289  for (row_index_t last_row = std::min<row_index_t>(std::ranges::distance(sequence2), row_size);
290  first_row_index < last_row;
291  ++first_row_index)
292  this->track_last_column_cell(*++alignment_column_it, *++cell_index_column_it);
293 
294  this->track_final_cell(*alignment_column_it, *cell_index_column_it);
295  }
296 
348  template <std::ranges::forward_range alignment_column_t,
349  std::ranges::input_range cell_index_column_t,
350  typename alphabet1_t,
351  std::ranges::input_range sequence2_t>
352  void compute_band_column(alignment_column_t && alignment_column,
353  cell_index_column_t && cell_index_column,
354  alphabet1_t const & alphabet1,
355  sequence2_t && sequence2)
356  {
357  // ---------------------------------------------------------------------
358  // Initial phase: prepare column and initialise first cell
359  // ---------------------------------------------------------------------
360 
361  auto current_alignment_column_it = alignment_column.begin();
362  auto cell_index_column_it = cell_index_column.begin();
363 
364  // Points to the last valid cell in the column.
365  decltype(current_alignment_column_it) next_alignment_column_it{current_alignment_column_it};
366  auto cell = *current_alignment_column_it;
367  cell = this->track_cell(
368  this->initialise_band_first_cell(cell.best_score(),
369  *++next_alignment_column_it,
370  this->scoring_scheme.score(alphabet1, *std::ranges::begin(sequence2))),
371  *cell_index_column_it);
372 
373  // ---------------------------------------------------------------------
374  // Iteration phase: iterate over column and compute each cell
375  // ---------------------------------------------------------------------
376 
377  for (auto && alphabet2 : sequence2 | std::views::drop(1))
378  {
379  current_alignment_column_it = next_alignment_column_it;
380  auto cell = *current_alignment_column_it;
381  cell = this->track_cell(
382  this->compute_inner_cell(cell.best_score(),
383  *++next_alignment_column_it,
384  this->scoring_scheme.score(alphabet1, alphabet2)),
385  *++cell_index_column_it);
386  }
387 
388  // ---------------------------------------------------------------------
389  // Final phase: track last cell
390  // ---------------------------------------------------------------------
391 
392  this->track_last_row_cell(*current_alignment_column_it, *cell_index_column_it);
393  }
394 };
395 
396 } // namespace seqan3::detail
The <concepts> header from C++20's standard library.
T forward(T... args)
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: traits.hpp:368
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: traits.hpp:388
constexpr auto slice
A view adaptor that returns a half-open interval on the underlying range.
Definition: slice.hpp:183
T min(T... args)
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:429
Provides seqan3::detail::pairwise_alignment_algorithm.
The <ranges> header from C++20's standard library.
Provides seqan3::views::slice.