有没有办法用更好的 "in place" 方法来做到这一点？

Question

这是 Biot-Savart law 的简单近似值。

我已经在函数 calc()、

中实现了积分（和）

如果空间点的数量很大，比如 10^7 或 10^8 -ish，是否可以编写 calc 以更有效地使用 NumPy 数组？感谢您的建议！

def calc(points, x_seg, idl_seg):

    r = points[:, None, :] - x_seg[None, :, :]             # START CALCULATION

    bottom = ((r**2).sum(axis=-1)**1.5)[...,None]     # 1/|r|**3 add axis for vector

    top = np.cross(idl_seg[None,:,:], r)                  # np.cross defaults to last axis

    db = (mu0 / four_pi) * top / bottom

    b = db.sum(axis=-2)               # sum over the segments of the current loop

    return b

编辑：例如，我可以做到这一点。现在只有两个大小为 nx * ny * nz * nseg * 3 的数组（r 和 hold）。也许我应该一次传递更小的 points 块，这样它就可以一次放入缓存？

def calc_alt(points, x_seg, idl_seg):

    r = points[:, None, :] - x_seg[None, :, :]             

    hold = np.ones_like(r)*((r**2).sum(axis=-1)**-1.5)[...,None]  # note **-1.5 neg

    b = (hold * np.cross(idl_seg[None,:,:], r)).sum(axis=-2)

    return b * (mu0 / four_pi)

发布其余代码以显示如何使用 calc。

import numpy as np
import matplotlib.pyplot as plt

pi, four_pi  = np.pi,  4. * np.pi
mu0          = four_pi * 1E-07 # Tesla m/A exact, defined
r0           = 0.05   # meters
I0           = 100.0  # amps
nx, ny, nz = 48, 49, 50

x,y,z = np.linspace(0,2*r0,nx), np.linspace(0,2*r0,ny), np.linspace(0,2*r0,nz) 
xg = np.zeros((nx, ny, nz, 3))  # 3D grid of position vectors
xg[...,0] = x[:, None, None]   # fill up the positions
xg[...,1] = y[None, :, None]
xg[...,2] = z[None, None, :]
xgv = xg.reshape(nx*ny*nz, 3)  # flattened view of spatial points

nseg = 32   # approximate the current loop as a set of discrete points I*dl 
theta = np.linspace(0, 2.*pi, nseg+1)[:-1]  # get rid of the repeat

xdl = np.zeros((nseg, 3))   # these are the position vectors
idl = np.zeros((nseg, 3))   # these are the current vectors

xdl[:,0],  xdl[:,1] = r0 * np.cos(theta),   r0 * np.sin(theta)
idl[:,0],  idl[:,1] = I0 * -np.sin(theta),  I0 * np.cos(theta)

b = calc(xgv, xdl, idl)           # HERE IS THE CALCULATION

bv = b.reshape(nx, ny, nz, 3)     # make a "3D view" again to use for plotting

bx, by, bz = bv[...,0], bv[...,1], bv[...,2]  # make component views

bperp = np.sqrt(bx**2 + by**2)  # new array for perp field

zround = np.round(z, 4)
iz = 5     # choose a transverse plane for a plot
fields    = [ bz,   bperp,   bx,   by]
names     = ['Bz', 'Bperp', 'Bx', 'By']
titles = ["approx " + name + " at z = " + str(zround[iz])
          for name in names]

plt.figure()
for i, field in enumerate(fields):
    print i
    plt.subplot(2, 2, i+1)
    plt.imshow(field[..., iz], origin='lower')  # fields at iz don't use Jet !!!
    plt.title(titles[i])
    plt.colorbar()
plt.show()

最后的绘图只是为了看看它是否有效。实际上，永远不要使用默认的颜色图。坏的，可怕的，调皮的杰特！在这种情况下，具有对称 vmin = -vmax 的发散 cmap 可能会很好。（参见 Jake VanderPlas 的 post, and the matplotlib documentation, and there's some lovely demos down here.

Answer 1

您可以压缩这些行：

b = db.sum(axis=-2)               # sum over the segments of the current loop

bv = b.reshape(nx, ny, nz, 3)     # make a "3D view" again to use for plotting

bx, by, bz = bv[...,0], bv[...,1], bv[...,2]

进入

bx, by, bz = np.split(db.sum(axis=-2).reshape(nx, ny, nz, 3), 3, -1)

我怀疑它是否对速度有任何影响。它是否使这一点更清楚是有争议的。

xdl = np.zeros((nseg, 3))   # these are the position vectors
idl = np.zeros((nseg, 3))   # these are the current vectors

xdl[:,0],  xdl[:,1] = r0 * np.cos(theta),   r0 * np.sin(theta)
idl[:,0],  idl[:,1] = I0 * -np.sin(theta),  I0 * np.cos(theta)

可以重写为（未测试）

xdl = r0 * np.array([np.cos(theta), np.sin(theta)]
idl = I0 * np.array([-np.sin(theta), np.cos(theta)]

虽然这些会使这些 (3,nseg)。请注意，split 的默认轴是 0。在第一个轴上组合和拆分通常更自然。另外 [None,...] 广播是自动的。

ng 构造也可能得到简化。

这些主要是外观上的变化，不会对性能产生大的影响。

Answer 2

我在 np.numexpr 中只有运行执行我在编辑中建议的（除其他外）- 将数组分成“块”，以便它们可以放入缓存，包括所有计算表达式所需的临时数组。

https://numexpr.readthedocs.io/projects/NumExpr3/en/latest/user_guide.html

有很好的解释here and especially in this wiki。

有没有办法用更好的 "in place" 方法来做到这一点？

Is there a way to do this with better "in place" methods?

python-2.7

numpy

physics