gdst_ukkonen_inserter.hpp

Go to the documentation of this file.

#ifndef GDST_UKKONEN_INSERTER_HPP_
#define GDST_UKKONEN_INSERTER_HPP_
#include "boost/assert.hpp"
#include "vdj_pipe/detail/disable_warnings.hpp"
#include "vdj_pipe/gdst/gdst.hpp"

namespace vdj_pipe{ namespace gdst{

class Ukkonen_inserter {
public:
   static void insert(Gdst& st, const Seq_id sid) {
      for( Ukkonen_inserter ui(st, sid); ui.has_next(); ui.next() ) ui.insert();
   }

   Ukkonen_inserter(Gdst& st, const Seq_id sid)
   : st_(st),
     sid_(sid),
     suff_(st_[sid].sequence()),
     active_node_(st_.root()),
     active_length_(st_[active_node_].n_),
     md_(0),
     prev_branch_(),
     prev_leaf_()
   {}

   void next() {suff_.remove_prefix(1);}
   bool has_next() const {return ! suff_.empty();}

   void insert() {
      for( ; ; ) {
         if( prev_branch_ && st_[prev_branch_].n_ == active_length_ + 1 ) {
            set_branch_sl(active_node_);
            prev_branch_ = Branch_id();
         }

         if( ! md_ && active_length_ == suff_.size() ) {
               add_to_leaf(active_node_);
               if( prev_leaf_ ) set_leaf_sl(active_node_);
               if( ! st_[active_node_].sl_ ) prev_leaf_ = active_node_;
               else prev_leaf_ = Branch_id();
               break;
         }

         BOOST_ASSERT(active_length_ < suff_.size());
         if( ! md_ ) {
            if( active_child() ) {
               md_ = 1;
               continue;
            }
            const Branch_id leaf = leaf_from_branch(active_node_, active_edge());
            if( prev_leaf_ ) set_leaf_sl(leaf);
            prev_leaf_ = leaf;
            prev_branch_ = Branch_id();
            break;
         }

         const unsigned l = edge_length();
         if( l <= md_ ) {
            md_ -= l;
            set_active_node(active_child());
            continue;
         }

         const Nucleotide n2 = letter(active_length_ + md_);
         if( active_length_ + md_ == suff_.size() ) {
            const Branch_id leaf = split_active_edge(active_length_ + md_, n2);
            st_.bm_[leaf].sb_ = leaf;
            add_to_leaf(leaf);
            if( prev_leaf_ ) set_leaf_sl(leaf);
            prev_leaf_ = leaf;
            prev_branch_ = Branch_id();
            break;
         }

         const Nucleotide n1 = nucleotide_index(suff_[active_length_ + md_]);
         if( n1 == n2 ) {
            ++md_;
            continue;
         }

         const Branch_id branch = split_active_edge(active_length_ + md_, n2);
         const Branch_id leaf = leaf_from_branch(branch, n1);
         st_.bm_[branch].sb_ = leaf;
         if( prev_branch_ ) set_branch_sl(branch);
         if( prev_leaf_ ) set_leaf_sl(leaf);
         if( st_[branch].n_ == 1) {
            st_.bm_[branch].sl_ = st_.root();
            prev_branch_ = Branch_id();
         }
         else prev_branch_ = branch;
         if( st_[leaf].n_ == 1) {
            st_.bm_[leaf].sl_ = st_.root();
            prev_leaf_ = Branch_id();
         }
         else prev_leaf_ = leaf;
         break;
      }

      if( active_node_ == st_.root() ) {
         if( md_ ) --md_;
      } else {
         if( const Branch_id bid = st_[active_node_].sl_ ) set_active_node(bid);
         else set_active_node(st_.root());
      }
   }

private:
   Gdst& st_;
   const Seq_id sid_;
   boost::string_ref suff_;
   Branch_id active_node_;
   unsigned active_length_;
   unsigned md_;
   Branch_id prev_branch_, prev_leaf_;

   Nucleotide active_edge() const {
      BOOST_ASSERT(active_length_ < suff_.size());
      return nucleotide_index(suff_[active_length_]);
   }

   Branch_id active_child() const {
      return st_.child(active_node_, active_edge());
   }

   void active_child(const Branch_id bid) {
      BOOST_ASSERT(active_length_ < suff_.size());
      const Nucleotide n = nucleotide_index(suff_[active_length_]);
      st_.child(active_node_, n, bid);
   }

   void set_branch_sl(const Branch_id bid) {
      BOOST_ASSERT(bid);
      BOOST_ASSERT(prev_branch_);
      BOOST_ASSERT( ! st_[prev_branch_].sl_);
      BOOST_ASSERT( st_[prev_branch_].n_ == st_[bid].n_ + 1);
      st_.bm_[prev_branch_].sl_ = bid;
   }

   void set_leaf_sl(const Branch_id bid) {
      BOOST_ASSERT(bid);
      BOOST_ASSERT( ! st_[prev_leaf_].sl_);
      BOOST_ASSERT( st_[prev_leaf_].n_ == st_[bid].n_ + 1);
      st_.bm_[prev_leaf_].sl_ = bid;
   }

   void add_to_leaf(const Branch_id bid) {
      Branch& b = st_.bm_[bid];
      if( ! b.leaf_ ) {
         b.leaf_ = st_.lm_.insert(Leaf());
      }
      st_.lm_[b.leaf_].insert(sid_);
   }

   Branch_id leaf_from_branch(
            const Branch_id bid,
            const Nucleotide n1
   ) {
      const Branch_id bid2 = st_.bm_.insert(Branch(suff_.size()));
      add_to_leaf(bid2);
      st_.bm_[bid2].sb_ = bid2;  //its own suffix branch
      st_.child(bid, n1, bid2);
      return bid2;
   }

   void set_active_node(const Branch_id bid) {
      active_node_ = bid;
      active_length_ = st_[active_node_].n_;
   }

   Branch_id split_active_edge(
            const unsigned i,
            const Nucleotide n2
   ) {
      BOOST_ASSERT(active_length_ < i && "new child index should be greater");
      const Branch_id bid2 = st_.bm_.insert(Branch(i));
      const Branch_id bid3 = active_child();
      BOOST_ASSERT(i < st_[bid3].n_ && "old child index should be greater");
      active_child(bid2);
      st_.child(bid2, n2, bid3);
      return bid2;
   }

   unsigned edge_length() const {
      const Branch_id c = active_child();
      BOOST_ASSERT(c);
      BOOST_ASSERT(st_[active_node_].n_ < st_[c].n_);
      return st_[c].n_ - st_[active_node_].n_;
   }

   Nucleotide letter(const unsigned i) const {
      BOOST_ASSERT(active_node_);
      BOOST_ASSERT(active_length_ < i);
      const Branch_id c = active_child();
      BOOST_ASSERT(c);
      BOOST_ASSERT(i < st_[c].n_);
      return st_.letter(c, i);
   }
};

}//namespace gdst
}//namespace vdj_pipe
#endif /* GDST_UKKONEN_INSERTER_HPP_ */