wavelet-matrix.hpp

Go to the documentation of this file.

#if !defined(HSDS_WAVELET_MATRIX_HPP_)
#define HSDS_WAVELET_MATRIX_HPP_

#include <vector>
#include <queue>
#include "hsds/scoped_ptr.hpp"
#include "hsds/bit-vector.hpp"
#include "hsds/vector.hpp"

namespace hsds {

bool uint2bit(uint64_t bits, uint64_t pos) {
    return ((bits >> (sizeof(uint64_t) * 8 - 1 - pos)) & 0x1ULL) == 0x1ULL;
}

// forward declaration
class Exception;

struct ListResult {
    ListResult(uint64_t c, uint64_t freq) :
            c(c), freq(freq) {
    }
    uint64_t c;     
    uint64_t freq;  

    bool operator <(const ListResult& lr) const {
        if (c != lr.c)
            return c < lr.c;
        return freq < lr.freq;
    }
};

class WaveletMatrix {
public:

    WaveletMatrix();

    virtual ~WaveletMatrix();

    void clear();

    void swap(WaveletMatrix& x);

    void build(std::vector<uint64_t>& src);

    inline uint64_t size() const {
        return size_;
    }

    inline uint64_t operator[](uint64_t pos) {
        return lookup(pos);
    }

    uint64_t lookup(uint64_t pos) const;

    uint64_t rank(uint64_t c, uint64_t pos) const;

    uint64_t rankLessThan(uint64_t c, uint64_t pos) const;

    uint64_t rankMoreThan(uint64_t c, uint64_t pos) const;

    void rankAll(uint64_t c, uint64_t begin_pos, uint64_t end_pos, uint64_t& rank, uint64_t& rank_less_than,
            uint64_t& rank_more_than) const;

    uint64_t select(uint64_t c, uint64_t rank) const;

    uint64_t selectFromPos(uint64_t c, uint64_t pos, uint64_t rank) const;

    uint64_t freq(uint64_t c) const;

    uint64_t freqSum(uint64_t min_c, uint64_t max_c) const;

    uint64_t freqRange(uint64_t min_c, uint64_t max_c, uint64_t begin_pos, uint64_t end_pos) const;

    void maxRange(uint64_t begin_pos, uint64_t end_pos, uint64_t& pos, uint64_t& val) const;

    void minRange(uint64_t begin_pos, uint64_t end_pos, uint64_t& pos, uint64_t& val) const;

    void quantileRange(uint64_t begin_pos, uint64_t end_pos, uint64_t k, uint64_t& pos, uint64_t& val) const;

    void listModeRange(uint64_t min_c, uint64_t max_c, uint64_t beg_pos, uint64_t end_pos, uint64_t num,
            std::vector<ListResult>& res) const;

    void listMinRange(uint64_t min_c, uint64_t max_c, uint64_t beg_pos, uint64_t end_pos, uint64_t num,
            std::vector<ListResult>& res) const;

    void listMaxRange(uint64_t min_c, uint64_t max_c, uint64_t beg_pos, uint64_t end_pos, uint64_t num,
            std::vector<ListResult>& res) const;

    void save(std::ostream& os) const throw (hsds::Exception);

    void load(std::istream& is) throw (hsds::Exception);

    uint64_t map(void* ptr, uint64_t mapSize) throw (hsds::Exception);

private:
    typedef hsds::Vector<hsds::BitVector> bv_type;
    typedef hsds::Vector<hsds::Vector<uint64_t> > range_type;
    typedef hsds::Vector<uint64_t> uint64_vector_type;

    uint64_t size_;
    const uint64_t bitSize_;
    uint64_t alphabetNum_;
    uint64_t alphabetBitNum_;
    bv_type bv_;
    range_type nodePos_;
    uint64_vector_type seps_;

    inline uint64_t bitSize() const {
        return bitSize_;
    }

    uint64_t getAlphabetNum(const std::vector<uint64_t>& array) const;
    uint64_t log2(uint64_t x) const;

    struct QueryOnNode {
        QueryOnNode(uint64_t beg_node, uint64_t end_node, uint64_t beg_pos, uint64_t end_pos, uint64_t depth,
                uint64_t prefix_char) :
                beg_node(beg_node),
                    end_node(end_node),
                    beg_pos(beg_pos),
                    end_pos(end_pos),
                    depth(depth),
                    prefix_char(prefix_char) {
        }
        uint64_t beg_node;
        uint64_t end_node;
        uint64_t beg_pos;
        uint64_t end_pos;
        uint64_t depth;
        uint64_t prefix_char;
        void print() {
            std::cout << beg_node << " " << end_node << " " << beg_pos << " " << end_pos << " " << depth << " ";
            for (uint64_t i = 0; i < depth; ++i) {
                std::cout << ((prefix_char >> (depth - (i + 1))) & 1LLU);
            }
            std::cout << std::endl;
        }
    };

    class ListModeComparator;
    class ListMinComparator;
    class ListMaxComparator;

    template<class Comparator>
    void listRange(uint64_t min_c, uint64_t max_c, uint64_t beg_pos, uint64_t end_pos, uint64_t num,
            std::vector<ListResult>& res) const {
        res.clear();
        if (end_pos > size_ || beg_pos >= end_pos)
            return;

        std::priority_queue<QueryOnNode, std::vector<QueryOnNode>, Comparator> qons;
        qons.push(QueryOnNode(0, size_, beg_pos, end_pos, 0, 0));

        while (res.size() < num && !qons.empty()) {
            QueryOnNode qon = qons.top();
            qons.pop();
            if (qon.depth >= alphabetBitNum_) {
                res.push_back(ListResult(qon.prefix_char, qon.end_pos - qon.beg_pos));
            } else {
                std::vector<QueryOnNode> next;
                expandNode(min_c, max_c, qon, next);
                for (size_t i = 0; i < next.size(); ++i) {
                    qons.push(next[i]);
                }
            }
        }
    }

    uint64_t prefixCode(uint64_t x, uint64_t len, uint64_t bit_num) const;
    bool checkPrefix(uint64_t prefix, uint64_t depth, uint64_t min_c, uint64_t max_c) const;
    void expandNode(uint64_t min_c, uint64_t max_c, const QueryOnNode& qon, std::vector<QueryOnNode>& next) const;
};

}

#endif /* !defined(HSDS_WAVELET_MATRIX_HPP_) */