alignment_configurator.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2020, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2020, 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 <functional>
#include <tuple>
#include <utility>
#include <vector>
 
#include <seqan3/alignment/configuration/all.hpp>
#include <seqan3/alignment/matrix/detail/alignment_score_matrix_one_column.hpp>
#include <seqan3/alignment/matrix/detail/alignment_score_matrix_one_column_banded.hpp>
#include <seqan3/alignment/matrix/detail/alignment_trace_matrix_full.hpp>
#include <seqan3/alignment/matrix/detail/alignment_trace_matrix_full_banded.hpp>
#include <seqan3/alignment/pairwise/policy/all.hpp>
#include <seqan3/alignment/pairwise/alignment_algorithm.hpp>
#include <seqan3/alignment/pairwise/align_result_selector.hpp>
#include <seqan3/alignment/pairwise/alignment_result.hpp>
#include <seqan3/alignment/pairwise/detail/type_traits.hpp>
#include <seqan3/alignment/pairwise/detail/concept.hpp>
#include <seqan3/alignment/pairwise/edit_distance_algorithm.hpp>
#include <seqan3/alignment/scoring/detail/simd_match_mismatch_scoring_scheme.hpp>
#include <seqan3/core/concept/tuple.hpp>
#include <seqan3/core/simd/simd.hpp>
#include <seqan3/core/type_traits/deferred_crtp_base.hpp>
#include <seqan3/core/type_traits/lazy.hpp>
#include <seqan3/core/type_traits/template_inspection.hpp>
#include <seqan3/core/type_list/type_list.hpp>
#include <seqan3/range/views/type_reduce.hpp>
#include <seqan3/range/views/zip.hpp>
 
namespace seqan3::detail
{
 
template <typename config_t>
    requires is_type_specialisation_of_v<config_t, configuration> && config_t::template exists<align_cfg::result>()
struct align_config_result_score
{
private:
    using result_config_t = std::remove_reference_t<decltype(seqan3::get<align_cfg::result>(std::declval<config_t>()))>;
 
public:
    using type = typename result_config_t::score_type;
};
 
template <typename range_type,
          typename alignment_config_type>
struct alignment_contract
{
private:
     using unref_range_type = std::remove_reference_t<range_type>;
    using first_seq_t  = std::tuple_element_t<0, value_type_t<std::ranges::iterator_t<unref_range_type>>>;
    using second_seq_t = std::tuple_element_t<1, value_type_t<std::ranges::iterator_t<unref_range_type>>>;
 
public:
    constexpr static bool expects_tuple_like_value_type()
    {
        return tuple_like<alignment_config_type> &&
               std::tuple_size_v<value_type_t<std::ranges::iterator_t<unref_range_type>>> == 2;
    }
 
    constexpr static bool expects_valid_scoring_scheme()
    {
        if constexpr (alignment_config_type::template exists<align_cfg::scoring>())
        {
            using scoring_type = std::remove_reference_t<
                                    decltype(get<align_cfg::scoring>(std::declval<alignment_config_type>()).value)
                                 >;
            return static_cast<bool>(scoring_scheme<scoring_type,
                                                   value_type_t<first_seq_t>,
                                                   value_type_t<second_seq_t>>);
        }
        else
        {
            return false;
        }
    }
 
    constexpr static bool expects_alignment_configuration()
    {
        return alignment_config_type::template exists<align_cfg::mode>();
    }
};
 
struct alignment_configurator
{
private:
 
    template <typename traits_t>
    struct select_matrix_policy
    {
    private:
        static constexpr bool only_coordinates = traits_t::result_type_rank < with_front_coordinate_type::rank;
 
        using score_matrix_t = std::conditional_t<traits_t::is_banded,
                                                  alignment_score_matrix_one_column_banded<typename traits_t::score_t>,
                                                  alignment_score_matrix_one_column<typename traits_t::score_t>>;
        using trace_matrix_t = std::conditional_t<traits_t::is_banded,
                                                  alignment_trace_matrix_full_banded<typename traits_t::trace_t,
                                                                                     only_coordinates>,
                                                  alignment_trace_matrix_full<typename traits_t::trace_t,
                                                                              only_coordinates>>;
 
    public:
        using type = deferred_crtp_base<alignment_matrix_policy, score_matrix_t, trace_matrix_t>;
    };
 
    template <typename traits_t>
    struct select_gap_policy
    {
    private:
        using score_t = typename traits_t::score_t;
        using is_local_t = std::bool_constant<traits_t::is_local>;
 
    public:
        using type = std::conditional_t<traits_t::is_vectorised,
                                        deferred_crtp_base<simd_affine_gap_policy, score_t, is_local_t>,
                                        deferred_crtp_base<affine_gap_policy, score_t, is_local_t>>;
    };
 
    template <typename traits_t, typename policy_traits_t>
    struct select_find_optimum_policy
    {
    private:
        using score_t = typename traits_t::score_t;
        static constexpr bool is_global_alignment = traits_t::is_global && !traits_t::is_aligned_ends;
 
    public:
        using type = std::conditional_t<traits_t::is_vectorised,
                                        deferred_crtp_base<simd_find_optimum_policy,
                                                           score_t,
                                                           std::bool_constant<is_global_alignment>,
                                                           policy_traits_t>,
                                        deferred_crtp_base<find_optimum_policy, policy_traits_t>>;
    };
 
public:
    template <align_pairwise_range_input sequences_t, typename config_t>
        requires is_type_specialisation_of_v<config_t, configuration>
    static constexpr auto configure(config_t const & cfg)
    {
 
        if constexpr (!config_t::template exists<align_cfg::result>())
        {
            // Set the default result value to be computed.
            return configure<sequences_t>(cfg | align_cfg::result{with_score});
        }
        else
        {
            // ----------------------------------------------------------------------------
            // Configure the type-erased alignment function.
            // ----------------------------------------------------------------------------
 
            using first_seq_t = std::tuple_element_t<0, value_type_t<sequences_t>>;
            using second_seq_t = std::tuple_element_t<1, value_type_t<sequences_t>>;
 
            using wrapped_first_t  = type_reduce_view<first_seq_t &>;
            using wrapped_second_t = type_reduce_view<second_seq_t &>;
 
            // The alignment executor passes a chunk over an indexed sequence pair range to the alignment algorithm.
            using indexed_sequence_pair_range_t = typename chunked_indexed_sequence_pairs<sequences_t>::type;
            using indexed_sequence_pair_chunk_t = std::ranges::range_value_t<indexed_sequence_pair_range_t>;
 
            // Select the result type based on the sequences and the configuration.
            using result_t = alignment_result<typename align_result_selector<std::remove_reference_t<wrapped_first_t>,
                                                                             std::remove_reference_t<wrapped_second_t>,
                                                                             config_t
                                                                            >::type
                                             >;
            using result_collection_t = std::vector<result_t>;  // Use a vector as return type.
            // Define the function wrapper type.
            using function_wrapper_t = std::function<result_collection_t(indexed_sequence_pair_chunk_t)>;
 
            // ----------------------------------------------------------------------------
            // Test some basic preconditions
            // ----------------------------------------------------------------------------
 
            using alignment_contract_t = alignment_contract<sequences_t, config_t>;
 
            static_assert(alignment_contract_t::expects_alignment_configuration(),
                          "Alignment configuration error: "
                          "The alignment can only be configured with alignment configurations.");
 
            static_assert(alignment_contract_t::expects_tuple_like_value_type(),
                          "Alignment configuration error: "
                          "The value type of the sequence ranges must model the seqan3::tuple_like "
                          "and must contain exactly 2 elements.");
 
            static_assert(alignment_contract_t::expects_valid_scoring_scheme(),
                          "Alignment configuration error: "
                          "Either the scoring scheme was not configured or the given scoring scheme cannot be invoked with "
                          "the value types of the passed sequences.");
 
            // ----------------------------------------------------------------------------
            // Configure the algorithm
            // ----------------------------------------------------------------------------
 
            // Use default edit distance if gaps are not set.
            auto const & gaps = cfg.template value_or<align_cfg::gap>(gap_scheme{gap_score{-1}});
            auto const & scoring_scheme = get<align_cfg::scoring>(cfg).value;
            auto align_ends_cfg = cfg.template value_or<align_cfg::aligned_ends>(free_ends_none);
 
            if constexpr (config_t::template exists<align_cfg::mode<detail::global_alignment_type>>())
            {
                // Only use edit distance if ...
                if (gaps.get_gap_open_score() == 0 &&  // gap open score is not set,
                    !(align_ends_cfg[2] || align_ends_cfg[3]) && // none of the free end gaps are set for second seq,
                    align_ends_cfg[0] == align_ends_cfg[1]) // free ends for leading and trailing gaps are equal in first seq.
                {
                    // TODO: Instead of relying on nucleotide scoring schemes we need to be able to determine the edit distance
                    //       option via the scheme.
                    if constexpr (is_type_specialisation_of_v<remove_cvref_t<decltype(scoring_scheme)>,
                                                              nucleotide_scoring_scheme>)
                    {
                        if ((scoring_scheme.score('A'_dna15, 'A'_dna15) == 0) &&
                            (scoring_scheme.score('A'_dna15, 'C'_dna15)) == -1)
                            return std::pair{configure_edit_distance<function_wrapper_t>(cfg), cfg};
                    }
                }
            }
 
            // ----------------------------------------------------------------------------
            // Check if invalid configuration was used.
            // ----------------------------------------------------------------------------
 
            // Do not allow max error configuration for alignments not computing the edit distance.
            if (config_t::template exists<align_cfg::max_error>())
                throw invalid_alignment_configuration{"The align_cfg::max_error configuration is only allowed for "
                                                      "the specific edit distance computation."};
            // Configure the alignment algorithm.
            return std::pair{configure_scoring_scheme<function_wrapper_t>(cfg), cfg};
        }
    }
 
private:
    template <typename function_wrapper_t, typename config_t>
    static constexpr function_wrapper_t configure_edit_distance(config_t const & cfg)
    {
        using traits_t = alignment_configuration_traits<config_t>;
 
        // ----------------------------------------------------------------------------
        // Unsupported configurations
        // ----------------------------------------------------------------------------
 
        if constexpr (traits_t::is_banded)
            throw invalid_alignment_configuration{"Banded alignments are yet not supported."};
 
        // ----------------------------------------------------------------------------
        // Configure semi-global alignment
        // ----------------------------------------------------------------------------
 
        // Get the value for the sequence ends configuration.
        auto align_ends_cfg = cfg.template value_or<align_cfg::aligned_ends>(free_ends_none);
        using align_ends_cfg_t = remove_cvref_t<decltype(align_ends_cfg)>;
 
        auto configure_edit_traits = [&] (auto is_semi_global)
        {
            struct edit_traits_type
            {
                using is_semi_global_type [[maybe_unused]] = remove_cvref_t<decltype(is_semi_global)>;
            };
 
            edit_distance_algorithm<remove_cvref_t<config_t>, edit_traits_type> algorithm{cfg};
            return function_wrapper_t{std::move(algorithm)};
        };
 
        // Check if it has free ends set for the first sequence trailing gaps.
        auto has_free_ends_trailing = [&] (auto first) constexpr
        {
            if constexpr (!decltype(first)::value)
            {
                return configure_edit_traits(std::false_type{});
            }
            else if constexpr (align_ends_cfg_t::template is_static<1>())
            {
#if defined(__GNUC__) && __GNUC__ >= 9
                constexpr bool free_ends_trailing = align_ends_cfg_t::template get_static<1>();
                return configure_edit_traits(std::integral_constant<bool, free_ends_trailing>{});
#else // ^^^ workaround / no workaround vvv
                return configure_edit_traits(std::integral_constant<bool, align_ends_cfg_t::template get_static<1>()>{});
#endif // defined(__GNUC__) && __GNUC__ >= 9
            }
            else // Resolve correct property at runtime.
            {
                if (align_ends_cfg[1])
                    return configure_edit_traits(std::true_type{});
                else
                    return configure_edit_traits(std::false_type{});
            }
        };
 
        // Check if it has free ends set for the first sequence leading gaps.
        // If possible use static information.
        if constexpr (align_ends_cfg_t::template is_static<0>())
            return has_free_ends_trailing(std::integral_constant<bool, align_ends_cfg_t::template get_static<0>()>{});
        else // Resolve correct property at runtime.
        {
            if (align_ends_cfg[0])
                return has_free_ends_trailing(std::true_type{});
            else
                return has_free_ends_trailing(std::false_type{});
        }
    }
 
    template <typename function_wrapper_t, typename ...policies_t, typename config_t>
    static constexpr function_wrapper_t configure_free_ends_initialisation(config_t const & cfg);
 
    template <typename function_wrapper_t, typename ...policies_t, typename config_t>
    static constexpr function_wrapper_t configure_free_ends_optimum_search(config_t const & cfg);
 
    template <typename function_wrapper_t, typename config_t>
    static constexpr function_wrapper_t configure_scoring_scheme(config_t const & cfg);
 
    template <typename function_wrapper_t, typename ...policies_t, typename config_t>
    static constexpr function_wrapper_t make_algorithm(config_t const & cfg)
    {
        using traits_t = alignment_configuration_traits<config_t>;
        using matrix_policy_t = typename select_matrix_policy<traits_t>::type;
        using gap_policy_t = typename select_gap_policy<traits_t>::type;
 
        return alignment_algorithm<config_t, matrix_policy_t, gap_policy_t, policies_t...>{cfg};
    }
};
 
template <typename function_wrapper_t, typename config_t>
constexpr function_wrapper_t alignment_configurator::configure_scoring_scheme(config_t const & cfg)
{
    using traits_t = alignment_configuration_traits<config_t>;
 
    using alignment_scoring_scheme_t =
        lazy_conditional_t<traits_t::is_vectorised,
                           lazy<simd_match_mismatch_scoring_scheme,
                                typename traits_t::score_t,
                                typename traits_t::scoring_scheme_alphabet_t,
                                typename traits_t::alignment_mode_t>,
                           typename traits_t::scoring_scheme_t>;
 
    using scoring_scheme_policy_t = deferred_crtp_base<scoring_scheme_policy, alignment_scoring_scheme_t>;
    return configure_free_ends_initialisation<function_wrapper_t, scoring_scheme_policy_t>(cfg);
}
 
// This function returns a std::function object which can capture runtime dependent alignment algorithm types through
// a fixed invocation interface which is already defined by the caller of this function.
template <typename function_wrapper_t, typename ...policies_t, typename config_t>
constexpr function_wrapper_t alignment_configurator::configure_free_ends_initialisation(config_t const & cfg)
{
    using traits_t = alignment_configuration_traits<config_t>;
    // Get the value for the sequence ends configuration.
    auto align_ends_cfg = cfg.template value_or<align_cfg::aligned_ends>(free_ends_none);
    using align_ends_cfg_t = decltype(align_ends_cfg);
 
    // This lambda augments the initialisation policy of the alignment algorithm
    // with the aligned_ends configuration from before.
    auto configure_leading_both = [&] (auto first_seq, auto second_seq) constexpr
    {
        // Define the trait for the initialisation policy
        struct policy_trait_type
        {
            using free_first_leading_t  [[maybe_unused]] = decltype(first_seq);
            using free_second_leading_t [[maybe_unused]] = decltype(second_seq);
        };
 
        // Make initialisation policy a deferred CRTP base and delegate to configure the find optimum policy.
        using gap_init_policy_t = deferred_crtp_base<affine_gap_init_policy, policy_trait_type>;
        return configure_free_ends_optimum_search<function_wrapper_t, policies_t..., gap_init_policy_t>(cfg);
    };
 
    if constexpr (traits_t::is_local)
    {
        return configure_leading_both(std::true_type{}, std::true_type{});
    }
    else
    {
        // This lambda determines the initialisation configuration for the second sequence given
        // the leading gap property for it.
        auto configure_leading_second = [&] (auto first) constexpr
        {
            // If possible use static information.
            if constexpr (align_ends_cfg_t::template is_static<2>())
            {
                using second_t = std::integral_constant<bool, align_ends_cfg_t::template get_static<2>()>;
                return configure_leading_both(first, second_t{});
            }
            else
            {   // Resolve correct property at runtime.
                if (align_ends_cfg[2])
                    return configure_leading_both(first, std::true_type{});
                else
                    return configure_leading_both(first, std::false_type{});
            }
        };
 
        // Here the initialisation configuration for the first sequence is determined given
        // the leading gap property for it.
        // If possible use static information.
        if constexpr (align_ends_cfg_t::template is_static<0>())
        {
            using first_t = std::integral_constant<bool, align_ends_cfg_t::template get_static<0>()>;
            return configure_leading_second(first_t{});
        }
        else
        {  // Resolve correct property at runtime.
            if (align_ends_cfg[0])
                return configure_leading_second(std::true_type{});
            else
                return configure_leading_second(std::false_type{});
        }
    }
}
 
// This function returns a std::function object which can capture runtime dependent alignment algorithm types through
// a fixed invocation interface which is already defined by the caller of this function.
template <typename function_wrapper_t, typename ...policies_t, typename config_t>
constexpr function_wrapper_t alignment_configurator::configure_free_ends_optimum_search(config_t const & cfg)
{
    using traits_t = alignment_configuration_traits<config_t>;
 
    // Get the value for the sequence ends configuration.
    auto align_ends_cfg = cfg.template value_or<align_cfg::aligned_ends>(free_ends_none);
    using align_ends_cfg_t = decltype(align_ends_cfg);
 
    // This lambda augments the find optimum policy of the alignment algorithm with the
    // respective aligned_ends configuration.
    auto configure_trailing_both = [&] (auto first_seq, auto second_seq) constexpr
    {
        struct policy_trait_type
        {
            using find_in_every_cell_type  [[maybe_unused]] = std::bool_constant<traits_t::is_local>;
            using find_in_last_row_type    [[maybe_unused]] = decltype(first_seq);
            using find_in_last_column_type [[maybe_unused]] = decltype(second_seq);
        };
 
        // We need to select the correct policy based on the configuration traits.
        using find_optimum_t = typename select_find_optimum_policy<traits_t, policy_trait_type>::type;
        return make_algorithm<function_wrapper_t, policies_t..., find_optimum_t>(cfg);
    };
 
    if constexpr (traits_t::is_local)
    {
        return configure_trailing_both(std::true_type{}, std::true_type{});
    }
    else
    {
        // This lambda determines the lookup configuration for the second sequence given
        // the trailing gap property for it.
        auto configure_trailing_second = [&] (auto first) constexpr
        {
            // If possible use static information.
            if constexpr (align_ends_cfg_t::template is_static<3>())
            {
                using second_t = std::integral_constant<bool, align_ends_cfg_t::template get_static<3>()>;
                return configure_trailing_both(first, second_t{});
            }
            else
            { // Resolve correct property at runtime.
                if (align_ends_cfg[3])
                    return configure_trailing_both(first, std::true_type{});
                else
                    return configure_trailing_both(first, std::false_type{});
            }
        };
 
        // Here the lookup configuration for the first sequence is determined given
        // the trailing gap property for it.
        // If possible use static information.
        if constexpr (align_ends_cfg_t::template is_static<1>())
        {
            using first_t = std::integral_constant<bool, align_ends_cfg_t::template get_static<1>()>;
            return configure_trailing_second(first_t{});
        }
        else
        { // Resolve correct property at runtime.
            if (align_ends_cfg[1])
                return configure_trailing_second(std::true_type{});
            else
                return configure_trailing_second(std::false_type{});
        }
    }
}
 
} // namespace seqan3::detail