Python + Alglib + numpy:如何避免将数组转换为列表?
-
12-11-2019 - |
题
上下文: 我最近发现了 alglib library (用于数值计算)似乎是我正在寻找的东西(强大的插值,数据分析......),无法在Numpy或Scipy中找到。 但是,我担心(例如,用于插值)的事实,它不接受numpy数组作为有效输入格式,但是只有常规python列表对象。
问题: 我已经挖掘了代码和文档,发现(按预期),此列表格式仅用于转换,因为库将无论如何将其转换为CTypes(CPython库只是底层C / C ++的接口图书馆)。
那是我担心的地方:在我的代码里面,我正在使用Numpy阵列,因为它是我正在执行的科学计算的大表现提升。因此,我担心必须将传递给 的任何数据转换为 例程到列表中(将转换为ctypes)将对性能产生巨大影响(我正在使用阵列可以在里面有数百万岁漂浮,千万阵列)。
问题: 您是否认为我确实会有绩效损失,或者您认为我应该开始修改 alglib 代码(只有python接口),以便它可以接受numpy arrays ,并只制作一个转换(从Numpy阵列到CTypes)?我甚至不知道这是可行的,因为它是一个很大的图书馆...... 也许你们有更好的想法或建议(即使在类似但不同的图书馆)...
编辑
似乎我的问题并没有得到很多兴趣,或者我的问题并不清楚/相关。或者也许没有人有解决方案或建议,但我怀疑周围有这么多专家:) 无论如何,我写了一个小型,快速肮脏的测试代码来说明问题...
#!/usr/bin/env python
import xalglib as al
import timeit
import numpy as np
def func(x):
return (3.14 *x**2.3 + x**3 -x**2.34 +x)/(1.+x)**2
def fa(x, y, val=3.14):
s = al.spline1dbuildakima(x, y)
return (al.spline1dcalc(s, val), func(val))
def fb(x, y, val=3.14):
_x = list(x)
_y = list(y)
s = al.spline1dbuildakima(_x, _y)
return (al.spline1dcalc(s, val), func(val))
ntot = 10000
maxi = 100
x = np.random.uniform(high=maxi, size=ntot)
y = func(x)
xl = list(x)
yl = list(y)
print "Test for len(x)=%d, and x between [0 and %.2f):" % (ntot, maxi)
print "Function: (3.14 *x**2.3 + x**3 -x**2.34 +x)/(1.+x)**2"
a, b = fa(xl, yl)
err = np.abs(a-b)/b * 100
print "(x=3.14) interpolated, exact =", (a, b)
print "(x=3.14) relative error should be <= 1e-2: %s (=%.2e)" % ((err <= 1e-2), err)
if __name__ == "__main__":
t = timeit.Timer(stmt="fa(xl, yl)", setup="from __main__ import fa, xl, yl, func")
tt = timeit.Timer(stmt="fb(x, y)", setup="from __main__ import fb, x, y, func")
v = 1000 * t.timeit(number=100)/100
vv = 1000 * tt.timeit(number=100)/100
print "%.2f usec/pass" % v
print "%.2f usec/pass" % vv
print "%.2f %% less performant using numpy arrays" % ((vv-v)/v*100.)
.
并运行,我得到了:
"""
Test for len(x)=10000, and x between [0 and 100.00):
Function: (3.14 *x**2.3 + x**3 -x**2.34 +x)/(1.+x)**2
(x=3.14) interpolated, exact = (3.686727834705164, 3.6867278531266905)
(x=3.14) relative error should be <= 1e-2: True (=5.00e-07)
25.85 usec/pass
28.46 usec/pass
10.09 % less performant using numpy arrays
"""
.
性能损失振荡在约8%和14%之间,这对我来说是巨大的...
解决方案
使c ++ alglib接受numpy阵列肯定是可行的:scipy这样做。问题真的是多么困难。您可能希望尝试一个半自动C ++→Python包装程序,如(从我开始用警告开始的那样:我没有专家):
- cython
- boost.python
- 编织 在一个不同的主题上:过去,我已经使用了Scipy中的插值样条曲线。我不确定这对您的需求来说是足够的,因为你没有找到你在Scipy中想要的一切。
其他提示
您可以创建自己的包装函数直接将numpy数组的数据缓冲区传递给向量的数据指针,这不会复制数据,并加快包装函数很多。以下代码传递给x_vector.ptr.p_ptr,其中x是numpy数组。
当您传递Numpy数组时,必须确保数组具有连续内存中的元素。以下代码不检查一下。
import xalglib as al
import numpy as np
import ctypes
def spline1dbuildakima(x, y):
n = len(x)
_error_msg = ctypes.c_char_p(0)
__c = ctypes.c_void_p(0)
__n = al.c_ptrint_t(n)
__x = al.x_vector(cnt=n, datatype=al.DT_REAL, owner=al.OWN_CALLER,
last_action=0,ptr=al.x_multiptr(p_ptr=x.ctypes.data))
__y = al.x_vector(cnt=n, datatype=al.DT_REAL, owner=al.OWN_CALLER,
last_action=0,ptr=al.x_multiptr(p_ptr=y.ctypes.data))
al._lib_alglib.alglib_spline1dbuildakima(
ctypes.byref(_error_msg),
ctypes.byref(__x),
ctypes.byref(__y),
ctypes.byref(__n),
ctypes.byref(__c))
__r__c = al.spline1dinterpolant(__c)
return __r__c
def func(x):
return (3.14 *x**2.3 + x**3 -x**2.34 +x)/(1.+x)**2
def fa(x, y, val=3.14):
s = spline1dbuildakima(x, y)
return al.spline1dcalc(s, val), func(val)
def fb(x, y, val=3.14):
s = al.spline1dbuildakima(x, y)
return al.spline1dcalc(s, val), func(val)
ntot = 10000
maxi = 100
x = np.random.uniform(high=maxi, size=ntot)
y = func(x)
xl = list(x)
yl = list(y)
import time
start = time.clock()
for i in xrange(100):
a, b = fa(x, y)
print time.clock()-start
err = np.abs(a-b)/b * 100
print a, b, err
start = time.clock()
for i in xrange(100):
a, b = fb(xl, yl)
print time.clock()-start
err = np.abs(a-b)/b * 100
print a, b, err
.
输出是:
0.722314760822 <- seconds of numpy array version
3.68672728107 3.68672785313 1.55166878281e-05
3.22011891502 <- seconds of list version
3.68672728107 3.68672785313 1.55166878281e-05
. 除了EOL的答案,您还可以尝试
- swig用numpy.i
为了生成处理numpy阵列的Python接口,但调用具有适当参数的底层C / C ++。
我发现文档足够明确地为一个小科学C库做到这一点,而不是在以前或拥有接口C和Python的巨大经验。