代码改变世界

求所有最大公共子序列的算法实现

2012-11-15 23:55  钱吉  阅读(10093)  评论(6)  编辑  收藏  举报

最近看了很多关于LCS(Longest common subsequence problem,最长公共子序列)的文章，大部分问题都只是求出最大公共子序列的长度，或者打印处其中的任意一个最大子序列即可，但是如何快速的打印出所有的最大长度子序列？这个问题好像看到的不多。本文给出了传统的DP(dynamic programming，动态规划)算法进行求解的过程，并用c语言实现。另外参考一篇论文实现了其中的一种打印所有最大公共子序列的算法，这个算法比起传统的算法而言，时间复杂度大大降低.

一：LCS解析

首先看下什么是子序列？定义就不写了，直接举例一目了然。如对于字符串：“student”，那么su，sud，sudt等都是它的子序列。它可以是连续的也可以不连续出现，如果是连续的出现，比如stud，一般称为子序列串，这里我们只讨论子序列。

什么是公共子序列？很简单，有两个字符串，如果包含共同的子序列，那么这个子序列就被称为公共子序列了。如“student”和“shade”的公共子序列就有“s”或者“sd”或者“sde”等。而其中最长的子序列就是所谓的最长公共子序列(LCS)。当然，最长公共子序列也许不止一个，比如：“ABCBDAB”和“BDCABA”，它们的LCS为“BCBA”，“BCAB”，“BDAB”。知道了这些概念以后就是如何求LCS的问题了。

通常的算法就是动态规划(DP)。假设现在有两个字符串序列：X={x1,x2,...xi...xm},Y={y1,y2,...yj...yn}。如果我们知道了X={x1,x2,...xi-1}和Y={y1,y2,...yj-1}的最大公共子序列L，那么接下来我们可以按递推的方法进行求解：

1)如果xi==yj,那么{L,xi(或yj)}就是新的LCS了，其长度也是len(L)+1。这个好理解，即序列{Xi,Yj}的最优解是由{Xi-1,Yj-1}求得的。

2)如果xiyj,那么可以转换为求两种情况下的LCS。

A: X={x1,x2,...xi}与Y={y1,y2,...yj-1}的LCS，假设为L1

B: X={x1,x2,...xi-1}与Y={y1,y2,...yj}的LCS，假设为L2

那么xiyj时的LCS=max{L1,L2},即取最大值。同样，实际上序列{Xi,Yj-1}和{Xi-1,Yj}都可以由{Xi-1,Yj-1}的最优解求得。

怎么样，是不是觉得这种方法很熟悉？当前问题的最优解总是包含了一个同样具有最优解的子问题，这就是典型的DP求解方法。好了，直接给出上面文字描述解法中求LCS长度的公式：

这里用一个二维数组存储LCS的长度信息，i，j分别表示两个字符串序列的下标值。这是求最大公共子序列长度的方法，如果要打印出最大公共子序列怎么办？我们还需要另外一个二维数组来保存求解过程中的路径信息，方便最后进行路径回溯，找到LCS。如果看着很含糊，我下面给出其实现过程。

 

二：DP实现

很多博文上面都有，基本上是用两个二维数组c[m][n]和b[m][n]，一个用来存储子字符串的LCS长度，一个用来存储路径方向，然后回溯。

其中二维数组b[i][j]的取值为1或2或3，其中取值为1时说明此时xi=yj，c[i][j]=c[i-1][j-1]+1。如果将二维数组看成一个矩阵，那么此时代表了一个从左上角到右下角的路径。如果取值为2，说明xi≠yj,且c[i][j]=c[i-1][j]，代表了一个从上到下的路径，同理取值为3代表一个从左到右的路径。

最后我们可以根据c[m][n]的值知道最大公共子序列的长度。然后根据b[i][j]回溯，可以打印一条LCS。其中b[i][j]=1的坐标点对应的字符同时在两个序列中出现，所以依次回溯这个二维数组就可以找到LCS了。这里给出实现代码：

 

View Code 
 1 #include <stdio.h>
 2 #include <stdlib.h>
 3 #include <string.h>
 4 #include "stack.h"
 5 #define MAX_LEN 1024
 6 typedef int **Matrix;
 7 
 8 void GetLCSLen(char *str1, char *str2, Matrix pc, Matrix pb, int nrow, int ncolumn);
 9 Matrix GreateMatrix(int nrow, int ncolumn);
10 void DeleteMatrix(Matrix p, int nrow, int ncolumn);
11 void TraceBack(char *str1, Matrix pb, int nrow, int ncolumn);
12 
13 void GetLCSLen(char *str1, char *str2, Matrix pc, Matrix pb, int nrow, int ncolumn)
14 {
15     int i,j;
16     /************initial the edge***************/
17     for(i=0; i<nrow; i++)
18     {
19         pc[i][0] = 0;
20         pb[i][0] = 0;
21     }
22     for(j=0; j<ncolumn; j++)
23     {
24         pc[0][j] = 0;
25         pb[0][j] = 0;
26     }
27     /************DP*****************************/
28     for(i=1; i<nrow; i++)
29     {
30         for(j=1; j<ncolumn; j++)
31         {
32             if(str1[i-1] == str2[j-1])
33             {
34                 pc[i][j] = pc[i-1][j-1] + 1;//由左上节点转移而来
35                 pb[i][j] = 1;//标记为1
36             }
37             else if(pc[i-1][j] >= pc[i][j-1])
38             {
39                 pc[i][j] = pc[i-1][j];//由上节点转移而来
40                 pb[i][j] = 2;//标记为2
41             }
42             else
43             {
44                 pc[i][j] = pc[i][j-1];//由左节点转移而来
45                 pb[i][j] = 3;//标记为2
46             }
47         }
48     }
49 }
50 void TraceBack(char *str1, Matrix pb, int nrow, int ncolumn)
51 {
52     int ntemp;
53     if(str1 == NULL || pb == NULL)
54         return;
55     if(nrow == 0 || ncolumn == 0)
56         return;
57     ntemp = pb[nrow-1][ncolumn-1];
58     switch(ntemp)
59     {
60     case 1:
61         printf("locate:(%d,%d),%4c\n", nrow-1, ncolumn-1, str1[nrow-2]);//打印公共字符，这里下标是nrow-2，因为矩阵的坐标值(i,j)比字符串的实际下标大1
62         TraceBack(str1, pb, nrow-1, ncolumn-1);//向左上角递归
63         break;
64     case 2:
65         TraceBack(str1, pb, nrow-1, ncolumn);//向上方向递归
66         break;
67     case 3:
68         TraceBack(str1, pb, nrow, ncolumn-1);//向左方向递归
69         break;
70     default:
71         break;
72     }
73 }

 

我们给出一个测试：

1     char str1[MAX_LEN] = "BADCDCBA";
2     char str2[MAX_LEN] = "ABCDCDAB";
3 
4     GetLCSLen(str1, str2, C, B, str1len+1, str2len+1);
5     TraceBack(str1, B, str1len+1, str2len+1);

详细的代码见文章结束处给出的链接。本测试output:BDCDB

问题：上面的方法中，我们没有单独考虑c[i-1][j]==c[i][j-1]的情况，所以在回溯的时候打印的字符只是其中一条最大公共子序列，如果存在多条公共子序列的情况下。怎么解决？我们对b[i][j]二维数组的取值添加一种可能，等于4，这代表了我们说的这种多支情况，那么回溯的时候我们可以根据这个信息打印更多可能的选择。这个过程就不写代码了，其实很简单，以下面的路径回溯图举例，你从（8,8)点开始按b[i][j]的值指示的方向回溯，把所有的路径遍历一遍，如果是能达到起点（0,0)的路径，就是LCS了，有多少条打印多少条。可是，

又出现问题了：你发现没有，在回溯路径的时候，如果采用一般的全搜索，会进行了很多无用功。即重复了很多，且会遍历了一些无效路径，因为这些路径最终不会到达终点(0,0)，比如节点(6,3),(7,2),(8,1)。因此加大算法复杂度和时间消耗。那么如何解决？看下面的这个方法，正式进入本文正题。

 路径回溯图：

 加入状态4后的状态图：

三：算法改进

 

上面提到路径回溯过程中，一般的方法就是遍历所有可能的路径，但是一些不可能构成最大公共子序列的跳跃点我们也会去计算。这里先解释下什么叫跳跃点，就是导致公共子序列长度发生变化的节点，即b[i][j]=1对应的节点(i,j)。Ok，接下来的问题是，如何不去考虑这些无效跳跃点，降低算法复杂度？参考论文里提出了这样一种方法：矩形搜索。

 

首先构造两个栈数据结构store和print。故名思议，一个用来储存节点，一个用来打印节点。栈的定义为：

 1 #define MAX_STACK_SIZE 1024
 2 typedef struct _Element
 3 {
 4     int nlcslen;
 5     int nrow;
 6     int ncolumn;
 7 }Element;
 8 typedef struct _Stack
 9 {
10     int top;
11     Element data[MAX_STACK_SIZE];
12 }Stack;

栈使用数组实现，并有一个指向顶点的下标值top。为了初始化需要，先构造了一个虚拟节点virtualnode，指向节点(m,n)的右下角，即(m+1,n+1),这个节点的LCS的长度假设为最大公共子序列长度+1。将虚拟节点压入栈store，然后然后执行下面的算法：

1)栈store为空吗？是的话退出。

2)否则从store栈顶弹出节点。

3)如果这个节点为边界元素(行或列的小标为1),则将此节点压入栈print，打印栈print里面的所有节点(除virtualnode外)。查看此时store栈顶节点的LCS长度，并以这个长度为参考，弹出print栈里面所有LCS长度小于等于这个参考值的节点，跳转到第1步。

4)如果不是边界元素，那么以该节点的左上节点(i-1,j-1)为出发点，沿着该出发点指示的方向，找到第一个跳跃点e1(即b[i][j]==1的点)。途中碰到分支节点(b[i][j]==4的点)时，沿着其上节点继续探索。

5)找到第一个跳跃点以后，重新回到第4步的出发点，沿着该节点指示的方向，找到第二个跳跃点e2。途中碰到分支节点(b[i][j]==4的点)时，沿着其左节点继续探索

6)如果e1和e2节点为同一节点，将该节点压入store栈，回到步骤1)。

7）如果不为同一节点，在e1和e2构成的矩形范围内，搜索出所有的跳跃点，并全部压入store栈，回到步骤1)。

 

不明白？不要紧，我们结合上面的矩阵图一步步按照算法来，看看到底是如何计算的：

第一步：压入虚拟节点6(9,9)到栈store，这里6表示这个节点的LCS长度，（9,9)表示坐标值。 

第二步：store栈不为空，则弹出store栈顶，压入print栈，这时候的两个栈的状态如下面的左图。沿出发点(8,8)出发，这是个分支节点，因为b[8][8]==4,所以选择向上走，搜索到e1跳跃点(7,8)，搜索路径为：(8,8)->(7,8)。然后回到(8,8)找e2点，这时选择向左走，找到e2跳跃点(8,7)。这两个跳跃点不同，所以以e1，e2为对角线方向围成的矩形内搜索所有跳跃点，这里只有e1,和e2本身两个节点，然后将它们压入栈store。此时两个栈的状态见下面的有图。蓝色底的节点表示有store栈弹出然后压到print栈，绿色底表示新压入到store栈的跳跃点，下面所有的图都这样表示。

   

第三步：弹出5(7,8)到print栈，搜索到新的两跳跃节点。

   

第四步：

    

第五步：

 

第六步：

 

第七步：关键步骤来了，因为此时从store栈弹出的节点是边界元素1(1,2)，所以我们打印print栈的所有元素(红色字体节点)，而这些元素恰好构成了一个最长公共子序列(自己揣摩一下)。打印完了以后，我们要对print栈进行清理，去除不需要的节点，按照步骤2，此时store栈顶节点的LCS为1，所以print栈中弹出的节点只有1(1,2)。弹完以后，print栈的状态如图所示。虚线框节点表示已弹出，下同。

  

第八步：继续弹出store栈顶，发现又是边界元素，继续压入print栈，打印print栈。清理print栈。

 

第九步：清理完后，继续步骤2.

 

好了，下面的过程就是重复进行上面的这些步骤了，自己动手画一下就一目了然了。

 

四：代码实现

用c语言实现关键代码：

View Code 
  1 void TraceBack2(char *str1, Matrix pc, Matrix pb, int nrow, int ncolumn);
  2 void PrintStack(Stack *ps, char *str1, int len1);
  3 void SearchE(Matrix pb, int curposx, int curposy, int *eposx, int *eposy, int ntype);
  4 
  5 void PrintStack(Stack *ps, char *str1, int len1)
  6 {
  7     if(ps == NULL || str1 == NULL)
  8         return;
  9     int ntemp = ps->top;
 10     int index = -1;
 11     while((index = ps->data[ntemp].nrow) <= len1)
 12     {
 13         ntemp--;
 14         printf("%2c",str1[index-1]);
 15     }
 16     printf("\n");
 17 }
 18 
 19 void TraceBack2(char *str1, Matrix pc, Matrix pb, int nrow, int ncolumn)
 20 {
 21     if(str1 == NULL || pc == NULL || pb == NULL)
 22         return;
 23     
 24     Stack store, print;//构造两个栈store，print
 25     Element storetop;//store栈的栈顶节点
 26     Element element;//临时变量
 27     Element virtualnode;//虚拟节点
 28     int ntoplen;//保存store栈顶节点的LCS长度
 29     int ex1,ey1,ex2,ey2;//矩形搜索的两个节点的坐标
 30     int i,j;
 31 
 32     InitialStack(&store);//初始化
 33     InitialStack(&print);
 34     virtualnode = CreateElement(pc[nrow-1][ncolumn-1]+1, nrow, ncolumn);
 35     Push(&store, &virtualnode);//压入虚拟节点到store
 36 
 37     while(!IsEmpty(&store))
 38     {
 39         Pop(&store, &storetop);//从栈顶取出一个节点
 40         if(storetop.nrow == 1 || storetop.ncolumn == 1)//如果是边界节点
 41         {
 42             Push(&print, &storetop);
 43             PrintStack(&print, str1, nrow-1);//打印print栈里面除虚拟节点之外的所有节点
 44             GetTop(&store, &element);
 45             ntoplen = element.nlcslen;//当前store的栈顶节点的LCS长度
 46 
 47             /**********弹出print栈中所有LCS长度小于等于ntoplen的节点**************/
 48             while(GetTop(&print, &element) && element.nlcslen<=ntoplen)
 49             {
 50                 Pop(&print, &element);
 51             }
 52         }
 53         else
 54         {
 55             Push(&print, &storetop);
 56             SearchE(pb, storetop.nrow-1, storetop.ncolumn-1, &ex1, &ey1, 0);
 57             SearchE(pb, storetop.nrow-1, storetop.ncolumn-1, &ex2, &ey2, 1/*also other value is ok*/);
 58 
 59             if(ex1 == ex2 && ey1 ==ey2)
 60             {
 61                 element = CreateElement(pc[ex1][ey1], ex1, ey1);
 62                 Push(&store,&element);//压入store栈，回到步骤2
 63             }
 64             else
 65             {
 66                 for(i=ey2; i<=ey1; i++)
 67                     for(j=ex1; j<=ex2; j++)
 68                     {
 69                         if(pb[i][j] == 1)
 70                         {
 71                             element = CreateElement(pc[i][j], i, j);
 72                             Push(&store, &element);
 73                         }
 74                     }
 75             }
 76         }
 77 
 78     }
 79 }
 80 void SearchE(Matrix pb, int curposx, int curposy, int *eposx, int *eposy, int ntype)
 81 {
 82     switch(pb[curposx][curposy])
 83     {
 84     case 1:
 85         *eposx = curposx;
 86         *eposy = curposy;
 87         return;
 88     case 2:
 89         SearchE(pb, curposx-1, curposy, eposx, eposy, ntype);
 90         break;
 91     case 3:
 92         SearchE(pb, curposx, curposy-1, eposx, eposy, ntype);
 93         break;
 94     case 4:
 95         if(ntype == 0)
 96             SearchE(pb, curposx-1, curposy, eposx, eposy, ntype);//搜索e1点，如过碰到分叉点，向上继续搜索
 97         else
 98             SearchE(pb, curposx, curposy-1, eposx, eposy, ntype);//搜索e2点，如过碰到分叉点，向左继续搜索
 99         break;
100     }
101 }

同样的测试，这种算法能打印出全部的最长公共子序列。

五：总结 

该算法能将传统算法的指数级复杂度降低到max{O(mn),O(ck)},k为最大公共子序列的个数。详细证明见论文。因为用两个栈存储了所有有效跳跃点，使得许多重复比较被忽略。栈的有序性也能巧妙的得到任意一条最大公共子序列。

本算法的全部实现代码下载路径：http://pan.baidu.com/share/link?shareid=125428&uk=285541510

欢迎讨论。

参考论文：

《利用矩阵搜索求所有最长公共子序列的算法》，宫洁卿，安徽工程科技学院学报，vol23,No.4,Dec.,2008