sparsification_mapper.hh

#ifndef SPARSIFICATION_MAPPER_HH
#define SPARSIFICATION_MAPPER_HH

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

#include <iostream>
#include <limits>

#include "aux.hh"
#include "basepairs.hh"
#include "ext_rna_data.hh"

// use for type safe index_t
// #include "type_wrapper.hh"


namespace LocARNA {

    // print a pair
    template <class T1, class T2>
    std::ostream& operator << (std::ostream& out, const std::pair<T1,T2>& pair) {
        return out << "(" << pair.first << "," << pair.second << ") ";
    }
    

class SparsificationMapper{

public:
    typedef BasePairs__Arc Arc; 
        typedef size_t ArcIdx; 
        typedef std::vector<ArcIdx> ArcIdxVec; 
        typedef pos_type matidx_t; 
        typedef pos_type seq_pos_t; 
    
        // note: the type safe index_t breaks current code, since
        // casts to size_t need to be explicite (using index_t's val()-method)
        // //! type-safe index type this is useful to distinguish index type
        // //! from other types that are defined as unsigned int
        // typedef type_wrapper<size_t> index_t;
    
        typedef size_t index_t; 

        struct info_for_pos{
                seq_pos_t seq_pos; 
                bool unpaired; 
                ArcIdxVec valid_arcs; 

                void reset(){
                        this->seq_pos=0;
                        this->valid_arcs.clear();
                        this->unpaired=false;
                }
        };

        typedef std::vector<info_for_pos> InfoForPosVec;

private:

        const BasePairs &bps; 
        const ExtRnaData &rnadata; 
        const double prob_unpaired_in_loop_threshold; 
        const double prob_basepair_in_loop_threshold; 
        size_type max_info_vec_size; 
        std::vector<InfoForPosVec> info_valid_seq_pos_vecs;

        std::vector<std::vector<matidx_t> > valid_mat_pos_vecs_before_eq;

        std::vector<std::vector<ArcIdxVec > > left_adj_vec;

        void compute_mapping_idx_arcs();

        void compute_mapping_idx_left_ends();

        void iterate_left_adj_list(pos_type cur_left_end,
                                   pos_type cur_pos,
                                   const Arc *inner_arc,
                                   info_for_pos &struct_pos);

        void valid_pos_external(pos_type cur_pos,const Arc *inner_arc, info_for_pos &struct_pos);



public:
        SparsificationMapper(const BasePairs &bps_,
                        const ExtRnaData &rnadata_,
                        double prob_unpaired_in_loop_threshold_,
                        double prob_basepair_in_loop_threshold_,
                        bool index_left_ends):
                                bps(bps_),
                                rnadata(rnadata_),
                                prob_unpaired_in_loop_threshold(prob_unpaired_in_loop_threshold_),
                                prob_basepair_in_loop_threshold(prob_basepair_in_loop_threshold_),
                                max_info_vec_size(0)
        {
                if(index_left_ends){
                        compute_mapping_idx_left_ends();
                }
                else{
                        compute_mapping_idx_arcs();
                }
        }

        size_type get_max_info_vec_size() const
        {
            return max_info_vec_size;
        }

        const InfoForPosVec &
        valid_seq_positions(index_t idx) const{
                return info_valid_seq_pos_vecs.at(idx);
        }

        const ArcIdxVec &
        valid_arcs_right_adj(index_t idx, matidx_t pos) const {
                return info_valid_seq_pos_vecs.at(idx).at(pos).valid_arcs;
        }

         matidx_t 
         first_valid_mat_pos_before_eq(index_t index, seq_pos_t pos,
                                       index_t left_end = 
                                       std::numeric_limits<index_t>::max()) const {
            if (left_end == std::numeric_limits<index_t>::max())
                left_end = index;
            assert (pos >= left_end); //tocheck
            return valid_mat_pos_vecs_before_eq.at(index).at(pos-left_end);
        }

        inline
        matidx_t first_valid_mat_pos_before(index_t index, seq_pos_t pos,
                                            index_t left_end =
                                            std::numeric_limits<index_t>::max()) const {
            // if (left_end == std::numeric_limits<index_t>::max()) assert (pos > index);
            return first_valid_mat_pos_before_eq(index, pos-1, left_end);
        }

        inline
        seq_pos_t get_pos_in_seq_new(index_t idx, matidx_t pos) const{
                assert(pos<number_of_valid_mat_pos(idx));
                return (info_valid_seq_pos_vecs.at(idx).at(pos).seq_pos);//+arc.left();
        }

        size_type number_of_valid_mat_pos(index_t idx) const{
                return info_valid_seq_pos_vecs.at(idx).size();
        }

        bool pos_unpaired(index_t idx,matidx_t pos)const{
                return info_valid_seq_pos_vecs.at(idx).at(pos).unpaired;
        }

        bool matching_without_gap_possible(const Arc &arc, seq_pos_t pos)const{
                const pos_type &mat_pos = first_valid_mat_pos_before(arc.idx(),pos,arc.left());
                return get_pos_in_seq_new(arc.idx(),mat_pos)==pos-1;
        }

        const ArcIdxVec &
        valid_arcs_left_adj(const Arc &arc, seq_pos_t pos) const{
                return left_adj_vec.at(arc.idx()).at(pos-arc.left());
        }

        virtual ~SparsificationMapper(){
        }

        matidx_t idx_geq(index_t index, seq_pos_t min_col, index_t left_end = std::numeric_limits<index_t>::max()) const{

                if (left_end == std::numeric_limits<index_t>::max())
                        left_end = index;

                size_t num_pos = number_of_valid_mat_pos(index);
                seq_pos_t last_mapped_seq_pos = get_pos_in_seq_new(index,num_pos-1);

                // if the position min_col is smaller than or equal to the left end (first mapped position), return 0
                if(min_col<=left_end) return 0;

                // if the position min_col is greater than the last mapped sequence position, the position does not exists
                // return the number of valid positions
                if(min_col>last_mapped_seq_pos) return num_pos;

                // the matrix position after the first valid position before the position min_col (first matrix position that
                // stores a sequence position that is greater than or equal to min_col-1)
                matidx_t idx_geq = first_valid_mat_pos_before_eq(index,min_col-1,left_end)+1;

                assert(get_pos_in_seq_new(index,idx_geq)>=min_col && !(get_pos_in_seq_new(index,idx_geq-1)>=min_col));

                return idx_geq;
        }

        matidx_t idx_after_leq(index_t index, seq_pos_t max_col, index_t left_end = std::numeric_limits<index_t>::max()) const{

                if (left_end == std::numeric_limits<index_t>::max())
                        left_end = index;

                size_t num_pos = number_of_valid_mat_pos(index);
                seq_pos_t last_mapped_seq_pos = get_pos_in_seq_new(index,num_pos-1);

                // if the position max_col is smaller than the left_end (first mapped position), return 0
                if(max_col<left_end) return 0;

                // if the position max_col is greater than or equal to the last mapped sequence position,
                // return the number of positions
                if(max_col>=last_mapped_seq_pos) return num_pos;

                // the matrix position after the first matrix position before or equal the sequence position max_col
                // (last matrix position that stores a sequence position that is lower than or equal max_col)
                matidx_t idx_after_leq = first_valid_mat_pos_before_eq(index,max_col,left_end)+1;

                assert(get_pos_in_seq_new(index,idx_after_leq-1)<=max_col && !(get_pos_in_seq_new(index,idx_after_leq)<=max_col));

                return idx_after_leq;
        }

private:

        bool is_valid_arc(const Arc &inner_arc,const Arc &arc)const{
                assert(inner_arc.left()>arc.left() && inner_arc.right()<arc.right());
                return rnadata.arc_in_loop_prob(inner_arc.left(),inner_arc.right(),arc.left(),arc.right())>=prob_basepair_in_loop_threshold;
        }

        bool is_valid_arc_external(const Arc &inner_arc)const{
                        return rnadata.arc_external_prob(inner_arc.left(),inner_arc.right())>=prob_basepair_in_loop_threshold;
        }

public:
        bool is_valid_pos(const Arc &arc,seq_pos_t pos) const{
            assert(arc.left()<pos && pos<arc.right());
            return rnadata.unpaired_in_loop_prob(pos,arc.left(),arc.right())>=prob_unpaired_in_loop_threshold; //todo: additional variable for external case
        }

        bool is_valid_pos_external(seq_pos_t pos) const{
            return rnadata.unpaired_external_prob(pos)>=prob_unpaired_in_loop_threshold; //todo: additional variable for external case
        }
};

template <class T>
std::ostream& operator<<(std::ostream& out, const std::vector<T>& vec){
        for(typename std::vector<T>::const_iterator it = vec.begin();it!=vec.end();it++){
                out << *it << " ";
        }
        return out;
}

std::ostream &operator << (std::ostream &out, const std::vector<SparsificationMapper::InfoForPosVec> &pos_vecs_);

std::ostream &operator << (std::ostream &out, const SparsificationMapper::InfoForPosVec &pos_vec_);

} //end namespace

#endif //  SPARSIFICATION_MAPPER_HH