gap_decorator.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2019, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2019, 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 <limits>
#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/range/container/concept.hpp>
#include <seqan3/range/detail/random_access_iterator.hpp>
#include <seqan3/range/view/view_all.hpp>
#include <seqan3/std/algorithm>
#include <seqan3/std/ranges>

namespace seqan3
{

template <std::ranges::ViewableRange inner_type>
    requires std::ranges::RandomAccessRange<inner_type> && std::ranges::SizedRange<inner_type> &&
             (std::is_const_v<std::remove_reference_t<inner_type>> || std::ranges::View<inner_type>)
class gap_decorator
{
private:
    class gap_decorator_iterator
    {
    private:
        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; // = reference
        using const_reference = reference;
        using pointer = value_type *;
        using iterator_category = std::bidirectional_iterator_tag;

        constexpr gap_decorator_iterator() = default;
        constexpr gap_decorator_iterator(gap_decorator_iterator const &) = default;
        constexpr gap_decorator_iterator & operator=(gap_decorator_iterator const &) = default;
        constexpr gap_decorator_iterator (gap_decorator_iterator &&) = default;
        constexpr gap_decorator_iterator & operator=(gap_decorator_iterator &&) = default;
        ~gap_decorator_iterator() = default;

        explicit constexpr gap_decorator_iterator(gap_decorator const & host_) noexcept :
            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;
            }
        }

        constexpr gap_decorator_iterator(gap_decorator const & host_,
                                                    typename gap_decorator::size_type const pos_) noexcept :
             host(&host_)
        {
            jump(pos_); // random access to pos
        }

        constexpr gap_decorator_iterator & operator++() noexcept
        {
            assert(host); // host is set
            ++pos;

            if (pos < left_gap_end)
            {   // we stay within the preceding gap stretch
                // no op but must precede other cases
            }
            else 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;
            }

            return *this;
        }

        constexpr gap_decorator_iterator & operator--() noexcept
        {
            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;
        }

        constexpr gap_decorator_iterator operator++(int) noexcept
        {
            gap_decorator_iterator cpy{*this};
            ++(*this);
            return cpy;
        }

        constexpr gap_decorator_iterator operator--(int) noexcept
        {
            gap_decorator_iterator cpy{*this};
            --(*this);
            return cpy;
        }

        constexpr reference operator*() const noexcept
        {
            return (is_at_gap) ? static_cast<reference>(gap{})
                               : static_cast<reference>(host->ungapped_view[ungapped_view_pos]);
        }

        constexpr friend bool operator==(gap_decorator_iterator const & lhs,
                                         gap_decorator_iterator const & rhs)
        {
            return lhs.pos == rhs.pos;
        }

        constexpr friend bool operator!=(gap_decorator_iterator const & lhs,
                                         gap_decorator_iterator const & rhs)
        {
            return lhs.pos != rhs.pos;
        }

        constexpr friend bool operator<(gap_decorator_iterator const & lhs,
                                        gap_decorator_iterator const & rhs)
        {
            return lhs.pos < rhs.pos;
        }

        constexpr friend bool operator>(gap_decorator_iterator const & lhs,
                                        gap_decorator_iterator const & rhs)
        {
            return lhs.pos > rhs.pos;
        }

        constexpr friend bool operator<=(gap_decorator_iterator const & lhs,
                                         gap_decorator_iterator const & rhs)
        {
            return lhs.pos <= rhs.pos;
        }

        constexpr friend bool operator>=(gap_decorator_iterator const & lhs,
                                         gap_decorator_iterator const & rhs)
        {
            return lhs.pos >= rhs.pos;
        }
    };

public:
    using value_type = gapped<value_type_t<inner_type>>;
    using reference = value_type;
    using const_reference = reference;
    using size_type = size_type_t<inner_type>;
    using difference_type = difference_type_t<inner_type>;
    using iterator = gap_decorator_iterator;
    using const_iterator = iterator;

    using unaligned_seq_type = inner_type;

    constexpr gap_decorator() = default;
    constexpr gap_decorator(gap_decorator const &) = default;
    constexpr gap_decorator & operator=(gap_decorator const &) = default;
    constexpr gap_decorator(gap_decorator && rhs) = default;
    constexpr gap_decorator & operator=(gap_decorator && rhs) = default;
    ~gap_decorator() = default;

    template <typename other_range_t>
         requires !std::Same<other_range_t, gap_decorator> &&
                  std::Same<remove_cvref_t<other_range_t>, remove_cvref_t<inner_type>> &&
                  std::ranges::ViewableRange<other_range_t> // at end, otherwise it competes with the move ctor
    gap_decorator(other_range_t && range) : ungapped_view{view::all(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 noexcept
    {
        if (anchors.size())
            return anchors.rbegin()->second + ungapped_view.size();

        return ungapped_view.size();
    }

    iterator insert_gap(iterator const it, size_type const count = 1)
    {
        assert(ungapped_view.size());

        if (!count) // [[unlikely]]
            return it;

        size_type const pos = std::distance(begin(), it);
        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(iterator const it)
    {
        assert(ungapped_view.size());

        // 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.");

        auto end_it = std::next(it);
        return erase_gap(it, end_it);
    }

    iterator erase_gap(iterator const first, iterator const last)
    {
        assert(ungapped_view.size());

        size_type const pos1 = std::distance(begin(), first);
        size_type const pos2 = std::distance(begin(), last);
        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_seq_t> // generic template to use forwarding reference
        requires std::Constructible<gap_decorator, unaligned_seq_t>
    friend void assign_unaligned(gap_decorator & dec, unaligned_seq_t && unaligned)
    {
        dec = unaligned;
    }

    iterator begin() const noexcept
    {
        return iterator{*this};
    }

    const_iterator cbegin() const noexcept
    {
        return const_iterator{*this};
    }

    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];
    }

    constexpr reference operator[](size_type const i) const noexcept
    {
        return *iterator{*this, i};
    }

    friend bool operator==(gap_decorator const & lhs, gap_decorator const & rhs) noexcept
    {
        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) noexcept
    {
        return !(lhs == rhs);
    }

    friend bool operator<(gap_decorator const & lhs, gap_decorator const & rhs) noexcept
    {
        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) noexcept
    {
        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) noexcept
    {
        return !(lhs <= rhs);
    }

    friend bool operator>=(gap_decorator const & lhs, gap_decorator const & rhs) noexcept
    {
        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()};

    constexpr size_type gap_length(set_iterator_type it) const noexcept
    {
        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;
        }
    }

    decltype(view::all(std::declval<inner_type &&>())) ungapped_view{};

    anchor_set_type anchors{};
};

template <std::ranges::ViewableRange 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>;

} // namespace seqan

namespace seqan3::detail
{

template <typename type>
constexpr int enable_view<seqan3::gap_decorator<type>> = 0;

template <typename type>
constexpr int enable_view<seqan3::gap_decorator<type> const> = 0;

} // namespace seqan3::detail