I'm quite desperate trying to implement an integration which result I wanna use as initial data for a PDE. The integral itself is scalar-valued and needs to be evaluated on the whole space $\mathbb{R}^3$. First of all, I do not expect a very fast super easy integration but hope that it can be done in a reasonable amount of time. Therefore I decided to use nquad in connection the LowLvelCallable option. I'm trying to do the implementation of my code in Python.
The LowLevelCallable option of scipy calls cfuncs with less overhead to save time in each individual integration step. Now I found a sample implementation for my purpose where quad is used instead of nquad, Sample. As far as I understood the idea of this sample is it to use a decorator to "trick" the signature of quad which demands an input of a cfunc in the form of sig = (types.CPointer(types.double), types.CPointer(types.void)) to bypass the problem that there isn't a possibility to extract data from a viod type when the code is written in Python.
That works pretty well, but now I need to adapt this approach for nquad nothing should change in principle, the only difference is that nquad demands a signature of the form sig=(types.intc, types.CPointer(types.double), types.CPointer(types.void)). Unfortunately, it turns out that just to edit the signature of the cfunc doesn't do the trick, and therefore without any background in high-performance computing and in C this turned into a very nasty try and error game. My hope is that someone can help me out and saying me what my mistake is.
Here is the sample of my code
import numpy as np
import numba as nb
from numba import types
from scipy import integrate, LowLevelCallable
import ctypes
from math import sqrt, exp
#define some constants
c = 137
sigma = 10
p_boost = 100
m = 1
gamma = sqrt(1 + p_boost ** 2 / c ** 2)
beta_x = p_boost / (gamma * m * c)
x = np.linspace(-4, 4, 32, endpoint = False)
y = np.linspace(-4, 4, 32, endpoint = False)
z = np.linspace(-4, 4, 32, endpoint = False)
#define arg data type
args_dtype = types.Record.make_c_struct([
('x_coord', types.float64),
('y_coord', types.float64),
('z_coord', types.float64),])
def create_jit_integrand_function(integrand_function,args,args_dtype):
jitted_function = nb.njit(integrand_function)
@nb.cfunc(types.intc,types.float64(types.float64,types.CPointer(args_dtype)))
def wrapped(n,x1,user_data_p):
#Array of structs
user_data = nb.carray(user_data_p, 1)
#Extract arg data
x_co=user_data[0].x_co
y_co=user_data[0].y_co
z_co=user_data[0].z_co
return jitted_function(n,x1, x_co, y_co, z_co)
return wrapped
def function_using_arrays(n ,p, x_co, y_co, z_co):
p_square = 0.0
for i in range(n):
p_square += p[i]**2
rel_energy = sqrt(p_square + 1)
px_new = beta_x * gamma * rel_energy / c + gamma * p[0]
p_square_new = px_new**2 + p[1]**2 + p[2]**2
rel_energy_new = sqrt(p_square_new +1)
return exp(-p_square_new/(2*sigma**2) - 1j * (p[0]*x_co + p[1]*y_co + p[2]*z_co)) * (rel_energy_new - beta_x * p[boost_index] + 0j)**(1/2) * rel_energy_new ** (-1)
def jit_with_dummy_data(args,args_dtype):
func=create_jit_integrand_function(function_using_arrays,args,args_dtype)
return func
def do_integrate_w_arrays(func,args,ranges):
integrand_func=LowLevelCallable(func.ctypes,user_data=args.ctypes.data_as(ctypes.c_void_p))
return integrate.nquad(integrand_func, ranges)
output = np.zeros((len(x),len(y),len(z)), dtype=complex)
for i in range(len(x)):
for j in range(len(y)):
for k in range(len(z)):
#creating args
x_co = x[i]
y_co = y[j]
z_co = z[k]
#hand over args
args=np.array((x_co, y_co, z_co),dtype=args_dtype)
func=jit_with_dummy_data(args,args_dtype)
#perform integration
ranges = ((-np.inf, np.inf), (-np.inf, np.inf), (-np.inf, np.inf))
output[i,j,k] = do_integrate_w_arrays(func,args,ranges)[0]
Unfortunately, it turns out that just to edit the signature of the cfunc doesn't do the trick, … My hope is that someone can help me out and saying me what my mistake is.
The mistake was to insert the additional integer parameter before the function signature instead of into the parameter list; correct:
@nb.cfunc(types.float64(types.intc, types.float64, types.CPointer(args_dtype)))
This correction reveals a naming mismatch in wrapped; correct:
x_co=user_data[0].x_coord
y_co=user_data[0].y_coord
z_co=user_data[0].z_coord
That correction allows the program to proceed to the next error:
NameError: name 'boost_index' is not defined
boost_index is used at the end of function_using_arrays; you have to see where this value shall come from.