pairprobmap.cpp

//  This program shows how you may extend the computation by a
//  own driver derived from PartitionDriverBase.
//  The driver PairProbMap requires the forward/backward sums 
//  to determine the pair probabilities for each letter in the 
//  sequences.

#include <partitionworkspace.hpp>
#include <weights.hpp>
#include <posterprob.hpp>

using namespace std;

// This class is a driver that computes the 
// pair probabilities.
class PairProbMap : public PartitionDriverBase
{
public:

  // the resulting pair probabilities are 
  // written to this matrix 
  vector<vector<double> > pair_probs;


  // tell the main algorithm that this driver 
  // needs calls of the method OnForwardBackward 
  virtual bool NeedForwardBackward()   { return true; } 

  // Before the actual computation we allocate the required memory
  virtual void BeforeComputation() 
  {
    long length1(ws->Length1());
    long length2(ws->Length2());
    pair_probs.resize(length1,vector<double>(length2,0));
  }

  // the results are recieved from the PartitionWorkspace object 
  // for each row and stored in our own matrix
  virtual void OnForwardBackward(PartitionWorkspace::ForwardBackwardInfo & 
                                 fb_info)
  {
    long j;
    if(fb_info.i >= 0) 
      {
        for(long j=0; j < ws->Length2(); j++) 
          pair_probs[fb_info.i][j] = 
            fb_info.PairProb(fb_info.i,j);
      }
  }
  
  // Write data to stdout which can be displayed by a 3d plot program 
  // like gnuplot
  void WriteResult() 
  {

    long length1 = ws->Length1();
    long length2 = ws->Length2();
    long i,j;
    for(i = 0; i < length1; i++) 
      {
        for(j=0; j < length2; j++) 
          {
            cout << i << " " << j << " " << pair_probs[i][j] << endl;
          }
      }
  }

};




int main()
{
  // the sequences 
  ProteinSequence seq1("MTLWMHLLPLLPQLPLWGSNLAPEFAEQQGCFSHLNRTVELDVEVEYELCGERT");
  ProteinSequence seq2("MELWVRLLPLAPQLALYGPEFIDQHTSIVEALAAYCGGHLGRAICSNLEVEVLYEVCEGRS");

  // alignment parameters 
  double temperature   = 1.0;
  string matrix        = "blosum62";
  int    gap_open      = 11;
  int    gap_extension = 1;
  bool   allow_DI      = true;
  bool   allow_ID      = false;
  bool   do_hmm        = true;

  // weights for the model (BoltzmannWeight or HmmWeight) 
  Weight * weight; 
  if(do_hmm) 
    {
      // pair hidden Markov model 
      weight = new HmmWeight(proteinAlphabet,
                             proteinAlphabet,
                             matrix,
                             gap_open,
                             gap_extension,
                             GlobalBounded,
                             allow_DI,
                             allow_ID);
    }
  else 
    {
      // construct the weights for the finite temperature model 
      weight = new BoltzmannWeight(proteinAlphabet,
                                   proteinAlphabet,
                                   matrix,
                                   temperature,
                                   gap_open,
                                   gap_extension,
                                   GlobalBounded,
                                   allow_DI,
                                   allow_ID);
    }

  // the workspace object 
  PartitionWorkspace pws(10000);

  // our own driver 
  PairProbMap        pair_prob;

  // add driver and perform the computation
  pws.AddDriver(pair_prob);
  pws.PartitionFunction(seq1,
                        seq2,
                        *weight);
  // write the result 
  pair_prob.WriteResult();
  return 0;
}