Loading...
Searching...
No Matches
edit_distance_unbanded.hpp
Go to the documentation of this file.
// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: BSD-3-Clause
#pragma once
#include <algorithm>
#include <bitset>
#include <ranges>
#include <utility>
#include <seqan3/alignment/matrix/detail/advanceable_alignment_coordinate.hpp>
#include <seqan3/alignment/matrix/detail/edit_distance_score_matrix_full.hpp>
#include <seqan3/alignment/matrix/detail/edit_distance_trace_matrix_full.hpp>
#include <seqan3/alignment/matrix/detail/matrix_concept.hpp>
#include <seqan3/alignment/pairwise/alignment_result.hpp>
#include <seqan3/alignment/pairwise/edit_distance_fwd.hpp>
#include <seqan3/core/configuration/configuration.hpp>
namespace seqan3::detail
{
template <typename derived_t, typename edit_traits>
class edit_distance_unbanded_max_errors_policy :
edit_traits
{
protected:
static_assert(edit_traits::use_max_errors, "This policy assumes that edit_traits::use_max_errors is true.");
friend derived_t;
edit_distance_unbanded_max_errors_policy() noexcept = default;
edit_distance_unbanded_max_errors_policy(edit_distance_unbanded_max_errors_policy const &) noexcept =
default;
edit_distance_unbanded_max_errors_policy(edit_distance_unbanded_max_errors_policy &&) noexcept =
default;
edit_distance_unbanded_max_errors_policy &
operator=(edit_distance_unbanded_max_errors_policy const &) noexcept = default;
edit_distance_unbanded_max_errors_policy &
operator=(edit_distance_unbanded_max_errors_policy &&) noexcept = default;
~edit_distance_unbanded_max_errors_policy() noexcept = default;
using typename edit_traits::score_type;
using typename edit_traits::word_type;
score_type max_errors{255};
size_t last_block{0u};
word_type last_score_mask{};
void max_errors_init(size_t block_count) noexcept
{
derived_t * self = static_cast<derived_t *>(this);
max_errors = -get<align_cfg::min_score>(self->config).score;
assert(max_errors >= score_type{0});
if (std::ranges::empty(self->query)) // [[unlikely]]
{
last_block = 0u;
self->score_mask = 0u;
last_score_mask = self->score_mask;
return;
}
last_block = block_count - 1u;
last_score_mask = self->score_mask;
// local_max_errors either stores the maximal number of _score (me.max_errors) or the needle size minus one. It
// is used for the mask computation and setting the initial score (the minus one is there because of the Ukkonen
// trick).
size_t const local_max_errors = std::min<size_t>(max_errors, std::ranges::size(self->query) - 1u);
self->score_mask = word_type{1u} << (local_max_errors % self->word_size);
last_block = std::min(local_max_errors / self->word_size, last_block);
self->_score = local_max_errors + 1u;
}
bool is_last_active_cell_within_last_row() const noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
return (self->score_mask == this->last_score_mask) && (this->last_block == self->vp.size() - 1u);
}
bool prev_last_active_cell() noexcept
{
derived_t * self = static_cast<derived_t *>(this);
self->score_mask >>= 1u;
if (self->score_mask != 0u)
return true;
if constexpr (edit_traits::is_global)
{
if (last_block == 0u) // [[unlikely]]
return false;
}
last_block--;
self->score_mask = word_type{1u} << (edit_traits::word_size - 1u);
return true;
}
void next_last_active_cell() noexcept
{
derived_t * self = static_cast<derived_t *>(this);
self->score_mask <<= 1u;
if (self->score_mask)
return;
self->score_mask = 1u;
last_block++;
}
bool update_last_active_cell() noexcept
{
derived_t * self = static_cast<derived_t *>(this);
// update the last active cell
while (!(self->_score <= max_errors))
{
self->advance_score(self->vn[last_block], self->vp[last_block], self->score_mask);
if (!prev_last_active_cell())
{
// prev_last_active_cell = false can only happen for global alignments
assert(edit_traits::is_global);
// we abort here if we don't need to compute a matrix, because the continued
// computation can't produce an alignment.
return !edit_traits::compute_matrix;
}
}
if (is_last_active_cell_within_last_row())
{
assert(self->_score <= max_errors);
if constexpr (edit_traits::is_semi_global)
self->update_best_score();
return self->on_hit();
}
else
{
next_last_active_cell();
self->advance_score(self->vp[last_block], self->vn[last_block], self->score_mask);
}
return false;
}
static size_t max_rows(word_type const score_mask,
unsigned const last_block,
score_type const score,
score_type const max_errors) noexcept
{
using score_matrix_type = typename edit_traits::score_matrix_type;
return score_matrix_type::max_rows(score_mask, last_block, score, max_errors);
}
};
template <typename derived_t, typename edit_traits>
class edit_distance_unbanded_global_policy :
edit_traits
{
protected:
static_assert(edit_traits::is_global || edit_traits::is_semi_global,
"This policy assumes that edit_traits::is_global or edit_traits::is_semi_global is true.");
friend derived_t;
edit_distance_unbanded_global_policy() noexcept = default;
edit_distance_unbanded_global_policy(edit_distance_unbanded_global_policy const &) noexcept =
default;
edit_distance_unbanded_global_policy(edit_distance_unbanded_global_policy &&) noexcept = default;
edit_distance_unbanded_global_policy &
operator=(edit_distance_unbanded_global_policy const &) noexcept = default;
edit_distance_unbanded_global_policy &
operator=(edit_distance_unbanded_global_policy &&) noexcept = default;
~edit_distance_unbanded_global_policy() noexcept = default;
using typename edit_traits::score_type;
score_type _best_score{};
void score_init() noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
_best_score = self->_score;
}
bool is_valid() const noexcept
{
[[maybe_unused]] derived_t const * self = static_cast<derived_t const *>(this);
// This condition uses the observation that after each computation of a column, _score has either the initial
// value of the first row (i.e. the entire column consist of INF's), has the value _score = max_errors + 1
// (there exists a cell within the column that has value <= max_errors, but is not on the last row) or _score <=
// max_errors (the score of the last active cell is <= max_errors)
if constexpr (edit_traits::use_max_errors)
return _best_score <= self->max_errors;
// When not using max_errors there is always a valid alignment, because the last row will always be updated and
// with it the score.
return true;
}
seqan3::detail::advanceable_alignment_coordinate<> invalid_coordinate() const noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
return {column_index_type{std::ranges::size(self->database)}, row_index_type{std::ranges::size(self->query)}};
}
void update_best_score() noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
_best_score = self->_score;
}
size_t end_positions_first() const noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
return std::ranges::size(self->database);
}
public:
std::optional<score_type> score() const noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
static_assert(edit_traits::compute_score,
"score() can only be computed if you specify the result type within "
"your alignment config.");
if (!self->is_valid())
return std::nullopt;
return -_best_score;
}
seqan3::detail::advanceable_alignment_coordinate<> end_positions() const noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
static_assert(edit_traits::compute_end_positions,
"end_positions() can only be computed if you specify the "
"result type within your alignment config.");
if (!self->is_valid())
return self->invalid_coordinate();
column_index_type const first{self->end_positions_first()};
row_index_type const second{std::ranges::size(self->query)};
return {first, second};
}
};
template <typename derived_t, typename edit_traits>
class edit_distance_unbanded_semi_global_policy : public edit_distance_unbanded_global_policy<derived_t, edit_traits>
{
protected:
static_assert(edit_traits::is_semi_global, "This policy assumes that edit_traits::is_semi_global is true.");
friend derived_t;
edit_distance_unbanded_semi_global_policy() noexcept = default;
edit_distance_unbanded_semi_global_policy(edit_distance_unbanded_semi_global_policy const &) noexcept =
default;
edit_distance_unbanded_semi_global_policy(edit_distance_unbanded_semi_global_policy &&) noexcept =
default;
edit_distance_unbanded_semi_global_policy &
operator=(edit_distance_unbanded_semi_global_policy const &) noexcept = default;
edit_distance_unbanded_semi_global_policy &
operator=(edit_distance_unbanded_semi_global_policy &&) noexcept = default;
~edit_distance_unbanded_semi_global_policy() noexcept = default;
using base_t = edit_distance_unbanded_global_policy<derived_t, edit_traits>;
using database_iterator = typename edit_traits::database_iterator;
using base_t::_best_score;
database_iterator _best_score_col{};
void score_init() noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
base_t::score_init();
_best_score_col = self->database_it_end;
}
void update_best_score() noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
// we have to make sure that update_best_score is only called after a score update within the last row.
if constexpr (edit_traits::use_max_errors)
{
assert(std::ranges::empty(self->query) || self->is_last_active_cell_within_last_row());
}
_best_score_col = (self->_score <= _best_score) ? self->database_it : _best_score_col;
_best_score = (self->_score <= _best_score) ? self->_score : _best_score;
}
size_t end_positions_first() const noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
// offset == 0u is a special case if database sequence is empty, because in this case the best column is zero.
size_t offset = std::ranges::empty(self->database) ? 0u : 1u;
return std::ranges::distance(std::ranges::begin(self->database), _best_score_col) + offset;
}
};
template <typename derived_t, typename edit_traits>
class edit_distance_unbanded_score_matrix_policy :
edit_traits
{
protected:
static_assert(edit_traits::compute_score_matrix,
"This policy assumes that edit_traits::compute_score_matrix is true.");
friend derived_t;
edit_distance_unbanded_score_matrix_policy() noexcept = default;
edit_distance_unbanded_score_matrix_policy(edit_distance_unbanded_score_matrix_policy const &) noexcept =
default;
edit_distance_unbanded_score_matrix_policy(edit_distance_unbanded_score_matrix_policy &&) noexcept =
default;
edit_distance_unbanded_score_matrix_policy &
operator=(edit_distance_unbanded_score_matrix_policy const &) noexcept = default;
edit_distance_unbanded_score_matrix_policy &
operator=(edit_distance_unbanded_score_matrix_policy &&) noexcept = default;
~edit_distance_unbanded_score_matrix_policy() noexcept = default;
using typename edit_traits::score_matrix_type;
score_matrix_type _score_matrix{};
void score_matrix_init()
{
derived_t const * self = static_cast<derived_t const *>(this);
_score_matrix = score_matrix_type{std::ranges::size(self->query) + 1u};
_score_matrix.reserve(std::ranges::size(self->database) + 1u);
}
public:
score_matrix_type const & score_matrix() const noexcept
{
static_assert(edit_traits::compute_score_matrix,
"score_matrix() can only be computed if you specify the "
"result type within your alignment config.");
return _score_matrix;
}
};
template <typename derived_t, typename edit_traits>
class edit_distance_unbanded_trace_matrix_policy :
edit_traits
{
protected:
static_assert(edit_traits::compute_trace_matrix,
"This policy assumes that edit_traits::compute_trace_matrix is true.");
friend derived_t;
edit_distance_unbanded_trace_matrix_policy() noexcept = default;
edit_distance_unbanded_trace_matrix_policy(edit_distance_unbanded_trace_matrix_policy const &) noexcept =
default;
edit_distance_unbanded_trace_matrix_policy(edit_distance_unbanded_trace_matrix_policy &&) noexcept =
default;
edit_distance_unbanded_trace_matrix_policy &
operator=(edit_distance_unbanded_trace_matrix_policy const &) noexcept = default;
edit_distance_unbanded_trace_matrix_policy &
operator=(edit_distance_unbanded_trace_matrix_policy &&) noexcept = default;
~edit_distance_unbanded_trace_matrix_policy() noexcept = default;
using typename edit_traits::alignment_result_type;
using typename edit_traits::trace_matrix_type;
using typename edit_traits::word_type;
std::vector<word_type> hp{};
std::vector<word_type> db{};
trace_matrix_type _trace_matrix{};
void trace_matrix_init(size_t block_count)
{
derived_t const * self = static_cast<derived_t const *>(this);
_trace_matrix = trace_matrix_type{std::ranges::size(self->query) + 1u};
_trace_matrix.reserve(std::ranges::size(self->database) + 1u);
hp.resize(block_count, 0u);
db.resize(block_count, 0u);
}
public:
trace_matrix_type const & trace_matrix() const noexcept
{
static_assert(edit_traits::compute_trace_matrix,
"trace_matrix() can only be computed if you specify the "
"result type within your alignment config.");
return _trace_matrix;
}
seqan3::detail::advanceable_alignment_coordinate<> begin_positions() const noexcept
{
derived_t const * self = static_cast<derived_t const *>(this);
static_assert(edit_traits::compute_begin_positions,
"begin_positions() can only be computed if you specify the "
"result type within your alignment config.");
if (!self->is_valid())
return self->invalid_coordinate();
auto [first, second] = self->end_positions();
detail::matrix_coordinate const end_positions{detail::row_index_type{second}, detail::column_index_type{first}};
auto trace_path = self->trace_matrix().trace_path(end_positions);
auto trace_path_it = std::ranges::begin(trace_path);
std::ranges::advance(trace_path_it, std::ranges::end(trace_path));
detail::matrix_coordinate const begin_positions = trace_path_it.coordinate();
return {detail::column_index_type{begin_positions.col}, detail::row_index_type{begin_positions.row}};
}
auto alignment() const noexcept
{
using alignment_t = std::remove_cvref_t<decltype(std::declval<alignment_result_type &>().alignment())>;
derived_t const * self = static_cast<derived_t const *>(this);
static_assert(edit_traits::compute_sequence_alignment,
"alignment() can only be computed if you specify the "
"result type within your alignment config.");
if (!self->is_valid())
return alignment_t{};
auto [first, second] = self->end_positions();
detail::matrix_coordinate const end_positions{detail::row_index_type{second}, detail::column_index_type{first}};
aligned_sequence_builder builder{self->database, self->query};
auto trace_path = self->trace_matrix().trace_path(end_positions);
return builder(trace_path).alignment;
}
};
template <typename value_t>
class proxy_reference
{
public:
proxy_reference() noexcept = default;
proxy_reference(proxy_reference const &) noexcept = default;
proxy_reference(proxy_reference &&) noexcept = default;
proxy_reference & operator=(proxy_reference const &) noexcept = default;
proxy_reference & operator=(proxy_reference &&) noexcept = default;
~proxy_reference() = default;
proxy_reference(value_t & t) noexcept : ptr(std::addressof(t))
{}
proxy_reference(value_t &&) = delete;
template <typename other_value_t>
requires std::convertible_to<other_value_t, value_t>
proxy_reference & operator=(other_value_t && u) noexcept
{
get() = std::forward<other_value_t>(u);
return *this;
}
value_t & get() const noexcept
{
assert(ptr != nullptr);
return *ptr;
}
operator value_t &() const noexcept
{
return get();
}
private:
value_t * ptr{nullptr};
};
template <std::ranges::viewable_range database_t,
std::ranges::viewable_range query_t,
typename align_config_t,
typename edit_traits>
class edit_distance_unbanded :
// Hide this section in doxygen, because it messes up the inheritance.
public edit_distance_base<edit_traits::use_max_errors,
edit_distance_unbanded_max_errors_policy,
edit_traits,
edit_distance_unbanded<database_t, query_t, align_config_t, edit_traits>>,
public edit_distance_base<edit_traits::is_global,
edit_distance_unbanded_global_policy,
edit_traits,
edit_distance_unbanded<database_t, query_t, align_config_t, edit_traits>>,
public edit_distance_base<edit_traits::is_semi_global,
edit_distance_unbanded_semi_global_policy,
edit_traits,
edit_distance_unbanded<database_t, query_t, align_config_t, edit_traits>>,
public edit_distance_base<edit_traits::compute_score_matrix,
edit_distance_unbanded_score_matrix_policy,
edit_traits,
edit_distance_unbanded<database_t, query_t, align_config_t, edit_traits>>,
public edit_distance_base<edit_traits::compute_trace_matrix,
edit_distance_unbanded_trace_matrix_policy,
edit_traits,
edit_distance_unbanded<database_t, query_t, align_config_t, edit_traits>>
{
public:
using edit_traits::word_size;
using typename edit_traits::align_config_type;
using typename edit_traits::database_type;
using typename edit_traits::query_type;
using typename edit_traits::score_type;
using typename edit_traits::word_type;
private:
template <typename other_derived_t, typename other_edit_traits>
friend class edit_distance_unbanded_max_errors_policy;
template <typename other_derived_t, typename other_edit_traits>
friend class edit_distance_unbanded_global_policy;
template <typename other_derived_t, typename other_edit_traits>
friend class edit_distance_unbanded_semi_global_policy;
template <typename other_derived_t, typename other_edit_traits>
friend class edit_distance_unbanded_score_matrix_policy;
template <typename other_derived_t, typename other_edit_traits>
friend class edit_distance_unbanded_trace_matrix_policy;
using edit_traits::compute_begin_positions;
using edit_traits::compute_end_positions;
using edit_traits::compute_matrix;
using edit_traits::compute_score;
using edit_traits::compute_score_matrix;
using edit_traits::compute_sequence_alignment;
using edit_traits::compute_trace_matrix;
using edit_traits::is_global;
using edit_traits::is_semi_global;
using edit_traits::use_max_errors;
using typename edit_traits::alignment_result_type;
using typename edit_traits::database_iterator;
using typename edit_traits::query_alphabet_type;
using typename edit_traits::score_matrix_type;
using typename edit_traits::trace_matrix_type;
database_t database;
query_t query;
align_config_t config;
static constexpr word_type hp0 = is_global ? 1u : 0u;
static constexpr word_type hn0 = 0u;
static constexpr word_type vp0 = ~word_type{0u};
static constexpr word_type vn0 = 0u;
score_type _score{};
word_type score_mask{0u};
std::vector<word_type> vp{};
std::vector<word_type> vn{};
std::vector<word_type> bit_masks{};
database_iterator database_it{};
database_iterator database_it_end{};
struct compute_state_trace_matrix
{
proxy_reference<word_type> db{};
};
struct compute_state : enable_state_t<edit_traits::compute_trace_matrix, compute_state_trace_matrix>
{
using hp_type = std::conditional_t<edit_traits::compute_trace_matrix, proxy_reference<word_type>, word_type>;
word_type b{};
word_type d0{};
hp_type hp{};
word_type hn{};
proxy_reference<word_type> vp{};
proxy_reference<word_type> vn{};
word_type carry_d0{};
word_type carry_hp{hp0};
word_type carry_hn{};
};
void add_state()
{
if constexpr (!use_max_errors && compute_score_matrix)
this->_score_matrix.add_column(vp, vn);
if constexpr (!use_max_errors && compute_trace_matrix)
this->_trace_matrix.add_column(this->hp, this->db, vp);
if constexpr (use_max_errors && compute_matrix)
{
size_t max_rows = this->max_rows(score_mask, this->last_block, _score, this->max_errors);
if constexpr (compute_score_matrix)
this->_score_matrix.add_column(vp, vn, max_rows);
if constexpr (compute_trace_matrix)
this->_trace_matrix.add_column(this->hp, this->db, vp, max_rows);
}
}
public:
edit_distance_unbanded() = delete;
edit_distance_unbanded(edit_distance_unbanded const &) = default;
edit_distance_unbanded(edit_distance_unbanded &&) = default;
edit_distance_unbanded & operator=(edit_distance_unbanded const &) = default;
edit_distance_unbanded & operator=(edit_distance_unbanded &&) = default;
~edit_distance_unbanded() = default;
edit_distance_unbanded(database_t _database,
query_t _query,
align_config_t _config,
edit_traits const & SEQAN3_DOXYGEN_ONLY(_traits)) :
database{std::forward<database_t>(_database)},
query{std::forward<query_t>(_query)},
config{std::forward<align_config_t>(_config)},
_score{static_cast<score_type>(std::ranges::size(query))},
database_it{std::ranges::begin(database)},
database_it_end{std::ranges::end(database)}
{
static constexpr size_t alphabet_size_ = alphabet_size<query_alphabet_type>;
size_t const block_count = (std::ranges::size(query) - 1u + word_size) / word_size;
score_mask = word_type{1u} << ((std::ranges::size(query) - 1u + word_size) % word_size);
this->score_init();
if constexpr (use_max_errors)
this->max_errors_init(block_count);
if constexpr (compute_score_matrix)
this->score_matrix_init();
if constexpr (compute_trace_matrix)
this->trace_matrix_init(block_count);
vp.resize(block_count, vp0);
vn.resize(block_count, vn0);
bit_masks.resize((alphabet_size_ + 1u) * block_count, 0u);
// encoding the letters as bit-vectors
for (size_t j = 0u; j < std::ranges::size(query); j++)
{
size_t const i = block_count * seqan3::to_rank(query[j]) + j / word_size;
bit_masks[i] |= word_type{1u} << (j % word_size);
}
add_state();
}
private:
template <bool with_carry>
static void compute_step(compute_state & state) noexcept
{
word_type x, t;
assert(state.carry_d0 <= 1u);
assert(state.carry_hp <= 1u);
assert(state.carry_hn <= 1u);
x = state.b | state.vn;
t = state.vp + (x & state.vp) + state.carry_d0;
state.d0 = (t ^ state.vp) | x;
state.hn = state.vp & state.d0;
state.hp = state.vn | ~(state.vp | state.d0);
if constexpr (with_carry)
state.carry_d0 = (state.carry_d0 != 0u) ? t <= state.vp : t < state.vp;
x = (state.hp << 1u) | state.carry_hp;
state.vn = x & state.d0;
state.vp = (state.hn << 1u) | ~(x | state.d0) | state.carry_hn;
if constexpr (with_carry)
{
state.carry_hp = state.hp >> (word_size - 1u);
state.carry_hn = state.hn >> (word_size - 1u);
}
}
template <bool with_carry>
void compute_kernel(compute_state & state, size_t const block_offset, size_t const current_block) noexcept
{
state.vp = proxy_reference<word_type>{this->vp[current_block]};
state.vn = proxy_reference<word_type>{this->vn[current_block]};
if constexpr (compute_trace_matrix)
{
state.hp = proxy_reference<word_type>{this->hp[current_block]};
state.db = proxy_reference<word_type>{this->db[current_block]};
}
state.b = bit_masks[block_offset + current_block];
compute_step<with_carry>(state);
if constexpr (compute_trace_matrix)
state.db = ~(state.b ^ state.d0);
}
void advance_score(word_type P, word_type N, word_type mask) noexcept
{
if ((P & mask) != word_type{0u})
_score++;
else if ((N & mask) != word_type{0u})
_score--;
}
bool on_hit() noexcept
{
// TODO: call external on_hit functor
return false;
}
inline bool small_patterns();
inline bool large_patterns();
inline void compute_empty_query_sequence()
{
assert(std::ranges::empty(query));
bool abort_computation = false;
for (; database_it != database_it_end; ++database_it)
{
if constexpr (is_global)
++_score;
else // is_semi_global
this->update_best_score();
// call on_hit
if constexpr (use_max_errors)
abort_computation = on_hit();
this->add_state();
if (abort_computation)
break;
}
}
void compute()
{
// limit search width for prefix search (if no matrix needs to be computed)
if constexpr (use_max_errors && is_global && !compute_matrix)
{
// Note: For global alignments we know that the database can only be max_length long to have a score less
// than or equal max_errors in the last cell.
//
// The following matrix shows a minimal value for each entry (because a diagonal must be +0 or +1, each
// diagonal is at least the value of the initial value in the first row)
// 0 1 2 3 4 5 6...
// 1 0 1 2 3 4 5...
// ...
// m ... 3 2 1 0 1 2 3 4 5
// Thus, after |query| + max_errors entries the score will always be higher than max_errors.
size_t const max_length = std::ranges::size(query) + this->max_errors + 1u;
size_t const haystack_length = std::min(std::ranges::size(database), max_length);
database_it_end -= std::ranges::size(database) - haystack_length;
}
// distinguish between the version for needles not longer than
// one machine word and the version for longer needles
// A special cases is if the second sequence is empty (vp.size() == 0u).
if (vp.size() == 0u) // [[unlikely]]
compute_empty_query_sequence();
else if (vp.size() == 1u)
small_patterns();
else
large_patterns();
if constexpr (is_global)
this->update_best_score();
}
public:
template <typename callback_t>
void operator()([[maybe_unused]] size_t const idx, callback_t && callback)
{
using traits_type = alignment_configuration_traits<align_config_t>;
using result_value_type = typename alignment_result_value_type_accessor<alignment_result_type>::type;
compute();
// First cache the begin and end positions if enabled by the edit distance traits.
// Note, that they might be activated even if the user did not configure them, but in order to
// compute for example the alignment both information are required and enabled internally.
auto cached_end_positions = this->invalid_coordinate();
auto cached_begin_positions = this->invalid_coordinate();
if constexpr (compute_end_positions)
cached_end_positions = this->end_positions();
if constexpr (compute_begin_positions && !compute_sequence_alignment)
{
static_assert(compute_end_positions, "End positions required to compute the begin positions.");
cached_begin_positions = this->begin_positions();
}
// Fill the result value type. Note we need to ask what was enabled on the user side in order to store
// the correct information in the alignment result type. This information is stored in the alignment
// configuration traits and not in the edit distance traits.
result_value_type res_vt{};
if constexpr (traits_type::output_sequence1_id)
res_vt.sequence1_id = idx;
if constexpr (traits_type::output_sequence2_id)
res_vt.sequence2_id = idx;
if constexpr (traits_type::compute_score)
res_vt.score = this->score().value_or(matrix_inf<score_type>);
if constexpr (traits_type::compute_sequence_alignment)
{
if (this->is_valid())
{
auto [first, second] = cached_end_positions;
detail::matrix_coordinate const end_positions{detail::row_index_type{second},
detail::column_index_type{first}};
aligned_sequence_builder builder{this->database, this->query};
auto trace_res = builder(this->trace_matrix().trace_path(end_positions));
res_vt.alignment = std::move(trace_res.alignment);
cached_begin_positions.first = trace_res.first_sequence_slice_positions.first;
cached_begin_positions.second = trace_res.second_sequence_slice_positions.first;
}
}
if constexpr (traits_type::compute_end_positions)
res_vt.end_positions = std::move(cached_end_positions);
if constexpr (traits_type::compute_begin_positions)
res_vt.begin_positions = std::move(cached_begin_positions);
callback(alignment_result_type{std::move(res_vt)});
}
};
template <std::ranges::viewable_range database_t,
std::ranges::viewable_range query_t,
typename align_config_t,
typename traits_t>
bool edit_distance_unbanded<database_t, query_t, align_config_t, traits_t>::small_patterns()
{
bool abort_computation = false;
// computing the blocks
while (database_it != database_it_end)
{
compute_state state{};
size_t const block_offset = seqan3::to_rank((query_alphabet_type)*database_it);
compute_kernel<false>(state, block_offset, 0u);
advance_score(state.hp, state.hn, score_mask);
// semi-global without max_errors guarantees that the score stays within the last row
if constexpr (is_semi_global && !use_max_errors)
this->update_best_score();
// updating the last active cell
if constexpr (use_max_errors)
abort_computation = this->update_last_active_cell();
add_state();
++database_it;
if (abort_computation)
return true;
}
return false;
}
template <std::ranges::viewable_range database_t,
std::ranges::viewable_range query_t,
typename align_config_t,
typename traits_t>
bool edit_distance_unbanded<database_t, query_t, align_config_t, traits_t>::large_patterns()
{
bool abort_computation = false;
while (database_it != database_it_end)
{
compute_state state{};
size_t const block_offset = vp.size() * seqan3::to_rank((query_alphabet_type)*database_it);
size_t block_count = vp.size();
if constexpr (use_max_errors)
block_count = this->last_block + 1;
// compute each block in the current column; carries between blocks will be propagated.
for (size_t current_block = 0u; current_block < block_count; current_block++)
compute_kernel<true>(state, block_offset, current_block);
advance_score(state.hp, state.hn, score_mask);
// semi-global without max_errors guarantees that the score stays within the last row
if constexpr (is_semi_global && !use_max_errors)
this->update_best_score();
if constexpr (use_max_errors)
{
// if the last active cell reached the end within the current block we have to compute the next block.
bool additional_block = score_mask >> (word_size - 1u);
bool reached_last_block = this->last_block + 1u == vp.size();
// If there is no next block we skip the computation.
if (reached_last_block)
additional_block = false;
if (additional_block)
{
size_t const current_block = this->last_block + 1u;
// this might not be necessary, but carry_d0 = 1u might have an influence on the result of vn and vp.
vp[current_block] = vp0;
vn[current_block] = vn0;
compute_kernel<false>(state, block_offset, current_block);
}
// updating the last active cell
abort_computation = this->update_last_active_cell();
}
add_state();
++database_it;
if (abort_computation)
return true;
}
return false;
}
template <typename database_t, typename query_t, typename config_t, typename traits_t>
edit_distance_unbanded(database_t && database, query_t && query, config_t config, traits_t)
-> edit_distance_unbanded<database_t, query_t, config_t, traits_t>;
} // namespace seqan3::detail
advanceable_alignment_coordinate.hpp
Provides seqan3::detail::advanceable_alignment_coordinate.
algorithm
alignment_result.hpp
Provides seqan3::alignment_result.
std::begin
T begin(T... args)
bitset
configuration.hpp
Provides seqan3::configuration and utility functions.
edit_distance_fwd.hpp
Forwards for seqan3::edit_distance_unbanded related types.
edit_distance_score_matrix_full.hpp
Provides seqan3::detail::edit_distance_score_matrix_full.
edit_distance_trace_matrix_full.hpp
Provides seqan3::detail::edit_distance_trace_matrix_full.
std::end
T end(T... args)
std::forward
T forward(T... args)
seqan3::alphabet_size
constexpr auto alphabet_size
A type trait that holds the size of a (semi-)alphabet.
Definition alphabet/concept.hpp:834
seqan3::to_rank
constexpr auto to_rank
Return the rank representation of a (semi-)alphabet object.
Definition alphabet/concept.hpp:152
seqan3::field::offset
@ offset
Sequence (seqan3::field::seq) relative start position (0-based), unsigned value.
Definition record.hpp:66
seqan3::pack_traits::size
constexpr size_t size
The size of a type pack.
Definition type_pack/traits.hpp:143
matrix_concept.hpp
Provides seqan3::detail::matrix.
std::min
T min(T... args)
seqan3::get
constexpr auto const & get(configuration< configs_t... > const &config) noexcept
This is an overloaded member function, provided for convenience. It differs from the above function o...
Definition configuration.hpp:412
std::vector::size
T size(T... args)
utility
Hide me