SeqAn3
The Modern C++ library for sequence analysis.
affine_gap_banded_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 <tuple>
16 #include <type_traits>
17 
21 #include <seqan3/std/concepts>
22 
23 namespace seqan3::detail
24 {
25 
26 // ----------------------------------------------------------------------------
27 // affine_gap_banded_policy
28 // ----------------------------------------------------------------------------
29 
36 template <typename derived_type, typename cell_type, typename align_local_t = std::false_type>
37 class affine_gap_banded_policy :
38  public affine_gap_policy<affine_gap_banded_policy<derived_type, cell_type, align_local_t>, cell_type, align_local_t>
39 {
40 private:
41 
43  using base_t = affine_gap_policy<affine_gap_banded_policy<derived_type, cell_type, align_local_t>,
44  cell_type,
45  align_local_t>;
46 
48  friend derived_type;
49 
51  using score_t = typename base_t::score_t;
52 
57  constexpr affine_gap_banded_policy() noexcept = default;
58  constexpr affine_gap_banded_policy(affine_gap_banded_policy const &) noexcept = default;
59  constexpr affine_gap_banded_policy(affine_gap_banded_policy &&) noexcept = default;
60  constexpr affine_gap_banded_policy & operator=(affine_gap_banded_policy const &) noexcept = default;
61  constexpr affine_gap_banded_policy & operator=(affine_gap_banded_policy &&) noexcept = default;
62  ~affine_gap_banded_policy() noexcept = default;
63 
72  template <typename score_cell_type, typename cache_t>
73  constexpr void compute_cell(score_cell_type && current_cell,
74  cache_t & cache,
75  score_t const score) const noexcept
76  {
77  using std::get;
78  // Unpack the cached variables.
79  auto & [score_entry, coordinate, trace_value] = current_cell;
80  auto & [current_entry, next_entry] = score_entry;
81  auto & [main_score, hz_score, hz_trace] = current_entry;
82  auto const & [prev_main_score, prev_hz_score, prev_hz_trace]= next_entry;
83  auto & [prev_cell, gap_open, gap_extend, opt] = cache;
84  auto & [tmp, vt_score, vt_trace] = prev_cell;
85 
86  (void) prev_main_score;
87 
88  //TODO For the local alignment we might be able to use GCC overlow builtin arithmetics, which
89  // allows us to check if overflow/underflow would happen. Not sure, if this helps with the performance though.
90  // (see https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html)
91  main_score += score;
92  // Compute where the max comes from.
93  if constexpr (decays_to_ignore_v<decltype(trace_value)>) // Don't compute a traceback
94  {
95  main_score = (main_score < vt_score) ? vt_score : main_score;
96  main_score = (main_score < prev_hz_score) ? prev_hz_score : main_score;
97  if constexpr (align_local_t::value)
98  main_score = (main_score < 0) ? 0 : main_score;
99  }
100  else // Compute any traceback
101  {
102  main_score = (main_score < vt_score) ? (trace_value = vt_trace, vt_score)
103  : (trace_value = trace_directions::diagonal | vt_trace, main_score);
104  main_score = (main_score < prev_hz_score) ? (trace_value = prev_hz_trace, prev_hz_score)
105  : (trace_value |= prev_hz_trace, main_score);
106  if constexpr (align_local_t::value)
107  main_score = (main_score < 0) ? (trace_value = trace_directions::none, 0) : main_score;
108  }
109 
110  // Check if this was the optimum. Possibly a noop.
111  static_cast<derived_type const &>(*this).check_score(
112  alignment_optimum{main_score, static_cast<alignment_coordinate>(coordinate)}, opt);
113 
114  tmp = main_score + gap_open;
115  vt_score += gap_extend;
116  hz_score = prev_hz_score + gap_extend;
117 
118  if constexpr (decays_to_ignore_v<decltype(trace_value)>)
119  {
120  vt_score = (vt_score < tmp) ? tmp : vt_score;
121  hz_score = (hz_score < tmp) ? tmp : hz_score;
122  }
123  else
124  {
125  vt_score = (vt_score < tmp) ? (vt_trace = trace_directions::up_open, tmp)
126  : (vt_trace = trace_directions::up, vt_score);
127  hz_score = (hz_score < tmp) ? (hz_trace = trace_directions::left_open, tmp)
128  : (hz_trace = trace_directions::left, hz_score);
129  }
130  }
131 
139  template <typename score_cell_type, typename cache_t>
140  constexpr void compute_first_band_cell(score_cell_type && current_cell,
141  cache_t & cache,
142  score_t const score) const noexcept
143  {
144  using std::get;
145  // Unpack the cached variables.
146  auto & [score_entry, coordinate, trace_value] = current_cell;
147  auto & [current_entry, next_entry] = score_entry;
148  auto & main_score = get<0>(current_entry);
149  auto const & prev_hz_score = get<1>(next_entry);
150  auto const & prev_hz_trace = get<2>(next_entry);
151  auto & [tmp, vt_score, vt_trace] = get<0>(cache);
152 
153  (void) tmp;
154 
155  // Compute the diagonal score and the compare with the previous horizontal value.
156  main_score += score;
157 
158  if constexpr (decays_to_ignore_v<decltype(trace_value)>) // Don't compute a traceback
159  {
160  main_score = (main_score < prev_hz_score) ? prev_hz_score : main_score;
161  if constexpr (align_local_t::value)
162  main_score = (main_score < 0) ? 0 : main_score;
163  }
164  else
165  {
166  main_score = (main_score < prev_hz_score) ? (trace_value = prev_hz_trace, prev_hz_score)
167  : (trace_value = trace_directions::diagonal, main_score);
168  if constexpr (align_local_t::value)
169  main_score = (main_score < 0) ? (trace_value = trace_directions::none, 0) : main_score;
170  }
171 
172  // Check if this was the optimum. Possibly a noop.
173  static_cast<derived_type const &>(*this).check_score(
174  alignment_optimum{main_score, static_cast<alignment_coordinate>(coordinate)}, get<3>(cache));
175  // At the top of the band we can not come from up but only diagonal or left, so the next vertical must be a
176  // gap open.
177  vt_score = main_score + get<1>(cache); // add gap open cost
178 
179  // Store vertical open as it is the only way to get into the first cell of the band from a vertical direction.
180  if constexpr (!decays_to_ignore_v<decltype(trace_value)>) // Traceback requested.
181  vt_trace = trace_directions::up_open;
182 
183  // Don't compute the horizontal value, since it will never be used.
184  }
185 
186  using base_t::make_cache;
187 };
188 
189 } // namespace seqan3::detail
The Concepts library.
Provides seqan3::detail::affine_gap_policy.
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 various type traits on generic types.
Provides various transformation traits used by the range module.