arc_matches.hh

#ifndef LOCARNA_ARC_MATCHES_HH
#define LOCARNA_ARC_MATCHES_HH

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

#include <algorithm>
#include <vector>

#include "scoring_fwd.hh"
#include "aux.hh"
#include "matrix.hh"
#include "basepairs.hh"

#include <assert.h>


namespace LocARNA {
    class Scoring;
    class Sequence;
    class RnaData;
    class AnchorConstraints;
    class TraceController;
    class MatchController;
    
    class ArcMatch {
    public:
        typedef std::vector<int>::size_type size_type; 
        typedef size_type idx_type; 
        typedef BasePairs__Arc Arc; 
    private:
        const Arc *arcA_; 
        const Arc *arcB_; 
        idx_type idx_;    

    public:
        
        ArcMatch(const Arc *arcA,const Arc *arcB, idx_type idx)
            : arcA_(arcA),
              arcB_(arcB),
              idx_(idx)
        {}
    
        const Arc &
        arcA() const {return *arcA_;}
    
        const Arc &
        arcB() const {return *arcB_;}
    
        idx_type
        idx()  const {return idx_;}
    };

    typedef std::vector<ArcMatch> ArcMatchVec;

    typedef std::vector<ArcMatch::idx_type> ArcMatchIdxVec;

    class ArcMatches {
    public:
        typedef std::vector<int>::size_type size_type; 
        typedef BasePairs__Arc Arc; 
    protected:
        
        size_type lenA; 
        size_type lenB; 
        
        BasePairs *bpsA; 
        BasePairs *bpsB; 
    
        /* Constraints and Heuristics */
    
        size_type max_length_diff; 
    
        size_type max_diff_at_am; 
    
        const MatchController &match_controller; 
    
        const AnchorConstraints &constraints; 
    
    
        bool is_valid_arcmatch(const Arc &arcA,const Arc &arcB) const;
    
        /* END constraints and heuristics */
    
        bool maintain_explicit_scores; 
    
        ArcMatchVec arc_matches_vec;

        size_type number_of_arcmatches;

        std::vector<score_t> scores;
    

        Matrix<ArcMatchIdxVec> common_right_end_lists;
       
    
        Matrix<ArcMatchIdxVec> common_left_end_lists;
    
    
        ArcMatchIdxVec inner_arcmatch_idxs;
    
        void
        init_inner_arc_matchs();
    
        class lex_greater_left_ends {
            const ArcMatches &arc_matches;
        public:
            
            explicit
            lex_greater_left_ends(const ArcMatches &arc_matches_)
                : arc_matches(arc_matches_)
            {}
            
            bool
            operator () (const ArcMatch::idx_type &i, const ArcMatch::idx_type &j) const {
                size_type ali = arc_matches.arcmatch(i).arcA().left();
                size_type bli = arc_matches.arcmatch(i).arcB().left();
                size_type alj = arc_matches.arcmatch(j).arcA().left();
                size_type blj = arc_matches.arcmatch(j).arcB().left();
            
                return (ali>alj) || (ali==alj && bli>blj);
            }
        };

        class tuple5 {
        public:
            typedef std::vector<int>::size_type size_type; 

            size_type i; 
            size_type j; 
            size_type k; 
            size_type l; 
            score_t score; 
        
            tuple5(size_type i_,size_type j_,size_type k_,size_type l_,score_t score_)
                :i(i_), j(j_), k(k_), l(l_), score(score_)
            {}
        };
            
    
    public:

        
        ArcMatches(const Sequence &seqA_, 
                   const Sequence &seqB_, 
                   const std::string &arcmatch_scores_file,
                   int probability_scale,
                   size_type max_length_diff,
                   size_type max_diff_at_am,
                   const MatchController &trace_controller,
                   const AnchorConstraints &constraints);
        
        
        
        ArcMatches(const RnaData &rnadataA, 
                   const RnaData &rnadataB,
                   double min_prob,
                   size_type max_length_diff,
                   size_type max_diff_at_am,
                   const MatchController &trace_controller,
                   const AnchorConstraints &constraints);
    
        ~ArcMatches();
    
        // for the mea probabilistic consistency transformation, support to read and write the arcmatch scores
        // this allows in general to have user defined arc-match scores
    
        void read_arcmatch_scores(const std::string &arcmatch_scores_file, int probability_scale);
    
    
        void write_arcmatch_scores(const std::string &arcmatch_scores_file, const Scoring &scoring) const;
    
    
        const BasePairs &
        get_base_pairsA() const {
            return *bpsA;
        }

        const BasePairs &
        get_base_pairsB() const {
            return *bpsB;
        }

        bool
        explicit_scores() const {
            return maintain_explicit_scores;
        }
    
        void
        make_scores_explicit(const Scoring &scoring);
        
        
        score_t
        get_score(const ArcMatch &am) const {
            assert(maintain_explicit_scores);
            return scores[am.idx()];
        }
        
        size_type num_arc_matches() const {
            return number_of_arcmatches;
        }
    
        const ArcMatch &arcmatch(size_type idx) const {
            assert(idx<number_of_arcmatches);
            return arc_matches_vec[idx];
        }
        
        // ============================================================
        // Iteration over arc matches
        //
    
        const ArcMatchIdxVec &
        common_right_end_list(size_type i, size_type j) const {
            return common_right_end_lists(i,j);
        }
    
        const ArcMatchIdxVec &
        common_left_end_list(size_type i, size_type j) const {
            return common_left_end_lists(i,j);
        }
    
    
        // ============================================================
    

        void get_max_right_ends(size_type al,size_type bl,size_type *max_ar,size_type *max_br, bool no_lonely_pairs) const; 

        
        void get_min_right_ends(size_type al,size_type bl,size_type *min_ar,size_type *min_br) const; 
    
        // ------------------------------------------------------------
        // inner arc matches

    
        bool
        exists_inner_arc_match(const ArcMatch &am) const {
            return inner_arcmatch_idxs[am.idx()] < num_arc_matches();
        }
        
        const ArcMatch &
        inner_arc_match(const ArcMatch &am) const {
            return arcmatch(inner_arcmatch_idxs[am.idx()]);
        }
        
        void
        sort_right_adjacency_lists();

        // ------------------------------------------------------------
        // iteration (in no specific order)

        typedef ArcMatchVec::const_iterator const_iterator;

        const_iterator begin() const {return arc_matches_vec.begin();}

        const_iterator end() const {return arc_matches_vec.begin()+number_of_arcmatches;}
    
    };

    
    class ArcMatchesIndexed : public ArcMatches {
    public:
        ArcMatchesIndexed(const Sequence &seqA_, 
                          const Sequence &seqB_, 
                          const std::string &arcmatch_scores_file,
                          int probability_scale,
                          size_type max_length_diff,
                          size_type max_diff_at_am,
                          const MatchController &trace_controller,
                          const AnchorConstraints &constraints)
            :ArcMatches(seqA_,seqB_,
                        arcmatch_scores_file,
                        probability_scale,
                        max_length_diff,
                        max_diff_at_am,
                        trace_controller,
                        constraints),
             am_index_()
        {
            build_arcmatch_index();
        }
        
        ArcMatchesIndexed(const RnaData &rnadataA,
                          const RnaData &rnadataB,
                          double min_prob,
                          size_type max_length_diff,
                          size_type max_diff_at_am,
                          const MatchController &trace_controller,
                          const AnchorConstraints &constraints)
            :ArcMatches(rnadataA,
                        rnadataB,
                        min_prob,
                        max_length_diff,
                        max_diff_at_am,
                        trace_controller,
                        constraints),
             am_index_()
        {
            build_arcmatch_index();
        }
        
    private:
        typedef std::pair<size_type,size_type> idx_pair_t;
        
        typedef unordered_map<idx_pair_t,ArcMatch::idx_type,pair_of_size_t_hash>::type am_index_type;
        
        am_index_type am_index_;
        
        void
        build_arcmatch_index();

        
    public:
        
        const ArcMatch::idx_type invalid_am_index() const {
            // The invalid arc match index is implemented to be the
            // maximum valid index + 1.
            // This allows an efficient implementation, where we push 
            // an invalid arc match to the end of vector arc_matches_vec.
            // Thus, we return size-1!

            return number_of_arcmatches;
        }
        
        const ArcMatch::idx_type
        am_index(const size_type &arcAIdx,const size_type &arcBIdx) const {
            am_index_type::const_iterator it = am_index_.find(idx_pair_t(arcAIdx,arcBIdx));
            if (am_index_.end() != it) {
                return it->second;
            } else {
                return invalid_am_index();
            }
        }

        const ArcMatch &
        am_index(const Arc &arcA,const Arc &arcB) const {
            return arc_matches_vec[am_index(arcA.idx(),arcB.idx())];
        }
    };

} // end namespace LocARNA

#endif // LOCARNA_ARC_MATCHES_HH