SeqAn3  3.2.0-rc.1
The Modern C++ library for sequence analysis.
edit_distance_trace_matrix_full.hpp
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------------------------------
2 // Copyright (c) 2006-2022, Knut Reinert & Freie Universität Berlin
3 // Copyright (c) 2016-2022, 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 <bitset>
16 
20 
21 namespace seqan3::detail
22 {
23 
31 template <typename word_t, bool is_semi_global, bool use_max_errors>
32 class edit_distance_trace_matrix_full
33 {
34 public:
36  template <std::ranges::viewable_range database_t,
37  std::ranges::viewable_range query_t,
38  typename align_config_t,
39  typename edit_traits>
40  friend class edit_distance_unbanded;
41 
45  edit_distance_trace_matrix_full() = default;
46  edit_distance_trace_matrix_full(edit_distance_trace_matrix_full const &) = default;
47  edit_distance_trace_matrix_full(edit_distance_trace_matrix_full &&) = default;
48  edit_distance_trace_matrix_full & operator=(edit_distance_trace_matrix_full const &) = default;
49  edit_distance_trace_matrix_full & operator=(edit_distance_trace_matrix_full &&) = default;
50  ~edit_distance_trace_matrix_full() = default;
51 
52 protected:
54  template <typename derived_t, typename edit_traits>
55  friend class edit_distance_unbanded_trace_matrix_policy;
56 
60  edit_distance_trace_matrix_full(size_t const rows_size) : rows_size{rows_size}, columns{}
61  {}
63 
64 private:
65  struct trace_path_iterator;
66 
67 public:
69  using word_type = word_t;
70 
72  static constexpr auto word_size = bits_of<word_type>;
73 
75  using value_type = detail::trace_directions;
76 
78  using reference = value_type;
79 
81  using size_type = size_t;
82 
91  void reserve(size_t const new_capacity)
92  {
93  columns.reserve(new_capacity);
94  }
95 
96 public:
98  reference at(matrix_coordinate const & coordinate) const noexcept
99  {
100  size_t row = coordinate.row;
101  size_t col = coordinate.col;
102 
103  assert(row < rows());
104  assert(col < cols());
105 
106  column_type const & column = columns[col];
107 
108  if constexpr (use_max_errors)
109  if (!(row < column.max_rows))
111 
112  if (row == 0u)
113  {
114  if constexpr (is_semi_global)
116 
117  if (col == 0u)
119 
120  return detail::trace_directions::left;
121  }
122 
123  size_t const idx = (row - 1u) / word_size;
124  size_t const offset = (row - 1u) % word_size;
125 
126  bool const left = std::bitset<word_size>(column.left[idx])[offset];
127  bool const diagonal = std::bitset<word_size>(column.diagonal[idx])[offset];
128  bool const up = std::bitset<word_size>(column.up[idx])[offset];
129 
130  auto const dir = (left ? detail::trace_directions::left : detail::trace_directions::none)
131  | (diagonal ? detail::trace_directions::diagonal : detail::trace_directions::none)
132  | (up ? detail::trace_directions::up : detail::trace_directions::none);
133 
134  return dir;
135  }
136 
138  size_t rows() const noexcept
139  {
140  return rows_size;
141  }
142 
144  size_t cols() const noexcept
145  {
146  return columns.size();
147  }
148 
155  auto trace_path(matrix_coordinate const & trace_begin) const
156  {
157  if (trace_begin.row >= rows() || trace_begin.col >= cols())
158  throw std::invalid_argument{"The given coordinate exceeds the matrix in vertical or horizontal direction."};
159 
160  using path_t = std::ranges::subrange<trace_path_iterator, std::default_sentinel_t>;
161  return path_t{trace_path_iterator{this, trace_begin}, std::default_sentinel};
162  }
163 
164 protected:
166  struct max_errors_state
167  {
170  size_t max_rows{};
171  };
172 
174  struct trace_matrix_state
175  {
179  std::vector<word_type> diagonal{};
182  };
183 
185  struct column_type : enable_state_t<true, trace_matrix_state>, enable_state_t<use_max_errors, max_errors_state>
186  {};
187 
193  void add_column(std::vector<word_type> left, std::vector<word_type> diagonal, std::vector<word_type> up)
194  requires (!use_max_errors)
195  {
196  column_type column{};
197  column.left = std::move(left);
198  column.diagonal = std::move(diagonal);
199  column.up = std::move(up);
200 
201  columns.push_back(std::move(column));
202  }
203 
210  void add_column(std::vector<word_type> left,
211  std::vector<word_type> diagonal,
213  size_t const max_rows)
214  requires use_max_errors
215  {
216  column_type column{};
217  column.left = std::move(left);
218  column.diagonal = std::move(diagonal);
219  column.up = std::move(up);
220  column.max_rows = max_rows;
221 
222  columns.push_back(std::move(column));
223  }
224 
225 private:
227  size_t rows_size{};
229  std::vector<column_type> columns{};
230 };
231 
246 template <typename word_t, bool is_semi_global, bool use_max_errors>
247 struct edit_distance_trace_matrix_full<word_t, is_semi_global, use_max_errors>::trace_path_iterator
248 {
253  using iterator_category = std::input_iterator_tag;
255  using value_type = detail::trace_directions;
257  using difference_type = std::ptrdiff_t;
259 
261  static constexpr value_type D = value_type::diagonal;
263  static constexpr value_type L = value_type::left;
265  static constexpr value_type U = value_type::up;
267  static constexpr value_type N = value_type::none;
268 
273  constexpr value_type operator*() const
274  {
275  value_type dir = parent->at(coordinate());
276 
277  if (dir == N)
278  return N;
279 
280  if ((dir & L) == L)
281  return L;
282  else if ((dir & U) == U)
283  return U;
284  else
285  return D;
286  }
287 
289  [[nodiscard]] constexpr matrix_coordinate const & coordinate() const
290  {
291  return coordinate_;
292  }
294 
299  constexpr trace_path_iterator & operator++()
300  {
301  value_type dir = *(*this);
302 
303  if ((dir & L) == L)
304  {
305  coordinate_.col = std::max<size_t>(coordinate_.col, 1) - 1;
306  }
307  else if ((dir & U) == U)
308  {
309  coordinate_.row = std::max<size_t>(coordinate_.row, 1) - 1;
310  }
311  else if ((dir & D) == D)
312  {
313  coordinate_.row = std::max<size_t>(coordinate_.row, 1) - 1;
314  coordinate_.col = std::max<size_t>(coordinate_.col, 1) - 1;
315  }
316 
317  // seqan3::trace_direction::none in an inner cell of the trace matrix is not possible in
318  // non-local alignments, e.g. global, semi-global, max-error, banded. On the other hand,
319  // this can happen in a cell of the first row or first colomn.
320  assert(dir != N || coordinate_.row == 0 || coordinate_.col == 0);
321 
322  return *this;
323  }
324 
326  constexpr void operator++(int)
327  {
328  ++(*this);
329  }
331 
336  friend bool operator==(trace_path_iterator const & it, std::default_sentinel_t)
337  {
338  return *it == value_type::none;
339  }
340 
342  friend bool operator==(std::default_sentinel_t, trace_path_iterator const & it)
343  {
344  return it == std::default_sentinel;
345  }
346 
348  friend bool operator!=(trace_path_iterator const & it, std::default_sentinel_t)
349  {
350  return !(it == std::default_sentinel);
351  }
352 
354  friend bool operator!=(std::default_sentinel_t, trace_path_iterator const & it)
355  {
356  return it != std::default_sentinel;
357  }
359 
361  edit_distance_trace_matrix_full const * parent{nullptr};
363  matrix_coordinate coordinate_{};
364 };
365 
366 } // namespace seqan3::detail
Provides utility functions for bit twiddling.
Forwards for seqan3::edit_distance_unbanded related types.
requires requires
The rank_type of the semi-alphabet; defined as the return type of seqan3::to_rank....
Definition: alphabet/concept.hpp:164
@ none
No flag is set.
Definition: debug_stream_type.hpp:32
@ offset
Sequence (seqan3::field::seq) relative start position (0-based), unsigned value.
T left(T... args)
matrix_index< size_t > matrix_coordinate
A coordinate type to access an element inside of a two-dimensional matrix.
Definition: matrix_coordinate.hpp:178
T operator!=(T... args)
T size(T... args)
Provides the declaration of seqan3::detail::trace_directions.