interleaved_bloom_filter.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 <seqan3/std/bit>
 
#include <sdsl/bit_vectors.hpp>
 
#include <seqan3/core/concept/cereal.hpp>
#include <seqan3/core/detail/strong_type.hpp>
 
namespace seqan3
{
 
enum data_layout : bool
{
    uncompressed, 
    compressed    
};
 
struct bin_count : public detail::strong_type<size_t, bin_count, detail::strong_type_skill::convert>
{
    using detail::strong_type<size_t, bin_count, detail::strong_type_skill::convert>::strong_type;
};
 
struct bin_size : public detail::strong_type<size_t, bin_size, detail::strong_type_skill::convert>
{
    using detail::strong_type<size_t, bin_size, detail::strong_type_skill::convert>::strong_type;
};
 
struct hash_function_count : public detail::strong_type<size_t, hash_function_count, detail::strong_type_skill::convert>
{
    using detail::strong_type<size_t, hash_function_count, detail::strong_type_skill::convert>::strong_type;
};
 
struct bin_index : public detail::strong_type<size_t, bin_index, detail::strong_type_skill::convert>
{
    using detail::strong_type<size_t, bin_index, detail::strong_type_skill::convert>::strong_type;
};
 
template <data_layout data_layout_mode_ = data_layout::uncompressed>
class interleaved_bloom_filter
{
private:
    template <data_layout data_layout_mode>
    friend class interleaved_bloom_filter;
 
    using data_type = std::conditional_t<data_layout_mode_ == data_layout::uncompressed,
                                         sdsl::bit_vector,
                                         sdsl::sd_vector<>>;
 
    size_t bins{};
    size_t technical_bins{};
    size_t bin_size_{};
    size_t hash_shift{};
    size_t bin_words{};
    size_t hash_funs{};
    data_type data{};
    static constexpr std::array<size_t, 5> hash_seeds{13572355802537770549ULL, // 2**64 / (e/2)
                                                      13043817825332782213ULL, // 2**64 / sqrt(2)
                                                      10650232656628343401ULL, // 2**64 / sqrt(3)
                                                      16499269484942379435ULL, // 2**64 / (sqrt(5)/2)
                                                      4893150838803335377ULL}; // 2**64 / (3*pi/5)
 
    inline constexpr size_t hash_and_fit(size_t h, size_t const seed) const
    {
        h *= seed;
        assert(hash_shift < 64);
        h ^= h >> hash_shift; // XOR and shift higher bits into lower bits
        h *= 11400714819323198485ULL; // = 2^64 / golden_ration, to expand h to 64 bit range
        // Use fastrange (integer modulo without division) if possible.
#ifdef __SIZEOF_INT128__
        h = static_cast<uint64_t>((static_cast<__uint128_t>(h) * static_cast<__uint128_t>(bin_size_)) >> 64);
#else
        h %= bin_size_;
#endif
        h *= technical_bins;
        return h;
    }
 
public:
    static constexpr data_layout data_layout_mode = data_layout_mode_;
 
    class membership_agent; // documented upon definition below
    template <std::integral value_t>
    class counting_agent_type; // documented upon definition below
 
    interleaved_bloom_filter() = default; 
    interleaved_bloom_filter(interleaved_bloom_filter const &) = default; 
    interleaved_bloom_filter & operator=(interleaved_bloom_filter const &) = default; 
    interleaved_bloom_filter(interleaved_bloom_filter &&) = default; 
    interleaved_bloom_filter & operator=(interleaved_bloom_filter &&) = default; 
    ~interleaved_bloom_filter() = default; 
 
    interleaved_bloom_filter(seqan3::bin_count bins_,
                             seqan3::bin_size size,
                             seqan3::hash_function_count funs = seqan3::hash_function_count{2u})
        requires (data_layout_mode == data_layout::uncompressed)
    {
        bins = bins_.get();
        bin_size_ = size.get();
        hash_funs = funs.get();
 
        if (bins == 0)
            throw std::logic_error{"The number of bins must be > 0."};
        if (hash_funs == 0 || hash_funs > 5)
            throw std::logic_error{"The number of hash functions must be > 0 and <= 5."};
        if (bin_size_ == 0)
            throw std::logic_error{"The size of a bin must be > 0."};
 
        hash_shift = std::countl_zero(bin_size_);
        bin_words = (bins + 63) >> 6; // = ceil(bins/64)
        technical_bins  = bin_words << 6; // = bin_words * 64
        data = sdsl::bit_vector(technical_bins * bin_size_);
    }
 
    interleaved_bloom_filter(interleaved_bloom_filter<data_layout::uncompressed> const & ibf)
        requires (data_layout_mode == data_layout::compressed)
    {
        std::tie(bins, technical_bins, bin_size_, hash_shift, bin_words, hash_funs) =
            std::tie(ibf.bins, ibf.technical_bins, ibf.bin_size_, ibf.hash_shift, ibf.bin_words, ibf.hash_funs);
 
        data = sdsl::sd_vector<>{ibf.data};
    }
 
    void emplace(size_t const value, bin_index const bin) noexcept
        requires (data_layout_mode == data_layout::uncompressed)
    {
        assert(bin.get() < bins);
        for (size_t i = 0; i < hash_funs; ++i)
        {
            size_t idx = hash_and_fit(value, hash_seeds[i]);
            idx += bin.get();
            assert(idx < data.size());
            data[idx] = 1;
        };
    }
 
    void clear(bin_index const bin) noexcept
        requires (data_layout_mode == data_layout::uncompressed)
    {
        assert(bin.get() < bins);
        for (size_t idx = bin.get(), i = 0; i < bin_size_; idx += technical_bins, ++i)
            data[idx] = 0;
    }
 
    template <typename rng_t>
        requires (data_layout_mode == data_layout::uncompressed)
    void clear(rng_t && bin_range) noexcept
    {
        static_assert(std::ranges::forward_range<rng_t>, "The range of bins to clear must model a forward_range.");
        static_assert(std::same_as<std::remove_cvref_t<std::ranges::range_reference_t<rng_t>>, bin_index>,
                      "The reference type of the range to clear must be seqan3::bin_index.");
#ifndef NDEBUG
        for (auto && bin : bin_range)
            assert(bin.get() < bins);
#endif // NDEBUG
 
        for (size_t offset = 0, i = 0; i < bin_size_; offset += technical_bins, ++i)
            for (auto && bin : bin_range)
                data[bin.get() + offset] = 0;
    }
 
    void increase_bin_number_to(bin_count const new_bins_)
        requires (data_layout_mode == data_layout::uncompressed)
    {
        size_t new_bins = new_bins_.get();
 
        if (new_bins < bins)
            throw std::invalid_argument{"The number of new bins must be >= the current number of bins."};
 
        // Equivalent to ceil(new_bins / 64)
        size_t new_bin_words = (new_bins + 63) >> 6;
 
        bins = new_bins;
 
        if (new_bin_words == bin_words) // No need for internal resize if bin_words does not change.
            return;
 
        size_t new_technical_bins = new_bin_words << 6;
        size_t new_bits = bin_size_ * new_technical_bins;
 
        size_t idx_{new_bits}, idx{data.size()};
        size_t delta = new_technical_bins - technical_bins + 64;
 
        data.resize(new_bits);
 
        for (size_t i = idx_, j = idx; j > 0; i -= new_technical_bins, j -= technical_bins)
        {
            size_t stop = i - new_technical_bins;
 
            for (size_t ii = i - delta, jj = j - 64; stop && ii >= stop; ii -= 64, jj -= 64)
            {
                uint64_t old = data.get_int(jj);
                data.set_int(jj, 0);
                data.set_int(ii, old);
            }
        }
 
        bin_words = new_bin_words;
        technical_bins = new_technical_bins;
    }
 
    membership_agent membership_agent() const
    {
        return typename interleaved_bloom_filter<data_layout_mode>::membership_agent{*this};
    }
 
    template <typename value_t = uint16_t>
    counting_agent_type<value_t> counting_agent() const
    {
        return counting_agent_type<value_t>{*this};
    }
 
    size_t hash_function_count() const noexcept
    {
        return hash_funs;
    }
 
    size_t bin_count() const noexcept
    {
        return bins;
    }
 
    size_t bin_size() const noexcept
    {
        return bin_size_;
    }
 
    size_t bit_size() const noexcept
    {
        return data.size();
    }
 
    friend bool operator==(interleaved_bloom_filter const & lhs, interleaved_bloom_filter const & rhs) noexcept
    {
        return std::tie(lhs.bins, lhs.technical_bins, lhs.bin_size_, lhs.hash_shift, lhs.bin_words, lhs.hash_funs,
                        lhs.data) ==
               std::tie(rhs.bins, rhs.technical_bins, rhs.bin_size_, rhs.hash_shift, rhs.bin_words, rhs.hash_funs,
                        rhs.data);
    }
 
    friend bool operator!=(interleaved_bloom_filter const & lhs, interleaved_bloom_filter const & rhs) noexcept
    {
        return !(lhs == rhs);
    }
 
    constexpr data_type & raw_data() noexcept
    {
        return data;
    }
 
    constexpr data_type const & raw_data() const noexcept
    {
        return data;
    }
 
    template <cereal_archive archive_t>
    void CEREAL_SERIALIZE_FUNCTION_NAME(archive_t & archive)
    {
        archive(bins);
        archive(technical_bins);
        archive(bin_size_);
        archive(hash_shift);
        archive(bin_words);
        archive(hash_funs);
        archive(data);
    }
};
 
template <data_layout data_layout_mode>
class interleaved_bloom_filter<data_layout_mode>::membership_agent
{
private:
    using ibf_t = interleaved_bloom_filter<data_layout_mode>;
 
    ibf_t const * ibf_ptr{nullptr};
 
public:
    class binning_bitvector;
 
    membership_agent() = default; 
    membership_agent(membership_agent const &) = default; 
    membership_agent & operator=(membership_agent const &) = default; 
    membership_agent(membership_agent &&) = default; 
    membership_agent & operator=(membership_agent &&) = default; 
    ~membership_agent() = default; 
 
    explicit membership_agent(ibf_t const & ibf) :
        ibf_ptr(std::addressof(ibf)), result_buffer(ibf.bin_count())
    {}
 
    binning_bitvector result_buffer;
 
    [[nodiscard]] binning_bitvector const & bulk_contains(size_t const value) & noexcept
    {
        assert(ibf_ptr != nullptr);
        assert(result_buffer.size() == ibf_ptr->bin_count());
 
        std::array<size_t, 5> bloom_filter_indices;
        std::memcpy(&bloom_filter_indices, &ibf_ptr->hash_seeds, sizeof(size_t) * ibf_ptr->hash_funs);
 
        for (size_t i = 0; i < ibf_ptr->hash_funs; ++i)
            bloom_filter_indices[i] = ibf_ptr->hash_and_fit(value, bloom_filter_indices[i]);
 
        for (size_t batch = 0; batch < ibf_ptr->bin_words; ++batch)
        {
           size_t tmp{-1ULL};
           for (size_t i = 0; i < ibf_ptr->hash_funs; ++i)
           {
               assert(bloom_filter_indices[i] < ibf_ptr->data.size());
               tmp &= ibf_ptr->data.get_int(bloom_filter_indices[i]);
               bloom_filter_indices[i] += 64;
           }
 
           result_buffer.data.set_int(batch << 6, tmp);
        }
 
        return result_buffer;
    }
 
    // `bulk_contains` cannot be called on a temporary, since the object the returned reference points to
    // is immediately destroyed.
    [[nodiscard]] binning_bitvector const & bulk_contains(size_t const value) && noexcept = delete;
 
};
 
template <data_layout data_layout_mode>
class interleaved_bloom_filter<data_layout_mode>::membership_agent::binning_bitvector
{
private:
    using data_type = sdsl::bit_vector;
    data_type data{};
 
    friend class membership_agent;
 
    template <std::integral value_t>
    friend class counting_vector;
 
public:
    binning_bitvector() = default; 
    binning_bitvector(binning_bitvector const &) = default; 
    binning_bitvector & operator=(binning_bitvector const &) = default; 
    binning_bitvector(binning_bitvector &&) = default; 
    binning_bitvector & operator=(binning_bitvector &&) = default; 
    ~binning_bitvector() = default; 
 
    explicit binning_bitvector(size_t const size) :
        data(size)
    {}
 
    size_t size() const noexcept
    {
        return data.size();
    }
 
    auto begin() noexcept
    {
        return data.begin();
    }
 
    auto begin() const noexcept
    {
        return data.begin();
    }
 
    auto end() noexcept
    {
        return data.end();
    }
 
    auto end() const noexcept
    {
        return data.end();
    }
 
    friend bool operator==(binning_bitvector const & lhs, binning_bitvector const & rhs) noexcept
    {
        return lhs.data == rhs.data;
    }
 
    friend bool operator!=(binning_bitvector const & lhs, binning_bitvector const & rhs) noexcept
    {
        return !(lhs == rhs);
    }
 
#ifdef SEQAN3_DEPRECATED_310
    SEQAN3_DEPRECATED_310 friend bool operator==(binning_bitvector const & lhs, sdsl::bit_vector const & rhs) noexcept
    {
        return lhs.data == rhs;
    }
 
    SEQAN3_DEPRECATED_310 friend bool operator==(sdsl::bit_vector const & lhs, binning_bitvector const & rhs) noexcept
    {
        return lhs == rhs.data;
    }
 
    SEQAN3_DEPRECATED_310 friend bool operator!=(binning_bitvector const & lhs, sdsl::bit_vector const & rhs) noexcept
    {
        return !(lhs.data == rhs);
    }
 
    SEQAN3_DEPRECATED_310 friend bool operator!=(sdsl::bit_vector const & lhs, binning_bitvector const & rhs) noexcept
    {
        return !(lhs == rhs.data);
    }
#endif // SEQAN3_DEPRECATED_310
 
    auto operator[](size_t const i) noexcept
    {
        assert(i < size());
        return data[i];
    }
 
    auto operator[](size_t const i) const noexcept
    {
        assert(i < size());
        return data[i];
    }
 
    constexpr data_type & raw_data() noexcept
    {
        return data;
    }
 
    constexpr data_type const & raw_data() const noexcept
    {
        return data;
    }
};
 
template <std::integral value_t>
class counting_vector : public std::vector<value_t>
{
private:
    using base_t = std::vector<value_t>;
public:
    counting_vector() = default; 
    counting_vector(counting_vector const &) = default; 
    counting_vector & operator=(counting_vector const &) = default; 
    counting_vector(counting_vector &&) = default; 
    counting_vector & operator=(counting_vector &&) = default; 
    ~counting_vector() = default; 
 
    using base_t::base_t;
 
    template <typename rhs_t>
        requires std::same_as<rhs_t,
                              interleaved_bloom_filter<data_layout::uncompressed>::membership_agent::binning_bitvector>
                 ||
                 std::same_as<rhs_t,
                              interleaved_bloom_filter<data_layout::compressed>::membership_agent::binning_bitvector>
    counting_vector & operator+=(rhs_t const & rhs)
    {
        assert(this->size() >= rhs.size()); // The counting vector may be bigger than what we need.
 
        // Each iteration can handle 64 bits, so we need to iterate `((rhs.size() + 63) >> 6` many times
        for (size_t batch = 0, bin = 0; batch < ((rhs.size() + 63) >> 6); bin = 64 * ++batch)
        {
            size_t tmp = rhs.data.get_int(batch * 64); // get 64 bits starting at position `batch * 64`
            if (tmp ^ (1ULL<<63)) // This is a special case, because we would shift by 64 (UB) in the while loop.
            {
                while (tmp > 0)
                {
                    // Jump to the next 1 and increment the corresponding vector entry.
                    uint8_t step = std::countr_zero(tmp);
                    bin += step++;
                    tmp >>= step;
                    ++(*this)[bin++];
                }
            }
            else
            {
                ++(*this)[bin + 63];
            }
        }
        return *this;
    }
 
    counting_vector & operator+=(counting_vector const & rhs)
    {
        assert(this->size() >= rhs.size()); // The counting vector may be bigger than what we need.
 
        std::transform(this->begin(), this->end(), rhs.begin(), this->begin(), std::plus<value_t>());
 
        return *this;
    }
};
 
template <data_layout data_layout_mode>
template <std::integral value_t>
class interleaved_bloom_filter<data_layout_mode>::counting_agent_type
{
private:
    using ibf_t = interleaved_bloom_filter<data_layout_mode>;
 
    ibf_t const * ibf_ptr{nullptr};
 
    typename ibf_t::membership_agent membership_agent;
 
public:
    counting_agent_type() = default; 
    counting_agent_type(counting_agent_type const &) = default; 
    counting_agent_type & operator=(counting_agent_type const &) = default; 
    counting_agent_type(counting_agent_type &&) = default; 
    counting_agent_type & operator=(counting_agent_type &&) = default; 
    ~counting_agent_type() = default; 
 
    explicit counting_agent_type(ibf_t const & ibf) :
        ibf_ptr(std::addressof(ibf)), membership_agent(ibf), result_buffer(ibf.bin_count())
    {}
 
    counting_vector<value_t> result_buffer;
 
    template <std::ranges::range value_range_t>
    [[nodiscard]] counting_vector<value_t> const & bulk_count(value_range_t && values) & noexcept
    {
        assert(ibf_ptr != nullptr);
        assert(result_buffer.size() == ibf_ptr->bin_count());
 
        static_assert(std::ranges::input_range<value_range_t>, "The values must model input_range.");
        static_assert(std::unsigned_integral<std::ranges::range_value_t<value_range_t>>,
                      "An individual value must be an unsigned integral.");
 
        std::ranges::fill(result_buffer, 0);
 
        for (auto && value : values)
            result_buffer += membership_agent.bulk_contains(value);
 
        return result_buffer;
    }
 
    // `bulk_count` cannot be called on a temporary, since the object the returned reference points to
    // is immediately destroyed.
    template <std::ranges::range value_range_t>
    [[nodiscard]] counting_vector<value_t> const & bulk_count(value_range_t && values) && noexcept = delete;
 
};
 
 
} // namespace seqan3