bitpacked_sequence.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2006-2026 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2026 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <concepts>
#include <iterator>
#include <ranges>
#include <type_traits>
 
#include <seqan3/alphabet/detail/alphabet_proxy.hpp>
#include <seqan3/alphabet/views/to_char.hpp>
#include <seqan3/alphabet/views/to_rank.hpp>
#include <seqan3/contrib/sdsl-lite.hpp>
#include <seqan3/core/concept/cereal.hpp>
#include <seqan3/core/range/detail/random_access_iterator.hpp>
#include <seqan3/utility/math.hpp>
#include <seqan3/utility/views/convert.hpp>
 
namespace seqan3
{
 
template <writable_semialphabet alphabet_type>
    requires std::regular<alphabet_type>
class bitpacked_sequence
{
private:
    static constexpr size_t bits_per_letter = detail::ceil_log2(alphabet_size<alphabet_type>);
 
    static_assert(bits_per_letter <= 64, "alphabet must be representable in at most 64bit.");
 
    using data_type = seqan3::contrib::sdsl::int_vector<bits_per_letter>;
 
    data_type data;
 
    class reference_proxy_type : public alphabet_proxy<reference_proxy_type, alphabet_type>
    {
    private:
        using base_t = alphabet_proxy<reference_proxy_type, alphabet_type>;
        friend base_t;
 
        std::ranges::range_reference_t<data_type> internal_proxy;
 
        constexpr void on_update() noexcept
        {
            // Alternative patch in the SDSL.
            // The problem is that, in the Fedora build type, there is a warning about accessing lo_set[65],
            // which is out of bounds (size is 64).
            // However, during runtime, lo_set[65] is never actually accessed (bogus warning).
            // Making `len > 64` undefined behaviour resolves the warning.
            // Note that the SDSL uses an older C++ standard and has not `std::unreachable()` available.
            // Also, `__builtin_unreachable()` is only valid for GCC and Clang. The SDSL also supports MSVC,
            // where there is an equivalent but different builtin for that.
            // The diff is for the amalgamated seqan3/contrib/sdsl-lite.hpp; but the patch should be applied upstream.
            // ```diff
            // diff --git a/include/seqan3/contrib/sdsl-lite.hpp b/include/seqan3/contrib/sdsl-lite.hpp
            // index a82da4ac2..b59628a56 100644
            // --- a/include/seqan3/contrib/sdsl-lite.hpp
            // +++ b/include/seqan3/contrib/sdsl-lite.hpp
            // @@ -510,6 +510,8 @@ constexpr uint32_t bits_impl<T>::sel11(uint64_t x, uint32_t i, uint32_t c)
            //  template <typename T>
            //  constexpr void bits_impl<T>::write_int(uint64_t * word, uint64_t x, uint8_t offset, const uint8_t len)
            //  {
            // +    if (len > 64)
            // +        __builtin_unreachable();
            //      x &= bits_impl<T>::lo_set[len];
            //      if (offset + len < 64)
            //      {
            // ```
            SEQAN3_WORKAROUND_GCC_BOGUS_MEMCPY_START(-Warray-bounds)
            internal_proxy = base_t::to_rank();
            SEQAN3_WORKAROUND_GCC_BOGUS_MEMCPY_STOP
        }
 
    public:
        // Import from base:
        using base_t::operator=;
 
        reference_proxy_type() = delete;
        constexpr reference_proxy_type(reference_proxy_type const &) noexcept = default;             
        constexpr reference_proxy_type(reference_proxy_type &&) noexcept = default;                  
        constexpr reference_proxy_type & operator=(reference_proxy_type const &) noexcept = default; 
        constexpr reference_proxy_type & operator=(reference_proxy_type &&) noexcept = default;      
        ~reference_proxy_type() noexcept = default;                                                  
 
        reference_proxy_type(std::ranges::range_reference_t<data_type> const internal) noexcept :
            internal_proxy{internal}
        {
            // Call alphabet_base's assign_rank to prevent calling on_update() during construction
            // which is not necessary, because internal_proxy is already correctly initialised!
            base_t::base_t::assign_rank(static_cast<alphabet_rank_t<alphabet_type>>(internal));
        }
    };
 
    static_assert(writable_alphabet<reference_proxy_type>);
    //NOTE(h-2): it is entirely unclear to me why we need this
    template <typename t>
        requires std::is_same_v<std::ranges::range_value_t<std::remove_cvref_t<t>>, alphabet_type>
    static constexpr bool has_same_value_type_v = true;
 
public:
    using value_type = alphabet_type;
    using reference = reference_proxy_type;
    using const_reference = alphabet_type;
    using iterator = detail::random_access_iterator<bitpacked_sequence>;
    using const_iterator = detail::random_access_iterator<bitpacked_sequence const>;
    using difference_type = std::ranges::range_difference_t<data_type>;
    using size_type = std::ranges::range_size_t<data_type>;
 
    bitpacked_sequence() = default;                                                 
    constexpr bitpacked_sequence(bitpacked_sequence const &) = default;             
    constexpr bitpacked_sequence(bitpacked_sequence &&) = default;                  
    constexpr bitpacked_sequence & operator=(bitpacked_sequence const &) = default; 
    constexpr bitpacked_sequence & operator=(bitpacked_sequence &&) = default;      
    ~bitpacked_sequence() = default;                                                
 
    template <typename other_range_t>
        requires (!std::same_as<bitpacked_sequence, std::remove_cvref_t<other_range_t>>)
              && std::ranges::input_range<other_range_t> && has_same_value_type_v<other_range_t>
    explicit bitpacked_sequence(other_range_t && range) :
        bitpacked_sequence{std::ranges::begin(range), std::ranges::end(range)}
    {}
 
    bitpacked_sequence(size_type const count, value_type const value) : data(count, to_rank(value))
    {}
 
    template <std::forward_iterator begin_iterator_type, typename end_iterator_type>
    bitpacked_sequence(begin_iterator_type begin_it, end_iterator_type end_it)
        requires std::sentinel_for<end_iterator_type, begin_iterator_type>
              && std::common_reference_with<std::iter_value_t<begin_iterator_type>, value_type>
    {
        insert(cend(), begin_it, end_it);
    }
 
    bitpacked_sequence(std::initializer_list<value_type> ilist) : bitpacked_sequence(std::begin(ilist), std::end(ilist))
    {}
 
    bitpacked_sequence & operator=(std::initializer_list<value_type> ilist)
    {
        assign(std::begin(ilist), std::end(ilist));
        return *this;
    }
 
    template <std::ranges::input_range other_range_t>
    void assign(other_range_t && range)
        requires std::common_reference_with<std::ranges::range_value_t<other_range_t>, value_type>
    {
        bitpacked_sequence rhs{std::forward<other_range_t>(range)};
        swap(rhs);
    }
 
    void assign(size_type const count, value_type const value)
    {
        bitpacked_sequence rhs{count, value};
        swap(rhs);
    }
 
    template <std::forward_iterator begin_iterator_type, typename end_iterator_type>
    void assign(begin_iterator_type begin_it, end_iterator_type end_it)
        requires std::sentinel_for<end_iterator_type, begin_iterator_type>
              && std::common_reference_with<std::iter_value_t<begin_iterator_type>, value_type>
    {
        bitpacked_sequence rhs{begin_it, end_it};
        swap(rhs);
    }
 
    void assign(std::initializer_list<value_type> ilist)
    {
        assign(std::begin(ilist), std::end(ilist));
    }
 
 
    iterator begin() noexcept
    {
        return iterator{*this};
    }
 
    const_iterator begin() const noexcept
    {
        return const_iterator{*this};
    }
 
    const_iterator cbegin() const noexcept
    {
        return const_iterator{*this};
    }
 
    iterator end() noexcept
    {
        return iterator{*this, size()};
    }
 
    const_iterator end() const noexcept
    {
        return const_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 bitpacked_sequence."};
        }
        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 bitpacked_sequence."};
        }
        return (*this)[i];
    }
 
    reference operator[](size_type const i) noexcept
    {
        assert(i < size());
        return data[i];
    }
 
    const_reference operator[](size_type const i) const noexcept
    {
        assert(i < size());
        return assign_rank_to(data[i], const_reference{});
    }
 
    reference front() noexcept
    {
        assert(size() > 0);
        return (*this)[0];
    }
 
    const_reference front() const noexcept
    {
        assert(size() > 0);
        return (*this)[0];
    }
 
    reference back() noexcept
    {
        assert(size() > 0);
        return (*this)[size() - 1];
    }
 
    const_reference back() const noexcept
    {
        assert(size() > 0);
        return (*this)[size() - 1];
    }
 
    constexpr data_type & raw_data() noexcept
    {
        return data;
    }
 
    constexpr data_type const & raw_data() const noexcept
    {
        return data;
    }
 
    bool empty() const noexcept
    {
        return size() == 0;
    }
 
    size_type size() const noexcept
    {
        return data.size();
    }
 
    size_type max_size() const noexcept
    {
        return data.max_size();
    }
 
    size_type capacity() const noexcept
    {
        return data.capacity();
    }
 
    void reserve(size_type const new_cap)
    {
        data.reserve(new_cap);
    }
 
    void shrink_to_fit()
    {
        data.shrink_to_fit();
    }
 
    void clear() noexcept
    {
        data.clear();
    }
 
    iterator insert(const_iterator pos, value_type const value)
    {
        return insert(pos, 1, value);
    }
 
    iterator insert(const_iterator pos, size_type const count, value_type const value)
    {
        auto const pos_as_num = std::distance(cbegin(), pos); // we want to insert BEFORE this position
 
        data.insert(data.begin() + pos_as_num, count, to_rank(value));
 
        return begin() + pos_as_num;
    }
 
    template <std::forward_iterator begin_iterator_type, typename end_iterator_type>
    iterator insert(const_iterator pos, begin_iterator_type begin_it, end_iterator_type end_it)
        requires std::sentinel_for<end_iterator_type, begin_iterator_type>
              && std::common_reference_with<std::iter_value_t<begin_iterator_type>, value_type>
    {
        auto const pos_as_num = std::distance(cbegin(), pos);
 
        auto v = std::ranges::subrange<begin_iterator_type, end_iterator_type>{begin_it, end_it}
               | views::convert<value_type> | views::to_rank;
        data.insert(data.begin() + pos_as_num, std::ranges::begin(v), std::ranges::end(v));
 
        return begin() + pos_as_num;
    }
 
    iterator insert(const_iterator pos, std::initializer_list<value_type> const & ilist)
    {
        return insert(pos, ilist.begin(), ilist.end());
    }
 
    iterator erase(const_iterator begin_it, const_iterator end_it)
    {
        if (begin_it >= end_it) // [[unlikely]]
            return begin() + std::distance(cbegin(), end_it);
 
        auto const begin_it_pos = std::distance(cbegin(), begin_it);
        auto const end_it_pos = std::distance(cbegin(), end_it);
 
        data.erase(data.cbegin() + begin_it_pos, data.cbegin() + end_it_pos);
 
        return begin() + begin_it_pos;
    }
 
    iterator erase(const_iterator pos)
    {
        return erase(pos, pos + 1);
    }
 
    void push_back(value_type const value)
    {
        data.push_back(to_rank(value));
    }
 
    void pop_back()
    {
        assert(size() > 0);
        data.pop_back();
    }
 
    void resize(size_type const count)
    {
        assert(count < max_size());
        data.resize(count);
    }
 
    void resize(size_type const count, value_type const value)
    {
        assert(count < max_size());
        data.resize(count, to_rank(value));
    }
 
    constexpr void swap(bitpacked_sequence & rhs) noexcept
    {
        std::swap(data, rhs.data);
    }
 
    constexpr void swap(bitpacked_sequence && rhs) noexcept
    {
        std::swap(data, rhs.data);
    }
 
    friend constexpr void swap(bitpacked_sequence & lhs, bitpacked_sequence & rhs) noexcept
    {
        std::swap(lhs, rhs);
    }
 
    friend constexpr void swap(bitpacked_sequence && lhs, bitpacked_sequence && rhs) noexcept
    {
        std::swap(lhs, rhs);
    }
 
    constexpr bool operator==(bitpacked_sequence const & rhs) const noexcept
    {
        return data == rhs.data;
    }
 
    constexpr bool operator!=(bitpacked_sequence const & rhs) const noexcept
    {
        return data != rhs.data;
    }
 
    constexpr bool operator<(bitpacked_sequence const & rhs) const noexcept
    {
        return data < rhs.data;
    }
 
    constexpr bool operator>(bitpacked_sequence const & rhs) const noexcept
    {
        return data > rhs.data;
    }
 
    constexpr bool operator<=(bitpacked_sequence const & rhs) const noexcept
    {
        return data <= rhs.data;
    }
 
    constexpr bool operator>=(bitpacked_sequence const & rhs) const noexcept
    {
        return data >= rhs.data;
    }
 
    template <cereal_archive archive_t>
    void CEREAL_SERIALIZE_FUNCTION_NAME(archive_t & archive)
    {
        archive(data);
    }
};
 
} // namespace seqan3