New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
support for scipy.special functions : feature request #3086
Comments
Agree that this needs to be done. As an interim, it's trivial to create these as needed following the docs http://numba.pydata.org/numba-doc/latest/extending/high-level.html. Here is an example for from numba.extending import get_cython_function_address
from numba import vectorize, njit
import ctypes
import numpy as np
addr = get_cython_function_address("scipy.special.cython_special", "beta")
functype = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double, ctypes.c_double)
beta_fn = functype(addr)
@vectorize('float64(float64, float64)')
def vec_beta(x, y):
return beta_fn(x, y)
@njit
def beta_in_njit(x, y):
return vec_beta(x, y)
d = np.linspace(0, 20, 50)
got = vec_beta(d, d[::-1])
got_njit_context = beta_in_njit(d, d[::-1])
from scipy.special import beta
expected = beta(d, d[::-1])
np.testing.assert_allclose(got, expected)
np.testing.assert_allclose(got_njit_context, expected) |
I also started a pure Numba port of special functions here: https://github.com/person142/special There is a (small) chance that what you want is there. |
@stuartarchibald Thanks a lot for the example. :) |
@stuartarchibald, Thank you very much for your example! Your example works very well for functions involving "double". ctypedef fused Dd_number_t: The numba manual said "check the extension module’s pyx_capi attribute". |
@guoyang0123 the
The fused types in e.g. from numba.extending import get_cython_function_address
from numba import vectorize, njit
import ctypes
import numpy as np
_PTR = ctypes.POINTER
_dble = ctypes.c_double
_ptr_dble = _PTR(_dble)
# signature is:
# void airy(Dd_number_t x0, Dd_number_t *y0, Dd_number_t *y1, Dd_number_t *y2, Dd_number_t *y3) nogil
# bind to the real space variant of the function
addr = get_cython_function_address("scipy.special.cython_special", "__pyx_fuse_1airy")
functype = ctypes.CFUNCTYPE(None, _dble, _ptr_dble, _ptr_dble, _ptr_dble, _ptr_dble)
airy_float64_fn = functype(addr)
@njit
def numba_airy_float64(x):
y0 = np.empty(1, dtype=np.float64)
y1 = np.empty(1, dtype=np.float64)
y2 = np.empty(1, dtype=np.float64)
y3 = np.empty(1, dtype=np.float64)
airy_float64_fn(x, y0.ctypes, y1.ctypes, y2.ctypes, y3.ctypes)
return y0[0], y1[0], y2[0], y3[0]
x = -15
got = numba_airy_float64(x)
from scipy.special import airy
expected = airy(x)
np.testing.assert_allclose(got, expected) |
Dear @stuartarchibald, |
@guoyang0123 no problem. Glad it's working for you! |
@stuartarchibald Thanks a lot! |
@max9111 yes, but be advised that because from numba.extending import get_cython_function_address
from numba import vectorize, njit
import ctypes
import numpy as np
_PTR = ctypes.POINTER
_dble = ctypes.c_double
_ptr_dble = _PTR(_dble)
addr = get_cython_function_address("scipy.special.cython_special", "__pyx_fuse_0airy")
# pass the complex value on the stack as two adjacent double values
functype = ctypes.CFUNCTYPE(None, _dble, _dble, _ptr_dble, _ptr_dble, _ptr_dble, _ptr_dble)
airy_complex128_fn = functype(addr)
@njit
def numba_airy_complex128(x):
y0 = np.empty(2, dtype=np.float64)
y1 = np.empty(2, dtype=np.float64)
y2 = np.empty(2, dtype=np.float64)
y3 = np.empty(2, dtype=np.float64)
airy_complex128_fn(np.real(x), np.imag(x), y0.ctypes, y1.ctypes, y2.ctypes, y3.ctypes)
def ascmplx(z):
return np.complex(z[0] + 1j * z[1])
return ascmplx(y0), ascmplx(y1), ascmplx(y2), ascmplx(y3)
x = -2 + 3j
got = numba_airy_complex128(x)
from scipy.special import airy
expected = airy(x)
np.testing.assert_allclose(got, expected) |
Would this work for functions that return a single complex value, e.g., I can get the real part without issue by using a |
Hi @mpetroff , |
@mpetroff I don't think this is possible right now. The returned complex number is too wide for any of the supported types and it is necessary for it to be declared as the returned item to match complex return behaviour in C99. In practice this is done however the compiler sees fit, but typically through putting the real/imag parts into
I'd expect this to be the case but purely down to luck of code generation. It so happens that the ctypes call places the real part in |
@mpetroff It is always a possibility to write a cython cdef function which exposes a simpler interface. Writing a Cyhon wrapperSpecial.pyx cimport scipy.special.cython_special
Setup.py
Wrapper
@stuartarchibald By commenting out the line |
@max9111 Thanks! I had started writing something with CFFI and |
I raised the issue of returning complex values from C-ABI functions at yesterday's core developer meeting. The outcome was that it's probably possible for Numba to be able to generate a C-ABI compliant wrapper (as clang and GCC are probably doing something sufficiently similar) to handle the case where a complex value is returned from a |
I stumbled upon this issue today which is kind of roadblock for introducing numba to my project. Is there any update or even an ETA for this feature? If not, I see some workarounds. I didn't fully understood them yet, can somebody explain it a little more in detail? I need the |
The key to find the right function is e.g. >>> {key: value for key, value in scipy.special.cython_special.__pyx_capi__.items() if 'erf' in key}
{'__pyx_fuse_0erf': <capsule object "__pyx_t_double_complex (__pyx_t_double_complex, int __pyx_skip_dispatch)" at 0x7fc717fd5780>,
'__pyx_fuse_1erf': <capsule object "double (double, int __pyx_skip_dispatch)" at 0x7fc717fd57b0>,
'__pyx_fuse_0erfc': <capsule object "__pyx_t_double_complex (__pyx_t_double_complex, int __pyx_skip_dispatch)" at 0x7fc717fd57e0>,
'__pyx_fuse_1erfc': <capsule object "double (double, int __pyx_skip_dispatch)" at 0x7fc717fd5810>,
'__pyx_fuse_0erfcx': <capsule object "__pyx_t_double_complex (__pyx_t_double_complex, int __pyx_skip_dispatch)" at 0x7fc717fd5840>,
'__pyx_fuse_1erfcx': <capsule object "double (double, int __pyx_skip_dispatch)" at 0x7fc717fd5870>,
'__pyx_fuse_0erfi': <capsule object "__pyx_t_double_complex (__pyx_t_double_complex, int __pyx_skip_dispatch)" at 0x7fc717fd58a0>,
'__pyx_fuse_1erfi': <capsule object "double (double, int __pyx_skip_dispatch)" at 0x7fc717fd58d0>} Then you can choose the one corresponds to the signature you need. |
Ok, so this gives me the info to construct the native call which numba cannot generate at the moment, if I understand correctly? |
Yes, following the airy example above together with the function you find out then you can construct a jittable function providing the same result. |
Ok, thanks for the clarification. Just as a comment: This functions are crucial for real science. Numba really needs to support this out of the box. |
Hm I don’t know that the burden should be entirely placed on Numba here because it’s unreasonable to expect it to reimplement the API of every major scientific Python package. At a certain point the community needs to figure out a better interop story. |
More concretely: I’m very interested in ideas as to how we could make |
Something like |
Making Numba able to compile SciPy math functions has been on our radar for a while (see #2109, @stuartarchibald even came up with the project name of "Scumba" 😄 ), but we've been trying to solidify our core functionality and NumPy support first. We would prefer to do it as a Numba extension (that we would be happy to host under the numba org). On the topic of minimizing code, we have wanted for a long time to figure out how to call Cython functions from Numba as easily as we can call CFFI and ctypes wrapped functions. If that would allow us to more directly use interfaces from SciPy and not have to reimplement stuff, that would be great. |
(Arrived here from SciPy-Dev) FWIW, this is an experiment I did to wrap directly the underlying C code with CFFI and call it from numba: https://github.com/poliastro/pycephes/tree/v0.1.0 (Explanatory article translated to English) @horta then followed a similar approach in https://github.com/limix/hcephes. |
Following the advises of @stuartarchibald and @ickc I solved my issues. A short question: how much likely is it that this solution will break with future updates, i.e. if scipy people make special more canonical? I really got nice speedups in my case. Without parallel it's on par with cython, with double as fast in my use-case. So I can throw out the cython dependency and have much easier maintenance. Thanks to Numba team! |
Not sure what the durability of this solution will be, but it has worked for some time now, and Cython extensions depend on this interface, so it is not likely to go away. We're having a discussion on the SciPy mailing list now to see if there is a more official way we could migrate toward. |
Followers of this ticket may also be interested in scumba discussion here: #4292 |
The Numba extension project for SciPy: https://github.com/numba/numba-scipy is now the place to request for Numba support of SciPy functionality. |
Closing this ticket, discussion can be moved to issues on: https://github.com/numba/numba-scipy |
If possible, it would be great if numba supports scipy special functions which are mostly universal functions and there is cython API for those functions as well.
The text was updated successfully, but these errors were encountered: