9

我正在对 3 个单词进行迭代，每个单词长约 500 万个字符，我想找到识别每个单词的 20 个字符的序列。也就是说，我想在一个单词中找到该单词唯一的所有长度为 20 的序列。我的问题是我编写的代码需要很长时间才能运行。我什至连一个单词都没有完成运行我的程序。

下面的函数获取一个包含字典的列表，其中每个字典包含 20 个可能的单词及其来自 500 万个长单词之一的位置。

如果有人知道如何优化它，我将非常感激，我不知道如何继续......

这是我的代码示例：

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
4

4 回答 4

10 
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)，只是一个很大的因素），可能会通过存储基本序列（字符串）的特殊“惰性切片”类以及开始和停止（或开始和长度）。

打印独特的切片：

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]
于 2009-12-21T18:42:06.653 回答
0

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

如果您可以提供有关输入数据的更多信息，则可能会提供特定的解决方案。

例如，英文文本（或任何其他书面语言）可能具有足够的重复性，以至于trie可以使用。然而，在最坏的情况下，它会耗尽构建所有 256^20 个键的内存。了解您的输入会有所不同。

编辑

我查看了一些基因组数据，看看这个想法是如何叠加起来的，使用硬编码的 [acgt]->[0123] 映射和每个 trie 节点 4 个子节点。

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

  2. enterobacteria phage lambda：48,502bp -> 40,921 个不同的 20 碱基序列，使用 529,384 个 trie 节点。

我没有遇到任何冲突，无论是在两个数据集中还是在两个数据集之间，尽管您的数据中可能存在更多冗余和/或重叠。你得试一试才能看到。

如果你确实得到了有用数量的碰撞，你可以尝试将三个输入一起走，构建一个树，记录每个叶子的起源，并在你去的时候修剪树中的碰撞。

如果您找不到修剪密钥的方法，您可以尝试使用更紧凑的表示。例如，您只需要2 位来存储 [acgt]/[0123]，这可能会以稍微复杂的代码为代价来节省空间。

我不认为你可以暴力破解 - 你需要找到一些方法来减少问题的规模，这取决于你的领域知识。

于 2009-12-22T14:56:16.647 回答
0

让我以罗杰·佩特的回答为基础。如果内存是一个问题，我建议不要使用字符串作为字典的键，而是可以使用字符串的哈希值。这将节省将字符串的额外副本存储为键的成本（最坏的情况是存储单个“单词”的 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 位数字，并且仍然没有冲突。所以切片将占用近四分之一的内存。在 Roger Pate 的回答中，这将节省大约 250 MB 的内存。

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

于 2009-12-22T21:53:46.433 回答
0

让我们尝试改进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)

或者，如果内存仍然太多，我们可以遍历每个可能的起始对（16 次而不是 DNA 的 4 次），或者每 3 次（64 次）等。

于 2009-12-28T22:39:47.090 回答