整理：最长公共子串及最长公共子序列

于 2012-06-05 22:03:37 发布
阅读量659 收藏 1
点赞数 1
文章标签： distance delete insert string c 存储

http://blog.youkuaiyun.com/bertzhang/article/details/7267772

http://www.cnblogs.com/zhangchaoyang/articles/2012070.html



最长公共子串

       找两个字符串的最长公共子串，这个子串要求在原字符串中是连续的。其实这又是一个序贯决策问题，可以用动态规划来求解。我们采用一个二维矩阵来记录中间的结果。这个二维矩阵怎么构造呢？直接举个例子吧："bab"和"caba"(当然我们现在一眼就可以看出来最长公共子串是"ba"或"ab")

　　 b　　a　　b

c　　0　　0　　0

a　　0　　1　　0

b　　1　　0　　1

a　　0　　1　　0

我们看矩阵的斜对角线最长的那个就能找出最长公共子串

不过在二维矩阵上找最长的由1组成的斜对角线也是件麻烦费时的事，下面改进：当要在矩阵是填1时让它等于其左上角元素加1。

　　 b　　a　　b

c　　0　　0　　0

a　　0　　1　　0

b　　1　　0　　2

a　　0　　2　　0

这样矩阵中的最大元素就是 最长公共子串的长度。

在构造这个二维矩阵的过程中由于得出矩阵的某一行后其上一行就没用了，所以实际上在程序中可以用两个一维数组(cur,temp)来代替这个矩阵。

 
#include<stdio.h>
#include<stdlib.h>
#include<string.h>

#define MAX 100
char* LCS(const char* s1,const char* s2);
int main()
{
	char s1[MAX],s2[MAX];

	scanf("%s",s1);
	scanf("%s",s2);
	printf("%s",LCS(s1,s2));
	return 0;
}

char* LCS(const char* s1,const char* s2)
{
	int len1=strlen(s1);
	int len2=strlen(s2);
	int i,j,maxlen=0,lastpos;

	int* cur=(int *)malloc(sizeof(int)*(len1+1));
	int* temp=(int*)malloc(sizeof(int)*(len1+1));
	char* result;

	for(i=0;i<len1;i++)
	{
		cur[i]=0;
		temp[i]=0;
	}
	for(j=0;j<len2;j++)
	{
		for(i=0;i<len1;i++)
		{
			if(s2[j]==s1[i])
			{
				if(i>0)
					cur[i]=temp[i-1]+1;
				else
					cur[i]=1;
			}

			if(cur[i]>maxlen)
			{
				maxlen=cur[i];
				lastpos=j;
			}
		}
		for(i=0;i<len1;i++)
			temp[i]=cur[i];
	}
	result=(char*)malloc(sizeof(char)*(maxlen+1));
	result[maxlen]='\0';
	i=maxlen-1;j=lastpos;
	while(i>=0)
	{
		result[i]=s2[j];
		i--;
		j--;
	}

	return result;

}

最长公共子序列

最长公共子序列与最长公共子串的区别在于最长公共子序列不要求在原字符串中是连续的，比如ADE和ABCDE的最长公共子序列是ADE。

我们用动态规划的方法来思考这个问题如是求解。首先要找到状态转移方程：

等号约定，C1是S1的最右侧字符，C2是S2的最右侧字符，S1‘是从S1中去除C1的部分，S2'是从S2中去除C2的部分。

LCS(S1,S2)等于下列3项的最大者：

（1）LCS（S1，S2’）

（2）LCS（S1’，S2）

（3）LCS（S1’，S2’）--如果C1不等于C2； LCS（S1'，S2'）+C1--如果C1等于C2；

边界终止条件：如果S1和S2都是空串，则结果也是空串。

下面我们同样要构建一个矩阵来存储动态规划过程中子问题的解。这个矩阵中的每个数字代表了该行和该列之前的LCS的长度。与上面刚刚分析出的状态转移议程相对应，矩阵中每个格子里的数字应该这么填，它等于以下3项的最大值：

（1）上面一个格子里的数字

（2）左边一个格子里的数字

（3）左上角那个格子里的数字（如果 C1不等于C2）； 左上角那个格子里的数字+1（ 如果C1等于C2）

举个例子：

 　　　　  G　　C　　T　　A

　　 0　　0　　0　　0　　0

G　 0　　1　　1　　1　　1

B　 0　　1　　1　　1　　1

T　 0　　1　　1　　2　　2

A    0　　1　　1　　2　　3

填写最后一个数字时，它应该是下面三个的最大者：

（1）上边的数字2

（2）左边的数字2

（3）左上角的数字2+1=3,因为此时C1==C2

所以最终结果是3。

在填写过程中我们还是记录下当前单元格的数字来自于哪个单元格，以方便最后我们回溯找出最长公共子串。有时候左上、左、上三者中有多个同时达到最大，那么任取其中之一，但是在整个过程中你必须遵循固定的优先标准。在我的代码中优先级别是左上>左>上。

下图给出了回溯法找出LCS的过程：

奉上代码：




#include<iostream>
#include<cstring>
#include<stack>
#include<utility>
#define LEFTUP 0
#define LEFT 1
#define UP 2
using  namespace  std;
int  Max( int  a, int  b, int  c, int  *max){             //找最大者时a的优先级别最高，c的最低.最大值保存在*max中
     int  res=0;           //res记录来自于哪个单元格
     *max=a;
     if (b>*max){
         *max=b;
         res=1;
     }
     if (c>*max){
         *max=c;
         res=2;
     }
     return  res;
}
//调用此函数时请注意把较长的字符串赋给str1,这主要是为了在回溯最长子序列时节省时间。如果没有把较长的字符串赋给str1不影响程序的正确执行。
string LCS( const  string &str1, const  string &str2){
     int  xlen=str1.size();                //横向长度
     int  ylen=str2.size();                //纵向长度
     if (xlen==0||ylen==0)                 //str1和str2中只要有一个为空，则返回空
         return  "" ;
     pair< int , int > arr[ylen+1][xlen+1];     //构造pair二维数组，first记录数据，second记录来源
     for ( int  i=0;i<=xlen;i++)          //首行清0
         arr[0][i].first=0;
     for ( int  j=0;j<=ylen;j++)          //首列清0
         arr[j][0].first=0;
     for ( int  i=1;i<=ylen;i++){
         char  s=str2.at(i-1);
         for ( int  j=1;j<=xlen;j++){
             int  leftup=arr[i-1][j-1].first;
             int  left=arr[i][j-1].first;
             int  up=arr[i-1][j].first;
             if (str1.at(j-1)==s)          //C1==C2
                 leftup++;
             int  max;
             arr[i][j].second=Max(leftup,left,up,&arr[i][j].first);
//          cout<<arr[i][j].first<<arr[i][j].second<<" ";
         }
//      cout<<endl;
     }        /*矩阵构造完毕*/
     //回溯找出最长公共子序列
     stack< int > st;
     int  i=ylen,j=xlen;
     while (!(i==0&&j==0)){
         if (arr[i][j].second==LEFTUP){
             if (arr[i][j].first==arr[i-1][j-1].first+1)
                 st.push(i);
             --i;
             --j;
         }
         else  if (arr[i][j].second==LEFT){
             --j;
         }
         else  if (arr[i][j].second==UP){
             --i;
         }
     }
     string res= "" ;
     while (!st.empty()){
         int  index=st.top()-1;
         res.append(str2.substr(index,1));
         st.pop();
     }
     return  res;
}
int  main(){
     string str1= "GCCCTAGCG" ;
     string str2= "GCGCAATG" ;
     string lcs=LCS(str1,str2);
     cout<<lcs<<endl;
     return  0;
}

扩展

对于一个字符串a可以通过增加一个字符、删除一个字符、修改一个字符，将字符串a变成字符串b，例如

a= abcddefg

b = abcefg

可以通过a字符串删除两个dd得到b字符串，也可以通过b字符串增加dd编程a字符串，从上面的分析可以知道，增加和删除的代价必须是相同的，这样a字符串变成b字符串的代价和b字符串变成a字符串的代价才会是相同的，否这可能产生代价不对称的情况。其实我们可以设定修改和增加（删除）的代价是不同的，当然也可以认为他们是一样的。

实际的计算过程可以如下进行：

1)比较a[i]和b[j];

2)如果a[i] == b[j]，那么distance = EditDistance(a[i + i], b[j + 1]) + 0;

3)如果a[i] != b[j]，那么可以经过如下操作使得a[i]等于b[j]

   a) a[i]前增加b[j]，那么distance = EditDistance(a[i], b[j + 1] + insert_cost

   b)b[j]前增加a[i]，那么distance = EditDistance(a[i + 1], b[j]) + insert_cost

   c)删除a[i]，那么distance = EditDistance(a[i + 1], b[j]) + delete_cost

   d)删除b[j]，那么distance = EditDistance(a[i], b[j + 1] + delete_cost

   e)a[i]变成b[j]，那么distance = EditDistance(a[i + 1], b[j + 1] + replace_cost

   f)b[j]变成a[i]，那么distance = EditDistance(a[i + 1], b[j + 1] + replace_cost

如果insert_cost == delete_cost，那么a添加字符变成b和b删除字符变成a是等价的，a[i]变成b[j]与b[j]变成a[i]也是等价的，因此实际需要考虑的代价就是下面3种情况：

i) distance = EditDistance(a[i], b[j + 1] + insert_cost(或delete_cost)

ii) distance = EditDistance(a[i + 1], b[j] + insert_cost(或delete_cost)

iii)distance = EditDistance(a[i + 1], b[j + 1] + replace_cost

如果我们回顾一下最长公共子序列问题（LCS），就会发现这个问题和LCS问题几乎是等价的。因为可以这样理解，找出a和b的LCS，保持LCS对齐不变，增加删除一些字符就完成了变换，而这样的代价应该是最小的（猜测的，没有证明）

这样一个问题，我们可以使用递归来解决。

当然的解决方案是用动态规划的方法解决，采用动态规划解决时，我们假设前面字符串都已经变换相同了，那么在a[i]变成b[j]的过程中需要对比如下代价：

1)如果a[i] == b[j]，那么前面的状态可能是a[i - 1] b[j - 1]或者 a[i - 1] b[j]或者a[i] b[j - 1]，我们要比较这些可能的转换过程中哪个代价更小；

2)如果a[i] != b[j]，那么：

   i)a[i] b[j]可能从a[i - 1] b[j - 1]状态通过replace a[i] to a[j] + replace的代价来实现；

  ii)a[i] b[j]可能从a[i ] b[j - 1]状态通过为a[i]前面添加一个b[j -1] + insert的代价来实现；

  iii)a[i] b[j]可能从a[i  - 1] b[j]状态通过删除a[i - 1] + delete的代价来实现；

而这些代价就通过一个向量cost向量来存储。

程序代码如下：


[cpp]  view plain copy
  1. #include <stdio.h>  
  2. #include <string>  
  3.   
  4. int Min(int a, int b, int c) {  
  5.   int tmp = a > b ? b : a;  
  6.   tmp = tmp > c ? c : tmp;  
  7.   return tmp;  
  8. }  
  9. int EditDistance(const std::string& a, int a_offset, const std::string& b, int b_offset) {  
  10.   if (a_offset == a.size() && b_offset < b.size()) {  
  11.     return EditDistance(a, a_offset, b, b_offset + 1) + 1 ;  
  12.   } else if (b_offset == b.size() && a_offset < a.size()) {  
  13.     return EditDistance(a, a_offset + 1, b, b_offset) + 1;  
  14.   } else if (a_offset == a.size() && b_offset == b.size()) {  
  15.     return 0;  
  16.   } else {  
  17.     if (a[a_offset] == b[b_offset]) {  
  18.       return EditDistance(a, a_offset + 1, b, b_offset + 1);  
  19.     } else {  
  20.       int distance1 = EditDistance(a, a_offset + 1, b, b_offset) + 1;  
  21.       int distance2 = EditDistance(a, a_offset, b, b_offset + 1) + 1;  
  22.       int distance3 = EditDistance(a, a_offset + 1, b, b_offset + 1) + 1;  
  23.       return Min(distance1, distance2, distance3);  
  24.     }  
  25.   }  
  26. }  
  27. int EditDistance_DP(const std::string& a, const std::string& b) {  
  28.   int** cost = new int*[a.size() + 1];  
  29.   for (int i = 0; i < a.size() + 1; ++i) {  
  30.     cost[i] = new int[b.size() + 1];  
  31.   }  
  32.   for (int i = 0; i < a.size() + 1; ++i) {  
  33.     for (int j = 0; j < b.size() + 1; ++j) {  
  34.       cost[i][j] = 0;  
  35.     }  
  36.   }  
  37.   for (int i = 0; i < a.size(); ++i) {  
  38.     for (int j = 0; j < b.size(); ++j) {  
  39.       if (a[i] == b[j]) {  
  40.         cost[i + 1][j + 1] = Min(cost[i][j], cost[i][j + 1], cost[i + 1][j]);  
  41.       } else {  
  42.         cost[i + 1][j + 1] = Min(cost[i][j] + 1, cost[i][j + 1] + 1, cost[i + 1][j] + 1);  
  43.       }  
  44.     }  
  45.   }  
  46.   for (int i = 0; i <= a.size(); ++i) {  
  47.     for (int j = 0; j <= b.size(); ++j) {  
  48.       printf("%d  ", cost[i][j]);  
  49.     }  
  50.     printf("\n");  
  51.   }  
  52.   int distance = cost[a.size()][b.size()];  
  53.   for (int i = 0; i < a.size() + 1; ++i) {  
  54.     delete[] cost[i];      
  55.   }  
  56.   delete[] cost;  
  57.   return distance;  
  58. }  
  59. int main(int argc, char** argv) {  
  60.   std::string a = "adefk";  
  61.   std::string b = "bdefg";  
  62.   printf("edit distance = %d\n", EditDistance(a, 0, b, 0));  
  63.   printf("edit distance = %d\n", EditDistance_DP(a, b));  
  64. }  
  65.         
sophie_wise8
关注
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值