动态规划——最长公共子序列LCS

【算法导论 ch15-4】

首先确定什么是最长公共子序列：

最长公共子串（longest common substring） ：其中的子串是原串的一个连续的部分；

最长公共子序列（longest common subsequence）：是不改变序列的顺序，从序列中得到任意的元素从而获得的新序列

二者区别=>子串中字符位置必须连续，子序列可以不必连续

（由于序列元素顺序的不同，得到的LCS可能不止一个）

用动态规划方法求解，需要先寻找最优子结构

定义c[i,j]为序列Xi和Yj的LCS长度，递归求解，将其保存到二维数组C中。可得递归式：

这样得到的表中最后一个值就是Xm和Yn的LCS；

还需要一个二维数组b，用于保存c[i,j]中的值是从c[i-1,j-1],c[i-1,j]还是c[i,j-1]中计算得到的，便于反向查找LCS元素时找到元素位置

伪代码：

便于反向查找LCS元素时找到元素位置

算导上的数组示意图：

构造LCS序列时，用b数组，从b[m,n]开始，延箭头反向追踪下去，每当遇到斜向尽头都证明该处对应的Xi和Yj是相同的，且是属于LCS的。一直追踪直到i或j到达最小。

为了保存下来每一个LCS，增加了一个数组用于保存结果，完整代码如下：

#include <stdio.h>
#include <string>
#include <math.h>
using namespace std;
/*可以考虑不适用b数组，直接从c[i-1,j-1],c[i-1,j],c[i,j-1]中计算得到c[i,j]运算方向*/
const string str_x = "BDCABA";//"ACCGGTCGAGTGCGCGGAAGCCGGCCGAA";//
const string str_y = "ABCBDAB";//"GTCGTTCGGAATGCCGTTGCTCTGTAAA";//

int **c;
int **b;
char *z; 
/*
建立的二维数组c和b的维数应当分别比序列x和y多1，用于存放初始的0值；
所以在c和b数组中存放的值，对应序列x和y要大1，下标容易弄错

代码输出一种LCS，但LCS不一定只有1种
*/
int count;
int lenLCS;
int getsublength(string x,string y,int m,int n)
{
	int i,j;
	//初始化
	for(i = 0;i<m;i++)
	{
		c[i][0] = 0;
		b[i][0] = 0;
	}
	for(i = 0;i<n;i++)
	{
		c[0][i] = 0;
		b[0][i] = 0;
	}
	//填写表格 ==> 表格c[i,j]是序列Xi和Yj的一个LCS的长度，用于得到LCS最终长度
	//b[i,j] ==>用来记录数组中箭头的方向，便于反向查找LCS元素时找到元素位置
	for (i=1;i<m ;i++ )
	{
		for (j = 1;j<n ;j++ )
		{
			if(x[i-1]==y[j-1])//斜角==3
			{
				c[i][j] = c[i-1][j-1] + 1;//LCS长度加上1
				b[i][j] = 3;
			}
			else if(c[i-1][j]==c[i][j-1])左和上相同==4 ==> 需要当做两种情况处理
			{
				c[i][j] = c[i-1][j];
				b[i][j] = 4;
			}
			else if(c[i-1][j]>c[i][j-1])//上==2
			{
				c[i][j] = c[i-1][j];
				b[i][j] = 2;
			}
			else if(c[i-1][j]<c[i][j-1])//左==1
			{
				c[i][j] = c[i][j-1];
				b[i][j] = 1;
			} 
		}
	}
	return c[m-1][n-1];
}
void getLCS(int x,int y,int k,char **ans)//增加一个数组保存结果,存下每一次运行的结果
{
	if( (x==0||y==0)&& k == 1 ) 
	{
		int t;
		count++;
		for (t = 1;t<lenLCS ;t++ )
		{
			//printf("%c",z[t]);
			ans[count][t] = z[t];
		}
		//printf("%d",k);
		//printf("\n");
		return;
	}
	if(b[x][y] == 4)
	{
		//
		getLCS(x-1,y,k,ans);
		getLCS(x,y-1,k,ans);
	}
	else if (b[x][y] == 3)
	{
		k--;
		z[k] = str_x[x-1];
		getLCS(x-1,y-1,k,ans);
	}
	else if(b[x][y] == 2) 
		getLCS(x-1,y,k,ans);
	else if(b[x][y] == 1) 
		getLCS(x,y-1,k,ans); 
	//return;
}
void printarray(int **A,int x,int y)
{
	int i,j;
	for (i = 0;i<x ;i++ )
	{
		for (j = 0;j<y ;j++ )
		{
			printf ("%3d",A[i][j]);
		}
		printf("\n");
	}
}
int main()
{
	int i,j;
	count = 0;
	int len_x = str_x.length();
	int len_y = str_y.length();

	c = new int* [len_x+1];
	b = new int *[len_x+1];
	for (i = 0;i<len_x+1 ;i++ )
	{
		c[i] = new int [len_y+1];
		b[i] = new int [len_y+1];
	}
	int sublen = getsublength(str_x,str_y,len_x+1,len_y+1);
	lenLCS = sublen + 1;
	printf("length of LCS = %d\n",sublen);

	printarray(c,len_x+1,len_y+1);
	printf("\n");
	printarray(b,len_x+1,len_y+1);

	z = new char [sublen+1];
	char** ans = new char *[len_x+1];
	for (i = 0;i<len_x ;i++ )
	{
		ans[i] = new char [sublen+1];
	}
	getLCS(len_x,len_y,sublen+1,ans);
	/*
	for (i = 1;i<sublen+1 ;i++ )
	{
		printf("%c",z[i]);
	}
	printf("\n");
	*/
	printf("\n%d kinds of LCS:\n",count);
	for (i = 1;i<=count ; i++)
	{
		for (j = 1;j<=sublen ;j++ )
		{
			printf("%c",ans[i][j]);
		}
		printf("\n");
	}
	/* 
	printf("%s\n",z);*/
	delete [] z;
	for (i=0;i<len_x+1 ;i++)
	{
		delete [] c[i];
		delete [] b[i];
		delete [] ans[i];
	}
	delete [] c;
	delete [] b;
	delete [] ans;
}