9

コンテキスト: 最近、alglibライブラリ(数値計算用) を発見しました。これは、私が探していたもの (ロバストな補間、データ分析...) のようで、実際には numpy や scipy で見つけることができませんでした。

ただし、(たとえば補間のために)有効な入力形式としてnumpy配列を受け入れず、通常のpythonリストオブジェクトのみを受け入れるという事実が心配です。

問題: コードとドキュメンテーションを少し掘り下げたところ、(予想通り) このリスト形式は移行用であることがわかりました。これは、ライブラリがとにかくそれを ctypes に変換するためです (cpython ライブラリは、基礎となる C の単なるインターフェイスです)。 /C++ ライブラリ)。

それが私の懸念です。コード内で、numpy 配列を使用しています。これは、実行している科学計算のパフォーマンスが大幅に向上するためです。したがって、alglibルーチンに渡されたデータをリスト (ctypes に変換される) に変換する必要があると、パフォーマンスに大きな影響を与えるのではないかと心配しています (内部に数十万の浮動小数点数を持つ可能性のある配列と、数千の浮動小数点数を持つ配列を扱っています)。アレイ)。

質問: 実際にパフォーマンスが低下すると思いますか?それとも、numpy 配列を受け入れて (numpy 配列から ctypes への) 変換を 1 つだけ行うように、 alglibコード (python インターフェイスのみ) の変更を開始する必要があると思いますか? )? これは非常に大きなライブラリであるため、これが実現可能かどうかさえわかりません...おそらく、より良いアイデアや提案があります（類似しているが異なるライブラリであっても）...

編集

私の問題があまり関心を集めていないか、私の質問が明確/関連性がないようです。または、誰も解決策やアドバイスを持っていないかもしれませんが、多くの専門家が周りにいるのではないかと思います:)とにかく、問題を説明するために、小さくて速くて汚いテストコードを書きました...

#!/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% の間で変動します。これは私にとって非常に大きな問題です...

4

3 に答える 3

5

numpy 配列のデータ バッファーをベクターのデータ ポインターに直接渡す独自のラップ関数を作成できます。これはデータをコピーせず、ラップ関数を大幅に高速化します。次のコードは、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
于 2012-02-27T01:39:11.580 に答える
3

C++ alglib が NumPy 配列を受け入れるようにすることは確かに可能です: SciPy はこれを行います。問題は、それがどれほど難しいかということです。次のような、半自動の C++ → Python ラッピング プログラムの 1 つを試してみるのもよいでしょう (最初に使用するものから始めます – 警告: 私は専門家ではありません)。

別の話題について: 私は過去に SciPy で補間スプラインを使用して成功しました。ただし、SciPy で必要なものがすべて見つかったわけではないため、これで十分かどうかはわかりません。

于 2012-02-26T00:47:45.333 に答える
1

EOLの回答に加えて、試すこともできます

NumPy 配列を処理するが、適切な引数で基礎となる C/C++ を呼び出す Python インターフェイスを生成するため。

以前にこれを行ったことがなく、C と Python のインターフェイスに関する豊富な経験がなくても、小さな科学的 C ライブラリに対してこれを行うには、ドキュメントが十分に明確であることがわかりました。

于 2012-02-26T03:25:22.187 に答える