I have written function with @cfunc wrapper (called integrand). This is for passing the function as C type callable for scipy.integrate.nquad integration routine. Right now all the arguments are doubles, but I would like to add a function as one of the arguments. Is there a way to do this? Or is there a alternative that would give me the same speed improvements that I got by switching to C type using numba?
hi @MadFisa could you share an example with code? Showing what you’d like to do. There’s a FunctionType that probably does what you want, but it’ll be easier to test in your own code.
Luk
My current code looks something like this (Which works).
from scipy import integrate,LowLevelCallable
from numba import cfunc, types
#Follows the signature as required per scipy documentation.
c_sig = types.double(types.intc, types.CPointer(types.double))
@cfunc(c_sig)
def integrand(n,args,pyFunc):
df = args[1]*args[2]*args[3]
return df
ctype_integrand = LowLevelCallable(integrand.ctypes)
#integrate it for some limits lim_1,lim_2.....
integrate.nquad(ctype_integrand, ((lim_1,lim_2),(lim_3,lim_4),(lim_5,lim_6)))
I want it to change it so that it takes an additional argument which is a function.
Here is an example of what I want to do. I want the function to be compatible with scipy quad as per the following link Integration (scipy.integrate) — SciPy v1.7.0 Manual. I am pretty sure I should be able to use *user_data to pass a function as an additional argument but not sure how.
from scipy import integrate,LowLevelCallable
from numba import cfunc, types
c_sig = types.double(types.intc, types.CPointer(types.double), types.someTypeForFunction)
@cfunc(c_sig)
def integrand(n,args,pyFunc):
xyz = args[1]*args[2]*args[3]
df = xyz * pyFunc(args[1],args[2])
return df
ctype_integrand = LowLevelCallable(integrand.ctypes)
integrate.nquad(ctype_integrand, ((lim_1,lim_2),(lim_3,lim_4),(lim_5,lim_6)))
I dont know what I should put as signature for pyFunc. Scipy documentation calls it as void *user_data which could be used to pass more data to the function. I want it to store the function that I want to pass as an argument
I have something for you but I cannot get it to work all the way.
from scipy import integrate,LowLevelCallable
from numba import cfunc, types
c_sig_foo =types.double(types.double, types.double)
@cfunc(c_sig_foo)
def foo(a, b):
return a+b
c_sig2 = types.double(types.intc, types.CPointer(types.double), types.FunctionType(types.double(types.double, types.double)))
@cfunc(c_sig2)
def integrand2(n, args, func):
df = args[1]*args[2]*args[3]
return df + func(1, 2)
print(integrand2(1, (1, 2, 3, 4), foo))
This works as expected. However, scipy requires the third argument to be void *user_data. That means that c_sig2 should be c_sig2 = types.double(types.intc, types.CPointer(types.double), types.voidptr).
What I don’t know is how to get that pointer and turn it into a function
@cfunc(c_sig2)
def integrand2(n, args, user_data):
df = args[1]*args[2]*args[3]
pyFunc = ....magic should happen here ...
return df + pyFunc(1, 2)
if user_data were an array then I would use the function carray. I suspect something similar is needed for functions but it might not exist.
I have spent a lot of time working with higher-order functions in Numba. I realized that many times higher-order functions are not needed. You could do
@cfunc(c_sig2)
def integrand_foo(n, args):
df = args[1]*args[2]*args[3]
return df + foo(1, 2)
@cfunc(c_sig2)
def integrand_bar(n, args):
df = args[1]*args[2]*args[3]
return df + bar(1, 2)
and you will probably get better performance.
For this kind of problems I usually make factory functions of the type
def make_integrand(pyFunc):
@cfunc(c_sig2)
def integrand(n, args):
df = args[1]*args[2]*args[3]
return df + pyFunc(1, 2)
return integrand
and then you do
# first function
integrand = make_integrand(foo)
ctype_integrand = LowLevelCallable(integrand.ctypes)
integrate.nquad(ctype_integrand, ((lim_1,lim_2),(lim_3,lim_4),(lim_5,lim_6)))
# another function
integrand = make_integrand(bar)
ctype_integrand = LowLevelCallable(integrand.ctypes)
integrate.nquad(ctype_integrand, ((lim_1,lim_2),(lim_3,lim_4),(lim_5,lim_6)))
You can even put them in a loop working that way.
Hope this helps,
Luk
@luk-f-a ,
You gave me some ideas that I could probably use. In particular, I think the make integrand (factory function) method should work for me. I had an idea similar to that in mind couldn’t quite put it into concrete form. Still it would be interesting to figure out a way to use the third null pointer mentioned in scipy documentation. I will update here if I find anything interesting. Thanks a lot.
@luk-f-a
Got the code to work using the function factory method. Here is my code.
from scipy import integrate,LowLevelCallable
from numba import cfunc, types
import numpy as np
c_sig = types.double(types.intc, types.CPointer(types.double))
def make_integrand(source_function):
@cfunc(c_sig)
def pyXspec_integrand(n,args):
theta = np.sqrt(2*299792458*args[2]/((1+args[6])*args[3]*3.086e+16))
x = (2*np.pi/(6.626e-34*299792458))*(1+args[6])*theta*args[1]*1.6022e-16*args[0]*1e-6
tau_a = np.power(args[1]*(1+args[6]),-args[4])*np.power(args[0]/0.1,4.-args[5])
temp = (np.sin(x)/x**2) - (np.cos(x)/x)
tau_theta = np.power(temp,2)
tau = tau_a*tau_theta
S = source_function(args[1]) # Where I use the source_function
dF = S*tau/args[2] # in erg/(kev cm^2 s)
dF *= 6.242e8 # to photons/(kev cm^2 s)
return dF
return pyXspec_integrand
c_sig_powerlaw = types.double(types.double)
@cfunc(c_sig_powerlaw)
def powerlaw(E):
return 1e-6*np.power(E, -1)
integrand = make_integrand(powerlaw)
ctype_integrand = LowLevelCallable(integrand.ctypes)
a_m = 0.025
a_M = 0.2
E_m = 0.3
E_M = 10.
t_m = 1e3
t_M = 8e3
s = 2
q = 4.5
z = 0
flux = integrate.nquad(ctype_integrand,((a_m,a_M), (E_m,E_M,)), args = (s,q,z))
Any general comments to improve the speed is welcome 
.
I would be calling the intgerate.nquad in a loop( for i in range(n) ). Thank you for your help.
you’re welcome! Will you call make_integrand inside the loop? How large is n?
No. Make_integrand will only be called once. This is the outline of how it would look like.
integrand = make_integrand(powerlaw)
ctype_integrand = LowLevelCallable(integrand.ctypes)
def dust_PL (engs,params,flux,flux_err,spec_num):
a_m, a_M, R_pc, s, q, tau0, z, S0, beta , nm= params
n = len(engs)
x1,x2=xp.AllData(spec_num).xflt #xp is a package
t_l=x1[1]
t_u=x2[1]
print(t_l)
for i in range(n - 1):
flux[i] = 0.1*tau0*integrate.nquad(ctype_integrand,
((a_m,a_M), (engs[i],engs[i+1]),(t_l,t_u)),
args = (R_pc, s, q, z, S0, beta))[0]
n would have size of around 400 typically. dust_PL would be given to a bayesian parameter estimator (a fitting routine), so I won’t know how many times it would be called beforehand. I expect it to be called 1000 - 10000 times. Let’s assume 5000 times for now.
This is a example using the *void user_data argument to pass numpy arrays as arguments Example. But I don’t know how to create a function inside the cfunc. The important part is the intrinsic function address_as_void_pointer to create a pointer from a int64 memory address, which could be passed to the cfunc using a structured numpy array.
Maybe that also helps when a function pointer is needed…
Can you explain what is going inside the function address_as_void_pointer? I don’t understand what is going on inside ( in particular codegen function).
From what I understand, this what the overall code is doing:
take a function with additional arguments as arrays to integrate (function_using_arrays) and pass it to create_jit_integrand_function to make a c type callable. wrapped function inside it acts kind of like a translator that takes the void pointer and find the arrays to pass to fucntion_using_arrays. This “translation” done using address_as_void_pointer. But I don’t quite understand how that function works.
-
All the additional data that has to be passed to the cfunc (pointer to the arrays, shape,… ) using
void *user_data. You could avoid all this using global variables, but than you have to recompile everything once you want to change the inputs. This takes longer than the whole integration without using any low-level callable, so this isn’t a real option here. -
This data (c_struct) can be created using a numpy array with the corresponding dtype
types.Record.make_c_struct([...]), but pointers are not allowed by Numba. Therefore pointers have to be passed as a memory address, but nb.carray needs a pointer. That’s whyaddress_as_void_pointeris needed to perform a reinterpret cast from an int64 value (the memory address) to a void pointer. In C this can be done usingvoid *pointer = &xwhere x is a memory address. If you want to generate a function pointer instead, this function has to be rewritten. I don’t know how this can be accomplished, but the answer should more or less be already somewhere in the Numba or LLVM-lite source code. -
create_jit_integrand_functionis just a wrapper for a Numba function. It would have been also possible to do all the calculation inside this function. -
The last question is now: How to create a low-level callable that has all arguments using the c-struct already available? This can be accomplished by ```
LowLevelCallable(func.ctypes,user_data=args.ctypes.data_as(ctypes.c_void_p))