unidirectional_search_algorithm.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 <ranges>
#include <type_traits>
 
#include <seqan3/alphabet/concept.hpp>
#include <seqan3/core/detail/test_accessor.hpp>
#include <seqan3/search/detail/search_common.hpp>
#include <seqan3/search/detail/search_traits.hpp>
#include <seqan3/search/fm_index/concept.hpp>
 
namespace seqan3::detail
{
 
enum class error_type : uint8_t
{
    deletion,  
    insertion, 
    matchmm,   
    none       
};
 
template <typename configuration_t, typename index_t, typename... policies_t>
class unidirectional_search_algorithm : protected policies_t...
{
private:
    using traits_t = search_traits<configuration_t>;
    using search_result_type = typename traits_t::search_result_type;
 
    static_assert(!std::same_as<search_result_type, empty_type>, "The search result type was not configured.");
 
public:
    unidirectional_search_algorithm() = default;                                                    
    unidirectional_search_algorithm(unidirectional_search_algorithm const &) = default;             
    unidirectional_search_algorithm(unidirectional_search_algorithm &&) = default;                  
    unidirectional_search_algorithm & operator=(unidirectional_search_algorithm const &) = default; 
    unidirectional_search_algorithm & operator=(unidirectional_search_algorithm &&) = default;      
    ~unidirectional_search_algorithm() = default;                                                   
 
    unidirectional_search_algorithm(configuration_t const & cfg, index_t const & index) : policies_t{cfg}...
    {
        stratum = cfg.get_or(search_cfg::hit_strata{0}).stratum;
        index_ptr = &index;
    }
 
    template <typename indexed_query_t, typename callback_t>
        requires (std::tuple_size_v<indexed_query_t> == 2)
              && std::ranges::forward_range<std::tuple_element_t<1, indexed_query_t>>
              && std::invocable<callback_t, search_result_type>
    void operator()(indexed_query_t && indexed_query, callback_t && callback)
    {
        auto && [query_idx, query] = indexed_query;
        auto error_state = this->max_error_counts(query); // see policy_max_error
 
        // construct internal delegate for collecting hits for later filtering (if necessary)
        std::vector<typename index_t::cursor_type> internal_hits{};
        delegate = [&internal_hits](auto const & it)
        {
            internal_hits.push_back(it);
        };
 
        perform_search_by_hit_strategy(internal_hits, query, error_state);
 
        this->make_results(std::move(internal_hits), query_idx, callback); // see policy_search_result_builder
    }
 
private:
    index_t const * index_ptr{nullptr};
 
    std::function<void(typename index_t::cursor_type const &)> delegate;
 
    uint8_t stratum{};
 
    // forward declaration
    template <bool abort_on_hit, typename query_t>
    bool search_trivial(typename index_t::cursor_type cur,
                        query_t & query,
                        typename index_t::cursor_type::size_type const query_pos,
                        search_param const error_left,
                        error_type const prev_error);
 
    template <typename query_t>
    void perform_search_by_hit_strategy(std::vector<typename index_t::cursor_type> & internal_hits,
                                        query_t & query,
                                        search_param error_state)
    {
        if constexpr (!traits_t::search_all_hits)
        {
            auto max_total = error_state.total;
            error_state.total = 0; // start search with less errors
 
            while (internal_hits.empty() && error_state.total <= max_total)
            {
                // * If you only want the best hit (traits_t::search_single_best_hit), you stop after finding the
                //   first hit, the hit with the least errors (`abort_on_hit` is true).
                // * If you are in strata mode (traits_t::search_strata_hits), you do the same as with best hits,
                //   but then do the extra step afterwards (`abort_on_hit` is true).
                // * If you want all best hits (traits_t::search_all_best_hits), you do not stop after the first
                //   hit but continue the current search algorithm/max_error pattern (`abort_on_hit` is true).
                constexpr bool abort_on_hit = !traits_t::search_all_best_hits;
                search_trivial<abort_on_hit>(index_ptr->cursor(), query, 0, error_state, error_type::none);
                ++error_state.total;
            }
 
            if constexpr (traits_t::search_strata_hits)
            {
                if (!internal_hits.empty())
                {
                    internal_hits.clear();
                    error_state.total += stratum - 1;
                    search_trivial<false>(index_ptr->cursor(), query, 0, error_state, error_type::none);
                }
            }
        }
        else // traits_t::search_all
        {
            // If you want to find all hits, you cannot stop once you found any hit (<false>)
            // since you have to find all paths in the search tree that satisfy the hit condition.
            search_trivial<false>(index_ptr->cursor(), query, 0, error_state, error_type::none);
        }
    }
};
 
template <typename configuration_t, typename index_t, typename... policies_t>
template <bool abort_on_hit, typename query_t>
inline bool unidirectional_search_algorithm<configuration_t, index_t, policies_t...>::search_trivial(
    typename index_t::cursor_type cur,
    query_t & query,
    typename index_t::cursor_type::size_type const query_pos,
    search_param const error_left,
    error_type const prev_error)
{
    // Exact case (end of query sequence or no errors left)
    if (query_pos == std::ranges::size(query) || error_left.total == 0)
    {
        // If not at end of query sequence, try searching the remaining suffix without any errors.
        using drop_size_t = std::ranges::range_difference_t<query_t>;
        if (query_pos == std::ranges::size(query)
            || cur.extend_right(std::views::drop(query, static_cast<drop_size_t>(query_pos))))
        {
            delegate(cur);
            return true;
        }
    }
    // Approximate case
    else
    {
        // Insertion
        // Only allow insertions if there is no match and we are not at the beginning of the query.
        bool const allow_insertion =
            (cur.query_length() > 0) ? cur.last_rank() != seqan3::to_rank(query[query_pos]) : true;
 
        if (allow_insertion && (prev_error != error_type::deletion || error_left.substitution == 0)
            && error_left.insertion > 0)
        {
            search_param error_left2{error_left};
            error_left2.insertion--;
            error_left2.total--;
 
            // Always perform a recursive call. Abort recursion if and only if recursive call found a hit and
            // abort_on_hit is set to true.
            if (search_trivial<abort_on_hit>(cur, query, query_pos + 1, error_left2, error_type::insertion)
                && abort_on_hit)
            {
                return true;
            }
        }
 
        // Do not allow deletions at the beginning of the query sequence
        if (((query_pos > 0 && error_left.deletion > 0) || error_left.substitution > 0) && cur.extend_right())
        {
            do
            {
                // Match (when error_left.substitution > 0) and Mismatch
                if (error_left.substitution > 0)
                {
                    bool delta = cur.last_rank() != seqan3::to_rank(query[query_pos]);
                    search_param error_left2{error_left};
                    error_left2.total -= delta;
                    error_left2.substitution -= delta;
 
                    if (search_trivial<abort_on_hit>(cur, query, query_pos + 1, error_left2, error_type::matchmm)
                        && abort_on_hit)
                    {
                        return true;
                    }
                }
 
                // Deletion (Do not allow deletions at the beginning of the query sequence.)
                if (query_pos > 0)
                {
                    // Match (when error_left.substitution == 0)
                    if (error_left.substitution == 0 && cur.last_rank() == seqan3::to_rank(query[query_pos]))
                    {
                        if (search_trivial<abort_on_hit>(cur, query, query_pos + 1, error_left, error_type::matchmm)
                            && abort_on_hit)
                        {
                            return true;
                        }
                    }
 
                    // Deletions at the end of the sequence are not allowed. When the algorithm
                    // arrives here, it cannot be at the end of the query and since deletions do not touch the query
                    // (i.e. increase query_pos) it won't be at the end of the query after the deletion.
                    // Do not allow deletions after an insertion.
                    if ((prev_error != error_type::insertion || error_left.substitution == 0)
                        && error_left.deletion > 0)
                    {
                        search_param error_left2{error_left};
                        error_left2.total--;
                        error_left2.deletion--;
                        // Only search for characters different from the corresponding query character.
                        // (Same character is covered by a match.)
                        if (cur.last_rank() != seqan3::to_rank(query[query_pos]))
                        {
                            if (search_trivial<abort_on_hit>(cur, query, query_pos, error_left2, error_type::deletion)
                                && abort_on_hit)
                            {
                                return true;
                            }
                        }
                    }
                }
            }
            while (cur.cycle_back());
        }
        else
        {
            // Match (when error_left.substitution == 0)
            if (cur.extend_right(query[query_pos]))
            {
                if (search_trivial<abort_on_hit>(cur, query, query_pos + 1, error_left, error_type::matchmm)
                    && abort_on_hit)
                {
                    return true;
                }
            }
        }
    }
 
    return false;
}
 
} // namespace seqan3::detail