求所有最大公共子序列的算法實現

作者：DarkHorse 2012-11-16 10:15:12

最近看了很多關于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了。這里給出實現代碼：

[[98078]]View Code

我們給出一個測試：

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棧的狀態如圖所示。虛線框節點表示已彈出，下同。