Python中成对相关的优化计算-Java 学习之路

给定一组离散位置（例如"sites"），它们以某些分类方式成对相关（例如一般接近度）并包含局部水平数据（例如种群大小），我希望 efficiently 计算成对位置上本地水平数据之间的平均相关系数具有相同关系的特征 .

例如，我假设100个站点并使用值1到25随机化它们的成对关系，产生三角矩阵 relations ：

import numpy as np

sites = 100
categ = 25

relations = np.random.randint(low=1, high=categ+1, size=(sites, sites))
relations = np.triu(relations) # set relation_ij = relation_ji
np.fill_diagonal(relations, 0) # ignore self-relation

我还在每个站点上复制了5000个模拟结果：

sims = 5000
res = np.round(np.random.rand(sites, sims),1)

为了计算每个特定关系类别的平均成对相关性，我首先计算每个关系类别 i 每个唯一站点对的模拟结果 res 之间的相关系数 res ，然后采用关系 i 的所有可能对的平均值：

rho_list = np.ones(categ)*99

for i in range(1, categ+1):
    idr = np.transpose(np.where(relations == i)) # pairwise site indices of the same relation category
    comp = np.vstack([res[idr[:,0]].ravel(), res[idr[:,1]].ravel()]) # pairwise comparisons of simulation results from the same relation category
    comp_uniq = np.reshape(comp.T, (len(idr), res.shape[1], -1)) # reshape above into pairwise comparisons of simulation results between unique site pairs

    rho = np.ones(len(idr))*99 # correlation coefficients of all unique site pairs of current relation category

    for j in range(len(idr)): # loop through unique site pairs
        comp_uniq_s = comp_uniq[j][np.all(comp_uniq!=0, axis=2)[j]].T # shorten comparisons by removing pairs with zero-valued result
        rho[j] = np.corrcoef(comp_uniq_s[0], comp_uniq_s[1])[0,1]

    rho_list[i-1] = np.nanmean(rho)

虽然这个脚本可以工作，但是一旦我增加了 sites = 400 ，那么整个计算可能需要6个多小时才能完成，这让我质疑我对数组函数的使用 . 这种糟糕表现的原因是什么？我该如何优化算法？

1 回答

我们可以使用 j 迭代器使用某些 masking 对最内层循环进行向量化，以处理在该循环的每次迭代中处理的数据的不规则性质 . 我们也可以放慢速度 np.corrcoef （灵感来自this post） . 此外，我们可以在外循环开始时优化几个步骤，特别是堆叠步骤，这可能是瓶颈 .

因此，完整的代码将减少到这样的东西 -

for i in range(1, categ+1):
    r,c = np.where(relations==i)

    A = res[r]
    B = res[c]

    mask0 = ~((A!=0) & (B!=0))
    A[mask0] = 0
    B[mask0] = 0

    count = mask0.shape[-1] - mask0.sum(-1,keepdims=1)
    A_mA = A - A.sum(-1, keepdims=1)/count
    B_mB = B - B.sum(-1, keepdims=1)/count

    A_mA[mask0] = 0
    B_mB[mask0] = 0

    ssA = np.einsum('ij,ij->i',A_mA, A_mA)
    ssB = np.einsum('ij,ij->i',B_mB, B_mB)
    rho = np.einsum('ij,ij->i',A_mA, B_mB)/np.sqrt(ssA*ssB)

    rho_list[i-1] = np.nanmean(rho)

Runtime tests

案例＃1：对于给定的样本数据，网站= 100

In [381]: %timeit loopy_app()
1 loop, best of 3: 7.45 s per loop

In [382]: %timeit vectorized_app()
1 loop, best of 3: 479 ms per loop

15x+ 加速 .

案例＃2： sites = 200

In [387]: %timeit loopy_app()
1 loop, best of 3: 1min 56s per loop

In [388]: %timeit vectorized_app()
1 loop, best of 3: 1.86 s per loop

In [390]: 116/1.86
Out[390]: 62.36559139784946

62x+ 加速 .

案例＃3：最后用 sites = 400

In [392]: %timeit vectorized_app()
1 loop, best of 3: 7.64 s per loop

这在OP结束时以循环方式获得了 6hrs+ .

从时间上看，很明显，对内部循环进行矢量化是获得大型 sites 显着加速的关键 .

回复于 2024-05-03T11:04:38+08:00

Python中成对相关的优化计算

1 回答

相关问题