动态规划解最长公共子序列问题

1 最长公共子序列问题概述

1.1 问题定义

序列A={a1,a2,...,am}A=\{a_1,a_2,...,a_m\}A={a1​,a2​,...,am​}，序列B={b1,b2,...,bn}B=\{b_1,b_2,...,b_n\}B={b1​,b2​,...,bn​}。如果存在一个序列C={c1,c2,...,ck}C=\{c_1,c_2,...,c_k\}C={c1​,c2​,...,ck​}，其中$c_i\in A,B$ $(i=1,2,...,k)$，且$c_1,c_2,...,c_k$在$A,B$中出现的先后顺序要保持一致，则称序列$C$是序列$A$和序列$B$的公共子序列。

要求找出序列$A$和序列$B$的最长的公共子序列。

1.2 举个栗子

序列A={1,2,3}A=\{1,2,3\}A={1,2,3}，序列B={2,3,4}B=\{2,3,4\}B={2,3,4}。

序列$A$的子序列有：{1},{2},{3},{1,2},{1,3},{2,3},{1,2,3}\{1\},\{2\},\{3\},\{1,2\},\{1,3\},\{2,3\},\{1,2,3\}{1},{2},{3},{1,2},{1,3},{2,3},{1,2,3}

序列$B$的子序列有：{2},{3},{4},{2,3},{2,4},{3,4},{2,3,4}\{2\},\{3\},\{4\},\{2,3\},\{2,4\},\{3,4\},\{2,3,4\}{2},{3},{4},{2,3},{2,4},{3,4},{2,3,4}

序列$A$和$B$的公共子序列有：{2},{3},{2,3}\{2\},\{3\},\{2,3\}{2},{3},{2,3}

所以序列$A$和$B$的最长公共子序列是：{2,3}\{2,3\}{2,3}

2 动态规划求解最长公共子序列问题

序列A={a1,a2,...,am}A=\{a_1,a_2,...,a_m\}A={a1​,a2​,...,am​}，序列B={b1,b2,...,bn}B=\{b_1,b_2,...,b_n\}B={b1​,b2​,...,bn​}。

假设序列C={c1,c2,...,ck}C=\{c_1,c_2,...,c_k\}C={c1​,c2​,...,ck​}是两个序列的最长公共子序列，则有以下推论：

• 如果$a_m=b_n$，则$c_k=a_m=b_n$，且序列{c1,c2,...,ck−1}\{c_1,c_2,...,c_{k-1}\}{c1​,c2​,...,ck−1​}是序列{a1,a2,...,am−1}\{a_1,a_2,...,a_{m-1}\}{a1​,a2​,...,am−1​}和序列{b1,b2,...,bn−1}\{b_1,b_2,...,b_{n-1}\}{b1​,b2​,...,bn−1​}的最长公共子序列
• 如果$a_m\ne b_n$，且$c_k\ne a_m$，则$c_k=b_n$，序列{c1,c2,...,ck}\{c_1,c_2,...,c_{k}\}{c1​,c2​,...,ck​}是序列{a1,a2,...,am−1}\{a_1,a_2,...,a_{m-1}\}{a1​,a2​,...,am−1​}和序列{b1,b2,...,bn}\{b_1,b_2,...,b_{n}\}{b1​,b2​,...,bn​}的最长公共子序列
• 如果$a_m\ne b_n$，且$c_k\ne b_n$，则$c_k=a_m$，序列{c1,c2,...,ck}\{c_1,c_2,...,c_{k}\}{c1​,c2​,...,ck​}是序列{a1,a2,...,am}\{a_1,a_2,...,a_{m}\}{a1​,a2​,...,am​}和序列{b1,b2,...,bn−1}\{b_1,b_2,...,b_{n-1}\}{b1​,b2​,...,bn−1​}的最长公共子序列

记$C[i,j]$为序列{a1,a2,...,ai}\{a_1,a_2,...,a_i\}{a1​,a2​,...,ai​}和序列{b1,b2,...,bj}\{b_1,b_2,...,b_j\}{b1​,b2​,...,bj​}的最长公共子序列(LCS)的长度。根据上述推论，可以得到$C[i,j]$的递推公式：

C[i,j]={0i=0或j=0C[i−1,j−1]+1i,j>0,xi=yimax{C[i,j−1],C[i−1,j]}x,j>0,xi≠yj C[i,j]=\left\{ \begin{array}{lcl} 0 & &{i=0 或 j=0}\\ C[i-1,j-1]+1 & &{i,j>0,x_i=y_i}\\ max\{C[i,j-1],C[i-1,j]\} & &{x,j>0,x_i\neq y_j} \end{array} \right. C[i,j]=⎩⎨⎧​0C[i−1,j−1]+1max{C[i,j−1],C[i−1,j]}​​i=0或j=0i,j>0,xi​=yi​x,j>0,xi​​=yj​​

在求最长公共子序列长度矩阵$C$的时候，需要借助一个辅助矩阵$B$来记录求得$C[i,j]$的过程：
$$B[i,j]=\{\begin{array}{lcl}0& & A[i]=B[i]\\ 1& & A[i]\ne B[i],C[i-1,j]>C[i,j-1]\\ 2& & A[i]\ne B[i],C[i-1,j]<C[i,j-1]\\ 3& & A[i]\ne B[i],C[i-1,j]=C[i,j-1]\end{array}$$

在求得最长公共子序列$C$以及辅助矩阵$B$之后，即可使用回溯法来得到所有的最长公共子序列。

3 实现代码

3.1 LCS.h

# ifndef _LCS_H
# define _LCS_H
# include <vector>
using namespace std;

/**
 *  序列类
 */
class Seq{
public:
    //数据指针
    int* data=nullptr;
    //数据长度
    int len=0;
    Seq(int* data,int len){
        this->data=data;
        this->len=len;
    }
};

// LCS类
class LCS{
private:
    //最长公共子序列长度矩阵
    int** C=nullptr;
    //辅助矩阵
    int** src=nullptr;

    //使用动态规划方法计算C矩阵和src矩阵
    void writeTable(Seq* seq1,Seq* seq2);
    //使用回溯方法来得到所有的最长公共子序列
    void createSubSeq(Seq* seq1,vector<int>lcs , int i,int j,int MaxLen,vector<vector<int>>&LcsSet);

public:

    /**
     * 对外的接口
     * 输入是两个序列
     * 输出是0-n条最长公共子序列
     */ 
    vector<vector<int>> getLCS(Seq* seq1,Seq* seq2);

};
#endif

3.2 LCS.cpp

# include "LCS.h"
#include <cstdio>
# include <vector>
#include <algorithm>

using namespace std;


/**
 * 动态规划法求最长公共子序列长度矩阵C以及辅助矩阵
 * 输入是两个序列
*/
void LCS::writeTable(Seq* seq1,Seq* seq2){
    // 此处有一个小技巧，扩展了一行一列，由于后面会出现C[i-1][j] C[i][j-1]，扩展一行一列可以避免下标出现-1
    int len1=seq1->len+1;
    int len2=seq2->len+1;
    //分配内存
    C=new int*[len1];
    src = new int*[len1];
    for(int i=0;i<len1;i++){
        C[i] = new int[len2];
        src[i] = new int[len2];
    }
    //计算C矩阵和辅助矩阵src
    for(int i=0;i<len1;i++){
        for(int j=0;j<len2;j++){
            if(i==0||j==0){ //自行扩展的一行一列，无意义
                C[i][j]=0;   //无实际用处
                src[i][j]=0; //无实际用处
            }else if(seq1->data[i-1]==seq2->data[j-1]){ 
                C[i][j]=C[i-1][j-1]+1;    
                src[i][j]=0; 
            }else{ 
                if(C[i-1][j]>C[i][j-1]){ 
                    C[i][j]=C[i-1][j];
                    src[i][j]=1;
                }else if(C[i-1][j]<C[i][j-1]){
                    C[i][j]=C[i][j-1];
                    src[i][j]=2;
                }else{
                    C[i][j]=C[i][j-1];
                    src[i][j]=3; 
                }
            }
        }
    }    
}

/**
 * 判断子序列是否和已有的子序列重复
 * 输入：LcsSet为已有的子序列集合，lcs为待加入的子序列
 * 输出：若重复返回true，否则false
 */
bool isRepeat(vector<vector<int>>LcsSet,vector<int>lcs){
    int flag=0;
    for(int i=0;i<LcsSet.size();i++){
        if(lcs.size()!=LcsSet.at(i).size()){
            return true;
        }
        for(int j=0;j<lcs.size();j++){
            if(lcs.at(j)!=LcsSet.at(i).at(j)){
               flag=1; 
            }
        }
        if(flag==0){
            return true;
        }
        flag=0;
    }
    return false;
}

/**
 * 使用回溯法来找出所有的最长公共子序列，此处使用辅助矩阵src完成，用到递归编程技巧
 * seq1是第一个序列，与C[i][j]的i对应；lcs是最长公共子序列；i,j是下标；Maxlen是最长公共子序列的长度；LcsSet是最长公共子序列的集合
 */
void LCS::createSubSeq(Seq* seq1,vector<int>lcs , int i,int j,int MaxLen,vector<vector<int>>&LcsSet){
    if(i==0||j==0){
        if(lcs.size()==MaxLen&&MaxLen>0){//已经回溯完成
            reverse(lcs.begin(),lcs.end()); //翻转一下
            if(!isRepeat(LcsSet,lcs)){ //去重
                LcsSet.push_back(lcs); //加入
            }
        }
        return;
    }
    switch (src[i][j]){
    case 0://左上角
        lcs.push_back(seq1->data[i-1]); //seq1对应i
        createSubSeq(seq1,lcs,i-1,j-1,MaxLen,LcsSet);
        break;
    case 1://向上
        createSubSeq(seq1,lcs,i-1,j,MaxLen,LcsSet);
        break;
    case 2://左
        createSubSeq(seq1,lcs,i,j-1,MaxLen,LcsSet);
        break;
    case 3://左右皆可
        createSubSeq(seq1,lcs,i-1,j,MaxLen,LcsSet);
        createSubSeq(seq1,lcs,i,j-1,MaxLen,LcsSet);
        break;
    default:
        break;
    }
}


vector<vector<int>> LCS::getLCS(Seq* seq1,Seq* seq2){
    if(seq1->len==0||seq2->len==0||seq1->data==nullptr||seq2->data==nullptr){
        return vector<vector<int>>();
    }
    writeTable(seq1,seq2);

    vector<int> lcs;
    vector<vector<int>> LcsSet;
    createSubSeq(seq1,lcs,seq1->len,seq2->len,C[seq1->len][seq2->len],LcsSet);

    return LcsSet;
}

3.3 test.cpp

#include "LCS.h"
#include <cstdio>
#include <vector>
using namespace std;


void test1(){
    LCS* lcs = new LCS();
    
    Seq seq1 = Seq(new int[7]{1,8,2,3,7,8,2},7);
    Seq seq2 = Seq(new int[9]{3,5,7,4,8,6,7,8,2},9);
    vector<vector<int>> LcsSet= lcs->getLCS(&seq1,&seq2);
    printf("nums of lcs : %d\n",LcsSet.size());
    for(int i=0;i<LcsSet.size();i++){
        for(int j=0;j<LcsSet.at(i).size();j++){
            printf("%d ",LcsSet.at(i).at(j));
        }
        printf("\n");
    }
}

int main(){
    test1();
    return 0;
}

3.4 Makefile

FLAG= -std=c++11 
LIBS= 

act:test.o  LCS.o
	g++ -o act test.o LCS.o $(FLAG) $(LIBS)


test.o:test.cpp 
	g++ -c test.cpp  $(FLAG)

LCS.o:LCS.cpp LCS.h
	g++ -c LCS.cpp $(FLAG)

4 测试结果