SeqAn3
The Modern C++ library for sequence analysis.
unbanded_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/reverse.hpp>
19 #include <range/v3/view/zip.hpp>
20 
25 #include <seqan3/std/span>
26 
27 namespace seqan3::detail
28 {
29 
33 struct gap_segment
34 {
36  size_t position;
38  size_t size;
39 };
40 
52 template <typename derived_t, typename score_allocator_t, typename trace_allocator_t>
53 class unbanded_score_trace_dp_matrix_policy :
54  public unbanded_score_dp_matrix_policy<unbanded_score_trace_dp_matrix_policy<derived_t,
55  score_allocator_t,
56  trace_allocator_t>,
57  score_allocator_t>
58 {
59 private:
60 
62  using base_t = unbanded_score_dp_matrix_policy<unbanded_score_trace_dp_matrix_policy<derived_t,
63  score_allocator_t,
64  trace_allocator_t>,
65  score_allocator_t>;
66 
68  friend derived_t;
70  using base_t::dimension_first_range;
72  using base_t::dimension_second_range;
74  using base_t::score_matrix;
76  using base_t::current_column_index;
77 
81  using cell_type = typename score_allocator_t::value_type;
85  using trace_type = typename trace_allocator_t::value_type;
87  using score_matrix_type = std::vector<cell_type, score_allocator_t>;
89  using trace_matrix_type = std::vector<trace_type, trace_allocator_t>;
91 
95  constexpr unbanded_score_trace_dp_matrix_policy() = default;
96  constexpr unbanded_score_trace_dp_matrix_policy(unbanded_score_trace_dp_matrix_policy const &) = default;
98  constexpr unbanded_score_trace_dp_matrix_policy(unbanded_score_trace_dp_matrix_policy &&) = default;
99  constexpr unbanded_score_trace_dp_matrix_policy & operator=(unbanded_score_trace_dp_matrix_policy const &) = default;
102  constexpr unbanded_score_trace_dp_matrix_policy & operator=(unbanded_score_trace_dp_matrix_policy &&) = default;
103  ~unbanded_score_trace_dp_matrix_policy() = default;
104 
112  template <typename first_range_t, typename second_range_t>
113  constexpr void allocate_matrix(first_range_t & first_range, second_range_t & second_range)
114  {
115  base_t::allocate_matrix(first_range, second_range);
116 
117  // We use the full matrix to store the trace direction.
118  trace_matrix.resize(dimension_first_range * dimension_second_range);
119  trace_matrix_iter = trace_matrix.begin();
120  }
121 
123  constexpr auto current_column()
124  {
125  advanceable_alignment_coordinate<advanceable_alignment_coordinate_state::row>
126  col_begin{column_index_type{current_column_index}, row_index_type{0u}};
127  advanceable_alignment_coordinate<advanceable_alignment_coordinate_state::row>
128  col_end{column_index_type{current_column_index}, row_index_type{dimension_second_range}};
129 
130  return std::view::zip(std::span{score_matrix},
131  std::view::iota(col_begin, col_end),
132  std::span{std::addressof(*trace_matrix_iter), dimension_second_range});
133  }
134 
136  constexpr void go_next_column() noexcept
137  {
138  base_t::go_next_column();
139  trace_matrix_iter += dimension_second_range;
140  }
141 
148  constexpr auto parse_traceback(alignment_coordinate const & back_coordinate)
149  {
150  // Store the trace segments.
151  std::deque<gap_segment> first_segments{};
152  std::deque<gap_segment> second_segments{};
153 
154  // Put the iterator to the position where the traceback starts.
155  auto direction_iter = std::ranges::begin(trace_matrix);
156  std::ranges::advance(direction_iter,
157  back_coordinate.first * dimension_second_range + back_coordinate.second);
158 
159  // Parse the trace until interrupt.
160  while (*direction_iter != trace_directions::none)
161  {
162  // parse until end of diagonal run
163  while (static_cast<bool>(*direction_iter & trace_directions::diagonal))
164  {
165  std::ranges::advance(direction_iter, -dimension_second_range - 1);
166  }
167 
168  // parse vertical gap -> record gap in first_segments (will be translated into gap of first sequence)
169  if (static_cast<bool>(*direction_iter & trace_directions::up) ||
170  static_cast<bool>(*direction_iter & trace_directions::up_open))
171  {
172  // Get the current column index (note the column based layout)
173  size_t pos = std::ranges::distance(std::ranges::begin(trace_matrix), direction_iter) /
174  dimension_second_range;
175  gap_segment gap{pos, 0u};
176 
177  // Follow gap until open signal is detected.
178  while (!static_cast<bool>(*direction_iter & trace_directions::up_open))
179  {
180  --direction_iter;
181  ++gap.size;
182  }
183  // explicitly follow opening gap
184  --direction_iter;
185  ++gap.size;
186  // record the gap
187  first_segments.push_front(std::move(gap));
188  continue;
189  }
190  // parse horizontal gap -> record gap in second_segments (will be translated into gap of second sequence)
191  if (static_cast<bool>(*direction_iter & trace_directions::left) ||
192  static_cast<bool>(*direction_iter & trace_directions::left_open))
193  {
194  // Get the current row index (note the column based layout)
195  size_t pos = std::ranges::distance(std::ranges::begin(trace_matrix), direction_iter) %
196  dimension_second_range;
197  gap_segment gap{pos, 0u};
198 
199  // Follow gap until open signal is detected.
200  while (!static_cast<bool>(*direction_iter & trace_directions::left_open))
201  {
202  std::ranges::advance(direction_iter, -dimension_second_range);
203  ++gap.size;
204  }
205  // explicitly follow opening gap
206  std::ranges::advance(direction_iter, -dimension_second_range);
207  ++gap.size;
208  second_segments.push_front(std::move(gap));
209  }
210  }
211 
212  // Get front coordinate.
213  auto c = column_index_type{std::ranges::distance(std::ranges::begin(trace_matrix), direction_iter) /
214  dimension_second_range};
215  auto r = row_index_type{std::ranges::distance(std::ranges::begin(trace_matrix), direction_iter) %
216  dimension_second_range};
217 
218  return std::tuple{alignment_coordinate{column_index_type{std::move(c)}, row_index_type{std::move(r)}},
219  first_segments,
220  second_segments};
221  }
222 
224  void print_trace_matrix()
225  {
226  auto printable = [](trace_directions dir)
227  {
228  std::string seq{};
229  if (dir == trace_directions::none)
230  seq.append("0");
231  if ((dir & trace_directions::diagonal) == trace_directions::diagonal)
232  seq.append("\\");
233  if ((dir & trace_directions::up) == trace_directions::up)
234  seq.append("|");
235  if ((dir & trace_directions::up_open) == trace_directions::up_open)
236  seq.append("^");
237  if ((dir & trace_directions::left) == trace_directions::left)
238  seq.append("-");
239  if ((dir & trace_directions::left_open) == trace_directions::left_open)
240  seq.append("<");
241  return seq;
242  };
243 
244  for (size_t row = 0; row < dimension_second_range; ++row)
245  {
246  for (size_t col = 0; col < dimension_first_range; ++col)
247  {
248  debug_stream << printable(trace_matrix[col * dimension_second_range + row]) << " ";
249  }
250  debug_stream << "\n";
251  }
252  }
253 
255  trace_matrix_type trace_matrix{};
257  std::ranges::iterator_t<trace_matrix_type> trace_matrix_iter{};
258 };
259 } // namespace seqan3::detail
::ranges::distance distance
Alias for ranges::distance. Returns the number of hops from first to last.
Definition: iterator:321
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
::ranges::size size
Alias for ranges::size. Obtains the size of a range whose size can be calculated in constant time...
Definition: ranges:189
Provides seqan3::detail::unbanded_score_dp_matrix_policy.
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 the declaration of seqan3::detail::trace_directions.
Definition: aligned_sequence_concept.hpp:35
Provides seqan3::view::persist.
Provides seqan3::detail::alignment_coordinate.