Python + alglib + NumPy：如何避免将数组转换为列表？

Posted 2023-02-23

【中文标题】Python + alglib + NumPy：如何避免将数组转换为列表？

【英文标题】：Python + alglib + NumPy: how to avoid converting arrays to lists?

【发布时间】：2012-02-25 23:05:46

【问题描述】：

背景： 我最近发现了alglib library（用于数值计算），这似乎是我一直在寻找的东西（稳健插值、数据分析......），但在 numpy 或 scipy 中找不到。

但是，我担心（例如，用于插值）它不接受 numpy 数组作为有效输入格式，而是仅常规 python 列表对象。

问题： 我对代码和文档进行了一些研究，发现（如预期的那样）这种列表格式仅用于转换，因为该库无论如何都会将其转换为 ctypes（cpython 库只是底层 C/C++ 的接口图书馆）。

这就是我担心的地方：在我的代码中，我正在使用 numpy 数组，因为它对我正在执行的科学计算来说是一个巨大的性能提升。因此，我担心必须将传递给 alglib 例程的任何数据转换为列表（将转换为 ctypes）会对性能产生巨大影响（我正在使用内部可能有数十万个浮点数和数千个数组的数组）。

问题： 你认为我确实会有性能损失，还是你认为我应该开始修改 alglib 代码（仅 python 接口），以便它可以接受 numpy 数组，并且只进行一次转换（从 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% 之间波动，这对我来说是巨大的......

【参考方案1】：

您可以创建自己的 wrap 函数，将 numpy 数组的数据缓冲区直接传递给向量的数据指针，这不会复制数据，并且大大加快了您的 wrap 函数。以下代码将 x.ctypes.data 传递给 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

【讨论】：

您可以致电x = np.ascontiguousarray(x) 以确保它在内存中是连续的。

哇，太好了！非常感谢海瑞。我以前从未使用过包装方法，所以我正在阅读 EOL 很好地指出的一些文档。你的例子对我来说来得正是时候，所以我现在有一些具体的东西可以训练。

@Sebastian 谢谢你的评论。这对我也有帮助。

【参考方案2】：

让 C++ alglib 接受 NumPy 数组当然是可行的：SciPy 就是这样做的。问题真的是它有多难。您可能想尝试一种半自动的 C++ → Python 包装程序，例如（从我要开始的那个开始——警告：我不是专家）：

Cython Boost.Python Weave

关于另一个主题：过去，我在 SciPy 中成功使用了插值样条曲线。不过，我不确定这是否足以满足您的需求，因为您没有在 SciPy 中找到所需的一切。

【讨论】：

感谢您的回答 EOL。 SciPy 包中的样条线按预期工作，但我缺少一些像 Akima 这样的样条线（我确实有一些数据需要的不仅仅是线性插值，但不够“表现良好”，无法使用三次或更高阶的样条线得到一些不需要的伪影，例如结果中的振荡）。此外，alglib 还具有 SciPy 所没有的其他一些数据分析和统计工具，因此对我来说使用这样的外部库是有意义的。对于包装方面，似乎 alglib 开发人员已经在使用一个用于 C->Python...

@mhavel：可以包装 C++ 代码，使其直接使用 NumPy 数组而不是 Python 列表；令人惊讶的是，alglib 的 Python 包装器不这样做/不允许这样做。包装您自己需要的函数可能比修改 alglib Python 包装器更容易：如果需要，您将能够更新包装器以更轻松地遵循 alglib 的更新。我列出的三个包装器应该对您有很大帮助。

哦，好吧:) 我不是 Python 方面的专家，更不用说包装了。所以我可以试一试。非常感谢。

@mhavel：谢谢。我不是为 Python 包装 C++ 代码的专家，但对于对快速科学计算感兴趣的人来说，研究一些列出的 C++ 包装器可能是一项不错的投资。

【参考方案3】：

除了EOL的答案你也可以试试

SWIG with numpy.i

为了生成处理 NumPy 数组但使用适当参数调用底层 C/C++ 的 Python 接口。

我发现文档足够清晰，可以为小型科学 C 库执行此操作，而之前从未这样做过，也没有丰富的 C 和 Python 接口经验。

【讨论】：

感谢您的回答，我会看看这个。