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

namespace LocARNA {


    //#define MEA_SCORING_OLD

    
    class RibosumFreq;
    class Ribofit;
    class Scoring;
    class Sequence;
    class BasePairs;
    class BasePairs__Arc;
    class ArcMatches;
    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;
    
    
        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(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); }


        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;


        score_t  
        arcDel(const BasePairs__Arc &arc, bool gapAorB, 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);
        }

        score_t gapX(size_type alignedToGap, bool gapInA) const {
            if (gapInA)
                return gapA(alignedToGap);
            else
                return gapB(alignedToGap);
        }

        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

} // end namespace LocARNA

#endif // LOCARNA_SCORING_HH