这是一种显然是二次时间的方法,但常数因子相对较低,因为除了长度为 1 的对象之外,它不构建任何子字符串对象。结果是一个 2 元组,
bestlen, list_of_results
其中bestlen
是重复相邻块的最长子串的长度,每个结果是一个3元组,
start_index, width, numreps
意味着被重复的子字符串是
the_string[start_index : start_index + width]
还有numreps
那些相邻的。永远都是这样
bestlen == width * numreps
问题描述留下了歧义。例如,考虑这个输出:
>>> crunch2("aaaaaabababa")
(6, [(0, 1, 6), (0, 2, 3), (5, 2, 3), (6, 2, 3), (0, 3, 2)])
因此,它找到了 5 种方法来将“最长”的拉伸视为长度为 6:
- 最初的“a”重复了 6 次。
- 最初的“aa”重复了 3 次。
- “ab”最左边的实例重复了 3 次。
- 最左边的“ba”重复了 3 次。
- 最初的“aaa”重复了 2 次。
它不会返回 intro 或 outro,因为从它返回的内容中推断这些是微不足道的:
- 简介是
the_string[: start_index]
。
- 结局是
the_string[start_index + bestlen :]
。
如果没有重复的相邻块,则返回
(0, [])
其他示例(来自您的帖子):
>>> crunch2("EEEFGAFFGAFFGAFCD")
(12, [(3, 4, 3)])
>>> crunch2("ACCCCCCCCCA")
(9, [(1, 1, 9), (1, 3, 3)])
>>> crunch2("ABCD")
(0, [])
它如何工作的关键:假设您W
每个都有相邻的重复宽度块。然后考虑将原始字符串与左移字符串进行比较时会发生什么W
:
... block1 block2 ... blockN-1 blockN ...
... block2 block3 ... blockN ... ...
然后你(N-1)*W
在相同的位置得到连续的相等字符。但这也适用于另一个方向:如果您向左移动W
并找到(N-1)*W
连续相等的字符,那么您可以推断:
block1 == block2
block2 == block3
...
blockN-1 == blockN
所以所有N
块都必须是块1的重复。
因此,代码反复将原始字符串向左移动(复制)一个字符,然后从左向右移动,识别最长的相等字符。这只需要一次比较一对字符。为了使“左移”高效(恒定时间),字符串的副本存储在collections.deque
.
编辑:update()
在许多情况下做了太多徒劳的工作;替换它。
def crunch2(s):
from collections import deque
# There are zcount equal characters starting
# at index starti.
def update(starti, zcount):
nonlocal bestlen
while zcount >= width:
numreps = 1 + zcount // width
count = width * numreps
if count >= bestlen:
if count > bestlen:
results.clear()
results.append((starti, width, numreps))
bestlen = count
else:
break
zcount -= 1
starti += 1
bestlen, results = 0, []
t = deque(s)
for width in range(1, len(s) // 2 + 1):
t.popleft()
zcount = 0
for i, (a, b) in enumerate(zip(s, t)):
if a == b:
if not zcount: # new run starts here
starti = i
zcount += 1
# else a != b, so equal run (if any) ended
elif zcount:
update(starti, zcount)
zcount = 0
if zcount:
update(starti, zcount)
return bestlen, results
使用正则表达式
[由于大小限制删除了这个]
使用后缀数组
这是迄今为止我发现的最快的,尽管仍然可以引发二次时间行为。
请注意,是否找到重叠字符串并不重要。正如crunch2()
上面程序所解释的(这里以次要方式详细阐述):
s
给定长度为的字符串n = len(s)
。
- 给定整数
i
和j
with 0 <= i < j < n
。
然后w = j-i
,如果 和是和c
之间共同的前导字符数,则(长度为)的块重复,从 开始,总共重复次数。s[i:]
s[j:]
s[i:j]
w
s[i]
1 + c // w
下面的程序直接查找所有重复的相邻块,并记住最大总长度的块。返回与 相同的结果crunch2()
,但有时顺序不同。
后缀数组简化了搜索,但几乎没有消除它。后缀数组直接查找<i, j>
具有最大值的对c
,但仅将搜索限制为最大化w * (1 + c // w)
。最坏的情况是形式的字符串letter * number
,例如"a" * 10000
.
我没有给出sa
下面模块的代码。它是冗长的,任何后缀数组的实现都会计算相同的东西。的输出suffix_array()
:
sa
是后缀数组,是这样的唯一排列range(n)
,对于所有i
in range(1, n)
, s[sa[i-1]:] < s[sa[i]:]
.
rank
这里不使用。
对于i
in range(1, n)
,给出从和lcp[i]
开始的后缀之间的最长公共前缀的长度。sa[i-1]
sa[i]
为什么会赢?部分原因是它永远不必搜索以相同字母开头的后缀(后缀数组通过构造使它们相邻),并检查重复的块,以及它是否是新的最好的,无论块有多大或重复了多少次。如上所述,这只是对c
and的简单算术w
。
免责声明:后缀数组/树对我来说就像连分数:我可以在必要时使用它们,并且可以惊叹于结果,但它们让我头疼。敏感的,敏感的,敏感的。
def crunch4(s):
from sa import suffix_array
sa, rank, lcp = suffix_array(s)
bestlen, results = 0, []
n = len(s)
for sai in range(n-1):
i = sa[sai]
c = n
for saj in range(sai + 1, n):
c = min(c, lcp[saj])
if not c:
break
j = sa[saj]
w = abs(i - j)
if c < w:
continue
numreps = 1 + c // w
assert numreps > 1
total = w * numreps
if total >= bestlen:
if total > bestlen:
results.clear()
bestlen = total
results.append((min(i, j), w, numreps))
return bestlen, results
一些时间
我将一个普通的英文单词文件读入一个字符串,xs
. 每行一个字:
>>> len(xs)
209755
>>> xs.count('\n')
25481
所以在大约 210K 字节中大约有 25K 字。这些是二次时间算法,所以我没想到它会跑得很快,但crunch2()
几个小时后它仍在运行——当我让它通宵运行时,它仍在运行。
这让我意识到它的update()
功能可以做大量徒劳的工作,使算法整体更像立方时间。所以我修复了那个。然后:
>>> crunch2(xs)
(44, [(63750, 22, 2)])
>>> xs[63750 : 63750+50]
'\nelectroencephalograph\nelectroencephalography\nelec'
这花了大约 38 分钟,这与我预期的差不多。
正则表达式版本crunch3()
只用了不到十分之一秒!
>>> crunch3(xs)
(8, [(19308, 4, 2), (47240, 4, 2)])
>>> xs[19308 : 19308+10]
'beriberi\nB'
>>> xs[47240 : 47240+10]
'couscous\nc'
如前所述,正则表达式版本可能找不到最佳答案,但这里还有其他东西在起作用:默认情况下,“。” 不匹配换行符,所以代码实际上做了很多微小的搜索。文件中约 25K 的换行符中的每一个都有效地结束了本地搜索范围。而是用标志编译正则表达式re.DOTALL
(因此不特殊处理换行符):
>>> crunch3(xs) # with DOTALL
(44, [(63750, 22, 2)])
14 分钟多一点。
最后,
>>> crunch4(xs)
(44, [(63750, 22, 2)])
不到 9 分钟。构建后缀数组的时间只是其中的一小部分(不到一秒)。这实际上令人印象深刻,因为尽管几乎完全“以 C 速度”运行,但并非总是正确的蛮力正则表达式版本较慢。
但这是相对意义上的。在绝对意义上,所有这些仍然是猪慢:-(
注意:下一节中的版本将其缩短到 5 秒以下(!)。
快得多
这个采用完全不同的方法。对于上面的大字典示例,它在不到 5 秒的时间内得到正确答案。
我对此感到相当自豪 ;-) 这是出乎意料的,而且我以前从未见过这种方法。它不进行任何字符串搜索,只是对索引集进行整数运算。
对于表单的输入,它仍然非常缓慢letter * largish_integer
。照原样,只要存在至少两个(不一定相邻,甚至不重叠!)子字符串(正在考虑的当前长度)的副本,它就会一直增加 1。所以,例如,在
'x' * 1000000
它将尝试从 1 到 999999 的所有子字符串大小。
但是,看起来可以通过重复将当前大小加倍(而不是仅仅加 1)、保存类、进行混合形式的二进制搜索来定位存在重复的最大子字符串大小来大大改善这种情况。
我将把它作为一个毫无疑问的乏味练习留给读者。我这的工作都干完了 ;-)
def crunch5(text):
from collections import namedtuple, defaultdict
# For all integers i and j in IxSet x.s,
# text[i : i + x.w] == text[j : j + x.w].
# That is, it's the set of all indices at which a specific
# substring of length x.w is found.
# In general, we only care about repeated substrings here,
# so weed out those that would otherwise have len(x.s) == 1.
IxSet = namedtuple("IxSet", "s w")
bestlen, results = 0, []
# Compute sets of indices for repeated (not necessarily
# adjacent!) substrings of length xs[0].w + ys[0].w, by looking
# at the cross product of the index sets in xs and ys.
def combine(xs, ys):
xw, yw = xs[0].w, ys[0].w
neww = xw + yw
result = []
for y in ys:
shifted = set(i - xw for i in y.s if i >= xw)
for x in xs:
ok = shifted & x.s
if len(ok) > 1:
result.append(IxSet(ok, neww))
return result
# Check an index set for _adjacent_ repeated substrings.
def check(s):
nonlocal bestlen
x, w = s.s.copy(), s.w
while x:
current = start = x.pop()
count = 1
while current + w in x:
count += 1
current += w
x.remove(current)
while start - w in x:
count += 1
start -= w
x.remove(start)
if count > 1:
total = count * w
if total >= bestlen:
if total > bestlen:
results.clear()
bestlen = total
results.append((start, w, count))
ch2ixs = defaultdict(set)
for i, ch in enumerate(text):
ch2ixs[ch].add(i)
size1 = [IxSet(s, 1)
for s in ch2ixs.values()
if len(s) > 1]
del ch2ixs
for x in size1:
check(x)
current_class = size1
# Repeatedly increase size by 1 until current_class becomes
# empty. At that point, there are no repeated substrings at all
# (adjacent or not) of the then-current size (or larger).
while current_class:
current_class = combine(current_class, size1)
for x in current_class:
check(x)
return bestlen, results
而且更快
crunch6()
在我的盒子上将较大的字典示例降至 2 秒以下。它结合了crunch4()
(后缀和 lcp 数组)和crunch5()
(在一组索引中找到所有具有给定步幅的算术级数)的想法。
像 一样crunch5()
,这也循环的次数等于重复的最长子串的长度(重叠与否)的多一倍。因为如果没有长度的重复n
,那么任何大于任何一个的大小都没有n
。这使得在不考虑重叠的情况下查找重复更容易,因为它是一个可利用的限制。当将“胜利”限制为相邻的重复时,就会崩溃。例如,“abcabc”中没有偶数长度为 1 的相邻重复,但有一个长度为 3。这似乎使任何形式的直接二分搜索都无效(大小的相邻重复的存在或不存在n
与存在任何其他大小的相邻重复)。
表单的输入'x' * n
仍然很糟糕。从 1 到 有所有长度的重复n-1
。
观察:我给出的所有程序都会生成所有可能的方法来分解最大长度的重复相邻块。例如,对于一个由 9 个“x”组成的字符串,它表示可以通过重复“x”9 次或重复“xxx”3 次得到。因此,令人惊讶的是,它们也都可以用作分解算法;-)
def crunch6(text):
from sa import suffix_array
sa, rank, lcp = suffix_array(text)
bestlen, results = 0, []
n = len(text)
# Generate maximal sets of indices s such that for all i and j
# in s the suffixes starting at s[i] and s[j] start with a
# common prefix of at least len minc.
def genixs(minc, sa=sa, lcp=lcp, n=n):
i = 1
while i < n:
c = lcp[i]
if c < minc:
i += 1
continue
ixs = {sa[i-1], sa[i]}
i += 1
while i < n:
c = min(c, lcp[i])
if c < minc:
yield ixs
i += 1
break
else:
ixs.add(sa[i])
i += 1
else: # ran off the end of lcp
yield ixs
# Check an index set for _adjacent_ repeated substrings
# w apart. CAUTION: this empties s.
def check(s, w):
nonlocal bestlen
while s:
current = start = s.pop()
count = 1
while current + w in s:
count += 1
current += w
s.remove(current)
while start - w in s:
count += 1
start -= w
s.remove(start)
if count > 1:
total = count * w
if total >= bestlen:
if total > bestlen:
results.clear()
bestlen = total
results.append((start, w, count))
c = 0
found = True
while found:
c += 1
found = False
for s in genixs(c):
found = True
check(s, c)
return bestlen, results
总是很快,并且已发布,但有时是错误的
在生物信息学中,事实证明这是在“串联重复”、“串联阵列”和“简单序列重复”(SSR)的名称下研究的。您可以搜索这些术语以找到相当多的学术论文,其中一些声称最坏情况的线性时间算法。
但这些似乎分为两个阵营:
- 要描述的那种线性时间算法,实际上是错误的:-(
- 算法如此复杂,甚至试图将它们变成可运行的代码都需要奉献精神:-(
在第一个阵营中,有几篇论文可以归结为crunch4()
上述内容,但没有其内部循环。我将在此之后使用代码,crunch4a()
. 这是一个例子:
“SA-SSR:一种基于后缀数组的算法,用于在大型基因序列中详尽有效地发现 SSR。”
皮克特等
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5013907/
crunch4a()
总是很快,但有时会出错。事实上,它为这里出现的每个示例找到至少一个最大重复拉伸,在几分之一秒内解决了较大的字典示例,并且对形式的字符串没有问题'x' * 1000000
。大部分时间都花在构建 suffix 和 lcp 数组上。但它可能会失败:
>>> x = "bcdabcdbcd"
>>> crunch4(x) # finds repeated bcd at end
(6, [(4, 3, 2)])
>>> crunch4a(x) # finds nothing
(0, [])
问题是不能保证相关的后缀在后缀数组中是相邻的。以“b”开头的后缀的顺序如下:
bcd
bcdabcdbcd
bcdbcd
要通过这种方法找到尾随重复块,需要将第一个与第三个进行比较。这就是为什么crunch4()
有一个内部循环,以尝试所有以一个共同字母开头的对。相关对可以由后缀数组中任意数量的其他后缀分隔。但这也使算法成为二次时间。
# only look at adjacent entries - fast, but sometimes wrong
def crunch4a(s):
from sa import suffix_array
sa, rank, lcp = suffix_array(s)
bestlen, results = 0, []
n = len(s)
for sai in range(1, n):
i, j = sa[sai - 1], sa[sai]
c = lcp[sai]
w = abs(i - j)
if c >= w:
numreps = 1 + c // w
total = w * numreps
if total >= bestlen:
if total > bestlen:
results.clear()
bestlen = total
results.append((min(i, j), w, numreps))
return bestlen, results
O(n log n)
这篇论文对我来说很合适,虽然我没有编码:
“使用后缀树简单灵活地检测连续重复”
延斯·斯托伊,丹·古斯菲尔德
https://csiflabs.cs.ucdavis.edu/~gusfield/tcs.pdf
但是,要使用次二次算法需要做出一些妥协。例如,"x" * n
有n-1
子串的形式"x"*2
,n-2
形式"x"*3
,...,所以只有O(n**2)
那些。因此,找到所有这些的任何算法也必然是最好的二次时间。
阅读论文以获取详细信息;-) 您正在寻找的一个概念是“原始”:我相信您只想要形式的重复,S*n
其中S
本身不能表示为较短字符串的重复。所以,例如,"x" * 10
是原始的,但"xx" * 5
不是。
迈出第一步O(n log n)
crunch9()
是我在评论中提到的“蛮力”算法的实现,来自:
“增强的后缀阵列及其在基因组分析中的应用”
易卜拉欣等
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.93.2217&rep=rep1&type=pdf
那里的实现草图只找到“分支串联”重复,我在这里添加了代码来推断任意数量的重复的重复,并包括非分支重复。虽然它仍然是最坏的情况,但对于您在评论中指出O(n**2)
的字符串来说,它比这里的其他任何东西都要快得多。seq
照原样,它复制(除了订单)与这里的大多数其他程序相同的详尽帐户。
该论文继续努力将最坏的情况减少到O(n log n)
,但这大大减慢了它的速度。因此,它会更加努力地战斗。我承认我失去了兴趣;-)
# Generate lcp intervals from the lcp array.
def genlcpi(lcp):
lcp.append(0)
stack = [(0, 0)]
for i in range(1, len(lcp)):
c = lcp[i]
lb = i - 1
while c < stack[-1][0]:
i_c, lb = stack.pop()
interval = i_c, lb, i - 1
yield interval
if c > stack[-1][0]:
stack.append((c, lb))
lcp.pop()
def crunch9(text):
from sa import suffix_array
sa, rank, lcp = suffix_array(text)
bestlen, results = 0, []
n = len(text)
# generate branching tandem repeats
def gen_btr(text=text, n=n, sa=sa):
for c, lb, rb in genlcpi(lcp):
i = sa[lb]
basic = text[i : i + c]
# Binary searches to find subrange beginning with
# basic+basic. A more gonzo implementation would do this
# character by character, never materialzing the common
# prefix in `basic`.
rb += 1
hi = rb
while lb < hi: # like bisect.bisect_left
mid = (lb + hi) // 2
i = sa[mid] + c
if text[i : i + c] < basic:
lb = mid + 1
else:
hi = mid
lo = lb
while lo < rb: # like bisect.bisect_right
mid = (lo + rb) // 2
i = sa[mid] + c
if basic < text[i : i + c]:
rb = mid
else:
lo = mid + 1
lead = basic[0]
for sai in range(lb, rb):
i = sa[sai]
j = i + 2*c
assert j <= n
if j < n and text[j] == lead:
continue # it's non-branching
yield (i, c, 2)
for start, c, _ in gen_btr():
# extend left
numreps = 2
for i in range(start - c, -1, -c):
if all(text[i+k] == text[start+k] for k in range(c)):
start = i
numreps += 1
else:
break
totallen = c * numreps
if totallen < bestlen:
continue
if totallen > bestlen:
bestlen = totallen
results.clear()
results.append((start, c, numreps))
# add non-branches
while start:
if text[start - 1] == text[start + c - 1]:
start -= 1
results.append((start, c, numreps))
else:
break
return bestlen, results
赚取奖励积分;-)
对于某些技术含义;-)crunch11()
是最坏情况 O(n log n)。除了 suffix 和 lcp 数组,这还需要rank
数组 ,sa
的逆:
assert all(rank[sa[i]] == sa[rank[i]] == i for i in range(len(sa)))
正如代码注释所指出的,它还依赖 Python 3 来提高速度(range()
行为)。这很浅,但重写会很乏味。
描述这个的论文有几个错误,所以如果这个代码与你读到的不完全匹配,不要翻出来。完全按照他们所说的去做,它会失败。
也就是说,代码变得非常复杂,我不能保证没有错误。它适用于我尝试过的一切。
表单的输入'x' * 1000000
仍然不快,但显然不再是二次时间。例如,将同一个字母重复一百万次的字符串在近 30 秒内完成。这里的大多数其他程序永远不会结束;-)
编辑:改为genlcpi()
使用半开放的 Python 范围;主要对外观进行了修改crunch11()
;'x' * 1000000
添加了“提前退出”,在最坏(如)情况下节省了大约三分之一的时间。
# Generate lcp intervals from the lcp array.
def genlcpi(lcp):
lcp.append(0)
stack = [(0, 0)]
for i in range(1, len(lcp)):
c = lcp[i]
lb = i - 1
while c < stack[-1][0]:
i_c, lb = stack.pop()
yield (i_c, lb, i)
if c > stack[-1][0]:
stack.append((c, lb))
lcp.pop()
def crunch11(text):
from sa import suffix_array
sa, rank, lcp = suffix_array(text)
bestlen, results = 0, []
n = len(text)
# Generate branching tandem repeats.
# (i, c, 2) is branching tandem iff
# i+c in interval with prefix text[i : i+c], and
# i+c not in subinterval with prefix text[i : i+c + 1]
# Caution: this pragmatically relies on that, in Python 3,
# `range()` returns a tiny object with O(1) membership testing.
# In Python 2 it returns a list - ahould still work, but very
# much slower.
def gen_btr(text=text, n=n, sa=sa, rank=rank):
from itertools import chain
for c, lb, rb in genlcpi(lcp):
origlb, origrb = lb, rb
origrange = range(lb, rb)
i = sa[lb]
lead = text[i]
# Binary searches to find subrange beginning with
# text[i : i+c+1]. Note we take slices of length 1
# rather than just index to avoid special-casing for
# i >= n.
# A more elaborate traversal of the lcp array could also
# give us a list of child intervals, and then we'd just
# need to pick the right one. But that would be even
# more hairy code, and unclear to me it would actually
# help the worst cases (yes, the interval can be large,
# but so can a list of child intervals).
hi = rb
while lb < hi: # like bisect.bisect_left
mid = (lb + hi) // 2
i = sa[mid] + c
if text[i : i+1] < lead:
lb = mid + 1
else:
hi = mid
lo = lb
while lo < rb: # like bisect.bisect_right
mid = (lo + rb) // 2
i = sa[mid] + c
if lead < text[i : i+1]:
rb = mid
else:
lo = mid + 1
subrange = range(lb, rb)
if 2 * len(subrange) <= len(origrange):
# Subrange is at most half the size.
# Iterate over it to find candidates i, starting
# with wa. If i+c is also in origrange, but not
# in subrange, good: then i is of the form wwx.
for sai in subrange:
i = sa[sai]
ic = i + c
if ic < n:
r = rank[ic]
if r in origrange and r not in subrange:
yield (i, c, 2, subrange)
else:
# Iterate over the parts outside subrange instead.
# Candidates i are then the trailing wx in the
# hoped-for wwx. We win if i-c is in subrange too
# (or, for that matter, if it's in origrange).
for sai in chain(range(origlb, lb),
range(rb, origrb)):
ic = sa[sai] - c
if ic >= 0 and rank[ic] in subrange:
yield (ic, c, 2, subrange)
for start, c, numreps, irange in gen_btr():
# extend left
crange = range(start - c, -1, -c)
if (numreps + len(crange)) * c < bestlen:
continue
for i in crange:
if rank[i] in irange:
start = i
numreps += 1
else:
break
# check for best
totallen = c * numreps
if totallen < bestlen:
continue
if totallen > bestlen:
bestlen = totallen
results.clear()
results.append((start, c, numreps))
# add non-branches
while start and text[start - 1] == text[start + c - 1]:
start -= 1
results.append((start, c, numreps))
return bestlen, results