Is there a way to use python's LowLevelCallable with interpolation function?

user_16528170 0
Viewed 180

I'm dealing with an integration involving lots of interpolation functions, e.g.

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)

where func_a, func_g1, func_g1p, func_g2, func_g2p and func_jl are all interpolation functions (done with cubic spline). Unfortunately the code performance is not satisfactory, so I'm thinking of accelerating it with scipy's LowLevelCallable.

From scipy's documents it seems that I should include in the C code just a function body of my integrand with certain function signature pattern, e.g.

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));
}

where x is the argument to integrate upon, and userdata contains some other parameters. But this means that there's no place for me to initialize my interpolation functions, unless I provide the table in userdata and generate the interpolation function every time I call the integrand (the performance won't be good I guess). So I'm wondering if there's any other ways to either enabling me to practically use the scipy.LowLevelCallable approach with interpolating functions, or accelerate the integration in another way.

Thanks!

scipy c python
1 Answers

Yes, but you have to rewrite the interpolating function

You mentioned, that you have cubic splines. Without considering special cases this can quite easily be implemented in C-code. The generation of the splines can be done in Python.

Example for the low-level callable

#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);
}

How to use the low-level callable from Python?

This can be accomplished using a ctypes-wrapper including some functions to generate the struct which can be passed using the void-pointer.

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)

Comparison of the implementations

Now compare the runtime of a pure Python implementation with the C-implementation.

Pure 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)

Using the low-level-callable

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)
user_40457740
Related