【发布时间】：2017-11-30 03:03:07
【问题描述】：
作为在 cython 中重写我的游戏引擎的一部分，我正在尝试提高我的 python+numpy 类在矩阵和向量数学方面的性能，因为这是我之前遇到的主要瓶颈之一。这组模块为Vector2/3/4、Matrix2/3/4 和Quaternion 等类型定义了类。
从glMatrix javascript library 中获取一页，我认为这次我可以做的一件事是从基于类的系统切换到只有一堆数学函数的模块，以减少更多开销。这样一来，我不必每次将两个向量相加时都返回一个新对象，而不必构造自定义对象。
为了测试这一点，我编写了一个基准演示，用于创建两个 Vec2 对象 a 和 b 将它们按组件相加得到 Vec2 对象 out。用于此的代码被分解为用于计时的main.py、用于 cython 代码的 vec2.pyx 和用于 python 代码的 pyvec2.py。以下是每个组件的代码：
main.py
import time
import array
import math3d.vec2 as vec2
import math3d.pyvec2 as pyvec2
def test(n, func, param_list):
start = time.time()
for i in range(n):
func(*param_list)
end = time.time()
print func, end-start
test(1000000, pyvec2.pyadd, [[1, 2], [3, 4]])
test(1000000, pyvec2.pyadd2, [[0, 0], [1, 2], [3, 4]])
test(1000000, vec2.add, [[1, 2], [3, 4]])
test(1000000, vec2.add2, [array.array("f", [1, 2]), array.array("f", [3, 4])])
test(1000000, vec2.add3, [array.array("f", [1, 2]), array.array("f", [3, 4])])
test(1000000, vec2.add4, [array.array("f", [1, 2]), array.array("f", [3, 4])])
test(1000000, vec2.add5, [[0, 0], [1, 2], [3, 4]])
test(1000000, vec2.add6, [array.array("f", [0, 0]), array.array("f", [1, 2]), array.array("f", [3, 4])])
test(1000000, vec2.add7, [array.array("f", [0, 0]), array.array("f", [1, 2]), array.array("f", [3, 4])])
test(1000000, vec2.add8, [array.array("f", [0, 0]), array.array("f", [1, 2]), array.array("f", [3, 4])])
test(1000000, vec2.add9, [[0, 0], [1, 2], [3, 4]])
vec2.pyx
from libc.stdlib cimport malloc, free
from cpython cimport array
import array
def add(list a, list b):
cdef float[2] out = [0, 0]
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
def add2(float[:] a, float[:] b):
cdef float[2] out = [0, 0]
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
def add3(array.array a, array.array b):
cdef float[2] out = [0, 0]
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
def add4(array.array a, array.array b):
cdef array.array out = array.array("f", [0, 0])
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
def add5(list out, list a, list b):
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
def add6(float[:] out, float[:] a, float[:] b):
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
def add7(array.array out, array.array a, array.array b):
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
def add8(array.array out, array.array a, array.array b):
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
def add9(out, a, b):
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
pyvec2.py
def pyadd(a, b):
out = [a[0] + b[0], a[1] + b[1]]
def pyadd2(out, a, b):
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
这里是运行main.py后的结果：
<function pyadd at 0x0000000003354828> 0.380000114441
<function pyadd2 at 0x0000000003354908> 0.31299996376
<built-in function add> 0.261000156403
<built-in function add2> 0.680999994278
<built-in function add3> 0.268000125885
<built-in function add4> 0.601000070572
<built-in function add5> 0.144999980927
<built-in function add6> 1.06299996376
<built-in function add7> 0.241000175476
<built-in function add8> 0.237999916077
<built-in function add9> 0.141000032425
由此可见，使用 python 列表似乎比类型化数组更快！此外，只是盲目地编译我的 python 函数而不输入 cython 会产生最好的结果！看起来执行代码所花费的大部分时间都用于在 python 类型之间进行转换。
因此，我很想知道是否有更快的方法在 cython 端执行数学运算，同时最大限度地减少传入参数的 python 开销。我对必须公开列表或 array.array 对象以直接访问我的向量或矩阵的内容并不真正感兴趣。将指向和来自 python 的指针传递到我的 cython 数学模块是理想的，但这似乎是不可能的，因为指针不是 python 对象。任何建议将不胜感激。
更新：
下面是我的Vec2 课程的代码。它由两个文件组成：第一个是基本_Vec 类，第二个从它继承为特定的Vec2 类。
vec.py
import math
import numpy as np
import random
class _Vec(object):
def __init__(self, *args):
try:
data, = args
data = np.array(data, dtype=np.float32)
cls_name = self.__class__.__name__
vec2_check = cls_name == "Vec2" and len(data) != 2
vec3_check = cls_name == "Vec3" and len(data) != 3
vec4_check = cls_name == "Vec4" and len(data) != 4
if any([vec2_check, vec3_check, vec4_check]) == True:
raise TypeError("{0} is not a valid {1}".format(data, cls_name))
except ValueError:
data = np.array(args, dtype=np.float32)
self._data = data
def __add__(self, other):
if isinstance(other, self.__class__):
return self.__class__(self._data + other._data)
return self.__class__(self._data + other)
def __radd__(self, other):
return self.__class__(other + self._data)
def __sub__(self, other):
if isinstance(other, self.__class__):
return self.__class__(self._data - other._data)
return self.__class__(self._data - other)
def __rsub__(self, other):
return self.__class__(other - self._data)
def __mul__(self, other):
if isinstance(other, self.__class__):
return self.__class__(self._data * other._data)
return self.__class__(self._data * other)
def __rmul__(self, other):
return self.__class__(other * self._data)
def __div__(self, other):
if isinstance(other, self.__class__):
return self.__class__(self._data / other._data)
return self.__class__(self._data / other)
def __rdiv__(self, other):
return self.__class__(other / self._data)
def __neg__(self):
return self.__class__(-self._data)
def __pos__(self):
return self.__class__(+self._data)
def __eq__(self, other):
return np.array_equal(self._data, other._data)
def __ne__(self, other):
return not self.__eq__(other)
def __lt__(self, other):
return self.square_length() < other.square_length()
def __le__(self, other):
return self.square_length() <= other.square_length()
def __gt__(self, other):
return self.square_length() > other.square_length()
def __ge__(self, other):
return self.square_length() >= other.square_length()
def __repr__(self):
return "{0}(data={1})".format(self.__class__.__name__, self.get_data())
def __str__(self):
return np.array_str(self._data)
def ceil(self):
return self.__class__(np.ceil(self._data))
def floor(self):
return self.__class__(np.floor(self._data))
def get_data(self):
return self._data.flatten().tolist()
def inverse(self):
return self.__class__(1.0/self._data)
def length(self):
return float(np.linalg.norm(self._data))
def negate(self):
return self.__class__(-self._data)
def normalize(self):
length = self.length()
if length == 0.0:
return self.__class__(np.zeros(self._data.shape()))
return self.__class__(self._data/length)
def round(self, decimal=0):
return self.__class__(np.round(self._data, decimal))
def square_length(self):
return float(np.sum(np.square(self._data)))
@classmethod
def distance(cls, a, b):
c = b - a
return c.length()
@classmethod
def dot(cls, a, b):
return float(np.dot(a._data, b._data))
@classmethod
def equals(cls, a, b, tolerance=0.0):
diffs = np.fabs((a - b)._data)
pairs = zip(list(np.fabs(a._data)), list(np.fabs(b._data)))
tolerance_calcs = [tolerance * max(1, a_val, b_val) for (a_val, b_val) in pairs]
tests = [d <= t for (d, t) in zip(diffs, tolerance_calcs)]
return all(tests)
@classmethod
def lerp(cls, a, b, t):
return a*(1-t) + b*t
@classmethod
def max_components(cls, a, b):
return cls(np.maximum(a._data, b._data))
@classmethod
def min_components(cls, a, b):
return cls(np.minimum(a._data, b._data))
@classmethod
def random(cls, n):
return cls(np.random.rand((n)))
@classmethod
def square_distance(cls, a, b):
c = b - a
return c.square_length()
vec2.py
from vec import _Vec
from vec3 import Vec3
from vec4 import Vec4
import math
import numpy as np
import random
class Vec2(_Vec):
@property
def x(self):
return float(self._data[0])
@x.setter
def x(self, value):
self._data[0] = float(value)
@property
def y(self):
return float(self._data[1])
@y.setter
def y(self, value):
self._data[1] = float(value)
def __repr__(self):
return "Vec2(x={0}, y={1})".format(self.x, self.y)
def transform_mat2(self, a):
prod = np.dot(a._data.T, self._data.T).T
return Vec2(prod)
def transform_mat3(self, a):
v3 = Vec3(self.get_data() + [1])
prod = np.dot(a._data.T, v3._data.T).T
return Vec2(prod[0:2])
def transform_mat4(self, a):
v4 = Vec4(self.get_data() + [0, 1])
prod = np.dot(a._data.T, v4._data.T).T
return Vec2(prod[0:2])
@classmethod
def random(cls):
return super(Vec2, cls).random(2)
更新 2：
正如 kazemakase 所指出的，main.py 中的一些值是整数。通过附加 .0 将所有内容定义为浮点数，我得到了以下时间：
<function pyadd at 0x0000000002FF4828> 0.384000062943
<function pyadd2 at 0x0000000002FF4908> 0.332000017166
<built-in function add> 0.227999925613
<built-in function add2> 0.640000104904
<built-in function add3> 0.258999824524
<built-in function add4> 0.556999921799
<built-in function add5> 0.145999908447
<built-in function add6> 0.983999967575
<built-in function add7> 0.217000007629
<built-in function add8> 0.236000061035
<built-in function add9> 0.131000041962
这些似乎与原始时间相似，在这种情况下 5 和 9 更快。
更新 3：
正如 BrenBarn 所指出的，解释我如何使用我原来的 python+numpy 类以及为什么我在它的性能上苦苦挣扎可能会很有用。最初，我的整个 3d 游戏库项目都是纯 Python 的，使用 PyOpenGL 来渲染图形。为了在我的 3d 世界中定位网格/模型，每一帧，我都需要计算一个 Mat4 变换矩阵，该矩阵定义该对象在世界中的位置、旋转和比例，然后上传到 GPU。当我定位许多对象（> 1000）时，我的应用程序中的帧速率会变慢。虽然暂时禁用 3d 确实在一定程度上提高了性能，但我的应用程序仍然落后于 60 fps 以下。那时我意识到在 python 中简单地计算Mat4 矩阵是罪魁祸首。当我删除这些计算并在原点处绘制所有 3d 对象时，性能又恢复了。
想知道如果我的所有 math3d 库都很慢，我想我会开始对所有这些类进行基准测试。我为 Vec2 类的 add 函数做了以下基准测试，因为它是最简单且计算成本最低的测试函数：
import time
from vec2 import Vec2
def add(a, b):
out = a + b
return out
def test(n, func, param_list):
start = time.time()
for i in range(n):
func(*param_list)
end = time.time()
print func, end-start
test(1000000, add, [Vec2([1.0, 2.0]), Vec2([3.0, 4.0])])
#<function add at 0x00000000041B1CF8> 2.81699991226 (using the real def add(a, b) function)
#<function add at 0x000000000362CCF8> 0.168999910355 (just passing in values to def add(a, b): pass)
那些 cmets 显示了我将两个 Vec2 加在一起所获得的时间。由此，我得出结论，数学很慢，python 函数调用 + 类开销很大，这是我的 3d 库中要考虑的性能瓶颈。我希望这能为这个问题提供一些理由。
更新 4：
Paul Cornelius 所说的关于从 python 到 cython 的翻译的惩罚只有一次让我思考：为什么不在我创建“Vec2 对象”时获取指向float * 对象的“指针”？然后，我可以为将来的数学运算传入指针，cython 可以取消引用这些指针以获取实际数据，然后可以执行数学运算。结果如下代码：
main.py
import time
import array
import math3d.vec2 as vec2
def make_list(a, b):
out = [a, b]
def test(n, func, param_list):
start = time.time()
for i in range(n):
func(*param_list)
end = time.time()
print func, end-start
test(1000000, vec2.create, [1, 2])
test(1000000, make_list, [1, 2])
a = vec2.create(1, 2)
#b = vec2.get_data(a)
b = vec2.create(3, 4)
c = vec2.create(0, 0)
test(1000000, vec2.add, [c, a, b])
test(1000000, vec2.add2, [[0, 0], [1, 2], [3, 4]])
print vec2.get_data(c)
vec2.pyx
def add(uintptr_t out, uintptr_t a, uintptr_t b):
cdef float *a_data = <float *>a
cdef float *b_data = <float *>b
cdef float *out_data = <float *>out
out_data[0] = a_data[0] + b_data[0]
out_data[1] = a_data[1] + b_data[1]
def create(float x, float y):
cdef float* a = <float *>malloc(sizeof(float))
a[:] = [x, y]
cdef uintptr_t u_ptr = <uintptr_t> a
return u_ptr
def get_data(uintptr_t u_ptr):
cdef float *b = <float *>u_ptr
return b[0], b[1]
def add2(out, a, b):
out[0] = a[0] + b[0]
out[1] = a[1] + b[1]
return out
计时结果
<built-in function create> 0.19000005722
<function make_list at 0x0000000002994828> 0.269999980927
<built-in function add> 0.111000061035
<built-in function add2> 0.141999959946
(4.0, 6.0)
当然，在 python 中处理 uintptr_t 整数是非常不安全的，因为它们可能会在 python 端使用 + 运算符错误地添加在一起。此外，尚不清楚这带来的轻微性能优势（100 万次操作只需 0.03 秒）是否真的值得。
【问题讨论】：
-
您不太可能从头开始编写比 numpy 更快的数学库。你可能会更好地展示你正在使用的“python+numpy 类”并尝试改进它们，假设它们是某种围绕 numpy 对象的包装器。
-
@BrenBarn 我为我的
Vec2类包含了我原来的“python+numpy 类”实现。希望这会有所帮助！ -
您将带有整数的列表传递给某些
addx函数。我猜情况 5 和 9 更快，因为它们执行整数加法而不是浮点运算。 -
我认为您应该问一个单独的问题，在该问题中展示自定义类的最简单示例，展示您如何尝试使用它们，并解释为什么您认为它们是性能瓶颈。猜测一下，从您发布的代码来看，您似乎正在 Python 类中设置便利方法/属性。如果你经常使用这些，它会比直接使用 numpy 对象慢。
-
@CodeSurgeon 很有趣。我想这些差异可能是由函数调用期间的类型检查/转换引起的。 （出于好奇，您可以尝试删除实际的
out=a+b代码和时间空函数。）请注意，无论如何，您的方法不太可能给您带来非常好的性能。add函数非常短，您仍然可以从 Python 单独调用它们，这会导致大量开销。为了充分利用 Cython，您应该尽可能长时间地留在编译后的代码中，而不要回调 Python。