动态规划之最长公共子序列问题 C++实现

原理

在之前的文章当中，作者论述了什么是动态规划，这次，我们来看看，如何用动态规划解决最长公共子序列问题。

这个问题经常运用在判断两种生物的相似度—-DNA比对上。对比俩串的方式有很多种，例如如果一个串是另一个的字串，那么可以说两个串是相似的：如果将一个串转换为另一个串的操作很少，那么也可以说这两个串是相似的。另一种衡量俩串 S1,S2 的相似度方式为：寻找第3个串 S3 ，它的所有元素也都出现在 S1 和 S2 中，且在三个串中出现的顺序都相同，但在 S1,S2 中不要求连续出现。我们将最后一种相似的概念描命名最长公共子序列问题。

其形式化定义如下：给定一个序列 X=⟨x1,x2…,xm⟩ ，另一个序列 Y=⟨y1,y2,…,yk⟩ 满足如下条件时成为 X 的子序列（sunsequence），即存在一严格递增的 X 的下标序列 ⟨i1,i2,…,ik⟩ ，对所有 j=1,2,…,k 满足 ij=zj 。例如， Z=⟨B,C,D,B⟩ 是 X=⟨A,B,C,B,D,A,B⟩ 的子序列，对应的下标为 ⟨2,3,5,7⟩ 。
给定一个序列 X=⟨x1,x2…,xm⟩ ，对 i=1,2,…,k ，定义 X 的第 i 前缀为 Xi=⟨x1,x2…,xi⟩ 。

刻画最长公共子序列的特征

LCS 的最优子结构：令 X=⟨x1,x2…,xm⟩ 和 Y=⟨y1,y2,…,yn⟩ 为两个序列， Z=⟨z1,z2,…,zk⟩ 为 X 和 Y 的任意 LCS 。
1.如果 xm=yn ，则 zk=xm=yn 且 Zk−1 是 Xm−1 和 Yn−1 的一个 LCS 。
2.如果 xm≠yn ，则 zk≠xm 且 Z 是 Xm−1 和 Y 的一个 LCS 。
3.如果 xm≠yn ，则 zk≠yn 且 Z 是 X 和 Yn−1 的一个 LCS 。

一个递归解

我们定义 c[i,j] 表示 Xi 和 Yj 的 LCS 的长度。则根据 LCS 问题的最优子结构性质，可得如下公式：

c[i,j]=⎧⎩⎨⎪⎪0c[i−1,j−1]+1max(c[i,j−1],c[i−1,j])if i=0 or j=0if i,j>0 and xi=yjif i,j>0 and xi≠ yj

源代码

#include <iostream>
#include <utility>
#include <vector>
#include <string>

using namespace  std;

//ACCGTCGAGTGCGCGGAAGCCGGCCGAA & CTCGTTCGGAATGCCGTTGCTCTGTAAA
string temp_strX = { "#ACCGTCGAGTGCGCGGAAGCCGGCCGAA" }, temp_strY = { "#CTCGTTCGGAATGCCGTTGCTCTGTAAA" }; 

//Memoized of Lcs
pair<vector<vector<int>>,vector<vector<int>>> Lcs_Length(const string &temp_strX, const string &strY) {
    auto temp_m = temp_strX.size() - 1, temp_n = temp_strY.size() - 1;
    vector<vector<int>> temp_VecB, temp_VecC;

    temp_VecB.resize(temp_m + 1);
    temp_VecC.resize(temp_m + 1);

    for(auto &i : temp_VecB) {
        i.resize(temp_n + 1);
    }
    for(auto &i : temp_VecC) {
        i.resize(temp_n + 1);
    }

    for(auto i = 1; i <= temp_m; ++i) {
        temp_VecC[i][0] = 0;
    }
    for(auto j = 0; j <= temp_n; ++j) {
        temp_VecC[0][j] = 0;
    }

    for(auto i = 1; i <= temp_m; ++i) {
        for(auto j = 1; j <= temp_n; ++j) {
            if(temp_strX[i] == temp_strY[j]) {
                temp_VecC[i][j] = temp_VecC[i - 1][j - 1] + 1;
                temp_VecB[i][j] = -1;
            }
            else if(temp_VecC[i - 1][j] >= temp_VecC[i][j - 1]) {
                temp_VecC[i][j] = temp_VecC[i - 1][j];
                temp_VecB[i][j] = -2;
            }
            else {
                temp_VecC[i][j] = temp_VecC[i][j - 1];
                temp_VecB[i][j] = -3;
            }
        }
    }

    return make_pair(temp_VecC, temp_VecB);
}

//Print
void Print_Lcs(const vector<vector<int>> & temp_VecB, const string &temp_strX, const size_t &i, const size_t &j) {
    if(i == 0 || j == 0) {
        return;
    }

    if(temp_VecB[i][j] == -1) {
        Print_Lcs(temp_VecB, temp_strX, i - 1, j - 1);
        cout << temp_strX[i];
    }
    else if(temp_VecB[i][j] == -2) {
        Print_Lcs(temp_VecB, temp_strX, i - 1, j);
    }
    else {
        Print_Lcs(temp_VecB, temp_strX, i, j - 1);
    }
}

int main() {
    auto temp_pair = Lcs_Length(temp_strX, temp_strY);
    Print_Lcs(temp_pair.second, temp_strX, temp_strX.size() - 1, temp_strY.size() - 1);

    return 0;
}