Loading...
Searching...
No Matches
edit_distance_trace_matrix_full.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 <bitset>
#include <seqan3/alignment/matrix/detail/trace_directions.hpp>
#include <seqan3/alignment/pairwise/edit_distance_fwd.hpp>
#include <seqan3/utility/detail/bits_of.hpp>
namespace seqan3::detail
{
template <typename word_t, bool is_semi_global, bool use_max_errors>
class edit_distance_trace_matrix_full
{
public:
template <std::ranges::viewable_range database_t,
std::ranges::viewable_range query_t,
typename align_config_t,
typename edit_traits>
friend class edit_distance_unbanded;
edit_distance_trace_matrix_full() = default;
edit_distance_trace_matrix_full(edit_distance_trace_matrix_full const &) = default;
edit_distance_trace_matrix_full(edit_distance_trace_matrix_full &&) = default;
edit_distance_trace_matrix_full & operator=(edit_distance_trace_matrix_full const &) = default;
edit_distance_trace_matrix_full & operator=(edit_distance_trace_matrix_full &&) = default;
~edit_distance_trace_matrix_full() = default;
protected:
template <typename derived_t, typename edit_traits>
friend class edit_distance_unbanded_trace_matrix_policy;
edit_distance_trace_matrix_full(size_t const rows_size) : rows_size{rows_size}, columns{}
{}
private:
struct trace_path_iterator;
public:
using word_type = word_t;
static constexpr auto word_size = bits_of<word_type>;
using value_type = detail::trace_directions;
using reference = value_type;
using size_type = size_t;
void reserve(size_t const new_capacity)
{
columns.reserve(new_capacity);
}
public:
reference at(matrix_coordinate const & coordinate) const noexcept
{
size_t row = coordinate.row;
size_t col = coordinate.col;
assert(row < rows());
assert(col < cols());
column_type const & column = columns[col];
if constexpr (use_max_errors)
if (!(row < column.max_rows))
return detail::trace_directions::none;
if (row == 0u)
{
if constexpr (is_semi_global)
return detail::trace_directions::none;
if (col == 0u)
return detail::trace_directions::none;
return detail::trace_directions::left;
}
size_t const idx = (row - 1u) / word_size;
size_t const offset = (row - 1u) % word_size;
bool const left = std::bitset<word_size>(column.left[idx])[offset];
bool const diagonal = std::bitset<word_size>(column.diagonal[idx])[offset];
bool const up = std::bitset<word_size>(column.up[idx])[offset];
auto const dir = (left ? detail::trace_directions::left : detail::trace_directions::none)
| (diagonal ? detail::trace_directions::diagonal : detail::trace_directions::none)
| (up ? detail::trace_directions::up : detail::trace_directions::none);
return dir;
}
size_t rows() const noexcept
{
return rows_size;
}
size_t cols() const noexcept
{
return columns.size();
}
auto trace_path(matrix_coordinate const & trace_begin) const
{
if (trace_begin.row >= rows() || trace_begin.col >= cols())
throw std::invalid_argument{"The given coordinate exceeds the matrix in vertical or horizontal direction."};
using path_t = std::ranges::subrange<trace_path_iterator, std::default_sentinel_t>;
return path_t{trace_path_iterator{this, trace_begin}, std::default_sentinel};
}
protected:
struct max_errors_state
{
size_t max_rows{};
};
struct trace_matrix_state
{
std::vector<word_type> left{};
std::vector<word_type> diagonal{};
std::vector<word_type> up{};
};
struct column_type : enable_state_t<true, trace_matrix_state>, enable_state_t<use_max_errors, max_errors_state>
{};
void add_column(std::vector<word_type> left, std::vector<word_type> diagonal, std::vector<word_type> up)
requires (!use_max_errors)
{
column_type column{};
column.left = std::move(left);
column.diagonal = std::move(diagonal);
column.up = std::move(up);
columns.push_back(std::move(column));
}
void add_column(std::vector<word_type> left,
std::vector<word_type> diagonal,
std::vector<word_type> up,
size_t const max_rows)
requires use_max_errors
{
column_type column{};
column.left = std::move(left);
column.diagonal = std::move(diagonal);
column.up = std::move(up);
column.max_rows = max_rows;
columns.push_back(std::move(column));
}
private:
size_t rows_size{};
std::vector<column_type> columns{};
};
template <typename word_t, bool is_semi_global, bool use_max_errors>
struct edit_distance_trace_matrix_full<word_t, is_semi_global, use_max_errors>::trace_path_iterator
{
using iterator_category = std::input_iterator_tag;
using value_type = detail::trace_directions;
using difference_type = std::ptrdiff_t;
static constexpr value_type D = value_type::diagonal;
static constexpr value_type L = value_type::left;
static constexpr value_type U = value_type::up;
static constexpr value_type N = value_type::none;
constexpr value_type operator*() const
{
value_type dir = parent->at(coordinate());
if (dir == N)
return N;
if ((dir & L) == L)
return L;
else if ((dir & U) == U)
return U;
else
return D;
}
[[nodiscard]] constexpr matrix_coordinate const & coordinate() const
{
return coordinate_;
}
constexpr trace_path_iterator & operator++()
{
value_type dir = *(*this);
if ((dir & L) == L)
{
coordinate_.col = std::max<size_t>(coordinate_.col, 1) - 1;
}
else if ((dir & U) == U)
{
coordinate_.row = std::max<size_t>(coordinate_.row, 1) - 1;
}
else if ((dir & D) == D)
{
coordinate_.row = std::max<size_t>(coordinate_.row, 1) - 1;
coordinate_.col = std::max<size_t>(coordinate_.col, 1) - 1;
}
// seqan3::trace_direction::none in an inner cell of the trace matrix is not possible in
// non-local alignments, e.g. global, semi-global, max-error, banded. On the other hand,
// this can happen in a cell of the first row or first colomn.
assert(dir != N || coordinate_.row == 0 || coordinate_.col == 0);
return *this;
}
constexpr void operator++(int)
{
++(*this);
}
friend bool operator==(trace_path_iterator const & it, std::default_sentinel_t)
{
return *it == value_type::none;
}
friend bool operator==(std::default_sentinel_t, trace_path_iterator const & it)
{
return it == std::default_sentinel;
}
friend bool operator!=(trace_path_iterator const & it, std::default_sentinel_t)
{
return !(it == std::default_sentinel);
}
friend bool operator!=(std::default_sentinel_t, trace_path_iterator const & it)
{
return it != std::default_sentinel;
}
edit_distance_trace_matrix_full const * parent{nullptr};
matrix_coordinate coordinate_{};
};
} // namespace seqan3::detail
bits_of.hpp
Provides utility functions for bit twiddling.
bitset
edit_distance_fwd.hpp
Forwards for seqan3::edit_distance_unbanded related types.
seqan3::field::offset
@ offset
Sequence (seqan3::field::seq) relative start position (0-based), unsigned value.
Definition record.hpp:66
std::left
T left(T... args)
std::max
T max(T... args)
std::rel_ops::operator!=
T operator!=(T... args)
trace_directions.hpp
Provides the declaration of seqan3::detail::trace_directions.
Hide me