SeqAn3  3.2.0-rc.1
The Modern C++ library for sequence analysis.
affine_gap_policy.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 
19 
20 namespace seqan3::detail
21 {
22 
23 // ----------------------------------------------------------------------------
24 // affine_gap_policy
25 // ----------------------------------------------------------------------------
26 
45 template <typename alignment_algorithm_t, typename score_t, typename align_local_t = std::false_type>
46 class affine_gap_policy
47 {
48 private:
50  friend alignment_algorithm_t;
51 
53  using alignment_state_t = alignment_algorithm_state<score_t>;
54 
59  constexpr affine_gap_policy() noexcept = default;
60  constexpr affine_gap_policy(affine_gap_policy const &) noexcept = default;
61  constexpr affine_gap_policy(affine_gap_policy &&) noexcept = default;
62  constexpr affine_gap_policy & operator=(affine_gap_policy const &) noexcept = default;
63  constexpr affine_gap_policy & operator=(affine_gap_policy &&) noexcept = default;
64  ~affine_gap_policy() noexcept = default;
65 
67  template <typename configuration_t>
68  affine_gap_policy(configuration_t const & /*config*/)
69  {}
71 
85  template <typename cell_t>
86  constexpr void
87  compute_cell(cell_t && current_cell, alignment_algorithm_state<score_t> & cache, score_t const score) const noexcept
88  {
89  // score_cell = seqan3::detail::alignment_score_matrix_proxy
90  // trace_cell = seqan3::detail::alignment_trace_matrix_proxy
91  auto & [score_cell, trace_cell] = current_cell;
92  constexpr bool with_trace = !decays_to_ignore_v<std::remove_reference_t<decltype(trace_cell.current)>>;
93  // Precompute the diagonal score.
94  score_t tmp = score_cell.diagonal + score;
95 
96  if constexpr (with_trace)
97  {
98  tmp = (tmp < score_cell.up) ? (trace_cell.current = trace_cell.up, score_cell.up)
99  : (trace_cell.current = trace_directions::diagonal | trace_cell.up, tmp);
100  tmp = (tmp < score_cell.r_left) ? (trace_cell.current = trace_cell.r_left, score_cell.r_left)
101  : (trace_cell.current |= trace_cell.r_left, tmp);
102  }
103  else
104  {
105  tmp = (tmp < score_cell.up) ? score_cell.up : tmp;
106  tmp = (tmp < score_cell.r_left) ? score_cell.r_left : tmp;
107  }
108 
109  if constexpr (align_local_t::value)
110  tmp = (tmp < 0) ? (trace_cell.current = trace_directions::none, 0) : tmp;
111 
112  // Store the current max score.
113  score_cell.current = tmp;
114  // Check if this was the optimum. Possibly a noop.
115  static_cast<alignment_algorithm_t const &>(*this).check_score_of_cell(current_cell, cache);
116 
117  // Prepare horizontal and vertical score for next column.
118  tmp += cache.gap_open_score;
119  score_cell.up += cache.gap_extension_score;
120  score_cell.w_left = score_cell.r_left + cache.gap_extension_score;
121 
122  score_cell.up = (score_cell.up < tmp) ? (trace_cell.up = trace_directions::up_open, tmp)
123  : (trace_cell.up = trace_directions::up, score_cell.up);
124  score_cell.w_left = (score_cell.w_left < tmp) ? (trace_cell.w_left = trace_directions::left_open, tmp)
125  : (trace_cell.w_left = trace_directions::left, score_cell.w_left);
126  }
127 
141  template <typename cell_t>
142  constexpr void compute_first_band_cell(cell_t && current_cell,
143  alignment_algorithm_state<score_t> & cache,
144  score_t const score) const noexcept
145  {
146  // Compute the diagonal score and the compare with the previous horizontal value.
147  // score_cell = seqan3::detail::alignment_score_matrix_proxy
148  // trace_cell = seqan3::detail::alignment_trace_matrix_proxy
149  auto & [score_cell, trace_cell] = current_cell;
150  score_cell.current = score_cell.diagonal + score;
151 
152  score_cell.current = (score_cell.current < score_cell.r_left)
153  ? (trace_cell.current = trace_cell.r_left, score_cell.r_left)
154  : (trace_cell.current = trace_directions::diagonal, score_cell.current);
155 
156  if constexpr (align_local_t::value)
157  {
158  score_cell.current =
159  (score_cell.current < 0) ? (trace_cell.current = trace_directions::none, 0) : score_cell.current;
160  }
161  // Check if this was the optimum. Possibly a noop.
162  static_cast<alignment_algorithm_t const &>(*this).check_score_of_cell(current_cell, cache);
163 
164  // At the top of the band we can not come from up but only diagonal or left, so the next vertical must be a
165  // gap open.
166  score_cell.up = score_cell.current + cache.gap_open_score; // add gap open cost
167  trace_cell.up = trace_directions::up_open;
168  }
169 
182  template <typename alignment_configuration_t>
183  constexpr void initialise_alignment_state(alignment_configuration_t const & config) noexcept
184  {
185  auto scheme =
186  config.get_or(align_cfg::gap_cost_affine{align_cfg::open_score{-10}, align_cfg::extension_score{-1}});
187 
188  alignment_state.gap_extension_score = static_cast<score_t>(scheme.extension_score);
189  alignment_state.gap_open_score = static_cast<score_t>(scheme.extension_score + scheme.open_score);
190  }
191 
192  alignment_state_t alignment_state{};
193 };
194 
195 } // namespace seqan3::detail
Provides seqan3::align_config::gap_cost_affine.
Provides seqan3::detail::alignment_algorithm_state.
Provides the declaration of seqan3::detail::trace_directions.
Provides concepts that do not have equivalents in C++20.