如何修正此Damerau-Levenshtein实现中的错误？
我又来了，带着一个比较长的问题。这次我尝试了几种基于Python的Damerau-Levenshtein编辑距离的实现，最后找到了下面提到的这个，叫做editdistance_reference()。这个方法似乎能给出正确的结果，而且实现效率也不错。
于是我决定把这个代码转换成Cython。在我的测试数据中，参考方法能处理11,000次比较（对于大约12个字母的单词对），而Cython版本每秒能处理超过200,000次比较。可惜的是，结果不对：当我打印出用于调试的变量thisrow时，我的版本无论输入什么数据，里面都是一堆1，而参考输出却完全不同。例如，测试'helo'和'world'的结果如下（ED表示我的函数，EDR是正确的参考）：
来自editdistance()的结果：
#ED A [0, 0, 0, 0, 0, 1]
#ED B [1, 0, 0, 0, 0, 1]
#ED B [1, 1, 0, 0, 0, 1]
#ED B [1, 1, 1, 0, 0, 1]
#ED B [1, 1, 1, 1, 0, 1]
#ED B [1, 1, 1, 1, 1, 1]
#ED A [0, 0, 0, 0, 0, 2]
#ED B [1, 0, 0, 0, 0, 2]
#ED B [1, 1, 0, 0, 0, 2]
#ED B [1, 1, 1, 0, 0, 2]
#ED B [1, 1, 1, 1, 0, 2]
#ED B [1, 1, 1, 1, 1, 2]
#ED A [0, 0, 0, 0, 0, 3]
#ED B [1, 0, 0, 0, 0, 3]
#ED B [1, 1, 0, 0, 0, 3]
#ED B [1, 1, 1, 0, 0, 3]
#ED B [1, 1, 1, 1, 0, 3]
#ED B [1, 1, 1, 1, 1, 3]
#ED A [0, 0, 0, 0, 0, 4]
#ED B [1, 0, 0, 0, 0, 4]
#ED B [1, 1, 0, 0, 0, 4]
#ED B [1, 1, 1, 0, 0, 4]
#ED B [1, 1, 1, 1, 0, 4]
#ED B [1, 1, 1, 1, 1, 4]
来自editdistance_reference()的结果：
#EDR A [0, 0, 0, 0, 0, 1]
#EDR B [1, 0, 0, 0, 0, 1]
#EDR B [1, 2, 0, 0, 0, 1]
#EDR B [1, 2, 3, 0, 0, 1]
#EDR B [1, 2, 3, 4, 0, 1]
#EDR B [1, 2, 3, 4, 5, 1]
#EDR A [0, 0, 0, 0, 0, 2]
#EDR B [2, 0, 0, 0, 0, 2]
#EDR B [2, 2, 0, 0, 0, 2]
#EDR B [2, 2, 3, 0, 0, 2]
#EDR B [2, 2, 3, 4, 0, 2]
#EDR B [2, 2, 3, 4, 5, 2]
#EDR A [0, 0, 0, 0, 0, 3]
#EDR B [3, 0, 0, 0, 0, 3]
#EDR B [3, 3, 0, 0, 0, 3]
#EDR B [3, 3, 3, 0, 0, 3]
#EDR B [3, 3, 3, 3, 0, 3]
#EDR B [3, 3, 3, 3, 4, 3]
#EDR A [0, 0, 0, 0, 0, 4]
#EDR B [4, 0, 0, 0, 0, 4]
#EDR B [4, 4, 0, 0, 0, 4]
#EDR B [4, 4, 4, 0, 0, 4]
#EDR B [4, 4, 4, 4, 0, 4]
#EDR B [4, 4, 4, 4, 4, 4]
我觉得我可能很笨，因为这个错误可能是非常明显的事情，但我就是找不到。
还有第二个问题：我用malloc为三个数组twoago、oneago和thisrow分配空间，然后它们以循环的方式互相交换。当我尝试free( twoago )等操作时，glibc会抱怨double free or corruption。我在网上查过这个问题；难道是指针交换的操作让glibc有点晕，所以它无法正确释放内存吗？
下面我先列出需要编译的setup.py（/path/to/python3.1 ./setup.py build_ext --inplace），然后是编辑距离的代码，这样感兴趣的人更容易复现。
还有一件事：这个是在Python3.1下运行的；有趣的是在*.pyx文件里，我们确实有裸的unicode字符串，但print仍然是一个语句，而不是一个函数。
我知道这里贴了很多代码，但问题是如果你删减太多，代码就无法运行。我相信除了editdistance()以外的所有方法都能正常工作，但如果你发现任何问题，请随时指出。
setup.py：
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
setup(
name = 'cython_dameraulevenshtein',
ext_modules = [
Extension( 'cython_dameraulevenshtein', [ 'cython_dameraulevenshtein.pyx', ] ), ],
cmdclass = {
'build_ext': build_ext }, )
cython_dameraulevenshtein.pyx（请滚动到最后查看有趣的部分）：
############################################################################################################
cdef extern from "stdlib.h":
ctypedef unsigned int size_t
void *malloc(size_t size)
void *realloc( void *ptr, size_t size )
void free(void *ptr)
#-----------------------------------------------------------------------------------------------------------
cdef inline unsigned int _minimum_of_two_uints( unsigned int a, unsigned int b ):
if a < b: return a
return b
#-----------------------------------------------------------------------------------------------------------
cdef inline unsigned int _minimum_of_three_uints( unsigned int a, unsigned int b, unsigned int c ):
if a < b:
if c < a:
return c
return a
if c < b:
return c
return b
#-----------------------------------------------------------------------------------------------------------
cdef inline int _warp( unsigned int limit, int value ):
return value if value >= 0 else limit + value
############################################################################################################
# ARRAYS THAT SAY SIZE ;-)
#-----------------------------------------------------------------------------------------------------------
cdef class Array_of_unsigned_int:
cdef unsigned int *data
cdef unsigned int length
#---------------------------------------------------------------------------------------------------------
def __cinit__( self, unsigned int length, fill_value = None ):
self.length = length
self.data = <unsigned int *>malloc( length * sizeof( unsigned int ) ) ###OBS### must check malloc doesn't return NULL pointer
if fill_value is not None:
self.fill( fill_value )
#---------------------------------------------------------------------------------------------------------
cdef fill( self, unsigned int value ):
cdef unsigned int idx
cdef unsigned int *d = self.data
for idx from 0 <= idx < self.length:
d[ idx ] = value
#---------------------------------------------------------------------------------------------------------
cdef resize( self, unsigned int length ):
self.data = <unsigned int *>realloc( self.data, length * sizeof( unsigned int ) ) ###OBS### must check realloc doesn't return NULL pointer
self.length = length
#---------------------------------------------------------------------------------------------------------
def free( self ):
"""Always remember the milk: Free up memory."""
free( self.data ) ###OBS### should free memory here
#---------------------------------------------------------------------------------------------------------
def as_list( self ):
"""Return the array as a Python list."""
R = []
cdef unsigned int idx
cdef unsigned int *d = self.data
for idx from 0 <= idx < self.length:
R.append( d[ idx ] )
return R
############################################################################################################
# CONVERTING UNICODE TO CHARACTER IDs (CIDs)
#---------------------------------------------------------------------------------------------------------
cdef unsigned int _UMX_surrogate_lower_bound = 0x10000
cdef unsigned int _UMX_surrogate_upper_bound = 0x10ffff
cdef unsigned int _UMX_surrogate_hi_lower_bound = 0xd800
cdef unsigned int _UMX_surrogate_hi_upper_bound = 0xdbff
cdef unsigned int _UMX_surrogate_lo_lower_bound = 0xdc00
cdef unsigned int _UMX_surrogate_lo_upper_bound = 0xdfff
cdef unsigned int _UMX_surrogate_foobar_factor = 0x400
#---------------------------------------------------------------------------------------------------------
cdef Array_of_unsigned_int _cids_from_text( text ):
"""Givn a ``text`` either as a Unicode string or as a ``bytes`` or ``bytearray``, return an instance of
``Array_of_unsigned_int`` that enumerates either the Unicode codepoints of each character or the value of
each byte. Surrogate pairs will be condensed into single values, so on narrow Python builds the length of
the array returned may be less than ``len( text )``."""
#.........................................................................................................
# Make sure ``text`` is either a Unicode string (``str``) or a ``bytes``-like thing:
is_bytes = isinstance( text, ( bytes, bytearray, ) )
assert is_bytes or isinstance( text, str ), '#121'
#.........................................................................................................
# Whether it is a ``str`` or a ``bytes``, we know the result can only have at most as many elements as
# there are characters in ``text``, so we can already reserve that much space (in the case of a Unicode
# text, there may be fewer CIDs if there happen to be surrogate characters):
cdef unsigned int length = <unsigned int>len( text )
cdef Array_of_unsigned_int R = Array_of_unsigned_int( length )
#.........................................................................................................
# If ``text`` is empty, we can return an empty array right away:
if length == 0: return R
#.........................................................................................................
# Otherwise, prepare to copy data:
cdef unsigned int idx = 0
#.........................................................................................................
# If ``text`` is a ``bytes``-like thing, use simplified processing; we just have to copy over all byte
# values and are done:
if is_bytes:
for idx from 0 <= idx < length:
R.data[ idx ] = <unsigned int>text[ idx ]
return R
#.........................................................................................................
cdef unsigned int cid = 0
cdef bool is_surrogate = False
cdef unsigned int hi = 0
cdef unsigned int lo = 0
cdef unsigned int chr_count = 0
#.........................................................................................................
# Iterate over all indexes in text:
for idx from 0 <= idx < length:
#.......................................................................................................
# If we met with a surrogate CID in the last cycle, then that was a high surrogate CID, and the
# corresponding low CID is on the current position. Having both, we can compute the intended CID
# and reset the flag:
if is_surrogate:
lo = <unsigned int>ord( text[ idx ] )
# IIRC, this formula was documented in Unicode 3:
cid = ( ( hi - _UMX_surrogate_hi_lower_bound ) * _UMX_surrogate_foobar_factor
+ ( lo - _UMX_surrogate_lo_lower_bound ) + _UMX_surrogate_lower_bound )
is_surrogate = False
#.......................................................................................................
else:
# Otherwise, we retrieve the CID from the current position:
cid = <unsigned int>ord( text[ idx ] )
#.....................................................................................................
if _UMX_surrogate_hi_lower_bound <= cid <= _UMX_surrogate_hi_upper_bound:
# If this CID is a high surrogate CID, set ``hi`` to this value and set a flag so we'll come back
# in the next cycle:
hi = cid
is_surrogate = True
continue
#.......................................................................................................
R.data[ chr_count ] = cid
chr_count += 1
#.........................................................................................................
# Surrogate CIDs take up two characters but end up as a single resultant CID, so the return value may
# have fewer elements than the naive string length indicated; in this case, we want to free some memory
# and correct array length data:
if chr_count != length:
R.resize( chr_count )
#.........................................................................................................
return R
#---------------------------------------------------------------------------------------------------------
def cids_from_text( text ):
cdef Array_of_unsigned_int c_R =_cids_from_text( text )
R = c_R.as_list()
c_R.free() ###OBS### should free memory here
return R
############################################################################################################
# SECOND-ORDER SIMILARITY
#-----------------------------------------------------------------------------------------------------------
cpdef float similarity( char *a, char *b ):
"""Given two byte strings ``a`` and ``b``, return their Damerau-Levenshtein similarity as a float between
0.0 and 1.1. Similarity is computed as ``1 - relative_editdistance( a, b )``, so a result of ``1.0``
indicates identity, while ``0.0`` indicates complete dissimilarity."""
return 1.0 - relative_editdistance( a, b )
#-----------------------------------------------------------------------------------------------------------
cpdef float relative_editdistance( char *a, char *b ):
"""Given two byte strings ``a`` and ``b``, return their relative Damerau-Levenshtein distance. The return
value is a float between 0.0 and 1.0; it is calculated as the absolute edit distance, divided by the
length of the longer string. Therefore, ``0.0`` indicates identity, while ``1.0`` indicates complete
dissimilarity."""
cdef int length = max( len( a ), len( b ) )
if length == 0: return 0.0
return editdistance( a, b ) / <float>length
############################################################################################################
# EDIT DISTANCE
#-----------------------------------------------------------------------------------------------------------
cpdef unsigned int editdistance( text_a, text_b ):
"""Given texts as Unicode strings or ``bytes`` / ``bytearray`` objects, return their absolute
Damerau-Levenshtein distance. Each deletion, insertion, substitution, and transposition is counted as one
difference, so the edit distance between ``abc`` and ``ab``, ``abcx``, ``abx``, ``acb``, respectively, is
``1``."""
#.........................................................................................................
# This should be fast in Python, as it can (and probably is) implemented by doing an identity check in
# the case of ``bytes`` and ``str`` objects:
if text_a == text_b: return 0
#.........................................................................................................
# Convert Unicode text to C array of unsigned integers:
cdef Array_of_unsigned_int a = _cids_from_text( text_a )
cdef Array_of_unsigned_int b = _cids_from_text( text_b )
R = c_editdistance( a, b )
#.........................................................................................................
# Always remember the milk:
a.free()
b.free()
#.........................................................................................................
return R
#-----------------------------------------------------------------------------------------------------------
cdef unsigned int c_editdistance( Array_of_unsigned_int cids_a, Array_of_unsigned_int cids_b ):
# Conceptually, this is based on a len(a) + 1 * len(b) + 1 matrix.
# However, only the current and two previous rows are needed at once,
# so we only store those.
#.........................................................................................................
# This shortcut is pretty useless if comparison is not very fast; therefore, it is done in the function
# that deals with the Python objects, q.v.
# if cids_a.equals( cids_b ): return 0
#.........................................................................................................
cdef unsigned int a_length = cids_a.length
cdef unsigned int b_length = cids_b.length
#.........................................................................................................
# Another shortcut: if one of the texts is empty, then the edit distance is trivially the length of the
# other text. This also works for two empty texts, but those have already been taken care of by the
# previous shortcut:
#.........................................................................................................
if a_length == 0: return b_length
if b_length == 0: return a_length
#.........................................................................................................
cdef unsigned int row_length = b_length + 1
cdef unsigned int row_length_1 = row_length - 1
cdef unsigned int row_bytecount = sizeof( unsigned int ) * row_length
cdef unsigned int *oneago = <unsigned int *>malloc( row_bytecount ) ###OBS### must check malloc doesn't return NULL pointer
cdef unsigned int *twoago = <unsigned int *>malloc( row_bytecount ) ###OBS### must check malloc doesn't return NULL pointer
cdef unsigned int *thisrow = <unsigned int *>malloc( row_bytecount ) ###OBS### must check malloc doesn't return NULL pointer
cdef unsigned int idx = 0
cdef unsigned int idx_a = 0
cdef unsigned int idx_b = 0
cdef int idx_a_1_text = 0
cdef int idx_b_1_row = 0
cdef int idx_b_2_row = 0
cdef int idx_b_1_text = 0
cdef unsigned int deletion_cost = 0
cdef unsigned int addition_cost = 0
cdef unsigned int substitution_cost = 0
#.........................................................................................................
# Equivalent of ``thisrow = list( range( 1, b_length + 1 ) ) + [ 0 ]``:
#print( '#305', cids_a.as_list(), cids_b.as_list(), a_length, b_length, row_length, row_length_1 )
for idx from 1 <= idx < row_length:
thisrow[ idx - 1 ] = idx
thisrow[ row_length - 1 ] = 0
#.........................................................................................................
for idx_a from 0 <= idx_a < a_length:
idx_a_1_text = _warp( a_length, idx_a - 1 )
twoago, oneago = oneago, thisrow
#.......................................................................................................
# Equivalent of ``thisrow = [ 0 ] * b_length + [ idx_a + 1 ]``:
for idx from 0 <= idx < row_length_1:
thisrow[ idx ] = 0
thisrow[ row_length - 1 ] = idx_a + 1
#.......................................................................................................
# some diagnostic output:
x = []
for idx from 0 <= idx < row_length: x.append( thisrow[ idx ] )
print
print '#ED A', x
#.......................................................................................................
for idx_b from 0 <= idx_b < b_length:
#.....................................................................................................
idx_b_1_row = _warp( row_length, idx_b - 1 )
idx_b_1_text = _warp( b_length, idx_b - 1 )
#.....................................................................................................
assert 0 <= idx_b_1_row < row_length, ( '#323', idx_b_1_row, )
assert 0 <= idx_a_1_text < a_length, ( '#324', idx_a_1_text, )
assert 0 <= idx_b_1_text < b_length, ( '#325', idx_b_1_text, )
#.....................................................................................................
deletion_cost = oneago[ idx_b ] + 1
addition_cost = thisrow[ idx_b_1_row ] + 1
substitution_cost = oneago[ idx_b_1_row ] + ( 1 if cids_a.data[ idx_a ]
!= cids_b.data[ idx_b ] else 0 )
thisrow[ idx_b ] = _minimum_of_three_uints( deletion_cost, addition_cost, substitution_cost )
#.....................................................................................................
# Transpositions:
if ( idx_a > 0
and idx_b > 0
and cids_a.data[ idx_a ] == cids_b.data[ idx_b_1_text ]
and cids_a.data[ idx_a_1_text ] == cids_b.data[ idx_b ]
and cids_a.data[ idx_a ] != cids_b.data[ idx_b ] ):
#...................................................................................................
idx_b_2_row = _warp( row_length, idx_b - 2 )
assert 0 <= idx_b_2_row < row_length, ( '#340', idx_b_2_row, )
thisrow[ idx_b ] = _minimum_of_two_uints( thisrow[ idx_b ], twoago[ idx_b_2_row ] + 1 )
#.....................................................................................................
# some diagnostic output:
x = []
for idx from 0 <= idx < row_length: x.append( thisrow[ idx ] )
print '#ED B', x
#.........................................................................................................
# Here, ``b_length - 1`` can't become negative, since we already tested for ``b_length == 0`` in the
# shortcut above:
cdef unsigned int R = thisrow[ b_length - 1 ]
#.........................................................................................................
# Always remember the milk:
# BUG: Activating below lines leads to glibc failing with ``double free or corruption``
#free( twoago )
#free( oneago )
#free( thisrow )e
#.........................................................................................................
return R
#-----------------------------------------------------------------------------------------------------------
def editdistance_reference( text_a, text_b ):
"""This method is believed to compute a correct Damerau-Levenshtein edit distance, with deletions,
insertions, substitutions, and transpositions. Do not touch it; it is here to validate results returned
from the above method. Code adapted from
http://mwh.geek.nz/2009/04/26/python-damerau-levenshtein-distance"""
# Conceptually, the implementation is based on a ``( len( seq1 ) + 1 ) * ( len( seq2 ) + 1 )`` matrix.
# However, only the current and two previous rows are needed at once, so we only store those. Python
# lists wrap around for negative indices, so we put the leftmost column at the *end* of the list. This
# matches with the zero-indexed strings and saves extra calculation.
b_length = len( text_b )
oneago = None
thisrow = list( range( 1, b_length + 1 ) ) + [ 0 ]
for idx_a in range( len( text_a ) ):
twoago, oneago, thisrow = oneago, thisrow, [ 0 ] * b_length + [ idx_a + 1 ]
#.......................................................................................................
# some diagnostic output:
print
print '#EDR A', thisrow
#.......................................................................................................
for idx_b in range( b_length ):
deletion_cost = oneago[ idx_b ] + 1
addition_cost = thisrow[ idx_b - 1 ] + 1
substitution_cost = oneago[ idx_b - 1 ] + ( text_a[ idx_a ] != text_b[ idx_b ] )
thisrow[ idx_b ] = min( deletion_cost, addition_cost, substitution_cost )
if ( idx_a > 0
and idx_b > 0
and text_a[ idx_a ] == text_b[ idx_b - 1 ]
and text_a[ idx_a - 1 ] == text_b[ idx_b ]
and text_a[ idx_a ] != text_b[ idx_b ] ):
thisrow[ idx_b ] = min( thisrow[ idx_b ], twoago[ idx_b - 2 ] + 1 )
#.....................................................................................................
# some diagnostic output:
print '#EDR B', thisrow
#.....................................................................................................
return thisrow[ len( text_b ) - 1 ]
编辑 我还把这段文字发到了pastebin和Cython邮件列表。
1 个回答
2
先做一些基础的调试。你知道在第二行输出中出错了，标记为 #ED B。错误的值似乎表明它在一开始找到了一个编辑，但之后再也没有找到更多。这可能是因为 min() 函数的某个参数被限制在了1。打印出 deletion_cost、substitution_cost、addition_cost……哪个是错的？为什么会错？打印一下输入的文本值。暂时禁用掉转置部分，看看这样问题是否消失。检查并重新检查 _warp 的用法（这真是个狡猾的小把戏）。如果你比较 "aaaaa" 和 "aaaaa" 会发生什么？"qwerty" 和 "qwerty" 呢？"xxxxx" 和 "yyyyy" 呢？这个问题在 bytes、bytearray 和 str 输入中都会出现吗？
关于内存管理的一些旁注：你可能想看看 这个，考虑使用 Python 特有的例程，而不是 malloc/free。如果有替代品的话，缩小你的数组似乎有点过了。
更新：我按照自己的建议进行了操作。删除成本出了问题。"oneago" 和 "thisrow" 是一样的。这个问题导致了错误的答案和重复的（-！没有损坏！-）释放：指针的循环洗牌并没有真正形成循环。
# twoago, oneago = oneago, thisrow ### BUG ###
twoago, oneago, thisrow = oneago, thisrow, twoago ### FIXED ###
更新 2： [评论容量太小] 只是普通的调试工作，没有什么特别的。 "专注于这个修复" 并不是 "超级易读"。参考代码确实为每次循环创建了一个新列表，这是可以的，因为 thisrow 不依赖于上一次循环的任何内容。它并不需要这样做，实际上，除了第一个和最后一个元素之外的初始化可以是随机数，这些数只是用来填充列表，以便可以索引，而不是像某些简单实现那样追加。因此，你可以完全模仿 "参考实现"，代价是多做一次（浪费的）malloc/free，或者你可以忽略 Python 特有的实现细节，仅将参考实现作为正确答案的来源。然后你可以接受我的修复，之后再通过减少 thisrow 数组的大部分初始化来节省时间。
更新 3： 这是一个替代的参考实现给你。它最初分配了3行，以避免在外部循环中创建列表的开销。它还避免了对 thisrow 除最后一个元素外的所有元素进行不必要的初始化。这使得转换为 C/Cython 更加容易。
def damlevref2(seq1, seq2):
# For Python 2.x as was the original.
# Appears to work on Python 1.5.2 as well :-)
seq2len = len(seq2)
twoago = [-777] * (seq2len + 1) # pseudo-malloc; any old rubbish will do
oneago = [-666] * (seq2len + 1) # ditto
thisrow = range(1, seq2len + 1) + [0]
for x in xrange(len(seq1)):
twoago, oneago, thisrow = oneago, thisrow, twoago # circular "pointer" shuffle
thisrow[-1] = x + 1
for y in xrange(seq2len):
delcost = oneago[y] + 1
addcost = thisrow[y - 1] + 1
subcost = oneago[y - 1] + (seq1[x] != seq2[y])
thisrow[y] = min(delcost, addcost, subcost)
if (x > 0 and y > 0 and seq1[x] == seq2[y - 1]
and seq1[x-1] == seq2[y] and seq1[x] != seq2[y]):
thisrow[y] = min(thisrow[y], twoago[y - 2] + 1)
return thisrow[seq2len - 1]