最长公共子序列LCS 问题的matlab实现代码

关于输出多个LCS（最长公共⼦序列）的简单技巧百度百科： ⼀个序列 S ，如果分别是两个或多个已知序列的⼦序列， 且是所有符合此条件序列中最长的，则 S 称为已知序列的最长公共⼦序列。

注意：S在已知序列中可以不连续；⽐如ABCBDAB和BDCABA的LCS为BCBA，BCBA不连续出现。

LCS通常利⽤动态规划算法求解，最简单的求解⽅式如《算法导论》所述： 1.设输⼊为X1X2…X m，Y1Y2…Y n，构造⼆维数组B、C； 其中B[i][j]存储C[i][j]的最优解指向，C[i][j]存储X[1…i]和Y[1…j]的LCS长度； 2.B、C数组同时更新： 2.1.令B[0][k]、C[0][k]（k=0…n）和B[k][0]、C[k][0]（k=0…m）为0； 2.2.按下标从⼩到⼤的顺序依次计算每⾏的B[i][j]和C[i][j]，规则如下： ①.若X[i]=Y[j]，则C[i][j]=C[i-1][j-1]+1，B[i][j]=0； ②.若X[i]！=Y[j]，则C[i][j]=Max(C[i][j-1],C[i-1][j])，B[i][j]=1或2； 3.更新完B、C后，从C[m][n]开始递归输出，根据B[i][j]的指向打印LCS；相信⼤家都明⽩这个原理，简单的⽤⼀张图来说明⼀下： 1.箭头代表B[i][j]的值，其中⽤指向左上⽅的箭头代表0，←和↑代表1和2(或2和1)； 2.箭头下⽅的数字为C[i][j]的值，即⼦问题X[1…i]和Y[1…j]的LCS长度；可是这样的构造只能输出⼀个LCS，但是我们发现BCBA、BCAB、BDAB都是该序列集的LCS。

为什么呢？因为该算法在C[i-1][j]=C[i][j-1]时，只记录了⼀种情况！即没有保存分⽀！如何简单的修改算法，以求得多个LCS呢？下⾯介绍⼀种改进的⽅法。

还是利⽤数组B、C，并且C的计算规则不变，下⾯修改B的计算规则： 1.若X[i]=Y[j]，B[i][j]=0； 2.若X[i]!=Y[j]，有以下判断： ①.若C[i-1][j]>C[i][j-1]，则C[i][j]=C[i-1][j]，B[i][j]=1； ②.若C[i-1][j]<C[i][j-1]，则C[i][j]=C[i][j-1]，B[i][j]=2； ③.若C[i-1][j]=C[i][j-1]，则C[i][j]的长度任取，B[i][j]=3；此时⽤另⼀张图表⽰该算法： 1.⽤四种箭头表⽰B[i][j]的值，其中←和↑都存在时，表⽰C[i-1][j]=C[i][j-1]； 2.输出的时候采⽤正则表达式的思想：遇到B[i][j]=3时，同时对两种情况进⾏递归； 形式为“（←+↑）”，其中←和↑分别表⽰两个⼦问题的LCS；补充：（A+B），在正则表达式中意义为A或B，即“+”表⽰⼆选⼀的关系。

最长公共子序列问题一．实验目的：1.加深对最长公共子序列问题算法的理解，实现最长公共子序列问题的求解算法；2.通过本次试验掌握将算法转换为上机操作；3.加深对动态规划思想的理解，并利用其解决生活中的问题。

二．实验内容：1.编写算法：实现两个字符串的最长公共子序列的求解；2.将输入与输出数据保存在文件之中，包括运行时间和运行结果；3.对实验结果进行分析。

三．实验操作：1.最长公共子序列求解：将两个字符串放到两个字符型数组中，characterString1和characterString2，当characterString1[m]= characterString2[m]时，找出这两个字符串m之前的最长公共子序列，然后在其尾部加上characterString1[m]，即可得到最长公共子序列。

当characterString1[m] ≠characterString2[m]时，需要解决两个子问题：即找出characterString1(m-1)和characterString2的一个最长公共子序列及characterString1和characterString2(m-1)的一个最长公共子序列，这两个公共子序列中较长者即为characterString1和characterString2的一个最长公共子序列。

2.动态规划算法的思想求解：动态规划算法是自底向上的计算最优值。

计算最长公共子序列长度的动态规划算法LCS-Length以characterString1和characterString2作为输入，输出两个数组result和judge1，其中result存储最长公共子序列的长度，judge1记录指示result的值是由那个子问题解答得到的，最后将最终的最长公共子序列的长度记录到result中。

以LCS-Length计算得到的数组judge1可用于快速构造序列最长公共子序列。

首先从judge1的最后开始，对judge1进行配对。

问题描述：一个给定序列的子序列是在该序列中删去若干元素后得到的序列。

确切地说，若给定序列X= { x1, x2,…, x m}，则另一序列Z= {z1, z2,…, z k}是X的子序列是指存在一个严格递增的下标序列{i1, i2,…, i k}，使得对于所有j=1,2,…,k有X ij=Z j。

例如，序列Z={B,C,D,B}是序列X={A,B,C,B,D,A,B}的子序列，相应的递增下标序列为{2,3,5,7}。

给定两个序列X和Y，当另一序列Z既是X的子序列又是Y的子序列时，称Z是序列X和Y的公共子序列。

例如，若X= { A, B, C, B, D, A, B}和Y= {B, D, C, A, B, A}，则序列{B,C,A}是X和Y的一个公共子序列，序列{B,C,B,A}也是X和Y的一个公共子序列。

而且，后者是X和Y的一个最长公共子序列，因为X和Y没有长度大于4的公共子序列。

给定两个序列X= {x1, x2, …, x m}和Y= {y1, y2, … , y n}，要求找出X和Y的一个最长公共子序列。

问题解析：设X= { A, B, C, B, D, A, B},Y= {B, D, C, A, B, A}。

求X，Y的最长公共子序列最容易想到的方法是穷举法。

对X的多有子序列，检查它是否也是Y的子序列，从而确定它是否为X和Y的公共子序列。

由集合的性质知，元素为m的集合共有2^m个不同子序列，因此，穷举法需要指数级别的运算时间。

进一步分解问题特性，最长公共子序列问题实际上具有最优子结构性质。

设序列X={x1,x2,……x m}和Y={y1,y2,……y n}的最长公共子序列为Z={z1,z2,……z k}。

则有：(1)若x m=y n,则z k=x m=y n,且z k-1是X m-1和Y n-1的最长公共子序列。

(2)若x m!=y n且z k!=x m，则Z是X m-1和Y的最长公共子序列。

// KSY.cpp : 定义控制台应用程序的入口点。

//#include "stdafx.h"#include <iostream>using namespace std;void LCSLength(intm,intn,char *x ,char *y, int **c, int **b) {inti ,j;for (i = 1; i<= m; i++) c[i][0] = 0;for (i = 1; i<= n; i++) c[0][i] = 0;for (i = 1; i<= m; i++)for (j = 1; j <= n; j++){if (x[i]==y[j]){c[i][j]=c[i-1][j-1]+1;b[i][j]=1;}else if (c[i-1][j]>=c[i][j-1]){c[i][j]=c[i-1][j];b[i][j]=2;}else{c[i][j]=c[i][j-1];b[i][j]=3;}}}void LCS(inti ,int j, char *x ,int **b){if (i ==0 || j==0) return;if (b[i][j]== 1){LCS(i-1,j-1,x,b);printf("%c",x[i]);}else if (b[i][j]== 2)LCS(i-1,j,x,b);else LCS(i,j-1,x,b);}constint M = 6;constint N = 5;void output(char *s,int n);void LCSLength(intm,intn,char *x,char *y,int * *c,int * *b); void LCS(inti,intj,char *x,int * *b);void main(){char x[] = {' ','B','C','E','F','G','T'};char y[] = {' ','C','D','F','J','G'};int **c = new int *[M+1];int **b = new int *[M+1];for(inti=0;i<=M;i++){c[i] = new int[N+1];b[i] = new int[N+1];}cout<<"序列X："<<endl;output(x,M);cout<<"序列Y："<<endl;output(y,N);LCSLength(M,N,x,y,c,b);cout<<"序列X、Y最长公共子序列长度为："<<c[M][N]<<endl; cout<<"序列X、Y最长公共子序列为："<<endl;LCS(M,N,x,b);cout<<endl;}void output(char *s,int n) {for(inti=1; i<=n; i++){cout<<s[i]<<" ";}cout<<endl;}。

最长公共子序列算法利用的算法最长公共子序列（LCS）算法利用了递归算法和动态规划算法。

递归算法的代码如下：```c#include<stdio.h>#include<string.h>char a[30], b[30];int m, n;int lcs(int i, int j) {if (i >= m || j >= n) //限定跳出条件return 0;if (a[i] == b[j]) //相等时都往后移动一位return lcs(i+1, j+1) + 1;else //比较得到较大的那一个return lcs(i +1, j) > lcs(i, j + 1)? lcs(i + 1, j) : lcs(i, j + 1);}int main() {strcpy(a, "CNBLOG");strcpy(b, "BELONG");//利用strcpy函数存储字符串m = strlen(a);n = strlen(b);printf("%d\n", lcs(0,0));return 0;}```递归算法的思路是利用一层层往上找，用栈存储，但随着数据规模的增加，算法的复杂度会以指数级增长，因此一般不采用递归算法。

动态规划算法采用二维数组来标识中间计算结果，避免重复的计算来提高效率。

设有字符串a(0…n)，b(0…m)，递推公式如下：字符串a对应的是二维数组num的行，字符串b对应的是二维数组num的列。

另外，采用二维数组flag来记录下标i和j的走向。

数字”1”表示斜向下，数字”2”表示水平向右，数字”3”表示竖直向下。

这样便于以后的求解最长公共子序列。

动态规划经典——最长公共⼦序列问题（LCS）和最长公共⼦串问题⼀.最长公共⼦序列问题(LCS问题)给定两个字符串A和B，长度分别为m和n，要求找出它们最长的公共⼦序列，并返回其长度。

例如： A = "Hel lo W o rld" B = "loo p"则A与B的最长公共⼦序列为 "loo",返回的长度为3。

此处只给出动态规划的解法：定义⼦问题dp[i][j]为字符串A的第⼀个字符到第 i 个字符串和字符串B 的第⼀个字符到第 j 个字符的最长公共⼦序列，如A为“app”,B为“apple”，dp[2][3]表⽰ “ap” 和 “app” 的最长公共字串。

注意到代码中 dp 的⼤⼩为 (n + 1) x (m + 1) ，这多出来的⼀⾏和⼀列是第 0 ⾏和第 0 列，初始化为 0，表⽰空字符串和另⼀字符串的⼦串的最长公共⼦序列，例如dp[0][3]表⽰ "" 和“app” 的最长公共⼦串。

当我们要求dp[i][j]，我们要先判断A的第i个元素B的第j个元素是否相同即判断A[i - 1]和 B[j -1]是否相同，如果相同它就是dp[i-1][j-1]+ 1，相当于在两个字符串都去掉⼀个字符时的最长公共⼦序列再加 1；否则最长公共⼦序列取dp[i][j - 1] 和dp[i - 1][j]中⼤者。

所以整个问题的初始状态为：dp[i][0]=0,dp[0][j]=0相应的状态转移⽅程为：dp[i][j]=max{dp[i−1][j],dp[i][j−1]},A[i−1]!=B[j−1] dp[i−1][j−1]+1,A[i−1]==B[j−1]代码的实现如下：class LCS{public:int findLCS(string A, int n, string B, int m){if(n == 0 || m == 0)//特殊输⼊return 0;int dp[n + 1][m + 1];//定义状态数组for(int i = 0 ; i <= n; i++)//初始状态dp[i][0] = 0;for(int i = 0; i <= m; i++)dp[0][i] = 0;for(int i = 1; i <= n; i++)for(int j = 1; j<= m; j++){if(A[i - 1] == B[j - 1])//判断A的第i个字符和B的第j个字符是否相同dp[i][j] = dp[i -1][j - 1] + 1;elsedp[i][j] = max(dp[i - 1][j],dp[i][j - 1]);}return dp[n][m];//最终的返回结果就是dp[n][m]}};该算法的时间复杂度为O(n*m)，空间复杂度为O(n*m)。

首先将要看到如何运用动态编程查找两个DNA 序列的最长公共子序列（longest common subsequence，LCS）。

发现了新的基因序列的生物学家通常想知道该基因序列与其他哪个序列最相似。

查找LCS 是计算两个序列相似程度的一种方法：LCS 越长，两个序列越相似。

子序列中的字符与子字符串中的字符不同，它们不需要是连续的。

例如，ACE是ABCDE的子序列，但不是它的子字符串。

请看下面两个DNA 序列：∙S1 = DE>GCCCTAGCGDE>∙S2 = DE>GCGCAATGDE>这两个序列的LCS 是GCCAG。

（请注意，这仅是一个LCS，而不是唯一的LCS，因为可能存在其他长度相同的公共子序列。

这种最优化问题和其他最优化问题的解可能不止一个。

）LCS 算法首先，考虑如何递归地计算LCS。

令：∙C1是S1最右侧的字符∙C2是S2最右侧的字符∙S1'是S1中“切掉” C1的部分∙S2'是S2中“切掉” C2的部分有三个递归子问题：∙L1 = LCS(S1', S2)∙L2 = LCS(S1, S2')∙L3 = LCS(S1', S2')结果表明（而且很容易使人相信）原始问题的解就是下面三个子序列中最长的一个：∙L1∙L2∙如果C1等于C2，则为L3后端加上C1，如果C1不等于C2，则为L3。

（基线条件（base case）是S1或S2为长度为零的字符串的情形。

在这种情况下，S1和S2的LCS 显然是长度为零的字符串。

）但是，就像计算斐波纳契数的递归过程一样，这个递归解需要多次计算相同的子问题。

可以证明，这种递归解法需要耗费指数级的时间。

相比之下，这一问题的动态编程解法的运行时间是Θ(mn)，其中m和n分别是两个序列的长度。

为了用动态编程有效地计算LCS，首先需要构建一个表格，用它保存部分结果。

沿着顶部列出一个序列，再沿着左侧从上到下列出另一个序列，如图2 所示：图2. 初始LCS 表格这种方法的思路是：将从上向下、从左到右填充表格，每个单元格包含一个数字，代表该行和该列之前的两个字符串的LCS 的长度。

求解所有最长公共子序列一、问题分析该部分思路同课件二、算法设计思路根据问题分析的结果，具体的算法设计思路如下：1)申明两个数组，用于保存比较的两个字符串；由于事先不知字符串大小，故动态的实现，这里用C++的容器。

2)申明全局变量，二维数组B和数组C。

数组C用于保存计算Xi和Yi的LCS值；数组B保存当前的C是从哪个子问题得来的。

为此，定义一个枚举类型，用于标识不同的方向，分别为对角线、向上、向左和向左向上四个方向。

3)根据动态规划，实现一个函数LCS_LENGTH，完成的功能是计算数组B和C。

具体过程是：先是动态申请二维数组B和C，他们的行列长度都增加1，目的就是方便计算。

将C的第0行和第0列都赋上0，即初始化。

开始计算C[i][j]，以行为主，一次计算C的每一个元素，即将两个数组逐一比较。

比较时就有两种情况，分别是若相等时，就将C[i][j]设置成C[i-1][j-1]，同时将B[i][j]设置成DIAGONAL。

若不相等时，比较C[i-1][j] 和C[i][j-1]的值，又有三种情况：一是C[i-1][j] 与C[i][j-1]相等，就随便把某一个赋给C[i][j]，比如C[i-1][j]，B[i][j]设置为UP_LEFT；二是若C[i-1][j] 大于C[i][j-1]，则将C[i-1][j]赋给C[i][j]，并且将B[i][j]设置成UP；最后是若C[i-1][j] 小于C[i][j-1]，则将C[i][j-1]赋给C[i][j]，并且将B[i][j]设置成LEFT。

4)根据第3）步骤的结果，就可以找出所有LCS了。

这里会用到回溯方法，具体实现可以用栈，也可以用递归。

本人使用的是递归，代码简单、易懂。

具体实现方法是：申请一个数组用于保存一个LCS，这个数组会反复使用，因此，一旦找到一个就会立即将它输出。

再设置一个变量curpos标识当前的数组下标，一个变量len 保存当前LCS数组的元素个数。