source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/extensions/mxscarna_src/BPPMatrix.hpp

//////////////////////////////////////////////////////////////
// BPPMatrix.hpp
//
// save the Base Pairing Probability Matrix for each sequences
//////////////////////////////////////////////////////////////

#ifndef __BBPMatrix_HPP__
#define __BBPMatrix_HPP__

#include <iostream>
#include <string>
#include "SparseMatrix.h"
#include "McCaskill.hpp"
#include "nrutil.h"



using namespace std;

namespace MXSCARNA{
class BPPMatrix {
private:

    int seqLength;       // sequence length;
    float cutOff;        // the threshold of probability 

    BPPMatrix() {}
public:
    SparseMatrix bppMat; // base pairing probability matrix 1-origin(1..seqLength)
    BPPMatrix(int bppmode, const string &seq, int seqLength, float inCutOff = 0.0001) : seqLength (seqLength), cutOff(inCutOff) {
      if( bppmode == 'r' ) // by katoh
        readBPPMatrix(seq);
      else if( bppmode == 'w' )
        setandwriteBPPMatrix(seq);
      else
        setBPPMatrix(seq);
    }
    BPPMatrix(int seqLength, float inCutOff, const Trimat<float> & tmpBppMat) : seqLength(seqLength), cutOff(inCutOff) {
      bppMat.SetSparseMatrix(seqLength, seqLength, tmpBppMat, cutOff);
    }

      
    void setBPPMatrix(const string &seq) {
        const char *tmpSeq = seq.c_str();
        McCaskill mc(seqLength, &tmpSeq[1]);
        mc.calcPartitionFunction();
        Trimat<float> tmpBppMat(seqLength + 1);

        for (int i = 0; i < seqLength; i++) {
          for(int j = i; j < seqLength; j++) {
            tmpBppMat.ref(i+1, j+1) = mc.getProb(i,j);
          }
        }
        bppMat.SetSparseMatrix(seqLength, seqLength, tmpBppMat, cutOff);
    }

    void setandwriteBPPMatrix(const string &seq) { // by katoh
        const char *tmpSeq = seq.c_str();
        McCaskill mc(seqLength, &tmpSeq[1]);
        mc.calcPartitionFunction();
        Trimat<float> tmpBppMat(seqLength + 1);
        float tmpbp;

        fprintf( stdout, ">\n" );
        for (int i = 0; i < seqLength; i++) {
          for(int j = i; j < seqLength; j++) {
        tmpbp = mc.getProb(i,j);
        if (tmpbp > 0.0001)
          fprintf( stdout, "%d %d %50.40f\n", i, j, tmpbp );
            tmpBppMat.ref(i+1, j+1) = tmpbp;
          }
        }
        bppMat.SetSparseMatrix(seqLength, seqLength, tmpBppMat, cutOff);
    }

    void readBPPMatrix(const string &seq) { // by katoh
        const char *tmpSeq = seq.c_str();
        char oneline[1000];
        int onechar;
        float prob;
        int posi, posj;
        static FILE *bppfp = NULL;


        if( bppfp == NULL )
        {
                bppfp = fopen( "_bpp", "r" );
                if( bppfp == NULL )
                {
                        fprintf( stderr, "Cannot open _bpp.\n" );
                        exit( 1 );
                }
        }

        Trimat<float> tmpBppMat(seqLength + 1);
        fgets( oneline, 999, bppfp );
        if( oneline[0] != '>' )
        {
                fprintf( stderr, "Format error\n" );
                exit( 1 );
        }
        while( 1 )
        {
                onechar = getc( bppfp ); 
                ungetc( onechar, bppfp );
                if( onechar == '>' || onechar == EOF ) break;

                fgets( oneline, 999, bppfp );
                sscanf( oneline, "%d %d %f", &posi, &posj, &prob );
//              fprintf( stderr, "%d %d -> %f\n", posi, posj, prob );
                tmpBppMat.ref(posi+1, posj+1) = prob;
        }

        bppMat.SetSparseMatrix(seqLength, seqLength, tmpBppMat, cutOff);
    }

    float GetProb(int i, int j) {
        return bppMat.GetValue(i,j);
    }

    float GetLength() const {
        return seqLength;
    }
    
    void updateBPPMatrix(const Trimat<float> &inbppMat) {
        bppMat.SetSparseMatrix(seqLength, seqLength, inbppMat, cutOff);
    }
};
}
#endif // __BPPMatrix_HPP__