src/LocARNA/scoring.hh

#ifndef LOCARNA_SCORING_HH
#define LOCARNA_SCORING_HH

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

#include <math.h>
#include <vector>

#include "aux.hh"

#include "scoring_fwd.hh"
#include "matrix.hh"

#ifndef NDEBUG
#include "sequence.hh"
#endif

#include "arc_matches.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 std::vector<pf_score_t> PFScoreVector;

    typedef Matrix<infty_score_t> ScoreMatrix;

    typedef Matrix<pf_score_t> PFScoreMatrix;

    typedef Matrix<double> ProbMatrix;

    class ScoringParams {
    public:
        const score_t basematch;

        const score_t basemismatch;

        const score_t indel;

        const score_t indel_loop;

        const score_t indel_opening;

        const score_t indel_opening_loop;

        const RibosumFreq *ribosum;

        const Ribofit *ribofit;

        const score_t unpaired_penalty;

        const score_t struct_weight;

        const score_t tau_factor;

        const score_t exclusion;

        const double exp_probA;

        const double exp_probB;

        const double temperature_alipf;

        const bool stacking;

        const bool new_stacking;

        const bool mea_scoring;

        const score_t alpha_factor;

        const score_t beta_factor;

        const score_t gamma_factor;

        const score_t probability_scale;

        /*
          The mea score is composed in the following way:

          params->probability_scale *
          (
          sum_basematchs (i,j)
          P(i~j)
          + alpha_factor/100 * (P(i unstructured) + P(j unstructured))
          +
          sum_arcmatchs (a,b)
          beta_factor/100 * (P(a)+P(b)) * P(al~bl) * P(ar~br) *
          ribosum_arcmatch(a,b)
          )
        */

    public:
        ScoringParams(score_t basematch_,
                      score_t basemismatch_,
                      score_t indel_,
                      score_t indel_loop_,
                      score_t indel_opening_,
                      score_t indel_opening_loop_,
                      RibosumFreq *ribosum_,
                      Ribofit *ribofit_,
                      score_t unpaired_penalty_,
                      score_t struct_weight_,
                      score_t tau_factor_,
                      score_t exclusion_,
                      double exp_probA_,
                      double exp_probB_,
                      double temp_,
                      bool stacking_,
                      bool new_stacking_,
                      bool mea_scoring_,
                      score_t alpha_factor_,
                      score_t beta_factor_,
                      score_t gamma_factor_,
                      score_t probability_scale_)
            : basematch(basematch_),
              basemismatch(basemismatch_),
              indel(indel_),
              indel_loop(indel_loop_),
              indel_opening(indel_opening_),
              indel_opening_loop(indel_opening_loop_),
              ribosum(ribosum_),
              ribofit(ribofit_),
              unpaired_penalty(unpaired_penalty_),
              struct_weight(struct_weight_),
              tau_factor(tau_factor_),
              exclusion(exclusion_),
              exp_probA(exp_probA_),
              exp_probB(exp_probB_),
              temperature_alipf(temp_),
              stacking(stacking_),
              new_stacking(new_stacking_),
              mea_scoring(mea_scoring_),
              alpha_factor(alpha_factor_),
              beta_factor(beta_factor_),
              gamma_factor(gamma_factor_),
              probability_scale(probability_scale_) {}
    };

    class Scoring {
    public:
        typedef BasePairs__Arc Arc; 

    private:
        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,
                bool exp_scores = false);

        void
        modify_by_parameter(score_t lambda);

        void
        apply_unpaired_penalty();

        score_t
        lambda() const {
            return lambda_;
        }

    private:
        // ------------------------------
        // 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

        // ------------------------------
        // 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; 

        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_exp_sigma();

        void
        precompute_gapcost();

        void
        precompute_exp_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;

        pf_score_t
        boltzmann_weight(score_t s) const {
            return exp(s / (pf_score_t)params->temperature_alipf);
        }

        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);
        }

        pf_score_t
        exp_basematch(size_type i, size_type j) const {
            return exp_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;

        template <bool gapAorB>
        score_t
        arcDel(const BasePairs__Arc &arc, bool stacked = false) const;

        pf_score_t
        exp_arcmatch(const ArcMatch &am) const {
            return boltzmann_weight(arcmatch(am));
        }

        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];
        }

        pf_score_t
        exp_gapA(size_type posA) const {
            assert(1 <= posA && posA <= seqA.length());
            return exp_gapcost_tabA[posA];
        }

        score_t
        gapB(size_type posB) const {
            assert(1 <= posB && posB <= seqB.length());

            return gapcost_tabB[posB];
        }

        pf_score_t
        exp_gapB(size_type posB) const {
            assert(1 <= posB && posB <= seqB.length());

            return exp_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;
        }

        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;
        }

        //
        // ------------------------------------------------------------

        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);
    }

} // end namespace LocARNA

#endif // LOCARNA_SCORING_HH