logo标签列表

Python numpy：在numpy 2-D数组的每一对列上执行函数？

pythonfunctionnumpy
6

我将尝试将一个函数应用于numpy数组中每一对列（每一列都是个体的基因型）。

例如：

[48]: g[0:10,0:10]

array([[ 1,  1,  1,  1,  1,  1,  1,  1,  1, -1],
      [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1],
      [ 1,  1,  1,  1,  1,  1, -1,  1,  1,  1],
      [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1],
      [ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
      [ 1,  1,  1,  1,  1,  1,  1,  1,  1, -1],
      [-1, -1,  0, -1, -1, -1, -1, -1, -1,  0],
      [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1],
      [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1],
      [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1]], dtype=int8)

我的目标是生成一个距离矩阵d，使得d中的每个元素都是比较g中每一列之间的成对距离。

d[0,1] = func(g[:,0], g[:,1])

任何想法都是很棒的！谢谢！
3
距离函数如何工作？ - Mazdak
比较两个numpy数组需要使用以下代码：def count_snp_diffs(x, y): return np.count_nonzero((x != y) & (x >= 0) & (y >= 0)) . 本代码来自http://alimanfoo.github.io/ - ksw
3个回答
2
您可以简单地定义该函数如下：
def count_snp_diffs(x, y): 
    return np.count_nonzero((x != y) & (x >= 0) & (y >= 0),axis=0)

然后使用由itertools.combinations生成的数组作为索引来调用它，以获取所有可能的列组合：

combinations = np.array(list(itertools.combinations(range(g.shape[1]),2)))
dist = count_snp_diffs(g[:,combinations[:,0]], g[:,combinations[:,1]])

此外，如果输出必须存储在一个矩阵中（对于大的g不推荐这样做，因为只有上三角部分会被填充，其余部分都是无用信息），可以使用相同的技巧来实现：

d = np.zeros((g.shape[1],g.shape[1]))
combinations = np.array(list(itertools.combinations(range(g.shape[1]),2)))
d[combinations[:,0],combinations[:,1]] = count_snp_diffs(g[:,combinations[:,0]], g[:,combinations[:,1]])

现在，d[i,j]返回列i和j之间的距离（而d[j,i]为零）。这种方法依赖于数组可以使用包含重复索引的列表或数组进行索引的事实：
a = np.arange(3)+4
a[[0,1,1,1,0,2,1,1]]
# Out
# [4, 5, 5, 5, 4, 6, 5, 5]

这里有一个逐步解释正在发生的事情。
调用g[:,combinations[:,0]]访问排列的第一列中的所有列，生成一个新数组，该数组与使用g[:,combinations[:,1]]生成的数组按列进行比较。因此，生成一个布尔数组diff。如果g有3列，则它看起来像这样，其中每列是列0,1，0,2和1,2的比较：
[[ True False False]
 [False  True False]
 [ True  True False]
 [False False False]
 [False  True False]
 [False False False]]

最后，对于每一列的值进行相加：
np.count_nonzero(diff,axis=0)
# Out
# [2 3 0]

此外，由于Python中的布尔类继承自整数类（大致上False==0和True==1，请参见此答案的"在Python中，False == 0 and True == 1是实现细节还是语言保证？"以获取更多信息）。np.count_nonzero为每个True位置添加1，这与使用np.sum获得的结果相同：
np.sum(diff,axis=0) 
# Out
# [2 3 0]

关于性能和内存的注意事项

对于大型数组，一次性处理整个数组可能需要过多的内存，导致出现内存错误，但是对于小型或中型数组，这往往是最快的方法。在某些情况下，按块处理可能会很有用：

combinations = np.array(list(itertools.combinations(range(g.shape[1]),2)))
n = len(combinations)
dist = np.empty(n)
# B = np.zeros((g.shape[1],g.shape[1]))
chunk = 200
for i in xrange(chunk,n,chunk):
    dist[i-chunk:i] = count_snp_diffs(g[:,combinations[i-chunk:i,0]], g[:,combinations[i-chunk:i,1]])
    # B[combinations[i-chunk:i,0],combinations[i-chunk:i,1]] = count_snp_diffs(g[:,combinations[i-chunk:i,0]], g[:,combinations[i-chunk:i,1]])
dist[i:] = count_snp_diffs(g[:,combinations[i:,0]], g[:,combinations[i:,1]])
# B[combinations[i:,0],combinations[i:,1]] = count_snp_diffs(g[:,combinations[i:,0]], g[:,combinations[i:,1]])

对于g.shape=(300,N)，在我的电脑上使用python 2.7、numpy 1.14.2和allel 1.1.10进行测试，%%timeit报告的执行时间如下：
  • 10列
    • numpy + 矩阵存储：107微秒
    • numpy + 1D存储：101微秒
    • allel：247微秒
  • 100列
    • numpy + 矩阵存储：15.7毫秒
    • numpy + 1D存储：16毫秒
    • allel：22.6毫秒
  • 1000列
    • numpy + 矩阵存储：1.54秒
    • numpy + 1D存储：1.53秒
    • allel：2.28秒
对于这些数组维度，纯numpy比allel模块略快，但应该针对您问题中的维度检查计算时间。
0

您可以使用np.dstack创建您期望的对，并使用np.apply_along_axis在第三个轴上应用该函数。

new = np.dstack((arr[:,:-1], arr[:, 1:]))
np.apply_along_axis(np.sum, 2, new)

例子：

In [86]: arr = np.array([[ 1,  1,  1,  1,  1,  1,  1,  1,  1, -1],
    ...:        [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1],
    ...:        [ 1,  1,  1,  1,  1,  1, -1,  1,  1,  1],
    ...:        [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1],
    ...:        [ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
    ...:        [ 1,  1,  1,  1,  1,  1,  1,  1,  1, -1],
    ...:        [-1, -1,  0, -1, -1, -1, -1, -1, -1,  0],
    ...:        [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1],
    ...:        [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1],
    ...:        [ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1]], dtype=np.int8)
    ...:        
    ...:        

In [87]: new = np.dstack((arr[:,:-1], arr[:, 1:]))

In [88]: new
Out[88]: 
array([[[ 1,  1],
        [ 1,  1],
        [ 1,  1],
        [ 1,  1],
        [ 1,  1],
        [ 1,  1],
        [ 1,  1],
        [ 1,  1],
        [ 1, -1]],

    ...

In [89]: 

In [89]: np.apply_along_axis(np.sum, 2, new)
Out[89]: 
array([[ 2,  2,  2,  2,  2,  2,  2,  2,  0],
       [ 2,  2,  2,  2,  2,  2,  2,  2,  2],
       [ 2,  2,  2,  2,  2,  0,  0,  2,  2],
       [ 2,  2,  2,  2,  2,  2,  2,  2,  2],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 2,  2,  2,  2,  2,  2,  2,  2,  0],
       [-2, -1, -1, -2, -2, -2, -2, -2, -1],
       [ 2,  2,  2,  2,  2,  2,  2,  2,  2],
       [ 2,  2,  2,  2,  2,  2,  2,  2,  2],
       [ 2,  2,  2,  2,  2,  2,  2,  2,  2]])
0
哦！现在我完全理解所需的输出了。我会编辑我的答案。从先验上看，我预计numpy会稍微快一点，但我不知道这对你来说是否很重要。 - OriolAbril
我编辑了答案。对于每个方法的计算时间进行了一些分析，但是你应该用你的 g 数组去验证，如果它真的很大，你可能会遇到 内存错误 或者使用 numpy 会变得更慢。 - OriolAbril
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接