在球壳内生成积分点

Question

给定一个球壳体积（即一个球体减去一个相同中心的较小球体），我需要计算该体积内的所有积分位置。

我可以通过计算外半径内的所有积分点（通过直到外径的所有整数的笛卡尔积）轻松地做到这一点，然后删除内半径内的点。

但是，在某些情况下，半径可能相当大，因此生成外半径内的所有点可能比我想要的更耗时。因此，如果可能的话，我想消除计算内半径内的点所浪费的处理时间（或至少减少必要的数量）。

这是一些当前的代码，但是半径可以更大。

import numpy as n
import numpy.linalg as la

def cartesian_product(*arrays):
    la = len(arrays)
    dtype = n.result_type(*arrays)
    arr = n.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(n.ix_(*arrays)):
        arr[..., i] = a
    return arr.reshape(-1, la)

def points_in_shell(inner_radius, outer_radius):
    cutoff = int(n.ceil(outer_radius))
    int_range = n.arange(-cutoff, cutoff+1)
    pos_initial = cartesian_product(int_range, int_range, int_range)
    pos_norm = la.norm(pos_initial, axis=1)

    valid_idx = n.where((pos_norm >= inner_radius) & (pos_norm <= outer_radius))
    pos_final = pos_initial[valid_idx]
    return pos_final

vals = points_in_shell(89.634, 95.254)

有什么想法吗？

Answer 1

更新：适应更新的问题，也使解决方案完全稀疏。

def f_sparser_oct(IR, OR):
    OR2, IR2 = OR*OR, IR*IR
    out = np.ones((int(OR)+1, 3), dtype=int)
    out[:, 0] = i = np.arange(int(OR)+1)
    i2 = i*i
    jh = np.sqrt(OR2-i2).astype(int)
    out = np.repeat(out, jh+1, axis=0)
    I, J, K = out.T
    J[0] = 0
    J[np.cumsum(jh[:-1]+1)] = -jh[:-1]
    J[:] = np.cumsum(J)
    kh = np.sqrt(OR2-i2[I]-i2[J]).astype(int)
    kl = np.ceil(np.sqrt(np.maximum(0, IR2-i2[I]-i2[J]))).astype(int)
    out = np.repeat(out, 1+kh-kl, axis=0)
    I, J, K = out.T
    K[0] = kl[0]
    K[np.cumsum(1 + kh[:-1] - kl[:-1])] = kl[1:] - kh[:-1]
    K[:] = np.cumsum(K)
    return out

def f_sparser_full(IR, OR):
    OR2, IR2 = OR*OR, IR*IR
    out = np.ones((2*int(OR)+1, 3), dtype=int)
    out[:, 0] = i = np.arange(-int(OR), int(OR)+1)
    i2 = i*i
    i2i = np.r_[i2[int(OR):], i2[:int(OR)]]
    jh = np.sqrt(OR2-i2).astype(int)
    out = np.repeat(out, 2*jh+1, axis=0)
    I, J, K = out.T
    J[0] = -jh[0]
    J[np.cumsum(2*jh[:-1]+1)] = -jh[1:] - jh[:-1]
    J[:] = np.cumsum(J)
    kh = np.sqrt(OR2-i2i[I]-i2i[J]).astype(int)
    kl = np.ceil(np.sqrt(np.maximum(0, IR2-i2i[I]-i2i[J]))).astype(int)
    szs = 2*(1 + kh - kl)-(kl==0)
    out = np.repeat(out, szs, axis=0)
    I, J, K = out.T
    K[0] = -kh[0]
    SZS = np.cumsum(szs)
    K[SZS[:-1]] = -kh[1:] - kh[:-1]
    K[SZS[kl!=0] - 1 - kh[kl!=0] + kl[kl!=0]] = 2 * kl[kl!=0]
    K[:] = np.cumsum(K)
    return out

演示：

>>> np.all(f_sparser_full(89.634, 95.254) == vals)
True
>>> t0 = time.perf_counter(); f_sparser_full(998, 1000).shape; t1 = time.perf_counter()
(25092914, 3)
>>> t1-t0
0.969647804973647

更新结束。

这里有两个快速的numpy 建议：第一个很短但占用大量内存。第二个有一个循环但使用较少的内存：

import numpy as np

def f_where(IR, OR):
    i, j, k = np.ogrid[:OR+1, :OR+1, :OR+1]
    R2 = i*i + j*j + k*k
    return np.argwhere((R2>=IR*IR) & (R2<=OR*OR))

def f_sparse(IR, OR):
    out = []
    for I in range(int(OR)+1):
        O = np.sqrt(OR*OR - I*I)
        i, j, k = np.ogrid[I:I+1, :O+1, :O+1]
        R2 = i*i + j*j + k*k
        out.append(np.argwhere((R2>=IR*IR) & (R2<=OR*OR))+(I, 0, 0))
    return np.concatenate(out, axis=0)

演示：

>>> np.all(f_sparse(89.634, 95.254) == vals)
True
>>> 
>>> import time
>>> 
>>> t0 = time.perf_counter(); f_where(298.5, 300.0).shape; t1 = time.perf_counter()
(211366, 3)
>>> t1-t0
0.20139808906242251
>>> t0 = time.perf_counter(); f_sparse(298.5, 300.0).shape; t1 = time.perf_counter()
(211366, 3)
>>> t1-t0
0.1110449800034985
>>> t0 = time.perf_counter(); f_sparse(888.513, 900.099).shape; t1 = time.perf_counter()
(14577694, 3)
>>> t1-t0
3.6043728749500588

Answer 2

它对numpy 不太友好，但显而易见的答案是计算所需的整数范围。对于每个 x，计算最大的 |y|这样 x^2+y^2R^2 （到目前为止只是外半径）；通过向下舍入 sqrt(R*R-x*x) 来做到这一点。

然后，对于该范围内的每个 y，计算 |z| 的范围这样 r^2x^2+y^2+z^2R^2，通过向下舍入 sqrt(R*R-x*x-y*y) 和向上舍入 sqrt(max(r*r-x*x-y*y,0))。然后在两个结果范围内生成所有带有 z 的三元组（但不是 0 两次！）。