【ZH奶酪】如何用Python計算最長公共子序列和最長公共子串

1. 什麼是最長公共子序列？什麼是最長公共子串？

1.1. 最長公共子序列（Longest-Common-Subsequences，LCS）

最長公共子序列（Longest-Common-Subsequences，LCS）是一個在一個序列集合中（一般爲兩個序列）用來查找全部序列中最長子序列的問題。這與查找最長公共子串的問題不一樣的地方是：子序列不須要在原序列中佔用連續的位置 。

最長公共子序列問題是一個經典的計算機科學問題，也是數據比較程序，好比Diff工具，和生物信息學應用的基礎。它也被普遍地應用在版本控制，好比Git用來調和文件之間的改變。

1.2 最長公共子串（Longest-Common-Substring，LCS）

最長公共子串（Longest-Common-Substring，LCS）問題是尋找兩個或多個已知字符串最長的子串。此問題與最長公共子序列問題的區別在於子序列沒必要是連續的，而子串卻必須是連續的。

2. 如何求解最長公共子序列?

例如序列str_a=world，str_b=wordl。序列wo是str_a和str_b的一個公共子序列，可是不是str_a和str_b的最長公共子序列，子序列word是str_a和str_b的一個LCS，序列worl也是。

暴力查找？
尋找LCS的一種方法是枚舉X全部的子序列，而後注意檢查是不是Y的子序列，並隨時記錄發現的最長子序列。假設X有m個元素，則X有2^m個子序列，指數級的時間，對長序列不實際。

分析問題，設str_a=<x1,x2,…,xm>和str_b=<y1,y2,…,yn>爲兩個序列，LCS(str_a,str_b)表示str_a和str_b的一個最長公共子序列，能夠看出

若是str_a[m] == str_b[n]，則LCS (str_a, str_b) = str_a[m] + LCS(str_a[1:m-1],str_b[1:n-1]) 若是str_a[m] != str_b[n]，則LCS(str_a,str_b)= max{LCS(str_a[1:m-1], str_b), LCS (str_a, str_b[n-1])}

LCS問題也具備重疊子問題性質：爲找出LCS(str_a,str_b)，可能須要找LCS(str_a[1:m-1], str_b)以及LCS (str_a, str_b[n-1])。但這兩個子問題都包含着LCS(str_a[1:m-1],str_b[1:n-1]).

2.1 基於遞歸的方法

根據上邊分析結果，能夠寫出簡潔易懂的遞歸方法。

def recursive_lcs(str_a, str_b):
  if len(str_a) == 0 or len(str_b) == 0:
    return 0
  if str_a[0] == str_b[0]:
    return recursive_lcs(str_a[1:], str_b[1:]) + 1
  else:
    return max([recursive_lcs(str_a[1:], str_b), recursive_lcs(str_a, str_b[1:])])
print recursive_lcs(str_a, str_b)

2.2 基於自底向上動態規劃的方法

根據上述分析問題，動態規劃遞推公式也很是明顯，能夠寫出動態規劃代碼：

def bottom_up_dp_lcs(str_a, str_b):
  """
  longest common subsequence of str_a and str_b
  """
  if len(str_a) == 0 or len(str_b) == 0:
    return 0
  dp = [[0 for _ in range(len(str_b) + 1)] for _ in range(len(str_a) + 1)]
  for i in range(1, len(str_a) + 1):
    for j in range(1, len(str_b) + 1):
      if str_a[i-1] == str_b[j-1]:
        dp[i][j] = dp[i-1][j-1] + 1
      else:
        dp[i][j] = max([dp[i-1][j], dp[i][j-1]])
  print "length of LCS is :",dp[len(str_a)][len(str_b)]
  # 輸出最長公共子序列
  i, j = len(str_a), len(str_b)
  LCS = ""
  while i > 0 and j > 0:
    if str_a[i-1] == str_b[j-1] \   # 這裏必定要比較a[i-1]和b[j-1]是否相等
        and dp[i][j] == dp[i-1][j-1] + 1:
      LCS = str_a[i - 1] + LCS
      i, j = i-1, j-1
      continue
    if dp[i][j] == dp[i-1][j]:
      i, j = i-1, j
      continue
    if dp[i][j] == dp[i][j-1]:
      i, j = i, j-1
      continue
  print "LCS is :", LCS
bottom_up_dp_lcs(str_a, str_b)

2.3 下降空間複雜度的動態規劃算法

根據上述問題分析以及2.2中的dp矩陣能夠看出，其實每一步的求解，只和三個元素有關：左邊的元素，上邊的元素，左上角的元素。所以咱們能夠進行空間優化，用一維數組代替二維矩陣。

def space_efficient_lcs(str_a, str_b):
  """
  longest common subsequence of str_a and str_b, with O(n) space complexity
  """
  if len(str_a) == 0 or len(str_b) == 0:
    return 0
  dp = [0 for _ in range(len(str_b) + 1)]
  for i in range(1, len(str_a) + 1):
    left_up = 0
    dp[0] = 0
    for j in range(1, len(str_b) + 1):
      left = dp[j-1]
      up = dp[j]
      if str_a[i-1] == str_b[j-1]:
        dp[j] = left_up + 1
      else:
        dp[j] = max([left, up])
      left_up = up
  print dp[len(str_b)]
space_efficient_lcs(str_a, str_b)

3. 如何求解最長公共子串？

最長公共子串比最長公共子序列的遞推公式要簡單一些。

dp[i][j]的含義也發生了變化：

• 在最長公共子序列中，dp[i][j]表示str_a[1:i]和str_b[1:j]的最長公共子序列，是從str_a的1和str_b的1開始計算的，即整個字符串的起始位置。
• 在最長公共子串中，dp[i][j]表示str_a[i':i]和str_b[j':j]的最長公共子串，由於str_a和str_b可能存在多個公共子串，因此i'和j'分別表示當前公共子串的起始位置。

也就是說：

• 當str_a[i] == str_b[j]時，dp[i][j] = dp[i-1][j-1]+ 1；
• 當str_a[i] != str_b[j]時，dp[i][j] = 0，即開始計算新的公共子串。

和最長公共子序列不一樣的是，在最長公共子串問題中，dp[m][n]不必定是最終結果，好比「abcdxy」和「abcfxy」，dp[m][n]存儲的是公共子串「xy」的長度，而不是公共子串「abc」的長度，因此須要一個變量單獨記錄最長子串的長度。

3.1 動態規劃算法

def bottom_up_dp_lcs(str_a, str_b):
  """
  longest common substring of str_a and str_b
  """
  if len(str_a) == 0 or len(str_b) == 0:
    return 0
  dp = [[0 for _ in range(len(str_b) + 1)] for _ in range(len(str_a) + 1)]
  max_len = 0
  lcs_str = ""
  for i in range(1, len(str_a) + 1):
    for j in range(1, len(str_b) + 1):
      if str_a[i-1] == str_b[j-1]:
        dp[i][j] = dp[i-1][j-1] + 1
        max_len = max([max_len, dp[i][j]])
        if max_len == dp[i][j]:
          lcs_str = str_a[i-max_len:i]
      else:
        dp[i][j] = 0
  print "length of LCS is :",max_len
  print "LCS :",lcs_str
bottom_up_dp_lcs(str_a, str_b)

3.2 優化空間複雜度的動態規劃算法

def space_efficient_lcs(str_a, str_b):
  """
  longest common substring of str_a and str_b, with O(n) space complexity
  """
  if len(str_a) == 0 or len(str_b) == 0:
    return 0
  max_len = 0
  dp = [0 for _ in range(len(str_b) + 1)]
  for i in range(1, len(str_a) + 1):
    left_up = 0
    for j in range(1, len(str_b) + 1):
      up = dp[j]
      if str_a[i-1] == str_b[j-1]:
        dp[j] = left_up + 1
        max_len = max([max_len, dp[j]])
      else:
        dp[j] = 0
      left_up = up
  print max_len
space_efficient_lcs(str_a, str_b)