Python，巨大迭代性能问题

Question

我正在遍历3个单词，每个单词大约有500万个字符，并且我想找到标识每个单词的20个字符序列。也就是说，我想在一个单词中找到所有长度为20的唯一序列。我的问题是，我编写的代码运行时间非常长。我甚至在整夜运行程序后都无法完成一个单词。

下面的函数接受一个包含字典的列表，每个字典都包含5百万长单词中每个可能的20个单词和它们的位置。

如果有人有优化这个问题的想法，我将非常感激，我不知道如何继续...

以下是我的代码示例：

def findUnique(list):
    # Takes a list with dictionaries and compairs each element in the dictionaries
    # with the others and puts all unique element in new dictionaries and finally
    # puts the new dictionaries in a list.
    # The result is a list with (in this case) 3 dictionaries containing all unique
    # sequences and their locations from each string.
    dicList=[]
    listlength=len(list)
    s=0
    valuelist=[]
    for i in list:
        j=i.values()
        valuelist.append(j)
    while s<listlength:
        currdic=list[s]
        dic={}
        for key in currdic:
            currval=currdic[key]
            test=True
            n=0
            while n<listlength:
                if n!=s:
                    if currval in valuelist[n]: #this is where it takes to much time
                        n=listlength
                        test=False
                    else:
                        n+=1
                else:
                    n+=1
            if test:
                dic[key]=currval
        dicList.append(dic)
        s+=1
    return dicList

- user236226

---

n 的平方乘以字典大小。难怪速度很慢。 - S.Lott

+1 非常感谢您实际发布代码，而不是要求我们读心术！ - PaulMcG

也许可以看一下这篇论文，它讨论了使用布隆过滤器来完成非常类似的任务：http://www.serpentine.com/bos/files/padl09.pdf。由于这篇论文是关于 Haskell 的，因此作为评论发布，希望对你有所帮助。 - Andrew Y

Answer 1

def slices(seq, length, prefer_last=False):
  unique = {}
  if prefer_last: # this doesn't have to be a parameter, just choose one
    for start in xrange(len(seq) - length + 1):
      unique[seq[start:start+length]] = start
  else: # prefer first
    for start in xrange(len(seq) - length, -1, -1):
      unique[seq[start:start+length]] = start
  return unique

# or find all locations for each slice:
import collections
def slices(seq, length):
  unique = collections.defaultdict(list)
  for start in xrange(len(seq) - length + 1):
    unique[seq[start:start+length]].append(start)
  return unique

这个函数（目前在我的iter_util模块中）的时间复杂度为O(n)（n为每个单词的长度），您可以使用set(slices(..))（与差分等集合运算一起使用）以获取所有单词中唯一的片段（下面有例子）。如果您不想跟踪位置，您也可以编写该函数来返回一个集合。内存使用量将很高（尽管仍然是O(n)，但因素很大），特别是当长度只有20时，可能会得到缓解，但不会太多，如果使用存储基本序列（字符串）加上开始和停止（或开始和长度）的特殊的“懒惰切片”类。

打印唯一的片段:

a = set(slices("aab", 2)) # {"aa", "ab"}
b = set(slices("abb", 2)) # {"ab", "bb"}
c = set(slices("abc", 2)) # {"ab", "bc"}
all = [a, b, c]
import operator
a_unique = reduce(operator.sub, (x for x in all if x is not a), a)
print a_unique # {"aa"}

包括以下位置：

a = slices("aab", 2)
b = slices("abb", 2)
c = slices("abc", 2)
all = [a, b, c]
import operator
a_unique = reduce(operator.sub, (set(x) for x in all if x is not a), set(a))
# a_unique is only the keys so far
a_unique = dict((k, a[k]) for k in a_unique)
# now it's a dict of slice -> location(s)
print a_unique # {"aa": 0} or {"aa": [0]}
               # (depending on which slices function used)

---

在一个更贴近你条件的测试脚本中，使用随机生成的5m个字符的单词和长度为20的切片，内存使用量非常高，我的测试脚本很快就达到了1G主内存限制，并开始把虚拟内存搞崩。此时，Python几乎没有花费任何时间在CPU上，我终止了它。将切片长度或单词长度（因为我使用的是完全随机的单词，这会减少重复并增加内存使用）缩小以适应主内存，它运行不到一分钟。这种情况再加上原始代码中的O(n ** 2)，将需要很长时间，这就是为什么算法的时间和空间复杂度都很重要。

import operator
import random
import string

def slices(seq, length):
  unique = {}
  for start in xrange(len(seq) - length, -1, -1):
    unique[seq[start:start+length]] = start
  return unique

def sample_with_repeat(population, length, choice=random.choice):
  return "".join(choice(population) for _ in xrange(length))

word_length = 5*1000*1000
words = [sample_with_repeat(string.lowercase, word_length) for _ in xrange(3)]
slice_length = 20
words_slices_sets = [set(slices(x, slice_length)) for x in words]
unique_words_slices = [reduce(operator.sub,
                              (x for x in words_slices_sets if x is not n),
                              n)
                       for n in words_slices_sets]
print [len(x) for x in unique_words_slices]

- Roger Pate

Answer 2

让我从罗杰·佩特的答案继续下去。如果内存是一个问题，我建议不要使用字符串作为字典的键，而应该使用字符串的哈希值。这将节省存储额外拷贝字符串作为键的成本（在最坏情况下，每个“单词”的存储量最多为20倍）。

import collections
def hashed_slices(seq, length, hasher=None):
  unique = collections.defaultdict(list)
  for start in xrange(len(seq) - length + 1):
    unique[hasher(seq[start:start+length])].append(start)
  return unique

（如果你真的想要变得高级一些，你可以使用滚动哈希，不过你需要改变函数。）

现在，我们可以将所有哈希值组合起来：

unique = []  # Unique words in first string

# create a dictionary of hash values -> word index -> start position
hashed_starts = [hashed_slices(word, 20, hashing_fcn) for word in words]
all_hashed = collections.defaultdict(dict)
for i, hashed in enumerate(hashed_starts) :
   for h, starts in hashed.iteritems() :
     # We only care about the first word
     if h in hashed_starts[0] :
       all_hashed[h][i]=starts

# Now check all hashes
for starts_by_word in all_hashed.itervalues() :
  if len(starts_by_word) == 1 :
    # if there's only one word for the hash, it's obviously valid
    unique.extend(words[0][i:i+20] for i in starts_by_word.values())
  else :
    # we might have a hash collision
    candidates = {}
    for word_idx, starts in starts_by_word.iteritems() :
      candidates[word_idx] = set(words[word_idx][j:j+20] for j in starts)
    # Now go that we have the candidate slices, find the unique ones
    valid = candidates[0]
    for word_idx, candidate_set in candidates.iteritems() :
      if word_idx != 0 :
        valid -= candidate_set
    unique.extend(valid)

我尝试扩展它以执行所有三个任务。虽然这是可能的，但复杂性会减弱算法效果。

请注意，我还没有对此进行测试。此外，代码可能有很多可以简化的地方，但是算法是有道理的。难点在于选择哈希。太多冲突会导致你毫无收获，而太少则会遇到内存问题。如果你只处理DNA基本代码，可以将20个字符的字符串哈希成一个40位数，仍然不会发生冲突。因此，每个切片将占用近四分之一的内存。这将节省大约250 MB的内存，比Roger Pate的答案节省更多。

代码仍然是O(N^2)，但常数应该要低得多。

- AFoglia

Answer 3

让我们试着改进Roger Pate的优秀回答。

首先，让我们保留集合而不是字典——它们本来就管理唯一性。

其次，由于我们很可能会比耗尽CPU时间（和耐心）更快地耗尽内存，因此我们可以为了内存效率而牺牲CPU效率。因此，也许只尝试以特定字母开头的20个字符串。对于DNA，这将要求减少75%。

seqlen = 20
maxlength = max([len(word) for word in words])
for startletter in letters:
    for letterid in range(maxlength):
        for wordid,word in words:
            if (letterid < len(word)):
                letter = word[letterid]
                if letter is startletter:
                    seq = word[letterid:letterid+seqlen]
                    if seq in seqtrie and not wordid in seqtrie[seq]:
                        seqtrie[seq].append(wordid)

或者，如果那样仍然太多内存，我们可以为每个可能的起始对进行遍历（与DNA的4次不同要进行16次），或者每3个（64次）等。

- Phil H

Answer 4

你说你有一个500万个字符长的“单词”，但我很难相信这在通常意义下是一个单词。

如果您可以提供更多关于输入数据的信息，可能会有一个具体的解决方案。

例如，英文文本（或其他任何书面语言）可能足够重复，以至于trie是可用的。但最坏的情况下，它将耗尽内存构造所有256^20个键。了解您的输入数据可以让我们做出更好的决策。

---

编辑

我查看了一些基因组数据，使用硬编码的 [acgt]->[0123] 映射和每个 trie 节点 4 个孩子节点来检查这个想法的效果。

1. 腺病毒2：35,937bp -> 使用469,339个 trie 节点得到35,899个不同的20个碱基序列

2. 肠道杆菌噬菌体 lambda：48,502bp -> 使用529,384个 trie 节点得到40,921种不同的20个碱基序列。

在这两个数据集中，我没有发现任何冲突，但你的数据可能存在更多的冗余和/或重叠。你需要亲自试一试才能知道。

如果你可以得到有用数量的冲突，你可以尝试将三个输入一起遍历，构建一个单独的 trie，记录每个叶子节点的来源，并在遍历时从 trie 中删除冲突。

如果你无法找到删除键的方法，可以尝试使用更紧凑的表示方式。例如，你只需要 2 位来存储 [acgt]/[0123]，这可能会节省空间，但代价是稍微更复杂的代码。

我认为你不能仅仅通过暴力方式解决这个问题 - 你需要找到一些方法来缩小问题的规模，这取决于你的领域知识。

- Useless

---

这个问题被标记为“生物信息学”，所以很可能那些不是英语单词，而是DNA序列。 - Roberto Bonvallet

就是这样。如果只涉及4个字符，那么它可能仍然有效... 4^20约等于10^12，因此只有在有许多常见的子树可折叠时才能使用。我对DNA了解不够，无法猜测。 - Useless