在具有不同数组形状的Python / Numpy中向量化Triple for循环

问题描述

我是Python / Numpy的新手，正在尝试将我的三重循环改进为更有效的计算，但是无法安静地找到方法。

I am new in Python/Numpy and is trying to improve my triple for loop into a more efficient calculation, but can't quiet figure out how to do it.

计算是在大小为（25,35）的网格上进行的，数组的形状为：

The calculations is carried out on a grid of the size (25,35) and the shapes of arrays is:

 A = (8760,25,35)
 B = (12,25,35)

A中的第一个维度对应于一年中的小时数（〜8760），B中的第一个维度对应于月数（12）。我想在第一个月使用B [0，：，：]中的值，在第二个月使用B [1，：，：]等。

The first dimensions in A corresponds to the number hours in a year (~8760), and the first dimension in B is the number of months(12). I want to use the values in B[0,:,:] for the first month, and B[1,:,:] for the second etc.

所以到目前为止，我以一种非常未经改进的方式创建了一个索引数组，该数组填充有1,1,1 ...，2,2,2 ...，12，以从B中提取值。我的代码带有一些随机数

So far I created, in a very unrefined way, a index array filled with 1,1,1...,2,2,2...,12 to extract the values from B. My code with some random numbers

    N,M = (25, 35)
    A = np.random.rand(8760,N,M)
    B = np.random.rand(12,N,M)
    q = len(A)/12
    index = np.hstack((np.full((1,q),1),np.full((1,q),2),np.full((1,q),3),np.full((1,q),4),np.full((1,q),5),np.full((1,q),6),np.full((1,q),7),np.full((1,q),8),np.full((1,q),9),np.full((1,q),10),np.full((1,q),11),np.full((1,q),12)))-1
    results = np.zeros((len(A),N,M))

    for t in xrange(len(A)):
        for i in xrange(N):
            for j in xrange(M):
                results[t][i][j] = some_function(A[t][i][j], B[index[(0,t)]][i][j],H = 80.)


    def some_function(A,B,H = 80.0):
        results = A*np.log(H/B)/np.log(10./B)
        return results

如何提高计算速度？

推荐答案

NumPy支持允许以高度优化的方式跨不同形状的数组执行元素操作。在您的情况下， A 和 B 中的行数和列数相同。但是，在第一维上，这两个数组的元素数量不同。从实现上看，似乎 B 的第一个维度元素是每个 q 编号重复的，直到我们转到第一个维度中的下一个元素。这与以下事实相吻合： B 的第一维元素数是 q 乘以first维的元素数 A 的维度。

NumPy suports broadcasting that allows elementwise operations to be performed across different shaped arrays in a highly optimized manner. In your case, you have the number of rows and columns in A and B the same. But, at the first dimension, the number of elements are different across these two arrays. Looking at the implementation, it seems B 's first dimension elements are repeated per q number until we go over to the next element in it's first dimension. This coincides with the fact that the number of elements in first dimension of B is q times the number of elements in first dimension of A.

现在，回到广播，解决方法是 split A 的第一个维度具有4D数组，这样我们在第一个维度中拥有元素数量此重塑后的4D数组的尺寸与B的第一维中的元素数量匹配。接下来，通过在第二个创建单例尺寸（无元素的尺寸），将 reshape B 转换为4D数组。 B [：，None，：，：] 的尺寸。然后，NumPy将使用广播魔术并执行广播的元素乘法，因为这就是我们在 some_function 中所做的工作。

Now, going back to broadcasting, the solution would be to split the first dimension of A to have a 4D array, such that we have the number of elements in the first dimension of this reshaped 4D array matching up with the number of elements in B's first dimension. Next up, reshape B to a 4D array as well by creating singleton dimension (dimension with no elements) at the second dimension with B[:,None,:,:]. Then, NumPy would use broadcasting magic and perform broadcasted elementwise multiplications, as that's what we are doing in our some_function.

这是使用能力-

Here's the vectorized implementation using NumPy's broadcasting capability -

H = 80.0
M,N,R = B.shape
B4D = B[:,None,:,:]
out = ((A.reshape(M,-1,N,R)*np.log(H/B4D))/np.log(10./B4D)).reshape(-1,N,R)

运行时测试和输出验证-

Runtime tests and output verification -

In [50]: N,M = (25, 35)
    ...: A = np.random.rand(8760,N,M)
    ...: B = np.random.rand(12,N,M)
    ...: H = 80.0
    ...: 

In [51]: def some_function(A,B,H = 80.0):
    ...:     return A*np.log(H/B)/np.log(10./B)
    ...: 

In [52]: def org_app(A,B,H):
    ...:    q = len(A)/len(B)
    ...:    index = np.repeat(np.arange(len(B))[:,None],q,axis=1).ravel()[None,:] # Simpler
    ...:    results = np.zeros((len(A),N,M))
    ...:    for t in xrange(len(A)):
    ...:        for i in xrange(N):
    ...:            for j in xrange(M):
    ...:                results[t][i][j] = some_function(A[t][i][j], B[index[(0,t)]][i][j])
    ...:    return results
    ...: 

In [53]: def vectorized_app(A,B,H):
    ...:    M,N,R = B.shape
    ...:    B4D = B[:,None,:,:]
    ...:    return ((A.reshape(M,-1,N,R)*np.log(H/B4D))/np.log(10./B4D)).reshape(-1,N,R)
    ...: 

In [54]: np.allclose(org_app(A,B,H), vectorized_app(A,B,H))
Out[54]: True

In [55]: %timeit org_app(A,B,H)
1 loops, best of 3: 1min 32s per loop

In [56]: %timeit vectorized_app(A,B,H)
10 loops, best of 3: 217 ms per loop

这篇关于在具有不同数组形状的Python / Numpy中向量化Triple for循环的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！