有没有一种方法可以使用python的LovelCallable with interpolation函数？

Question

我正在处理一个涉及很多插值函数的积分，例如

import scipy.integrate as integrate

integrate.quad(lambda t: 1/func_a(t)*(func_g2(tp)*func_g1p(t)-func_g1(tp)*func_g2p(t))
               *func_jl(t)/((tau0-t)**2 * klist[84]**2), tp, tau0, limit=100)

其中func_a，func_g1，func_g1p，func_g2，func_g2p和func_jl都是插值函数（使用三次样条曲线完成）。不幸的是，代码性能并不令人满意，所以我正在考虑使用scipy的LowLevelCallable来加速它

从scipy的文档来看，似乎我应该在C代码中只包含我的被积函数的一个函数体，具有特定的函数签名模式，例如

double f(double * x, void * userdata)
{
// integrand here, e.g.
// return 1/func_a(t)*(func_g2(tp)*func_g1p(t)-func_g1(tp)*func_g2p(t))
//               *func_jl(t)/(pow(tau0-t, 2) * pow(klist[84], 2));
}

其中，x是要进行积分的参数，userdata包含一些其他参数。但这意味着我没有地方初始化插值函数，除非我在userdata中提供表，并在每次调用被积函数时生成插值函数（我想性能不会很好）。因此，我想知道是否还有其他方法可以让我实际使用scipy.lowlevellable方法和插值函数，或者以另一种方式加速集成

谢谢

Answer 1

是的，但是你必须重写插值函数

你提到过，你有三次样条曲线。不考虑特殊情况，这可以很容易地用C代码实现。样条曲线的生成可以在Python中完成

低级可调用的示例

#ifdef _WIN32
#    define API __declspec(dllexport)
#else
#  define API
#endif

#include <math.h>
#include <stdint.h>
#include <stddef.h>
#include <stdio.h>

/*
intervalls have to be increasing, no checks are performed.
*/

double eval_spline(double xp,size_t k,size_t m,double *c,double *x){
    size_t interv=0;
    
    for (size_t i=1;i<m;i++){
        if (xp>x[i]){interv+=1;}
    }
    
    k=k-1;
    double acc=0;
    
    for (size_t i=0;i<k+1;i++){
        acc+=c[i*m+interv]*pow(xp-x[interv],k-i);
    }
    return acc;
}

struct data{
   double *cs_1_c;
   double *cs_1_x;
   size_t  cs_1_k;
   size_t  cs_1_m;
   double *cs_2_c;
   double *cs_2_x;
   size_t  cs_2_k;
   size_t  cs_2_m;
};

API double Integrand(double t, void *userdata){
    struct data input=*(struct data*)userdata;
    
    return eval_spline(t,input.cs_1_k,input.cs_1_m,input.cs_1_c,input.cs_1_x)
      *pow(eval_spline(t+0.5,input.cs_2_k,input.cs_2_m,input.cs_2_c,input.cs_2_x),2);
}

如何使用Python中的低级可调用函数？

这可以通过使用ctypes包装器来实现，该包装器包括一些生成结构的函数，该结构可以使用void指针来传递

import ctypes
import os
from scipy import integrate, LowLevelCallable
import numpy as np

lib = ctypes.cdll.LoadLibrary("integrand.dll")

def get_integrand():
    lib.Integrand.restype =  ctypes.c_double
    lib.Integrand.argtypes = (ctypes.c_double, ctypes.c_void_p)
    return lib.Integrand

dble_p=ctypes.POINTER(ctypes.c_double)
class args_struct(ctypes.Structure):
    _fields_ = [('cs_1_c',   dble_p),
                ('cs_1_x',   dble_p),
                ('cs_1_k',   ctypes.c_size_t),
                ('cs_1_m',   ctypes.c_size_t),
                ('cs_2_c',   dble_p),
                ('cs_2_x',   dble_p),
                ('cs_2_k',   ctypes.c_size_t),
                ('cs_2_m',   ctypes.c_size_t)]

def create_struct(cs_1,cs_2):
    #The arrays must be c-contigous, if not make them c-contiguous
    cs_1_c = np.ascontiguousarray(cs_1.c)
    cs_1_x = np.ascontiguousarray(cs_1.x)
    cs_2_c = np.ascontiguousarray(cs_2.c)
    cs_2_x = np.ascontiguousarray(cs_2.x)
    
    args=args_struct(ctypes.cast(cs_1_c.ctypes.data,dble_p),
                     ctypes.cast(cs_1_x.ctypes.data,dble_p),
                     cs_1_c.shape[0],
                     cs_1_c.shape[1],
                     ctypes.cast(cs_2_c.ctypes.data,dble_p),
                     ctypes.cast(cs_2_x.ctypes.data,dble_p),
                     cs_2_c.shape[0],
                     cs_2_c.shape[1],)
    return ctypes.cast(ctypes.pointer(args),ctypes.c_void_p)

def do_integrate_w_arrays(func,args,lolim=0, hilim=1):
    integrand_func=LowLevelCallable(func,user_data=args)
    return integrate.quad(integrand_func, lolim, hilim)

实现的比较

现在比较纯Python实现和C实现的运行时

纯Python

import numpy as np
from scipy.interpolate import CubicSpline
import scipy.integrate as integrate

x = np.arange(10)
y1 = np.sin(x)
y2 = np.cos(x)**2

cs_1 = CubicSpline(x, y1)
cs_2 = CubicSpline(x, y2)


def Integrand_orig(t):
    return cs_1(t)*cs_2(t+0.5)**2

%timeit res=integrate.quad(Integrand_orig,a=-0.5,b=10.5)
#15.2 ms ± 237 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

使用低级可调用

import wrapper

func=wrapper.get_integrand()

args=wrapper.create_struct(cs_1,cs_2)
wrapper.do_integrate_w_arrays(func,args,lolim=-0.5, hilim=10.5)
#136 µs ± 445 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)