Python + alglib + NumPy: как избежать преобразования массивов в списки? - PullRequest
9 голосов
/ 26 февраля 2012

Контекст: Я недавно обнаружил alglib библиотека (для численных расчетов), которая, кажется, мне нужнаискал (надежная интерполяция, анализ данных ...) и не мог найти ничего в numpy или scipy.

Однако меня беспокоит тот факт, что (например, для интерполяции) он не принимает numpyмассив в качестве допустимого формата ввода, но только обычные объекты списка Python.

Проблема: Я немного покопался в коде и документации и нашел (как и ожидалось)) что этот формат списка предназначен только для перехода, так как библиотека в любом случае преобразует его в ctypes (библиотека cpython является просто интерфейсом для базовой библиотеки C / C ++).

Вот что меня беспокоит: внутримой код, я работаю с пустыми массивами, потому что это значительно повышает производительность для научных расчетов, которые я выполняю на нем.Поэтому я боюсь, что необходимость преобразования любых данных, переданных в подпрограммы alglib в список (который будет преобразован в ctypes), окажет огромное влияние на производительность (я работаю с массивами, которыеможет иметь сотни тысяч плавающих внутри и тысячи массивов).

Вопрос: Как вы думаете, у меня действительно будет потеря производительности, или вы думаете, что я должен начать изменять alglib код (только интерфейс python), чтобы он мог принимать числовые массивы и выполнять только одно преобразование (из числовых массивов в 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%, что для меня огромно ...

Ответы [ 3 ]

5 голосов
/ 27 февраля 2012

Вы можете создать свою собственную функцию обтекания, которая передает буфер данных массива 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
3 голосов
/ 26 февраля 2012

Заставить C ++ alglib принимать массивы NumPy, безусловно, выполнимо: SciPy делает это.Вопрос действительно в том, насколько это сложно.Возможно, вы захотите попробовать одну из полуавтоматических программ обертывания C ++ → Python, например (начиная с той, с которой я бы начал - предупреждение: я не эксперт):

По другой теме: я использовал интерполяционные сплайныв SciPy с успехом, в прошлом.Я не уверен, что этого будет достаточно для ваших нужд, поскольку вы не нашли в SciPy всего, что хотели.

1 голос
/ 26 февраля 2012

В дополнение к ответам EOL вы также можете попробовать

для генерации интерфейса Python, который работает с массивами NumPy, но вызывает соответствующий C / C ++ с соответствующими аргументами.

Я обнаружил, что документы достаточно ясны, чтобы сделать это для небольшой научной библиотеки C, даже не делая этого раньше или не имея большого опыта взаимодействия с C и Python.

...