在 C 中使用函数时将参数传递给 scipy.LowLevelCallable答案

【问题标题】：Passing arguments to scipy.LowLevelCallable while using functions in C在 C 中使用函数时将参数传递给 scipy.LowLevelCallable
【发布时间】：2020-08-19 17:05:52
【问题描述】：

我正在尝试使用 C 定义的函数在 SciPy 中进行数值积分。给出的示例案例 here (SciPy documentation) 工作正常。

就我而言，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    */
}

并且，下面是调用编译后得到的testlib.so文件的test.py函数，

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 个这样的参数，在调用 py 文件中以 np.float64 的形式提供。

我想知道是否可以将参数作为数组作为 user_data 传递给函数 f2 。

从documentation到nquad，发现参数要作为数组传递，C函数中int n是传递的参数个数。
另外，我愿意尝试其他选项，例如 cython、pyCapsule，但没有经验。使用 numba 和 jit 发现非常相似的问题，没有传递其他参数。 Using numba and jit for integration: SE

用于编译testlib.c：$ gcc -shared -fPIC -o testlib.so testlib.c

【问题讨论】：

顺便说一句，您可能应该编译 testlib.so 并优化（即 -O2 或 -O3 或任何最佳）以获得最佳性能。
只要你只传递标量 dtype np.float64 这里有一个例子stackoverflow.com/a/60619681/4045774，它展示了如何使用 Numba 做到这一点。您可以像在纯 Python 中那样将附加参数作为 args 传递。

标签： python scipy cython numerical-integration

【解决方案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-array 中（正如我最初理解你的问题）：

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)

但仍需要检查形状/步幅。

【讨论】：