分类： C/C++
2012-01-27 16:39:58
动态规划(1)：最长公共子序列
[注：本文原来是写在另外一个网志上的一个技术专题，为了方便查找，都挪过来。其中全部程序都打包到dynamicProg.tar文件中，可以到“资源中心”下载。]
生物学应用中我们经常需要比较两个或更多个DNA序列。有机体的DNA序列是由AGCT四种碱基排列构成。对了，原文提到了AGCT所代表的原词： adenine，腺嘌呤脱氧核苷酸，guanine，鸟嘌呤脱氧核苷酸，cytosine，胞核嘧啶脱氧核苷酸，thymine，胸腺嘧啶脱氧核苷酸。 靠，世界上怎么有这么生僻的词，如果有个土生土长的中国人谈起DNA的时候脱口而出A strand of deoxyribonucleic acid consists of a string of the bases adenine, guanine, cytosine and thymine：太帅了。扯远了，还说正题：例如两个可能的DNA序列是
S1= ACCGGTCGAGTGCGCGGAAGCCGGCCGAA
S2 = GTCGTTCGGAATGCCGTTGCTCTGTAAA
比较两个DNA的目标之一，就是要知道它们到底有"多么"的相似。相似度的衡量可以有很多标准，例如我们可以说，如果其中一个DNA是另外一个DNA的子 序列，他们就是相似的。上面的例子种S1, S2都不是另外一个的子序列。或者我们也可以这样定义相似度：如果只用很少的修正操作（比如替换或者插入、删除）就能让一个DNA变成另一个，就说它们是 相似的。还可以这样定义：对两个DNA，S1和S2，找到一个碱基序列S3，S3中出现的全部碱基都既出现在S1中又出现在S2中，而且出现的顺序相同， 但是可以不连续：在这个前提下能找到的最长的S3越长，就说S1和S2越相似。在我们上面的例子中，能找到的最长的S3为 GTCGTCGGAAGCCGGCCGAA：
S1 ACCGGTCGAGTGCGCGGAAGCCGGCCGAA
S3 ____GTCG__T_CG__GAAGCCGGCCGAA
S2 GTCGTTCGGAATGCCGTTGCTCTGTAAA
S3 GTCGT_CGGA__GCCG__GC_C_G_AA_
上面对相似性的最后一种定义，可以形式化为传说中的"最长公共子序列"问题。某个序列的子序列定义为原序列中的0个或多个元素被去掉之后剩下的元素序列。形式化定义为：给定一个序列
X=x1, x2, x3, ..., xm
和另外一个序列
Z=z1, z2, z3, ..., zk
Z是X的子序列当且仅当存在一个严格递增序列i1, i2, …, ik，其中元素都是1-m之间的整数（X的元素的下标），而且对全部j=1, 2, …, k都有x[ij]=zj。例如Z=BCDB是X=ABCBDAB的子序列，对应的序列I为：i1=2, i2=3, i3=5, i4=7。
给定序列X和Y，如果Z既是X的子序列又是Y的子序列，就说序列Z是X, Y的公共子序列。例如若X=ABCBDAB，Y=BDCABA，那么它们有公共子序列BCA。这个例子中BCA不是XY的最长公共子序列，因为BCBA也 是XY的公共子序列而且比BCA更长。BDAB也是上述XY的公共子序列。事实上BABA和BDAB都是XY的最长公共子序列，因为XY没有长度为5或更 长的公共子序列。
在最长公共子序列问题中，我们拿到的是两个序列X=x1, x2, ..., xm和Y=y1, y2, ..., yn，并且希望找到XY的最长公共子序列。下面就是这个问题的动态规划解法。
先定义"前缀（prefix）"：给定一个序列X=x1, x2, ..., xm，对i=0, 1, ..., m，Xi=x1, x2, ..., xi都是X的前缀。
定理：对序列X=x1, x2, ..., xm和Y=y1, y2, ..., yn，Z=z1, z2, ..., zk是XY的最长公共子序列
1- 如果xm=yn，那么zk=xm=yn而且Z[k-1]是X[m-1]和Y[n-1]的最长公共子序列。
2- 如果xm!=yn
a. zk!=xm，则Z是X[m-1]和Y的最长公共子序列。
b. zk!=yn，则Z是X和Y[n-1]的最长公共子序列。
证明：
1- 如果zk!=xm，可以在序列Z的最后加上xm=yn，得到一个长度为k+1的公共子序列，这与假设Z是最长公共子序列矛盾。因此一定有zk =xm=yn。这样Z[k-1]是X[m-1]和Y[n-1]的长度为k-1的公共子序列，我们想要证明Z[k-1]也是最长公共子序列。假定X[m- 1]和Y[n-1]有公共子序列长度超过k-1，那么在X[m-1]和Y[n-1]后面分别添加xm=yn，即可得到XY的一个长度超过k的子序列，这与 假设Z是最长公共子序列矛盾。
2- 如果对X[m-1]和Y存在公共子序列W，长度大于k，那么W也会是Xm和Y的公共子序列，这和假设Z（长度k）是XY的最长公共子序列矛盾。2b同理可证。
上述定理说明，试图寻找任何X=x1, x2, ..., xm和Y=y1, y2, ..., yn的最长公共子序列的过程，只须解决至多两个子问题。如果xm=yn，就找X[m-1]和Y[n-1]的最长公共子序列，并在其后添上xm=yn就得到 了XY的最长公共子序列。如果xm!=yn，就分别求X[m-1]和Y的最长公共子序列，以及X和Y[n-1]的最长公共子序列，这两个哪一个比较长，哪 一个就是XY的最长公共子序列。
为了求X和Y[n-1]的，以及X[m-1]和Y的最长公共子序列，可能需要求出下列各组的最长公共子序列：
X[m-1] Y[n-2]
X[m-1] Y[n-1]
X Y[n-2]
X[m-2] Y[n-1]
X[m-1] Y[n-1]
X[m-2] Y
其中求X[m-1] Y[n-1]的最长公共子序列，是重复出现的子问题。算法进一步进行，还会有很多更细分的子问题重复出现，应该看到这是应用动态规划解题的一个好机会。定 义c[i,j]为序列Xi和Yj的最长公共子序列的长度，如果i=0或者j=0显然有c[i,j]=0，从上面的定理容易得到下面的递归公式：
c[i,j]=0 如果i=0或者j=0
c[i,j]=c[i-1,j-1]+1 如果xi=yj
c[i,j]=min(c[i,j-1], c[i-1,j]) 如果xi!=yj
在这个公式的基础上不难得到一个指数复杂度的递归算法求最长公共子序列，但是考虑到对长度分别为m, n的序列X, Y，其前缀分别有m+1, n+1个，不同的子问题个数的最多有(m+1)(n+1)个，可以用动态规划法自底向上求解。算法是：
对给定的Xm, Yn
1- 如果m=0或者n=0，返回0，算法结束。
2- 建立数组c[m+1][n+1]，将c[0][0..n], c[0..m][0]初始化为0。[注：c[i][j]在i=0或j=0时取值为0。]
3- 按照从左到右，从上到下的顺序填表。对c[i][j]来说，如果xi=yj，那么令delta=1, c[i][j]的值是以下三者中最大的一个：c[i-1][j], c[i][j-1], c[i-1][j-1]+delta。
4- 表c填满后，c[m][n]的值就是最大公共子序列的长度。
5- 填表过程结束后，按照第3步的规则倒推，即可知道哪些i, j的值是出现在最长公共子序列中的下标。
程序在清单lcs.c中给出，标准C。
动态规划(2)：编辑距离
先说明什么是编辑。给定两个字符串x[1..m]和y[1..n]，可以对他们进行很多编辑操作。目标是找到一系列的编辑操应用于x使它变为y。用一个中 间变量z来保存中间结果，算法开始的时候z是一个空字符串，算法结束的时候对所有j=1, 2, ..., n有zj=yj。用变量i标记当前正被处理的字符在x中的下标，用j标记z的当前下标，我们的编辑操作就是作用于i, j和z上面的。开始时有i=j=1，要求是在转换过程中一一检查x中的每一个字符，也就是说在转换结束的时候必须有i=m+1。
可供选用的编辑操作包括以下六种：
拷贝：把一个字符从x拷贝到z，即令zj=xi，之后令i=i+1, j=j+1。这个操作检查了xi。
替换：将x中的一个字符替换为另外一个，即令zj=c，之后令i=i+1, j=j+1。这个操作检查了xi。
删除：从x中删除一个字符，即令i=i+1，j不变。这个操作检查了xi。
插入：（呵呵呵，这是我最喜欢的操作！）向z中插入一个字符，即令zj=c，之后令j=j+1，i不变。这个操作没有检查x中的任何字符。
换位：把x中将要处理的下两个字符拷贝到x中，但是要颠倒顺序，即令zj=x[i+1], z[j+1]=xi之后令i=i+2, j=j+2。这个操作检查了xi和x[i+1]。
截断：删掉x中剩下的全部字符，即令i=m+1。这个操作检查了x中当前尚未被检查到的全部字符，但是这个操作只能作为最后一个操作出现在编辑操作序列中。
下面举例说明。设原字符串为'algorithm'，目标为'altruistic'。下面的编辑操作序列就可以把原字符串变换为目标字符串（注意这只是可能的编辑操作序列之一，事实上还有其他编辑序列能完成同样的功能，*标出的是i的位置）：
Operation # x # z
- # *algorithm # _
copy # a*lgorithm # a_
copy # al*gorithm # al_
replace by t # alg*orithm # alt_
delete # algo*rithm # alt_
copy # algor*ithm # altr_
insert u # algor*ithm # altru_
insert i # algor*ithm # altrui_
insert s # algor*ithm # altruis_
twiddle # algorit*hm # altruisti_
insert c # algorit*hm # altruistic_
kill # algorithm* # altruistic_
每一个变换操作都有相应的代价。假设(1)所有编辑操作的代价都是常数并且是已知的，(2)拷贝和替换操作的代价，小于等效的删除、插入操作的代价的和。定义给定的编辑操作序列的代价是其中每一个编辑操作的代价的总和，对上面的例子来说，编辑序列的代价是
(3*cost(copy)) +cost(replace) +cost(delete) +(4*cost(insert)) +cost(twiddle) +cost(kill).
a. 给定两个序列x[1..m]和y[1..n]以及各编辑操作的代价，x到y的编辑距离定义为：能把序列x变为y的编辑操作序列的最小代价。描 述一个用来求给定序列x[1..m]和y[1..n]的编辑距离的动态规划算法，打印最佳的编辑序列，分析该算法的时间和空间复杂度。
编辑距离问题，是DNA序列对齐问题的泛化。可以有很多方式对齐DNA序列，以衡量它们的相似程度。有一种对齐方式是，可以在两个DNA序列的任何位置 （包括两头）插入空格，但是要求通过这个过程得到的序列x1和y1长度相同，而且在两个序列的同一位置上不能都是空格。做到这一步之后我们为每一个位置分 配一个"得分"，例如位置j按照如下规则得分：(1)如果x1[j]==y1[j]，1分；(2)如果x1[i]!=y1[j]并且两者都不是空格，-1 分；(3)如果x1[j]或者y1[j]是空格，-2分。某种特定对齐结果的得分，是全部位置的得分的总和。例如对序列x=GATCGGCAT 和 y=CAATGTGAATC，一个可能的对齐如下：
G ATCG GCAT
CAAT GTGAATC
-*++*+*+-++*
其中+对应1分，-对应-1分，*对应-2分。这个对齐的总分是-1-2+2-2+1-2+1-1+2-2=-1-2+1-2=-4。
b. 解释如何从给定的编辑操作集选择一个子集，从而把寻找总分最高的对齐结果的问题，转化为求序列编辑距离的问题的一个特例。
用动态规划解题。
记给定的两个序列是X[1..m]和Y[1..n]。对给定序列X=x1, x2, ..., xm，以及i=0, 1, ..., m，定义序列Xi=x1, x2, ..., xi都为X的前缀。记c[i][j]为前缀Xi和Yj的编辑距离。由题意可以知道c[0][0]=0，对k=1, 2, .., n都有c[0][k]=k*cost(insert)，（即当X都没有被检查，而Y的第1, 2, .., n个元素已经被插入的时候，总开销就是k个插入操作的代价）对k2=1, 2, .., m都有c[k2][0]=k2*cost(delete)。（即当Y完全没有被构造，而Y的第1, 2, .., n个元素已经被删除的时候，总开销就是k个删除操作的代价。）
两个前缀Xi和Yj之间的编辑距离c[i][j]，是以下各值中最小的一个：
c[i-1][j] + cost(insert)
c[i][j-1] + cost(delete)
c[i-1][j-1] + cost(replace) [if xi!=yj]
c[i-1][j-1] + cost(copy) [if xi==yj]
c[i-2][j-2] + cost(twiddle) [ if i>1 && i
这个递推公式中没有kill操作。由于kill操作可以一次删除一个子串，我们对这个操作特别处理如下。首先按照上述过程，忽略kill操作的存在，填表求出只用另外五个操作能达到的编辑距离。
对所有i=1..m-1, j=1..n，当X被检查到xi，Y被生成到yj的时候，用一个kill操作可以删除X中剩下的元素，用(n-i)个insert操作可以插入Y中尚未出 现的全部元素。这样即可把X中的元素全部检查到，并且生成完整的Y。这种方式把X变成Y的总开销为c[i][j]+cost(kill)+(n-i) *cost(insert)。对所有c[i][j]求这个值，并且和上述忽略kill操作的前提下求出的编辑距离相比，取其中最小值，即求得编辑距离。
按照上面算法解题的程序清单：ed.c
复杂度。该算法除去kill操作不考虑，剩下的部分是对一个(m+1)(n+1)的表格进行填写的过程，填写每一个格c[i][j]花费的开销，只是根据 c[i-1][j], c[i][j-1], c[i-1][j-1], c[i-2][j-1]进行简单的加减运算，得到五个数字，选择其中最小的一个。所以这部分的时间复杂度为O(mn)。考虑kill操作的步骤，是对上述 表格的每个格，在常数的时间，内求出用一个kill操作和n个插入操作完成转换的开销，复杂度也是O(mn)。这两个部分相加，总的时间复杂度仍然是mn 的常数倍，O(mn)。
做过了最长公共子序列和编辑距离这两道题，解开了让我迷惑了很久的一个问题--UNIX的diff工具是用什么算法比较文件的。
动态规划这个系列已经写了两篇了，这个系列的第三篇--求数列的最长单调递增子序列--没有这两题难度大，不过以我并不灵光的头脑，还是想了许久才想出来的。
动态规划(3)：最长单调递增子序列
先回答编辑距离问题中的第二问。DNA对齐问题中，最高得分的对齐方法可以通过选择下面四个操作并相应赋给开销，来转换成一个编辑距离问题：
insert, -2
delete, -2
copy, +1
replace, -1
(1)如果x1[j]==y1[j]，1分；对应的是copy操作；
(2)如果x1[i]!=y1[j]并且两者都不是空格，-1分；对应的是replace操作；
(3)如果x1[j]或者y1[j]是空格，-2分，对应的是insert和delete操作。
最长单调子序列问题也是可以用动态规划解决的典型问题。问题的描述是对一个整数数列X=x1, x2, .., xn，求其子序列xi1, xi2, .., xim，要求1 <= i1 < i2 < .. < im <= n，且xi1, xi2, .., xim单调递增。
解决该问题的算法描述如下。把xi之前的最长单调递增子序列的长度记为ci，为了找到最长单调序列中的元素，另外设p1..pn记录，pi的值就是pi结尾的子序列中位于pi之前的元素，在原序列中的下标。如果xi作为一个子序列的第一个元素出现，令pi=0。
显然有c1=1，p1=0。对i=2..n，首先记下x1..x[i-1]中小于xi的全部元素的下标m1, m2, .., mp。之后取cm1, cm2, .., cmp中的最大值cx，令ci=cx+1，pi=x。如果x1..x[i-1]中没有比xi小的元素，令ci=1，pi=0。
上述过程结束后，检察c1..cn，其中最大值ck就是原数列的最长单调子数列的长度，确定其中元素的过程为，令i=k，count = 0:
result[count]=xi
如果pi=0结束
count++
i=pi
返回开头
分析复杂度。对全部i=1..n，需要找到小于xi的全部元素，及他们对应的ci的最大值。这可以通过对x1..xi的一次扫描完成，扫描一个元素所需的 时间(t0)是一个常数。寻找该子数列所含元素的操作，因为我们用了标记数组p，所需的时间就是检查ck个标记并且把相应的xi拷贝出来的时间，这是ck 的一次函数。因此总的运行时间是
S=t0*[1+2+...+(n-1)] + t1*ck <= t0*(n*n-n)/2 + t1*n
总的时间复杂度是O(n*n)。
程序清单： lmis.c
《算法导论》的习题中提到，最长单调子数列问题有时间复杂度O(nlgn)的算法，不过我比较笨，看了书上的提示，想了半天还是没想出来怎么能把时间复杂度减小到O(nlgn)，在这里记一笔，等以后看到了或者水平提高了再补吧。
