SeqAn3
The Modern C++ library for sequence analysis.
banded_score_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 <limits>
16 #include <memory>
17 
18 #include <range/v3/view/repeat_n.hpp>
19 
25 #include <seqan3/std/ranges>
26 #include <seqan3/std/span>
27 
28 namespace seqan3::detail
29 {
30 
36 template <typename derived_t, typename allocator_type>
37 class banded_score_dp_matrix_policy :
38  public unbanded_score_dp_matrix_policy<banded_score_dp_matrix_policy<derived_t, allocator_type>, allocator_type>
39 {
40 private:
41 
43  using base_t = unbanded_score_dp_matrix_policy<banded_score_dp_matrix_policy<derived_t, allocator_type>, allocator_type>;
44 
46  friend derived_t;
47 
51  using cell_type = typename base_t::cell_type;
54  using score_matrix_type = typename base_t::score_matrix_type;
56 
58  static constexpr std::tuple_element_t<0, cell_type> INF =
60 
64  constexpr banded_score_dp_matrix_policy() = default;
65  constexpr banded_score_dp_matrix_policy(banded_score_dp_matrix_policy const &) = default;
66  constexpr banded_score_dp_matrix_policy(banded_score_dp_matrix_policy &&) = default;
67  constexpr banded_score_dp_matrix_policy & operator=(banded_score_dp_matrix_policy const &) = default;
68  constexpr banded_score_dp_matrix_policy & operator=(banded_score_dp_matrix_policy &&) = default;
69  ~banded_score_dp_matrix_policy() = default;
70 
72 public:
73 
82  template <typename first_range_t, typename second_range_t, typename band_t>
83  constexpr void allocate_matrix(first_range_t & first_range, second_range_t & second_range, band_t const & band)
84  {
85  dimension_first_range = std::ranges::distance(std::ranges::begin(first_range), std::ranges::end(first_range)) + 1;
86  dimension_second_range = std::ranges::distance(std::ranges::begin(second_range), std::ranges::end(second_range)) + 1;
87 
88  // If upper_bound is negative, set it to 0 and trim the second sequences accordingly.
89  band_column_index = std::max(static_cast<uint_fast32_t>(band.upper_bound), static_cast<uint_fast32_t>(0));
90  // If lower_bound is positive, set it to 0 and trim the first sequences accordingly.
91  band_row_index = std::abs(std::min(static_cast<int_fast32_t>(band.lower_bound),
92  static_cast<int_fast32_t>(0)));
93 
94  // If the band is wider than the sequence length, limit the band width.
95  band_column_index = std::min(band_column_index, static_cast<uint_fast32_t>(dimension_second_range - 1));
96  band_row_index = std::min(band_row_index, static_cast<uint_fast32_t>(dimension_first_range - 1));
97 
98  band_size = band_column_index + band_row_index + 1;
99 
100  // Reserve one more cell to deal with last cell in the banded column which needs only the diagonal and up cell.
101  // TODO: introduce specific named cell types with initialisation values.
102  score_matrix.resize(band_size + 1);
103 
104  using std::get;
105  get<0>(score_matrix.back()) = INF;
106  get<1>(score_matrix.back()) = INF;
107  get<2>(score_matrix.back()) = trace_directions::none;
108  current_column_index = 0;
109  // Position the iterator to the right offset within the band.
110  current_matrix_iter = std::ranges::begin(score_matrix) + band_column_index;
111  }
112 
114  constexpr auto current_column() noexcept
115  {
116  auto span = current_band_size();
117 
118  assert(span > 0u); // The span must always be greater than 0.
119 
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  ranges::view::repeat_n(std::ignore, span) | std::view::common);
133  }
134 
136  constexpr void go_next_column() noexcept
137  {
138  // Update the current_column_index.
139  base_t::go_next_column();
140  // Still in the initialisation phase and need to update the current matrix iterator until begin is reached.
141  if (current_matrix_iter != std::ranges::begin(score_matrix))
142  --current_matrix_iter;
143  }
144 
146  constexpr uint_fast32_t current_band_size() const noexcept
147  {
148  // Distance from begin of band until end of the entire column (not end of the band).
149  int_fast32_t remaining_column_size =
150  static_cast<int_fast32_t>(dimension_second_range) -
151  std::max(static_cast<int_fast32_t>(0), static_cast<int_fast32_t>(current_column_index - band_column_index));
152 
153  // The matrix was trimmed to fit the band exactly, thus the second term cannot be greater or equal
154  // than the first term in the equation above.
155  assert(remaining_column_size > 0);
156 
157  // using const_iter = typename score_matrix_type::const_iterator;
158  // The current band size is the min of the remaining column size and the size of the current span of the band.
159  return std::min(static_cast<uint_fast32_t>(remaining_column_size),
160  static_cast<uint_fast32_t>(
161  std::ranges::distance(current_matrix_iter, std::ranges::end(score_matrix)) - 1));
162  }
163 
172  constexpr uint_fast32_t second_range_begin_offset() const noexcept
173  {
174  assert(current_column_index > band_column_index);
175 
176  return current_column_index - band_column_index - 1;
177  }
178 
180  constexpr bool band_touches_last_row() const noexcept
181  {
182  if (current_column_index > band_column_index) // TODO [[likely]]
183  return (second_range_begin_offset() + current_band_size() + 1) == dimension_second_range;
184  else
185  return current_band_size() >= dimension_second_range;
186  }
187 
201  template <typename first_range_t, typename second_range_t, typename band_t>
202  constexpr auto trim_sequences(first_range_t & first_range,
203  second_range_t & second_range,
204  band_t const & band) const noexcept
205  {
206  using band_type = decltype(band.lower_bound);
207 
208  band_type dimension_first = std::ranges::distance(std::ranges::begin(first_range), std::ranges::end(first_range));
209  band_type dimension_second = std::ranges::distance(std::ranges::begin(second_range), std::ranges::end(second_range));
210 
211  auto trim_first_range = [&]() constexpr
212  {
213  size_t begin_pos = std::max(band.lower_bound - 1, static_cast<band_type>(0));
214  size_t end_pos = std::min(band.upper_bound + dimension_second, dimension_first);
215  return first_range | view::slice(begin_pos, end_pos);
216  };
217 
218  auto trim_second_range = [&]() constexpr
219  {
220  size_t begin_pos = std::abs(std::min(band.upper_bound + 1, static_cast<band_type>(0)));
221  size_t end_pos = std::min(dimension_first - band.lower_bound, dimension_second);
222  return second_range | view::slice(begin_pos, end_pos);
223  };
224 
225  return std::tuple{trim_first_range(), trim_second_range()};
226  }
227 
231  constexpr auto map_banded_coordinate_to_range_position(alignment_coordinate coordinate) const noexcept
232  {
233  using as_int_t = std::make_signed_t<decltype(coordinate.first)>;
234  // Refine the row coordinate to match the original sequence coordinates since the first position of the
235  // trace matrix is shifted by the value of the band_column_index, i.e. the upper bound of the band.
236  //
237  // case 1: ends in column before the band_column_index: subtract the offset from the actual row coordinate.
238  // case 2: ends in column after the band_column_index: add the offset to the actual row coordinate.
239  coordinate.second += static_cast<as_int_t>(coordinate.first - band_column_index);
240  return coordinate;
241  }
242 
243 private:
244 
245  using base_t::score_matrix;
246  using base_t::dimension_first_range;
247  using base_t::dimension_second_range;
248  using base_t::current_column_index;
249 
251  typename score_matrix_type::iterator current_matrix_iter;
253  uint_fast32_t band_column_index{};
255  uint_fast32_t band_row_index{};
257  uint_fast32_t band_size{};
258 };
259 
260 } // namespace seqan3::detail
::ranges::distance distance
Alias for ranges::distance. Returns the number of hops from first to last.
Definition: iterator:321
constexpr auto zip
A range adaptor that transforms a tuple of range into a range of tuples.
Definition: ranges:948
Provides various shortcuts for common std::ranges functions.
constexpr auto slice
A view adaptor that returns a half-open interval on the underlying range.
Definition: slice.hpp:144
T min(T... args)
constexpr auto common
A range adaptor that makes any range model std::ranges::CommonRange (at the expense of some performan...
Definition: ranges:447
Provides seqan3::detail::unbanded_score_dp_matrix_policy.
Adaptations of concepts from the Ranges TS.
::ranges::begin begin
Alias for ranges::begin. Returns an iterator to the beginning of a range.
Definition: ranges:174
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.
T max(T... args)
Definition: aligned_sequence_concept.hpp:35
auto const get
A view calling std::get on each element in a range.
Definition: get.hpp:66
Provides seqan3::view::slice.
Provides seqan3::detail::alignment_coordinate.
constexpr auto repeat_n
A view factory that repeats a given value n times.
Definition: repeat_n.hpp:97
::ranges::end end
Alias for ranges::end. Returns an iterator to the end of a range.
Definition: ranges:179