python - 使用叠加模拟球体周围的扬声器-需要提高速度

Question

注意：自发布以来，速度有了极大提高，请参见底部的编辑。

我在利用循环方面有一些有效的代码，我敢肯定应该有一种更快的方法。输出数组的大小最终变得非常大，因此当我尝试使其他数组的大小与输出大小相同时，我很快就会耗尽内存。

我正在模拟围绕球形指向中心的许多扬声器。我有一个扬声器的仿真，我想通过使用叠加原理来利用此单个仿真。基本上，我想总结单个换能器仿真的旋转副本，以获得最终结果。

我在圆柱坐标（“ polar_coord_r”，“ polar_coord_z”）中进行了声压数据的轴对称模拟。来自模拟的压力场在每个R和Z值处都是唯一的，并由2D数组（“ P_real_RZ”）完整描述。

我想将这个压力场的多个旋转副本汇总到3D笛卡尔输出数组中。每个副本都旋转到球体上的不同位置。当前，我指定一个带有x，y，z点的旋转，因为它允许我进行矢量数学运算（球面坐标不允许我如此优雅地执行此操作）。输出阵列相当大（770×770×804）。

我有一些工作代码可以从扬声器的单个副本（“换能器”）获取输出。每个片段大约需要12秒，因此添加新扬声器需要两个多小时！我想要一打扬声器的副本，所以需要很长时间。

该代码获取具有常量X的切片，并计算该切片中每个位置的R和Z位置。 “ r_distance”是一个2D数组，其中包含距在原点和一个点（“点”）之间经过的线的径向距离。类似地，“ z_distance”是一个2D数组，其中包含沿着同一条线的距离。

为了获得切片的压力，我找到与计算出的R距离和Z距离最匹配的“ polar_coord_r”和“ polar_coord_z”的索引。我使用这些索引来查找要在输出中的每个值上放置什么压力值（来自P_real_RZ）。

一些定义：


xx，yy和zz是一维数组，用于描述通过输出量的切片
XXX，YYY和ZZZ是numpy.meshgrid生成的3D数组
point是定义扬声器旋转方向的点。基本上，这只是扬声器中心的位置向量。
P_real_RZ是一个2D数组，用于指定每个唯一R和Z值的实际压力。
polar_coord_r和polar_coord_z是一维数组，用于定义R和Z的唯一值，并在这些值上定义了P_real_RZ。
current_transducer（此代码到目前为止仅代表一个）是当前点的压力值计算机。
输出是将多个扬声器/换能器加在一起的结果。


任何加快该代码的建议将不胜感激。

工作循环：

for i, x in enumerate(xx):
    # Creates a unit vector from origin to a point
    vector = normalize(point)

    # Create a slice of the cartesian space with constant X
    xyz_slice = np.array([x*np.ones_like(XXX[i,:,:]), YYY[i,:,:], ZZZ[i,:,:]])

    # Projects the position vector of each point of the slice onto the unit vector.
    projection = np.array(list(map(np.dot, xyz_slice, vector )))

    # Normalizes the projection which results in the Z distance along the line passing through the point
    #z_distance = np.apply_along_axis(np.linalg.norm, 0, projection) # this is the slow bit
    z_distance = np.linalg.norm(projection, axis=0) # I'm an idiot

    # Uses vector math to determine the distance from the line
    # Each point in the XYZ slice is the sum of vector along the line and the vector away from the line (radial vector).
    # By extension the position of the xyz point minus the projection of the point against the unit vector, results in the radial vector
    # Norm the radial vector to get the R value for everywhere in the slice
    #r_distance = np.apply_along_axis(np.linalg.norm, 0, xyz_slice - projection) # this is the slow bit
    r_distance = np.linalg.norm(xyz_slice - projection, axis=0) # I'm an idiot

    # Map the pressure data to each point in the slice using the R and Z distance with the RZ pressure slice.
    # look for a more efficient way to do this perhaps. currently takes about 12 seconds per slice
    r_indices = r_map_v(r_distance) # 1.3 seconds by itself
    z_indices = z_map_v(z_distance)
    r_indices[r_indices>384] = 384 # find and remove indicies above the max for r_distance
    z_indices[r_indices>803] = 803
    current_transducer[i,:,:] = P_real_RZ[z_indices, r_indices]

    # Sum the mapped pressure data into the output.
    output += current_transducer

我还尝试使用3D笛卡尔数组形式的模拟数据。这就是来自模拟的所有x，y和z值的压力数据都与输出大小相同。我可以沿一个方向旋转3D阵列（布置在球体上的扬声器不需要两个旋转）。 waaaay占用了太多内存，并且仍然很慢。我最终通过这种方法收到内存错误。

编辑：我发现了一种稍微简单的方法，但是仍然很慢。我已经更新了上面的代码，以便不再有嵌套循环。

我运行了一个行探查器，到目前为止最慢的行是包含np.apply_along_axis（）的两行。恐怕我可能不得不重新考虑如何完全做到这一点。

最终编辑：我最初有一个嵌套循环，我以为是问题所在。我不知道是什么让我觉得我需要对linalg.norm使用apply_along_axis。无论如何，这就是问题所在。

Answer 1

我还没有找到优化该代码的所有方法，但是这个问题突然出现了：“我运行了一个行探查器，到目前为止，最慢的行是两个包含np.apply_along_axis（）的行。” np.linalg.norm接受axis参数。您可以更换线

z_distance = np.apply_along_axis(np.linalg.norm, 0, projection)

与

z_distance = np.linalg.norm(projection, axis=0)

（对于np.apply_along_axis和np.linalg.norm的其他用法也是如此）。
那应该可以提高性能。