sdsl/include/sdsl/cst3_sct3p.hpp

/* sdsl - succinct data structures library
    Copyright (C) 2010 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_CST_SCT3P
#define INCLUDED_SDSL_CST_SCT3P

#include "int_vector.hpp"
#include "algorithms.hpp"
#include "iterators.hpp"
#include "lcp_support_tree.hpp"
#include "lcp_wt.hpp"
#include "bp_support.hpp"
#include "csa_wt.hpp" // for std initialization of cst_sct3p 
#include "csa_uncompressed.hpp"
#include "cst_iterators.hpp"
#include "cst_sct.hpp" // for lcp_interval
#include "rank_support.hpp"
#include "testutils.hpp"
#include "util.hpp"
#include <iostream>
#include <algorithm>
#include <cassert>
#include <cstring> // for strlen
#include <iomanip>
#include <iterator>
#include <stack> // for the calculation of the balanced parantheses sequence
#include <ostream>

namespace sdsl
{

struct cst_tag; // forward declaration

template<class Int = int_vector<>::size_type /*size_t*/>
struct bp_interval_p
{
    Int ipos;  // position of the i+1th opening parenthesis in the balanced parenthese sequence
    Int cipos; // positon of the matching closing parenthesis of the i+1th openeing parenthesis in the balanced parentheses sequence
    Int jp1pos; // position of the j+2th opening parentheis in the balanced parentheses sequence


    bp_interval_p(Int ipos=0, Int cipos=0, Int jp1pos=0):ipos(ipos),cipos(cipos),jp1pos(jp1pos) {};


    bool operator<(const bp_interval_p& interval)const {
        if (ipos!=interval.ipos)
            return ipos<interval.ipos;
        return jp1pos<interval.jp1pos;
    }


    bool operator==(const bp_interval_p& interval)const {
        return ipos==interval.ipos and
               jp1pos==interval.jp1pos;
    }


    bool operator!=(const bp_interval_p& interval)const {
        return !(*this==interval);
    }


    bp_interval_p& operator=(const bp_interval_p& interval) {
        ipos = interval.ipos;
        cipos = interval.cipos;
        jp1pos = interval.jp1pos;
        return *this;
    }
};


template<class Int>
inline std::ostream& operator<<(std::ostream& os, const bp_interval_p<Int>& interval)
{
    os<<"-("<<interval.ipos<<","<<interval.cipos<<","<<interval.jp1pos<<")";
    return os;
}



template<class Csa = csa_wt<>, class Lcp = lcp_support_tree<lcp_wt<> >, class Bp_support = bp_support_sada<>, class Rank_support = rank_support_v5<> >
class cst_sct3p
{
    public:
        typedef typename Csa::value_type                                                         value_type;    // STL Container requirement/TODO: ist das nicht gleich node type???
        typedef cst_dfs_const_forward_iterator<cst_sct3p>                        const_iterator;// STL Container requirement
//      typedef const_iterator                                                                           iterator;              // STL Container requirement
        typedef cst_bottom_up_const_forward_iterator<cst_sct3p>          const_bottom_up_iterator;
//      typedef  const_bottom_up_iterator                                                        bottom_up_iterator;
        typedef const value_type                                                                         const_reference;
        typedef const_reference                                                                          reference;
        typedef const_reference*                                                                         pointer;
        typedef const pointer                                                                            const_pointer;
        typedef typename Csa::size_type                                                          size_type;             // STL Container requirement
        typedef size_type                                                                                        cst_size_type;
        typedef ptrdiff_t                                                                                        difference_type; // STL Container requirement
        typedef Csa                                                                                                      csa_type;
        typedef typename Lcp::template type<cst_sct3p>::lcp_type                 lcp_type;
        typedef Bp_support                                                                                       bp_support_type;
        typedef typename Csa::pattern_type                                                       pattern_type;
        typedef typename Csa::char_type                                                          char_type;
        typedef bp_interval_p<size_type>                                                         node_type; 
        typedef Rank_support                                                                             fc_rank_support_type;

        typedef cst_tag                                                                                          index_category;
    private:
        Csa                                     m_csa;
        lcp_type                                m_lcp;
        bit_vector                              m_bp;
        bp_support_type                 m_bp_support;
        bit_vector                              m_first_child; // implementation note: no rank structure is needed for the first_child bit_vector, except for id()
        fc_rank_support_type    m_first_child_rank;
        uint8_t                                 m_sigma;
        size_type                               m_nodes;

        void copy(const cst_sct3p& cst) {
            m_csa                       = cst.m_csa;
            copy_lcp(m_lcp, cst.m_lcp, *this);
            m_bp                        = cst.m_bp;
            m_bp_support        = cst.m_bp_support;
            m_bp_support.set_vector(&m_bp);
            m_first_child       = cst.m_first_child;
            m_first_child_rank = cst.m_first_child_rank;
            m_first_child_rank.set_vector(&m_first_child);
            m_sigma                     = cst.m_sigma;
            m_nodes                     = cst.m_nodes;
        }

        // Get the first l index of a [i,j] interval.
        /* I.e. given an interval [i,j], the function returns the position of the smallest entry lcp[k] with \f$ i<k\leq j \f$
         * \par Time complexity
         *       \f$ \Order{1} \f$
         */
        inline size_type get_first_l_index(const node_type& node, size_type& kpos, size_type& ckpos)const {
            if (node.cipos > node.jp1pos) { // corresponds to m_lcp[i] <= m_lcp[j+1]
                ckpos   = node.jp1pos-1;
            } else { // corresponds to m_lcp[i] > m_lcp[j+1]
                ckpos   = node.cipos-1;
            }
            assert(m_bp[ckpos]==0);
            kpos        = m_bp_support.find_open(ckpos);
            return m_bp_support.rank(kpos)-1;
        }

        // Get the ith l index of a node
        // if there exists no ith l-index return node.j+1
        size_type get_ith_l_index(const node_type& node, size_type i, size_type& kpos, size_type& ckpos)const {
            uint64_t children = 0;

            assert(i > 0);
            if (node.cipos > node.jp1pos) { // corresponds to m_lcp[i] <= m_lcp[j+1]
                ckpos   = node.jp1pos-1;
            } else { // corresponds to m_lcp[i] > m_lcp[j+1]
                ckpos   = node.cipos-1;
            }
            assert(m_bp[ckpos] == 0);   // at least the first l-index should be present, i.e. node is not leaf
            if (1 == i) {
                kpos    = m_bp_support.find_open(ckpos);
                return m_bp_support.rank(kpos)-1;
            } else { // i > 1
                size_type r = ckpos - m_bp_support.rank(ckpos); // numbers of closing parentheses - 1 = index of first child in m_first_child
                if (r+1 >= i) { // if there exist more than i l-indices
                    // check if m_first_child[r-i+1..r-1] consists of zeros
                    const uint64_t* p = m_first_child.data() + (r>>6);
                    uint8_t offset = r&0x3F;

                    uint64_t w = (*p) & bit_magic::Li1Mask[offset];
                    if (w) {
                        children = offset - bit_magic::l1BP(w) + 1;
                    } else if (m_first_child.data() == p) { // w==0 and we are in the first word
                        children = offset + 2; // da bit_magic::l1BP(w)=-1 sein muesste
                    } else {
                        children = offset + 2;
                        while (p > m_first_child.data()) {
                            w = *(--p);
                            if (0==w)
                                children += 64;
                            else {
                                children += (63-bit_magic::l1BP(w));
                                break;
                            }
                        }
                        children += (w==0);
                    }
                    if (i < children) {  // there exists an ith l-index
                        ckpos -= (i-1);
                        assert(m_bp[ckpos] == 0);
                        kpos   = m_bp_support.find_open(ckpos);
                        return m_bp_support.rank(kpos)-1;
                    }
                }
                // if i >= children
                kpos = node.jp1pos;
                if (kpos < m_bp.size())
                    ckpos = m_bp_support.find_close(kpos);
                else
                    ckpos = m_bp.size();
                return node.j+1;
            }
        }

        // Get the postion of the first l index of a l-[i,j] interval in the balanced parentheses sequence.
        /* \par Time complexity
         *              \f$ \Order{1} \f$
         */
        inline size_type get_pos_of_closing_parenthesis_of_the_first_l_index(const node_type& node)const {
            if (node.cipos > node.jp1pos) { // corresponds to m_lcp[i] <= m_lcp[j+1]
                return node.jp1pos-1;
            } else { // corresponds to m_lcp[i] > m_lcp[j+1]
                return node.cipos-1;
            }
        }

        // Get the next smaller value.
        /* \par Time complexity
         *      \f$ \Order{1} \f$
         */
        // Returns n if there is no next smaller value in [i+1..n-1]
        inline size_type nsv(size_type i, size_type ipos)const { // possible optimization: calculate also position of nsv, i.e. next ( following position cipos
            size_type cipos = m_bp_support.find_close(ipos);
            size_type result = m_bp_support.rank(cipos);
            return result;
        }

        // Get the previous smaller value.
        /*
         * \par Time complexity
         *    \f$ \Order{\frac{\sigma}{w}} \f$, where w=64 is the word size, can be implemented in \f$\Order{1}\f$ with rank and select
         */
        inline size_type psv(size_type i, size_type ipos, size_type cipos, size_type& psvpos, size_type& psvcpos)const {
            if (cipos >= m_bp.size() - m_sigma) { // if lcp[i]==0 => psv is the 0th index by definition
                psvpos = 0;
                psvcpos = m_bp.size()-1;
                return 0;
            }
            size_type r0 = cipos - m_bp_support.rank(cipos); // index of clothing parenthesis in first_child bp
            size_type next_first_child = 0;
            const uint64_t* p = m_first_child.data() + (r0>>6);
            uint64_t w = (*p) >> (r0&0x3F);
            if (w) { // if w!=0
                next_first_child = cipos + bit_magic::r1BP(w);
                if (cipos == next_first_child and m_bp[next_first_child+1]) {
                    psvpos = m_bp_support.enclose(ipos);
                    psvcpos = m_bp_support.find_close(psvpos);
                    return m_bp_support.rank(psvpos)-1;
                }
            } else {
                cipos += 64-(r0&0x3F);
                ++p;
                while (!(w=*p)) { // while w==0
                    ++p;
                    cipos += 64;
                }
                next_first_child = cipos + bit_magic::r1BP(w);
            }
            if (!m_bp[next_first_child+1]) { // if next parenthesis is a closing one
                psvcpos = next_first_child+1;
                psvpos = m_bp_support.find_open(psvcpos);
                return m_bp_support.rank(psvpos)-1;
            } else {
                psvpos = m_bp_support.enclose(m_bp_support.find_open(next_first_child));
                psvcpos = m_bp_support.find_close(psvpos);
                return m_bp_support.rank(psvpos)-1;
            }
        }

        // Range minimum query based on the rr_enclose method.
        /* \par Time complexity
         *     \f$ \Order{\rrenclose} \f$
         */
        inline size_type rmq(size_type l, size_type r)const {
            size_type i         = m_bp_support.select(l+1);
            size_type j         = m_bp_support.select(r+1);
            size_type fc_i      = m_bp_support.find_close(i);
            if (j < fc_i) { // i < j < find_close(j) < find_close(i)
                return l;
            } else { // i < find_close(i) < j < find_close(j)
                size_type ec = m_bp_support.rr_enclose(i,j);
                if (ec == m_bp_support.size()) {// no restricted enclosing pair found
                    return r;
                } else { // found range restriced enclosing pair
                    return m_bp_support.rank(ec)-1; // subtract 1, as the index is 0 based
                }
            }
        }

    public:
        const Csa& csa;                         
        const lcp_type& lcp;                            
        const bit_vector& bp;                   
        const bp_support_type& bp_support;      

        const bit_vector& first_child_bv;
        const fc_rank_support_type& first_child_rank;

        /* @{ */

        cst_sct3p(): csa(m_csa), lcp(m_lcp), bp(m_bp), bp_support(m_bp_support), first_child_bv(m_first_child), first_child_rank(m_first_child_rank) {}

        // Constructor for the cst_sct3p taking a string for that the compressed suffix tree should be calculated.
        /*
         * \param str Text for which the \f$ \CST \f$ should be constructed. The text should be terminated by a zero byte.
         * \pre The text has to be terminated by a zero byte.
         */
//              cst_sct3p(const unsigned char *str);

        template<uint8_t int_width, class size_type_class, uint8_t int_width_1, class size_type_class_1, uint8_t int_width_2, class size_type_class_2>
        cst_sct3p(const std::string& cst_file_name,
                  int_vector_file_buffer<int_width, size_type_class>& lcp_buf,
                  int_vector_file_buffer<int_width_1, size_type_class_1>& sa_buf,
                  int_vector_file_buffer<int_width_2, size_type_class_2>& isa_buf,
                  std::string dir,
                  bool build_ony_bps);

        cst_sct3p(tMSS& file_map, const std::string& dir, const std::string& id, bool build_only_bps);


        cst_sct3p(const cst_sct3p& cst):csa(m_csa),lcp(m_lcp),bp(m_bp),bp_support(m_bp_support),first_child_bv(m_first_child), first_child_rank(m_first_child_rank) {
            copy(cst);
        }

        /* @} */

        ~cst_sct3p() {}

        void construct(tMSS& file_map, const std::string& dir, const std::string& id, bool build_only_bps);


        size_type size()const {
            return m_bp.size()>>1;
        }


        static size_type max_size() {
            return Csa::max_size();
        }


        bool empty()const {
            return m_csa.empty();
        }


        void swap(cst_sct3p& cst);


        const_iterator begin()const {
            if (0 == m_bp.size())  // special case: tree is uninitialized
                return end();
            else if (2 == m_bp.size()) { // special case: the root is a leaf
                return const_iterator(this, root(), true, true);
            }
            return const_iterator(this, root(), false, true);
        };


        const_iterator end()const {
            return const_iterator(this, root(), true, false);
        }

        const_bottom_up_iterator begin_bottom_up()const {
            if (0 == m_bp.size())  // special case: tree is uninitialized
                return end_bottom_up();
            return const_bottom_up_iterator(this, leftmost_leaf_in_the_subtree(root()));
        }

        const_bottom_up_iterator end_bottom_up()const {
            return const_bottom_up_iterator(this, root(), false);
        }


        cst_sct3p& operator=(const cst_sct3p& cst);


        bool operator==(const cst_sct3p& cst)const;


        bool operator!=(const cst_sct3p& cst)const;


        size_type serialize(std::ostream& out, structure_tree_node* v=NULL, std::string name="")const;


        void load(std::istream& in);

        /* @{ */


        node_type root() const {
            return node_type(0, m_bp.size()-1, m_bp.size());
        }


        bool is_leaf(const node_type& v)const {
            return v.i==v.j;
        }


        node_type ith_leaf(size_type i)const {
            assert(i > 0 and i <= size());
            size_type ipos = m_bp_support.select(i);
            size_type jp1pos;
            if (i == size())
                jp1pos = m_bp.size();
            else if (m_bp[ipos+1])
                jp1pos = ipos+1;
            else
                jp1pos = m_bp_support.select(i+1);
            return node_type(ipos, m_bp_support.find_close(ipos), jp1pos);
        }


        size_type leaves_in_the_subtree(const node_type& v)const {
            return v.j-v.i+1;
        }


        node_type leftmost_leaf_in_the_subtree(const node_type& v)const {
            return ith_leaf(v.i+1);
        }


        node_type rightmost_leaf_in_the_subtree(const node_type& v)const {
            return ith_leaf(v.j+1);
        }


        size_type lb(const node_type& v)const {
            return v.i;
        }


        size_type rb(const node_type& v)const {
            return v.j;
        }
        /*
                        bool lcp_value_equals_zero(size_type i)const{
                                size_type ipos = m_bp_support.find_close(m_bp_support.select(i+1));
                                return ipos >= (m_bp.size() - m_sigma);
                        }
        */


        node_type parent(const node_type& v) const {
//std::cout<<"parent "<<depth(v)<<v<<std::endl;
#ifndef NDEBUG
//                      std::cout<<"parent interval of "<<v.l<<"-["<<v.left<<" "<<v.right<<"]"<<std::endl;
#endif
            if (v.cipos > v.jp1pos) { // LCP[i] <= LCP[j+1]
                size_type psv_pos, psv_cpos, psv_v, nsv_v, nsv_p1pos;
                psv_v = psv(v.j+1, v.jp1pos, m_bp_support.find_close(v.jp1pos), psv_pos, psv_cpos);
                nsv_v = nsv(v.j+1, v.jp1pos)-1;
                if (nsv_v == size()-1) {
                    nsv_p1pos = m_bp.size();
                } else { // nsv_v < size()-1
                    nsv_p1pos = m_bp_support.select(nsv_v+2);
                }
                return node_type(psv_v, nsv_v, psv_pos, psv_cpos, nsv_p1pos);
            } else { // LCP[i] > LCP[j+1]
                size_type psv_pos, psv_cpos, psv_v;
                psv_v = psv(v.i, v.ipos, v.cipos, psv_pos, psv_cpos);
                return node_type(psv_v, v.j, psv_pos, psv_cpos, v.jp1pos);
            }
        }


        node_type sibling(const node_type& v)const {
//std::cout<<"sibling "<<depth(v)<<v<<std::endl;
// Vorgehen:(1) Bestimmen, ob v noch einen rechten Bruder hat. Entsp. parent hat gleicht rechte Grenze wie v. Speziallfall rechter Rand ist hier schon behandelt!
            if (v.cipos < v.jp1pos) { // LCP[i] > LCP[j+1] => v hat die selbe rechte Grenze wie parent(v) => kein rechter Bruder
                return root();
            }
//          (2) Es existiert ein rechter Bruder, LCP[j+1] >= LCP[i] und j>i
            // es gilt nun: v.cipos > v.jp1pos
            size_type cjp1posm1 = m_bp_support.find_close(v.jp1pos)-1; // v.cipos-2 ???
//std::cerr<<"cjp1posm1="<<cjp1posm1<<std::endl;
//                      size_type cjp1posm1 = v.cipos-1;
//                      if( cjp1posm1+1 != v.cipos ){
//                              std::cout<<v<<" "<<cjp1posm1+1<<" "<<v.cipos<<std::endl;
//                      }
//                      size_type cjp1posm1 = v.cipos-2;
            // an der stelle von cjp1posm1 steht 1, falls v das letzte Kind ist
            bool last_child = m_bp[cjp1posm1];
            // an der stelle von cjp1posm1 steht 0 -> entweder nicht letztes Kind oder letztes Kind
            if (!last_child) {
                size_type first_child_idx = cjp1posm1 - m_bp_support.rank(cjp1posm1);
                last_child = m_first_child[first_child_idx]; // if first_child indicator is true => the new sibling is the rightmost sibling
            }
            if (last_child) {
                size_type nsv_v = nsv(v.j+1, v.jp1pos)-1, nsv_p1pos;
                if (nsv_v == size()-1) {
                    nsv_p1pos = m_bp.size();
                } else {
                    nsv_p1pos = m_bp_support.select(nsv_v+2);
                }
                return node_type(v.j+1, nsv_v, v.jp1pos, m_bp_support.find_close(v.jp1pos), nsv_p1pos);
            } else {
//std::cerr<<"not last child: "<<std::endl;
                /*
                                        if( cjp1posm1+1 != v.cipos-1 ){
                                                std::cout<<v<<" "<<cjp1posm1+1<<" "<<v.cipos<<std::endl;
                                                std::cout<<"LCP["<<v.i<<"]="<<m_lcp[v.i]<<" "<<"LCP["<<v.j<<"]="<<m_lcp[v.j];
                                                if(v.j<size())
                                                        std::cout<<" LCP[j+1]="<<m_lcp[v.j+1]<<std::endl;
                                        }
                */
//std::cerr<<"find_open(cjp1posm1) = "<< m_bp_support.find_open(cjp1posm1) << std::endl;
                size_type new_j = m_bp_support.rank(m_bp_support.find_open(cjp1posm1))-2;
                return node_type(v.j+1, new_j, v.jp1pos, m_bp_support.find_close(v.jp1pos), m_bp_support.select(new_j+2));
            }
        }


        node_type ith_child(const node_type& v, size_type i)const {
//                      std::cerr<<i<<"th child "<<depth(v)<<v<<" i="<<i<<std::endl;
            assert(i > 0);
            if (is_leaf(v))  // if v is a leave, v has no child
                return root();
            if (1 == i) {
                size_type k = 0, kpos = 0, k_find_close = 0;
                // v is not a leave: v has at least two children
                k = get_first_l_index(v, kpos, k_find_close);// get first l-index k and the position of k
                return node_type(v.i, k-1, v.ipos, v.cipos, kpos);
            } else { // i > 1
                size_type k1, kpos1, k_find_close1;
                k1 = get_ith_l_index(v, i-1, kpos1, k_find_close1);
                if (k1 == v.j+1)
                    return root();
                size_type k2, kpos2, k_find_close2;
                k2 = get_ith_l_index(v, i, kpos2, k_find_close2);
                return node_type(k1, k2-1, kpos1, k_find_close1, kpos2);
            }
        }


        size_type degree(const node_type& v)const {
            if (is_leaf(v))  // if v is a leave, v has no child
                return 0;
            // v is not a leave: v has at least two children
            size_type r = get_pos_of_closing_parenthesis_of_the_first_l_index(v);
            /*                  if( m_bp[r-1] ){// if there exists no next l-index
                                        return 2;
                                }
            */
            size_type r0 = r - m_bp_support.rank(r);
            const uint64_t* p = m_first_child.data() + (r0>>6);
            uint8_t offset = r0&0x3F;

            uint64_t w = (*p) & bit_magic::Li1Mask[offset];
            if (w) {
                return offset-bit_magic::l1BP(w)+1;
            } else if (m_first_child.data() == p) { // w==0 and we are in the first word
                return offset+2; // da bit_magic::l1BP(w)=-1 sein muesste
            } else {
                size_type res = offset+2;
                while (p > m_first_child.data()) {
                    w = *(--p);
                    if (0==w)
                        res += 64;
                    else {
                        return res + (63-bit_magic::l1BP(w));
                    }
                }
                return res + (w==0);
            }
        }



        // Returns the next sibling of node v.
        // Only for tests.
        node_type sibling_naive(const node_type& v)const {
            if (v==root())
                return root();
            node_type parent = this->parent(v);
            assert(parent != v);
            size_type nr = degree(parent);
            for (size_type i=1; i <= nr; ++i)
                if (ith_child(parent, i) == v and i!=nr)
                    return ith_child(parent, i+1);
            return root();
        }



        // TODO const unsigned char c durch char_type ersetzen
        node_type child(const node_type& v, const unsigned char c, size_type& char_pos)const {
            if (is_leaf(v))  // if v is a leaf = (), v has no child
                return root();
            // else v = ( (     ))
            uint16_t cc = m_csa.char2comp[c];
            if (cc==0 and c!=0) // TODO: aendere char2comp so ab, dass man diesen sonderfall nicht braucht
                return root();
            size_type char_ex_max_pos = m_csa.C[cc+1], char_inc_min_pos = m_csa.C[cc];

            size_type d                 = depth(v);

//                      (1) check the first child
//                      char_pos = m_csa(m_csa[v.i]+d);  // replaced by the next line
            char_pos = get_char_pos(v.i, d, m_csa);
            if (char_pos >= char_ex_max_pos) {// the first character of the first child interval is lex. greater than c
                // => all other first characters of the child intervals are also greater than c => no solution
                return root();
            } else if (char_pos >= char_inc_min_pos) { // i.e. char_pos < char_ex_max_pos and char_pos >= char_inc_min_pos
                return ith_child(v, 1);
            }

            size_type child_cnt         = degree(v);

//                      (2) check the last child
//                      char_pos = m_csa(m_csa[v.j]+d); // replaced by the next line
            char_pos = get_char_pos(v.j, d, m_csa);
            if (char_pos < char_inc_min_pos) {// the first character of the last child interval is lex. smaller than c
                // =>   all other first characters of the child intervals are also smaller than c => no solution
                return root();
            } else if (char_pos < char_ex_max_pos) { // i.e. char_pos < char_ex_max_pos and char_pos >= char_inc_min_pos
                return ith_child(v, child_cnt);
            }

//                      (3) binary search for c in the children [2..children)
            size_type l_bound = 2, r_bound = child_cnt, mid, kpos, ckpos, l_index;
            while (l_bound < r_bound) {
                mid = (l_bound + r_bound) >> 1;

                l_index = get_ith_l_index(v, mid-1, kpos, ckpos);
//                              char_pos = m_csa(m_csa[l_index]+d); // replaced by the next line
                char_pos = get_char_pos(l_index, d, m_csa);

                if (char_inc_min_pos > char_pos) {
                    l_bound = mid+1;
                } else if (char_ex_max_pos <= char_pos) {
                    r_bound = mid;
                } else { // char_inc_min_pos <= char_pos < char_ex_max_pos => found child
                    // we know that the child is not the last child, see (2)
                    // find next l_index: we know that a new l_index exists: i.e. assert( 0 == m_bp[ckpos-1]);
                    size_type lp1_index = m_bp_support.rank(m_bp_support.find_open(ckpos-1))-1;
                    size_type jp1pos = m_bp.size();
                    if (lp1_index-1 < size()-1) {
                        jp1pos = m_bp_support.select(lp1_index+1);
                    }
                    return node_type(l_index, lp1_index-1, kpos, ckpos, jp1pos);
                }
            }
            return root();
        }

        // \sa child(node_type v, const unsigned char c, size_type &char_pos)
        node_type child(const node_type& v, const unsigned char c) {
            size_type char_pos;
            return child(v, c, char_pos);
        }


        unsigned char edge(const node_type& v, size_type d)const {
#ifndef NDEBUG
            if (d < 1 or d > depth(v)) {
                throw std::out_of_range("OUT_OF_RANGE_ERROR: "+util::demangle(typeid(this).name())+"cst_sct3p<>::edge(node_type v, size_type d). d == 0 or d > depth(v)!");
            }
#endif
//                      size_type order = m_csa(m_csa[v.i]+d-1); replaced by the next line
            size_type order = get_char_pos(v.i, d-1, m_csa);
            uint16_t c_begin = 1, c_end = m_sigma+1, mid;
            while (c_begin < c_end) {
                mid = (c_begin+c_end)>>1;
                if (m_csa.C[mid] <= order) {
                    c_begin = mid+1;
                } else {
                    c_end = mid;
                }
            }
            return m_csa.comp2char[c_begin-1];
        }


        node_type lca(node_type v, node_type w)const {
            if (v.i > w.i or (v.i == w.i and v.j < w.j)) {
                std::swap(v, w);
            }
            assert(v.i < w.i or (v.i==w.i and v.j >= j));
            assert(!(v.i < w.i and v.j > w.i and v.j < w.j));   // assert that v and w do not overlapp
            if (v.j >= w.j) { // v encloses w or v==w
                return v;
            } else { // v.i < v.j < w.i < w.j
                size_type min_index = rmq(v.i+1, w.j);
                size_type min_index_pos         = m_bp_support.select(min_index+1);
                size_type min_index_cpos        = m_bp_support.find_close(min_index_pos);

                if (min_index_cpos >= m_bp.size() - m_sigma) { // if lcp[min_index]==0 => return root
                    return root();
                }
                size_type new_j = nsv(min_index, min_index_pos)-1;
                size_type new_ipos, new_icpos;
                size_type new_i = psv(min_index, min_index_pos, min_index_cpos, new_ipos, new_icpos);
                size_type jp1pos = m_bp.size();
                if (new_j < size()-1) {
                    jp1pos = m_bp_support.select(new_j+2);
                }
                return node_type(new_i, new_j, new_ipos, new_icpos, jp1pos);
            }
        }


        size_type depth(const node_type& v)const {
            if (v.i == v.j) {
                return size()-m_csa[v.i];
            } else if (v == root()) {
                return 0;
            } else {
                size_type kpos, ckpos;
                size_type l = get_first_l_index(v, kpos, ckpos);
                return m_lcp[l];
            }
        }


        // TODO: can be implemented in O(1) with o(n) space. See Jansson, Sadakane, Sung, SODA 2007, "Ultra-succinct Representation of Ordered Trees"
        size_type node_depth(node_type v)const {
            size_type d = 0;
            while (v != root()) {
                ++d;
                v = parent(v);
            }
            return d;
        }


        node_type sl(const node_type& v)const {
            if (v == root())
                return root();
            // get interval with first char deleted
            size_type i  = m_csa.psi[v.i];
            if (is_leaf(v)) {
                if (v.i==0 and v.j==0) // if( v.l==1 )
                    return root();
                else
                    return ith_leaf(i+1);
            }
            size_type j  = m_csa.psi[v.j];
            assert(i < j);
            size_type min_index = rmq(i+1, j); // rmq
            size_type min_index_pos     = m_bp_support.select(min_index+1);
            size_type min_index_cpos    = m_bp_support.find_close(min_index_pos);
            if (min_index_cpos >= m_bp.size() - m_sigma) { // if lcp[min_index]==0 => return root
                return root();
            }
            size_type new_j = nsv(min_index, min_index_pos)-1;
            size_type new_ipos, new_icpos;
            size_type new_i = psv(min_index, min_index_pos, min_index_cpos, new_ipos, new_icpos);
            size_type jp1pos = m_bp.size();
            if (new_j < size()-1) {
                jp1pos = m_bp_support.select(new_j+2);
            }
            return node_type(new_i, new_j, new_ipos, new_icpos, jp1pos);
        }


        /*
         *  \param v A valid not of a cst_sct3p.
         *  \param c The character which should be prepended to the string of the current node.
         *      \return root() if the Weiner link of (v, c) does not exist, otherwise the Weiner link is returned.
         *  \par Time complexity
         *              \f$ \Order{ t_{rank\_bwt} } \f$
         *
         */
        node_type wl(const node_type& v, const unsigned char c) const {
            size_type c_left    = m_csa.rank_bwt(v.i, c);
            size_type c_right   = m_csa.rank_bwt(v.j+1, c);
            if (c_left == c_right)  // there exists no Weiner link
                return root();
            if (c_left+1 == c_right)
                return ith_leaf(m_csa.C[m_csa.char2comp[c]] + c_left + 1);
            else {
                size_type left  = m_csa.char2comp[c] + c_left;
                size_type right = m_csa.char2comp[c] + c_right - 1;
                assert(left < right);

                size_type ipos = m_bp_support.select(left+1);
                size_type jp1pos = m_bp.size();
                if (right < size()-1) {
                    jp1pos = m_bp_support.select(right+2);
                }
                return node_type(left, right, ipos, m_bp_support.find_close(ipos), jp1pos);
            }
        }


        size_type sn(const node_type& v)const {
            assert(is_leaf(v));
            return m_csa[v.i];
        }


        size_type id(const node_type& v)const {
            if (is_leaf(v)) { // return i in the range from 0..csa.size()-1
                return v.i;
            }
            size_type clpos; // closing parentheses of the l-index
            if (v.cipos > v.jp1pos) { // corresponds to m_lcp[i] <= m_lcp[j+1]
                clpos   = v.jp1pos-1;
            } else { // corresponds to m_lcp[i] > m_lcp[j+1]
                clpos   = v.cipos-1;
            }
            assert(m_bp[clpos]==0);
            size_type r0clpos = clpos-m_bp_support.rank(clpos);

            return size()+m_first_child_rank(r0clpos);
        }

        size_type nodes()const {
            return m_nodes;
        }

        /* \param lb Left bound of the lcp-interval [lb..rb] (inclusive).
         * \param rb Right bound of the lcp-interval [lb..rb] (inclusive).
         * \param l  Depth of the lcp-interval.
         *\ return The node in the suffix tree corresponding lcp-interval [lb..rb]
         */
        node_type node(size_type lb, size_type rb, size_type l=0) const {
            size_type ipos = m_bp_support.select(lb+1);
            size_type jp1pos;
            if (rb == size()-1) {
                jp1pos = m_bp.size();
            } else {
                jp1pos = m_bp_support.select(rb+2);
            }
            return node_type(ipos, m_bp_support.find_close(ipos), jp1pos);
        }


        size_type tlcp_idx(size_type i) const {
            size_type ipos      = m_bp_support.select(i+1);
            size_type cipos = m_bp_support.find_close(ipos);
            return m_first_child_rank.rank(((ipos+cipos-1)>>1)-i);
        }
        /* @} */

};

// == template functions ==


template<class Csa, class Lcp, class Bp_support, class Rank_support>
template<uint8_t int_width, class size_type_class, uint8_t int_width_1, class size_type_class_1, uint8_t int_width_2, class size_type_class_2>
cst_sct3p<Csa, Lcp, Bp_support, Rank_support>::cst_sct3p(const std::string& csa_file_name,
        int_vector_file_buffer<int_width, size_type_class>& lcp_buf,
        int_vector_file_buffer<int_width_1, size_type_class_1>& sa_buf,
        int_vector_file_buffer<int_width_2, size_type_class_2>& isa_buf,
        std::string dir="./",
        bool build_only_bps=false
                                                        ):csa(m_csa), lcp(m_lcp), bp(m_bp), bp_support(m_bp_support), first_child_bv(m_first_child), first_child_rank(m_first_child_rank)
{
    std::string id =  util::to_string(util::get_pid())+"_"+util::to_string(util::get_id()).c_str();

    write_R_output("cst", "construct BPS", "begin", 1, 0);
    m_nodes = algorithm::construct_supercartesian_tree_bp_succinct_and_first_child(lcp_buf, m_bp, m_first_child) + m_bp.size()/2;
    write_R_output("cst", "construct BPS", "end", 1, 0);

    if (!build_only_bps) {
        util::load_from_file(m_csa, csa_file_name.c_str());

        write_R_output("cst", "construct CLCP", "begin", 1, 0);
        construct_lcp(m_lcp, *this, lcp_buf, isa_buf);
        write_R_output("cst", "construct CLCP", "end", 1, 0);
    }
    write_R_output("cst", "construct BPSS", "begin", 1, 0);
    m_bp_support = Bp_support(&m_bp);
    util::init_support(m_first_child_rank, &m_first_child);
    write_R_output("cst", "construct BPSS", "end", 1, 0);
    m_sigma = degree(root());
}


template<class Csa, class Lcp, class Bp_support, class Rank_support>
cst_sct3p<Csa, Lcp, Bp_support, Rank_support>::cst_sct3p(tMSS& file_map, const std::string& dir, const std::string& id, bool build_only_bps = false):csa(m_csa), lcp(m_lcp), bp(m_bp), bp_support(m_bp_support), first_child_bv(m_first_child), first_child_rank(m_first_child_rank)
{
    construct(file_map, dir, id, build_only_bps);
}


template<class Csa, class Lcp, class Bp_support, class Rank_support>
void cst_sct3p<Csa, Lcp, Bp_support, Rank_support>::construct(tMSS& file_map, const std::string& dir, const std::string& id, bool build_only_bps = false)
{
    write_R_output("cst", "construct BPS", "begin", 1, 0);
    int_vector_file_buffer<> lcp_buf(file_map["lcp"].c_str());
    m_nodes = algorithm::construct_supercartesian_tree_bp_succinct_and_first_child(lcp_buf, m_bp, m_first_child) + m_bp.size()/2;
    write_R_output("cst", "construct BPS", "end", 1, 0);
    write_R_output("cst", "construct BPSS", "begin", 1, 0);
    m_bp_support = Bp_support(&m_bp);
    util::init_support(m_first_child_rank, &m_first_child);
    write_R_output("cst", "construct BPSS", "end", 1, 0);

    if (!build_only_bps) {
        write_R_output("cst", "construct CLCP", "begin", 1, 0);
        construct_lcp(m_lcp, *this, file_map, dir, id);
        write_R_output("cst", "construct CLCP", "end", 1, 0);
    }
    if (!build_only_bps) {
        util::load_from_file(m_csa, file_map["csa"].c_str());
    }
    m_sigma = degree(root());
}

template<class Csa, class Lcp, class Bp_support, class Rank_support>
typename cst_sct3p<Csa, Lcp, Bp_support, Rank_support>::size_type cst_sct3p<Csa, Lcp, Bp_support, Rank_support>::serialize(std::ostream& out, structure_tree_node* v, 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 += m_csa.serialize(out, child, "csa");
    written_bytes += m_lcp.serialize(out, child, "lcp");
    written_bytes += m_bp.serialize(out, child, "bp");
    written_bytes += m_bp_support.serialize(out, child, "bp_support");
    written_bytes += m_first_child.serialize(out, child, "mark_child");
    written_bytes += m_first_child_rank.serialize(out, child, "mark_child_rank");
    written_bytes += util::write_member(m_sigma, out, child, "sigma");
    written_bytes += util::write_member(m_nodes, out, child, "node_cnt");
    structure_tree::add_size(child, written_bytes);
    return written_bytes;
}

template<class Csa, class Lcp, class Bp_support, class Rank_support>
void cst_sct3p<Csa, Lcp, Bp_support, Rank_support>::load(std::istream& in)
{
    m_csa.load(in);
    load_lcp(m_lcp, in, *this);
    m_bp.load(in);
    m_bp_support.load(in, &m_bp);
    m_first_child.load(in);
    m_first_child_rank.load(in, &m_first_child);
    util::read_member(m_sigma, in);
    util::read_member(m_nodes, in);
#ifdef SDSL_DEBUG
    assert(algorithm::check_bp_support(m_bp, m_bp_support));
    std::cerr<<"checked bp_support"<<std::endl;
#endif
}

template<class Csa, class Lcp, class Bp_support, class Rank_support>
cst_sct3p<Csa, Lcp, Bp_support, Rank_support>& cst_sct3p<Csa, Lcp, Bp_support, Rank_support>::operator=(const cst_sct3p& cst)
{
    if (this != &cst) {
        copy(cst);
    }
    return *this;
}





} // end namespace sdsl


#endif