policy_affine_gap_recursion.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2021, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2021, Knut Reinert & MPI für molekulare Genetik
// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
// shipped with this file and also available at: https://github.com/seqan/seqan3/blob/master/LICENSE.md
// -----------------------------------------------------------------------------------------------------
 
#pragma once
 
#include <tuple>
 
#include <seqan3/alignment/configuration/align_config_gap_cost_affine.hpp>
#include <seqan3/alignment/matrix/detail/affine_cell_proxy.hpp>
#include <seqan3/alignment/pairwise/detail/type_traits.hpp>
#include <seqan3/core/detail/template_inspection.hpp>
#include <seqan3/utility/simd/algorithm.hpp>
#include <seqan3/utility/simd/concept.hpp>
 
namespace seqan3::detail
{
 
template <typename alignment_configuration_t>
class policy_affine_gap_recursion
{
protected:
    using traits_type = alignment_configuration_traits<alignment_configuration_t>;
    using original_score_type = typename traits_type::original_score_type;
    using score_type = typename traits_type::score_type;
    using affine_score_tuple_t = std::tuple<score_type, score_type, score_type>;
    using affine_cell_type = affine_cell_proxy<affine_score_tuple_t>;
 
    score_type gap_extension_score{};
    score_type gap_open_score{};
 
    bool first_row_is_free{};
    bool first_column_is_free{};
 
    policy_affine_gap_recursion() = default; 
    policy_affine_gap_recursion(policy_affine_gap_recursion const &) = default; 
    policy_affine_gap_recursion(policy_affine_gap_recursion &&) = default; 
    policy_affine_gap_recursion & operator=(policy_affine_gap_recursion const &) = default; 
    policy_affine_gap_recursion & operator=(policy_affine_gap_recursion &&) = default; 
    ~policy_affine_gap_recursion() = default; 
 
    explicit policy_affine_gap_recursion(alignment_configuration_t const & config)
    {
        // Get the gap scheme from the config or choose -1 and -10 as default.
        auto const & selected_gap_scheme = config.get_or(align_cfg::gap_cost_affine{align_cfg::open_score{-10},
                                                                                    align_cfg::extension_score{-1}});
 
        gap_extension_score = maybe_convert_to_simd(selected_gap_scheme.extension_score);
        gap_open_score = maybe_convert_to_simd(selected_gap_scheme.open_score) + gap_extension_score;
 
        auto method_global_config = config.get_or(align_cfg::method_global{});
        first_row_is_free = method_global_config.free_end_gaps_sequence1_leading;
        first_column_is_free = method_global_config.free_end_gaps_sequence2_leading;
    }
 
    template <typename affine_cell_t>
    affine_cell_type compute_inner_cell(score_type diagonal_score,
                                        affine_cell_t previous_cell,
                                        score_type const sequence_score) const noexcept
    {
        diagonal_score += sequence_score;
        score_type horizontal_score = previous_cell.horizontal_score();
        score_type vertical_score = previous_cell.vertical_score();
 
        diagonal_score = (diagonal_score < vertical_score) ? vertical_score : diagonal_score;
        diagonal_score = (diagonal_score < horizontal_score) ? horizontal_score : diagonal_score;
 
        score_type tmp = diagonal_score + gap_open_score;
        vertical_score += gap_extension_score;
        horizontal_score += gap_extension_score;
 
        // store the vertical_score and horizontal_score value in the next path
        vertical_score = (vertical_score < tmp) ? tmp : vertical_score;
        horizontal_score = (horizontal_score < tmp) ? tmp : horizontal_score;
 
        return {diagonal_score, horizontal_score, vertical_score};
    }
 
    affine_cell_type initialise_origin_cell() const noexcept
    {
        return {score_type{},
                first_row_is_free ? score_type{} : gap_open_score,
                first_column_is_free ? score_type{} : gap_open_score};
    }
 
    template <typename affine_cell_t>
    affine_cell_type initialise_first_column_cell(affine_cell_t previous_cell) const noexcept
    {
        score_type new_vertical = previous_cell.vertical_score() + gap_extension_score;
        return {previous_cell.vertical_score(),
                previous_cell.vertical_score() + gap_open_score,
                first_column_is_free ? previous_cell.vertical_score() : new_vertical};
    }
 
    template <typename affine_cell_t>
    affine_cell_type initialise_first_row_cell(affine_cell_t previous_cell) const noexcept
    {
        score_type new_horizontal_score = previous_cell.horizontal_score() + gap_extension_score;
        return {previous_cell.horizontal_score(),
                first_row_is_free ? previous_cell.horizontal_score() : new_horizontal_score,
                previous_cell.horizontal_score() + gap_open_score};
    }
 
    score_type lowest_viable_score() const noexcept
    {
        if constexpr (simd_concept<score_type>)
            assert(gap_open_score[0] <= 0 && gap_extension_score[0] <= 0);
        else
            assert(gap_open_score <= 0 && gap_extension_score <= 0);
 
        return maybe_convert_to_simd(std::numeric_limits<original_score_type>::lowest()) -
               (gap_open_score + gap_extension_score);
    }
 
    template <typename score_t>
        requires arithmetic<std::remove_cvref_t<score_t>>
    constexpr auto maybe_convert_to_simd(score_t && score) const noexcept
    {
        if constexpr (simd_concept<score_type>)
            return simd::fill<score_type>(std::forward<score_t>(score));
        else // Return unmodified.
            return std::forward<score_t>(score);
    }
};
} // namespace seqan3::detail