在使用C语言的函数时向scipy.LowLevelCallable传递参数。

问题描述 投票:1回答:1

我试图在SciPy中使用C定义的函数进行数值积分。所给出的例子是 在这里(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

python scipy cython numerical-integration
1个回答
1
投票

它们是不止一种剥皮的方式,但如果你使用的是 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-array的内存是否连续,关于如何获得这些信息的例子,可以在以下地方查找到 numpy.ctypeslib.as_ctypes.

也许一个不那么低级的获取指针的方式是

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

但仍然需要对形状特征进行检查。

© www.soinside.com 2019 - 2024. All rights reserved.