算法导论程序38--最长公共子序列（Python）

最长公共子序列问题：

S1=ACCGGTCGAGTGCGCGGAAGCCGGCCGAA

S2=GTCGTTCGGAATGCCGTTGCTCTGTAAA

相似度的衡量：

寻找第三个串S3,它的所有碱基也都出现在S1和S2中，且在三个串中出现的顺序都相同，但在S1和S2中不要求连续出现。可以找到的S3越长，就可以认为S1和S2的相似度越高。

S3=GTCGTCGGAAGCCGGCCGAA

一个给定序列的子序列，就是将给定的序列中零个或多个元素去掉之后得到的结果。

其形式化定义如下：给定一个序列X=<x1，x2,....,xm>，另一个序列Z=<z1,z2,...zk>满足如下条件时称为X的子序列。即存在严格递增的X的下标序列

给定两个序列X和Y，如果Z既是X的子序列，也是Y的子序列，我们称它是X和Y的公共子序列。

最长公共子序列问题（longest-common-subsequence problem）：

给定两个序列X=<x1,x2,....,xm>和Y=<y1,y2,....,yn>，求X和Y长度最长的公共子序列。本节将展示如何用动态规划方法高效地求解LCS问题。

与矩阵链乘法问题相似：设计LCS问题的递归算法首先要建立最优解的递归式。

我们定义C[i, j]表示Xi和Yj的LCS的长度。

计算LCS的长度：

过程LCS-LENGTH接受两个序列X=<x1,x2,...,xm>和Y=<y1,y2,...,yn>为输入，它将c[i, j]的值保存在表c[0..m 0..n]中，并按行主次序计算表项。过程还维护一个表b[1..m 1..n]帮助构造最优解。b[i, j]指向的表项对应计算c[i, j]时所选择的子问题最优解。过程返回表b和表c，c[m , n]保存了X和Y的LCS的长度。

def lcs_length(X,Y):
    m=len(X)
    n=len(Y)
    c=[[0 for j in range(n+1)]for i in range(m+1)]
    b=[[" " for j in range(n+1)]for i in range(m+1)]
    for i in range(1,m+1):
        for j in range(1,n+1):
            if X[i-1]==Y[j-1]:
                c[i][j]=c[i-1][j-1]+1
                b[i][j]="left up"
            elif c[i-1][j]>=c[i][j-1]:
                c[i][j]=c[i-1][j]
                b[i][j]="up"
            else:
                c[i][j]=c[i][j-1]
                b[i][j]="left"
    return c,b 

def print_lcs(b,X,i,j):
    if i==0 or j==0:
        return
    if b[i][j]=="left up":
        print_lcs(b,X,i-1,j-1)
        print(X[i-1],end='')
    elif b[i][j]=="up":
        print_lcs(b,X,i-1,j)
    else:
        print_lcs(b,X,i,j-1)

if __name__=="__main__":
    X=["A","B","C","B","D","A","B"]
    Y=["B","D","C","A","B","A"]
    c,b=lcs_length(X,Y)
    for i in range(len(X)+1):
        for j in range(len(Y)+1):
            print(c[i][j],' ',end='')
        print()
    for i in range(len(X)+1):
        for j in range(len(Y)+1):
            print(b[i][j],' ',end='')
        print()  
    print_lcs(b,X,len(X),len(Y))

运行：

0  0  0  0  0  0  0  
0  0  0  0  1  1  1  
0  1  1  1  1  2  2  
0  1  1  2  2  2  2  
0  1  1  2  2  3  3  
0  1  2  2  2  3  3  
0  1  2  2  3  3  4  
0  1  2  2  3  4  4  
                     
   up  up  up  left up  left  left up  
   left up  left  left  up  left up  left  
   up  up  left up  left  up  up  
   left up  up  up  up  left up  left  
   up  left up  up  up  up  up  
   up  up  up  left up  up  left up  
   left up  up  up  up  left up  up  
BCBA

过程的运行时间为O(m+n)，因为每次递归调用i和j至少有一个会减少1。