edge_probs.hh

#ifndef LOCARNA_EDGE_PROBS
#define LOCARNA_EDGE_PROBS
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
#include <string>
#include <cmath>
 
#include "aux.hh"
#include "matrix.hh"
 
namespace LocARNA {
 
    class StralScore;
    class RnaData;
    template <class T, size_t N> class Alphabet;
    class Sequence;
    class TraceController;
    class FreeEndgaps;
 
    class EdgeProbs {
    public:
        using size_type=size_t; 
 
        EdgeProbs(std::istream &in, size_type lenA, size_type lenB) {
            read_sparse(in, lenA, lenB);
        }
 
        std::ostream &
        write_sparse(std::ostream &out, double threshold) const;
 
        size_type
        lenA() const {
            return probs_.sizes().first;
        }
 
        size_type
        lenB() const {
            return probs_.sizes().second;
        }
 
        double
        prob(size_t i, size_t j) const {
            assert(i < probs_.sizes().first);
            assert(j < probs_.sizes().second);
 
            return probs_(i, j);
        }
 
    protected:
        Matrix<double> probs_; 
 
        std::istream &
        read_sparse(std::istream &in, size_type lenA, size_type lenB);
 
        void
        read_sparse(const std::string &filename,
                    size_type lenA,
                    size_type lenB);
        EdgeProbs() {}
    };
 
    class MatchProbs : public EdgeProbs {
    public:
        MatchProbs(std::istream &in, size_type lenA, size_type lenB)
            : EdgeProbs(in, lenA, lenB)
        {
        }
    protected:
        MatchProbs() {}
    };
 
    template <class pf_score_t>
    class PFGotoh {
    public:
        using size_type=size_t; 
 
        PFGotoh(const RnaData &rnaA,
                const RnaData &rnaB,
                const TraceController &trace_controller,
                const Matrix<double> &sim_mat,
                const Alphabet<char, 4> &alphabet,
                double gap_opening,
                double gap_extension,
                double pf_struct_weight,
                double temp,
                const FreeEndgaps &free_endgaps,
                bool flag_local);
        const pf_score_t&
        z() const {return z_;}
 
    protected:
 
        size_type lenA_;
        size_type lenB_;
 
        double temp_;
 
        bool flag_local_;
 
        pf_score_t z_;
 
        Matrix<pf_score_t> zM_;
 
        Matrix<pf_score_t> zA_;
 
        Matrix<pf_score_t> zB_;
 
        Matrix<pf_score_t> zMr_; 
        Matrix<pf_score_t> zAr_; 
        Matrix<pf_score_t> zBr_; 
 
        void
        pf_gotoh(Matrix<pf_score_t> &zM,
                 Matrix<pf_score_t> &zA,
                 Matrix<pf_score_t> &zB,
                 const TraceController &trace_controller,
                 const StralScore &score,
                 const FreeEndgaps &free_endgaps);
 
        bool fail() const {return z_==(pf_score_t)0 || std::isnan(z_) || std::isinf(z_);}
 
    };
 
    class TraceProbs : public EdgeProbs {
    public:
        TraceProbs(std::istream &in, size_type lenA, size_type lenB)
            : EdgeProbs(in, lenA, lenB)
        {}
    protected:
        TraceProbs() {}
    };
 
    template <class  pf_score_t>
    class PFMatchProbs : public MatchProbs, PFGotoh<pf_score_t> {
    public:
        using size_type=EdgeProbs::size_type; 
 
        PFMatchProbs(const RnaData &rnaA,
                     const RnaData &rnaB,
                     const TraceController &trace_controller,
                     const Matrix<double> &sim_mat,
                     const Alphabet<char, 4> &alphabet,
                     double gap_opening,
                     double gap_extension,
                     double pf_struct_weight,
                     double temp,
                     const FreeEndgaps &free_endgaps,
                     bool flag_local);
 
        bool fail() const {return PFGotoh<pf_score_t>::fail();}
 
        const pf_score_t &
        z() const {return PFGotoh<pf_score_t>::z();}
    };
 
 
    template <class pf_score_t>
    class PFTraceProbs : public TraceProbs, PFGotoh<pf_score_t> {
    public:
        using size_type=TraceProbs::size_type; 
 
        PFTraceProbs(const RnaData &rnaA,
                     const RnaData &rnaB,
                     const TraceController &trace_controller,
                     const Matrix<double> &sim_mat,
                     const Alphabet<char, 4> &alphabet,
                     double gap_opening,
                     double gap_extension,
                     double pf_struct_weight,
                     double temp,
                     const FreeEndgaps &free_endgaps,
                     bool flag_local);
 
        bool fail() const {return PFGotoh<pf_score_t>::fail();}
 
        const pf_score_t &
        z() const {return PFGotoh<pf_score_t>::z();}
    };
 
    /*
     * @brief Provide match probabilities calculated by pair HMM (deprecated)
     *
     * For computing the probabilities the class
     * uses a pairHMM analogously to PROBCONS.
     * Also, the class reads transition probabilities
     * from a file in the format of Probcons, which
     * allows to use their parameter files
     *
     * @note the pairHMM approach does not work very well yet (parametrization
     * is not obvious, correctness of implementation is still critical, uses
     * only first sequence of multiple alignments,
     * numerical problems without log transformation). Instead of fixing
     * the issues of the pairHMM implementation, the probalign-like partition
     * function approach was implemented.
     */
    class PairHMMMatchProbs : public MatchProbs {
    public:
 
        class PairHMMParams {
        public:
            // ------------------------------------------------------------
            // transition probabilities
            // there are three states M, X, Y (plus implicitely start and end
            // state)
            double initM;       
            double initX;       
            double initY;       
            double startX;      
            double startY;      
            double extendM;     
            double extendX;     
            double extendY;     
            double startMFromX; 
            double startMFromY; 
 
            std::string basenames; 
 
            Matrix<double> emmission; 
            std::vector<double>
                background; 
 
            explicit PairHMMParams(const std::string &filename);
        };
 
        PairHMMMatchProbs(const Sequence &seqA,
                          const Sequence &seqB,
                          const PairHMMParams &params) {
            pairHMM_probs(seqA, seqB, params);
        }
 
    protected:
        void
        pairHMM_probs(const Sequence &seqA,
                      const Sequence &seqB,
                      const PairHMMParams &params);
    };
}
 
#include "edge_probs.icc"
 
 
#endif // LOCARNA_EDGE_PROBS