scoring.hh

#ifndef LOCARNA_SCORING_HH
#define LOCARNA_SCORING_HH
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
#include <cmath>
#include <vector>
 
#include "aux.hh"
 
#include "scoring_fwd.hh"
#include "matrix.hh"
#include "sequence.hh"
#include "arc_matches.hh"
#include "named_arguments.hh"
 
namespace LocARNA {
 
    //#define MEA_SCORING_OLD
 
    class RibosumFreq;
    class Ribofit;
    class Scoring;
    class Sequence;
    class BasePairs;
    class BasePairs__Arc;
    class ArcMatch;
    class MatchProbs;
    class RnaData;
 
    typedef std::vector<infty_score_t> ScoreVector;
 
    typedef Matrix<infty_score_t> ScoreMatrix;
 
    typedef Matrix<double> ProbMatrix;
 
    class ScoringParams {
    public:
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(match, score_t, 50);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(mismatch, score_t, 0);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(indel, score_t, -150);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(indel_opening, score_t, -750);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(indel_loop, score_t, -200);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(indel_opening_loop, score_t, -800);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(ribosum, const RibosumFreq *, nullptr);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(ribofit, const Ribofit *, nullptr);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(unpaired_penalty, score_t, 0);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(struct_weight, score_t, 200);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(tau_factor, score_t, 50);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(exclusion, score_t, 0);
 
        DEFINE_NAMED_ARG_FEATURE(exp_probA, double);
 
        DEFINE_NAMED_ARG_FEATURE(exp_probB, double);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(temperature_alipf, double, 150);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(stacking, bool, false);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(new_stacking, bool, false);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(mea_scoring, bool, false);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(mea_alpha, score_t, 0);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(mea_beta, score_t, 200);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(mea_gamma, score_t, 100);
 
        DEFINE_NAMED_ARG_DEFAULT_FEATURE(probability_scale, score_t, 10000);
 
        using valid_args = std::tuple<
            match,
            mismatch,
            indel,
            indel_loop,
            indel_opening,
            indel_opening_loop,
            ribosum,
            ribofit,
            unpaired_penalty,
            struct_weight,
            tau_factor,
            exclusion,
            exp_probA,
            exp_probB,
            temperature_alipf,
            stacking,
            new_stacking,
            mea_scoring,
            mea_alpha,
            mea_beta,
            mea_gamma,
            probability_scale
            >;
 
        /*
          The mea score is composed in the following way:
 
          params->probability_scale_ *
          (
          sum_basematchs (i,j)
          P(i~j)
          + mea_alpha/100 * (P(i unstructured) + P(j unstructured))
          +
          sum_arcmatchs (a,b)
          mea_beta/100 * (P(a)+P(b)) * P(al~bl) * P(ar~br) *
          ribosum_arcmatch(a,b)
          )
        */
 
    public:
        template <typename... Args>
        ScoringParams(Args... argpack) {
            static_assert( type_subset_of<
                           std::tuple<Args...>,
                           valid_args>::value,
                           "Invalid type in named arguments pack." );
 
            auto args = std::make_tuple(argpack...);
 
            //mandatory
            exp_probA_ = get_named_arg<exp_probA>(args);
            exp_probB_ = get_named_arg<exp_probB>(args);
 
            //optional
            match_ = get_named_arg_opt<match>(args);
            mismatch_ = get_named_arg_opt<mismatch>(args);
            indel_ = get_named_arg_opt<indel>(args);
            indel_loop_ = get_named_arg_opt<indel_loop>(args);
            indel_opening_ = get_named_arg_opt<indel_opening>(args);
            indel_opening_loop_ = get_named_arg_opt<indel_opening_loop>(args);
            ribosum_ = get_named_arg_opt<ribosum>(args);
            ribofit_ = get_named_arg_opt<ribofit>(args);
            unpaired_penalty_ = get_named_arg_opt<unpaired_penalty>(args);
            struct_weight_ = get_named_arg_opt<struct_weight>(args);
            tau_factor_ = get_named_arg_opt<tau_factor>(args);
            exclusion_ = get_named_arg_opt<exclusion>(args);
            temperature_alipf_ = get_named_arg_opt<temperature_alipf>(args);
            stacking_ = get_named_arg_opt<stacking>(args);
            new_stacking_ = get_named_arg_opt<new_stacking>(args);
            mea_scoring_ = get_named_arg_opt<mea_scoring>(args);
            mea_alpha_ = get_named_arg_opt<mea_alpha>(args);
            mea_beta_ = get_named_arg_opt<mea_beta>(args);
            mea_gamma_ = get_named_arg_opt<mea_gamma>(args);
            probability_scale_ = get_named_arg_opt<probability_scale>(args);
        }
    };
 
    class Scoring {
    public:
        typedef BasePairs__Arc Arc; 
 
    protected:
        const ScoringParams *params; 
 
        const ArcMatches *arc_matches_; 
 
        const MatchProbs *match_probs; 
 
        const RnaData &rna_dataA; 
        const RnaData &rna_dataB; 
        const Sequence &seqA;     
        const Sequence &seqB;     
 
        score_t lambda_;
 
    public:
        Scoring(const Sequence &seqA,
                const Sequence &seqB,
                const RnaData &rna_dataA,
                const RnaData &rna_dataB,
                const ArcMatches &arc_matches,
                const MatchProbs *match_probs,
                const ScoringParams &params);
 
        void
        modify_by_parameter(score_t lambda);
 
        void
        apply_unpaired_penalty();
 
        score_t
        lambda() const {
            return lambda_;
        }
 
    protected:
        // ------------------------------
        // tables for precomputed score contributions
        //
        Matrix<score_t>
            sigma_tab; 
 
        std::vector<score_t> gapcost_tabA; 
        std::vector<score_t> gapcost_tabB; 
 
        std::vector<score_t> weightsA; //<! weights of base pairs in A
        std::vector<score_t> weightsB; //<! weights of base pairs in B
 
        std::vector<score_t> stack_weightsA; //<! weights of stacked base
        //<! pairs in A
        std::vector<score_t> stack_weightsB; //<! weights of stacked base
        //<! pairs in B
 
        Matrix<size_t> identity; 
 
        void
        precompute_sequence_identities();
 
        score_t
        round2score(double d) const {
            return (score_t)((d < 0) ? (d - 0.5) : (d + 0.5));
        }
 
        score_t
        sigma_(int i, int j) const;
 
        void
        precompute_sigma();
 
        void
        precompute_gapcost();
 
        void
        precompute_weights();
 
        void
        precompute_weights(const RnaData &rna_data,
                           const BasePairs &bps,
                           double exp_prob,
                           std::vector<score_t> &weights,
                           std::vector<score_t> &stack_weights);
 
        score_t
        probToWeight(double p, double prob_exp) const;
 
        double
        ribosum_arcmatch_prob(const Arc &arcA, const Arc &arcB) const;
 
        score_t
        riboX_arcmatch_score(const Arc &arcA, const Arc &arcB) const;
 
        void
        subtract(std::vector<score_t> &v, score_t x) const;
 
        void
        subtract(Matrix<score_t> &m, score_t x) const;
 
    public:
        // ------------------------------------------------------------
        // SCORE CONTRIBUTIONS
 
        score_t
        basematch(size_type i, size_type j) const {
            return sigma_tab(i, j);
        }
 
        score_t
        arcmatch(const ArcMatch &am, bool stacked = false) const;
 
        score_t
        arcmatch(const BasePairs__Arc &arcA,
                 const BasePairs__Arc &arcB,
                 bool stacked = false) const;
 
        const ArcMatches *
        arc_matches() const {
            return arc_matches_;
        }
 
        template <bool gapAorB>
        score_t
        arcDel(const BasePairs__Arc &arc, bool stacked = false) const;
 
        score_t
        arcmatch_stacked(const ArcMatch &am) const {
            return arcmatch(am, true);
        }
 
        template <bool gapInA>
        inline
        score_t
        gapX(size_type alignedToGap) const;
 
        score_t
        gapA(size_type posA) const {
            assert(1 <= posA && posA <= seqA.length());
 
            return gapcost_tabA[posA];
        }
 
        score_t
        gapB(size_type posB) const {
            assert(1 <= posB && posB <= seqB.length());
 
            return gapcost_tabB[posB];
        }
 
        score_t
        exclusion() const {
            return params->exclusion_;
        }
 
        score_t
        indel_opening() const {
            return params->indel_opening_;
        }
 
        score_t
        loop_indel_score(const score_t score) const {
            return round2score(score * params->indel_loop_ / params->indel_);
        }
        score_t
        indel_opening_loop() const {
            return params->indel_opening_loop_;
        }
 
        //
        // ------------------------------------------------------------
 
        double
        prob_exp(size_type len) const;
 
        bool
        stacking() const {
            return params->stacking_ || params->new_stacking_;
        }
 
        bool
        is_stackable_arcA(const Arc &a) const;
 
        bool
        is_stackable_arcB(const Arc &a) const;
 
        bool
        is_stackable_am(const ArcMatch &am) const;
 
    }; // end class Scoring
 
    template <bool isA>
    score_t
    Scoring::arcDel(const Arc &arcX,
                    bool stacked) const { // TODO Important Scoring scheme for
                                          // aligning an arc to a gap is not
                                          // defined and implemented!
 
        if (arc_matches_
                ->explicit_scores()) { // will not take stacking into account!!!
            std::cerr
                << "ERROR sparse explicit scores is not supported!"
                << std::endl; // TODO: Supporting explicit scores for arcgap
            assert(!arc_matches_->explicit_scores());
        }
 
        if (!params->mea_scoring_) {
            return
                // base pair weights
                loop_indel_score(
                    gapX<isA>(arcX.left()) +
                    gapX<isA>(arcX.right())) + // score of aligining base-pair
                                               // ends wo gap, such that it is
                                               // multiply by gamma_loop ratio
                (stacked ? (isA ? stack_weightsA[arcX.idx()]
                                : stack_weightsB[arcX.idx()])
                         : (isA ? weightsA[arcX.idx()] : weightsB[arcX.idx()]));
        } else {
            std::cerr << "ERROR sparse mea_scoring is not supported!"
                      << std::endl; // TODO: Supporting mea_scoring for arcgap
            assert(!params->mea_scoring_);
            return 0;
        }
    }
 
    template <>
    inline
    score_t
    Scoring::gapX<true>(size_type alignedToGap) const {
        return gapA(alignedToGap);
    }
 
    template <>
    inline
    score_t
    Scoring::gapX<false>(size_type alignedToGap) const {
        return gapB(alignedToGap);
    }
 
 
    template<typename T>
    class PFScoring : public Scoring {
        public:
        using pf_score_t = T;
 
        using PFScoreVector = std::vector<pf_score_t>;
 
        using PFScoreMatrix = Matrix<pf_score_t>;
 
 
        PFScoring(const Sequence &seqA,
                const Sequence &seqB,
                const RnaData &rna_dataA,
                const RnaData &rna_dataB,
                const ArcMatches &arc_matches,
                const MatchProbs *match_probs,
                const ScoringParams &params);
 
        pf_score_t
        exp_basematch(size_type i, size_type j) const {
            return exp_sigma_tab(i, j);
        }
 
        pf_score_t
        exp_arcmatch(const ArcMatch &am) const {
            return boltzmann_weight(arcmatch(am));
        }
 
        pf_score_t
        exp_gapA(size_type posA) const {
            assert(1 <= posA && posA <= seqA.length());
            return exp_gapcost_tabA[posA];
        }
 
        pf_score_t
        exp_gapB(size_type posB) const {
            assert(1 <= posB && posB <= seqB.length());
 
            return exp_gapcost_tabB[posB];
        }
 
        pf_score_t
        exp_indel_opening() const {
            return exp_indel_opening_score;
        }
 
        pf_score_t
        exp_indel_opening_loop() const {
            return exp_indel_opening_loop_score;
        }
 
 
    protected:
        void
        precompute_exp_gapcost();
 
 
        void
        precompute_exp_sigma();
 
 
        pf_score_t
        boltzmann_weight(score_t s) const {
            return std::exp(s / (pf_score_t)params->temperature_alipf_);
        }
 
 
    private:
        // ------------------------------
        // tables for precomputed exp score contributions for partition function
        //
        Matrix<pf_score_t>
            exp_sigma_tab; 
        pf_score_t exp_indel_opening_score; 
        pf_score_t exp_indel_opening_loop_score; 
        std::vector<pf_score_t>
            exp_gapcost_tabA; 
        std::vector<pf_score_t>
            exp_gapcost_tabB; 
 
    };
 
    template <typename T>
    PFScoring<T>::PFScoring(const Sequence &seqA,
                         const Sequence &seqB,
                         const RnaData &rna_dataA,
                         const RnaData &rna_dataB,
                         const ArcMatches &arc_matches,
                         const MatchProbs *match_probs,
                         const ScoringParams &params)
        : Scoring(seqA,seqB,rna_dataA,rna_dataB,
                  arc_matches,match_probs,params) {
 
        exp_indel_opening_score = boltzmann_weight(params.indel_opening_);
        exp_indel_opening_loop_score =
            boltzmann_weight(params.indel_opening_loop_);
        precompute_exp_sigma();
        precompute_exp_gapcost();
    }
 
    template <typename T>
    void
    PFScoring<T>::precompute_exp_sigma() {
        size_type lenA = seqA.length();
        size_type lenB = seqB.length();
 
        exp_sigma_tab.resize(lenA + 1, lenB + 1);
 
        for (size_type i = 1; i <= lenA; ++i) {
            for (size_type j = 1; j <= lenB; ++j) {
                exp_sigma_tab(i, j) = boltzmann_weight(sigma_tab(i, j));
            }
        }
    }
 
    template <typename T>
    void
    PFScoring<T>::precompute_exp_gapcost() {
        size_type lenA = seqA.length();
        size_type lenB = seqB.length();
 
        // resize and create tables
        exp_gapcost_tabA.resize(lenA + 1);
        exp_gapcost_tabB.resize(lenB + 1);
 
        for (size_type i = 1; i < lenA + 1; i++) {
            exp_gapcost_tabA[i] = boltzmann_weight(gapcost_tabA[i]);
        }
 
        for (size_type i = 1; i < lenB + 1; i++) {
            exp_gapcost_tabB[i] = boltzmann_weight(gapcost_tabB[i]);
        }
    }
 
} // end namespace LocARNA
 
#endif // LOCARNA_SCORING_HH