sdsl/include/sdsl/wt_int.hpp

Go to the documentation of this file.

/* sdsl - succinct data structures library
    Copyright (C) 2009 Simon Gog

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see http://www.gnu.org/licenses/ .
*/
#ifndef INCLUDED_SDSL_INT_WAVELET_TREE
#define INCLUDED_SDSL_INT_WAVELET_TREE

#include "int_vector.hpp"
#include "rank_support_v.hpp"
#include "select_support_mcl.hpp"
#include "bitmagic.hpp"
#include "testutils.hpp"
#include "temp_write_read_buffer.hpp"
#include "util.hpp"
#include <set> // for calculating the alphabet size
#include <map> // for mapping a symbol to its lexicographical index
#include <algorithm> // for std::swap
#include <stdexcept>
#include <vector>

namespace sdsl
{


template<class RandomAccessContainer=int_vector<>,
         class BitVector                 = bit_vector,
         class RankSupport               = typename BitVector::rank_1_type,
         class SelectSupport     = typename BitVector::select_1_type,
         class SelectSupportZero = typename BitVector::select_0_type>
class wt_int
{
    public:
        typedef typename RandomAccessContainer::size_type               size_type;
        typedef typename RandomAccessContainer::value_type              value_type;
        typedef BitVector                                                                               bit_vector_type;
        typedef RankSupport                                                                             rank_1_type;
        typedef SelectSupport                                                                   select_1_type;
        typedef SelectSupportZero                                                               select_0_type;
    protected:
        size_type                       m_size;
        size_type                       m_sigma;                //<- \f$ |\Sigma| \f$
        bit_vector_type         m_tree;                 // bit vector to store the wavelet tree
        rank_1_type                     m_tree_rank;    // rank support for the wavelet tree bit vector
        select_1_type           m_tree_select1; // select support for the wavelet tree bit vector
        select_0_type           m_tree_select0;
        uint32_t                        m_logn;

        void copy(const wt_int& wt) {
            m_size                      = wt.m_size;
            m_sigma             = wt.m_sigma;
            m_tree                      = wt.m_tree;
            m_tree_rank         = wt.m_tree_rank;
            m_tree_rank.set_vector(&m_tree);
            m_tree_select1      = wt.m_tree_select1;
            m_tree_select1.set_vector(&m_tree);
            m_tree_select0      = wt.m_tree_select0;
            m_tree_select0.set_vector(&m_tree);
            m_logn                      = wt.m_logn;
        }

    public:

        const size_type& sigma; 
        const bit_vector_type& tree; 

        wt_int():m_size(0),m_sigma(0), m_logn(0), sigma(m_sigma), tree(m_tree) {};



        wt_int(const RandomAccessContainer& rac, uint32_t logn=0):m_size(rac.size()), m_sigma(0), sigma(m_sigma), tree(m_tree) {
#ifdef SDSL_DEBUG_INT_WAVELET_TREE
            std::cerr<<"wt_int construct: size="<<m_size<<std::endl;
#endif
            size_type n =       m_size;
            if (logn == 0) { // detect the biggest value in rac and set logn properly
                value_type x = 1;
                for (size_type i=0; i < n; ++i)
                    if (rac[i] > x)
                        x = rac[i];
                m_logn  = bit_magic::l1BP(x)+1; // we need logn bits to represent all values in the range [0..x]
            } else {
                m_logn = logn;
            }
            bit_vector tree(n*m_logn, 0);  // initialize the tree

            int_vector<>        perms[2];
            perms[0].set_int_width(m_logn); perms[1].set_int_width(m_logn);
            perms[0].resize(n); perms[1].resize(n);

            for (size_type i=0; i < n; ++i) { // copy original array to perms[0]
                perms[0][i] = rac[i];
            }

            size_type tree_pos = 0;

            uint64_t            mask_old = 1ULL<<(m_logn);
            for (uint32_t k=0; k<m_logn; ++k) {
                int_vector<>&   act_perm         = perms[k%2];
                int_vector<>&   next_perm        = perms[1-(k%2)];
                size_type               start            = 0;
                const uint64_t          mask_new = 1ULL<<(m_logn-k-1);
                do {
                    size_type   i               = start;
                    size_type   cnt1    =       0;
                    uint64_t    start_value = (act_perm[i]&mask_old);
                    uint64_t    x;
                    while (i < n and((x=act_perm[i])&mask_old)==start_value) {
                        if (x&mask_new) {
                            ++cnt1; tree[tree_pos] = 1;
                        }
                        ++tree_pos;
                        ++i;
                    }
                    size_type cnt0 = i-start-cnt1;
                    if (k+1 < m_logn) { // inner node
                        size_type idxs[2] = {start, start+cnt0};
                        for (i=start, start = start+cnt0+cnt1; i < start; ++i) {
                            next_perm[idxs[(act_perm[i]&mask_new)>0 ]++] = act_perm[i];
                        }
                    } else { // leaf node
                        start += cnt0+cnt1;
                        ++m_sigma;  // increase sigma for each leaf
                    }
                } while (start < n);
                mask_old += mask_new;
            }
#ifdef SDSL_DEBUG_INT_WAVELET_TREE
            if (tree.size()<100) {
                std::cerr<<"tree="<<tree<<std::endl;
            }
#endif
            util::assign(m_tree, tree);
            util::init_support(m_tree_rank, &m_tree);
            util::init_support(m_tree_select0, &m_tree);
            util::init_support(m_tree_select1, &m_tree);
        }


        template<uint8_t int_width>
        wt_int(int_vector_file_buffer<int_width>& buf, uint32_t logn=0, std::string dir="./")
            : m_size(0),m_sigma(0), m_logn(0), sigma(m_sigma), tree(m_tree) {
            construct(buf, logn, dir);
        }


        // TODO external construction method which occupies only additional constant memory
        template<uint8_t int_width>
        void construct(int_vector_file_buffer<int_width>& buf, uint32_t logn=0, std::string dir="./") {
            buf.reset();
            size_type n = buf.int_vector_size;  // set n
            m_size = n;                         // set sigma and size
            m_sigma = 0;
            temp_write_read_buffer<> buf1(5000000, buf.int_width, dir);   // buffer for elements in the rigth node
            int_vector<int_width> rac(n, 0, buf.int_width);                               // initialze rac

            value_type x = 1;  // variable for the biggest value in rac
            for (size_type i=0,r=0,r_sum=0; i<n;) { // detect the biggest value in rac
                for (; i < r+r_sum; ++i) {
                    if (buf[i-r_sum] > x)
                        x = buf[i-r_sum];
                    rac[i] = buf[i-r_sum];
                }
                r_sum += r; r = buf.load_next_block();
            }

            if (logn == 0) {
                m_logn  = bit_magic::l1BP(x)+1; // we need logn bits to represent all values in the range [0..x]
            } else {
                m_logn = logn;
            }
            std::string tree_out_buf_file_name = (dir+"m_tree"+util::to_string(util::get_pid())+"_"+util::to_string(util::get_id()));
            std::ofstream tree_out_buf(tree_out_buf_file_name.c_str(),
                                       std::ios::binary | std::ios::trunc | std::ios::out);   // open buffer for tree
            size_type bit_size = n*m_logn;
            tree_out_buf.write((char*) &bit_size, sizeof(bit_size));    // write size of bit_vector

            size_type tree_pos = 0;
            uint64_t tree_word = 0;

            uint64_t            mask_old = 1ULL<<(m_logn);
            for (uint32_t k=0; k<m_logn; ++k) {
                size_type               start            = 0;
                const uint64_t          mask_new = 1ULL<<(m_logn-k-1);
                do {
                    buf1.reset();
                    size_type   i               = start;
                    size_type   cnt0    =       0;
                    uint64_t    start_value = (rac[i]&mask_old);
                    uint64_t    x;
                    while (i < n and((x=rac[i])&mask_old)==start_value) {
                        if (x&mask_new) {
                            tree_word |= (1ULL << (tree_pos&0x3FULL));
                            buf1 << x;
                        } else {
                            rac[start + cnt0++ ] = x;
                        }
                        ++tree_pos;
                        if ((tree_pos & 0x3FULL) == 0) { // if tree_pos % 64 == 0 write old word
                            tree_out_buf.write((char*) &tree_word, sizeof(tree_word));
                            tree_word = 0;
                        }
                        ++i;
                    }
                    buf1.write_close();
                    size_type cnt1 = i-start-cnt0;
                    if (k+1 < m_logn) { // inner node
                        for (i=start + cnt0, start = start+cnt0+cnt1; i < start; ++i) {
                            buf1 >> x;
                            rac[ i ] = x;
                        }
                    } else { // leaf node
                        start += cnt0+cnt1;
                        ++m_sigma; // increase sigma for each leaf
                    }

                } while (start < n);
                mask_old += mask_new;
            }
            if ((tree_pos & 0x3FULL) != 0) { // if tree_pos % 64 > 0 => there are remaining entries we have to write
                tree_out_buf.write((char*) &tree_word, sizeof(tree_word));
            }
            tree_out_buf.close();
            rac.resize(0);
            bit_vector tree;
            util::load_from_file(tree, tree_out_buf_file_name.c_str());
            std::remove(tree_out_buf_file_name.c_str());
#ifdef SDSL_DEBUG_INT_WAVELET_TREE
            if (tree.size()<100) {
                std::cerr<<"tree="<<tree<<std::endl;
            }
#endif
            util::assign(m_tree, tree);
            util::init_support(m_tree_rank, &m_tree);
            util::init_support(m_tree_select0, &m_tree);
            util::init_support(m_tree_select1, &m_tree);
        }

        wt_int(const wt_int& wt):sigma(m_sigma) {
            copy(wt);
        }

        wt_int& operator=(const wt_int& wt) {
            if (this != &wt) {
                copy(wt);
            }
            return *this;
        }

        void swap(wt_int& wt) {
            if (this != &wt) {
                std::swap(m_size, wt.m_size);
                std::swap(m_sigma,  wt.m_sigma);
                m_tree.swap(wt.m_tree);
                m_tree_rank.swap(wt.m_tree_rank); // rank swap after the swap of the bit vector m_tree
                m_tree_rank.set_vector(&m_tree);
                m_tree_select1.swap(wt.m_tree_select1); // select1 swap after the swap of the bit vector m_tree
                m_tree_select1.set_vector(&m_tree);
                m_tree_select0.swap(wt.m_tree_select0); // select0 swap after the swap of the bit vector m_tree
                m_tree_select0.set_vector(&m_tree);
                std::swap(m_logn,  wt.m_logn);
            }
        }

        size_type size()const {
            return m_size;
        }

        bool empty()const {
            return m_size == 0;
        }


        value_type operator[](size_type i)const {
//              size_type zeros_left_of_pos = 0;
            size_type offset = 0;
            value_type res = 0;
            size_type node_size = m_size;
            for (uint32_t k=0; k < m_logn; ++k) {
                res <<= 1;
                size_type ones_before_o   = m_tree_rank(offset);
                size_type ones_before_i   = m_tree_rank(offset + i) - ones_before_o;
                size_type ones_before_end = m_tree_rank(offset + node_size) - ones_before_o;
//                      std::cerr<<"k="<<k<<" i="<<i<<" offset="<<offset<<" node_size="<<node_size<<" obo="<<ones_before_o<<" obi="<<ones_before_i<<" obe="<<ones_before_end<<std::endl;
                if (m_tree[offset+i]) { // one at position i => follow right child
//                              zeros_left_of_pos += (node_size - ones_before_end);
                    offset += (node_size - ones_before_end);
                    node_size = ones_before_end;
                    i = ones_before_i;
                    res |= 1;
                } else { // zero at position i => follow left child
                    node_size = (node_size - ones_before_end);
                    i = (i-ones_before_i);
                }
                //offset = (k+1)*m_size + zeros_left_of_pos;
                offset += m_size;
            }
//              std::cerr<<" offset="<<offset<<" node_size="<<node_size<<std::endl;
            return res;
        };

        //TODO: buchstaben, die nicht im orginaltext vorkommen behandeln!!!

        size_type rank(size_type i, value_type c)const {
            size_type offset = 0;
            uint64_t mask        = (1ULL) << (m_logn-1);
            size_type node_size = m_size;
            for (uint32_t k=0; k < m_logn and i; ++k) {
                size_type ones_before_o   = m_tree_rank(offset);
                size_type ones_before_i   = m_tree_rank(offset + i) - ones_before_o;
                size_type ones_before_end = m_tree_rank(offset + node_size) - ones_before_o;
                if (c & mask) { // search for a one at this level
                    offset += (node_size - ones_before_end);
                    node_size = ones_before_end;
                    i = ones_before_i;
                } else { // search for a zero at this level
                    node_size = (node_size - ones_before_end);
                    i = (i-ones_before_i);
                }
                offset += m_size;
                mask >>= 1;
            }
            return i;
        };


        // TODO: was ist wenn c gar nicht vorkommt, oder es keine i Vorkommen gibt?
        size_type select(size_type i, value_type c)const {
            // possible optimization: if the array is a permutation we can start at the bottom of the tree
            size_type offset = 0;
            uint64_t mask        = (1ULL) << (m_logn-1);
            size_type node_size = m_size;
            size_type offsets[m_logn+1];                        // offsets down the path
            size_type ones_before_os[m_logn+1];   // ones before the offsets
            offsets[0] = ones_before_os[0] = 0;

            for (uint32_t k=0; k < m_logn and node_size; ++k) {
                size_type ones_before_o   = m_tree_rank(offset);
                ones_before_os[k] = ones_before_o;
                size_type ones_before_end = m_tree_rank(offset + node_size) - ones_before_o;
                if (c & mask) { // search for a one at this level
                    offset += (node_size - ones_before_end);
                    node_size = ones_before_end;
                } else { // search for a zero at this level
                    node_size = (node_size - ones_before_end);
                }
                offset += m_size;
                offsets[k+1] = offset;
                mask >>= 1;
            }
            if (node_size < i) {
                std::cerr<<"c="<<c<<" does not occure "<<i<<" times in the WT"<<std::endl;
                return m_size;
            }
            mask = 1ULL;
            for (uint32_t k=m_logn; k>0; --k) {
                offset = offsets[k-1];
                size_type ones_before_o = ones_before_os[k-1];
                if (c & mask) { // right child => search i'th
                    i = m_tree_select1(ones_before_o + i) - offset + 1;
                } else { // left child => search i'th zero
                    i = m_tree_select0(offset - ones_before_o + i) - offset + 1;
                }
                mask <<= 1;
            }
            return i-1;
        };


        size_type range_search_2d(size_type lb, size_type rb, value_type vlb, value_type vrb,
                                  std::vector<size_type>* idx_result=NULL,
                                  std::vector<value_type>* val_result=NULL
                                 ) const {
            size_type offsets[m_logn+1];
            size_type ones_before_os[m_logn+1];
            offsets[0] = 0;
            if (vrb > (1ULL << m_logn))
                vrb = (1ULL << m_logn);
            if (vlb > vrb)
                return 0;
            size_type cnt_answers = 0;
            _range_search_2d(lb, rb, vlb, vrb, 0, 0, m_size, offsets, ones_before_os, 0, idx_result, val_result, cnt_answers);
            return cnt_answers;
        }

        // add parameter path
        // ilb interval left bound
        // irb interval right bound
        void _range_search_2d(size_type lb, size_type rb, value_type vlb, value_type vrb, size_type depth,
                              size_type ilb, size_type node_size, size_type offsets[], size_type ones_before_os[], size_type path,
                              std::vector<size_type>* idx_result, std::vector<size_type>* val_result, size_type& cnt_answers)
        const {
            if (lb > rb)
                return;
            if (depth == m_logn) {
                if (idx_result != NULL) {
                    for (size_type j=1; j <= node_size; ++j) {
                        size_type i = j;
                        size_type c = path;
                        for (uint32_t k=m_logn; k>0; --k) {
                            size_type offset = offsets[k-1];
                            size_type ones_before_o = ones_before_os[k-1];
                            if (c&1) {
                                i = m_tree_select1(ones_before_o + i) - offset + 1;
                            } else {
                                i = m_tree_select0(offset - ones_before_o + i) - offset + 1;
                            }
                            c >>= 1;
                        }
                        idx_result->push_back(i-1); // add resulting index; -1 cause of 0 based indexing
                    }
                }
                if (val_result != NULL) {
                    for (size_type j=1; j <= node_size; ++j) {
                        val_result->push_back(path);
                    }
                }
                cnt_answers += node_size;
                return;
            }
            size_type irb = ilb + (1ULL << (m_logn-depth));
            size_type mid = (irb + ilb)>>1;

            size_type offset = offsets[depth];

            size_type ones_before_o             = m_tree_rank(offset);
            ones_before_os[depth]               = ones_before_o;
            size_type ones_before_lb    = m_tree_rank(offset + lb);
            size_type ones_before_rb    = m_tree_rank(offset + rb + 1);
            size_type ones_before_end   = m_tree_rank(offset + node_size);
            size_type zeros_before_o    = offset - ones_before_o;
            size_type zeros_before_lb   = offset + lb - ones_before_lb;
            size_type zeros_before_rb   = offset + rb + 1 - ones_before_rb;
            size_type zeros_before_end  = offset + node_size - ones_before_end;
            if (vlb < mid and mid) {
                size_type nlb                   = zeros_before_lb - zeros_before_o;
                size_type nrb                   = zeros_before_rb - zeros_before_o;
                offsets[depth+1]                = offset + m_size;
                if (nrb)
                    _range_search_2d(nlb, nrb-1, vlb, std::min(vrb,mid-1), depth+1, ilb, zeros_before_end - zeros_before_o, offsets, ones_before_os, path<<1, idx_result, val_result, cnt_answers);
            }
            if (vrb >= mid) {
                size_type nlb                   = ones_before_lb - ones_before_o;
                size_type nrb                   = ones_before_rb - ones_before_o;
                offsets[depth+1]                = offset + m_size + (zeros_before_end - zeros_before_o);
                if (nrb)
                    _range_search_2d(nlb, nrb-1, std::max(mid, vlb), vrb, depth+1, mid, ones_before_end - ones_before_o, offsets, ones_before_os, (path<<1)+1 ,idx_result, val_result, cnt_answers);
            }
        }

        size_type serialize(std::ostream& out, structure_tree_node* v=NULL, std::string name="")const {
            structure_tree_node* child = structure_tree::add_child(v, name, util::class_name(*this));
            size_type written_bytes = 0;
            written_bytes += util::write_member(m_size, out, child, "size");
            written_bytes += util::write_member(m_sigma, out, child, "sigma");
            written_bytes += m_tree.serialize(out, child, "tree");
            written_bytes += m_tree_rank.serialize(out, child, "tree_rank");
            written_bytes += m_tree_select1.serialize(out, child, "tree_select_1");
            written_bytes += m_tree_select0.serialize(out, child, "tree_select_0");
            written_bytes += util::write_member(m_logn, out, child, "log_n");
            structure_tree::add_size(child, written_bytes);
            return written_bytes;
        }

        void load(std::istream& in) {
            util::read_member(m_size, in);
            util::read_member(m_sigma, in);
            m_tree.load(in);
            m_tree_rank.load(in, &m_tree);
            m_tree_select1.load(in, &m_tree);
            m_tree_select0.load(in, &m_tree);
            util::read_member(m_logn, in);
        }
        /*
                void print_info()const{
                        size_type rle_ed = 0;
                        for(size_type i=0; i < m_tree.size(); ++i){

                        }
                }
        */

};

}// end namespace sds

#endif // end file