Python numpy：对 numpy 二维数组中的每对列执行函数？

我试图将一个函数应用于 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])

任何想法都会很棒！谢谢你！

您可以简单地将函数定义为：

def count_snp_diffs(x, y): 
    return np.count_nonzero((x != y) & (x >= 0) & (y >= 0),axis=0)

然后使用生成的数组作为索引来调用它itertools.combinations， 为了得到all可能的列组合：

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

另外，如果输出must存储在矩阵中（对于大的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]])

Now, d[i,j]返回列之间的距离i and 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]

这是对正在发生的事情的一步一步的解释。

Calling g[:,combinations[:,0]]访问排列第一列中的所有列，生成一个新数组，将其与生成的数组逐列进行比较g[:,combinations[:,1]]。因此，一个布尔数组diff被生成。如果g有 3 列，看起来像这样，其中每一列都是列的比较0,1, 0,2 and 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中的boolean类继承自integer类（大致False==0和和True==1，看这个answer of “Python 中的 False == 0 和 True == 1 是实现细节还是由语言保证的？”了解更多信息）。这np.count_nonzero每项加 1True位置，这与获得的结果相同np.sum:

np.sum(diff,axis=0) 
# Out
# [2 3 0]

关于性能和内存的注释

对于大型数组，一次处理整个数组可能需要太多内存，您可以获得Memory Error然而，对于小型或中型阵列，它往往是最快的方法。在某些情况下，按块工作可能很有用：

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]])

For g.shape=(300,N)，报告的执行时间%%timeit在我的电脑上，python 2.7、numpy 1.14.2 和 allel 1.1.10 是：

• 10 columns
  • numpy + 矩阵存储：107 µs
  • numpy + 一维存储：101 µs
  • 等位基因：247 µs
• 100 columns
  • numpy + 矩阵存储：15.7 毫秒
  • numpy + 一维存储：16 毫秒
  • 等位基因：22.6 毫秒
• 1000 columns
  • numpy + 矩阵存储：1.54 s
  • numpy + 一维存储：1.53 秒
  • 等位基因：2.28 s

使用这些数组维度，纯 numpy 比 allen 模块快一点，但应检查问题中的维度的计算时间。