pairwise_alignment_algorithm_banded.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 <seqan3/std/concepts>
#include <seqan3/std/ranges>
 
#include <seqan3/alignment/pairwise/detail/pairwise_alignment_algorithm.hpp>
#include <seqan3/range/views/drop.hpp>
#include <seqan3/range/views/take.hpp>
#include <seqan3/range/views/slice.hpp>
 
namespace seqan3::detail
{
 
template <typename alignment_configuration_t, typename ...policies_t>
    requires is_type_specialisation_of_v<alignment_configuration_t, configuration>
class pairwise_alignment_algorithm_banded :
    protected pairwise_alignment_algorithm<alignment_configuration_t, policies_t...>
{
protected:
    using base_algorithm_t = pairwise_alignment_algorithm<alignment_configuration_t, policies_t...>;
 
    // import types from base class.
    using typename base_algorithm_t::traits_type;
    using typename base_algorithm_t::alignment_result_type;
    using typename base_algorithm_t::score_type;
 
    static_assert(!std::same_as<alignment_result_type, empty_type>, "Alignment result type was not configured.");
    static_assert(traits_type::is_banded, "Alignment configuration must have band configured.");
 
public:
    pairwise_alignment_algorithm_banded() = default; 
    pairwise_alignment_algorithm_banded(pairwise_alignment_algorithm_banded const &) = default; 
    pairwise_alignment_algorithm_banded(pairwise_alignment_algorithm_banded &&) = default; 
    pairwise_alignment_algorithm_banded & operator=(pairwise_alignment_algorithm_banded const &) = default; 
    pairwise_alignment_algorithm_banded & operator=(pairwise_alignment_algorithm_banded &&) = default; 
    ~pairwise_alignment_algorithm_banded() = default; 
 
    pairwise_alignment_algorithm_banded(alignment_configuration_t const & config) : base_algorithm_t(config)
    {}
 
    template <indexed_sequence_pair_range indexed_sequence_pairs_t, typename callback_t>
        requires std::invocable<callback_t, alignment_result_type>
    void operator()(indexed_sequence_pairs_t && indexed_sequence_pairs, callback_t && callback)
    {
        using std::get;
 
        for (auto && [sequence_pair, idx] : indexed_sequence_pairs)
        {
            size_t const sequence1_size = std::ranges::distance(get<0>(sequence_pair));
            size_t const sequence2_size = std::ranges::distance(get<1>(sequence_pair));
 
            auto && [alignment_matrix, index_matrix] = this->acquire_matrices(sequence1_size,
                                                                              sequence2_size,
                                                                              this->lowest_viable_score());
 
            compute_matrix(get<0>(sequence_pair), get<1>(sequence_pair), alignment_matrix, index_matrix);
            this->make_result_and_invoke(std::forward<decltype(sequence_pair)>(sequence_pair),
                                         std::move(idx),
                                         this->optimal_score,
                                         this->optimal_coordinate,
                                         callback);
        }
    }
 
protected:
    template <std::ranges::forward_range sequence1_t,
              std::ranges::forward_range sequence2_t,
              std::ranges::input_range alignment_matrix_t,
              std::ranges::input_range index_matrix_t>
        requires std::ranges::forward_range<std::ranges::range_reference_t<alignment_matrix_t>> &&
                 std::ranges::forward_range<std::ranges::range_reference_t<index_matrix_t>>
    void compute_matrix(sequence1_t && sequence1,
                        sequence2_t && sequence2,
                        alignment_matrix_t && alignment_matrix,
                        index_matrix_t && index_matrix)
    {
        // ---------------------------------------------------------------------
        // Initialisation phase: allocate memory and initialise first column.
        // ---------------------------------------------------------------------
 
        this->reset_optimum(); // Reset the tracker for the new alignment computation.
 
        auto alignment_matrix_it = alignment_matrix.begin();
        auto indexed_matrix_it = index_matrix.begin();
 
        size_t row_size = std::max<int32_t>(0, -this->lower_diagonal);
        size_t const column_size = std::max<int32_t>(0, this->upper_diagonal);
        this->initialise_column(*alignment_matrix_it, *indexed_matrix_it, sequence2 | views::take(row_size));
 
        // ---------------------------------------------------------------------
        // 1st recursion phase: band intersects with the first row.
        // ---------------------------------------------------------------------
 
        for (auto sequence1_value : sequence1 | views::take(column_size))
        {
            this->compute_column(*++alignment_matrix_it,
                                 *++indexed_matrix_it,
                                 sequence1_value,
                                 sequence2 | views::take(++row_size));
        }
 
        // ---------------------------------------------------------------------
        // 2nd recursion phase: iterate until the end of the matrix.
        // ---------------------------------------------------------------------
 
        size_t first_row_index = 0u;
        for (auto sequence1_value : sequence1 | views::drop(column_size))
        {
            compute_band_column(*++alignment_matrix_it,
                                *++indexed_matrix_it | views::drop(first_row_index + 1),
                                sequence1_value,
                                sequence2 | views::slice(first_row_index, ++row_size));
            ++first_row_index;
        }
 
        // ---------------------------------------------------------------------
        // Final phase: track score of last column
        // ---------------------------------------------------------------------
 
        auto alignment_column = *alignment_matrix_it;
        auto cell_index_column = *indexed_matrix_it | views::drop(first_row_index);
 
        auto alignment_column_it = alignment_column.begin();
        auto cell_index_column_it = cell_index_column.begin();
 
        this->track_last_column_cell(*alignment_column_it, *cell_index_column_it);
 
        for (size_t last_row = std::min<size_t>(std::ranges::distance(sequence2), row_size);
             first_row_index < last_row;
             ++first_row_index)
            this->track_last_column_cell(*++alignment_column_it, *++cell_index_column_it);
 
        this->track_final_cell(*alignment_column_it, *cell_index_column_it);
    }
 
    template <std::ranges::input_range alignment_column_t,
              std::ranges::input_range cell_index_column_t,
              semialphabet sequence1_value_t,
              std::ranges::input_range sequence2_t>
    void compute_band_column(alignment_column_t && alignment_column,
                             cell_index_column_t && cell_index_column,
                             sequence1_value_t const & sequence1_value,
                             sequence2_t && sequence2)
    {
        // ---------------------------------------------------------------------
        // Initial phase: prepare column and initialise first cell
        // ---------------------------------------------------------------------
 
        auto alignment_column_it = alignment_column.begin();
        auto cell_index_column_it = cell_index_column.begin();
 
        auto cell = *alignment_column_it;
        cell = this->track_cell(
                this->initialise_band_first_cell(cell.best_score(),
                                                 *++alignment_column_it,
                                                 this->scoring_scheme.score(sequence1_value,
                                                                            *std::ranges::begin(sequence2))),
                *cell_index_column_it);
 
        // ---------------------------------------------------------------------
        // Iteration phase: iterate over column and compute each cell
        // ---------------------------------------------------------------------
 
        for (auto && sequence2_value : sequence2 | views::drop(1))
        {
            auto cell = *alignment_column_it;
            cell = this->track_cell(
                this->compute_inner_cell(cell.best_score(),
                                         *++alignment_column_it,
                                         this->scoring_scheme.score(sequence1_value, sequence2_value)),
                *++cell_index_column_it);
        }
 
        // ---------------------------------------------------------------------
        // Final phase: track last cell
        // ---------------------------------------------------------------------
 
        this->track_last_row_cell(*alignment_column_it, *cell_index_column_it);
    }
};
 
} // namespace seqan3::detail