src/LocARNA/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,
                        const 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(const 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 verbose_);

        ~ExactMatcher();

        void
        compute_arcmatch_score();

        void
        test_arcmatch_score();

        void
        trace_EPMs(bool suboptimal);
    };

} // end namespace

#endif //  EXACT_MATCHER_HH