interleaved_bloom_filter.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2023, 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 <algorithm>
#include <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_type; // 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::compressed> const & ibf)
        requires (data_layout_mode == data_layout::uncompressed)
    {
        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::bit_vector{ibf.data.begin(), ibf.data.end()};
    }
 
    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_type membership_agent() const
    {
        return membership_agent_type{*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_type
{
private:
    using ibf_t = interleaved_bloom_filter<data_layout_mode>;
 
    ibf_t const * ibf_ptr{nullptr};
 
public:
    class binning_bitvector;
 
    membership_agent_type() = default;                                          
    membership_agent_type(membership_agent_type const &) = default;             
    membership_agent_type & operator=(membership_agent_type const &) = default; 
    membership_agent_type(membership_agent_type &&) = default;                  
    membership_agent_type & operator=(membership_agent_type &&) = default;      
    ~membership_agent_type() = default;                                         
 
    explicit membership_agent_type(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_type::binning_bitvector
{
private:
    using data_type = sdsl::bit_vector;
    data_type data{};
 
    friend class membership_agent_type;
 
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);
    }
 
    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>;
 
    template <typename binning_bitvector_t>
    static constexpr bool is_binning_bitvector =
        std::same_as<binning_bitvector_t,
                     interleaved_bloom_filter<data_layout::uncompressed>::membership_agent_type::binning_bitvector>
        || std::same_as<binning_bitvector_t,
                        interleaved_bloom_filter<data_layout::compressed>::membership_agent_type::binning_bitvector>;
 
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 binning_bitvector_t>
        requires is_binning_bitvector<binning_bitvector_t>
    counting_vector & operator+=(binning_bitvector_t const & binning_bitvector)
    {
        for_each_set_bin(binning_bitvector,
                         [this](size_t const bin)
                         {
                             ++(*this)[bin];
                         });
        return *this;
    }
 
    template <typename binning_bitvector_t>
        requires is_binning_bitvector<binning_bitvector_t>
    counting_vector & operator-=(binning_bitvector_t const & binning_bitvector)
    {
        for_each_set_bin(binning_bitvector,
                         [this](size_t const bin)
                         {
                             assert((*this)[bin] > 0);
                             --(*this)[bin];
                         });
        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;
    }
 
    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(),
                       [](auto a, auto b)
                       {
                           assert(a >= b);
                           return a - b;
                       });
 
        return *this;
    }
 
private:
    template <typename binning_bitvector_t, typename on_bin_fn_t>
    void for_each_set_bin(binning_bitvector_t && binning_bitvector, on_bin_fn_t && on_bin_fn)
    {
        assert(this->size() >= binning_bitvector.size()); // The counting vector may be bigger than what we need.
 
        // Jump to the next 1 and return the number of places jumped in the bit_sequence
        auto jump_to_next_1bit = [](size_t & x)
        {
            auto const zeros = std::countr_zero(x);
            x >>= zeros; // skip number of zeros
            return zeros;
        };
 
        // Each iteration can handle 64 bits
        for (size_t bit_pos = 0; bit_pos < binning_bitvector.size(); bit_pos += 64)
        {
            // get 64 bits starting at position `bit_pos`
            size_t bit_sequence = binning_bitvector.raw_data().get_int(bit_pos);
 
            // process each relative bin inside the bit_sequence
            for (size_t bin = bit_pos; bit_sequence != 0u; ++bin, bit_sequence >>= 1)
            {
                // Jump to the next 1 and
                bin += jump_to_next_1bit(bit_sequence);
 
                on_bin_fn(bin);
            }
        }
    }
};
 
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};
 
    membership_agent_type 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