bit_vector.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
 
// Modified by Enrico Seiler <enrico.seiler AT fu-berlin.de>:
// * To be a single header and non-template (was CRTP)
// * Changed default allocator to 64-byte aligned allocator
// * Changed `bit_vector operator~()` to be auto-vectorizable
// * Changed `binary_transform_impl` to be auto-vectorizable
// * Replaced `seqan3::detail::bits_of<chunk_type>` with `sizeof(chunk_type) * CHAR_BIT`
// * Changed `difference_type` to `size_type` in `operator[]`
// * Changed `1` to `1ULL` in `(1ULL << to_local_chunk_position(size()))` of `resize()`
// * Changed `bit_reference & operator=` to be in accordance with STL
// * Changed `resize` to do nothing when reducing size besides setting new size.
 
#pragma once
 
#include <algorithm>   // for initializer_list, __for_each, for_each, __fill, copy, fill
#include <bit>         // for popcount, countr_zero
#include <cassert>     // for assert
#include <climits>     // for CHAR_BIT
#include <compare>     // for strong_ordering, operator==
#include <cstddef>     // for size_t, ptrdiff_t
#include <cstdint>     // for uint64_t
#include <iterator>    // for iter_reference_t, __next, distance, iter_difference_t, iter_va...
#include <memory>      // for allocator, allocator_traits, assume_aligned, _LIBCPP_COMPRESSE...
#include <ranges>      // for __fn, begin, end
#include <stdexcept>   // for out_of_range
#include <string>      // for char_traits, operator+, operator""s, basic_string, to_string
#include <type_traits> // for conditional_t
#include <utility>     // for swap
#include <vector>      // for vector
 
#include <cereal/cereal.hpp>     // for binary_data
#include <cereal/macros.hpp>     // for CEREAL_LOAD_FUNCTION_NAME, CEREAL_SAVE_FUNCTION_NAME
#include <cereal/specialize.hpp> // for specialization, specialize
 
#include <hibf/cereal/concepts.hpp>           // for cereal_archive, cereal_text_archive
#include <hibf/contrib/aligned_allocator.hpp> // for aligned_allocator
#include <hibf/platform.hpp>                  // for _LIBCPP_VERSION
 
namespace seqan::hibf
{
 
class bit_vector :
    public std::vector<uint64_t,
                       typename std::allocator_traits<
                           seqan::hibf::contrib::aligned_allocator<bool, 64u>>::template rebind_alloc<uint64_t>>
{
private:
    using allocator_t = seqan::hibf::contrib::aligned_allocator<bool, 64u>;
 
    using base_t = std::vector<uint64_t, typename std::allocator_traits<allocator_t>::template rebind_alloc<uint64_t>>;
 
    using chunk_type = uint64_t;
 
    template <bool is_const>
    class bit_reference
    {
    private:
        friend class bit_vector;
 
        using maybe_const_chunk_type = std::conditional_t<is_const, chunk_type const, chunk_type>;
 
        maybe_const_chunk_type * _chunk{}; 
        chunk_type _chunk_mask{};          
 
        constexpr bit_reference(maybe_const_chunk_type * chunk, size_type const local_chunk_position) noexcept :
            _chunk{chunk},
            _chunk_mask{static_cast<chunk_type>(1) << to_local_chunk_position(local_chunk_position)}
        {}
 
    public:
        bit_reference() = delete; 
        bit_reference(bit_reference const & other) = default;
        bit_reference(bit_reference && other) = default;
        bit_reference & operator=(bit_reference const & other)
        {
            *this = bool(other);
            return *this;
        }
        bit_reference & operator=(bit_reference && other) noexcept
        {
            *this = bool(other);
            return *this;
        }
 
        constexpr bit_reference & operator=(bool const bit) noexcept
        {
            bit ? set() : clear();
            return *this;
        }
 
        // Needed to model std::output_iterator<bit_iterator, bool>, which requires the assignment to an const && version
        // of the proxy.
        // NOLINTNEXTLINE(misc-unconventional-assign-operator)
        constexpr bit_reference const & operator=(bool const bit) const noexcept
            requires (!is_const)
        {
            bit ? set() : clear();
            return *this;
        }
 
        constexpr operator bool() const noexcept
        {
            return *_chunk & _chunk_mask;
        }
 
        constexpr bit_reference & flip() noexcept
        {
            (*this) ? clear() : set();
            return *this;
        }
 
    private:
        constexpr void set() noexcept
        {
            *_chunk |= _chunk_mask;
        }
 
        constexpr void clear() noexcept
        {
            *_chunk &= ~_chunk_mask;
        }
    };
 
    template <bool is_const>
    class bit_iterator
    {
    private:
        template <bool>
        friend class bit_iterator;
 
        using maybe_const_chunk_type = std::conditional_t<is_const, chunk_type const, chunk_type>;
 
        maybe_const_chunk_type * _chunk{}; 
        size_type _chunk_position{};       
 
    public:
        using value_type = bool;                                   
        using reference = bit_reference<is_const>;                 
        using pointer = void;                                      
        using difference_type = std::ptrdiff_t;                    
        using iterator_category = std::random_access_iterator_tag; 
        using iterator_concept = std::random_access_iterator_tag;  
 
        bit_iterator() = default; 
 
        explicit constexpr bit_iterator(maybe_const_chunk_type * chunk) noexcept : _chunk{chunk}, _chunk_position{0}
        {}
 
        constexpr bit_iterator(bit_iterator<!is_const> const & other) noexcept
            requires (is_const)
            : _chunk{other._chunk}, _chunk_position{other._chunk_position}
        {}
 
        constexpr reference operator*() const noexcept
        {
            return reference{_chunk, _chunk_position};
        }
 
        constexpr reference operator[](difference_type const count) const noexcept
        {
            return *((*this) + count);
        }
 
        constexpr bit_iterator & operator++() noexcept
        {
            _chunk += !static_cast<bool>(to_local_chunk_position(++_chunk_position));
            return *this;
        }
 
        constexpr bit_iterator operator++(int) noexcept
        {
            bit_iterator tmp{*this};
            ++(*this);
            return tmp;
        }
 
        constexpr bit_iterator & operator+=(difference_type const count) noexcept
        {
            //           chunk:|    0   |    1   |    2   |    3   |    4   |    5   |
            //                 |--------|--------|--------|--------|-x------|--------|
            //  chunk_position:|01234567|01234567|01234567|01234567|01234567|01234567|
            // global position:|01234567|89012345|67890123|45678901|23456789|01234567|
            //                 |0       |  1     |    2   |      3 |        |4       |
            if (count < 0)
            {
                size_type updated_count = modulo_mask - to_local_chunk_position(_chunk_position) - count;
                _chunk_position = modulo_mask - to_local_chunk_position(updated_count);
                _chunk -= to_chunk_position(updated_count);
                //(to_chunk_position(-count) + (old_chunk_position < _chunk_position));
            }
            else
            {
                _chunk += to_chunk_position(to_local_chunk_position(_chunk_position) + count);
                _chunk_position = to_local_chunk_position(_chunk_position + count);
            }
 
            return *this;
        }
 
        constexpr bit_iterator operator+(difference_type const count) const noexcept
        {
            bit_iterator tmp{*this};
            return tmp += count;
        }
 
        friend constexpr bit_iterator operator+(difference_type const count, bit_iterator rhs) noexcept
        {
            return rhs + count;
        }
 
        constexpr bit_iterator & operator--() noexcept
        {
            _chunk -= !static_cast<bool>(to_local_chunk_position(--_chunk_position));
            return *this;
        }
 
        constexpr bit_iterator operator--(int) noexcept
        {
            bit_iterator tmp{*this};
            --(*this);
            return tmp;
        }
 
        constexpr bit_iterator & operator-=(difference_type const count) noexcept
        {
            return *this += -count;
        }
 
        constexpr bit_iterator operator-(difference_type const count) const noexcept
        {
            bit_iterator tmp{*this};
            return tmp -= count;
        }
 
        template <bool is_const_other>
        constexpr difference_type operator-(bit_iterator<is_const_other> rhs) const noexcept
        {
            return ((_chunk - rhs._chunk) << division_mask) -     // number of bits between chunks.
                   to_local_chunk_position(rhs._chunk_position) + // minus the first bits in rhs.
                   to_local_chunk_position(_chunk_position);      // plus the first bits of the lhs
        }
 
        template <bool is_const_other>
        bool operator==(bit_iterator<is_const_other> const & rhs) const
        {
            return _chunk == rhs._chunk
                && (to_local_chunk_position(_chunk_position) == to_local_chunk_position(rhs._chunk_position));
        }
 
        template <bool is_const_other>
        std::strong_ordering operator<=>(bit_iterator<is_const_other> const & rhs) const
        {
            if (std::strong_ordering order = _chunk <=> rhs._chunk; order == std::strong_ordering::equivalent)
                return to_local_chunk_position(_chunk_position) <=> to_local_chunk_position(rhs._chunk_position);
            else
                return order;
        }
    };
 
public:
    using iterator = bit_iterator<false>;
    using const_iterator = bit_iterator<true>;
    using value_type = std::iter_value_t<iterator>;
    using reference = std::iter_reference_t<iterator>;
    using const_reference = std::iter_reference_t<const_iterator>;
    using size_type = size_t;
    using difference_type = std::iter_difference_t<iterator>;
    using allocator_type = allocator_t;
 
private:
    static constexpr size_type chunk_size = sizeof(chunk_type) * CHAR_BIT;
    static constexpr size_type modulo_mask = chunk_size - 1u;
    static constexpr size_type division_mask = std::countr_zero(chunk_size);
 
    size_type _size{};
 
public:
    constexpr bit_vector(allocator_type const & alloc = allocator_type{}) : base_t{alloc}
    {}
    constexpr bit_vector(bit_vector const &) = default;             
    constexpr bit_vector(bit_vector &&) = default;                  
    constexpr bit_vector & operator=(bit_vector const &) = default; 
    constexpr bit_vector & operator=(bit_vector &&) = default;      
    constexpr ~bit_vector() = default;                              
 
    constexpr bit_vector(size_type const count, bool const bit, allocator_type const & alloc = allocator_type{}) :
        base_t{alloc}
    {
        assign(count, bit);
    }
 
    constexpr bit_vector(std::initializer_list<bool> list, allocator_type const & alloc = allocator_type{}) :
        base_t{alloc}
    {
        assign(list);
    }
 
    constexpr bit_vector(size_type const count, allocator_type const & alloc = allocator_type{}) :
        bit_vector{count, bool{}, alloc}
    {}
 
    template <std::input_iterator iterator_t, std::sentinel_for<iterator_t> sentinel_t>
        requires std::assignable_from<bool &, std::iter_reference_t<iterator_t>>
    constexpr void assign(iterator_t first, sentinel_t last)
    {
        bit_vector tmp{}; // To ensure strong exception guarantee.
        if constexpr (std::sized_sentinel_for<sentinel_t, iterator_t>)
            tmp.reserve(std::ranges::distance(first, last));
 
        std::ranges::copy(first, last, std::back_inserter(tmp));
 
        // ----- no exception after this.
        swap(tmp);
        set_new_size(std::ranges::distance(begin(), end()));
    }
 
    constexpr void assign(std::initializer_list<bool> const & ilist)
    {
        assign(std::ranges::begin(ilist), std::ranges::end(ilist));
    }
 
    constexpr void assign(size_type const count, bool const bit)
    {
        resize(count, bit);
        std::ranges::for_each(*as_base(),
                              [value = fill_chunk(bit)](chunk_type & chunk)
                              {
                                  chunk = value;
                              });
    }
 
    constexpr reference operator[](size_type const position) noexcept
    {
        assert(position < size());
 
        return *std::ranges::next(begin(), position);
    }
 
    constexpr const_reference operator[](size_type const position) const noexcept
    {
        assert(position < size());
 
        return *std::ranges::next(begin(), position);
    }
 
    constexpr reference back() noexcept
    {
        assert(!empty()); // Calling on empty container is undefined behaviour.
 
        return (*this)[size() - 1u];
    }
 
    constexpr const_reference back() const noexcept
    {
        assert(!empty()); // Calling on empty container is undefined behaviour.
 
        return (*this)[size() - 1u];
    }
 
    constexpr bool all() const noexcept
    {
        if (empty())
            return true;
 
        constexpr chunk_type mask = ~chunk_type{};
        chunk_type const * const ptr = data();
        size_type const last_chunk = host_size_impl(size()) - 1u;
 
        for (size_t i = 0; i < last_chunk; ++i)
            if (ptr[i] != mask)
                return false;
 
        // The last chunk may contain unset bits beyond `size()`. Those must not be considered.
        return (ptr[last_chunk] | ~last_chunk_mask()) == mask;
    }
 
    constexpr bool any() const noexcept
    {
        if (empty())
            return false;
 
        chunk_type const * const ptr = data();
        size_type const last_chunk = host_size_impl(size()) - 1u;
 
        for (size_t i = 0; i < last_chunk; ++i)
            if (ptr[i])
                return true;
 
        // The last chunk may contain set bits beyond `size()`. Those must not be considered.
        return (ptr[last_chunk] & last_chunk_mask()) != chunk_type{};
    }
 
    constexpr bool none() const noexcept
    {
        return !any();
    }
 
    constexpr size_t count() const noexcept
    {
        if (empty())
            return 0u;
 
        chunk_type const * const ptr = data();
        size_type const size = host_size_impl(this->size());
        size_t result{};
 
        for (size_t i = 0; i < size; ++i)
            result += std::popcount(ptr[i]);
 
        // The last chunk may contain bits beyond `size()`, e.g., after `flip()`. Those must not be counted.
        return result - std::popcount(ptr[size - 1] & ~last_chunk_mask());
    }
 
    constexpr size_type size() const noexcept
    {
        return _size;
    }
 
    constexpr bool empty() const noexcept
    {
        return _size == 0u;
    }
 
    constexpr size_type capacity() const noexcept
    {
        return base_t::capacity() * chunk_size;
    }
 
    constexpr void reserve(size_type const new_capacity)
    {
        base_t::reserve(host_size_impl(new_capacity));
    }
 
    constexpr void push_back(bool bit)
    {
        size_t const new_size = size() + 1u;
        resize(new_size);
        // ---- no exception after this point.
        set_new_size(new_size);
        back() = bit; // set the bit.
    }
 
    constexpr void resize(size_type const count, bool const bit = {})
    {
        base_t::resize(host_size_impl(count));
 
        size_t const old_size = size();
        set_new_size(count);
 
        // If bit is true and we increase the size.
        if (bit && size() > old_size)
            std::ranges::fill(begin() + old_size, end(), bit);
    }
 
    constexpr void clear() noexcept
    {
        base_t::clear();
        _size = 0u;
    }
 
    constexpr bit_vector & operator&=(bit_vector const & rhs) noexcept
    {
        assert(rhs.size() == size());
 
        return binary_transform_impl(rhs,
                                     [](auto const & left_chunk, auto const & right_chunk)
                                     {
                                         return left_chunk & right_chunk;
                                     });
    }
 
    constexpr bit_vector & operator|=(bit_vector const & rhs) noexcept
    {
        assert(rhs.size() == size());
 
        return binary_transform_impl(rhs,
                                     [](auto const & left_chunk, auto const & right_chunk)
                                     {
                                         return left_chunk | right_chunk;
                                     });
    }
 
    constexpr bit_vector & operator^=(bit_vector const & rhs) noexcept
    {
        assert(rhs.size() == size());
 
        return binary_transform_impl(rhs,
                                     [](auto const & left_chunk, auto const & right_chunk)
                                     {
                                         return left_chunk ^ right_chunk;
                                     });
    }
 
    constexpr bit_vector operator~() const noexcept
    {
        bit_vector tmp(size());
 
        tmp.binary_transform_impl(*this,
                                  [](auto const &, auto const & right_chunk)
                                  {
                                      return ~right_chunk;
                                  });
 
        return tmp;
    }
 
    constexpr friend bit_vector operator&(bit_vector lhs, bit_vector const & rhs) noexcept
    {
        return lhs &= rhs;
    }
 
    constexpr friend bit_vector operator|(bit_vector lhs, bit_vector const & rhs) noexcept
    {
        return lhs |= rhs;
    }
 
    constexpr friend bit_vector operator^(bit_vector lhs, bit_vector const & rhs) noexcept
    {
        return lhs ^= rhs;
    }
 
    constexpr bit_vector & and_not(bit_vector const & rhs) noexcept
    {
        assert(rhs.size() == size());
 
        return binary_transform_impl(rhs,
                                     [](auto const & left_chunk, auto const & right_chunk)
                                     {
                                         return left_chunk & ~right_chunk;
                                     });
    }
 
    constexpr bit_vector & flip() noexcept
    {
        std::ranges::for_each(*as_base(),
                              [](chunk_type & chunk)
                              {
                                  chunk = ~chunk;
                              });
        return *this;
    }
 
    constexpr bit_vector & flip(size_type position)
    {
        using namespace std::literals;
 
        if (position >= size())
            throw std::out_of_range{"The given posisiton "s + std::to_string(position) + " is out of the range [0, "s
                                    + std::to_string(size()) + ")!"s};
 
        (*this)[position].flip();
        return *this;
    }
 
    constexpr void swap(bit_vector & other) noexcept
    {
        base_t::swap(*other.as_base());
        std::swap(_size, other._size);
    }
 
    constexpr iterator begin() noexcept
    {
        return iterator{base_t::data()};
    }
 
    constexpr const_iterator begin() const noexcept
    {
        return const_iterator{base_t::data()};
    }
 
    constexpr const_iterator cbegin() const noexcept
    {
        return begin();
    }
 
    constexpr iterator end() noexcept
    {
        return begin() + size();
    }
 
    constexpr const_iterator end() const noexcept
    {
        return begin() + size();
    }
 
    constexpr const_iterator cend() const noexcept
    {
        return end();
    }
 
    [[gnu::always_inline]] inline constexpr chunk_type * data() noexcept
    {
        return std::assume_aligned<allocator_type::alignment>(base_t::data());
    }
 
    [[gnu::always_inline]] inline constexpr chunk_type const * data() const noexcept
    {
        return std::assume_aligned<allocator_type::alignment>(base_t::data());
    }
 
    template <cereal_archive archive_t>
    void CEREAL_LOAD_FUNCTION_NAME(archive_t & archive)
    {
        // Not using `cereal::make_size_tag(_size)`, because the size tag is inferred for text (XML/JSON) archives.
        // For text archives, `cereal::make_size_tag(_size)` would be the number of elements serialised in the for-loop.
        // E.g., `_size == 100` would store `2` (`== host_size_impl(_size)`).
        archive(_size);
        size_t const vector_size = host_size_impl(_size);
 
        resize_for_overwrite(vector_size);
 
        if constexpr (cereal_text_archive<archive_t>)
        {
            for (auto && v : *as_base())
                archive(v);
        }
        else
        {
            archive(cereal::binary_data(data(), vector_size * sizeof(chunk_type)));
        }
    }
 
    template <cereal_archive archive_t>
    void CEREAL_SAVE_FUNCTION_NAME(archive_t & archive) const
    {
        // Not using `cereal::make_size_tag(_size)`, because the size tag is inferred for text (XML/JSON) archives.
        // For text archives, `cereal::make_size_tag(_size)` would be the number of elements serialised in the for-loop.
        // E.g., `_size == 100` would store `2` (`== host_size_impl(_size)`).
        archive(_size);
 
        if constexpr (cereal_text_archive<archive_t>)
        {
            for (auto && v : *as_base())
                archive(v);
        }
        else
        {
            archive(cereal::binary_data(data(), base_t::size() * sizeof(chunk_type)));
        }
    }
 
private:
// If nothing else works: Just use `resize`.
#ifndef HIBF_UNINITIALISED_RESIZE
    inline constexpr void resize_for_overwrite(size_t const size)
    {
        base_t::resize(size);
    }
#else
// libc++: reinterpret_cast to a struct that has the same layout as std::vector.
// All internal members are private.
#    if defined(_LIBCPP_VERSION)
    inline void resize_for_overwrite(size_t const size)
    {
        struct fake_vector
        {
            using allocator_t = typename base_t::allocator_type;
            using pointer = typename std::allocator_traits<allocator_t>::pointer;
 
            pointer begin;
            pointer end;
#        pragma GCC diagnostic push
#        pragma GCC diagnostic ignored "-Wgnu-anonymous-struct" // Warning since clang-21
            _LIBCPP_COMPRESSED_PAIR(pointer, end_cap = nullptr, allocator_type, allocator);
#        pragma GCC diagnostic pop
        };
 
        static_assert(sizeof(fake_vector) == sizeof(base_t));
        static_assert(alignof(fake_vector) == alignof(base_t));
 
        if (size > base_t::capacity())
            base_t::reserve(size);
// Annotate the new memory as contiguous container for llvm's address sanitizer.
#        if __has_feature(address_sanitizer)
        std::__annotate_contiguous_container<allocator_t>(base_t::data(),
                                                          base_t::data() + base_t::capacity(),
                                                          base_t::data() + base_t::size(),
                                                          base_t::data() + size);
#        endif
        fake_vector & vec = reinterpret_cast<fake_vector &>(*this);
        vec.end = vec.begin + size;
    }
// libstdc++: The internal members are protected, so we can access them.
#    else
    inline constexpr void resize_for_overwrite(size_t const size)
    {
        if (size > base_t::capacity())
            base_t::reserve(size);
 
        this->_M_impl._M_finish = this->_M_impl._M_start + size;
    }
#    endif
#endif
 
    template <typename binary_operator_t>
    constexpr bit_vector & binary_transform_impl(bit_vector const & rhs, binary_operator_t && op) noexcept
    {
        chunk_type * const lhs_data = data();
        chunk_type const * const rhs_data = rhs.data();
        size_type const size = host_size_impl(this->size());
 
        for (size_t i = 0; i < size; ++i)
            lhs_data[i] = op(lhs_data[i], rhs_data[i]);
 
        return *this;
    }
 
    constexpr size_type host_size_impl(size_type const count) const noexcept
    {
        return chunks_needed(count);
    }
 
    constexpr void set_new_size(size_type const new_size) noexcept
    {
        _size = new_size;
    }
 
    constexpr base_t const * as_base() const noexcept
    {
        return static_cast<base_t const *>(this);
    }
 
    constexpr base_t * as_base() noexcept
    {
        return static_cast<base_t *>(this);
    }
 
    constexpr size_type chunks_needed(size_type const count) const noexcept
    {
        return (count + 63u) >> 6; // ceil(count/64)
    }
 
    constexpr chunk_type last_chunk_mask() const noexcept
    {
        assert(!empty());
 
        size_type const bits_in_last_chunk = to_local_chunk_position(size());
        return bits_in_last_chunk ? (chunk_type{1u} << bits_in_last_chunk) - 1u : ~chunk_type{};
    }
 
    constexpr chunk_type fill_chunk(bool const bit) const noexcept
    {
        return (bit) ? ~chunk_type{} : chunk_type{};
    }
 
    static constexpr size_type to_local_chunk_position(size_type const position) noexcept
    {
        return position & modulo_mask; // e.g. position % 64
    }
 
    static constexpr size_type to_chunk_position(size_type const position) noexcept
    {
        return position >> division_mask; // e.g. position / 64
    }
};
 
} // namespace seqan::hibf
 
// See https://uscilab.github.io/cereal/serialization_functions.html#inheritance
// seqan::hibf::bit_vector's base class is std::vector
// If we include <cereal/types/vector.hpp> for std::vector serialisation (e.g., HIBF),
// cereal provides these as non-member load/save functions.
// Since both load/save non-member functions (std::vector) and load/save member functions (bit_vector) are available,
// cereal needs to be told which one to use.
namespace cereal
{
 
template <typename archive_t>
struct specialize<archive_t, seqan::hibf::bit_vector, cereal::specialization::member_load_save>
{};
 
} // namespace cereal