La documentación discute el uso de los cfunc de cfunc
como argumento scipy.integrate.quad
de scipy.integrate.quad
. Necesito lo mismo con parámetro adicional.
Básicamente estoy tratando de hacer algo como esto:
import numpy as np from numba import cfunc import numba.types voidp = numba.types.voidptr def integrand(t, params): a = params[0] # this is additional parameter return np.exp(-t/a) / t**2 nb_integrand = cfunc(numba.float32(numba.float32, voidp))(integrand)
Sin embargo, no funciona, porque se supone que los voidptr
son voidptr
/ void*
y no pueden transformarse al double
. Tengo el siguiente mensaje de error:
TypingError: Failed at nopython (nopython frontend) Invalid usage of getitem with parameters (void*, int64) * parameterized
No encontré ninguna información sobre cómo extraer valores de void*
en Numba. En C, debería ser algo como a = *((double*) params)
. ¿Es posible hacer lo mismo en Numba?
1. Pasando argumentos adicionales a través de scipy.integrate.quad
Los documentos quad
dicen:
Si el usuario desea mejorar el rendimiento de la integración, entonces
f
puede ser unscipy.LowLevelCallable
con una de las firmas:
double func(double x)
double func(double x, void *user_data)
double func(int n, double *xx)
double func(int n, double *xx, void *user_data)
El
user_data
es la información contenida en elscipy.LowLevelCallable
. En las formas de llamada conxx
,n
es la longitud de la matrizxx
que contienexx[0] == x
y el rest de los elementos son números contenidos en el argumentoargs
dequad
.
Por lo tanto, para pasar un argumento adicional a integrand
través de quad
, es mejor usar la firma double func(int n, double *xx)
.
Puede escribir un decorador a su función integrand para transformarlo en un LowLevelCallable
como:
import numpy as np import scipy.integrate as si import numba from numba import cfunc from numba.types import intc, CPointer, float64 from scipy import LowLevelCallable def jit_integrand_function(integrand_function): jitted_function = numba.jit(integrand_function, nopython=True) @cfunc(float64(intc, CPointer(float64))) def wrapped(n, xx): return jitted_function(xx[0], xx[1]) return LowLevelCallable(wrapped.ctypes) @jit_integrand_function def integrand(t, *args): a = args[0] return np.exp(-t/a) / t**2 def do_integrate(func, a): """ Integrate the given function from 1.0 to +inf with additional argument a. """ return si.quad(func, 1, np.inf, args=(a,)) print(do_integrate(integrand, 2.)) >>>(0.326643862324553, 1.936891932288535e-10)
O si no desea el decorador, cree el LowLevelCallable
manualmente y páselo a quad
.
2. Envolviendo la función integrand
No estoy seguro de que lo siguiente cumpla con sus requisitos, pero también podría ajustar su función integrand
para lograr el mismo resultado:
import numpy as np from numba import cfunc import numba.types def get_integrand(*args): a = args[0] def integrand(t): return np.exp(-t/a) / t**2 return integrand nb_integrand = cfunc(numba.float64(numba.float64))(get_integrand(2.)) import scipy.integrate as si def do_integrate(func): """ Integrate the given function from 1.0 to +inf. """ return si.quad(func, 1, np.inf) print(do_integrate(get_integrand(2))) >>>(0.326643862324553, 1.936891932288535e-10) print(do_integrate(nb_integrand.ctypes)) >>>(0.326643862324553, 1.936891932288535e-10)
3. Casting desde voidptr
a un tipo python
No creo que esto sea posible todavía. De esta discusión en 2016, parece que voidptr
solo está aquí para pasar un contexto a una callback de C.
El caso del puntero void * sería para las API en las que el código C extranjero no intenta desreferir el puntero, sino que simplemente lo devuelve a la callback como una forma para que la callback retenga el estado entre las llamadas. No creo que sea particularmente importante en este momento, pero quería plantear el problema.
Y probando lo siguiente:
numba.types.RawPointer('p').can_convert_to( numba.typing.context.Context(), CPointer(numba.types.Any))) >>>None
¡Tampoco parece alentador!