Я боролся с проблемой, из-за которой у меня проблемы с намотанием головы, и поэтому я не знаю, как начать ее решать. Мой опыт программирования на C очень ограничен, и я думаю, что это причина, по которой я не могу добиться прогресса.
У меня есть функция, которая использует numpy.interp и scipy.integrate.quad для выполнения определенного интеграла. Так как я использую quad для интеграции, и согласно его документации:
Функция или метод Python для интеграции. Если func принимает много аргументов, он интегрируется вдоль оси, соответствующей первому аргументу.
Если пользователь желает улучшить производительность интеграции, тогда f может быть scipy.LowLevelCallable с одной из сигнатур:
double func(double x)
double func(double x, void *user_data)
double func(int n, double *xx)
double func(int n, double *xx, void *user_data)
user_data - это данные, содержащиеся в scipy.LowLevelCallable. В формах вызова с xx, n - это длина массива xx, который содержит xx[0] == x, а остальные элементы являются числами, содержащимися в аргументе args аргумента quad.
Кроме того,некоторые сигнатуры вызовов ctypes поддерживаются для обратной совместимости, но их не следует использовать в новом коде.
Мне нужно использовать объекты scipy.LowLevelCallable для ускорения моего кода, и мне нужно придерживатьсямоя функция дизайна одной из вышеперечисленных подписей. Более того, поскольку я не хочу усложнять все это с библиотеками C и компиляторами, я хочу решить эту проблему "на лету" с помощью инструментов, доступных из numba, в частности numba.cfunc, что позволяет мнепередать Python C. API.
Я смог решить эту проблему для интегральной функции, которая принимает в качестве входных данных переменную интегрирования и произвольное число скалярных параметров:
from scipy import integrate, LowLevelCallable
from numba import njit, cfunc
from numba.types import intc, float64, CPointer
def jit_integrand_function(integrand_function):
jitted_function = njit(integrand_function)
@cfunc(float64(intc, CPointer(float64)))
def wrapped(n, xx):
return jitted_function(xx[0], xx[1], xx[2], xx[3])
return LowLevelCallable(wrapped.ctypes)
@jit_integrand_function
def regular_function(x1, x2, x3, x4):
return x1 + x2 + x3 + x4
def do_integrate_wo_arrays(a, b, c, lolim=0, hilim=1):
return integrate.quad(regular_function, lolim, hilim, (a, b, c))
>>> print(do_integrate_wo_arrays(1,2,3,lolim=2, hilim=10))
(96.0, 1.0658141036401503e-12)
Этот код работаетпросто хорошо. Я в состоянии выполнить функцию jitnd и вернуть функцию jited как объект LowLevelCallable. Однако на самом деле мне нужно передать моему интегранту два numpy.array с, и описанные выше разрывы конструкции:
from scipy import integrate, LowLevelCallable
from numba import njit, cfunc
from numba.types import intc, float64, CPointer
def jit_integrand_function(integrand_function):
jitted_function = njit(integrand_function)
@cfunc(float64(intc, CPointer(float64)))
def wrapped(n, xx):
return jitted_function(xx[0], xx[1], xx[2], xx[3])
return LowLevelCallable(wrapped.ctypes)
@jit_integrand_function
def function_using_arrays(x1, x2, array1, array2):
res1 = np.interp(x1, array1[0], array1[1])
res2 = np.interp(x2, array2[0], array2[1])
return res1 + res2
def do_integrate_w_arrays(a, lolim=0, hilim=1):
foo = np.arange(20, dtype=np.float).reshape(2, -1)
bar = np.arange(60, dtype=np.float).reshape(2, -1)
return integrate.quad(function_using_arrays, lolim, hilim, (a, foo, bar))
>>> print(do_integrate_w_arrays(3, lolim=2, hilim=10))
Traceback (most recent call last):
File "C:\ProgramData\Miniconda3\lib\site-packages\IPython\core\interactiveshell.py", line 3267, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "<ipython-input-63-69c0074d4936>", line 1, in <module>
runfile('C:/Users/mosegui/Desktop/fos4x_pkg_develop/python-packages/fos4x_tec/fos4x_tec/test_scipy_numba.py', wdir='C:/Users/mosegui/Desktop/fos4x_pkg_develop/python-packages/fos4x_tec/fos4x_tec')
File "C:\Program Files\JetBrains\PyCharm Community Edition 2018.3.4\helpers\pydev\_pydev_bundle\pydev_umd.py", line 197, in runfile
pydev_imports.execfile(filename, global_vars, local_vars) # execute the script
File "C:\Program Files\JetBrains\PyCharm Community Edition 2018.3.4\helpers\pydev\_pydev_imps\_pydev_execfile.py", line 18, in execfile
exec(compile(contents+"\n", file, 'exec'), glob, loc)
File "C:/Users/mosegui/Desktop/fos4x_pkg_develop/python-packages/fos4x_tec/fos4x_tec/test_scipy_numba.py", line 29, in <module>
@jit_integrand_function
File "C:/Users/mosegui/Desktop/fos4x_pkg_develop/python-packages/fos4x_tec/fos4x_tec/test_scipy_numba.py", line 13, in jit_integrand_function
@cfunc(float64(intc, CPointer(float64)))
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\decorators.py", line 260, in wrapper
res.compile()
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler_lock.py", line 32, in _acquire_compile_lock
return func(*args, **kwargs)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\ccallback.py", line 69, in compile
cres = self._compile_uncached()
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\ccallback.py", line 82, in _compile_uncached
cres = self._compiler.compile(sig.args, sig.return_type)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\dispatcher.py", line 81, in compile
raise retval
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\dispatcher.py", line 91, in _compile_cached
retval = self._compile_core(args, return_type)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\dispatcher.py", line 109, in _compile_core
pipeline_class=self.pipeline_class)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler.py", line 528, in compile_extra
return pipeline.compile_extra(func)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler.py", line 326, in compile_extra
return self._compile_bytecode()
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler.py", line 385, in _compile_bytecode
return self._compile_core()
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler.py", line 365, in _compile_core
raise e
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler.py", line 356, in _compile_core
pm.run(self.state)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler_machinery.py", line 328, in run
raise patched_exception
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler_machinery.py", line 319, in run
self._runPass(idx, pass_inst, state)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler_lock.py", line 32, in _acquire_compile_lock
return func(*args, **kwargs)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler_machinery.py", line 281, in _runPass
mutated |= check(pss.run_pass, internal_state)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\compiler_machinery.py", line 268, in check
mangled = func(compiler_state)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\typed_passes.py", line 94, in run_pass
state.locals)
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\typed_passes.py", line 66, in type_inference_stage
infer.propagate()
File "C:\Users\mosegui\AppData\Roaming\Python\Python36\site-packages\numba\typeinfer.py", line 951, in propagate
raise errors[0]
numba.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Failed in nopython mode pipeline (step: nopython frontend)
Invalid use of Function(<built-in function getitem>) with argument(s) of type(s): (float64, Literal[int](0))
* parameterized
In definition 0:
All templates rejected with literals.
In definition 1:
All templates rejected without literals.
In definition 2:
All templates rejected with literals.
In definition 3:
All templates rejected without literals.
In definition 4:
All templates rejected with literals.
In definition 5:
All templates rejected without literals.
In definition 6:
All templates rejected with literals.
In definition 7:
All templates rejected without literals.
In definition 8:
All templates rejected with literals.
In definition 9:
All templates rejected without literals.
In definition 10:
All templates rejected with literals.
In definition 11:
All templates rejected without literals.
In definition 12:
All templates rejected with literals.
In definition 13:
All templates rejected without literals.
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of intrinsic-call at C:/Users/mosegui/Desktop/fos4x_pkg_develop/python-packages/fos4x_tec/fos4x_tec/test_scipy_numba.py (32)
[2] During: typing of static-get-item at C:/Users/mosegui/Desktop/fos4x_pkg_develop/python-packages/fos4x_tec/fos4x_tec/test_scipy_numba.py (32)
File "test_scipy_numba.py", line 32:
def diff_moment_edge(radius, alpha, chord_df, aerodyn_df):
<source elided>
# # calculate blade twist for radius
# sensor_twist = np.arctan((2 * rated_wind_speed) / (3 * rated_rotor_speed * (sensor_radius / 30.0) * radius)) * (180.0 / np.pi)
^
[1] During: resolving callee type: type(CPUDispatcher(<function function_using_arrays at 0x0000020C811827B8>))
[2] During: typing of call at C:/Users/mosegui/Desktop/fos4x_pkg_develop/python-packages/fos4x_tec/fos4x_tec/test_scipy_numba.py (15)
[3] During: resolving callee type: type(CPUDispatcher(<function function_using_arrays at 0x0000020C811827B8>))
[4] During: typing of call at C:/Users/mosegui/Desktop/fos4x_pkg_develop/python-packages/fos4x_tec/fos4x_tec/test_scipy_numba.py (15)
File "test_scipy_numba.py", line 15:
def jit_integrand_function(integrand_function):
<source elided>
jitted_function = njit(integrand_function)
^
Ну, теперь, очевидно, это не работает, потому что в этом примере я не модифицировалдизайн декоратора. Но это как раз и есть суть моего вопроса: я не до конца понимаю эту ситуацию и поэтому не знаю, как изменить аргументы cfunc для передачи массива в качестве параметра и при этом соответствовать требованиям scipy.integrate.quad для подписи. В документации numba, которая представляет CPointers, есть пример того, как передать массив в numba.cfunc:
ABI собственной платформы, используемые C или C ++, не имеютПонятие о форме массива, как в Numpy. Одним из распространенных решений является передача необработанного указателя данных и одного или нескольких аргументов размера (в зависимости от размерности). Numba должен предоставить способ перестроить представление массива этих данных внутри обратного вызова.
from numba import cfunc, carray
from numba.types import float64, CPointer, void, intp
# A callback with the C signature `void(double *, double *, size_t)`
@cfunc(void(CPointer(float64), CPointer(float64), intp))
def invert(in_ptr, out_ptr, n):
in_ = carray(in_ptr, (n,))
out = carray(out_ptr, (n,))
for i in range(n):
out[i] = 1 / in_[i] ```
Я так или иначе понимаю, что CPointer используется для построения массива в C, как в сигнатуремоего примера декоратора CPointer(float64) собирает все числа с плавающей точкой, переданные в качестве аргументов, и помещает их в массив. Тем не менее, я все еще не могу собрать его и посмотреть, как я могу использовать его для передачи массива, а не для создания массива из набора float аргументов, которые я передаю.
РЕДАКТИРОВАТЬ:
Ответ по @max9111 сработал в том смысле, что смог перейти к scipy.integrate.quad a LowLevelCallable, что улучшило временную эффективность вычислений. Это очень ценно, так как теперь стало намного понятнее, как управление адресами памяти работает в C. Несмотря на то, что концепция структурированного массива не существует в native C, я могу создать структурированный массив в Python с данными, которые C собирается хранить. в непрерывной области памяти и доступ к ней / указать на него через уникальный адрес памяти. Отображение, предоставляемое структурированным массивом, позволяет идентифицировать различные компоненты этой области памяти.
Несмотря на то, что решение @max9111 работает и решает вопрос, который я опубликовал изначально, с точки зрения Python этот подход вносит определенные накладные расходы, которые при определенных условиях могут занимать больше времени, чем время, полученное при вызове * 1067. * функция интеграции через LowLevelCallable.
В моем реальном случае я использую интеграцию как шаг задачи двумерной оптимизации. Каждый шаг оптимизации требует двойной интеграции, а интеграция требует девяти скалярных параметров и двух массивов. Пока я не смог решить интеграцию с помощью LowLevelCallable, единственное, что я мог сделать для ускорения кода, было просто njit функция интегрирования. И это работало прилично, хотя интеграция все еще выполнялась через API Python.
В моем случае реализация решения @max9111 значительно повысила эффективность во времени интеграции (примерно с 0,0009 с за шагпримерно до 0,0005). Тем не менее, шаг создания структурированного массива, C-распаковка данных, передача их объединенному подынтегральному выражению и возврат LowLevelCallable добавили дополнительно 0,3 с в среднем на одну итерацию, что ухудшило мою ситуацию.
Здесь есть некоторый игрушечный код для демонстрации того, как подход LowLevelCallable становится худшим вариантом, чем больше он подходит для итеративного процесса:
import ctypes
import timeit
from tqdm import tqdm
import numpy as np
from scipy import integrate, LowLevelCallable
import numba as nb
from numba import types
import matplotlib.pyplot as plt
##################################################
# creating some sample data and parameters
a = 3
foo = np.arange(200, dtype=np.float64).reshape(2, -1)
bar = np.arange(600, dtype=np.float64).reshape(2, -1)
lim1 = 0
lim2 = 1
@nb.njit
def function_using_arrays(x1, x2, array1, array2):
res1 = np.interp(x1, array1[0], array1[1])
res2 = np.interp(x2, array2[0], array2[1])
return res1 + res2
##################################################
# JIT INTEGRAND
def do_integrate_w_arrays_jit(a, array1, array2, lolim=0, hilim=1):
return integrate.quad(function_using_arrays, lolim, hilim, (a, array1, array2))
def process_jit_integrand():
do_integrate_w_arrays_jit(a, foo, bar, lolim=lim1, hilim=lim2)
##################################################
# LOWLEV CALLABLE
def create_jit_integrand_function(integrand_function,args,args_dtype):
@nb.cfunc(types.float64(types.float64,types.CPointer(args_dtype)))
def wrapped(x1,user_data_p):
#Array of structs
user_data = nb.carray(user_data_p, 1)
#Extract the data
x2=user_data[0].a
array1=user_data[0].foo
array2=user_data[0].bar
return integrand_function(x1, x2, array1, array2)
return wrapped
def do_integrate_w_arrays_lowlev(func,args,lolim=0, hilim=1):
integrand_func = LowLevelCallable(func.ctypes,user_data=args.ctypes.data_as(ctypes.c_void_p))
return integrate.quad(integrand_func, lolim, hilim)
def process_lowlev_callable():
args_dtype = types.Record.make_c_struct([
('a', types.float64),
('foo', types.NestedArray(dtype=types.float64, shape=foo.shape)),
('bar', types.NestedArray(dtype=types.float64, shape=bar.shape)),])
args=np.array((a, foo, bar), dtype=args_dtype)
func = create_jit_integrand_function(function_using_arrays,args,args_dtype)
do_integrate_w_arrays_lowlev(func, args, lolim=0, hilim=1)
##################################################
repetitions = range(100)
jit_integrand_delays = [timeit.timeit(stmt=process_jit_integrand,
number=repetition) for repetition in tqdm(repetitions)]
lowlev_callable_delays = [timeit.timeit(stmt=process_lowlev_callable,
number=repetition) for repetition in tqdm(repetitions)]
fig, ax = plt.subplots()
ax.plot(repetitions, jit_integrand_delays, label="jit_integrand")
ax.plot(repetitions, lowlev_callable_delays, label="lowlev_callable")
ax.set_xlabel('number of repetitions')
ax.set_ylabel('calculation time (s)')
ax.set_title("Comparison calculation time")
plt.tight_layout()
plt.legend()
plt.savefig(f'calculation_time_comparison_{repetitions[-1]}_reps.png')
Здесь два варианта (только сравнение с интегралом против @max9111Решение) сравниваются. В модифицированной версии решения @max9111 я навсегда связал функцию интегрирования (function_using_arrays) и удалил этот шаг из create_jit_integrand_function, что сокращает время «накладных расходов» на 20%. Более того, ради скорости я также отключил функцию jit_with_dummy_data и включил ее функциональность в тело process_lowlev_callable, в основном для избежания ненужного вызова функции. Ниже найдите время расчета для обоих решений для серии до 100 циклов:
Как видите, если вы находитесь витеративный процесс, время, сэкономленное в каждом отдельном вычислении (30+% !!), не окупает накладные расходы, переносимые парой дополнительных функций, которые необходимо реализовать для построения LowLevelCallable (функции, которые также итеративно вызываются и выполняются черезPython C API).
Итог: это решение очень хорошо для сокращения времени вычисления в одном очень тяжелом интеграле, но простое совпадение с подынтегральным выражением кажется более выгодным при решении средних интегралов в итерационном процессе, посколькудополнительные функции, необходимые для LowlevelCallable, которые нужно вызывать так же часто, как и сама интеграция, берут свое.
В любом случае, большое спасибо. Хотя это решение не сработает для меня, я узнал ценные вещи и считаю, что мой вопрос решен.
РЕДАКТИРОВАТЬ 2:
Я неправильно понял части @max9111 'Решение и роль, которую сыграла функция create_jit_integrand_function, и я ошибочно компилировал LowLevelCallable на каждом шаге своей оптимизации (чего мне не нужно делать, потому что параметры и массивы, передаваемые в интегральное изменение, каждыйитерации, их форма и, следовательно, структура C остается неизменной).
Реорганизованная версия кода из приведенного выше EDIT , которая имеет смысл:
import ctypes
import timeit
from tqdm import tqdm
import numpy as np
from scipy import integrate, LowLevelCallable
import numba as nb
from numba import types
import matplotlib.pyplot as plt
##################################################
# creating some sample data and parameters
a = 3
foo = np.arange(200, dtype=np.float64).reshape(2, -1)
bar = np.arange(600, dtype=np.float64).reshape(2, -1)
lim1 = 0
lim2 = 1
def function_using_arrays(x1, x2, array1, array2):
res1 = np.interp(x1, array1[0], array1[1])
res2 = np.interp(x2, array2[0], array2[1])
return res1 + res2
##################################################
# JIT INTEGRAND
def do_integrate_w_arrays_jit(a, array1, array2, lolim=0, hilim=1):
return integrate.quad(nb.njit(function_using_arrays), lolim, hilim, (a, array1, array2))
def process_jit_integrand():
do_integrate_w_arrays_jit(a, foo, bar, lolim=lim1, hilim=lim2)
##################################################
# LOWLEV CALLABLE
def create_jit_integrand_function(integrand_function, args_dtype):
jitted_function = nb.njit(integrand_function)
@nb.cfunc(types.float64(types.float64,types.CPointer(args_dtype)))
def wrapped(x1,user_data_p):
#Array of structs
user_data = nb.carray(user_data_p, 1)
#Extract the data
x2=user_data[0].a
array1=user_data[0].foo
array2=user_data[0].bar
return jitted_function(x1, x2, array1, array2)
return wrapped
def do_integrate_w_arrays_lowlev(func,args,lolim=0, hilim=1):
integrand_func=LowLevelCallable(func.ctypes,user_data=args.ctypes.data_as(ctypes.c_void_p))
return integrate.quad(integrand_func, lolim, hilim)
def process_lowlev_callable():
do_integrate_w_arrays_lowlev(func, np.array((a, foo, bar), dtype=args_dtype), lolim=0, hilim=1)
##################################################
repetitions = range(100)
jit_integrand_delays = [timeit.timeit(stmt=process_jit_integrand, number=repetition) for repetition in tqdm(repetitions)]
args_dtype = types.Record.make_c_struct([
('a', types.float64),
('foo', types.NestedArray(dtype=types.float64, shape=foo.shape)),
('bar', types.NestedArray(dtype=types.float64, shape=bar.shape)),])
func = create_jit_integrand_function(function_using_arrays, args_dtype)
lowlev_callable_delays = [timeit.timeit(stmt=process_lowlev_callable, number=repetition) for repetition in tqdm(repetitions)]
fig, ax = plt.subplots()
ax.plot(repetitions, jit_integrand_delays, label="jit_integrand")
ax.plot(repetitions, lowlev_callable_delays, label="lowlev_callable")
ax.set_xlabel('number of repetitions')
ax.set_ylabel('calculation time (s)')
ax.set_title("Comparison calculation time")
plt.tight_layout()
plt.legend()
plt.savefig(f'calculation_time_comparison_{repetitions[-1]}_reps.png')
В этой конфигурации построение LowLevelCallable (которое действительно стоит немного времени) должно быть выполнено только один раз, и весь процесс на несколько порядков быстрее:
и крупный план для lowlev_callable:
