BitString.cuh

Go to the documentation of this file.

 
#pragma once
 
#include <array>
#include <bit>       // for std::countl_zero, std::bit_ceil
#include <compare>   // for std::strong_ordering and <=> operator
#include <cstddef>   // for size_t
#include <cstdint>   // for uintmax_t
#include <numeric>   // Provides std::min (for two arguments)
#include <vector>
#include <algorithm> // Include for std::min/max
 
// For file I/O
#include <fstream>
#include <stdexcept> // for std::runtime_error
#include <string>
 
#include <iostream>
 
// Include CUDA headers for parallel operations
#include "msw.cuh"
#include "cuda_utils.cuh"
 
static constexpr size_t MIN_PAR_COMPARE_WORD_SIZE = 1;
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS = sizeof(CHAR_T) * 8>
class BitString
{
public:
    static constexpr unsigned WORD_SIZE_BITS = sizeof(uintmax_t) * 8;
    static constexpr unsigned ALPHA = WORD_SIZE_BITS / CHAR_SIZE_BITS;
    static constexpr size_t MIN_PAR_COMPARE_CHAR_SIZE = MIN_PAR_COMPARE_WORD_SIZE * ALPHA;
 
    explicit BitString(const CHAR_T* data, size_t size);
 
    template<typename Container>
    explicit BitString(const Container& data);
 
    explicit BitString(const std::string& data)
        : BitString(reinterpret_cast<const CHAR_T*>(data.c_str()), data.size()) // Reinterpret cast assumes CHAR_T layout matches char
    {
        // Body intentionally empty
    }
 
    BitString() = default;
 
    void push_back(CHAR_T c) noexcept;
 
    CHAR_T operator[](size_t index) const;
 
    size_t size() const noexcept { return m_size; }
 
    size_t word_count() const noexcept { return m_data.size(); }
 
    void clear() noexcept
    {
        m_data.clear();
        m_size = 0;
    }
 
    struct ResultLCP
    {
        std::strong_ordering result; 
        size_t lcp;                  
    };
 
    ResultLCP seq_k_compare(const BitString& other, size_t lcp, size_t max_compare) const noexcept;
 
    ResultLCP seq_k_compare(const BitString& other, size_t lcp) const noexcept
    {
        return seq_k_compare(other, lcp, size());
    }
 
    ResultLCP par_k_compare(const BitString& other, size_t lcp, size_t max_compare, uintmax_t* d_a, uintmax_t* d_largeblock, size_t& max_copied) const noexcept;
 
    ResultLCP par_k_compare(const BitString& other, size_t lcp, uintmax_t* d_a, uintmax_t* d_largeblock, size_t& max_copied) const noexcept
    {
        return par_k_compare(other, lcp, size(), d_a, d_largeblock, max_copied);
    }
 
    const uintmax_t* data() const noexcept { return m_data.data(); }
 
private:
    std::vector<uintmax_t> m_data;
    size_t m_size = 0;
 
    static constexpr std::array<unsigned, ALPHA> GET_SHIFTS();
    static constexpr std::array<unsigned, ALPHA> m_shifts = GET_SHIFTS();
 
    static constexpr std::array<uintmax_t, ALPHA> GET_MASKS();
    static constexpr std::array<uintmax_t, ALPHA> m_masks = GET_MASKS();
 
    static constexpr CHAR_T GET_CHAR(uintmax_t word, unsigned bit_index);
 
public:
    // Forward declaration for the friend function template
    template<typename CT, unsigned CSB>
    friend std::ostream& operator<<(std::ostream& os, const BitString<CT, CSB>& bs);
 
    void print_bytes(std::ostream& os) const;
 
    void to_file(const std::string& filename) const;
 
    void to_file(std::ofstream& file) const;
 
    static BitString from_file(const std::string& filename);
 
    static BitString from_file(std::ifstream& file);
 
    std::string to_string() const noexcept;
 
public:
    class Iterator
    {
    public:
        using iterator_category = std::forward_iterator_tag;
        using value_type = CHAR_T;
        using difference_type = std::ptrdiff_t;
        using pointer = const CHAR_T*;   // Or void if value is returned by copy
        using reference = CHAR_T; // Return by value
 
        Iterator(const BitString& bs, size_t index)
            : m_bs(bs), m_index(index)
        {
            // Body intentionally empty
        }
 
        CHAR_T operator*() const noexcept
        {
            return m_bs[m_index];
        }
 
        Iterator& operator++() noexcept
        {
            ++m_index;
            return *this;
        }
 
        Iterator operator++(int) noexcept
        {
            Iterator temp = *this;
            ++(*this);
            return temp;
        }
 
 
        bool operator==(const Iterator& other) const noexcept
        {
            return m_index == other.m_index;
        }
 
    private:
        const BitString& m_bs;
        size_t m_index;
    };
 
    Iterator begin() const noexcept
    {
        return Iterator(*this, 0);
    }
 
    Iterator end() const noexcept
    {
        return Iterator(*this, size());
    }
};
 
// == Implementation of Template Methods ==
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
constexpr std::array<unsigned, BitString<CHAR_T, CHAR_SIZE_BITS>::ALPHA> BitString<CHAR_T, CHAR_SIZE_BITS>::GET_SHIFTS()
{
    std::array<unsigned, ALPHA> shifts;
    for (unsigned i = 0; i < ALPHA; ++i)
    {
        shifts[ALPHA - i - 1] = i * CHAR_SIZE_BITS;
    }
    return shifts;
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
constexpr std::array<uintmax_t, BitString<CHAR_T, CHAR_SIZE_BITS>::ALPHA> BitString<CHAR_T, CHAR_SIZE_BITS>::GET_MASKS()
{
    uintmax_t mask = 0;
    for (unsigned i = 0; i < CHAR_SIZE_BITS; ++i)
    {
        mask |= static_cast<uintmax_t>(1) << i;
    }
 
    std::array<uintmax_t, ALPHA> masks;
    for (unsigned i = 0; i < ALPHA; ++i)
    {
        masks[i] = mask << m_shifts[i];
    }
    return masks;
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
constexpr CHAR_T BitString<CHAR_T, CHAR_SIZE_BITS>::GET_CHAR(uintmax_t word, unsigned bit_index)
{
    const uintmax_t masked_word = word & m_masks[bit_index];
    const uintmax_t shifted_word = masked_word >> m_shifts[bit_index];
    return static_cast<CHAR_T>(shifted_word);
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
BitString<CHAR_T, CHAR_SIZE_BITS>::BitString(const CHAR_T* data, size_t data_size)
    : m_size(data_size)
{
    const size_t word_count = (size() + ALPHA - 1) / ALPHA;
    m_data.resize(word_count, 0);
 
    for (size_t i = 0; i < size(); ++i)
    {
        const size_t word_index = i / ALPHA;
        const size_t bit_index = i % ALPHA;
        const uintmax_t mask = static_cast<uintmax_t>(data[i]) << m_shifts[bit_index];
        m_data[word_index] |= mask;
    }
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
template<typename Container>
BitString<CHAR_T, CHAR_SIZE_BITS>::BitString(const Container& data)
    : m_size(data.size()) // Get size from the container
{
    const size_t word_count = (size() + ALPHA - 1) / ALPHA;
    m_data.resize(word_count, 0);
 
    for (size_t i = 0; i < size(); ++i)
    {
        const size_t word_index = i / ALPHA;
        const size_t bit_index = i % ALPHA;
        const uintmax_t mask = static_cast<uintmax_t>(data[i]) << m_shifts[bit_index];
        m_data[word_index] |= mask;
    }
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
void BitString<CHAR_T, CHAR_SIZE_BITS>::push_back(CHAR_T c) noexcept
{
    const size_t word_index = size() / ALPHA;
    const size_t bit_index = size() % ALPHA;
    const uintmax_t mask = static_cast<uintmax_t>(c) << m_shifts[bit_index];
 
    if (word_index >= m_data.size())
    {
        m_data.push_back(0);
    }
 
    m_data[word_index] |= mask;
    ++m_size;
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
CHAR_T BitString<CHAR_T, CHAR_SIZE_BITS>::operator[](size_t index) const
{
    return GET_CHAR(m_data[index / ALPHA], index % ALPHA);
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
auto BitString<CHAR_T, CHAR_SIZE_BITS>::seq_k_compare(const BitString& other, size_t lcp, size_t max_compare) const noexcept -> ResultLCP
{
    const size_t min_size = std::min(size(), other.size());
    const size_t max_compare_until = lcp + max_compare;
    size_t word_index = lcp / ALPHA;
 
    lcp -= lcp % ALPHA; // Reset to word boundary
 
    while (lcp < min_size && lcp < max_compare_until)
    {
        const uintmax_t XOR = m_data[word_index] ^ other.m_data[word_index];
 
        if (XOR == 0)
        {
            lcp += ALPHA;
            ++word_index;
            continue;
        }
 
        const size_t l_zero = std::countl_zero(XOR);
 
        lcp += l_zero / CHAR_SIZE_BITS;
 
        const size_t bit_index = lcp % ALPHA;
 
        // Mismatch found within limit
        const CHAR_T lhs = GET_CHAR(m_data[word_index], bit_index);
        const CHAR_T rhs = GET_CHAR(other.m_data[word_index], bit_index);
 
        return { lhs <=> rhs, lcp };
    }
 
    if (lcp >= min_size)
    {
        return { size() <=> other.size(), min_size };
    }
 
    return { std::strong_ordering::equal, lcp };
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
auto BitString<CHAR_T, CHAR_SIZE_BITS>::par_k_compare(const BitString& other, size_t lcp, size_t max_compare, uintmax_t* d_a, uintmax_t* d_largeblock, size_t& max_copied) const noexcept -> ResultLCP
{
    // Fallback for small comparisons
    if (max_compare <= MIN_PAR_COMPARE_CHAR_SIZE)
    {
        return seq_k_compare(other, lcp, max_compare);
    }
 
    const size_t min_size = std::min(size(), other.size());
    const size_t min_word_count = std::min(word_count(), other.word_count());
    size_t word_index = lcp / ALPHA;
 
    // Calculate how many full words remain to potentially compare
    const size_t num_remaining_words = (min_word_count > word_index) ? (min_word_count - word_index) : 0;
 
    // Calculate max words to compare based on max_compare characters
    // Need to compare up to the word containing character (lcp + max_compare - 1)
    const size_t end_char_index = lcp + max_compare;
    const size_t end_word_index = (end_char_index + ALPHA - 1) / ALPHA; // Word index containing the last char to check
    const size_t max_compare_words = (end_word_index > word_index) ? (end_word_index - word_index) : 0;
 
 
    const size_t actual_max_compare_words = std::min(max_compare_words, num_remaining_words);
 
    // Fallback if the actual number of words to compare is too small
    if (actual_max_compare_words <= MIN_PAR_COMPARE_WORD_SIZE)
    {
        return seq_k_compare(other, lcp, max_compare);
    }
 
    // Ensure data is on the device
    if (max_copied < word_index + actual_max_compare_words)
    {
        size_t end_copy_index = word_index + actual_max_compare_words;
        
        // get next power of 2 greater than or equal to end_copy_index
        size_t end_copy = std::min(std::bit_ceil(end_copy_index), word_count());
        if (end_copy > max_copied) {
            copy_to_device(d_a + max_copied, data() + max_copied, end_copy - max_copied);
            max_copied = end_copy;
        }
    }
 
    lcp -= lcp % ALPHA; // Reset LCP to the start of the current word boundary for parallel comparison
 
    // Perform parallel mismatch search
    auto msw = par_find_mismatch(d_a + word_index, other.data() + word_index, d_largeblock, actual_max_compare_words);
 
    // Calculate LCP based on mismatch word index (msw)
    lcp += msw * ALPHA;
 
    // If no mismatch found within the compared words
    if (msw == actual_max_compare_words)
    {
        // Check if the comparison reached the end of either string
        if (lcp >= min_size)
        {
            return { size() <=> other.size(), min_size };
        }
 
        return { std::strong_ordering::equal, lcp };
    }
 
    // Mismatch found by parallel kernel. Refine LCP within the mismatch word.
    word_index = lcp / ALPHA; // Update word_index to the mismatch word
 
    // Find exact mismatch character within the word
    const uintmax_t XOR = m_data[word_index] ^ other.m_data[word_index];
    const size_t l_zero = std::countl_zero(XOR);
 
    lcp += l_zero / CHAR_SIZE_BITS;
 
    const size_t bit_index = lcp % ALPHA;
 
    // Mismatch is within the limit
    const CHAR_T lhs = GET_CHAR(m_data[word_index], bit_index);
    const CHAR_T rhs = GET_CHAR(other.m_data[word_index], bit_index);
 
    return { lhs <=> rhs, lcp };
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
void BitString<CHAR_T, CHAR_SIZE_BITS>::print_bytes(std::ostream& os) const
{
    static constexpr unsigned shift = WORD_SIZE_BITS - 8;
 
    for (size_t i = 0, a = 0; i + ALPHA <= size(); i += ALPHA, ++a)
    {
        uintmax_t word = m_data[a];
        
        for (unsigned c = 0; c < WORD_SIZE_BITS / 8; c++, word <<= 8)
        {
            os << static_cast<uint8_t>(word >> shift);
        }
    }
 
    // Handle last partial word if it exists
    if (size() % ALPHA != 0 && !m_data.empty())
    {
        uintmax_t word = m_data.back();
        unsigned remaining_chars = size() % ALPHA;
        unsigned remaining_bits = remaining_chars * CHAR_SIZE_BITS;
        unsigned bytes_to_print = (remaining_bits + 7) / 8; // Round up bytes
 
        // Process the required bytes from the most significant side
        for (unsigned byte_idx = 0; byte_idx < bytes_to_print; ++byte_idx)
        {
            os << static_cast<uint8_t>(word >> (shift - byte_idx * 8));
        }
    }
}
 
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
std::ostream& operator<<(std::ostream& os, const BitString<CHAR_T, CHAR_SIZE_BITS>& bs)
{
    for (CHAR_T c : bs)
    {
        os << c;
    }
    return os;
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
void BitString<CHAR_T, CHAR_SIZE_BITS>::to_file(std::ofstream& file) const
{
    file.write(reinterpret_cast<const char*>(&m_size), sizeof(size_t));
    file.write(reinterpret_cast<const char*>(m_data.data()), m_data.size() * sizeof(uintmax_t));
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
void BitString<CHAR_T, CHAR_SIZE_BITS>::to_file(const std::string& filename) const
{
    std::ofstream file(filename, std::ios::binary);
    if (!file)
    {
        throw std::runtime_error("Failed to open file for writing: " + filename);
    }
    to_file(file);
    // file automatically closed by destructor
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
BitString<CHAR_T, CHAR_SIZE_BITS> BitString<CHAR_T, CHAR_SIZE_BITS>::from_file(std::ifstream& file)
{
    BitString<CHAR_T, CHAR_SIZE_BITS> bs;
 
    file.read(reinterpret_cast<char*>(&bs.m_size), sizeof(size_t));
    if (!file) { // Check read success
        throw std::runtime_error("Failed to read size from file stream.");
    }
 
    const size_t word_count = (bs.size() + ALPHA - 1) / ALPHA;
    bs.m_data.resize(word_count, 0); // Resize before reading into it
 
    if (word_count > 0) { // Only read if there's data expected
        file.read(reinterpret_cast<char*>(bs.m_data.data()), word_count * sizeof(uintmax_t));
        if (!file) { // Check read success
            throw std::runtime_error("Failed to read data from file stream.");
        }
    }
 
 
    return bs;
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
BitString<CHAR_T, CHAR_SIZE_BITS> BitString<CHAR_T, CHAR_SIZE_BITS>::from_file(const std::string& filename)
{
    std::ifstream file(filename, std::ios::binary);
    if (!file)
    {
        throw std::runtime_error("Failed to open file for reading: " + filename);
    }
    return from_file(file);
    // file automatically closed by destructor
}
 
template<typename CHAR_T, unsigned CHAR_SIZE_BITS>
std::string BitString<CHAR_T, CHAR_SIZE_BITS>::to_string() const noexcept
{
    std::string str;
    str.reserve(size());
 
    for (CHAR_T c : *this)
    {
        // Assumes CHAR_T is convertible to char, potential truncation if wider.
        str.push_back(static_cast<char>(c));
    }
 
    return str;
}