我是Python / Numpy的新手，正在尝试将我的三重循环改进为更有效的计算，但无法静静地弄清楚该怎么做。
计算是在大小为（25,35）的网格上进行的，数组的形状为：
A = (8760,25,35)
B = (12,25,35)
A中的第一维度对应于一年中的小时数（〜8760），B中的第一维度对应于月数（12）。我想在第一个月使用B [0，：，：]中的值，在第二个月使用B [1，:::]等。
到目前为止，我以一种非常未经改进的方式创建了一个索引数组，该索引数组填充有1,1,1 ...，2,2,2 ...，12，以从B中提取值。我的代码带有一些随机数
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支持broadcasting允许以高度优化的方式跨不同形状的数组执行元素操作。你的情况，你有行数和列数A和B相同的。但是，在第一维上，这两个数组的元素数量不同。看一下实现，似乎B每个q数字都重复了第一个维度元素，直到我们转到第一个维度中的下一个元素为止。这与一个事实，即在第一维中的元素的数量相一致B是q倍的第一尺寸元件的数量A。
现在，回到broadcasting，解决方案是将的第一个维度拆分为A一个4D数组，这样我们可以使此重塑后的4D数组的第一个维度中的元素数量与B的第一个维度中的元素数量匹配。接下来，reshape B通过使用，在第二维上创建单例维（无元素的维），从而生成4D数组B[:,None,:,:]。然后，NumPy将使用广播魔术并执行广播的元素乘法，因为这就是我们在中所做的some_function。
这是使用NumPy's broadcasting功能的矢量化实现-
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)
运行时测试和输出验证-
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
本文收集自互联网，转载请注明来源。
如有侵权，请联系 [email protected] 删除。
我来说两句