alignment.t.cpp

ViennaRNA-1.6.1

/*
  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 <[e-mail]>
*/

#ifndef _ALIGNMENT_T_CPP_
#define _ALIGNMENT_T_CPP_

#include <algorithm>
#include <cassert>

#include "alignment.h"
#include "debug.h"

#include "misc.t.cpp"
#include "ppforest.t.cpp"

/* ****************************************** */
/*    Constructor and Destruktor functions    */
/* ****************************************** */

template<class R,class L,class AL>
Alignment<R,L,AL>::Alignment(const PPForest<L> *ppfx, const PPForest<L> *ppfy, const Algebra<R,L> &alg, bool local, bool noSpeedup)
: m_suboptimalsPercent(100)
{
  if(local)
    calculateLocal(ppfx,ppfy,alg,noSpeedup);
  else
    calculateGlobal(ppfx,ppfy,alg,noSpeedup);
}

template<class R,class L,class AL>
Alignment<R,L,AL>::Alignment(const PPForest<L> *ppfx, const PPForest<L> *ppfy, const RNA_Algebra<R,L> &rnaAlg)
: m_suboptimalsPercent(100)
{
	assert(ppfx != NULL);
	assert(ppfy != NULL);

	Ulong m,n,h,cols;
	long i,k;
	Uint j,l,r;

	R score,h_score;

	// alloc space for the score matrix, backtrack structure  and , if wanted, for the calculation-order-matrix
	m_mtrxSize=ppfx->getNumCSFs()*ppfy->getNumCSFs();

	m_mtrx=new R[m_mtrxSize];
	m_rowStart=new Ulong[ppfx->getNumCSFs()];
	m_ppfx = new PPForest<L>(*ppfx);				// copy the ppforests
	m_ppfy = new PPForest<L>(*ppfy);
	m_rnaAlg=&rnaAlg;
	m_alg=(const Algebra<R,L>*)&rnaAlg;
	m_localOptimum=rnaAlg.worst_score();

	// initialize variables
	m=ppfx->size();
	n=ppfy->size();
	cols=ppfy->getNumCSFs();

	m_rowStart[0]=0;
	for(h=1;h<ppfx->getNumCSFs();h++)
		m_rowStart[h]=m_rowStart[h-1]+cols;

	// align forests fx and fy

	// the easiest case .. 
	setMtrxVal(0,0,rnaAlg.empty());

	// align fx to the empty forest (fill first row of array)
	for(i=m-1;i>=0;i--)  // for all nodes in fx
	{
		for(j=1;j<=ppfx->getMaxLength(i);j++)  // for all non empty csfs induced by i
		{
			score = rnaAlg.del(ppfx->label(i),
				               getMtrxVal(ppfx->down(i),0),
				               getMtrxVal(ppfx->over(i,j),0));	  

			setMtrxVal(ppfx->indexpos(i,j),0,score);
		}
	}

	// align fy to the empty forest (fill first column of array)
	for(k=n-1;k>=0;k--)  // for all nodes in fx
	{
		for(l=1;l<=ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
		{
			score = rnaAlg.insert(getMtrxVal(0,ppfy->down(k)),
					              ppfy->label(k),
					              getMtrxVal(0,ppfy->over(k,l)));

			setMtrxVal(0,ppfy->indexpos(k,l),score);
		}
	}

	// align the rest
	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<=ppfx->getMaxLength(i);j++)    // for all non empty csfs induced by i
				for(l=1;l<=ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
				{
					// basepair replacement
					if(ppfx->down(i) && ppfy->down(k))
					{
						// must be two P nodes !!
						h_score = rnaAlg.replacepair(ppfx->label(i+1),
							                         ppfy->label(k+1),
							                         getMtrxVal(ppfx->mdown(i),ppfy->mdown(k)),
							                         ppfx->label(ppfx->getRightmostBrotherIndex(i+1)),
							                         ppfy->label(ppfy->getRightmostBrotherIndex(k+1)),
							                         getMtrxVal(ppfx->over(i,j),ppfy->over(k,l)));
					}
					else
					{
						h_score = rnaAlg.replace(ppfx->label(i),
							                     getMtrxVal(ppfx->down(i),ppfy->down(k)),
							                     ppfy->label(k),
							                     getMtrxVal(ppfx->over(i,j),ppfy->over(k,l)));
					}

					score=h_score;

					// delete
					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 = rnaAlg.del(ppfx->label(i),
							                 getMtrxVal(ppfx->down(i),ppfy->indexpos(k,r)),
							                 getMtrxVal(ppfx->over(i,j),ppfy->indexpos(h,l-r)));

						score=rnaAlg.choice(score,h_score);
						h=ppfy->rb(h);
					}

					// insert
					h=i;
					for(r=0;r<=j;r++) // for all splits of fx
					{
						h_score = rnaAlg.insert(getMtrxVal(ppfx->indexpos(i,r),ppfy->down(k)),
							                    ppfy->label(k),
							                    getMtrxVal(ppfx->indexpos(h,j-r),ppfy->over(k,l)));

						score=rnaAlg.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>
Alignment<R,L,AL>::~Alignment()
{
    delete[] m_mtrx;
    delete[] m_rowStart;
    delete m_ppfx;
    delete m_ppfy;
}



/* ****************************************** */
/*            Private functions               */
/* ****************************************** */

template<class R,class L,class AL>
void Alignment<R,L,AL>::calculateLocal(const PPForest<L> *ppfx, const PPForest<L> *ppfy, const Algebra<R,L> &alg, bool noSpeedup)
{
	assert(ppfx != NULL);
	assert(ppfy != NULL);

	Ulong m,n,h,cols;
	long i,k;
	Uint j,l,r;

	R score,h_score;

	// alloc space for the score matrix, backtrack structure  and , if wanted, for the calculation-order-matrix
	m_mtrxSize=ppfx->getNumCSFs()*ppfy->getNumCSFs();

	m_mtrx=new R[m_mtrxSize];
	m_rowStart=new Ulong[ppfx->getNumCSFs()];
	m_ppfx = new PPForest<L>(*ppfx);				// copy the ppforests
	m_ppfy = new PPForest<L>(*ppfy);
	m_alg=&alg;
	m_rnaAlg=NULL;
	m_localOptimum=alg.worst_score();


	// initialize variables
	m=ppfx->size();
	n=ppfy->size();
	cols=ppfy->getNumCSFs();

	m_rowStart[0]=0;
	for(h=1;h<ppfx->getNumCSFs();h++)
		m_rowStart[h]=m_rowStart[h-1]+cols;

	// align forests fx and fy

	// the easiest case .. 
	setMtrxVal(0,0,alg.empty());

	// align fx to the empty forest (fill first row of array)
	for(i=m-1;i>=0;i--)  // for all nodes in fx
	{
		for(j=1;j<=ppfx->getMaxLength(i);j++)  // for all non empty csfs induced by i
		{
			score = alg.del(ppfx->label(i),
				            getMtrxVal(ppfx->down(i),0),
				            getMtrxVal(ppfx->over(i,j),0));	  

			setMtrxVal(ppfx->indexpos(i,j),0,score);
		}
	}

	// align fy to the empty forest (fill first column of array)
	for(k=n-1;k>=0;k--)  // for all nodes in fx
	{
		for(l=1;l<=ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
		{
			score = alg.insert(getMtrxVal(0,ppfy->down(k)),
					            ppfy->label(k),
					            getMtrxVal(0,ppfy->over(k,l)));

			setMtrxVal(0,ppfy->indexpos(k,l),score);
		}
	}

	// align the rest
	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<=ppfx->getMaxLength(i);j++)    // for all non empty csfs induced by i
				for(l=1;l<=ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
				{
					// 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);
				}

	/*
					// delete
					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
					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>
void Alignment<R,L,AL>::calculateGlobal(const PPForest<L> *ppfx, const PPForest<L> *ppfy, const Algebra<R,L> &alg, bool noSpeedup)
{
	assert(ppfx != NULL);
	assert(ppfy != NULL);

	Ulong m,n,h,cols;
	long i,k;
	Uint j,l,r;

	R score,h_score;

	// alloc space for the score matrix, backtrack structure  and , if wanted, for the calculation-order-matrix
	m_mtrxSize=ppfx->getNumCSFs()*ppfy->getNumCSFs();

	m_mtrx=new R[m_mtrxSize];
	m_rowStart=new Ulong[ppfx->getNumCSFs()];
	m_ppfx = new PPForest<L>(*ppfx);				// copy the ppforests
	m_ppfy = new PPForest<L>(*ppfy);
	m_alg=&alg;
	m_rnaAlg=NULL;
	m_localOptimum=alg.worst_score();


	// initialize variables
	m=ppfx->size();
	n=ppfy->size();
	cols=ppfy->getNumCSFs();

	m_rowStart[0]=0;
	for(h=1;h<ppfx->getNumCSFs();h++)
		m_rowStart[h]=m_rowStart[h-1]+cols;

	// align forests fx and fy

	// the easiest case .. 
	setMtrxVal(0,0,alg.empty());

	// align fx to the empty forest (fill first row of array)
	for(i=m-1;i>=0;i--)  // for all nodes in fx
	{
		for(j=1;j<=ppfx->getMaxLength(i);j++)  // for all non empty csfs induced by i
		{
			score = alg.del(ppfx->label(i),
				            getMtrxVal(ppfx->down(i),0),
				            getMtrxVal(ppfx->over(i,j),0));	  

			setMtrxVal(ppfx->indexpos(i,j),0,score);
		}
	}

	// align fy to the empty forest (fill first column of array)
	for(k=n-1;k>=0;k--)  // for all nodes in fx
	{
		for(l=1;l<=ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
		{
			score = alg.insert(getMtrxVal(0,ppfy->down(k)),
					            ppfy->label(k),
					            getMtrxVal(0,ppfy->over(k,l)));

			setMtrxVal(0,ppfy->indexpos(k,l),score);
		}
	}

	// align the rest
	for(i=m-1;i>=0;i--)  // for all nodes in fx  
		for(k=n-1;k>=0;k--)  // for all nodes in fx
		        {
			        j=ppfx->getMaxLength(i);
				for(l=1;l<=ppfy->getMaxLength(k);l++)  // for all non empty csfs induced by k
				{
					// 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);
						  }
					      }