text: /*
    Copyright by Matthias Hoechsmann (C) 2002-2004
    =====================================                                   
    You may use, copy and distribute this file freely as long as you
    - do not change the file,
    - leave this copyright notice in the file,
    - do not make any profit with the distribution of this file
    - give credit where credit is due
    You are not allowed to copy or distribute this file otherwise
    The commercial usage and distribution of this file is prohibited
    Please report bugs and suggestions to <mhoechsm@TechFak.Uni-Bielefeld.DE>
  */
  
  #ifndef WIN32
  #include "config.h"
  #endif
  
  #include <cmath>
  
  #ifdef HAVE_LIBG2
  #include <g2.h>
  #include <g2_PS.h>
  #endif
  
  #include <fstream>
  #include <iostream>
  #include <iomanip>
  #include <sstream>
  #include <string.h>
  
  #ifdef HAVE_LIBXMLPLUSPLUS
  #include <libxml++/libxml++.h>
  #endif
  
  #include "matrix.h"
  #include "rna_profile_alignment.h"
  #include "rnafuncs.h"
  #include "utils.h"
  
  #ifdef HAVE_LIBRNA  // This features require the ViennaRNA library
  #include "part_func.h"
  #include "fold_vars.h"
  
  extern "C"{
  #include "fold.h"
  }
  
  #endif
  
  #include "ppforest.t.cpp"
  
  ostream& operator <<(ostream &s,RNA_Alphabet_Profile p)
  {
      
    s.precision(2);
    s  << p.p[ALPHA_PRO_BASE_A] << "\\n";
    s.precision(2);
    s << p.p[ALPHA_PRO_BASE_C] << "\\n";
    s.precision(2);
    s << p.p[ALPHA_PRO_BASE_G] << "\\n";
    s.precision(2);
    s << p.p[ALPHA_PRO_BASE_U] << "\\n";
    s.precision(2);
    s << p.p[ALPHA_PRO_BASEPAIR] << "\\n";
    s.precision(2);
    s << p.p[ALPHA_PRO_GAP];
      s.precision(2);
      s << p.p[ALPHA_PRO_BASE];
  
    //  for(Uint i=0;i<p.columnStr.length();i++)
    //    s << p.columnStr[i] << "\\n";
  
    return s;
  }
  
  /* ****************************************** */
  /*    Constructor and Destructor functions    */
  /* ****************************************** */
  
  RNAProfileAlignment::RNAProfileAlignment(const string &baseStr, const string &viennaStr, const string &name)
  : PPForestAli<RNA_Alphabet_Profile,RNA_Alphabet_Profile>(RNAFuncs::treeSize(viennaStr)),
    m_name(name),
    m_numStructures(1)
  {
  	buildForest(baseStr,viennaStr);
  	addStrName(name);
  };
  
  #ifdef HAVE_LIBRNA  // This features require the ViennaRNA library
  RNAProfileAlignment::RNAProfileAlignment(const string &baseStr, const string &name, const string &constraint, double t)
    : PPForestAli<RNA_Alphabet_Profile,RNA_Alphabet_Profile>(2*baseStr.length()),
      m_name(name),
      m_numStructures(1)
  {
    char *viennaStr=NULL;
    
    // calculate partition function for the sequence
    do_backtrack=1;
    init_pf_fold(baseStr.length());
  
    //if(constraint.length()>0)
    //pf_fold((char*)baseStr.c_str(),(char*)constraint.c_str());    // expicit conversion to non-const value, but pf_fold does not alter baseStr
    //else
    pf_fold((char*)baseStr.c_str(),NULL);    // expicit conversion to non-const value, but pf_fold does not alter baseStr
  
    viennaStr=new char[baseStr.length()+1];
    dangles=2;
    fold((char*)baseStr.c_str(),viennaStr);
  
    setSize(RNAFuncs::treeSize(viennaStr));
    buildForest(baseStr,viennaStr,true);
    
    free_pf_arrays();
    delete[] viennaStr;
    
    //  hasSequence=true;
    addStrName(name);
  }
  #endif
  
  RNAProfileAlignment::RNAProfileAlignment(const string &filename)
  {
    // read ppforest from file
    size_type size;
    ifstream s;
    char *colStr;
  
    s.open(filename.c_str());
  
    // first of all save the size
    s.read((char*)&size,sizeof(size_type));
    s.read((char*)&m_numStructures,sizeof(Uint));
    
    colStr=new char[m_numStructures+1];
    colStr[m_numStructures]=0;   // terminate string
    
    // initialize structures
    //  ::PPForestBase(m_size);  
    //  m_lb=new L[ppf.size()];
    initialize(size);
  
    // save the arrays
      for(Uint i=0;i<size;i++)
        {
  	label_type l;
  	
  	for(int r=0;r<RNA_ALPHABET_SIZE;r++)
  	  {
  	    double d;
  	    s.read(reinterpret_cast<char*>(&d),sizeof(double));
  	    l.p[r]=d;
  	  }
  
  	s.read(colStr,sizeof(char)*m_numStructures);
  
  	l.columnStr=colStr;
  	m_lb[i]=l;
      }
  
    s.read(reinterpret_cast<char*>(m_rb),sizeof(size_type)*size);
    s.read(reinterpret_cast<char*>(m_noc),sizeof(size_type)*size);
    //  s.read((char*)m_sumUpCSF,sizeof(size_type)*size);
    //  s.read((char*)m_rmb,sizeof(size_type)*size);    
  
    for(Uint r=0;r<m_numStructures;r++)
      {
        char str[21];
        str[20]=0;
        
        s.read(str,sizeof(char)*20);
        m_strNames.push_back(string(str));
      }
  
  
    calcSumUpCSF();
    calcRMB();  
  }
  
  /*
  RNAProfileForest::RNAProfileForest(const Profile_RNA_Alignment_Forest &ppf)
    : RNAForestBase<RNA_Alphabet_Profile>(ppf)
  {
    m_numStructures=ppf.m_numStructuresX+ppf.m_numStructuresY;
  }*/
  
  /* ****************************************** */
  /*            Private functions               */
  /* ****************************************** */
  
  void RNAProfileAlignment::makeRepLabel(size_type node, RNA_Alphabet_Profile a,  RNA_Alphabet_Profile b)
      {
  	  double p,q;
  	  double m;
  
  	  p=m_numStructuresX;	// weight number of structures in profile
  	  q=m_numStructuresY;
  	  m=p+q;		
  	  p/=m;
  	  q/=m;
  
  	  // profile
        for(int i=0;i<RNA_ALPHABET_SIZE;i++)
  		m_lb[node].p[i]=p*a.p[i]+q*b.p[i];
  
  	  // alignment column
  	  m_lb[node].columnStr=a.columnStr + b.columnStr;
      };
    
  void RNAProfileAlignment::makeDelLabel(size_type node, RNA_Alphabet_Profile a)
      {
  	  double p,q;
  	  double m;
  
  	  p=m_numStructuresX;	// weight number of structures in profile
  	  q=1;
  	  m=p+q;		
  	  p/=m;
  	  q/=m;
  
  	  // profile
  	  m_lb[node]=a;
  	  //m_lb[node].p[ALPHA_GAP]=(1+m_lb[node].p[ALPHA_GAP])/2.0;	
  
  	  for(int i=0;i<RNA_ALPHABET_SIZE;i++)
  	    {
  	      m_lb[node].p[i]*=p;
  	    }
  	  m_lb[node].p[ALPHA_PRO_GAP]+=q;
  
  	  // alignment column
  	  m_lb[node].columnStr=a.columnStr + string(m_numStructuresY,ALPHA_GAP);
      };
    
  void RNAProfileAlignment::makeInsLabel(size_type node, RNA_Alphabet_Profile b)
      {
  	  double p,q;
  	  double m;
  
  	  p=1;						// weight number of structures in profile
  	  q=m_numStructuresY;
  	  m=p+q;		
  	  p/=m;
  	  q/=m;
  
  	  // profile
  	  m_lb[node]=b;
  	  //m_lb[node].p[ALPHA_GAP]=(1+m_lb[node].p[ALPHA_GAP])/2.0;	
  
  	  for(int i=0;i<RNA_ALPHABET_SIZE;i++)
  	    {
  	      m_lb[node].p[i]*=q;
  	    }
  	  m_lb[node].p[ALPHA_PRO_GAP]+=p;
  	  
  	  //	  m_lb[node].p[ALPHA_PRO_GAP]=(p+q*m_lb[node].p[ALPHA_PRO_GAP])/2.0;
  
  	  // alignment column
  	  m_lb[node].columnStr=string(m_numStructuresX,ALPHA_GAP) + b.columnStr;
      };
    
  void RNAProfileAlignment::showLabel(ostream &s,RNA_Alphabet_Profile p) const
    {
      s << p;
    };
  
  
  void RNAProfileAlignment::makeLabel(RNA_Alphabet_Profile &p,char c)
  { 
    int i;
  
    // initialize profile entries to zero
    for(i=0;i<RNA_ALPHABET_SIZE;i++)
      p.p[i]=0.0;
  
    i=alpha2RNA_Alpha(c);
    p.p[i]=1.0;
  
    // if it is acgu it is also a base 
    if(i<=ALPHA_PRO_BASE_U)
      p.p[ALPHA_PRO_BASE]=1.0;
  
    // set column string
    p.columnStr=c;
  }
  
  void RNAProfileAlignment::buildForest(const string &baseStr, const string &viennaStr, bool use_bp_prob)
  {
  	Ulong basePairCount=0,maxDepth=0,stackPtr;
  	Uint baseStrLen,viennaStrLen,node,*nodeStack, *numChildrenStack, *baseposStack;
  
  	assert(RNAFuncs::isRNAString(baseStr));
  
  	RNAFuncs::isViennaString(viennaStr,basePairCount,maxDepth);
  
  	baseStrLen=baseStr.length();
  	viennaStrLen=viennaStr.length();
  
  	if(baseStrLen)
  		hasSequence=true;
  	else
  		hasSequence=false;	
  
  	// check if base string and vienna string have the same length
  	//  if(baseStr.length() > 0)
  	//    if(baseStr.length() != viennaStrLen)
  	//      throw RNAForestExceptionInput(RNAForestExceptionInput::Error_BaseStringAndViennaStringIncompatible);
  
  	nodeStack=new Uint[maxDepth+1];
  	numChildrenStack=new Uint[maxDepth+1];
  	baseposStack=new Uint[maxDepth+1];
  	memset(nodeStack,0,sizeof(Uint)*maxDepth+1);
  	memset(numChildrenStack,0,sizeof(Uint)*maxDepth+1);
  	memset(baseposStack,0,sizeof(Uint)*maxDepth+1);
  
  	// fill PPForest structure
  	stackPtr=0;
  	node=0;
  	for(Uint i=0;i<viennaStrLen;i++)
  	{		
  		switch(viennaStr[i])
  		{
  		case '.':
  			// set label
  			if(baseStrLen)
  				makeLabel(m_lb[node],baseStr[i]);
  			else
  				makeLabel(m_lb[node],'B');				
  
  			// set right brother
  			if(node==size()-1)
  				setRightBrotherIndex(node,0);
  			else
  				setRightBrotherIndex(node,node+1);
  
  			// set num children
  			setNumChildren(node,0);
  
  			// increase stack values
  			numChildrenStack[stackPtr]++;
  
  			node++;
  			break;
  
  		case '(':
  			// set label
  			makeLabel(m_lb[node],'P');			
  
  			// increase stack values
  			numChildrenStack[stackPtr]++;
  			baseposStack[stackPtr]=i+1;
  
  			// push 
  			stackPtr++;
  			nodeStack[stackPtr]=node;
  			numChildrenStack[stackPtr]=1;
  
  			node++;
  
  			// set label
  			if(baseStrLen)
  				makeLabel(m_lb[node],baseStr[i]);
  			else
  				makeLabel(m_lb[node],'B');
  
  			// set right brother
  			setRightBrotherIndex(node,node+1);
  
  			// set num children
  			setNumChildren(node,0);
  
  			node++;				
  			break;
  		case ')':
  			// set label
  			if(baseStrLen)
  				makeLabel(m_lb[node],baseStr[i]);
  			else
  				makeLabel(m_lb[node],'B');  
  
  			// set right brother
  			setRightBrotherIndex(node,0);    
  
  			// set num children
  			setNumChildren(node,0);  
  
  			// pop
  			if(node==size()-1)
  				setRightBrotherIndex(nodeStack[stackPtr],0);
  			else
  				setRightBrotherIndex(nodeStack[stackPtr],node+1);  
  
  			setNumChildren(nodeStack[stackPtr],numChildrenStack[stackPtr]+1);
  			
  			stackPtr--;
  
  #ifdef HAVE_LIBRNA			
  			// set basepair probability
  			if(use_bp_prob)
  			  {			    
  			    m_lb[nodeStack[stackPtr]].p[ALPHA_PRO_BASEPAIR]=pr[iindx[baseposStack[stackPtr]]-(i+1)];
  			    m_lb[nodeStack[stackPtr]].p[ALPHA_PRO_GAP]=1-m_lb[nodeStack[stackPtr]].p[ALPHA_PRO_BASEPAIR];
  			  }			
  #endif
  
  			
  			node++;
  			break;							
  		}		
  	}	
  	
  	delete[] nodeStack;
  	delete[] numChildrenStack;
  	delete[] baseposStack;
  	calcSumUpCSF();
  
  	calcRMB();
  	
  	assert (m_noc[m_size-1]==0);
  }
  
  void RNAProfileAlignment::getStructureAlignmentFromCSF(string &s, deque<double> &pairprob, double t,size_type i, size_type j) const
  {
    size_type h;
    double bestPairScore;
    //  list<Uint> leftPairList;
    //list<Uint> rightPairList;
    //Uint bestLeftIndex,bestRightIndex;
    //Uint lastLeftIndex,lastRightIndex;
    
    //  QWATCH(i);
    //  QWATCH(j);	
  
    if(j==0)
      return;
    
    if(isPair(i))
      {
        // TRACE(DBG_GET_PROFILE_STRUCTURE,"Profile_RNA_Alignment_Forest::getStructureAlignmentFromCSF","basepair");
        // WATCH(DBG_GET_PROFILE_STRUCTURE,"Profile_RNA_Alignment_Forest::getStructureAlignmentFromCSF",m_lb[i].p[ALPHA_BASEPAIR]);
  
        bestPairScore=bestPairs(i);
  
        // backtrack best pairs
        if(bestPairScore == bestPairs(i+1) + bestPairs(getRightmostBrotherIndex(i+1)) || m_lb[i].p[ALPHA_PRO_BASEPAIR] < t)
        {
  	// i pairs not
  	//	cout << "unpaired" << endl;
  	getStructureAlignmentFromCSF(s,pairprob,t,i+1,noc(i));
  	//	cout << "back to:" << endl;
  	//	QWATCH(i);
  	//	QWATCH(j);
        }
        else
        {
  	//	cout << "paired" << endl;
  	// i pairs
  	s += '(';
  	pairprob.push_back(m_lb[i].p[ALPHA_PRO_BASEPAIR]);
  	
  	// left path - righthand best pairs
  	h=i+1;
  	while(h < size() && isPair(h))
  	{
  	  //	  cout << "left" << endl;
  	  //	  QWATCH(h);
  
  	  assert((int)noc(h)-1>=0);
  	  getStructureAlignmentFromCSF(s,pairprob,t,rb(h+1),noc(h)-1);
  	  h=h+1;
  	}
  
  
  	assert((int)noc(i)-2>=0);
  	getStructureAlignmentFromCSF(s,pairprob,t,rb(i+1),noc(i)-2);
  	//	cout << "back to:" << endl;
  	//	QWATCH(i);
  	//	QWATCH(j);
  
  	// right path - lefthand best pairs
  	h=getRightmostBrotherIndex(i+1);
  	while(h < size() && isPair(h))
  	{
  	  //	  cout << "right" << endl;
  	  //	  QWATCH(h);
  
  	  assert((int)noc(h)-1>=0);
  	  getStructureAlignmentFromCSF(s,pairprob,t,h+1,noc(h)-1);
  	  //h=h+1;
  	  h=getRightmostBrotherIndex(h+1);
  	}      
   
  
  	s += ')';
  	pairprob.push_back(m_lb[i].p[ALPHA_PRO_BASEPAIR]);
        }
      }
    else
      {
        s+= '.';
        pairprob.push_back(m_lb[i].p[ALPHA_PRO_BASE]);
      }
    
    // right forest
    getStructureAlignmentFromCSF(s,pairprob,t,rb(i),j-1);
  }
  
  double RNAProfileAlignment::bestPairs(size_type node) const
  {
    size_type i=node;
    double d1,d2;
  
    WATCH(DBG_GET_PROFILE_STRUCTURE,"RNAProfileForest::getStructureAlignment",node);
  
    if(isBase(node))
      return 0;
    else
      {
        // node pairs not
        d1=bestPairs(node+1) + bestPairs(getRightmostBrotherIndex(node+1));
  
        // node pairs
        d2=label(node).p[ALPHA_PRO_BASEPAIR];
  
        // left path - righthand best pairs
        i=node+1;
        while(i < size() && isPair(i))
  	{
  	  d2+=bestPairs(getRightmostBrotherIndex(i+1));
  	  i=i+1;
  	}
  
        // right path - lefthand best pairs
        i=getRightmostBrotherIndex(node+1);
        while(isPair(i) && i < size())
  	{
  	  d2+=bestPairs(i+1);
  	  i=getRightmostBrotherIndex(i+1);
  	}
  
        return max(d1,d2);
      }
  } 
  
  #ifdef HAVE_LIBG2  // This features require the g2 library
  void RNAProfileAlignment::drawBaseCircles(int device_id,const BaseProbs &bp,double center_x,double center_y) const
  {
    const double box_size=6.0;
    const double max_radius=3.0;
  
    double xpos,ypos;
    int color;
    //  double dashes=0.5;
  
    // draw the base probabilities as circles on the edges of the square and the gap probability
    // as the center square
  
    // upper left corner = a
    xpos=center_x-box_size/2.0;
    ypos=center_y+box_size/2.0;
    
    color=g2_ink(device_id,1,0,0);
    g2_pen(device_id,color);
    g2_filled_circle(device_id,xpos,ypos,max_radius*bp.a);
  
    // upper right corner = c
    xpos=center_x+box_size/2.0;
    ypos=center_y+box_size/2.0;
    
    color=g2_ink(device_id,0,1,0);
    g2_pen(device_id,color);
    g2_filled_circle(device_id,xpos,ypos,max_radius*bp.c);
  
    // lower right corner = g
    xpos=center_x+box_size/2.0;
    ypos=center_y-box_size/2.0;
    
    color=g2_ink(device_id,0,0,1);
    g2_pen(device_id,color);
    g2_filled_circle(device_id,xpos,ypos,max_radius*bp.g);
  
    // lower right corner = u
    xpos=center_x-box_size/2.0;
    ypos=center_y-box_size/2.0;
    
    color=g2_ink(device_id,1,0,1);
    g2_pen(device_id,color);
    g2_filled_circle(device_id,xpos,ypos,max_radius*bp.u);
  
    // gap in the center  
    color=g2_ink(device_id,0,0,0);
    g2_pen(device_id,color);
    g2_filled_circle(device_id,center_x,center_y,max_radius*bp.gap);  
  
    // draw rectangle for orientation
    //g2_set_dash(device_id,1,&dashes);
    color=g2_ink(device_id,0.75,0.75,0.75);
    g2_pen(device_id,color);
    g2_rectangle(device_id,center_x-box_size/2,center_y-box_size/2,center_x+box_size/2,center_y+box_size/2);
  }
  #endif
  
  inline double RNAProfileAlignment::getMlBaseFreq(const BaseProbs &bp) const
  {
    return max(max(max(bp.a,bp.c),bp.g),bp.u);
  }
  
  char RNAProfileAlignment::getMlBase(const BaseProbs &bp) const
  {
    double p = getMlBaseFreq(bp);
    if(bp.a==p)
      return ALPHA_BASE_A;
    if(bp.c==p)
      return ALPHA_BASE_C;
    if(bp.g==p)
      return ALPHA_BASE_G;
    if(bp.u==p)
      return ALPHA_BASE_U;
  
    assert(false); // you should never get so far ...
    return 'x';  
  }
  
  void RNAProfileAlignment::getSeqAli(string &seq,Uint row, Uint i, Uint j) const
  {
  	if(i==0 && j==getMaxLength(0))
  		seq="";
  
  	if(j==0)
  		return;
  
  	// basepair=internal node
  	if(isPair(i))
  	{
  		getSeqAli(seq,row,i+1,noc(i));
  		getSeqAli(seq,row,rb(i),j-1);		
  	}
  	else
  	{
  		// base=leaf
  		seq+=m_lb[i].columnStr[row];
  		getSeqAli(seq,row,rb(i),j-1);
  	}
  }
  
  void RNAProfileAlignment::getStructAli(string &s,Uint row) const
  {
    s="";
  
    map<Uint,Uint> pairs;
    makePairTable(pairs, row);
  
    // iterate through leaves nodes and use information of pairs
    for(Uint i=0;i<m_size;i++)
      {
  	  RNA_Alphabet c;
  	  
  	  c=m_lb[i].columnStr[row];
  
        if(isBase(i))
  	  {
  		if(c != '-')
  		{
  			if(pairs.find(i) != pairs.end())	// is base paired ?
  			{
  				Uint j=pairs[i];
  				if(i<j)
  					s+='(';
  				else
  					s+=')';
  			}
  			else
  				s+='.';
  		}
  		else
  			s+='-';
  	  }
      }
  }
  
  void RNAProfileAlignment::makePairTable(map<Uint,Uint> &pairs, Uint row) const
  {
  	pair<int,int> *baseIndex;		// left and right base
  	RNA_Alphabet c;
  
  	assert(row<m_numStructures);
  	baseIndex=new pair<int,int>[m_size];
  
  	// initialize pairs, all gaps
  	for(int i=size()-1;i>=0;i--)
  	  {
  	    baseIndex[i].first=-1;
  	    baseIndex[i].second=-1;
  	  }
  
  	for(int i=size()-1;i>=0;i--)
  	{
  	        c=m_lb[i].columnStr[row];
  
  		if(isLeave(i))
  		  {
  			if(c!=ALPHA_GAP)
  			{
  			    baseIndex[i].first=i;
  			    baseIndex[i].second=i;
  			}
  		  }
  		else
  		{
  			// internal node
  			// leftmost and rightmost base
  			bool lmBaseFound=false;
  			for(size_type r=0,h=i+1;r<getMaxLength(i+1);r++,h=rb(h))
  			{					  
  				// leftmost base
  				if(!lmBaseFound && baseIndex[h].first != -1)
  				{
  					baseIndex[i].first=baseIndex[h].first;
  					lmBaseFound=true;
  				}
  
  				// rightmost base
  				if(baseIndex[h].second != -1)
  				{
  					baseIndex[i].second=baseIndex[h].second;
  				}
  			}
  
  			// report pairing bases if P node
  			if(c==ALPHA_BASEPAIR)
  			{
  			        assert(baseIndex[i].first != -1);
  			        assert(baseIndex[i].second != -1);
  				assert(baseIndex[i].first < baseIndex[i].second);
  
  				pairs[baseIndex[i].first]=baseIndex[i].second;
  				pairs[baseIndex[i].second]=baseIndex[i].first;
  			}
  		}
  	}
  
  	delete[] baseIndex;       
  }
  
  void RNAProfileAlignment::filterConsensus(string &structure, deque<double> &pairprob, deque<BaseProbs> &baseprobs, double minFreq) const
  {
  	deque<double>::iterator itPP;
  	//	deque<BaseProbs>::iterator itBP;
  	string::iterator it;
  
  	for(itPP=pairprob.begin(),it=structure.begin();itPP!=pairprob.end();itPP++,it++)
  	{
  		if(*itPP<minFreq)
  		{
  			*itPP=0.0;
  			*it='.';
  		}
  	}
  
  	/*
  	for(it=structure.begin(),itPP=pairprob.begin(),itBP=baseprobs.begin();it!=structure.end();)
  	{
  		if(itBP->base<minFreq)
  		{				
  			structure.erase(it);
  			baseprobs.erase(itBP);
  			pairprob.erase(itPP);
  
  			it=structure.begin();
  			itPP=pairprob.begin();
  			itBP=baseprobs.begin();			
  		}
  		else
  		{
  			it++;
  			itPP++;
  			itBP++;
  		}
  	}
  	*/
  }
  
  /* ****************************************** */
  /*             Public functions               */
  /* ****************************************** */
  
  void RNAProfileAlignment::getSequenceAlignment(deque<BaseProbs> &baseprob) const
  {
  	BaseProbs bp;
  
  	// generate base strings
  	for(size_type i=0;i<size();i++)
  	{
  		if(isBase(i))
  		{	
  			bp.a=label(i).p[ALPHA_PRO_BASE_A];
  			bp.c=label(i).p[ALPHA_PRO_BASE_C];
  			bp.g=label(i).p[ALPHA_PRO_BASE_G];
  			bp.u=label(i).p[ALPHA_PRO_BASE_U];
  			bp.gap=label(i).p[ALPHA_PRO_GAP];
  			bp.base=bp.a+bp.c+bp.g+bp.u;
  
  			baseprob.push_back(bp);
  		}
  	} 
  }
  
  void RNAProfileAlignment::getStructureAlignment(double t,string &s, deque<double> &pairprob) const
  {
  	WATCH(DBG_GET_PROFILE_STRUCTURE,"Profile_RNA_Alignment_Forest::getStructureAlignment",size());
  
  	s="";
  	getStructureAlignmentFromCSF(s,pairprob,t,0,getMaxLength(0));
  }
  
  #ifdef HAVE_LIBG2  // This features require the g2 library
  void RNAProfileAlignment::squigglePlot(const string &filename, SquigglePlotOptions &options) const
  {
  	const double base_fontsize=8;
  	const Uint num_grey_colors=100;
  	const double min_grey_color=1.0;
  
  	string seq,structure;
  	string base,structname;
  	float *X,*Y,min_X=0,max_X=0,min_Y=0,max_Y=0;
  	Uint i;
  	short *pair_table;
  	int id_PS,id;
  	int ps_grey_colors[num_grey_colors];
  	int ps_color_red;
  	int ps_color_black;
  	double xpos,ypos;
  
  	deque<double> pairprob;
  	deque<BaseProbs> baseprobs;
  
  	getStructureAlignment(options.minPairProb,structure,pairprob);
  	getSequenceAlignment(baseprobs);
  
  	//  filterConsensus(structure,pairprob,baseprobs,0.5);
  
  	//assert(baseprobs.size() == structure.size());
  	if(baseprobs.size() != structure.size())
  		cerr <<  "Error in resolving consensus structure!" << endl;
  
  	X = new float[structure.size()];
  	Y = new float[structure.size()];
  
  	pair_table = make_pair_table(structure.c_str());
  	i = naview_xy_coordinates(pair_table, X, Y);
  	if(i!=structure.size())
  		cerr << "strange things happening in squigglePlot ..." << endl;
  
  	// calculate image dimesions
  	for(i=0;i<structure.size();i++)
  	{
  		min_X=min(min_X,X[i]);
  		max_X=max(max_X,X[i]);
  		min_Y=min(min_Y,Y[i]);
  		max_Y=max(max_Y,Y[i]);
  	}
  
  	//  id_PS  = g2_open_PS("ali.ps", g2_A4, g2_PS_port);
  	id_PS  = g2_open_EPSF((char*)filename.c_str());
  	id     = g2_open_vd();
  	g2_attach(id,id_PS);
  
  	//  cout << "min_X: " << min_X <<",max_X: " << max_X << ",min_Y: " << min_Y << "max_Y: " << max_Y << endl; 
  	g2_set_coordinate_system(id_PS,595/2.0,842/2.0,0.5,0.5);
  	g2_set_line_width(id,0.2);
  
  
  	// set colors
  	double intv=min_grey_color/(double)num_grey_colors;
  	for(i=0;i<num_grey_colors;i++)
  	{
  		double grey_color=min_grey_color-i*intv;
  		ps_grey_colors[i]=g2_ink(id_PS,grey_color,grey_color,grey_color);
  	}
  
  	ps_color_black=g2_ink(id_PS,0,0,0);
  	if(options.greyColors)
  		ps_color_red=g2_ink(id_PS,0,0,0);
  	else
  		ps_color_red=g2_ink(id_PS,1,0,0);
  
  	// draw sequence
  	g2_set_font_size(id,base_fontsize);
  	for(i=0;i<structure.size();i++)
  	{
  
  		if(options.mostLikelySequence)
  		{
  			double p=getMlBaseFreq(baseprobs[i]);
  
  			//base color
  			if(p==1)
  				g2_pen(id,ps_color_red);
  			else
  				g2_pen(id,ps_grey_colors[(int)floor(p*num_grey_colors-1)]);
  			  
  			base=getMlBase(baseprobs[i]);
  
  			xpos=X[i]-base.length()*base_fontsize/2.0;
  			ypos=Y[i]-4;
  			g2_string(id,xpos,ypos,(char*)base.c_str());     
  		}
  		else
  		{
  			drawBaseCircles(id_PS,baseprobs[i],X[i],Y[i]);
  		}
  
  		// connection to next base
  		if(i<structure.size()-1)
  		{
  			if((1-baseprobs[i].gap)*(1-baseprobs[i+1].gap)==1)
  				g2_pen(id,ps_color_red);
  			else
  				g2_pen(id,ps_grey_colors[(int)floor((1-baseprobs[i].gap)*(1-baseprobs[i+1].gap)*num_grey_colors-1)]);
  
  			g2_line(id,X[i],Y[i],X[i+1],Y[i+1]);
  		}
  	}
  
  	// draw pairings
  	// !!! pair_table indexing begins at 1 !!!
  	for(i=0;i<structure.size();i++)
  	{
  		if((unsigned short)pair_table[i+1]>i+1)
  		{	    	    
  			// pairs in both structures
  			if(pairprob[i]==1)
  				g2_pen(id,ps_color_red);
  			else
  				g2_pen(id,ps_grey_colors[(int)floor(pairprob[i]*num_grey_colors-1)]);	   	    
  
  			g2_line(id,X[i],Y[i],X[pair_table[i+1]-1],Y[pair_table[i+1]-1]);
  		}
  	}
  
  	g2_flush(id);
  	g2_close(id);
  
  	free(pair_table);
  	DELETE(X);
  	DELETE(Y);
  }
  #endif
  
  void RNAProfileAlignment::printSeqAli() const
  {
  	Uint i,l;
  	deque<string> seqs;
  	string seq,info;
  
  	// get alignment rows
  	for(i=0;i<m_numStructures;i++)
  	{
  	  
  		getSeqAli(seq,i,0,getMaxLength(0));
  		seqs.push_back(seq);
  	}
  
  	l=seq.length();
  
  	// sequence 
  	//	if(hasSequence)
  	//	{
  		// calculate info line
  		for(i=0;i<l;i++)
  		{
  				bool equal=true;
  			
  				for(Uint r=1;r<m_numStructures;r++)
  			{
  				if((seqs[r])[i]!=(seqs[r-1])[i])
  				{
  				equal=false;
  				break;
  				}
  			}
  
  			info += equal ? "*" : " ";		  
  		}  
  
  		// print it
  		// sequences
  		for(i=0;i<l;i+=55)
  		{
  			for(Uint r=0;r<m_numStructures;r++)
  				cout << setw(20) << setfill(' ') << left << m_strNames[r].substr(0,20) << setw(5) << " " << RNAFuncs::UpperCase(seqs[r].substr(i,55)) << endl;
  
  			cout << setw(25) << " " << info.substr(i,55) << endl;
  			cout << endl;
  		}
  		//	}
  }
  
  deque<pair<string,string> > RNAProfileAlignment::getSeqAli() const
  {
  	Uint i;
  	deque<pair<string,string> > seqAli;
  	pair<string,string> p;
  	string seq;
  
  	 // get alignment rows^M
           for(i=0;i<m_numStructures;i++)
           {
             getSeqAli(seq,i,0,getMaxLength(0));
  	   p.first = m_strNames[i];
  	   p.second = seq;
             seqAli.push_back(p);
           }
  
  	 return seqAli; 
  }
  
  void RNAProfileAlignment::printStrAli() const
  {
          Uint i,l;
  	deque<string> strs;
  	string str,info;
  
  	// get alignment rows
  	for(i=0;i<m_numStructures;i++)
  	{
  		getStructAli(str,i);
  
  		strs.push_back(str);
  	}
  
  	l=str.length();
  
  
  	for(i=0;i<l;i++)
  	  {
  	    bool equal=true;
  	    
  	    for(Uint r=1;r<m_numStructures;r++)
  	      {
  		if((strs[r])[i]!=(strs[r-1])[i])
  		  {
  		    equal=false;
  		    break;
  		  }
                }
  	    
  	    info += equal ? "*" : " ";		  
  	  }  
  	
  	// structures
  	for(i=0;i<l;i+=55)
  	{
  		for(Uint r=0;r<m_numStructures;r++)
  			cout << setw(20) << setfill(' ') << left << m_strNames[r].substr(0,20) << setw(5) << " " << strs[r].substr(i,55) << endl;
  
  		cout << setw(25) << " " <<  info.substr(i,55) << endl;
  		cout << endl;
  	}
   
    
  }
  
  deque<string> RNAProfileAlignment::getStrAli() const
  {
  	Uint i;
  	deque<string> strs;
          string str;
  
          // get alignment rows^M
          for(i=0;i<m_numStructures;i++)
          {
            getStructAli(str,i);
            strs.push_back(str);
          }
  
  	return strs;
  }
  
  string RNAProfileAlignment::getConsSeq() const
  {
  	string seq;
  	deque<BaseProbs> baseprobs;
  	getSequenceAlignment(baseprobs);
  
  	 // build sequence
  	 deque<BaseProbs>::const_iterator it;
  	 for(it=baseprobs.begin();it!=baseprobs.end();it++)
  	 {
  	   seq += getMlBase(*it);
  	 }
  
  	 return seq;
  }
  
  string RNAProfileAlignment::getConsStr(double minPairProb) const
  {
  	string str;
  	deque<double> pairprob;
  	getStructureAlignment(minPairProb,str,pairprob);
  
  	return str;
  }
  
  deque<double> RNAProfileAlignment::getBaseProb() const
  {
  	deque<BaseProbs> baseprobs;
  	deque<double> bestBaseprob;
  
  	getSequenceAlignment(baseprobs);
          // build sequence
          deque<BaseProbs>::const_iterator it;
          for(it=baseprobs.begin();it!=baseprobs.end();it++)
          {
             bestBaseprob.push_back(getMlBaseFreq(*it));
          }
  
  	return bestBaseprob;
  }
  
  //deque<pair<pair<int,int>,double> > RNAProfileAlignment::getPairProb(double minPairProb)
  void RNAProfileAlignment::getPairProb(double &minPairProb, deque<pair<pair<int,int>,double> > &pairprobs)
  {
    //deque<pair<pair<int,int>,double> > structureProbs;
  	deque<double> pairprob;
  	string str;
  	
  	getStructureAlignment(minPairProb,str,pairprob);
  
  	// store the position of open brackets in an integer deque 
  	deque<int> openBrackets;
  	// iterator for structureProbs
  	pair<int,int> tmpPair;
  	pair<pair<int,int>,double> tmpPair2;
  	
  	// iterate through structure str	
  	for(int i=0;i<=str.length();i++)
  	{
  	  if(str[i]=='(') {
  	    openBrackets.push_back(i);
  	  } else if(str[i]==')') {
  	    // insert the int value of the corresponding opening bracket into structureProbs
  	    // insert position of the actual closing bracket
  	    // insert probability according to the actual bases
  	    tmpPair = pair<int,int> (openBrackets.back(),i);
  	    tmpPair2 = pair<pair<int,int>,double> (tmpPair,pairprob[i]);
  	    pairprobs.push_back(tmpPair2);
  	    // delete position of the last opening bracket
  	    openBrackets.pop_back();
  	  }
  	}
  }
  
  void RNAProfileAlignment::printFastaAli(bool noStructure) const
  {
  	string str,seq;
  
  	// get alignment rows
  	for(Uint i=0;i<m_numStructures;i++)
  	{
  		getStructAli(str,i);	  
  		getSeqAli(seq,i,0,getMaxLength(0));
  
  		cout << ">" << m_strNames[i] << endl;
  		cout << seq << endl;
  		if(!noStructure)
  		  cout << str << endl;
  	}
  }
  
  void RNAProfileAlignment::printConsensus(double minPairProb) const
  {
  	const int mountain_high=10;
  	int j;
  
  	string seq,structure;
  
  	deque<BaseProbs> baseprobs;
  	deque<double> pairprob;
  
  	getSequenceAlignment(baseprobs);
  	getStructureAlignment(minPairProb,structure,pairprob);
  
  	// build sequence
  	deque<BaseProbs>::const_iterator it;
  	for(it=baseprobs.begin();it!=baseprobs.end();it++)
  	{
  		seq += getMlBase(*it);
  	}
  
  //	assert(seq.size() == structure.size());
  	if(seq.size() != structure.size())
  		cerr <<  "Error in resolving consensus structure!" << endl;
  
  	for(Uint i=0;i<seq.length();i+=55)
  	{
  		// show structure frequency mountain
  		for(j=0;j<mountain_high;j++)
  		{
  			cout << setw(20) << setfill(' ') << " ";
  			cout << setw(3) << right << 100 - (100 / mountain_high) * j << "% ";
  
  			for(int k=0;k<55;k++)
  			{
  				if(i+k<seq.length())
  				{
  					if(getMlBaseFreq(baseprobs[i+k])>=1-(1/(double)mountain_high)*(double)j)
  						cout << "*";
  					else
  						cout << " ";
  				}
  			}
  			cout << endl;
  		}
  
  		cout << setw(25) << " " << RNAFuncs::UpperCase(seq.substr(i,55)) << endl;
  		cout << setw(25) << " " << structure.substr(i,55) << endl;
  
  		// show structure frequency mountain
  		for(j=0;j<mountain_high;j++)
  		{
  			cout << setw(20) << setfill(' ') << " ";
  			cout << setw(3) << right << (100 / mountain_high) * (j+1) << "% ";
  
  			for(int k=0;k<55;k++)
  			{		
  				if(i+k<seq.length())
  				{
  					if(pairprob[i+k]>=(1/(double)mountain_high)*(double)j)
  						cout << "*";
  					else
  						cout << " ";
  				}
  			}
  			cout << endl;
  		}
  
  		cout << endl;
  	}
  }
  
  void RNAProfileAlignment::addStrNames(const deque<string>& strNames)
  {
  	deque<string>::const_iterator it;
  	for(it=strNames.begin();it!=strNames.end();it++)
  		m_strNames.push_back(*it);
  }
  
  void RNAProfileAlignment::save(const string &filename)
  {
    ofstream s(filename.c_str());
  
    // save the pforest to stream in binary format
    
    // first of all save the size
    s.write(reinterpret_cast<const char*>(&m_size),sizeof(size_type));
    s.write(reinterpret_cast<const char*>(&m_numStructures),sizeof(Uint));
    
    // save the arrays
    for(Uint i=0;i<m_size;i++)
      {
      for(int r=0;r<RNA_ALPHABET_SIZE;r++)
        s.write(reinterpret_cast<const char*>(&m_lb[i].p[r]),sizeof(double));
  
      s.write(reinterpret_cast<const char*>(m_lb[i].columnStr.c_str()),sizeof(char)*m_numStructures);
      }
    
    
    //  for(Uint i=0;i<m_size;i++)
    //  {
    //    s.write(reinterpret_cast<char*>(&m_lb[i]),sizeof(label_type)*m_size);
    //  }
  
    //  s.write(reinterpret_cast<char*>(m_lb),sizeof(label_type)*m_size);
    s.write(reinterpret_cast<const char*>(m_rb),sizeof(size_type)*m_size);
    s.write(reinterpret_cast<const char*>(m_noc),sizeof(size_type)*m_size);
    
    for(Uint r=0;r<m_numStructures;r++)
      {
        ostringstream ss;
  
        ss << setw(20) << setfill(' ') << left <<  m_strNames[r];
        s.write(reinterpret_cast<const char*>(ss.str().c_str()),sizeof(char)*20);
      }
    
  
    //  s.write(reinterpret_cast<char*>(m_sumUpCSF),sizeof(size_type)*m_size);
    //  s.write(reinterpret_cast<char*>(m_rmb),sizeof(size_type)*m_size);
    
    //  s.write(m_name;
    //Uint m_numStructures;
    //Uint m_numStructuresX;
    //Uint m_numStructuresY;
    //deque<string> m_strNames;
    //bool hasSequence;
  }
  