首页 > 代码库 > Min Edit Distance

Min Edit Distance

2024-07-15 11:07:37   220人阅读

Min Edit Distance

————两字符串之间的最小距离

PPT原稿参见Stanford；http://www.stanford.edu/class/cs124/lec/med.pdf

Tips：由于本人水平有限，对MED的个人理解可能有纰漏之处，请勿尽信。

Edit：个人理解指编辑之意，也即对于两个字符串，对其中的一个进行各种编辑操作（插入、删除、替换）使其变为另一个字符串。要解决的问题是求出最小的编辑操作次数是多少。

clip_image002[4]

基因系列比对

定义距离：

X，Y是大小分别为n，m的字符串。

定义D(i,j)表示X[1..i],Y[1..j]两子字符串的距离。则：

X与Y的距离为D(n,m)

应用动态规划的方法：

clip_image012[4]

对于D(i,j)的计算结果做成表格（矩阵），D(i,j)的运算结果可以建立在之前的结果之上。

对本算法的一点个人理解：

设 i：输出子字符串长度

    j：输入子字符串长度

D(0,0)=0

D(i,0)：输入0个字符转换到i个字符的输出，也即要插入i个字符，代价为insertCost*i

D(0,j)：目标长度为0，输入长度为j，所以代价为：deletionCost*j

D(i,j)的上一步可以分为三种情况：

1、上一步输入长度为j，输出长度为i-1，那么现在这一步肯定至少还要插入一个字符才能达到i长度输出：

D(i-1,j)+inseartCost*1

2、上一步输入长度为j-1,输出长度为i-1，那么现在这一步第j个输入只需要做替换处理（如果第j个输入与第i个输出不相等）或者保持不变（如果第j个输入与第i个输出相等）：

D(i-1,j-1)+substituteCost*(source[j]==target[i] ? 0 : 1)

3、上一步输入长度为j-1,输出长度为i，那么现在这一步由于又多了一个字符，所以要把多的这个字符删除：

D(i,j-1)+deletionCost*1

由于我们要求的是最小Edit Distance，自然，就是上述三种情况中最小值做为D(i,j)的值。具体算法如下：

最小编辑距离动态算法（Levenshtein）：

clip_image014[4]

D(n,m)即为最小距离。

 

字符串对齐：

    每一次计算D(i,j)时记录当前数据是由哪个数据得的。当D(n,m)的计算完成后，即可从D(n,m)出发进行回溯(backtrace).得到和条从D(n,m)到D(0,0)的路径:

clip_image019[4]

clip_image021[4]

上图中从(0,0)到(n,m)的每一个非递减路径即为两字符串的对齐关系。

附加回溯过程的MED：

clip_image023[4]

Weighted Edit Distance:（加权编辑距离）

为什么要加权？

因为有些字符被写错的概率要大些（如搜索引擎中经常能自动搜索到相近的词句）

clip_image025[4]

算法：

clip_image027[4]

 

Levenshtein算法python实现：

#===============================================================================# Using dynamic programming to realize the MED algorithm(Levenshtein)# MED: short for Minimum Edit Distance#===============================================================================import typesimport numpy as npclass MED:    def __init__(self):        self.insCost = 1    #insertion cost        self.delCost = 1    #deletion cost        self.subsCost = 1   #substitution cost        self.D = 0    def mDistance(self, word1, word2):        assert type(word1) is types.StringType        assert type(word2) is types.StringType        size1 = len(word1)        size2 = len(word2)        if size1==0 or size2 ==0:            return size1*self.delCost+size2*self.insCost        D_mat = np.zeros((size1+1,size2+1))        D_mat[:,0] = range(size1+1)        D_mat[0,:] = range(size2+1)        for i in range(1,size1+1):            for j in range(1,size2+1):                insert_cost = D_mat[i-1, j] + self.insCost*1                delete_cost = D_mat[i, j-1] + self.delCost*1                if word1[i-1]==word2[j-1]:                    substitue_cost = D_mat[i-1, j-1]                else:                    substitue_cost = D_mat[i-1, j-1] + self.subsCost*1                D_mat[i,j] = min(insert_cost, delete_cost, substitue_cost)        self.D = D_mat        return D_mat[size1, size2]if __name__ == "__main__":    word1 = "Function"    word2 = "fanctional"    med = MED()    print "MED distance is :" ,med.mDistance(word1, word2)

输出结果：

MED distance is : 4.0

 

扩展应用1：利用回溯做字符对齐

对原ED计算函数做点更改（每次得到MED时记录该值的来源,然后从D(n,m)开始利用记录的来源回溯至起始位置，得到该MED的完整路径）：

def computingAlignment(self, word1,  word2):        assert type(word1) is types.StringType        assert type(word2) is types.StringType        size1 = len(word1)        size2 = len(word2)        if size1==0 or size2 ==0:            return size1*self.delCost+size2*self.insCost        D_mat = np.zeros((size1+1,size2+1))        D_rec = np.zeros((size1+1, size2+1))        D_mat[:,0] = range(size1+1)        D_mat[0,:] = range(size2+1)        for i in range(1,size1+1):            for j in range(1,size2+1):                insert_cost = D_mat[i-1, j] + self.insCost*1                delete_cost = D_mat[i, j-1] + self.delCost*1                if word1[i-1]==word2[j-1]:                    substitue_cost = D_mat[i-1, j-1]                else:                    substitue_cost = D_mat[i-1, j-1] + self.subsCost*1                D_mat[i,j] = min(insert_cost, delete_cost, substitue_cost)                if D_mat[i,j] == insert_cost:#Record Where the min val comes from                    D_rec[i,j] = 1                elif D_mat[i,j]== delete_cost:                    D_rec[i,j] = 2                else:                    D_rec[i,j] = 3        self.D = D_mat        self.D_rec = D_rec        # BackTrace        alignRevPath=[]#Get the reverse path        j = size2        i = size1        while i!=0 or j!=0:#Be carefull of this row            alignRevPath.append([i,j,D_rec[i,j]])            if D_rec[i,j]==1:                i -=1            elif D_rec[i,j]==2:                j -=1            elif D_rec[i,j]==3:                i -=1                j-=1            elif D_rec[i,j]==0:                if i>0:                    i -= 1                if j>0:                    j -= 1        alignStr1 =[]        alignStr2 =[]        if alignRevPath[-1][0]!=0:#process the first postion of the path            alignStr1.append(word1[alignRevPath[-1][0]-1])        else:            alignStr1.append(‘*‘)        if alignRevPath[-1][1]!=0:            alignStr2.append(word2[alignRevPath[-1][1]-1])        else:            alignStr2.append(‘*‘)                for i in range(len(alignRevPath)-1, 0, -1): #process the rest of the path            k = np.subtract(alignRevPath[i-1], alignRevPath[i])            bK = k>0            if bK[0]!=0:                alignStr1.append(word1[alignRevPath[i-1][0]-1])            else:                alignStr1.append(‘*‘)                            if bK[1]!=0:                alignStr2.append(word2[alignRevPath[i-1][1]-1])            else:                alignStr2.append(‘*‘)        return (alignStr1, alignStr2)

上面的代码中alignRevPath用来保存路径上的每一个位置，其每个元素都为3元列表，前两维为路上的坐标，第3维取0、1、2、3四种值，0表示路径到达边界了，1表示当前的ED结果由前一个insert操作得到，2表示当前ED结果由前一个Delete得到，3表示当前ED结果由前一个substitue得到。

在主程序中加入如下测试代码：

word1 = "efnction"word2 = "faunctional"med = MED()w1, w2 = med.computingAlignment(word1, word2)print ‘ ‘.join(w1)print ‘| ‘*len(w1)print ‘ ‘.join(w2)

得到的输出：

clip_image002[6]

 

扩展应用2：匹配最长子字符串(不一定连续）

原理：

clip_image002[8]

clip_image004

图示：

clip_image006

从上图可以看出，匹配的并不一定是连续的子串.这是因为我们的惩罚项设置为：

clip_image008

也即迭代过程clip_image010中的s(xi,yj)取-1（不匹配）或1（匹配）。当子串中有不匹配的字符出现时，将会对之前的匹配计数减1。

如果想匹配最长连续子串，可以令惩罚项F(i,j)为0（不匹配）或F(i-1,j-1)+1（匹配）

Python 实现（对computingAlignment做少许更改）：

def longestSubstr(self, word1,  word2):        assert type(word1) is types.StringType        assert type(word2) is types.StringType        size1 = len(word1)        size2 = len(word2)        if size1==0 or size2 ==0:            return size1*self.delCost+size2*self.insCost        D_mat = np.zeros((size1+1,size2+1))        D_rec = np.zeros((size1+1, size2+1))        D_mat[:,0] = 0        D_mat[0,:] = 0        for i in range(1,size1+1):            for j in range(1,size2+1):                insert_cost = D_mat[i-1, j] - self.insCost*1                delete_cost = D_mat[i, j-1] - self.delCost*1                if word1[i-1]==word2[j-1]:                    substitue_cost = D_mat[i-1, j-1] +1                else:                    substitue_cost = D_mat[i-1, j-1] - self.subsCost*1#                     substitue_cost = 0                D_mat[i,j] = max(0,insert_cost, delete_cost, substitue_cost)                if D_mat[i,j] == insert_cost:#Record Where the min val comes from                    D_rec[i,j] = 1                elif D_mat[i,j]== delete_cost:                    D_rec[i,j] = 2                elif D_mat[i,j]==substitue_cost:                    D_rec[i,j] = 3                            self.D = D_mat        self.D_rec = D_rec        maxVal = np.max(D_mat)        maxBool = D_mat == maxVal        numMax = np.sum(maxBool)        alignRevPaths=[]        for i in range(numMax):            alignRevPaths.append([])        pathStarts=[]        for i in range(size1+1):            for j in range(size2+1):                if maxBool[i,j]:                    pathStarts.append([i,j])                for m in range(numMax):            i = pathStarts[m][0]            j = pathStarts[m][1]            while i!=0 and j!=0:                alignRevPaths[m].append([i,j,D_rec[i,j]])                if D_rec[i,j]==1:                    i -=1                elif D_rec[i,j]==2:                    j -=1                elif D_rec[i,j]==3:                    i -=1                    j-=1                elif D_rec[i,j]==0:                    break                if D_mat[i,j]==0:                    break        str1=[]        str2=[]        for m in range(numMax):            str1.append([])            str2.append([])            p = alignRevPaths[m]            for i in range(len(p)-1, -1,-1):                str1[m].append(word1[p[i][0]-1])                str2[m].append(word2[p[i][1]-1])        return ([‘‘.join(i) for i in str1],[‘‘.join(i) for i in str2])

代码测试：

word1 = "ATCAT"word2 = "ATTATC"med = MED()str1,str2=  med.longestSubstr(word1, word2)print "Longest match substr:"print "match substr in word1:", str1print "match substr in word2:", str2

输出结果：

Longest match substr:

match substr in word1: [‘ATC‘, ‘ATCAT‘]

match substr in word2: [‘ATC‘, ‘ATTAT‘]

如果把longestSubstr中的

substitue_cost = D_mat[i-1, j-1] - self.subsCost*1

改为：

substitue_cost = 0

即可用来求最大连续子字符串，同样运行上述测试代码，可得：

Longest match substr:

match substr in word1: [‘ATC‘]

match substr in word2: [‘ATC‘]

改用另外一组字符串进行实验进行进一步验证：

word1 = "fefnction"word2 = "faunctional"med = MED()str1,str2=  med.longestSubstr(word1, word2)print "Longest match substr:"print "match substr in word1:", str1print "match substr in word2:", str2

当求解的是最长子字符串（非连续）时，输出：

Longest match substr:

match substr in word1: [‘nction‘]

match substr in word2: [‘nction‘]

当求解的是最长连续子字符串时，输出：

Longest match substr:

match substr in word1: [‘nction‘]

match substr in word2: [‘nction‘]

 

在以上的算法中，MED的一个最大的特点就是利用了矩阵保存之前处理的结果数据，以做为下一次的输入。对于最长子字符串的查找，同样套用了MED的框架，但仔细一想我们会发现，其实最长子串的查找并不一定需要记录路径alignRevPaths，有了alignRevPaths这个路径，我们编程时是根据这个路径处理字符的。路径的设置主要是为了在computingAlignment中对字符进行对齐，因为字符对齐的情况下，字符不一定是连续，比如会有如下对齐的形式中的“＊”号：

clip_image002[10]

所以我们才想到要用路径做记录。

但在最长子串中，子串肯定是连续的，自然路径也就不需要了。

下面对最长子串程序做简化：

def longestSubstr2(self,word1, word2):        assert type(word1) is types.StringType        assert type(word2) is types.StringType        size1 = len(word1)        size2 = len(word2)        if size1==0 or size2 ==0:            return size1*self.delCost+size2*self.insCost        D_mat = np.zeros((size1+1,size2+1))        D_mat[:,0] = 0        D_mat[0,:] = 0        for i in range(1,size1+1):            for j in range(1,size2+1):                insert_cost = D_mat[i-1, j] - self.insCost*1                delete_cost = D_mat[i, j-1] - self.delCost*1                if word1[i-1]==word2[j-1]:                    substitue_cost = D_mat[i-1, j-1] +1                else:                    substitue_cost = D_mat[i-1, j-1] - self.subsCost*1#                     substitue_cost = 0                D_mat[i,j] = max(0,insert_cost, delete_cost, substitue_cost)        self.D = D_mat        maxVal = np.max(D_mat)        maxBool = D_mat == maxVal        numMax = np.sum(maxBool)        alignRevPaths=[]        for i in range(numMax):            alignRevPaths.append([])        pathStarts=[]        for i in range(size1+1):            for j in range(size2+1):                if maxBool[i,j]:                    pathStarts.append([i,j])                str1=[]        str2=[]        for m in range(numMax):            str1.append([])            str2.append([])            i = pathStarts[m][0]            j = pathStarts[m][1]            s1Tmp = []            s2Tmp = []            while D_mat[i,j]!=0:                s1Tmp.append(word1[i-1])                s2Tmp.append(word2[j-1])                i -= 1                j -= 1            str1[m]=[s1Tmp[len(s1Tmp)-i-1] for i in range(len(s1Tmp))]            str2[m]=[s2Tmp[len(s2Tmp)-i-1] for i in range(len(s2Tmp))]        return ([‘‘.join(i) for i in str1],[‘‘.join(i) for i in str2])

使用以下字符做实验：

word1 = "ATCAT"

word2 = "ATTATC"

输出结果同样为：

Longest match substr:

match substr in word1: [‘ATC‘, ‘ATCAT‘]

match substr in word2: [‘ATC‘, ‘ATTAT‘]

一点个人感悟：

虽然这里的动态规划算法看上去有点不可思议，但联想起线性代数中的矩阵运算也就不难理解了。

就拿最长子串的程序来说，其实际过程仍可看做：对每word1中的每一个字符在word2中进行查找，当匹配第一个字符后继续匹配第二个字符，然后第三、第四……个，直到有字符不匹配时，记录该匹配成功的串长度及字串始末位置。

而这里的Smith-waterman算法将这个匹配记录在一个2维矩阵（数组）中。我们都知道，线性代数中的矩阵乘法同样是可以展开为各元素的乘与累加操作的。而这里，我认为，发明者正是利用了这一思想。


联系
我们