ribosum.hh

#ifndef LOCARNA_RIBOSUM_HH
#define LOCARNA_RIBOSUM_HH
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
#include <array>
#include <cstdlib>
#include <fstream>
#include <sstream>
#include <math.h>
#include <string>
 
#include "alphabet.hh"
#include "matrix.hh"
 
namespace LocARNA {
 
    class Ribosum {
    public:
        typedef Matrix<double> matrix_t; 
        typedef Alphabet<std::string, 4> base_alphabet_type; 
        typedef Alphabet<char, 4> char_alphabet_type; 
        typedef Alphabet<std::string, 16> arc_alphabet_type; 
 
        explicit Ribosum(const std::string &filename);
 
        virtual ~Ribosum();
 
        const matrix_t &
        get_basematch_scores() const {
            return bm_;
        }
 
        const matrix_t &
        get_arcmatch_scores() const {
            return am_;
        }
 
        const base_alphabet_type &
        string_alphabet() const {
            return basename_alphabet_;
        }
 
        const char_alphabet_type &
        alphabet() const {
            return char_basename_alphabet_;
        }
 
        const std::string &
        get_name() const {
            return name_;
        }
 
        double
        basematch_score(char i, char j) const {
            return bm_(alphabet().idx(i), alphabet().idx(j));
        }
 
        double
        arcmatch_score(char i, char j, char k, char l) const {
            return am_(alphabet().idx(i) * 4 + alphabet().idx(j),
                       alphabet().idx(k) * 4 + alphabet().idx(l));
        }
 
        friend std::ostream &
        operator<<(std::ostream &out, const Ribosum &ribosum);
    protected:
        std::string name_; 
        matrix_t bm_;      
        matrix_t am_;      
 
        base_alphabet_type basename_alphabet_; 
        arc_alphabet_type arcname_alphabet_;  
 
        char_alphabet_type
            char_basename_alphabet_; 
 
        template<size_t N>
        std::istream &
        read_matrix(std::istream &in,
                    matrix_t &mat,
                    const Alphabet<std::string, N> &names) const;
 
        template<size_t N>
        std::ostream &
        write_matrix(std::ostream &out,
                     const matrix_t &mat,
                     const Alphabet<std::string, N> &alph) const;
 
        Ribosum();
 
        void
        read_ribosum(std::istream &in);
 
        char_alphabet_type
        make_char_alphabet() const;
 
        void
        set_basename_alphabet(const std::array<std::string,4> &a) {
            basename_alphabet_ = base_alphabet_type(a);
            char_basename_alphabet_ = make_char_alphabet();
        }
 
        void
        set_arcname_alphabet(const std::array<std::string,16> &a) {
            arcname_alphabet_ = arc_alphabet_type(a);
        }
    };
 
    class RibosumFreq : public Ribosum {
    public:
        explicit RibosumFreq(const std::string &filename);
 
        double
        base_prob(char i) const {
            return base_probs_(alphabet().idx(i), 0);
        }
 
        double
        base_nonstruct_prob(char i) const {
            return base_nonstruct_probs_(alphabet().idx(i), 0);
        }
 
        const matrix_t &
        get_base_probs() const {
            return base_probs_;
        }
 
        const matrix_t &
        get_base_nonstruct_probs() const {
            return base_nonstruct_probs_;
        }
 
        double
        basepair_prob(char i, char j) const {
            return basepair_probs_(alphabet().idx(i), alphabet().idx(j));
        }
 
        const matrix_t &
        get_basepair_probs() const {
            return basepair_probs_;
        }
 
        double
        basematch_prob(char i, char j) const {
            return basematch_probs_(alphabet().idx(i), alphabet().idx(j));
        }
 
        const matrix_t &
        get_basematch_probs() const {
            return basematch_probs_;
        }
 
        double
        arcmatch_prob(char i, char j, char k, char l) const {
            return arcmatch_probs_(alphabet().idx(i) * 4 + alphabet().idx(j),
                                   alphabet().idx(k) * 4 + alphabet().idx(l));
        }
 
        const matrix_t &
        get_arcmatch_probs() const {
            return arcmatch_probs_;
        }
 
        double
        base_unpaired_prob(char i) const;
 
        double
        basematch_score_corrected(char i, char j) const;
 
        void
        print_basematch_scores_corrected(std::ostream &out) const;
 
        void
        read_matrix(std::istream &in,
                    const std::string &header,
                    matrix_t &mat,
                    size_t xdim,
                    size_t ydim);
 
        void
        write_ICC_code(std::ostream &out, const std::string &ribname) const;
 
        std::ostream &
        write_matrix(std::ostream &out,
                     const std::string &name,
                     const Matrix<double> &mat) const;
 
        friend std::ostream &
        operator<<(std::ostream &out, const RibosumFreq &ribosum);
 
    protected:
        RibosumFreq();
 
        matrix_t base_probs_;           
        matrix_t base_nonstruct_probs_; 
        matrix_t basepair_probs_;       
        matrix_t basematch_probs_;      
        matrix_t arcmatch_probs_;       
 
        void
        write_CC_matrix(std::ostream &out,
                        const std::string &ribname,
                        const std::string &matname,
                        int x,
                        int y,
                        const Ribosum::matrix_t &m) const;
 
    private:
        void
        read_frequencies(std::istream &in);
    };
 
 
    template <size_t N>
    std::istream &
    Ribosum::read_matrix(std::istream &in,
                         Matrix<double> &mat,
                         const Alphabet<std::string, N> &alph) const {
        auto siz = alph.size();
 
        std::string line;
 
        while (std::getline(in, line) && line == "")
            ;
 
        {
            std::istringstream linestream(line);
 
            for (size_t i = 0; i < siz; i++) {
                std::string name;
                linestream >> name;
 
                if (name != alph[i])
                    throw(std::ifstream::failure(
                        "Expecting correct table header. Found: " + line));
            }
        }
 
        mat.resize(siz, siz);
        for (size_t i = 0; i < siz; i++) {
            std::getline(in, line);
            std::istringstream linestream(line);
            std::string base;
            linestream >> base;
            if (base != alph[i])
                throw(std::ifstream::failure("Expecting base name " +
                                             alph[i] + " as row header"));
 
            for (size_t j = 0; j <= i; j++) {
                double number;
                linestream >> number;
                mat(i, j) = mat(j, i) = number;
            }
        }
        return in;
    }
 
    template<size_t N>
    std::ostream &
    Ribosum::write_matrix(std::ostream &out,
                          const Matrix<double> &mat,
                          const Alphabet<std::string, N> &alph) const {
        out << alph << std::endl;
        out << mat << std::endl;
        return out;
    }
 
 
} // end namespace LocARNA
 
#endif // LOCARNA_RIBOSUM_HH