SeqAn3
The Modern C++ library for sequence analysis.
banded_score_trace_dp_matrix_policy.hpp
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------------------------------
2 // Copyright (c) 2006-2019, Knut Reinert & Freie Universität Berlin
3 // Copyright (c) 2016-2019, 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 <deque>
16 
17 #include <range/v3/view/iota.hpp>
18 #include <range/v3/view/zip.hpp>
19 
23 #include <seqan3/std/span>
24 
25 namespace seqan3::detail
26 {
33 template <typename derived_t, typename score_allocator_t, typename trace_allocator_t>
34 class banded_score_trace_dp_matrix_policy :
35  public banded_score_dp_matrix_policy<banded_score_trace_dp_matrix_policy<derived_t,
36  score_allocator_t,
37  trace_allocator_t>,
38  score_allocator_t>
39 {
40 private:
41 
43  using base_t = banded_score_dp_matrix_policy<banded_score_trace_dp_matrix_policy<derived_t,
44  score_allocator_t,
45  trace_allocator_t>,
46  score_allocator_t>;
47 
49  friend derived_t;
50 
51  // Import members from base class.
52  using base_t::score_matrix;
53  using base_t::dimension_first_range;
54  using base_t::dimension_second_range;
55  using base_t::current_column_index;
56  using base_t::current_matrix_iter;
57  using base_t::band_column_index;
58  using base_t::band_row_index;
59  using base_t::band_size;
60 
61  // Import member functions from base class
62  using base_t::current_band_size;
63  using base_t::second_range_begin_offset;
64  using base_t::band_touches_last_row;
65  using base_t::trim_sequences;
66  using base_t::map_banded_coordinate_to_range_position;
67 
71  using cell_type = typename score_allocator_t::value_type;
74  using trace_type = typename trace_allocator_t::value_type;
76  using score_matrix_type = std::vector<cell_type, score_allocator_t>;
78  using trace_matrix_type = std::vector<trace_type, trace_allocator_t>;
80 
84  constexpr banded_score_trace_dp_matrix_policy() = default;
85  constexpr banded_score_trace_dp_matrix_policy(banded_score_trace_dp_matrix_policy const &) = default;
86  constexpr banded_score_trace_dp_matrix_policy(banded_score_trace_dp_matrix_policy &&) = default;
87  constexpr banded_score_trace_dp_matrix_policy & operator=(banded_score_trace_dp_matrix_policy const &) = default;
90  constexpr banded_score_trace_dp_matrix_policy & operator=(banded_score_trace_dp_matrix_policy &&) = default;
91  ~banded_score_trace_dp_matrix_policy() = default;
92 
102  template <typename first_range_t, typename second_range_t, typename band_t>
103  constexpr void allocate_matrix(first_range_t & first_range, second_range_t & second_range, band_t const & band)
104  {
105  base_t::allocate_matrix(first_range, second_range, band);
106 
107  trace_matrix.resize(band_size * dimension_first_range);
108  trace_matrix_iter = std::ranges::begin(trace_matrix) + band_column_index;
109  }
110 
112  constexpr auto current_column() noexcept
113  {
114  auto span = base_t::current_band_size();
115 
116  assert(span > 0u); // The span must always be greater than 0.
117 
118  // The front coordinate in the current column begins at it - begin(matrix).
119  // The back coordinate ends at it - begin(matrix) + current_band_size
120  advanceable_alignment_coordinate<advanceable_alignment_coordinate_state::row>
121  col_begin{column_index_type{current_column_index},
122  row_index_type{static_cast<size_t>(std::ranges::distance(std::ranges::begin(score_matrix),
123  current_matrix_iter))}};
124  advanceable_alignment_coordinate<advanceable_alignment_coordinate_state::row>
125  col_end{column_index_type{current_column_index}, row_index_type{col_begin.second + span}};
126 
127  // Return zip view over current column and current column shifted by one to access the previous horizontal.
128  auto zip_score = std::view::zip(std::span{std::addressof(*current_matrix_iter), span},
129  std::span{std::addressof(*(current_matrix_iter + 1)), span});
130  return std::view::zip(std::move(zip_score),
131  std::view::iota(col_begin, col_end),
132  std::span{std::addressof(*trace_matrix_iter), span});
133  }
134 
136  constexpr void go_next_column() noexcept
137  {
138  base_t::go_next_column();
139  // Move trace back pointer to begin of next column
140  trace_matrix_iter = std::ranges::begin(trace_matrix) + band_size * current_column_index;
141  // As long as we shift the band towards the band_column_index, jump to the correct begin.
142  if (current_column_index < band_column_index)
143  trace_matrix_iter += band_column_index - current_column_index;
144  }
145 
152  constexpr auto parse_traceback(alignment_coordinate const & back_coordinate) const
153  {
154  // Store the trace segments.
155  std::deque<gap_segment> first_segments{};
156  std::deque<gap_segment> second_segments{};
157 
158  // Put the iterator to the position where the traceback starts.
159  auto direction_iter = std::ranges::begin(trace_matrix);
160  std::ranges::advance(direction_iter, back_coordinate.first * band_size +
161  back_coordinate.second);
162 
163  // Parse the trace until interrupt.
164  while (*direction_iter != trace_directions::none)
165  {
166  // parse until end of diagonal run
167  while (static_cast<bool>(*direction_iter & trace_directions::diagonal))
168  {
169  std::ranges::advance(direction_iter, -static_cast<int_fast32_t>(band_size));
170  }
171 
172  size_t col_pos = std::ranges::distance(std::ranges::begin(trace_matrix), direction_iter) / band_size;
173  // parse vertical gap -> record gap in first_segments (will be translated into gap of first sequence)
174  if (static_cast<bool>(*direction_iter & trace_directions::up) ||
175  static_cast<bool>(*direction_iter & trace_directions::up_open))
176  {
177  // Get the current column index (note the column based layout)
178 
179  gap_segment gap{col_pos, 0u};
180 
181  // Follow gap until open signal is detected.
182  while (!static_cast<bool>(*direction_iter & trace_directions::up_open))
183  {
184  --direction_iter;
185  ++gap.size;
186  }
187  // explicitly follow opening gap
188  --direction_iter;
189  ++gap.size;
190  // record the gap
191  first_segments.push_front(std::move(gap));
192  continue;
193  }
194  // parse horizontal gap -> record gap in second_segments (will be translated into gap of second sequence)
195  if (static_cast<bool>(*direction_iter & trace_directions::left) ||
196  static_cast<bool>(*direction_iter & trace_directions::left_open))
197  {
198  // Get the current row index (note the column based layout)
199  size_t pos = std::ranges::distance(std::ranges::begin(trace_matrix), direction_iter) % band_size +
200  static_cast<int_fast32_t>(col_pos - band_column_index);
201  gap_segment gap{pos, 0u};
202 
203  // Follow gap until open signal is detected.
204  while (!static_cast<bool>(*direction_iter & trace_directions::left_open))
205  {
206  std::ranges::advance(direction_iter, -static_cast<int_fast32_t>(band_size) + 1);
207  ++gap.size;
208  }
209  // explicitly follow opening gap
210  std::ranges::advance(direction_iter, -static_cast<int_fast32_t>(band_size) + 1);
211  ++gap.size;
212  second_segments.push_front(std::move(gap));
213  }
214  }
215 
216  // Get front coordinate.
217  auto c = column_index_type{
218  static_cast<uint_fast32_t>(std::ranges::distance(std::ranges::begin(trace_matrix), direction_iter) /
219  band_size)};
220  auto r = row_index_type{
221  static_cast<uint_fast32_t>(std::ranges::distance(std::ranges::begin(trace_matrix), direction_iter) %
222  band_size)};
223 
224  // Validate correct coordinates.
225  auto front_coordinate = map_banded_coordinate_to_range_position(
226  alignment_coordinate{column_index_type{c}, row_index_type{r}});
227  assert(front_coordinate.first >= 0u);
228  assert(front_coordinate.first <= back_coordinate.first);
229 
230  assert(front_coordinate.second >= 0u);
231  assert(front_coordinate.second <= map_banded_coordinate_to_range_position(back_coordinate).second);
232 
233  return std::tuple{front_coordinate, first_segments, second_segments};
234  }
235 
237  constexpr void print_trace_matrix() const
238  {
239  auto printable = [](trace_directions dir)
240  {
241  std::string seq{};
242  if (dir == trace_directions::none)
243  seq.append("0");
244  if ((dir & trace_directions::diagonal) == trace_directions::diagonal)
245  seq.append("\\");
246  if ((dir & trace_directions::up) == trace_directions::up)
247  seq.append("|");
248  if ((dir & trace_directions::up_open) == trace_directions::up_open)
249  seq.append("^");
250  if ((dir & trace_directions::left) == trace_directions::left)
251  seq.append("-");
252  if ((dir & trace_directions::left_open) == trace_directions::left_open)
253  seq.append("<");
254  return seq;
255  };
256 
257  // First part: moving band right
258  for (size_t col = 0; col < band_column_index; ++col)
259  {
260  auto it = std::ranges::begin(trace_matrix) + (band_size * col) + (band_column_index - col);
261  for (size_t row = 0; row <= std::min(dimension_second_range - 1, band_row_index + col); ++row, ++it)
262  {
263  debug_stream << /*std::ranges::distance(std::ranges::begin(trace_matrix), it)*/printable(*it) << ',';
264  }
265  debug_stream << '\n';
266  }
267 
268  // Second part: moving band down.
269  for (size_t col = band_column_index; col < dimension_first_range; ++col)
270  {
271  for (size_t padding = 0; padding < col - band_column_index; ++padding)
272  debug_stream << " ,";
273 
274  auto it = std::ranges::begin(trace_matrix) + (band_size * col);
275  for (size_t row = 0; row < band_size; ++row, ++it)
276  {
277  // If the band moves out of the matrix do not try to print the characters.
278  if (col - band_column_index + row >= dimension_second_range)
279  continue;
280  debug_stream << printable(*it) << ',';
281  }
282  debug_stream << '\n';
283  }
284  }
285 
287  trace_matrix_type trace_matrix{};
289  typename std::ranges::iterator_t<trace_matrix_type> trace_matrix_iter{};
290 };
291 } // namespace seqan3::detail
::ranges::distance distance
Alias for ranges::distance. Returns the number of hops from first to last.
Definition: iterator:321
Provides seqan3::detail::unbanded_score_trace_dp_matrix.
constexpr sequenced_policy seq
Global execution policy object for sequenced execution policy.
Definition: execution.hpp:54
constexpr auto zip
A range adaptor that transforms a tuple of range into a range of tuples.
Definition: ranges:948
T min(T... args)
debug_stream_type debug_stream
A global instance of seqan3::debug_stream_type.
Definition: debug_stream.hpp:249
::ranges::iterator_t iterator_t
Alias for ranges::iterator_t. Obtains the iterator type of a range.
Definition: ranges:204
::ranges::begin begin
Alias for ranges::begin. Returns an iterator to the beginning of a range.
Definition: ranges:174
::ranges::advance advance
Alias for ranges::advance. Advances the iterator by the given distance.
Definition: iterator:316
T addressof(T... args)
Provides std::span from the C++20 standard library.
constexpr auto iota
Generates a sequence of elements by repeatedly incrementing an initial value.
Definition: ranges:647
Provides seqan3::detail::banded_score_dp_matrix_policy.
Definition: aligned_sequence_concept.hpp:35
Provides seqan3::detail::alignment_coordinate.