提问人:Stand with Gaza 提问时间:5/13/2020 更新时间:5/25/2020 访问量:1205
查找最长相邻重复非重叠子字符串
Find longest adjacent repeating non-overlapping substring
问:
(这个问题不是关于音乐的,但我以音乐为例 一个用例。
在音乐中,构建乐句的常用方法是作为音符序列 其中中间部分重复一次或多次。因此,这句话 由引子、循环部分和结尾组成。这是一个 例:
[ E E E F G A F F G A F F G A F C D ]
我们可以“看到”前奏是 [ E E E] 重复部分是 [ F G A F ] 和结尾是 [ C D ]。因此,拆分列表的方法是
[ [ E E E ] 3 [ F G A F ] [ C D ] ]
其中第一项是介绍,第二项是 重复部分是重复的,第三部分是结尾。
我需要一种算法来执行这样的拆分。
但有一点需要注意,那就是可能有多种方法可以 拆分列表。例如,上面的列表可以分为:
[ [ E E E F G A ] 2 [ F F G A ] [ F C D ] ]
但这是一个更糟糕的分裂,因为前奏和结尾更长。所以 该算法的标准是找到最大化 循环部分的长度,并最小化 前奏和结尾。这意味着正确的拆分
[ A C C C C C C C C C A ]
是
[ [ A ] 9 [ C ] [ A ] ]
因为前奏和结尾的总长度是 2,而长度 循环部分是 9。
此外,虽然前奏和结尾可以是空的,但只有“真实”重复是空的 允许。因此,以下拆分将被禁止:
[ [ ] 1 [ E E E F G A F F G A F F G A F C D ] [ ] ]
可以把它看作是找到最佳的“压缩”, 序列。请注意,某些序列中可能没有任何重复:
[ A B C D ]
对于这些堕落的情况,任何合理的结果是允许的。
这是我的算法实现:
def find_longest_repeating_non_overlapping_subseq(seq):
candidates = []
for i in range(len(seq)):
candidate_max = len(seq[i + 1:]) // 2
for j in range(1, candidate_max + 1):
candidate, remaining = seq[i:i + j], seq[i + j:]
n_reps = 1
len_candidate = len(candidate)
while remaining[:len_candidate] == candidate:
n_reps += 1
remaining = remaining[len_candidate:]
if n_reps > 1:
candidates.append((seq[:i], n_reps,
candidate, remaining))
if not candidates:
return (type(seq)(), 1, seq, type(seq)())
def score_candidate(candidate):
intro, reps, loop, outro = candidate
return reps - len(intro) - len(outro)
return sorted(candidates, key = score_candidate)[-1]
我不确定它是否正确,但它通过了我做过的简单测试 描述。它的问题在于它的速度很慢。我看过了 在后缀树,但它们似乎不适合我的用例,因为 我所追求的子字符串应该是不重叠和相邻的。
答:
看起来您要做的几乎是 LZ77 压缩算法。您可以根据我链接到的 Wikipedia 文章中的参考实现检查您的代码。
评论
这是我对你所说的内容的实现。它与你的非常相似,但它跳过了子字符串,这些子字符串已被检查为先前子字符串的重复。
from collections import namedtuple
SubSequence = namedtuple('SubSequence', ['start', 'length', 'reps'])
def longest_repeating_subseq(original: str):
winner = SubSequence(start=0, length=0, reps=0)
checked = set()
subsequences = ( # Evaluates lazily during iteration
SubSequence(start=start, length=length, reps=1)
for start in range(len(original))
for length in range(1, len(original) - start)
if (start, length) not in checked)
for s in subsequences:
subseq = original[s.start : s.start + s.length]
for reps, next_start in enumerate(
range(s.start + s.length, len(original), s.length),
start=1):
if subseq != original[next_start : next_start + s.length]:
break
else:
checked.add((next_start, s.length))
s = s._replace(reps=reps)
if s.reps > 1 and (
(s.length * s.reps > winner.length * winner.reps)
or ( # When total lengths are equal, prefer the shorter substring
s.length * s.reps == winner.length * winner.reps
and s.reps > winner.reps)):
winner = s
# Check for default case with no repetitions
if winner.reps == 0:
winner = SubSequence(start=0, length=len(original), reps=1)
return (
original[ : winner.start],
winner.reps,
original[winner.start : winner.start + winner.length],
original[winner.start + winner.length * winner.reps : ])
def test(seq, *, expect):
print(f'Testing longest_repeating_subseq for {seq}')
result = longest_repeating_subseq(seq)
print(f'Expected {expect}, got {result}')
print(f'Test {"passed" if result == expect else "failed"}')
print()
if __name__ == '__main__':
test('EEEFGAFFGAFFGAFCD', expect=('EEE', 3, 'FGAF', 'CD'))
test('ACCCCCCCCCA', expect=('A', 9, 'C', 'A'))
test('ABCD', expect=('', 1, 'ABCD', ''))
为我传递了你的所有三个示例。这似乎是可能有很多奇怪的边缘情况的事情,但考虑到它是一种优化的蛮力,它可能更像是更新规范的问题,而不是修复代码本身中的错误。
这是一种明显的二次时间方法,但常数因子相对较低,因为它不会构建任何子字符串对象,除了长度为 1 的子字符串对象。结果是一个 2 元组,
bestlen, list_of_results
其中 是重复相邻块的最长子串的长度,每个结果是一个 3 元组,bestlen
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 次。
它不会返回前奏或结尾,因为从它返回的内容中推断出这些是微不足道的:
- 介绍是 .
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
因此,所有块都必须是 block1 的重复项。N
因此,代码会反复移动(复制)原始字符串留下一个字符,然后从左到右移动两个字符,以标识相等字符的最长段。这只需要一次比较一对字符。为了使“左移”有效(恒定时间),字符串的副本存储在 .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)
- 给定 ints 和 .
i
j
0 <= i < j < n
然后,如果 和 是 和 之间共有的前导字符数,则重复(长度)块,从 开始,总共 1 次。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
是后缀数组,使得 对于所有 in , 的唯一排列。range(n)
i
range(1, n)
s[sa[i-1]:] < s[sa[i]:]
rank
此处不使用。对于 in ,给出从 和 开始的后缀之间最长公共前缀的长度。
i
range(1, n)
lcp[i]
sa[i-1]
sa[i]
为什么会赢?部分原因是它永远不必搜索以相同字母开头的后缀(后缀数组,通过构造,使它们相邻),并且检查重复的块,以及它是否是新的最佳块,无论块有多大或重复多少次,都需要很小的恒定时间。如上所述,这只是 和 上的微不足道的算术运算。c
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
所以大约 25K 字,大约 210K 字节。这些是二次时间算法,所以我没想到它会跑得很快,但下班后仍在运行 - 当我让它过夜时仍在运行。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 秒的时间内获得正确答案。
我为此感到自豪;-)这是出乎意料的,我以前从未见过这种方法。它不做任何字符串搜索,只对索引集进行整数算术。
对于表单的输入,它仍然非常慢。照原样,只要子字符串(当前考虑的长度)至少存在两个(不一定是相邻的,甚至不重叠的)副本,它就会继续上升 1。所以,例如,在letter * largish_integer
'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()
将 Largish 字典示例放在我的盒子上不到 2 秒。它结合了(后缀和 lcp 数组)和(在一组索引中查找具有给定步幅的所有算术级数)的想法。crunch4()
crunch5()
像 一样,这也循环了多次,等于重复的最长子字符串的长度(重叠与否)。因为如果没有长度的重复,则没有大于任何一个大小的重复。这使得在不考虑重叠的情况下查找重复变得更加容易,因为这是一个可利用的限制。当将“胜利”限制为相邻的重复时,这就会被打破。例如,“abcabc”中没有偶数长度为 1 的相邻重复项,但有一个长度为 3 的重复项。这似乎使任何形式的直接二叉搜索都是徒劳的(大小的相邻重复的存在与否并不能说明任何其他大小的相邻重复的存在)。crunch5()
n
n
n
表单的输入仍然很悲惨。从 1 到 的所有长度都有重复。'x' * n
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发现。”
Pickett et alia
crunch4a()
总是很快,但有时是错误的。事实上,它为此处出现的每个示例找到至少一个最大重复拉伸,在几分之一秒内解决了大型字典示例,并且对 .大部分时间都花在构建后缀和 lcp 数组上。但它可能会失败:'x' * 1000000
>>> 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)
这篇论文在我看来是正确的,尽管我没有对它进行编码:
“使用后缀树简单灵活地检测连续重复”
延斯·斯托耶,丹·古斯菲尔德
不过,要获得亚二次算法,需要做出一些妥协。例如,具有 、 形式的子字符串 、 ...,因此只有这些子字符串。因此,任何找到所有这些的算法也必然是充其量的二次时间。"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
赚取奖励积分 ;-)
对于一些技术意义;-) 是最坏情况 O(n log n)。除了后缀和 lcp 数组之外,这还需要数组 的 inverse:crunch11()
rank
sa
assert all(rank[sa[i]] == sa[rank[i]] == i for i in range(len(sa)))
正如代码注释所指出的,它还依赖于 Python 3 来提高速度(行为)。这很肤浅,但重写起来会很乏味。range()
描述这一点的论文有几个错误,所以如果这段代码与你读到的内容不完全匹配,请不要翻出来。完全按照他们所说的去做,它就会失败。
也就是说,代码变得越来越复杂,我不能保证没有错误。它适用于我尝试过的所有东西。
表单的输入仍然不快,但显然不再是二次时间。例如,一个字符串重复同一个字母一百万次,在将近 30 秒内完成。这里的大多数其他程序永远不会结束;-)'x' * 1000000
编辑:更改为使用半开放的Python范围;主要对外观进行了更改;添加了“提前出局”,在最坏(类似)的情况下节省了大约三分之一的时间。genlcpi()
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
评论
aaaaaa
[ [] 6 a [] ]
[ [] 3 aa [] ]
[ [] 3 [ [] 2 a [] ] [] ]
crunch6
crunch6()
seq[434 : 434 + 2160] == seq[2596 : 2596 + 2160]
O(n log n)
n
O(log n)
O(n log n)
S+S
O(n**2)
O(n log n)
上一个:在间隔不均匀的时间序列中检测峰值
评论