alignment.t.cpp

ViennaRNA-1.6.1

					  }

					// set value
					setMtrxVal(ppfx->indexpos(i,j),ppfy->indexpos(k,l),score);
				}

				l=ppfy->getMaxLength(k);
				for(j=1;j<=ppfx->getMaxLength(i);j++)    // for all non empty csfs induced by i
				{
					// replace
					score = alg.replace(ppfx->label(i),
							            getMtrxVal(ppfx->down(i),ppfy->down(k)),
							            ppfy->label(k),
							            getMtrxVal(ppfx->over(i,j),ppfy->over(k,l)));

					// delete				       
					if(ppfx->noc(i)==0 && !noSpeedup)  // no child
					  {
					    h_score = alg.del(ppfx->label(i),0,
							              getMtrxVal(ppfx->over(i,j),ppfy->indexpos(k,l)));

						score=alg.choice(score,h_score);
					  }
					else					  
					  {	
					    if(ppfx->rb(i)==0 && !noSpeedup) // no right brother
					      {
						h_score = alg.del(ppfx->label(i),
								  getMtrxVal(ppfx->down(i),ppfy->indexpos(k,l)),
								  0);
						
						score=alg.choice(score,h_score);
					      }
					    else
					      {
						h=k;               // h is the node where the suffix of the split begins		   					
						for(r=0;r<=l;r++)  // for all splits of fy
						  {
						    h_score = alg.del(ppfx->label(i),
							              getMtrxVal(ppfx->down(i),ppfy->indexpos(k,r)),
							              getMtrxVal(ppfx->over(i,j),ppfy->indexpos(h,l-r)));
						    
						    score=alg.choice(score,h_score);
						    h=ppfy->rb(h);
						  }
					      }
					  }

					// insert
					if(ppfy->noc(k)==0 && !noSpeedup) // no child
					  {
					    h_score = alg.insert(0,
							                 ppfy->label(k),
							                 getMtrxVal(ppfx->indexpos(i,j),ppfy->over(k,l)));

						score=alg.choice(score,h_score);
					  }
					else
					  {
					    if(ppfy->rb(k)==0 && !noSpeedup) // no right brother
					      {
						h_score = alg.insert(getMtrxVal(ppfx->indexpos(i,j),ppfy->down(k)),
								     ppfy->label(k),
								     0);
						
						score=alg.choice(score,h_score);
					      }					 
					    else
					      {
						h=i;
						for(r=0;r<=j;r++) // for all splits of fx
						  {
						    h_score = alg.insert(getMtrxVal(ppfx->indexpos(i,r),ppfy->down(k)),
							                 ppfy->label(k),
									 getMtrxVal(ppfx->indexpos(h,j-r),ppfy->over(k,l)));
						    
						    score=alg.choice(score,h_score);
						    h=ppfx->rb(h);
						  }
					      }
					  }

					// set value
					setMtrxVal(ppfx->indexpos(i,j),ppfy->indexpos(k,l),score);
				}
			}

	//      showArray(m_mtrx,ppfx->getNumCSFs(),ppfy->getNumCSFs());
	resetOptLocalAlignment(100);
}

template<class R,class L,class AL>
Uint Alignment<R,L,AL>::backtrack(PPForestAli<L,AL> &ppf,Uint i, Uint j, Uint k, Uint l, Uint &node)
{
	R score, h_score;
	Uint p_node,rb_node,noc,rbs,r,h;

	// empty alignment
	if(j==0 && l==0)
	{
		return 0;
	}

	score = getMtrxVal(m_ppfx->indexpos(i,j),m_ppfy->indexpos(k,l));
	p_node=node;
	node++;

	// could it be a replacement
	if(j>0 && l>0)
	{
		// check for basepair replacement only if Algebra is of type RNA_Algebra
		if(m_rnaAlg && m_ppfx->down(i) && m_ppfy->down(k))
		{
			h_score = m_rnaAlg->replacepair(m_ppfx->label(i+1),
				                            m_ppfy->label(k+1),
				                            getMtrxVal(m_ppfx->mdown(i),m_ppfy->mdown(k)),
				                            m_ppfx->label(m_ppfx->getRightmostBrotherIndex2(i+1)),
				                            m_ppfy->label(m_ppfy->getRightmostBrotherIndex2(k+1)),
				                            getMtrxVal(m_ppfx->over(i,j),m_ppfy->over(k,l)));      

			if(score == h_score)
			{
				// it is a basepair replacement
				ppf.makeRepLabel(p_node,m_ppfx->label(i),m_ppfy->label(k));   // P

				// set labels of leftmost child
				ppf.makeRepLabel(p_node+1,m_ppfx->label(i+1),m_ppfy->label(k+1));
				ppf.setRightBrotherIndex(p_node+1,p_node+1+1);
				ppf.setNumChildren(p_node+1,0);  // base node has no children
				node++;

				// down alignment
				assert(m_ppfx->noc(i)>=2);
				assert(m_ppfy->noc(k)>=2);
				noc=backtrack(ppf,i+1+1,m_ppfx->noc(i)-2,k+1+1,m_ppfy->noc(k)-2,node);  // !! mdown !!
				ppf.setNumChildren(p_node,noc+2);

				if(noc==0)
				{
					ppf.setRightBrotherIndex(p_node+1,p_node+1+1);
					ppf.setRightBrotherIndex(p_node+1+1,0);
				}
				else
				{
					ppf.setRightBrotherIndex(ppf.getRightmostBrotherIndex2(p_node+1+1),node);
				}


				// set labels of leftmost child
				ppf.makeRepLabel(node,m_ppfx->label(m_ppfx->getRightmostBrotherIndex2(i+1+1)),m_ppfy->label(m_ppfy->getRightmostBrotherIndex2(k+1+1)));
				ppf.setRightBrotherIndex(node,0);
				ppf.setNumChildren(node,0);  // base node has no children	  
				node++;
				rb_node=node; // !!	  

				// right alignment
				rbs=backtrack(ppf,m_ppfx->rb(i),j-1,m_ppfy->rb(k),l-1,node);
				if(rbs)
					ppf.setRightBrotherIndex(p_node,rb_node);
				else
					ppf.setRightBrotherIndex(p_node,0);

				return rbs+1;
			}
		}
		else
		{
			h_score = m_alg->replace(m_ppfx->label(i),
				                     getMtrxVal(m_ppfx->down(i),m_ppfy->down(k)),
				                     m_ppfy->label(k),
				                     getMtrxVal(m_ppfx->over(i,j),m_ppfy->over(k,l)));

			if(score == h_score)
			{
				// it is a replacement
				ppf.makeRepLabel(p_node,m_ppfx->label(i),m_ppfy->label(k));

				// down alignment
				noc=backtrack(ppf,i+1,m_ppfx->noc(i),k+1,m_ppfy->noc(k),node);
				ppf.setNumChildren(p_node,noc);
				rb_node=node;

				// right alignment
				rbs=backtrack(ppf,m_ppfx->rb(i),j-1,m_ppfy->rb(k),l-1,node);
				if(rbs)
					ppf.setRightBrotherIndex(p_node,rb_node);
				else
					ppf.setRightBrotherIndex(p_node,0);

				return rbs+1;
			}
		}
	}

	// could it be a deletion 
	if(j>0)
	{
		h=k;               // h is the node where the suffix of the split begins
		for(r=0;r<=l;r++)  // for all splits of fy
		{
			h_score = m_alg->del(m_ppfx->label(i),
				                 getMtrxVal(m_ppfx->down(i),m_ppfy->indexpos(k,r)),
				                 getMtrxVal(m_ppfx->over(i,j),m_ppfy->indexpos(h,l-r)));

			if(score == h_score)
			{
				// it is a deletion
				ppf.makeDelLabel(p_node,m_ppfx->label(i));

				// down alignment
				noc=backtrack(ppf,i+1,m_ppfx->noc(i),k,r,node);
				ppf.setNumChildren(p_node,noc);
				rb_node=node;

				// right alignment
				rbs=backtrack(ppf,m_ppfx->rb(i),j-1,h,l-r,node);
				if(rbs)
					ppf.setRightBrotherIndex(p_node,rb_node);
				else
					ppf.setRightBrotherIndex(p_node,0);

				return rbs+1;
			}

			//	  if(r<l)       // do not calculate rightbrother of h=0=no-rightbrother
			h=m_ppfy->rb(h);
		}
	}

	// could it be an insertion
	if(l>0)
	{
		h=i;
		for(r=0;r<=j;r++) // for all splits of fx
		{
			h_score = m_alg->insert(getMtrxVal(m_ppfx->indexpos(i,r),m_ppfy->down(k)),
				                    m_ppfy->label(k),
				                    getMtrxVal(m_ppfx->indexpos(h,j-r),m_ppfy->over(k,l)));

			if(score == h_score)
			{
				// it is an insertion
				ppf.makeInsLabel(p_node,m_ppfy->label(k));

				// down alignment
				noc=backtrack(ppf,i,r,k+1,m_ppfy->noc(k),node);
				ppf.setNumChildren(p_node,noc);
				rb_node=node;

				// right alignment
				rbs=backtrack(ppf,h,j-r,m_ppfy->rb(k),l-1,node);
				if(rbs)
					ppf.setRightBrotherIndex(p_node,rb_node);
				else
					ppf.setRightBrotherIndex(p_node,0);

				return rbs+1;      
			}

			//	  if(r<j)
			h=m_ppfx->rb(h);
		}
	}

	// you should never be here
	cerr << "Strange things happening in backtrack" << endl;
	exit(EXIT_FAILURE);
} 



/* ****************************************** */
/*             Public functions               */
/* ****************************************** */

template<class R,class L,class AL>
R Alignment<R,L,AL>::getGlobalOptimum()
{
  return getMtrxVal(m_ppfx->indexpos(0,m_ppfx->getMaxLength(0)),m_ppfy->indexpos(0,m_ppfy->getMaxLength(0)));
};

template<class R,class L,class AL>
double Alignment<R,L,AL>::getGlobalOptimumRelative()
{
  double opt;
  double max_x,max_y;

  opt=(double)getGlobalOptimum();

  if(m_rnaAlg)
    {
      max_x=(double)m_ppfx->maxScore(*m_rnaAlg);
      max_y=(double)m_ppfy->maxScore(*m_rnaAlg);
    }
  else
    {
      max_x=(double)m_ppfx->maxScore(*m_alg);
      max_y=(double)m_ppfy->maxScore(*m_alg);
    }

  assert(max_x+max_y>0);	
  opt=2*opt/(max_x+max_y);  

  return opt;
};

template<class R,class L,class AL>
R Alignment<R,L,AL>::getSILOptimum()
{
  R silOptimum;
  int k=0;
  Uint j,l=0;
  Ulong n;
  
  n=m_ppfy->size();

  if(m_rnaAlg)
    silOptimum=m_rnaAlg->worst_score();
  else
    silOptimum=m_alg->worst_score();

  j=m_ppfx->getMaxLength(0);

  // find the best match
  for(k=n-1;k>=0;k--)  // for all nodes in fx
    for(l=1;l<=m_ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
      {
	silOptimum=m_alg->choice(silOptimum,getMtrxVal(m_ppfx->indexpos(0,j),m_ppfy->indexpos(k,l)));
      }

  return silOptimum;
}

template<class R,class L,class AL>
void Alignment<R,L,AL>::getOptGlobalAlignment(PPForestAli<L,AL> &ppfali)
{
  Uint node=0;

  // allocate a forest of the maximal size that a forest alignment can have
  ppfali.initialize(m_ppfx->size()+m_ppfy->size());
  backtrack(ppfali,0,m_ppfx->getMaxLength(0),0,m_ppfy->getMaxLength(0),node);
  ppfali.setSize(node);
  ppfali.calcSumUpCSF();
  ppfali.calcRMB();
}

template<class R,class L,class AL>
void Alignment<R,L,AL>::resetOptLocalAlignment(int suboptimalsPercent)
{
	m_localAlis.clear();
	m_suboptimalsPercent=suboptimalsPercent/100.0;
	m_localSubOptimum=m_localOptimum;
};

// calculate the score of the next best local alignment that is not "included" 
// in a local alignment returned by getOptLocalAlignment before 
template<class R,class L,class AL>
bool Alignment<R,L,AL>::nextLocalSuboptimum()
{
	Ulong m,n;
	int i=0,k=0;
	Uint j=0,l=0;
	
  m_localSubOptimum=m_alg->worst_score();
  m=m_ppfx->size();
  n=m_ppfy->size();

  // find a matrix element that is optimal
  for(i=m-1;i>=0;i--)  // for all nodes in fx  
    for(k=n-1;k>=0;k--)  // for all nodes in fx
      for(j=1;j<=m_ppfx->getMaxLength(i);j++)    // for all non empty csfs induced by i
		for(l=1;l<=m_ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
		{
			bool disjoint=true;

			// check if i,j,k,l is included 
			typename list<CSFPair>::const_iterator it;
			for(it=m_localAlis.begin();it!=m_localAlis.end();it++)
			{
				if(!m_ppfx->isDisjoint(it->i,it->j,i,j) || !m_ppfy->isDisjoint(k,l,it->k,it->l))
				{
					disjoint=false;
					break;
				}
				
			}

			if(disjoint)
			{
				if(getMtrxVal(m_ppfx->indexpos(i,j),m_ppfy->indexpos(k,l))>=m_suboptimalsPercent*m_localOptimum)
					m_localSubOptimum=m_alg->choice(m_localSubOptimum,getMtrxVal(m_ppfx->indexpos(i,j),m_ppfy->indexpos(k,l)));
			}
		}

	return m_localSubOptimum!=m_alg->worst_score();	
}

template<class R,class L,class AL>
void Alignment<R,L,AL>::getOptLocalAlignment(PPForestAli<L,AL> &ppfali,Uint &xbasepos, Uint &ybasepos)
{
  Ulong m,n;
  int i=0,k=0;
  Uint j=0,l=0,node=0;

  m=m_ppfx->size();
  n=m_ppfy->size();
  // allocate a forest of the maximal size that a forest alignment can have
  ppfali.initialize(m_ppfx->size()+m_ppfy->size());

  // find a matrix element that is optimal
  for(i=m-1;i>=0;i--)  // for all nodes in fx  
    for(k=n-1;k>=0;k--)  // for all nodes in fx
      for(j=1;j<=m_ppfx->getMaxLength(i);j++)    // for all non empty csfs induced by i
		for(l=1;l<=m_ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
		{
			bool disjoint=true;

			// check if i,j,k,l is included 
			typename list<CSFPair>::const_iterator it;
			for(it=m_localAlis.begin();it!=m_localAlis.end();it++)
			{
				if(!m_ppfx->isDisjoint(it->i,it->j,i,j) || !m_ppfy->isDisjoint(k,l,it->k,it->l))
				{
					disjoint=false;
					break;
				}
				
			}

			if(disjoint && getMtrxVal(m_ppfx->indexpos(i,j),m_ppfy->indexpos(k,l)) == m_localSubOptimum)
				goto found;
		}

 found: 
  backtrack(ppfali,i,j,k,l,node);
  ppfali.setSize(node);
  ppfali.calcSumUpCSF();
  ppfali.calcRMB();  

  m_localAlis.push_back(CSFPair(i,j,k,l));
  xbasepos=m_ppfx->countLeaves(i);
  ybasepos=m_ppfy->countLeaves(k);
}

template<class R,class L,class AL>
void Alignment<R,L,AL>::getOptSILAlignment(PPForestAli<L,AL> &ppfali,Uint &ybasepos)
{
  R silOptimum;
  Ulong m,n;
  int k=0;
  Uint j=0,l=0,node=0;

  m=m_ppfx->size();
  n=m_ppfy->size();
  // allocate a forest of the maximal size that a forest alignment can have
  ppfali.initialize(m_ppfx->size()+m_ppfy->size());

  silOptimum=getSILOptimum();
  j=m_ppfx->getMaxLength(0);
  
  // find a matrix element that is optimal
  for(k=n-1;k>=0;k--)  // for all nodes in fx
    for(l=1;l<=m_ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
      {
	if(silOptimum == getMtrxVal(m_ppfx->indexpos(0,j),m_ppfy->indexpos(k,l)))
	  goto found;
      }

 found: 
  backtrack(ppfali,0,j,k,l,node);
  ppfali.setSize(node);
  ppfali.calcSumUpCSF();
  ppfali.calcRMB();  

  ybasepos=m_ppfy->countLeaves(k);
}


#endif