exact_matcher.hh

#ifndef EXACT_MATCHER_HH
#define EXACT_MATCHER_HH

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

#include <iostream>
#include <sstream>
#include <list>
#include <algorithm>
#include <limits>
#include <iterator>
#include "trace_controller.hh"
#include "sparsification_mapper.hh"
#include "tuples.hh"
#include "scoring.hh"
#include "ext_rna_data.hh"
#include "aux.hh"

extern "C" {
#   include <ViennaRNA/fold_vars.h>
#   include <ViennaRNA/utils.h>
#   include <ViennaRNA/PS_dot.h>
#   include <ViennaRNA/fold.h>
    int    PS_rna_plot(char *string, char *structure, char *file);
    int    PS_rna_plot_a(char *string, char *structure, char *file, char *pre, char *post);
    float  fold(const char *sequence, char *structure);
}

namespace LocARNA {

    typedef size_t                                      size_type;
    typedef std::vector<unsigned int>                   intVec;
    typedef std::pair<unsigned int,unsigned int>        intPair;
    typedef std::pair<intPair, intPair>                         intPPair;
    typedef const intPPair*                     intPPairPTR;
    typedef std::vector<intPPair>::const_iterator       IntPPairCITER;


    class SinglePattern
    {
    public:
        SinglePattern(){};

        SinglePattern(const std::string& myId_,const std::string& seqId_,const intVec& mySinglePattern_)
            :myId(myId_),seqId(seqId_),pattern(mySinglePattern_)
        {};

        virtual ~SinglePattern() { pattern.clear(); };

        const std::string&        getmyId()  const { return myId; };

        const std::string&      getseqId() const {return seqId; };

        const intVec&        getPat() const { return pattern; };

    private:

        std::string         myId; 
        std::string          seqId; 
        intVec         pattern; 
    };

    class PatternPair
    {
    public:
        PatternPair(){};

        PatternPair(const std::string& myId,const SinglePattern& myFirstPat,const SinglePattern& mySecPat, const std::string& structure_, int& score_)
            : id(myId),first(myFirstPat),second(mySecPat), structure(structure_), EPMscore(score_)
        {
            if (first.getPat().size() != second.getPat().size()){
                std::cerr << "Error! PatternPair cannot be constructed due to different sizes of SinglePatterns!" << std::endl;
            }
            score = EPMscore;
            size = first.getPat().size();
        };

        virtual ~PatternPair()
        {
            insideBounds.clear();
        };

        const std::string&              getId() const { return id; };

        const int&              getSize() const { return size; };

        const SinglePattern&    getFirstPat() const { return first; };

        const SinglePattern&    getSecPat() const { return second;};

        void            resetBounds();

        void            setOutsideBounds(intPPair myPPair);

        const   intPPair                getOutsideBounds() const { return outsideBounds; };

        void            addInsideBounds(intPPair myPPair);

        const   std::vector<intPPair>& getInsideBounds() const { return insideBounds; };

        void                    setEPMScore(int myScore);

        const           int                     getScore() const { return score;  };

        const           int                     getEPMScore() const { return EPMscore; };

        const std::string& get_struct() const{return structure;};

    private:
        std::string             id; 
        int             size; 
        SinglePattern   first; 
        SinglePattern   second; 

        std::string             structure; 
        int                             score; 
        int                             EPMscore; 
        std::vector<intPPair>   insideBounds; 
        intPPair           outsideBounds; 
    };

    class PatternPairMap
    {
    public:
        typedef  PatternPair                                  selfValueTYPE; 
        typedef  PatternPair*                                                   SelfValuePTR; 

        typedef  std::multimap<int,SelfValuePTR,std::greater<int> >     orderedMapTYPE; 
        typedef  orderedMapTYPE::const_iterator               orderedMapCITER; 
        typedef  orderedMapTYPE::iterator                     orderedMapITER; 
        typedef  std::list<SelfValuePTR>                           patListTYPE; 
        typedef  patListTYPE::iterator                        patListITER; 
        typedef  patListTYPE::const_iterator                  patListCITER; 
        typedef unordered_map<std::string,SelfValuePTR>::type PatternIdMapTYPE; 


        PatternPairMap();

        PatternPairMap(const PatternPairMap& myPairMap)
            :patternList(myPairMap.patternList),
             patternOrderedMap(myPairMap.patternOrderedMap),
             idMap(myPairMap.idMap)  { minPatternSize = 100000;};

        virtual ~PatternPairMap();

        void              add( const std::string& id,
                               const SinglePattern& first,
                               const SinglePattern& second,
                               const std::string& structure,
                               int score
                               );

        void              add(const SelfValuePTR value);

        void              makeOrderedMap();

        void              updateFromMap();

        const    PatternPair&      getPatternPair(const std::string& id)const;

        const    SelfValuePTR      getPatternPairPTR(const std::string& id)const;

        const    patListTYPE&      getList() const;

        const    orderedMapTYPE&   getOrderedMap() const;

        orderedMapTYPE&   getOrderedMap2();

        const    int               size()   const;

        int              getMapBases();

        int              getMapEPMScore();

        const    int             getMinPatternSize() const { return minPatternSize; };

    private:

        patListTYPE        patternList; 
        orderedMapTYPE     patternOrderedMap; 
        PatternIdMapTYPE   idMap; 
        int minPatternSize; 
    };

    std::ostream &operator << (std::ostream &out, const PatternPairMap::patListTYPE &pat_pair_map);

    class LCSEPM
    {
    public:

        LCSEPM(const    Sequence&               seqA_,
               const    Sequence&               seqB_,
               const    PatternPairMap& myPatterns,
               PatternPairMap& myLCSEPM
               )

            :seqA(seqA_),
             seqB(seqB_),
             matchedEPMs(myLCSEPM),
             patterns(myPatterns)
        {};

        virtual         ~LCSEPM();

        void            MapToPS(const std::string& sequenceA, const std::string& sequenceB, PatternPairMap& myMap, const std::string& file1, const std::string& file2);

        void            calculateLCSEPM(bool quiet);
    
        std::pair<SequenceAnnotation,SequenceAnnotation>
        anchor_annotation();
    
        void            output_locarna(const std::string& sequenceA, const std::string& sequenceB, const std::string& outfile);

        void            output_clustal(const std::string& outfile_name);

    private:

        struct HoleCompare2 {
            bool operator()(const intPPairPTR & h1, const intPPairPTR & h2) const {
                // first compare size of holes
                if (h1->first.second - h1->first.first-1 < h2->first.second - h2->first.first-1){
                    return true; }
                // compare if holes are identical in both structures
                if (h1->first.second - h1->first.first-1 == h2->first.second - h2->first.first-1){
                    if ((h1->first.first == h2->first.first) && (h1->first.second == h2->first.second) &&
                        (h1->second.first==h2->second.first) && (h1->second.second==h2->second.second))
                        { return true; }
                }

                return false;
            }
        };

        typedef     std::multimap<intPPairPTR,PatternPairMap::SelfValuePTR,HoleCompare2>        HoleOrderingMapTYPE2;
        typedef     HoleOrderingMapTYPE2::const_iterator HoleMapCITER2;


        void    preProcessing                   ();
        void    calculateHoles3                 (bool quiet);
        void    calculatePatternBoundaries      (PatternPair* myPair);
        void    calculateTraceback2             (const int i,const int j,const int k,const int l,std::vector < std::vector<int> > holeVec);
        int     D_rec2                          (const int& i,const  int& j,const int& k,const int& l,std::vector < std::vector<int> >& D_h,const bool debug);

        int     max3                            (int a, int b, int c)
        {
            int tmp = a>b? a:b;
            return (tmp>c? tmp:c);
        };

        char* getStructure(PatternPairMap& myMap, bool firstSeq, int length);

        std::string intvec2str(const std::vector<unsigned int>& V, const std::string &delim){
            std::stringstream oss;
            copy(V.begin(), V.end(), std::ostream_iterator<unsigned int>(oss, delim.c_str()));
            std::string tmpstr;
            tmpstr = oss.str();
            if (tmpstr.length()>0) tmpstr.erase(tmpstr.end()-1);
            return tmpstr;
        }

        std::string upperCase(const std::string &seq){
            std::string s= "";
            for(unsigned int i= 0; i<seq.length(); i++)
                s+= toupper(seq[i]);
            return s;
        }

        std::vector< std::vector <std::vector<PatternPairMap::SelfValuePTR> > > EPM_Table2;
        HoleOrderingMapTYPE2                    holeOrdering2;
        const   Sequence&                               seqA;
        const   Sequence&                               seqB;
        PatternPairMap&                         matchedEPMs;
        const   PatternPairMap&                         patterns;
    };


    class SparseTraceController: public TraceController{

    private:

        typedef SparsificationMapper::matidx_t matidx_t; 
        typedef SparsificationMapper::seq_pos_t seqpos_t; 
        typedef SparsificationMapper::index_t index_t; 

    public:

        typedef std::pair<matidx_t,matidx_t> matpos_t; 
        typedef std::pair<seqpos_t,seqpos_t> pair_seqpos_t; 

    private:
        const SparsificationMapper &sparse_mapperA; 
        const SparsificationMapper &sparse_mapperB; 

    public:
        SparseTraceController(const SparsificationMapper &sparse_mapperA_,const SparsificationMapper &sparse_mapperB_,const TraceController &trace_controller_):
            TraceController::TraceController(trace_controller_),
            sparse_mapperA(sparse_mapperA_),
            sparse_mapperB(sparse_mapperB_)

        {

        }

        virtual ~SparseTraceController(){}; 

        const SparsificationMapper& get_sparse_mapperA() const{
            return sparse_mapperA;
        }

        const SparsificationMapper& get_sparse_mapperB() const{
            return sparse_mapperB;
        }


        matidx_t min_col_idx(index_t indexA, index_t indexB, matidx_t idx_i,
                             index_t left_endB =  std::numeric_limits<index_t>::max()) const{

            seqpos_t i = sparse_mapperA.get_pos_in_seq_new(indexA,idx_i);
            return sparse_mapperB.idx_geq(indexB,min_col(i),left_endB);
        }

        matidx_t idx_after_max_col_idx(index_t indexA, index_t indexB, matidx_t idx_i,
                                       index_t left_endB =  std::numeric_limits<index_t>::max()) const{

            seqpos_t i = sparse_mapperA.get_pos_in_seq_new(indexA,idx_i);
            return sparse_mapperB.idx_after_leq(indexB,max_col(i),left_endB);
        }

        matpos_t diag_pos_bef(index_t indexA, index_t indexB,pair_seqpos_t cur_pos_seq,
                              index_t left_endA =  std::numeric_limits<index_t>::max(), index_t left_endB =  std::numeric_limits<index_t>::max())const{

            bool debug_valid_mat_pos = false;

            if(debug_valid_mat_pos) std::cout << "first valid mat pos before with tc " << std::endl;

            seqpos_t i = cur_pos_seq.first;
            seqpos_t j = cur_pos_seq.second;
            
            matidx_t idx_after_max_col;

            // find valid matrix position based on the SparsificationMapper
            matidx_t cur_row = sparse_mapperA.first_valid_mat_pos_before(indexA,i,left_endA);
            matidx_t col_before = sparse_mapperB.first_valid_mat_pos_before(indexB,j,left_endB);

            //bool valid_pos_found = false;

            // find a valid position that is valid also based on the TraceController
            // go through the rows and find an interval that includes the column col_before or lies
            // before the column col_before
            for(;;--cur_row){

                matidx_t min_col = min_col_idx(indexA,indexB,cur_row,left_endB);
                idx_after_max_col = idx_after_max_col_idx(indexA,indexB,cur_row,left_endB);

                if(debug_valid_mat_pos) std::cout << "interval " << min_col << "," << idx_after_max_col << std::endl;

                // valid interval found
                if(min_col<idx_after_max_col && min_col<=col_before){
                    //valid_pos_found=true;
                    break;
                }

                if(cur_row==0){
                    break;
                }

            }

            //assert(valid_pos_found);
            assert(idx_after_max_col>0);

            matidx_t max_col = idx_after_max_col-1;

            // the column of the new position is the col_before or lies before it
            matpos_t result = matpos_t(cur_row,std::min(max_col,col_before));

            assert(is_valid_idx_pos(indexA,indexB,result));

            return result;

        }

        pair_seqpos_t pos_in_seq(index_t idxA, index_t idxB,//const Arc &a, const Arc &b,
                                 const matpos_t &cur_pos) const{

            return pair_seqpos_t(sparse_mapperA.get_pos_in_seq_new(idxA,cur_pos.first),
                                 sparse_mapperB.get_pos_in_seq_new(idxB,cur_pos.second));
        }

        bool matching_wo_gap(index_t idxA, index_t idxB,
                             matpos_t idx_pos_diag, pair_seqpos_t seq_pos_to_be_matched) const{
            pair_seqpos_t pos_diag = pos_in_seq(idxA,idxB,idx_pos_diag);
            return (pos_diag.first+1 == seq_pos_to_be_matched.first) && (pos_diag.second+1==seq_pos_to_be_matched.second);
        }

        bool pos_unpaired(index_t idxA, index_t idxB,
                          matpos_t pos) const{
            return sparse_mapperA.pos_unpaired(idxA,pos.first)
                && sparse_mapperB.pos_unpaired(idxB,pos.second);
        }

        bool is_valid_idx_pos(index_t idxA, index_t idxB,
                              matpos_t mat_pos) const{
            pair_seqpos_t seq_pos = pos_in_seq(idxA,idxB,mat_pos);
            return is_valid(seq_pos.first,seq_pos.second);
        }

    };


    class EPM{

    public:
        typedef BasePairs__Arc Arc; 

        typedef SparsificationMapper::seq_pos_t seqpos_t; 
        typedef SparseTraceController::matpos_t matpos_t; 
        typedef SparsificationMapper::ArcIdx ArcIdx; 
        typedef SparseTraceController::pair_seqpos_t pair_seqpos_t; 
        typedef std::pair<ArcIdx,ArcIdx> PairArcIdx; 
        typedef std::vector<PairArcIdx> PairArcIdxVec; 

        typedef triple<seqpos_t,seqpos_t,char> el_pat_vec;

        typedef std::vector<el_pat_vec> pat_vec_t; 

    private:

        pat_vec_t pat_vec; 

        score_t score; 
        int state; 
        matpos_t cur_pos; 
        score_t max_tol_left; 
        bool first_insertion; 
        bool invalid; 

        PairArcIdxVec am_to_do;

        struct compare_el_pat_vec {
        public:
            bool
            operator () (const EPM::el_pat_vec &el1,const EPM::el_pat_vec &el2)const {
                seqpos_t el1_pos1 = el1.first;
                seqpos_t el1_pos2 = el1.second;
                seqpos_t el2_pos1 = el2.first;
                seqpos_t el2_pos2 = el2.second;
                char el1_struc = el1.third;
                char el2_struc = el2.third;
                return (el1_pos1 < el2_pos1) || (el1_pos1==el2_pos1 && el1_pos2 < el2_pos2)
                    || (el1_pos1==el2_pos1 && el1_pos2 == el2_pos2 && el1_struc<el2_struc);
            }

        };

        struct compare_el_am_to_do {
        public:
            bool
            operator () (const EPM::PairArcIdx &el1, const EPM::PairArcIdx &el2)const {
                return (el1.first<el2.first) || ((el1.first==el2.first) && el1.second < el2.second);
            }
        };

    public:

        EPM():
            score(0),
            state(0),
            cur_pos(matpos_t(0,0)),
            max_tol_left(0),
            first_insertion(true),
            invalid(false)
        {}

        virtual ~EPM(){} 

        //-----------------------------------------------------------------------
        // getter methods
        //-----------------------------------------------------------------------

        score_t get_score() const{return score;}

        int get_state() const{return state;}

        const matpos_t & get_cur_pos() const{return cur_pos;}

        const score_t & get_max_tol_left() const{return max_tol_left;}

        bool get_first_insertion() const{return first_insertion;}

        bool is_invalid() const{
            return invalid;
        }

        //-----------------------------------------------------------------------
        // setter methods
        //-----------------------------------------------------------------------

        void set_score(score_t score_){ score=score_;}

        void set_state(int state_){     state=state_;}

        void set_cur_pos(matpos_t cur_pos_){cur_pos = cur_pos_;}

        void set_max_tol_left(score_t tol){max_tol_left=tol;}

        void set_first_insertion(bool first_insertion_){first_insertion=first_insertion_;}

        void set_invalid(){
            invalid = true;
        }

        const PairArcIdx& get_am(PairArcIdxVec::size_type idx) const{
            assert(idx<am_to_do.size());
            return am_to_do[idx];
        }

        PairArcIdxVec::size_type number_of_am(){return am_to_do.size();}

        void clear_am_to_do(){am_to_do.clear();}

        PairArcIdxVec::const_iterator am_begin() const {return am_to_do.begin();}

        PairArcIdxVec::const_iterator am_end() const {return am_to_do.end();}

        el_pat_vec pat_vec_at(pat_vec_t::size_type idx) const{
            assert(idx<pat_vec.size());
            return pat_vec[idx];
        }

        pat_vec_t::size_type pat_vec_size() const{return pat_vec.size();}

        pat_vec_t::const_iterator begin() const{return pat_vec.begin();}

        pat_vec_t::const_iterator end() const{return pat_vec.end();}


        pair_seqpos_t last_matched_pos(){
            assert(!pat_vec.empty());
            return pair_seqpos_t(pat_vec.back().first,pat_vec.back().second);
        }

        void add(seqpos_t posA,seqpos_t posB,char c){pat_vec.push_back(el_pat_vec(posA,posB,c));}

        void overwrite(seqpos_t posA,seqpos_t posB,char c,pat_vec_t::size_type pos){
            if(pat_vec.size()<=pos){pat_vec.push_back(el_pat_vec(posA,posB,c));}
            pat_vec.at(pos)=el_pat_vec(posA,posB,c);
        }

        void add_am(const Arc &a, const Arc &b){
            add(a.right(),b.right(),')');
            add(a.left(),b.left(),'(');
        }

        void store_am(const Arc &a, const Arc &b){
            const PairArcIdx &pair_arc_idx = PairArcIdx(a.idx(),b.idx());
            //store the pair of arc indices in the am_to_do datastructure
            am_to_do.push_back(pair_arc_idx);
        }

        PairArcIdx next_arcmatch(){
            PairArcIdx arc_idx = am_to_do.back();
            am_to_do.pop_back();
            return arc_idx;
        }

        void sort_patVec(){sort(pat_vec.begin(), pat_vec.end(),compare_el_pat_vec());}

        void sort_am_to_do(){
            sort(am_to_do.begin(), am_to_do.end(), compare_el_am_to_do());
        }

        void insert_epm(const EPM &epm_to_insert){
            pat_vec.insert(pat_vec.end(),epm_to_insert.begin(),epm_to_insert.end());
        }

        bool includes(const EPM &epm_to_test) const{
            assert(pat_vec_size()>=epm_to_test.pat_vec_size());
            return std::includes(this->begin(),this->end(),
                                 epm_to_test.begin(),epm_to_test.end(),
                                 compare_el_pat_vec());
        }

        bool includes_am(const EPM &epm_to_test) const{
            return std::includes(am_begin(),am_end(),
                                 epm_to_test.am_begin(),epm_to_test.am_end(),
                                 compare_el_am_to_do());
        }

        void print_epm(std::ostream &out, bool verbose) const{
            out << "_________________________________________________" << std::endl;
            out << "epm with score " << this->score << std::endl;
            out << " ";
            for(pat_vec_t::const_iterator it=pat_vec.begin();it!=pat_vec.end();++it){
                out << it->first << ":" << it->second << " ";
            }
            out << std::endl;
            out << " ";
            for(pat_vec_t::const_iterator it=pat_vec.begin();it!=pat_vec.end();++it){
                out << it->third;
            }
            out << std::endl;
            out << "am_to_do " << am_to_do << std::endl;
            out << "tolerance left " << this->max_tol_left << std::endl;
            if(verbose){
                out << "score " << score << std::endl;
                out << "pos " << this->cur_pos.first << "," << this->cur_pos.second << std::endl;
                out << "state " << this->state << std::endl;
            }
            out << "______________________________________________________" << std::endl;
        }
    };

    inline bool operator< (const EPM &epm1, const EPM &epm2) {
        return epm1.get_max_tol_left()>epm2.get_max_tol_left();
    }

    inline std::ostream & operator << (std::ostream &out, const EPM &epm){
        epm.print_epm(out,false);
        return out;
    }

    template <class T1>
    T1 max3(const T1 &first,const T1 &second, const T1 &third){
        return max(max(first,second),third);
    }

    template <class T1>
    T1 max4(const T1 &first,const T1 &second, const T1 &third, const T1 &fourth){
        return max(max3(first,second,third),fourth);
    }


    class ExactMatcher {
        typedef BasePairs__Arc Arc; 

        typedef SparsificationMapper::ArcIdx ArcIdx; 
        typedef SparsificationMapper::ArcIdxVec ArcIdxVec; 
        typedef SparsificationMapper::matidx_t matidx_t; 
        typedef SparsificationMapper::seq_pos_t seqpos_t; 
        typedef SparsificationMapper::index_t index_t; 
        typedef SparseTraceController::matpos_t matpos_t; 
        typedef SparseTraceController::pair_seqpos_t pair_seqpos_t; 

        typedef EPM::PairArcIdx PairArcIdx; 
        typedef EPM::PairArcIdxVec PairArcIdxVec; 

        typedef std::list<EPM> epm_cont_t; 
        typedef epm_cont_t::iterator epm_it_t; 
        typedef std::pair<score_t,epm_cont_t > el_map_am_to_do_t; 
        
        typedef unordered_map<PairArcIdx,el_map_am_to_do_t,pair_of_size_t_hash>::type map_am_to_do_t;
    private:

        //<state, max_tol, current matrix position, potential arcMatch, sequence position to be matched>
        typedef quintuple<int,infty_score_t,matpos_t,PairArcIdx,pair_seqpos_t> poss_L_LR;

        //infty_score_t because of the check_poss, change to score_t!!!
        typedef triple<int,infty_score_t,matpos_t> poss_in_G;//<state,max_tol,current matrix position>


        const Sequence &seqA; 
        const Sequence &seqB; 

        const RnaData &rna_dataA; 
        const RnaData &rna_dataB; 

        const ArcMatches &arc_matches; 

        const BasePairs &bpsA; 
        const BasePairs &bpsB; 
        const SparseTraceController &sparse_trace_controller; 
        // (valid positions in the matrix)
        const SparsificationMapper &sparse_mapperA; 
        const SparsificationMapper &sparse_mapperB; 
        PatternPairMap &foundEPMs; 

        ScoreMatrix L; 
        ScoreMatrix G_A; 
        ScoreMatrix G_AB; 
        ScoreMatrix LR; 
        ScoreMatrix F; 
        ScoreMatrix Dmat; 

        int alpha_1; 
        int alpha_2; 
        int alpha_3; 

        int difference_to_opt_score; 
        // worse than the optimal score are considered
        int min_score; 
        long int max_number_of_EPMs; 
        long int cur_number_of_EPMs; 

        bool inexact_struct_match; 
        score_t struct_mismatch_score; 
        bool add_filter; 

        bool verbose; 

        pair_seqpos_t pos_of_max; 

        enum{in_LR,in_G_A,in_G_AB,in_L,in_F}; 

        const Arc pseudo_arcA; 
        const Arc pseudo_arcB; 

        const
        infty_score_t &D(const ArcMatch &am) const {
            return D(am.arcA(),am.arcB());
        }

        infty_score_t &D(const ArcMatch &am) {
            return D(am.arcA(),am.arcB());
        }

        const
        infty_score_t &D(const Arc &a, const Arc &b) const {
            return Dmat(a.idx(),b.idx());
        }


        infty_score_t &D(const Arc &a, const Arc &b){
            return Dmat(a.idx(),b.idx());
        }

        bool nucleotide_match(seqpos_t pos_seqA, seqpos_t pos_seqB) const {
            assert(pos_seqA>=1 && pos_seqA<=seqA.length() &&
                   pos_seqB>=1 && pos_seqB<=seqB.length()); //seqA and seqB are 1-based
            return (seqA[pos_seqA]==seqB[pos_seqB]);
        }

        bool seq_matching(ArcIdx idxA, ArcIdx idxB, matpos_t cur_mat_pos, pair_seqpos_t cur_seq_pos) const {
            seqpos_t i = cur_seq_pos.first;
            seqpos_t j = cur_seq_pos.second;

            return sparse_trace_controller.pos_unpaired(idxA,idxB,cur_mat_pos) &&
                nucleotide_match(i,j);
        }

        // ----------------------------------------
        // fill matrices


        void initialize_gap_matrices();

        void init_Fmat();

        void init_mat(ScoreMatrix &mat, const Arc &a, const Arc &b,
                      infty_score_t first_entry,infty_score_t first_col, infty_score_t first_row);

        pair_seqpos_t compute_LGLR(const Arc &a, const Arc &b, bool suboptimal);

        infty_score_t compute_matrix_entry(const Arc &a, const Arc &b,matpos_t mat_pos,
                                           matpos_t mat_pos_diag, bool matrixLR, bool suboptimal);

        infty_score_t seq_str_matching(const Arc &a, const Arc &b, matpos_t mat_pos_diag,
                                       pair_seqpos_t seq_pos_to_be_matched, score_t add_score,bool matrixLR,
                                       bool suboptimal);

        void compute_F();

        // --------------------------------------------
        // helper functions

        score_t score_for_seq_match();

        infty_score_t score_for_am(const Arc &a, const Arc &b) const;

        score_t score_for_stacking(const Arc &a, const Arc &b,
                                   const Arc &inner_a,const Arc &inner_b);

        void add_foundEPM(EPM &cur_epm, bool count_EPMs);

        bool check_PPM(){
            if(this->difference_to_opt_score != -1) return true; //when we use the parameter difference_to_opt_score,
            //we enumerate all EPMs regardless of whether the number extends the max_number_of_EPMs
            if(cur_number_of_EPMs>=max_number_of_EPMs+1) return false;
            else return true;
        }

        void find_start_pos_for_tb(bool suboptimal, score_t difference_to_opt_score=-1, bool count_EPMs=false);

        bool check_num_EPMs(){
            double valid_deviation = 0.8;
            return (cur_number_of_EPMs>=max_number_of_EPMs*valid_deviation && cur_number_of_EPMs<=max_number_of_EPMs);
        }


        // --------------------------------------------
        // heuristic traceback

        void trace_F_heuristic(pos_type i, pos_type j,EPM &cur_epm);

        void trace_LGLR_heuristic(const Arc &a, const Arc &b, EPM &cur_epm);

        bool trace_seq_str_matching_heuristic(const Arc &a, const Arc &b,int &state,
                                              matpos_t &cur_mat_pos, matpos_t mat_pos_diag,
                                              pair_seqpos_t seq_pos_to_be_matched, score_t add_score);


        // --------------------------------------------
        // suboptimal traceback

        void trace_F_suboptimal(pos_type i,pos_type j,score_t max_tol, bool recurse, bool count_EPMs);

        void apply_filter(epm_cont_t &found_epms);

        void trace_LGLR_suboptimal(const Arc &a, const Arc &b, score_t max_tol,
                                   epm_cont_t &found_epms, bool recurse, bool count_EPMs);


        void trace_seq_str_matching_subopt(const Arc &a, const Arc &b,
                                           score_t score_contr, matpos_t mat_pos_diag, pair_seqpos_t seq_pos_to_be_matched,
                                           const PairArcIdx &am, poss_L_LR &poss, epm_it_t cur_epm, epm_cont_t &found_epms,
                                           map_am_to_do_t &map_am_to_do, bool count_EPMs);

        bool check_poss(const Arc &a, const Arc &b, const poss_L_LR &pot_new_poss,
                        poss_L_LR &poss, epm_it_t cur_epm, epm_cont_t &found_epms,
                        map_am_to_do_t &am_to_do_for_cur_am, bool count_EPMs);


        void store_new_poss(const Arc &a, const Arc &b, bool last_poss,
                            const poss_L_LR &new_poss, poss_L_LR &poss, epm_it_t cur_epm,
                            epm_cont_t &found_epms, map_am_to_do_t &am_to_do_for_cur_am, bool count_EPMs);

        void trace_G_suboptimal(const Arc &a, const Arc &b, const poss_L_LR &pot_new_poss,
                                poss_L_LR &poss, epm_it_t cur_epm, epm_cont_t &found_epms,
                                map_am_to_do_t &map_am_to_do, bool count_EPMs);

        bool is_valid_gap(const Arc &a, const Arc &b,
                          const poss_L_LR &pot_new_poss);

        void preproc_fill_epm(map_am_to_do_t &am_to_do, epm_it_t cur_epm,
                              epm_cont_t &found_epms, bool count_EPMs, score_t min_allowed_score=-1);

        void fill_epm(const map_am_to_do_t &map_am_to_do, size_type vec_idx,
                      std::vector<score_t> &max_tol_left_up_to_pos, std::vector<const EPM*> &epms_to_insert,
                      score_t min_score, epm_it_t cur_epm, epm_cont_t &found_epms,bool count_EPMs);

        // --------------------------------------------
        // debugging/testing

        // print the matrices in the condensed form
        void print_matrices(const Arc &a, const Arc &b, size_type offset_A, size_type offset_B, bool suboptimal,bool add_info);

        // checks whether an epm is valid, i.e. only one gap per arc match etc.
        bool validate_epm(const EPM &epm_to_test) const;

        // checks the validity of the epm list, i.e. that no epm is contained in another epm (all
        // epms are maximally extended)
        bool validate_epm_list(epm_cont_t &found_epms) const;

    public:


        ExactMatcher(const Sequence &seqA_,
                     const Sequence &seqB_,
                     const RnaData &rna_dataA_,
                     const RnaData &rna_dataB_,
                     const ArcMatches &arc_matches_,
                     const SparseTraceController &sparse_trace_controller_,
                     PatternPairMap &foundEPMs_,
                     int alpha_1_,
                     int alpha_2_,
                     int alpha_3_,
                     score_t difference_to_opt_score_,
                     score_t min_score_,
                     long int max_number_of_EPMs_,
                     bool inexact_struct_match_,
                     score_t struct_mismatch_score_,
                     bool apply_filter_,
                     bool opt_verbose_
                     );

        ~ExactMatcher();

        void compute_arcmatch_score();

        void test_arcmatch_score();

        void trace_EPMs(bool suboptimal);


    };


} //end namespace

#endif //  EXACT_MATCHER_HH