policy_alignment_matrix.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <tuple>
 
#include <seqan3/alignment/exception.hpp>
#include <seqan3/alignment/matrix/detail/coordinate_matrix.hpp>
#include <seqan3/alignment/pairwise/detail/type_traits.hpp>
#include <seqan3/core/configuration/configuration.hpp>
#include <seqan3/core/detail/template_inspection.hpp>
 
namespace seqan3::detail
{
 
template <typename traits_t, typename alignment_matrix_t>
    requires (is_type_specialisation_of_v<traits_t, alignment_configuration_traits> &&
              requires (alignment_matrix_t & matrix, typename traits_t::score_type const initial_score)
              {
                  { matrix.resize(column_index_type{size_t{}}, row_index_type{size_t{}}, initial_score) };
              })
class policy_alignment_matrix
{
protected:
    using score_type = typename traits_t::score_type;
    using matrix_index_type = typename traits_t::matrix_index_type;
 
    int32_t lower_diagonal{};
    int32_t upper_diagonal{};
    bool last_column_is_free{};
    bool last_row_is_free{};
 
    policy_alignment_matrix() = default;                                            
    policy_alignment_matrix(policy_alignment_matrix const &) = default;             
    policy_alignment_matrix(policy_alignment_matrix &&) = default;                  
    policy_alignment_matrix & operator=(policy_alignment_matrix const &) = default; 
    policy_alignment_matrix & operator=(policy_alignment_matrix &&) = default;      
    ~policy_alignment_matrix() = default;                                           
 
    template <typename alignment_configuration_t>
        requires (is_type_specialisation_of_v<alignment_configuration_t, configuration>)
    policy_alignment_matrix(alignment_configuration_t const & config)
    {
        using seqan3::get;
 
        auto band = config.get_or(seqan3::align_cfg::band_fixed_size{});
 
        lower_diagonal = band.lower_diagonal;
        upper_diagonal = band.upper_diagonal;
 
        bool invalid_band = upper_diagonal < lower_diagonal;
        std::string error_cause = (invalid_band) ? " The upper diagonal is smaller than the lower diagonal." : "";
 
        if constexpr (traits_t::is_global)
        {
            auto method_global_config = get<seqan3::align_cfg::method_global>(config);
 
            bool first_row_is_free = method_global_config.free_end_gaps_sequence1_leading;
            bool first_column_is_free = method_global_config.free_end_gaps_sequence2_leading;
 
            last_row_is_free = method_global_config.free_end_gaps_sequence1_trailing;
            last_column_is_free = method_global_config.free_end_gaps_sequence2_trailing;
            // band starts in first column without free gaps or band starts in first row without free gaps.
            invalid_band |= (upper_diagonal < 0 && !first_column_is_free) || (lower_diagonal > 0 && !first_row_is_free);
            error_cause += " The band starts in a region without free gaps.";
        }
 
        if (invalid_band)
            throw invalid_alignment_configuration{"The selected band [" + std::to_string(lower_diagonal) + ":"
                                                  + std::to_string(upper_diagonal)
                                                  + "] cannot be used with the current "
                                                    "alignment configuration:"
                                                  + error_cause};
    }
 
    auto acquire_matrices(size_t const sequence1_size,
                          size_t const sequence2_size,
                          score_type initial_score = score_type{}) const
    {
        assert(sequence1_size < static_cast<uint64_t>(std::numeric_limits<int64_t>::max()));
        assert(sequence2_size < static_cast<uint64_t>(std::numeric_limits<int64_t>::max()));
 
        if constexpr (traits_t::is_banded)
            check_valid_band_configuration(sequence1_size, sequence2_size);
 
        static thread_local alignment_matrix_t alignment_matrix{};
        static thread_local coordinate_matrix<matrix_index_type> index_matrix{};
 
        // Increase dimension by one for the initialisation of the matrix.
        size_t const column_count = sequence1_size + 1;
        size_t row_count = sequence2_size + 1;
 
        index_matrix.resize(column_index_type{column_count}, row_index_type{row_count});
 
        if constexpr (traits_t::is_banded)
        {
            assert(upper_diagonal - lower_diagonal + 1 > 0); // Band size is a positive integer.
            // Allocate one more cell to compute the last cell of the band with standard recursion function.
            row_count = std::min<int64_t>(upper_diagonal - lower_diagonal + 2, row_count);
        }
 
        alignment_matrix.resize(column_index_type{column_count}, row_index_type{row_count}, initial_score);
 
        return std::tie(alignment_matrix, index_matrix);
    }
 
    void check_valid_band_configuration(size_t const sequence1_size, size_t const sequence2_size) const
    {
        bool const upper_diagonal_ends_before_last_cell = (upper_diagonal + sequence2_size) < sequence1_size;
        bool const lower_diagonal_ends_behind_last_cell = (-lower_diagonal + sequence1_size) < sequence2_size;
 
        bool invalid_band = false;
        std::string error_cause{};
 
        if constexpr (traits_t::is_global)
        {
            // band ends in last column without free gaps or band ends in last row without free gaps.
            invalid_band |= (lower_diagonal_ends_behind_last_cell && !last_column_is_free)
                         || (upper_diagonal_ends_before_last_cell && !last_row_is_free);
            error_cause = "The band ends in a region without free gaps.";
        }
 
        if (invalid_band)
            throw invalid_alignment_configuration{"The selected band [" + std::to_string(lower_diagonal) + ":"
                                                  + std::to_string(upper_diagonal)
                                                  + "] cannot be used with the current "
                                                    "alignment configuration: "
                                                  + error_cause};
    }
};
} // namespace seqan3::detail