Передача аргументов в scipy.LowLevelCallable при использовании функций в C - PullRequest
Новая
       33

Передача аргументов в scipy.LowLevelCallable при использовании функций в C

1 голос
/ 05 мая 2020

Я пытаюсь использовать C определенную функцию для численного интегрирования в SciPy. Пример приведенного здесь (документация SciPy) отлично работает.

В моем случае файл testlib.c - это 

/* testlib.c */

#include <math.h>
#define PI 3.14159265358979323846

double factor(double phi, double r) {
    double val = (2*PI)/(pow(r, 2) + 3*cos(phi));
    return val;
}

//------------------------------------

double f2(int n, double *x, void *user_data) {
    double c = *(double *)user_data;
    double v1 = factor(c, 0.25); // value of phi defined inline but it is an argument
    return v1 + x[0] - x[1] ; /* corresponds to v1 + x - y    */
}

И функция test.py вызов файла testlib.so, полученного после компиляции, приведен ниже: 

import os, ctypes
from scipy import integrate, LowLevelCallable

lib = ctypes.CDLL(os.path.abspath('testlib.so'))

# define return type in .restype
lib.f2.restype = ctypes.c_double

# define argument type in .argtypes
lib.f2.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_void_p)

# additional argument, here a constant, casting needed
c = ctypes.c_double(1.0)
user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)

# pass extra argument
func = LowLevelCallable(lib.f2, user_data)

# quadrature in 3-dimensions
out=integrate.nquad(func, [[0, 10], [-10, 0]])

print(out)

# -----------------------------------------------
phi = 0.25   #  How to pass these to the C function 
esv = 1.25   
cv1 = 0.03   
cv2 = -0.15  
cv3 = 3.0

Мой вопрос : Как передать дополнительные аргументы, такие как c, в функцию f2. В моем случае у меня есть 5 таких аргументов, доступных как np.float64 в вызывающем файле py.

Интересно, могу ли я передать аргументы в виде массива user_data функции f2.

  • Из документации для nquad обнаружено, что аргументы должны передаваться как массив, а int n в функции C - это количество переданных аргументов .

  • Кроме того, я открыт, чтобы попробовать другие варианты, такие как cython, pyCapsule, но без опыта в них. Обнаружен очень похожий вопрос с использованием numba и jit, где не передаются дополнительные аргументы. Использование numba и jit для интеграции: SE

Для компиляции testlib.c: $ gcc -shared -fPIC -o testlib.so testlib.c

1 Ответ

1 голос
/ 05 мая 2020

Это несколько способов снять шкуру с кошки, но если вы используете ctypes, есть возможность придерживаться ctypes, например:

Вы можете создать массив и инициализировать его с помощью значения, например: 

ptr_to_buffer=(ctypes.c_double*5)(phi,esv,cv1,cv2,cv3)
user_data = ctypes.cast(ptr_to_buffer, ctypes.c_void_p)

или если данные уже находятся в массиве numpy (как я изначально понял ваш вопрос): 

import numpy as np
a=np.array([1.0, 2.0, 3.0, 4.0, 5.0], np.float64)

import ctypes
ptr_to_buffer=(ctypes.c_double*5).from_buffer(a)
user_data = ctypes.cast(ptr_to_buffer, ctypes.c_void_p)

в данном случае user_data является копией a и не разделяет память, что иногда хорошо, а иногда нет.

Для больших массивов можно также разрешить user_data также делить память с numpy - array: 

user_data = ctypes.c_void_p(a.__array_interface__['data'][0])

, что можно проверить с помощью: 

ctypes.cast(user_data, ctypes.POINTER(ctypes.c_double))[0] = 42.0
print(a[0])
# 42.0 and not 1.0

Для этого варианта вам действительно нужно проверить, что память numpy -массив является непрерывной, для Например, как получить эту информацию, можно посмотреть, например, в numpy.ctypeslib.as_ctypes.

Возможно, менее низкоуровневым способом получения указателя является 

user_data =  ctypes.cast(np.ctypeslib.as_ctypes(a), ctypes.c_void_p)

но все же необходима проверка формы / шагов.

...