Numpy:向量化一个集成二维数组的函数
Numpy: vectorizing a function that integrate 2D array
我需要对二维数组执行以下积分:
也就是说,网格中的每个点都得到值 RC,它是整个场与场值 U 在特定点 (x,y) 之间的差的 2D 积分,乘以归一化内核,在一维版本中为:
到目前为止我所做的是对索引的低效迭代:
def normalized_bimodal_kernel_2D(x,y,a,x0=0.0,y0=0.0):
""" Gives a kernel that is zero in x=0, and its integral from -infty to
+infty is 1.0. The parameter a is a length scale where the peaks of the
function are."""
dist = (x-x0)**2 + (y-y0)**2
return (dist*np.exp(-(dist/a)))/(np.pi*a**2)
def RC_2D(U,a,dx):
nx,ny=U.shape
x,y = np.meshgrid(np.arange(0,nx, dx),np.arange(0,ny,dx), sparse=True)
UB = np.zeros_like(U)
for i in xrange(0,nx):
for j in xrange(0,ny):
field=(U-U[i,j])*normalized_bimodal_kernel_2D(x,y,a,x0=i*dx,y0=j*dx)
UB[i,j]=np.sum(field)*dx**2
return UB
def centerlizing_2D(U,a,dx):
nx,ny=U.shape
x,y = np.meshgrid(np.arange(0,nx, dx),np.arange(0,ny,dx), sparse=True)
UB = np.zeros((nx,ny,nx,ny))
for i in xrange(0,nx):
for j in xrange(0,ny):
UB[i,j]=normalized_bimodal_kernel_2D(x,y,a,x0=i*dx,y0=j*dx)
return UB
您可以在此处查看 centeralizing
函数的结果:
U=np.eye(20)
plt.imshow(centerlizing(U,10,1)[10,10])
我确定我还有其他错误,因此热烈欢迎任何反馈,但我真正感兴趣的是了解如何以矢量化方式更快地完成此操作。
假设 dx=1
,因为我不确定您要用该离散化做什么:
def normalized_bimodal_kernel_2D(x, y, a):
#generating a 4-d tensor instead of 1d vector
dist = (x[:,None,None,None] - x[None,None,:,None])**2 +\
(y[None,:,None,None] - y[None,None,None,:])**2
return (dist * np.exp(-(dist / a))) / (np.pi * a**2)
def RC_2D(U, a):
nx, ny = U.shape
x, y = np.arange(nx), np.arange(ny)
U4 = U[:, :, None, None] - U[None, None, :, :] #Another 4d
k = normalized_bimodal_kernel_2D(x, y, a)
return np.einsum('ijkl,ijkl->ij', U4, k)
def centerlizing_2D(U, a):
nx, ny = U.shape
x, y = np.arange(nx), np.arange(ny)
return normalized_bimodal_kernel_2D(x, y, a)
基本上,numpy
中的向量化 for
循环是添加更多维度的问题。您在 2D U
向量上进行了两次循环,因此要对其进行向量化,只需将其转换为 4D。
normalized_bimodal_kernel_2D
在两个嵌套循环中被调用,每个循环只偏移一小步。这会重复很多计算。
centerlizing_2D
的优化是为更大的范围计算一次内核,然后定义 UB
以将转移的观点纳入其中。这可以使用 stride_tricks
,不幸的是它是相当高级的 numpy。
def centerlizing_2D_opt(U,a,dx):
nx,ny=U.shape
x,y = np.meshgrid(np.arange(-nx//2, nx+nx//2, dx),
np.arange(-nx//2, ny+ny//2, dx), # note the increased range
sparse=True)
k = normalized_bimodal_kernel_2D(x, y, a, x0=nx//2, y0=ny//2)
sx, sy = k.strides
UB = as_strided(k, shape=(nx, ny, nx*2, ny*2), strides=(sy, sx, sx, sy))
return UB[:, :, nx:0:-1, ny:0:-1]
assert np.allclose(centerlizing_2D(U,10,1), centerlizing_2D_opt(U,10,1)) # verify it's correct
是的,速度更快:
%timeit centerlizing_2D(U,10,1) # 100 loops, best of 3: 9.88 ms per loop
%timeit centerlizing_2D_opt(U,10,1) # 10000 loops, best of 3: 85.9 µs per loop
接下来我们优化RC_2D
,用优化后的centerlizing_2D
例程表示:
def RC_2D_opt(U,a,dx):
UB_tmp = centerlizing_2D_opt(U, a, dx)
U_tmp = U[:, :, None, None] - U[None, None, :, :]
UB = np.sum(U_tmp * UB_tmp, axis=(0, 1))
return UB
assert np.allclose(RC_2D(U,10,1), RC_2D_opt(U,10,1))
%timeit RC_2D(U,10, 1)
的表现:
#original: 100 loops, best of 3: 13.8 ms per loop
#@DanielF's: 100 loops, best of 3: 6.98 ms per loop
#mine: 1000 loops, best of 3: 1.83 ms per loop
为了适合您的公式,让 U
成为一个函数。
然后你只需要将 x,y,x',y'
和 np.ix_
放在四个不同的维度上,然后仔细翻译你的公式。 Numpy 广播将完成剩下的工作。
a=20
x,y,xp,yp=np.ix_(*[np.linspace(0,1,a)]*4)
def U(x,y) : return np.float32(x == y) # function "eye"
def f(x,y,xp,yp,a):
r2=(x-xp)**2+(y-yp)**2
return r2*np.exp(-r2/a)*(U(xp,yp) - U(x,y))/np.pi/a/a
#f(x,y,xp,yp,a).shape is (20, 20, 20, 20)
RC=f(x,y,xp,yp,a).sum(axis=(2,3))
#RC.shape is (20, 20)
我需要对二维数组执行以下积分:
到目前为止我所做的是对索引的低效迭代:
def normalized_bimodal_kernel_2D(x,y,a,x0=0.0,y0=0.0):
""" Gives a kernel that is zero in x=0, and its integral from -infty to
+infty is 1.0. The parameter a is a length scale where the peaks of the
function are."""
dist = (x-x0)**2 + (y-y0)**2
return (dist*np.exp(-(dist/a)))/(np.pi*a**2)
def RC_2D(U,a,dx):
nx,ny=U.shape
x,y = np.meshgrid(np.arange(0,nx, dx),np.arange(0,ny,dx), sparse=True)
UB = np.zeros_like(U)
for i in xrange(0,nx):
for j in xrange(0,ny):
field=(U-U[i,j])*normalized_bimodal_kernel_2D(x,y,a,x0=i*dx,y0=j*dx)
UB[i,j]=np.sum(field)*dx**2
return UB
def centerlizing_2D(U,a,dx):
nx,ny=U.shape
x,y = np.meshgrid(np.arange(0,nx, dx),np.arange(0,ny,dx), sparse=True)
UB = np.zeros((nx,ny,nx,ny))
for i in xrange(0,nx):
for j in xrange(0,ny):
UB[i,j]=normalized_bimodal_kernel_2D(x,y,a,x0=i*dx,y0=j*dx)
return UB
您可以在此处查看 centeralizing
函数的结果:
U=np.eye(20)
plt.imshow(centerlizing(U,10,1)[10,10])
假设 dx=1
,因为我不确定您要用该离散化做什么:
def normalized_bimodal_kernel_2D(x, y, a):
#generating a 4-d tensor instead of 1d vector
dist = (x[:,None,None,None] - x[None,None,:,None])**2 +\
(y[None,:,None,None] - y[None,None,None,:])**2
return (dist * np.exp(-(dist / a))) / (np.pi * a**2)
def RC_2D(U, a):
nx, ny = U.shape
x, y = np.arange(nx), np.arange(ny)
U4 = U[:, :, None, None] - U[None, None, :, :] #Another 4d
k = normalized_bimodal_kernel_2D(x, y, a)
return np.einsum('ijkl,ijkl->ij', U4, k)
def centerlizing_2D(U, a):
nx, ny = U.shape
x, y = np.arange(nx), np.arange(ny)
return normalized_bimodal_kernel_2D(x, y, a)
基本上,numpy
中的向量化 for
循环是添加更多维度的问题。您在 2D U
向量上进行了两次循环,因此要对其进行向量化,只需将其转换为 4D。
normalized_bimodal_kernel_2D
在两个嵌套循环中被调用,每个循环只偏移一小步。这会重复很多计算。
centerlizing_2D
的优化是为更大的范围计算一次内核,然后定义 UB
以将转移的观点纳入其中。这可以使用 stride_tricks
,不幸的是它是相当高级的 numpy。
def centerlizing_2D_opt(U,a,dx):
nx,ny=U.shape
x,y = np.meshgrid(np.arange(-nx//2, nx+nx//2, dx),
np.arange(-nx//2, ny+ny//2, dx), # note the increased range
sparse=True)
k = normalized_bimodal_kernel_2D(x, y, a, x0=nx//2, y0=ny//2)
sx, sy = k.strides
UB = as_strided(k, shape=(nx, ny, nx*2, ny*2), strides=(sy, sx, sx, sy))
return UB[:, :, nx:0:-1, ny:0:-1]
assert np.allclose(centerlizing_2D(U,10,1), centerlizing_2D_opt(U,10,1)) # verify it's correct
是的,速度更快:
%timeit centerlizing_2D(U,10,1) # 100 loops, best of 3: 9.88 ms per loop
%timeit centerlizing_2D_opt(U,10,1) # 10000 loops, best of 3: 85.9 µs per loop
接下来我们优化RC_2D
,用优化后的centerlizing_2D
例程表示:
def RC_2D_opt(U,a,dx):
UB_tmp = centerlizing_2D_opt(U, a, dx)
U_tmp = U[:, :, None, None] - U[None, None, :, :]
UB = np.sum(U_tmp * UB_tmp, axis=(0, 1))
return UB
assert np.allclose(RC_2D(U,10,1), RC_2D_opt(U,10,1))
%timeit RC_2D(U,10, 1)
的表现:
#original: 100 loops, best of 3: 13.8 ms per loop
#@DanielF's: 100 loops, best of 3: 6.98 ms per loop
#mine: 1000 loops, best of 3: 1.83 ms per loop
为了适合您的公式,让 U
成为一个函数。
然后你只需要将 x,y,x',y'
和 np.ix_
放在四个不同的维度上,然后仔细翻译你的公式。 Numpy 广播将完成剩下的工作。
a=20
x,y,xp,yp=np.ix_(*[np.linspace(0,1,a)]*4)
def U(x,y) : return np.float32(x == y) # function "eye"
def f(x,y,xp,yp,a):
r2=(x-xp)**2+(y-yp)**2
return r2*np.exp(-r2/a)*(U(xp,yp) - U(x,y))/np.pi/a/a
#f(x,y,xp,yp,a).shape is (20, 20, 20, 20)
RC=f(x,y,xp,yp,a).sum(axis=(2,3))
#RC.shape is (20, 20)