3

我试图通过对二维数组应用操作来提高 numpy 性能，问题是数组中每个元素的值取决于该元素的 i,j 位置。

显然，最简单的方法是使用嵌套的 for 循环，但我想知道是否有更好的方法通过引用 np.indices 或类似的东西？这是我的“愚蠢”代码：

for J in range(1025):
    for I in range(1025):
        PSI[I][J] = A*math.sin((float(I+1)-.5)*DI)*math.sin((float(J+1)-.5)*DJ)
        P[I][J] = PCF*(math.cos(2.*float(I)*DI)+math.cos(2.*float(J)*DJ))+50000.
4

3 回答 3

7

由于您在两个数组之间进行乘法运算，因此可以在使用后使用外部arange函数来获取 sin/cos 的数组。

像这样的东西（使用 numpy 的三角函数，因为它们是矢量化的）

PSI_i = numpy.sin((arange(1,1026)-0.5)*DI)
PSI_j = numpy.sin((arange(1,1026)-0.5)*DJ)
PSI = A*outer(PSI_i, PSI_j)

P_i = numpy.cos(2.*arange(1,1026)*DI)
P_j = numpy.cos(2.*arange(1,1026)*DJ)
P = PCF*outer(P_i, P_j) + 50000

如果您的环境是使用from numpy import *or设置的from pylab import *，那么在您的 trig 函数之前不需要这些numpy.前缀。我将它们保留在其中以将它们与其他方法区分开来math，这不适用于这种方法。

于 2012-05-30T17:55:02.650 回答
1

您可以使用索引获取索引值的网格：

I,J=np.indices(PSI.shape)
#All constants set to one
PSI2=np.sin(I+1-.5)*np.sin(J+1-.5)
print PSI-PSI2 # should be zero.

我用 ipython 做了一些计时：

import numpy as np
import math
A = 1
P = 1
DI = 1
DJ = 1

def a():
    PSI=np.zeros((1025,1025))
    for J in range(1025):
        for I in range(1025):
            PSI[I][J] = A*math.sin((float(I+1)-.5)*DI)*math.sin((float(J+1)-.5)*DJ)
%timeit a()

def b():
    PSI=np.zeros((1025,1025))
    for I,J in np.ndindex(*PSI.shape):
        PSI[I,J] = A*math.sin((float(I+1)-.5)*DI)*math.sin((float(J+1)-.5)*DJ)        
%timeit b()

def c():
    I,J=np.indices((1025, 1025))
    P2=A*np.sin((I+1-.5)*DI)*np.sin((J+1-.5)*DJ)    
%timeit c()

def d():
    PSI_i = np.sin((np.arange(1,1026)-0.5)*DI)
    PSI_j = np.sin((np.arange(1,1026)-0.5)*DJ)
    PSI = A*np.outer(PSI_i, PSI_j)    
%timeit d()

结果在我的机器上一点也不奇怪：

1 loops, best of 3: 1.75 s per loop
1 loops, best of 3: 3.51 s per loop
10 loops, best of 3: 77.1 ms per loop
100 loops, best of 3: 7.16 ms per loop
于 2012-05-30T17:55:38.977 回答
0

试试 numpy 的 ndenumerate 函数，它返回值和索引：

>>> a
array([[5, 5, 5],
       [1, 2, 3]])


>>> for index, value in numpy.ndenumerate(a):
...     print index, value

(0, 0) 5
(0, 1) 5
(0, 2) 5
(1, 0) 1
(1, 1) 2
(1, 2) 3
于 2012-05-30T17:52:27.250 回答