关键词:
【中文标题】我可以使用 Numba、矢量化或多处理加速这种空气动力学计算吗?【英文标题】:Can I speed up this aerodynamics calculation with Numba, vectorization, or multiprocessing? 【发布时间】:2021-06-19 08:52:59 【问题描述】:问题:
我正在尝试提高 Python 中空气动力学函数的速度。
功能集:
import numpy as np
from numba import njit
def calculate_velocity_induced_by_line_vortices(
points, origins, terminations, strengths, collapse=True
):
# Expand the dimensionality of the points input. It is now of shape (N x 1 x 3).
# This will allow NumPy to broadcast the upcoming subtractions.
points = np.expand_dims(points, axis=1)
# Define the vectors from the vortex to the points. r_1 and r_2 now both are of
# shape (N x M x 3). Each row/column pair holds the vector associated with each
# point/vortex pair.
r_1 = points - origins
r_2 = points - terminations
r_0 = r_1 - r_2
r_1_cross_r_2 = nb_2d_explicit_cross(r_1, r_2)
r_1_cross_r_2_absolute_magnitude = (
r_1_cross_r_2[:, :, 0] ** 2
+ r_1_cross_r_2[:, :, 1] ** 2
+ r_1_cross_r_2[:, :, 2] ** 2
)
r_1_length = nb_2d_explicit_norm(r_1)
r_2_length = nb_2d_explicit_norm(r_2)
# Define the radius of the line vortices. This is used to get rid of any
# singularities.
radius = 3.0e-16
# Set the lengths and the absolute magnitudes to zero, at the places where the
# lengths and absolute magnitudes are less than the vortex radius.
r_1_length[r_1_length < radius] = 0
r_2_length[r_2_length < radius] = 0
r_1_cross_r_2_absolute_magnitude[r_1_cross_r_2_absolute_magnitude < radius] = 0
# Calculate the vector dot products.
r_0_dot_r_1 = np.einsum("ijk,ijk->ij", r_0, r_1)
r_0_dot_r_2 = np.einsum("ijk,ijk->ij", r_0, r_2)
# Calculate k and then the induced velocity, ignoring any divide-by-zero or nan
# errors. k is of shape (N x M)
with np.errstate(divide="ignore", invalid="ignore"):
k = (
strengths
/ (4 * np.pi * r_1_cross_r_2_absolute_magnitude)
* (r_0_dot_r_1 / r_1_length - r_0_dot_r_2 / r_2_length)
)
# Set the shape of k to be (N x M x 1) to support numpy broadcasting in the
# subsequent multiplication.
k = np.expand_dims(k, axis=2)
induced_velocities = k * r_1_cross_r_2
# Set the values of the induced velocity to zero where there are singularities.
induced_velocities[np.isinf(induced_velocities)] = 0
induced_velocities[np.isnan(induced_velocities)] = 0
if collapse:
induced_velocities = np.sum(induced_velocities, axis=1)
return induced_velocities
@njit
def nb_2d_explicit_norm(vectors):
return np.sqrt(
(vectors[:, :, 0]) ** 2 + (vectors[:, :, 1]) ** 2 + (vectors[:, :, 2]) ** 2
)
@njit
def nb_2d_explicit_cross(a, b):
e = np.zeros_like(a)
e[:, :, 0] = a[:, :, 1] * b[:, :, 2] - a[:, :, 2] * b[:, :, 1]
e[:, :, 1] = a[:, :, 2] * b[:, :, 0] - a[:, :, 0] * b[:, :, 2]
e[:, :, 2] = a[:, :, 0] * b[:, :, 1] - a[:, :, 1] * b[:, :, 0]
return e
上下文:
此函数由Ptera Software 使用,这是一个用于扑翼空气动力学的开源求解器。如下面的配置文件输出所示,它是迄今为止 Ptera Software 运行时间的最大贡献者。
目前,Ptera Software 运行一个典型案例只需 3 多分钟,我的目标是在 1 分钟内完成。
该函数接受一组点、起点、终点和强度。在每个点上,它都会找到由线涡流引起的感应速度,这些线涡流的特征在于起点、终点和强度的组。如果塌陷为真,则输出是由于涡流在每个点处引起的累积速度。如果为 false,则函数输出每个涡旋对每个点的速度的贡献。
在典型的运行过程中,速度函数被调用大约 2000 次。起初,调用涉及具有相对较小输入参数(大约 200 个点、起点、终点和强度)的向量。后来的调用涉及大量输入参数(大约 400 个点和大约 6,000 个起源、终止和强度)。理想的解决方案对于所有大小的输入都应该是快速的,但提高大输入调用的速度更为重要。
为了测试,我建议使用您自己的函数实现运行以下脚本:
import timeit
import matplotlib.pyplot as plt
import numpy as np
n_repeat = 2
n_execute = 10 ** 3
min_oom = 0
max_oom = 3
times_py = []
for i in range(max_oom - min_oom + 1):
n_elem = 10 ** i
n_elem_pretty = np.format_float_scientific(n_elem, 0)
print("Number of elements: " + n_elem_pretty)
# Benchmark Python.
print("\tBenchmarking Python...")
setup = '''
import numpy as np
these_points = np.random.random((''' + str(n_elem) + ''', 3))
these_origins = np.random.random((''' + str(n_elem) + ''', 3))
these_terminations = np.random.random((''' + str(n_elem) + ''', 3))
these_strengths = np.random.random(''' + str(n_elem) + ''')
def calculate_velocity_induced_by_line_vortices(points, origins, terminations,
strengths, collapse=True):
pass
'''
statement = '''
results_orig = calculate_velocity_induced_by_line_vortices(these_points, these_origins,
these_terminations,
these_strengths)
'''
times = timeit.repeat(repeat=n_repeat, stmt=statement, setup=setup, number=n_execute)
time_py = min(times)/n_execute
time_py_pretty = np.format_float_scientific(time_py, 2)
print("\t\tAverage Time per Loop: " + time_py_pretty + " s")
# Record the times.
times_py.append(time_py)
sizes = [10 ** i for i in range(max_oom - min_oom + 1)]
fig, ax = plt.subplots()
ax.plot(sizes, times_py, label='Python')
ax.set_xscale("log")
ax.set_xlabel("Size of List or Array (elements)")
ax.set_ylabel("Average Time per Loop (s)")
ax.set_title(
"Comparison of Different Optimization Methods\nBest of "
+ str(n_repeat)
+ " Runs, each with "
+ str(n_execute)
+ " Loops"
)
ax.legend()
plt.show()
以前的尝试:
我之前为加速此函数所做的尝试包括对其进行矢量化(效果很好,因此我保留了这些更改)并尝试了 Numba 的 JIT 编译器。我对 Numba 的结果好坏参半。当我尝试在整个速度函数的修改版本上使用 Numba 时,我的结果比以前慢得多。但是,我发现 Numba 显着加快了我在上面实现的叉积和范数函数。
更新:
更新 1:
根据 Mercury 的评论(已被删除),我替换了
points = np.expand_dims(points, axis=1)
r_1 = points - origins
r_2 = points - terminations
两次调用以下函数:
@njit
def subtract(a, b):
c = np.empty((a.shape[0], b.shape[0], 3))
for i in range(a.shape[0]):
for j in range(b.shape[0]):
for k in range(3):
c[i, j, k] = a[i, k] - b[j, k]
return c
这导致速度从 227 秒增加到 220 秒。这个更好!但是还是不够快。
我还尝试将 njit fastmath 标志设置为 true,并使用 numba 函数而不是调用 np.einsum。都没有提高速度。
更新 2:
有了 Jérôme Richard 的回答,运行时间现在是 156 秒,减少了 29%!我很满意接受这个答案,但如果您认为可以改进他们的工作,请随时提出其他建议!
【问题讨论】:
出色的工作矢量化了你所做的。这看起来不错的样子。我不是 numba 专家,但在某些情况下,我认为 numba 可以更好地 处理非矢量化代码。尽管可能很痛苦,但值得用 numba 恢复到普通 python 中的 for 循环,看看是否有帮助 很遗憾,我不知道答案。 简要地看一下 repo,你好像按顺序调用了这个函数 3 次,你是否考虑过并行化这些调用本身,即在单独的线程/进程中运行它们? github.com/camUrban/PteraSoftware/blob/… @wingedNorthropi 请注意,第一次调用 Numba 函数非常慢,因为必须编译代码。但是,您可以将编译后的代码放在缓存中以降低成本。或者,您可以为 Numba 函数提供类型,以便可以提前完成编译。最后,Numba 有时对代码进行矢量化的效率低于原生预编译 Numpy 调用。 @wingedNorthropi 答案已经使用了多个(Numba)线程,因此多处理不会帮助您的程序更快(至少对于此功能而言)。建议的解决方案仍然高度受内存限制。所以我认为在普通 CPU 上进一步改进代码的唯一方法是分解代码以便动态计算。 【参考方案1】:首先,如果您主要使用parallel=True
和prange
手动请求它,Numba 可以执行并行计算,从而产生更快的代码。这对大数组很有用(但对小数组没有用)。
此外,您的计算主要是内存限制。因此,当它们没有被多次重用时,或者更普遍地,当它们不能被重新计算时(以相对便宜的方式),你应该避免创建大数组。例如 r_0
就是这种情况。
此外,内存访问模式很重要:当访问在内存中连续并且缓存/RAM的使用效率更高时,向量化会更有效。因此,arr[0, :, :] = 0
应该比arr[:, :, 0] = 0
更快。同样,arr[:, :, 0] = arr[:, :, 1] = 0
应该比arr[:, :, 0:2] = 0
慢,因为前者执行非连续的内存传递,而后者只执行一个更连续的内存传递。有时,转置您的数据可能会有所帮助,以便以下计算更快。
此外,Numpy 倾向于创建许多分配成本高昂的临时数组。当输入数组很小时,这是一个大问题。在大多数情况下,Numba jit 可以避免这种情况。
最后,关于您的计算,将 GPU 用于大型数组(绝对不用于小型数组)可能是个好主意。你可以看看 cupy 或 clpy 来很容易地做到这一点。
这是一个在 CPU 上工作的优化实现:
import numpy as np
from numba import njit, prange
@njit(parallel=True)
def subtract(a, b):
c = np.empty((a.shape[0], b.shape[0], 3))
for i in prange(c.shape[0]):
for j in range(c.shape[1]):
for k in range(3):
c[i, j, k] = a[i, k] - b[j, k]
return c
@njit(parallel=True)
def nb_2d_explicit_norm(vectors):
res = np.empty((vectors.shape[0], vectors.shape[1]))
for i in prange(res.shape[0]):
for j in range(res.shape[1]):
res[i, j] = np.sqrt(vectors[i, j, 0] ** 2 + vectors[i, j, 1] ** 2 + vectors[i, j, 2] ** 2)
return res
# NOTE: better memory access pattern
@njit(parallel=True)
def nb_2d_explicit_cross(a, b):
e = np.empty(a.shape)
for i in prange(e.shape[0]):
for j in range(e.shape[1]):
e[i, j, 0] = a[i, j, 1] * b[i, j, 2] - a[i, j, 2] * b[i, j, 1]
e[i, j, 1] = a[i, j, 2] * b[i, j, 0] - a[i, j, 0] * b[i, j, 2]
e[i, j, 2] = a[i, j, 0] * b[i, j, 1] - a[i, j, 1] * b[i, j, 0]
return e
# NOTE: avoid the slow building of temporary arrays
@njit(parallel=True)
def cross_absolute_magnitude(cross):
return cross[:, :, 0] ** 2 + cross[:, :, 1] ** 2 + cross[:, :, 2] ** 2
# NOTE: avoid the slow building of temporary arrays again and multiple pass in memory
# Warning: do the work in-place
@njit(parallel=True)
def discard_singularities(arr):
for i in prange(arr.shape[0]):
for j in range(arr.shape[1]):
for k in range(3):
if np.isinf(arr[i, j, k]) or np.isnan(arr[i, j, k]):
arr[i, j, k] = 0.0
@njit(parallel=True)
def compute_k(strengths, r_1_cross_r_2_absolute_magnitude, r_0_dot_r_1, r_1_length, r_0_dot_r_2, r_2_length):
return (strengths
/ (4 * np.pi * r_1_cross_r_2_absolute_magnitude)
* (r_0_dot_r_1 / r_1_length - r_0_dot_r_2 / r_2_length)
)
@njit(parallel=True)
def rDotProducts(b, c):
assert b.shape == c.shape and b.shape[2] == 3
n, m = b.shape[0], b.shape[1]
ab = np.empty((n, m))
ac = np.empty((n, m))
for i in prange(n):
for j in range(m):
ab[i, j] = 0.0
ac[i, j] = 0.0
for k in range(3):
a = b[i, j, k] - c[i, j, k]
ab[i, j] += a * b[i, j, k]
ac[i, j] += a * c[i, j, k]
return (ab, ac)
# Compute `np.sum(arr, axis=1)` in parallel.
@njit(parallel=True)
def collapseArr(arr):
assert arr.shape[2] == 3
n, m = arr.shape[0], arr.shape[1]
res = np.empty((n, 3))
for i in prange(n):
res[i, 0] = np.sum(arr[i, :, 0])
res[i, 1] = np.sum(arr[i, :, 1])
res[i, 2] = np.sum(arr[i, :, 2])
return res
def calculate_velocity_induced_by_line_vortices(points, origins, terminations, strengths, collapse=True):
r_1 = subtract(points, origins)
r_2 = subtract(points, terminations)
# NOTE: r_0 is computed on the fly by rDotProducts
r_1_cross_r_2 = nb_2d_explicit_cross(r_1, r_2)
r_1_cross_r_2_absolute_magnitude = cross_absolute_magnitude(r_1_cross_r_2)
r_1_length = nb_2d_explicit_norm(r_1)
r_2_length = nb_2d_explicit_norm(r_2)
radius = 3.0e-16
r_1_length[r_1_length < radius] = 0
r_2_length[r_2_length < radius] = 0
r_1_cross_r_2_absolute_magnitude[r_1_cross_r_2_absolute_magnitude < radius] = 0
r_0_dot_r_1, r_0_dot_r_2 = rDotProducts(r_1, r_2)
with np.errstate(divide="ignore", invalid="ignore"):
k = compute_k(strengths, r_1_cross_r_2_absolute_magnitude, r_0_dot_r_1, r_1_length, r_0_dot_r_2, r_2_length)
k = np.expand_dims(k, axis=2)
induced_velocities = k * r_1_cross_r_2
discard_singularities(induced_velocities)
if collapse:
induced_velocities = collapseArr(induced_velocities)
return induced_velocities
在我的机器上,此代码比在大小为 10**3
的数组上的初始实现快 2.5 倍。它还使用了一点更少的内存。
【讨论】:
哇,这是一个巨大的改进。运行时间现在为 156 秒,增加了 29%。这是有道理的,因为您将大约 60% 的代码速度提高了 2.5 倍!我会将其添加为更新。使用多处理时避免重新编译 numba 代码
】使用多处理时避免重新编译numba代码【英文标题】:Avoidrecompilationofnumbacodewhenusingmultiprocessing【发布时间】:2020-12-3116:05:28【问题描述】:我一直在使用numba进行多处理。唯一的问题-numba会分别重新编译每个进程的代码。(当进... 查看详情
如何加速这种双 for 循环?
...imizationalgorithm。为了加快计算速度,我想对这个瓶颈进行矢量化处理。我知道N大约是k的一百倍。MyLoglik=0for(iinc(1:N))for(jinc(1:k))MyLoglik=MyLog 查看详情
Python:numba 可以在 nopython 模式下处理字符串数组吗?
】Python:numba可以在nopython模式下处理字符串数组吗?【英文标题】:Python:cannumbaworkwitharraysofstringsinnopythonmode?【发布时间】:2015-11-1010:35:27【问题描述】:我正在使用pandas0.16.2、numpy1.9.2和numba0.20。有没有办法让numba在nopython模式... 查看详情
使用numba加速opencvpython视频流代码。提升6.5倍性能(代码片段)
使用Numba对OpenCVPython视频处理代码加速。性能提升6.5倍目标问题:在OpenCVPython中视频处理是比较耗资源的,从而造成画面卡顿,如果跳帧处理可能造成丢失关键数据。用Numba对OpenCV代码加速是1个较好的改进方法。只须... 查看详情
我可以在 iOS 上进行这种代码向量化,还都有哪些替代方法?
...述】:我遇到了一个有趣的blogpost,他在谈论某种通过“矢量化代码”加快处理速度的精湛技术。很科学。他正在使用一种名为SS 查看详情
用于 SciPy 集成和插值的 Numba
...样条)中。我做了几百次这些集成,所以我认为这是Numba可以提升的东西。看起来Numb 查看详情
Numba 中的稀疏矩阵
...nNumba【发布时间】:2013-10-2513:21:20【问题描述】:我希望使用Numba(http://numba.pydata.org/)加速我的机器学习算法(用Python编写)。请注意,该算法将稀疏矩阵作为其输入数据。在我的纯Python实现中,我使用了来自Scipy的csr_matrix和相... 查看详情
numba安装和使用
numba是针对python加速的包,类似cython,pypy,优势是代码改动少首先要安装llvmliteapt-getinstallllvm-3.8LLVM_CONFIG=/usr/local/llvm38/3.8.1/lib/llvm-3.8/bin/llvm-configpipinstallllvmlite#看自己路径在哪程序里importnumba@numba.jit()装饰要加速的函数 查看详情
Python numpy:无法将 datetime64[ns] 转换为 datetime64[D](与 Numba 一起使用)
...我想将一个日期时间数组传递给一个Numba函数(它不能被矢量化,否则会很慢)。我了解Numba支持nump 查看详情
使用 Numba 处理 pandas DataFrame 时间序列的有效方法
】使用Numba处理pandasDataFrame时间序列的有效方法【英文标题】:EfficientwaytoprocesspandasDataFrametimeserieswithNumba【发布时间】:2014-07-0101:30:34【问题描述】:我有一个包含1,500,000行的DataFrame。这是我从QuantQuote.com购买的一分钟级股市数... 查看详情
使用带有 numba njit 功能的字典
】使用带有numbanjit功能的字典【英文标题】:UsingDictionarieswithnumbanjitfunction【发布时间】:2019-07-3109:30:11【问题描述】:当输入和返回是字典时,如何使用numba加速函数?我熟悉将numba用于接受数字并返回数组的函数,如下所示:... 查看详情
即使对于巨型矩阵,NUMBA CUDA 也比并行 CPU 慢
...,我发现它们都比并行CPU方法慢。带有CUDA目标和模板的矢量化甚至更糟,所以我尝试创建一个自定义内核。您随处可见的一篇博文是https://gist.gith 查看详情
用python做策略回测,耗时很长,有啥加速办法
...能。有几个numpy的加速包,比如numexpr.安装IntelMKL.最后,可以讲关键部分用c/c++实现。如果无法避开python的for,建议使用Numba来提速,理想情况下可以达到和numpy向量化差不多的速度。 查看详情
加速python中的元素数组乘法
...发布时间】:2013-10-1607:56:31【问题描述】:我一直在尝试使用numba和numexpr来加快简单的逐元素矩阵乘法。我一直没能得到更好的结果,它们基本上(速度方面)都相当于numpys乘法函数。有没有人在这方面有运气?我是否使用了num... 查看详情
对于纯 numpy 代码,使用 numba 的收益在哪里?
...加速纯numpy代码时的收益来自哪里。是否有任何分析工具可以让您查看jitted函数?演示代码(如下)只是使用非常基本的矩阵乘法来为计算机提供工作。观察到的收益来自:更 查看详情
使用矢量化通过反向乘法加速 matlab 代码
】使用矢量化通过反向乘法加速matlab代码【英文标题】:Speedupmatlabcodewithbackwardmultiplicationusingvectorization【发布时间】:2020-09-3001:04:49【问题描述】:我需要减少以下用Matlab编写的代码的运行时间:dt=0.001;dt05=dt^0.5;length_t=1.0e6;%a:ar... 查看详情
为啥 np.hypot 和 np.subtract.outer 与香草广播相比非常快?使用 Numba 并行加速 numpy 进行距离矩阵计算
】为啥np.hypot和np.subtract.outer与香草广播相比非常快?使用Numba并行加速numpy进行距离矩阵计算【英文标题】:Whynp.hypotandnp.subtract.outerveryfastcomparedtovanillabroadcast?UsingNumbaforspeedupnumpyinparallelfordistancematrixcalculation为什么np.hypot和np.subtr 查看详情
如何使用 numba 在 GPU 上泛化快速矩阵乘法
...程在他们的网站上阅读它,目前我坚持使用他们的示例,可以在这里找到:https://numba.pydata.org/numb 查看详情