--- /dev/null
+/*
+ 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 _ALIGNMENT_T_CPP_
+#define _ALIGNMENT_T_CPP_
+
+#include <algorithm>
+#include <cassert>
+#include <new>
+#include <exception>
+#include <limits.h>
+
+#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
+ // check for an overflow
+ if (ppfx->getNumCSFs() > ULONG_MAX / ppfy->getNumCSFs()) {
+ cerr << "Error: Overflow in calculation matrix multiplication. Calculation terminated." << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ m_mtrxSize=ppfx->getNumCSFs()*ppfy->getNumCSFs();
+
+ // maximum array size is 2GB
+ if(m_mtrxSize > 2000000000 || ppfx->getNumCSFs() > 2000000000) {
+ cerr << "Error: Maximum array size of 2GB exceeded due to large input data. Calculation terminated." << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ 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);
+ }
+
+ 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
+ if (ppfx->getNumCSFs() > ULONG_MAX / ppfy->getNumCSFs()) {
+ cerr << "Error: Overflow in calculation matrix multiplication. Calculation terminated." << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ m_mtrxSize=ppfx->getNumCSFs()*ppfy->getNumCSFs();
+
+ if(m_mtrxSize > 2000000000 || ppfx->getNumCSFs() > 2000000000) {
+ cerr << "Error: Maximum array size of 2GB exceeded due to large input data. Calculation terminated." << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ 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
+ if (ppfx->getNumCSFs() > ULONG_MAX / ppfy->getNumCSFs()) {
+ cerr << "Error: Overflow in calculation matrix multiplication. Calculation terminated." << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ m_mtrxSize=ppfx->getNumCSFs()*ppfy->getNumCSFs();
+
+ if(m_mtrxSize > 2000000000 || ppfx->getNumCSFs() > 2000000000) {
+ cerr << "Error: Maximum array size of 2GB exceeded due to large input data. Calculation terminated." << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ 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);
+ }
+ }
+ }
+
+ // 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
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
git://source.jalview.org / jabaws.git / blobdiff