src/LocARNA/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