new mafft v 6.857 with extensions

[jabaws.git] / binaries / src / mafft / extensions / mxscarna_src / AlifoldMEA.cpp

diff --git a/binaries/src/mafft/extensions/mxscarna_src/AlifoldMEA.cpp b/binaries/src/mafft/extensions/mxscarna_src/AlifoldMEA.cpp
new file mode 100644 (file)
index 0000000..340353f
--- /dev/null
+++ b/binaries/src/mafft/extensions/mxscarna_src/AlifoldMEA.cpp
@@ -0,0 +1,211 @@
+#include "AlifoldMEA.h"
+
+namespace MXSCARNA{
+
+const int AlifoldMEA::TURN = 3;
+
+void
+AlifoldMEA::
+Run()
+{
+    makeProfileBPPMatrix(alignment);
+    Initialization();
+    DP();
+    TraceBack();
+}
+
+void
+AlifoldMEA::
+makeProfileBPPMatrix(const MultiSequence *Sequences)
+{
+    int length = Sequences->GetSequence(0)->GetLength();
+
+    Trimat<float> *consBppMat = new Trimat<float>(length + 1);
+    fill(consBppMat->begin(), consBppMat->end(), 0);
+
+    for(int i = 1; i <= length; i++) 
+       for (int j = i; j <= length; j++) 
+           bppMat.ref(i, j) = 0;
+
+
+    int number = Sequences->GetNumSequences();
+    for(int seqNum = 0; seqNum < number; seqNum++) {
+       SafeVector<int> *tmpMap = Sequences->GetSequence(seqNum)->GetMappingNumber();
+       int label = Sequences->GetSequence(seqNum)->GetLabel();
+       BPPMatrix *tmpBppMatrix = BPPMatrices[label];
+       
+       for(int i = 1; i <= length ; i++) {
+           int originI = tmpMap->at(i);
+           for(int j = i; j <= length; j++) {
+               int originJ = tmpMap->at(j);
+               if(originI != 0 && originJ != 0) {
+                   float tmpProb = tmpBppMatrix->GetProb(originI, originJ);
+                   bppMat.ref(i, j) += tmpProb;
+               }
+           }
+       }
+    }
+
+       /* compute the mean of base pairing probability  */
+    for(int i = 1; i <= length; i++) {
+       for(int j = i; j <= length; j++) {
+           bppMat.ref(i,j) = bppMat.ref(i,j)/(float)number;
+       }
+    }
+
+    for (int i = 1; i <= length; i++) {
+       float sum = 0;
+       for (int j = i; j <= length; j++) {
+           sum += bppMat.ref(i,j);
+       }
+       Qi[i] = 1 - sum;
+    }
+
+    for (int i = 1; i <= length; i++) {
+       float sum = 0;
+       for (int j = i; j >= 1; j--) {
+           sum += bppMat.ref(j, i);
+       }
+       Qj[i] = 1 - sum;
+    }
+}
+
+void
+AlifoldMEA::
+Initialization()
+{
+    int length = alignment->GetSequence(0)->GetLength();
+
+    for (int i = 1; i <= length; i++) {
+       for (int j = i; j <= length; j++) {
+           M.ref(i,j) = 0;
+           traceI.ref(i,j) = 0;
+           traceJ.ref(i,j) = 0;
+       }
+    }
+
+    for (int i = 1; i <= length; i++) {
+       M.ref(i,i)   = Qi[i]; 
+       traceI.ref(i,i) = 0;
+       traceJ.ref(i,i) = 0;
+    }
+
+    for (int i = 1; i <= length - 1; i++) {
+       M.ref(i, i+1) =  Qi[i+1];
+       traceI.ref(i,i + 1) = 0;
+       traceJ.ref(i,i + 1) = 0;
+    }
+
+    for (int i = 0; i <= length; i++) {
+       ssCons[i] = '.';
+    }
+}
+
+void
+AlifoldMEA::
+DP()
+{
+    float g    = BasePairConst; // see scarna.hpp
+    int length = alignment->GetSequence(0)->GetLength();
+    
+    for (int i = length - 1; i >= 1; i--) {
+       for (int j = i + TURN + 1; j <= length; j++) {
+           float qi       = Qi[i];
+           float qj       = Qj[j];
+           float p        = bppMat.ref(i,j);
+
+           
+           float maxScore = qi + M.ref(i+1, j);
+           int tmpI = i+1;
+           int tmpJ = j;
+           
+           float tmpScore = qj + M.ref(i, j-1);
+           if (tmpScore > maxScore) {
+               maxScore = tmpScore;
+               tmpI     = i;
+               tmpJ     = j - 1;
+           }
+           
+           tmpScore = g*2*p + M.ref(i+1, j-1);
+           if (tmpScore > maxScore) {
+               maxScore = tmpScore;
+               tmpI     = i + 1;
+               tmpJ     = j - 1;
+           }
+           
+           for (int k = i + 1; k < j - 1; k++) {
+               tmpScore = M.ref(i,k) + M.ref(k+1,j);
+               if (tmpScore > maxScore) {
+                   maxScore = tmpScore;
+                   tmpI = i;
+                   tmpJ = j;
+               }
+           }
+           M.ref(i,j)       = maxScore;
+           traceI.ref(i, j) = tmpI;
+           traceJ.ref(i, j) = tmpJ;
+       }
+    }
+}
+
+void
+AlifoldMEA::
+TraceBack()
+{
+
+    int length = alignment->GetSequence(0)->GetLength();
+    SafeVector<int> stackI((length + 1)*(length+1));
+    SafeVector<int> stackJ((length + 1)*(length+1));
+    int pt = 0;
+
+    stackI[pt] = traceI.ref(1, length);
+    stackJ[pt] = traceJ.ref(1, length);
+    ++pt;
+    
+    while(pt != 0) {
+       --pt;
+       int tmpI = stackI[pt];
+       int tmpJ = stackJ[pt];
+       int nextI = traceI.ref(tmpI, tmpJ);
+       int nextJ = traceJ.ref(tmpI, tmpJ);
+
+       if (tmpI < tmpJ) {
+           if (tmpI + 1  == nextI && tmpJ == nextJ) {
+               stackI[pt] = nextI;
+               stackJ[pt] = nextJ;
+               ++pt;
+           }
+           else if (tmpI == nextI && tmpJ - 1 == nextJ) {
+               stackI[pt] = nextI;
+               stackJ[pt] = nextJ;
+               ++pt;
+           }
+           else if (tmpI + 1 == nextI && tmpJ - 1== nextJ) {
+               stackI[pt] = nextI;
+               stackJ[pt] = nextJ;
+               ++pt;
+               ssCons[tmpI] = '(';
+               ssCons[tmpJ] = ')';
+           }
+           else if (tmpI == nextI && tmpJ == nextJ) {
+               float maxScore = IMPOSSIBLE;
+               int maxK = 0;
+
+               for (int k = tmpI + 1; k < tmpJ - 1; k++) {
+                   float tmpScore = M.ref(tmpI,k) + M.ref(k+1,tmpJ);
+                   if (tmpScore > maxScore) {
+                       maxScore = tmpScore;
+                       maxK = k;
+                   }
+               }
+               stackI[pt] = traceI.ref(tmpI, maxK);
+               stackJ[pt] = traceJ.ref(tmpI, maxK);
+               ++pt;
+               stackI[pt] = traceI.ref(maxK+1, tmpJ);
+               stackJ[pt] = traceJ.ref(maxK+1, tmpJ);
+               ++pt;
+           }
+       }
+    }
+}
+}