Huang et al. Algorithm_BMLCS(Dynamic programming, 2008)

Main features
• searching phase in O(rnδ) time complexity;
Abstract of the paper [Huang 2008]

The longest common subsequence and sequence alignment problems have been studied extensively and they can be regarded as the relationship measurement between sequences. However, most of them treat sequences evenly or consider only two sequences. Recently, with the rise of whole-genome duplication research, the doubly conserved synteny relationship among three sequences should be considered. It is a brand new model to ﬁnd a merging way for understanding the interleaving relationship of sequences. Here, we deﬁne the merged LCS problem for measuring the interleaving relationship among three sequences. An O (n3 ) algorithm is ﬁrst proposed for solving the problem, where n is the sequence length. We further discuss the variant version of this problem with the block information. For the blocked merged LCS problem, we propose an algorithm with time complexity O (n2 m), where m is the number of blocks. An improved O (n2 + nm2 ) algorithm is further proposed for the same blocked problem.

C++ source code
```int BMLCS_DP(vector< int> T, vector< int> A, vector< int> B, int Tlen, int Alen, int Blen, vector< int> EOBA, vector< int> EOBB) {

vector< vector< int>> vat(Alen + 1, vector< int>(Tlen + 1));
vector< vector< int>> vbt(Blen + 1, vector< int>(Tlen + 1));
vector< vector< int>> vabt1(Alen + 1, vector< int>(Blen + 1));
vector< vector< int>> vabt2(Alen + 1, vector< int>(Blen + 1));
//use LCS algorithm to calcluate similiar between sequence A and sequence T
vat[0][0] = 0;
for (int i = 1; i <= Alen; ++i) {
vat[i][0] = 0;
}
for (int j = 1; j <= Tlen; ++j) {
vat[0][j] = 0;
}
for (int i = 1; i <= Alen; ++i) {
for (int j = 1; j <= Tlen; ++j) {
if (A[i] != T[j]) {
vat[i][j] = max(vat[i - 1][j], vat[i][j - 1]);
}
else {
vat[i][j] = vat[i - 1][j - 1] + 1;
}
}
}

//use LCS algorithm to calcluate similiar between sequence B and sequence T
vbt[0][0] = 0;
for (int i = 1; i <= Blen; ++i) {
vbt[i][0] = 0;
}
for (int j = 1; j <= Tlen; ++j) {
vbt[0][j] = 0;
}
for (int i = 1; i <= Blen; ++i) {
for (int j = 1; j <= Tlen; ++j) {
if (B[i] != T[j]) {
vbt[i][j] = max(vbt[i - 1][j], vbt[i][j - 1]);
}
else {
vbt[i][j] = vbt[i - 1][j - 1] + 1;
}
}
}
for (int i = 0; i <= Alen; ++i) {
for (int j = 0; j <= Blen; ++j) {
vabt1[i][j] = 0;
}
}
for (int i = 0; i <= Alen; ++i) {
for (int j = 0; j <= Blen; ++j) {
vabt2[i][j] = 0;
}
}
//use BMLCS-DP algorithm

int max = 0;
int cali = 0, calj = 0;
int tmpb = 0;
int layer = 1;
while (1) {

vabt1[0][0] = 0;
vabt2[0][0] = 0;
if (layer % 2 == 1) {
cali = 0, calj = 0;
for (int i = 1; i <= Alen; ++i) {
vabt1[i][0] = vat[i][layer];
if (i < EOBA[cali]) {   // ai!=EOB amd bj = EOB
for (int j = 0; j < EOBB.size(); ++j) {
max = 0;
tmpb = EOBB[j];
if (A[i] == T[layer]) {
if ((vabt2[i - 1][tmpb] + 1) > max)
max = vabt2[i - 1][tmpb] + 1;
}
else {
if (vabt1[i - 1][tmpb] > max)
max = vabt1[i - 1][tmpb];
if (vabt2[i][tmpb] > max)
max = vabt2[i][tmpb];
}
vabt1[i][tmpb] = max;
}
}
else {// i = EOBA[cali]
calj = 0;
for (int j = 1; j <= Blen; ++j) {
vabt1[0][j] = vbt[j][layer];
max = 0;
if (j < EOBB[calj]) {// ai=EOB amd bj!= EOB
if (B[j] == T[layer]) {
if ((vabt2[i][j - 1] + 1) > max)
max = vabt2[i][j - 1] + 1;
}
else {
if (vabt1[i][j - 1] > max)
max = vabt1[i][j - 1];
if (vabt2[i][j] > max)
max = max = vabt2[i][j];
}
vabt1[i][j] = max;
}
else {  // ai=EOB amd bj= EOB
if (A[i] == T[layer]) {
if ((vabt2[i - 1][j] + 1) > max)
max = vabt2[i - 1][j] + 1;
}
else {
if (vabt1[i - 1][j] > max)
max = vabt1[i - 1][j];
}
if (B[j] == T[layer]) {
if ((vabt2[i][j - 1] + 1) > max)
max = vabt2[i][j - 1] + 1;
}
else {
if (vabt1[i][j - 1] > max)
max = vabt1[i][j - 1];
}
if (A[i] != T[layer] || B[j] != T[layer]) {
if (vabt2[i][j] > max)
max = vabt2[i][j];
}
vabt1[i][j] = max;
++calj;
}
}
++cali;
}
}
if (layer == Tlen)
return vabt1[Alen][Blen];
++layer;
}
else {
cali = 0, calj = 0;
for (int i = 1; i <= Alen; ++i) {
vabt2[i][0] = vat[i][layer];
if (i < EOBA[cali]) {   // ai!=EOB amd bj = EOB
for (int j = 0; j < EOBB.size(); ++j) {
max = 0;
tmpb = EOBB[j];
if (A[i] == T[layer]) {
if ((vabt1[i - 1][tmpb] + 1) > max)
max = vabt1[i - 1][tmpb] + 1;
}
else {
if (vabt2[i - 1][tmpb] > max)
max = vabt2[i - 1][tmpb];
if (vabt1[i][tmpb] > max)
max = vabt1[i][tmpb];
}
vabt2[i][tmpb] = max;
}
}
else {// i = EOBA[cali]
calj = 0;
for (int j = 1; j <= Blen; ++j) {
vabt2[0][j] = vbt[j][layer];
max = 0;
if (j < EOBB[calj]) {// ai=EOB amd bj!= EOB
if (B[j] == T[layer]) {
if ((vabt1[i][j - 1] + 1) > max)
max = vabt1[i][j - 1] + 1;
}
else {
if (vabt2[i][j - 1] > max)
max = vabt2[i][j - 1];
if (vabt1[i][j] > max)
max = max = vabt1[i][j];
}
vabt2[i][j] = max;
}
else {  // ai=EOB amd bj= EOB
if (A[i] == T[layer]) {
if ((vabt1[i - 1][j] + 1) > max)
max = vabt1[i - 1][j] + 1;
}
else {
if (vabt2[i - 1][j] > max)
max = vabt2[i - 1][j];
}
if (B[j] == T[layer]) {
if ((vabt1[i][j - 1] + 1) > max)
max = vabt1[i][j - 1] + 1;
}
else {
if (vabt2[i][j - 1] > max)
max = vabt2[i][j - 1];
}
if (A[i] != T[layer] || B[j] != T[layer]) {
if (vabt1[i][j] > max)
max = vabt1[i][j];
}
vabt2[i][j] = max;
++calj;
}
}
++cali;
}
}
if (layer == Tlen)
return vabt2[Alen][Blen];
++layer;
}
}
}

```
The files
All_BMLCS.cpp

Reference
Huang, K.-S., Yang, C.-B., Tseng, K.-T., Ann, H.-Y., Peng, Y.-H., 2008. Eﬃcient algorithms for ﬁnding interleaving relationship between sequences. Information Processing Letters 105, 188–193.