二维numpy.take快速吗？

Question

numpy.take可应用于2维

np.take(np.take(T,ix,axis=0), iy,axis=1 )

我测试了离散二维拉普拉斯的模板

ΔT = T[ix-1,iy] + T[ix+1, iy] + T[ix,iy-1] + T[ix,iy+1] - 4 * T[ix,iy]

采用2种捕获方案和常用的numpy.array方案.文中引入了p和q两种函数，实现了一种精简的代码编写，并按不同的顺序对轴0和1进行了进位。这是代码：

nx = 300; ny= 300
T  = np.arange(nx*ny).reshape(nx, ny)
ix = np.linspace(1,nx-2,nx-2,dtype=int) 
iy = np.linspace(1,ny-2,ny-2,dtype=int)
#------------------------------------------------------------
def p(Φ,kx,ky):
    return np.take(np.take(Φ,ky,axis=1), kx,axis=0 )
#------------------------------------------------------------
def q(Φ,kx,ky):
    return np.take(np.take(Φ,kx,axis=0), ky,axis=1 )
#------------------------------------------------------------
%timeit ΔT_n = T[0:nx-2,1:ny-1] + T[2:nx,1:ny-1] + T[1:nx-1,0:ny-2]  + T[1:nx-1,2:ny] - 4.0 * T[1:nx-1,1:ny-1] 
%timeit ΔT_t = p(T,ix-1,iy)  + p(T,ix+1,iy)  + p(T,ix,iy-1)  + p(T,ix,iy+1)  - 4.0 * p(T,ix,iy)
%timeit ΔT_t = q(T,ix-1,iy)  + q(T,ix+1,iy)  + q(T,ix,iy-1)  + q(T,ix,iy+1)  - 4.0 * q(T,ix,iy)
.
1000 loops, best of 3: 944 µs per loop
100 loops, best of 3: 3.11 ms per loop
100 loops, best of 3: 2.02 ms per loop

结果似乎是显而易见的：

1. 通常的numpy指数算术数最快
2. 取方案q花费100%的时间(= C-订购?)
3. 采取方案p需要200%的时间(= Fortran-订购?)

甚至连1-dimensional 枕骨手册的例子都没有指出numpy.take是快速的：

a = np.array([4, 3, 5, 7, 6, 8])
indices = [0, 1, 4]
%timeit np.take(a, indices)
%timeit a[indices]
.
The slowest run took 6.58 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 4.32 µs per loop
The slowest run took 7.34 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 3.87 µs per loop

有没有人有过如何快速制作numpy.take的经验？这将是一种灵活而有吸引力的精益代码编写方法，它在编码和

被告知要迅速执行也是。感谢您的一些提示，以改善我的做法！

Answer 1

索引版本可能会使用这样的片对象来清理：

T[0:nx-2,1:ny-1] + T[2:nx,1:ny-1] + T[1:nx-1,0:ny-2]  + T[1:nx-1,2:ny] - 4.0 * T[1:nx-1,1:ny-1]

sy1 = slice(1,ny-1)
sx1 = slice(1,nx-1)
sy2 = slice(2,ny)
sy_2 = slice(0,ny-2)
T[0:nx-2,sy1] + T[2:nx,sy1] + T[sx1,xy_2]  + T[sx1,sy2] - 4.0 * T[sx1,sy1]

Answer 2

谢谢@Divakar和hpaulj！是的，使用slice也是可行的。对所有4种方法进行比较得出：

1. 最快前aequo: t(usual np)和t(slice)
2. t(take) =2* t(slice)
3. t(ix_) =3* t(slice)

在这里，代码和结果：

import numpy as np
from numpy import ix_ as r
nx = 500;    ny = 500
T = np.arange(nx*ny).reshape(nx, ny)

ix = np.arange(1,nx-1); 
iy = np.arange(1,ny-1);

jx = slice(1,nx-1); jxm = slice(0,nx-2); jxp = slice(2,nx)
jy = slice(1,ny-1); jym = slice(0,ny-2); jyp = slice(2,ny)

#------------------------------------------------------------
def p(U,kx,ky):
    return np.take(np.take(U,kx, axis=0), ky,axis=1)
#------------------------------------------------------------

%timeit ΔT_slice= -T[jxm,jy]     + T[jxp,jy]     - T[jx,jym]     + T[jx,jyp]     - 0.0 * T[jx,jy]
%timeit ΔT_npy  = -T[0:nx-2,1:ny-1] + T[2:nx,1:ny-1] - T[1:nx-1,0:ny-2]  + T[1:nx-1,2:ny] - 0.0 * T[1:nx-1,1:ny-1]
%timeit ΔT_take = -p(T,ix-1,iy)  + p(T,ix+1,iy)  - p(T,ix,iy-1)  + p(T,ix,iy+1)  - 0.0 * p(T,ix,iy)
%timeit ΔT_ix_  = -T[r(ix-1,iy)] + T[r(ix+1,iy)] - T[r(ix,iy-1)] + T[r(ix,iy+1)] - 0.0 * T[r(ix,iy)]
.
100 loops, best of 3: 3.14 ms per loop
100 loops, best of 3: 3.13 ms per loop
100 loops, best of 3: 7.03 ms per loop
100 loops, best of 3: 9.58 ms per loop

关于对视图和复制的讨论，以下几点可能具有指导意义：

print("if False --> a view ;   if True --> a copy"  )
print("_slice_ :", T[jx,jy].base is None)
print("_npy_   :", T[1:nx-1,1:ny-1].base is None)
print("_take_  :", p(T,ix,iy).base is None)
print("_ix_    :", T[r(ix,iy)].base is None)
.
if False --> a view ;   if True --> a copy
_slice_ : False
_npy_   : False
_take_  : True
_ix_    : True