gap_decorator.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 <seqan3/std/algorithm>
#include <limits>
#include <seqan3/std/ranges>
#include <set>
#include <tuple>
#include <type_traits>
 
#include <seqan3/alignment/exception.hpp>
#include <seqan3/alphabet/concept.hpp>
#include <seqan3/alphabet/gap/gap.hpp>
#include <seqan3/alphabet/gap/gapped.hpp>
#include <seqan3/utility/views/type_reduce.hpp>
 
namespace seqan3
{
 
template <std::ranges::viewable_range inner_type>
    requires std::ranges::random_access_range<inner_type> && std::ranges::sized_range<inner_type> &&
             (std::is_const_v<std::remove_reference_t<inner_type>> || std::ranges::view<inner_type>)
class gap_decorator
{
private:
    // Declaration of class's iterator types; for the definition see below.
    class gap_decorator_iterator;
 
    using iterator = gap_decorator_iterator;
    using const_iterator = iterator;
 
    using ungapped_view_type = decltype(views::type_reduce(std::declval<inner_type &&>()));
 
public:
    using value_type = gapped<std::ranges::range_value_t<inner_type>>;
 
    using reference = value_type;
 
    using const_reference = reference;
 
    using size_type = std::ranges::range_size_t<inner_type>;
 
    using difference_type = std::ranges::range_difference_t<inner_type>;
 
    using unaligned_sequence_type = inner_type;
 
    gap_decorator() = default;
    gap_decorator(gap_decorator const &) = default; 
    gap_decorator & operator=(gap_decorator const &) = default; 
    gap_decorator(gap_decorator && rhs) = default; 
    gap_decorator & operator=(gap_decorator && rhs) = default; 
    ~gap_decorator() = default; 
 
    template <typename other_range_t>
         requires (!std::same_as<other_range_t, gap_decorator>) &&
                  std::same_as<std::remove_cvref_t<other_range_t>, std::remove_cvref_t<inner_type>> &&
                  std::ranges::viewable_range<other_range_t> // at end, otherwise it competes with the move ctor
    gap_decorator(other_range_t && range) : ungapped_view{views::type_reduce(std::forward<inner_type>(range))}
    {} // TODO (@smehringer) only works for copyable views. Has to be changed once views are not required to be copyable anymore.
    // !\}
 
    size_type size() const
    {
        if (anchors.size())
            return anchors.rbegin()->second + ungapped_view.size();
 
        return ungapped_view.size();
    }
 
    iterator insert_gap(const_iterator const it, size_type const count = 1)
    {
        if (!count) // [[unlikely]]
            return it;
 
        size_type const pos = it - begin();
        assert(pos <= size());
 
        set_iterator_type it_set = anchors.upper_bound(anchor_gap_t{pos, bound_dummy});
 
        if (it_set == anchors.begin()) // will also catch if anchors is empty since begin() == end()
        {
            anchors.emplace_hint(anchors.begin(), anchor_gap_t{pos, count});
        }
        else // there are gaps before pos
        {
            --it_set;
            auto gap_len{it_set->second};
            if (it_set != anchors.begin())
                gap_len -= (*(std::prev(it_set))).second;
 
            if (it_set->first + gap_len >= pos) // extend existing gap
            {
                anchor_gap_t gap{it_set->first, it_set->second + count};
                it_set = anchors.erase(it_set);
                anchors.insert(it_set, gap);
            }
            else                                  // insert new gap
            {
                anchor_gap_t gap{pos, it_set->second + count};
                ++it_set;
                anchors.insert(it_set, gap);
            }
        }
 
        // post-processing: reverse update of succeeding gaps
        rupdate(pos, count);
        return iterator{*this, pos};
    }
 
    iterator erase_gap(const_iterator const it)
    {
        // check if [it, it+gap_len[ covers [first, last[
        if ((*it) != gap{}) // [[unlikely]]
            throw gap_erase_failure("The range to be erased does not correspond to a consecutive gap.");
 
        return erase_gap(it, std::next(it));
    }
 
    iterator erase_gap(const_iterator const first, const_iterator const last)
    {
        size_type const pos1 = first - begin();
        size_type const pos2 = last - begin();
        set_iterator_type it = anchors.upper_bound(anchor_gap_t{pos1, bound_dummy}); // first element greater than pos1
 
        if (it == anchors.begin())
            throw gap_erase_failure{"There is no gap to erase in range [" + std::to_string(pos1) + "," +
                                    std::to_string(pos2) + "]."};
 
        --it;
        size_type const gap_len = gap_length(it);
 
        // check if [it, it+gap_len[ covers [first, last[
        if ((it->first + gap_len) < pos2) // [[unlikely]]
        {
            throw gap_erase_failure{"The range to be erased does not correspond to a consecutive gap."};
        }
        // case 1: complete gap is deleted
        else if (gap_len == pos2 - pos1)
        {
            it = anchors.erase(it);
        }
        // case 2: gap to be deleted in tail or larger than 1 (equiv. to shift tail left, i.e. pos remains unchanged)
        else
        {
            anchor_gap_t gap{it->first, it->second - pos2 + pos1};
            it = anchors.erase(it);
            it = anchors.insert(it, gap); // amortized constant because of hint
            ++it; // update node after the current
        }
 
        // post-processing: forward update of succeeding gaps
        update(it, pos2 - pos1);
 
        return iterator{*this, pos1};
    }
 
    template <typename unaligned_sequence_t> // generic template to use forwarding reference
        requires std::assignable_from<gap_decorator &, unaligned_sequence_t>
    friend void assign_unaligned(gap_decorator & dec, unaligned_sequence_t && unaligned)
    {
        dec = unaligned;
    }
 
    const_iterator begin() const noexcept
    {
        return iterator{*this};
    }
 
    const_iterator cbegin() const noexcept
    {
        return const_iterator{*this};
    }
 
    const_iterator end() const noexcept
    {
        return iterator{*this, size()};
    }
 
    const_iterator cend() const noexcept
    {
        return const_iterator{*this, size()};
    }
 
    reference at(size_type const i)
    {
        if (i >= size()) // [[unlikely]]
            throw std::out_of_range{"Trying to access element behind the last in gap_decorator."};
        return (*this)[i];
    }
 
    const_reference at(size_type const i) const
    {
        if (i >= size()) // [[unlikely]]
            throw std::out_of_range{"Trying to access element behind the last in gap_decorator."};
        return (*this)[i];
    }
 
    reference operator[](size_type const i) const
    {
        return *iterator{*this, i};
    }
 
    friend bool operator==(gap_decorator const & lhs, gap_decorator const & rhs)
    {
        if (lhs.size()  == rhs.size()  &&
            lhs.anchors == rhs.anchors &&
            std::ranges::equal(lhs.ungapped_view, rhs.ungapped_view))
        {
            return true;
        }
 
        return false;
    }
 
    friend bool operator!=(gap_decorator const & lhs, gap_decorator const & rhs)
    {
        return !(lhs == rhs);
    }
 
    friend bool operator<(gap_decorator const & lhs, gap_decorator const & rhs)
    {
        auto lit = lhs.begin();
        auto rit = rhs.begin();
 
        while (lit != lhs.end() && rit != rhs.end() && *lit == *rit)
            ++lit, ++rit;
 
        if (rit == rhs.end())
            return false;           //  lhs == rhs, or rhs prefix of lhs
        else if (lit == lhs.end())
            return true;            // lhs prefix of rhs
 
        return *lit < *rit;
    }
 
    friend bool operator<=(gap_decorator const & lhs, gap_decorator const & rhs)
    {
        auto lit = lhs.begin();
        auto rit = rhs.begin();
 
        while (lit != lhs.end() && rit != rhs.end() && *lit == *rit)
            ++lit, ++rit;
 
        if (lit == lhs.end())
            return true;            // lhs == rhs, or lhs prefix of rhs
        else if (rit == rhs.end())
            return false;           // rhs prefix of lhs
 
        return *lit < *rit;
    }
 
    friend bool operator>(gap_decorator const & lhs, gap_decorator const & rhs)
    {
        return !(lhs <= rhs);
    }
 
    friend bool operator>=(gap_decorator const & lhs, gap_decorator const & rhs)
    {
        return !(lhs < rhs);
    }
 
private:
    using anchor_gap_t = typename std::pair<size_t, size_t>;
 
    using anchor_set_type = std::set<anchor_gap_t>;
 
    using set_iterator_type = typename anchor_set_type::iterator;
 
    constexpr static size_t bound_dummy{std::numeric_limits<size_t>::max()};
 
    size_type gap_length(set_iterator_type it) const
    {
        return (it == anchors.begin()) ? it->second : it->second - (*std::prev(it)).second;
    }
 
    void rupdate(size_type const pos, size_type const offset)
    {
        for (auto it = std::prev(anchors.end(), 1); it->first > pos;)
        {
            anchors.emplace_hint(it, anchor_gap_t{it->first + offset, it->second + offset});
            anchors.erase(*it--);
        }
    }
 
    void update(set_iterator_type it, size_type const offset)
    {
        while (it != anchors.end())
        {
            anchor_gap_t gap{it->first - offset, it->second - offset};
            it = anchors.erase(it);
            it = anchors.insert(it, gap);
            ++it;
        }
    }
 
    ungapped_view_type ungapped_view{};
 
    anchor_set_type anchors{};
};
 
template <std::ranges::viewable_range urng_t>
    requires (!std::ranges::view<std::remove_reference_t<urng_t>>)
gap_decorator(urng_t && range) -> gap_decorator<std::remove_reference_t<urng_t> const &>;
 
template <std::ranges::view urng_t>
gap_decorator(urng_t range) -> gap_decorator<urng_t>;
 
template <std::ranges::viewable_range inner_type>
    requires std::ranges::random_access_range<inner_type> && std::ranges::sized_range<inner_type> &&
             (std::is_const_v<std::remove_reference_t<inner_type>> || std::ranges::view<inner_type>)
class gap_decorator<inner_type>::gap_decorator_iterator
{
protected:
    typename std::add_pointer_t<gap_decorator const> host{nullptr};
    typename gap_decorator::size_type pos{0u};
    int64_t ungapped_view_pos{0}; // must be signed because we need this value to be -1 in case of leading gaps.
    typename gap_decorator::size_type left_gap_end{0};
    typename gap_decorator::set_iterator_type anchor_set_it{};
    bool is_at_gap{true};
 
    void jump(typename gap_decorator::size_type const new_pos)
    {
        assert(new_pos <= host->size());
        pos = new_pos;
 
        anchor_set_it = host->anchors.upper_bound(anchor_gap_t{pos, host->bound_dummy});
        ungapped_view_pos = pos;
 
        if (anchor_set_it != host->anchors.begin())
        {
            typename gap_decorator::set_iterator_type prev{std::prev(anchor_set_it)};
            size_type gap_len{prev->second};
 
            if (prev != host->anchors.begin())
                gap_len -= std::prev(prev)->second;
 
            ungapped_view_pos -= prev->second;
            left_gap_end = prev->first + gap_len;
        }
 
        if (ungapped_view_pos != static_cast<int64_t>(host->ungapped_view.size()) &&
            pos >= left_gap_end && (anchor_set_it == host->anchors.end() || pos < anchor_set_it->first))
            is_at_gap = false;
        else
            is_at_gap = true;
    }
 
public:
    using difference_type = typename gap_decorator::difference_type;
    using value_type = typename gap_decorator::value_type;
    using reference = typename gap_decorator::const_reference;
    using pointer = value_type *;
    using iterator_category = std::bidirectional_iterator_tag;
 
    gap_decorator_iterator() = default; 
    gap_decorator_iterator(gap_decorator_iterator const &) = default; 
    gap_decorator_iterator & operator=(gap_decorator_iterator const &) = default; 
    gap_decorator_iterator(gap_decorator_iterator &&) = default; 
    gap_decorator_iterator & operator=(gap_decorator_iterator &&) = default; 
    ~gap_decorator_iterator() = default; 
 
    explicit gap_decorator_iterator(gap_decorator const & host_) :
        host(&host_), anchor_set_it{host_.anchors.begin()}
    {
        if (host_.anchors.size() && (*host_.anchors.begin()).first == 0) // there are gaps at the very front
        {
            --ungapped_view_pos; // set ungapped_view_pos to -1 so operator++ works without an extra if-branch.
            left_gap_end = anchor_set_it->second;
            ++anchor_set_it;
        }
        else
        {
            is_at_gap = false;
        }
    }
 
    gap_decorator_iterator(gap_decorator const & host_, typename gap_decorator::size_type const pos_) : host(&host_)
    {
        jump(pos_); // random access to pos
    }
 
    gap_decorator_iterator & operator++()
    {
        assert(host); // host is set
        ++pos;
 
        if (pos < left_gap_end) // we stay within the preceding gap stretch
            return *this;
 
        if (anchor_set_it == host->anchors.end() || pos < anchor_set_it->first)
        {   // proceed within the view since we are right of the previous gap but didn't arrive at the right gap yet
            ++ungapped_view_pos;
            if (ungapped_view_pos != static_cast<int64_t>(host->ungapped_view.size()))
                is_at_gap = false;
        }
        else
        {   // we arrived at the right gap and have to update the variables. ungapped_view_pos remains unchanged.
            left_gap_end = anchor_set_it->first + anchor_set_it->second -
                            ((anchor_set_it != host->anchors.begin()) ? (std::prev(anchor_set_it))->second : 0);
            ++anchor_set_it;
            is_at_gap = true;
 
            if (left_gap_end == host->size()) // very last gap
                ++ungapped_view_pos;
        }
 
        return *this;
    }
 
    gap_decorator_iterator operator++(int)
    {
        gap_decorator_iterator cpy{*this};
        ++(*this);
        return cpy;
    }
 
    gap_decorator_iterator & operator+=(difference_type const skip)
    {
        this->jump(this->pos + skip);
        return *this;
    }
 
    gap_decorator_iterator operator+(difference_type const skip) const
    {
        return gap_decorator_iterator{*(this->host), this->pos + skip};
    }
 
    friend gap_decorator_iterator operator+(difference_type const skip, gap_decorator_iterator const & it)
    {
        return it + skip;
    }
 
    gap_decorator_iterator & operator--()
    {
        assert(host); // host is set
        --pos;
 
        if (pos < left_gap_end)
        {   // there was no gap before but we arrive at the left gap and have to update the variables.
            (anchor_set_it != host->anchors.begin()) ? --anchor_set_it : anchor_set_it;
 
            if (anchor_set_it != host->anchors.begin())
            {
                auto prev = std::prev(anchor_set_it);
                left_gap_end = prev->first + prev->second -
                               ((prev != host->anchors.begin()) ? std::prev(prev)->second : 0);
            }
            else // [[unlikely]]
            {
                left_gap_end = 0;
            }
            is_at_gap = true;
        }
        else if (anchor_set_it == host->anchors.end() || pos < anchor_set_it->first)
        {   // we are neither at the left nor right gap
            --ungapped_view_pos;
            is_at_gap = false;
        }
        // else -> no op (we are still within the right gap stretch)
 
        return *this;
    }
 
    gap_decorator_iterator operator--(int)
    {
        gap_decorator_iterator cpy{*this};
        --(*this);
        return cpy;
    }
 
    gap_decorator_iterator & operator-=(difference_type const skip)
    {
        this->jump(this->pos - skip);
        return *this;
    }
 
    gap_decorator_iterator operator-(difference_type const skip) const
    {
        return gap_decorator_iterator{*(this->host), this->pos - skip};
    }
 
    friend gap_decorator_iterator operator-(difference_type const skip, gap_decorator_iterator const & it)
    {
        return it - skip;
    }
 
    difference_type operator-(gap_decorator_iterator const lhs) const noexcept
    {
        return static_cast<difference_type>(this->pos - lhs.pos);
    }
 
    reference operator*() const
    {
        return (is_at_gap) ? reference{gap{}} : reference{host->ungapped_view[ungapped_view_pos]};
    }
 
    reference operator[](difference_type const n) const
    {
        return *(*this + n);
    }
 
    friend bool operator==(gap_decorator_iterator const & lhs, gap_decorator_iterator const & rhs) noexcept
    {
        return lhs.pos == rhs.pos;
    }
 
    friend bool operator!=(gap_decorator_iterator const & lhs, gap_decorator_iterator const & rhs) noexcept
    {
        return lhs.pos != rhs.pos;
    }
 
    friend bool operator<(gap_decorator_iterator const & lhs, gap_decorator_iterator const & rhs) noexcept
    {
        return lhs.pos < rhs.pos;
    }
 
    friend bool operator>(gap_decorator_iterator const & lhs, gap_decorator_iterator const & rhs) noexcept
    {
        return lhs.pos > rhs.pos;
    }
 
    friend bool operator<=(gap_decorator_iterator const & lhs, gap_decorator_iterator const & rhs) noexcept
    {
        return lhs.pos <= rhs.pos;
    }
 
    friend bool operator>=(gap_decorator_iterator const & lhs, gap_decorator_iterator const & rhs) noexcept
    {
        return lhs.pos >= rhs.pos;
    }
};
 
} // namespace seqan