如何加速 itertools 组合？答案

【问题标题】：How to speed up itertools combinations?如何加速 itertools 组合？
【发布时间】：2015-07-15 14:50:43
【问题描述】：

作为我自己对基因型网络研究的一部分，我有一段代码，我试图在其中从字符串列表构建一个单差网络。流程是这样的：

遍历字符串的所有成对组合。
如果字符串相差一个位置，则在网络中在它们之间绘制一条边。
如果它们没有一个位置差异，则通过。

我现在的代码块是这样的：

from itertools import combinations
import Levenshtein as lev # a package that wraps a C-implemented levenshtein distance
import networkx as nx
strings = [...list of strings...]

G = nx.Graph()

for (s1, s2) in combinations(strings, 2):
    if s1 not in G.nodes():
        G.add_node(s1)
    if s2 not in G.nodes():
        G.add_node(s2)

    if lev.distance(s1, s2) == 1:
        G.add_edge(s1, s2)

显然没有办法提高图构造过程的计算复杂度——它总是O(n**2)。至少，在我有限的知识范围内，我是这么认为的——也许我错了？

也就是说，考虑到我需要进行的比较次数的正常规模（约 2000-5000 次），如果总体上我能获得几个数量级的加速，那么实际计算时间将仍然可以接受 - 使用当前的 Python 实现，构建图形需要〜几天。使用正确的导入（下文未说明），我尝试了下面的 Cython 实现，但无法弄清楚如何使其更快：

cpdef cython_genotype_network(sequences):

    G = nx.Graph()
    cdef:
        unicode s1
        unicode s2

    for (s1, s2) in combinations(sequences, 2):
        if lev.distance(s1, s2) == 1:
            G.add_edge(s1, s2)

    return G

特别是，Cython 期望 bytes，而不是 str 用于 s1 和 s2。该代码块引发错误。

所以...我来回答我的两个问题：

Q1：Cython 实施会有帮助吗？以及如何修复 bytes vs. str 错误？
Q2：是否可以改用numpy 来解决这个问题？从numpy 矩阵转换为NetworkX 很容易；但是，我似乎无法弄清楚如何在每行和每列对应一个字符串的 n×n 空矩阵中广播 Levenshtein 距离函数。

更新 1：如何生成示例数据

要生成字符串：

from random import choice

def create_random_nucleotides_python(num_nucleotides):
    """
    Creates random nucleotides of length num_nucleotides.
    """

    sequence = ''

    letters = ['A', 'T', 'G', 'C']

    for i in range(num_nucleotides):
        sequence = sequence + choice(letters)

    return sequence


def mutate_random_position(string):
    """
    Mutates one position in the nucleotide sequence at random.
    """

    positions = [i for i in range(len(string))]
    pos_to_mut = choice(positions)

    letters = ['A', 'T', 'G', 'C']

    new_string = ''
    for i, letter in enumerate(string):
        if i == pos_to_mut:
            new_string = new_string + choice(letters)
        else:
            new_string = new_string + letter

    return new_string


# Create 100 Python srings by mutating a first sequence, then randomly choosing stuff to mutate a single position.
base_sequence = create_random_nucleotides_python(1000)

sequences = [base_sequence]

for i in range(99):
    sequence = choice(sequences)
    mutseq = mutate_random_position(sequence)
    sequences.append(mutseq)

【问题讨论】：

什么是.nodes()？它不是一个集合，它应该是
基于networkx文档，Graph.nodes()返回一个列表，但底层实现是一个字典。
那是O(n) 每次查找，我认为节点是键？
是的，我认为这是正确的。
好吧，if s1 not in G 或 if G.has_node(s1) 应该给你 0(1) 次查找，我想 G.has_node 是 if s1 not in G 的语法糖，但我没有检查

标签： python numpy combinations cython levenshtein-distance

【解决方案1】：

关于复杂性：

您正在考虑每对字符串。你不需要那个。您可以考虑每个字符串的所有 1 距离字符串：

# I would start by populating the whole graph:
for s1 in strings:
    if s1 not in G.nodes():
        G.add_node(s1)
# Then compute the leven-1:
for s1 in strings:
    for s2 in variations(s1):
        if s2 in G.nodes():
            G.add_edge(s1, s2)

现在你只需要一个比O(n)短的variations(string)：

这将返回距离为 1 的所有变体。（只有 1 个编辑|删除|插入）

def variations(string):
    for i in range(len(string)):
        # delete
        yield string[:i] + string[i+1:]
        # edit
        for l in 'ATGC':
            yield string[:i] + l + string[i+1:]
        # insert
        for l in 'ATGC':
            yield string[:i] + l + string[i:]

    # insert at the end
    for l in 'ATGC':
        yield string + l

现在，它的复杂度是 O(m^2)（因为字符串 concat），其中 m 是最长序列的大小。如果它是已知的，它就是一个常数，而现在就是O(1)

如果序列的大小都相同，则只能计算编辑。

否则，您可以将序列从大到小排序，并且只计算编辑和删除。

或者，您可以按大小对字符串进行排序，而不是将所有字符串与所有其他字符串进行比较，而是比较大小差异为

【讨论】：

有趣的方法！这些序列实际上是可变长度的氨基酸序列（但长度差异不大）。
在字符串大小这一点上，如果我们知道最大字符串长度，复杂度是O(1)？
@ericmjl 因为它变成了一个常数。如果您的字符串大小达到最大值，例如 10（但与 10000 相同），那么您知道字符串有 100 种变体，无论池中有多少字符串。这意味着整个算法是O(n)，其中 n 是字符串的数量，因为如果你有两倍的字符串，算法将占用两倍的大小。但是，如果字符串的大小未知，则表示计算中的变量
现在，从您的 mutate_random_position 函数来看，字符串的大小似乎没有改变，显然是 1000。
这是正确的，它只适用于合成数据。现实世界的数据要混乱得多 :) 我只在合成数据上测试不同的算法。

【解决方案2】：

关于cython 代码可能更快的一些想法

for (s1, s2) in combinations(sequences, 2):
    if lev.distance(s1, s2) == 1:
        G.add_edge(s1, s2)

我认为itertools.combinations 已编译，因此它本身很快。你调用它一次，所以for ... 部分可能和它来的一样快。尽管如此，在cython 中循环访问sequences 并不难。

lev.distance 似乎使用已编译的代码。是否可以直接导入和调用该代码？查看lev 来源。

我认为lev.distance 最终会调用带有签名的 c 函数：

distance_py(PyObject *self, PyObject *args)

基本上是

lev_edit_distance(size_t len1, const lev_byte *string1,
              size_t len2, const lev_byte *string2,
              int xcost)

https://github.com/ztane/python-Levenshtein/blob/master/Levenshtein/_levenshtein.c

G.add_edge 是做什么的？是否有更简单的数据结构可以收集这些“边缘”。也许作为元组列表？

networkx 是否提供批量添加节点和边的方法？您是否总是必须使用add_node 和add_edge？或者你能给它一个节点列表和节点对（元组）列表吗？

看起来 graph 是 Python 字典或字典字典。

networkx 允许您通过元组列表添加一堆边：

>>> G.add_edges_from([(1,2),(1,3)])

http://networkx.github.io/documentation/latest/tutorial/tutorial.html#edges

在python代码中，这一段可能是不必要的：

if s1 not in G.nodes():
    G.add_node(s1)
if s2 not in G.nodes():
    G.add_node(s2)

只要做：

G.add_nodes_from(strings)

http://networkx.github.io/documentation/latest/tutorial/tutorial.html#nodes

【讨论】：