Problems using numba and math functions from scipy.special.cython_special

Numba Community Support
1

Maybe someone could help me:

Python 3.7
scipy 1.5.0
numpy 1.18.1
numba 0.49.1
numba-scipy 0.1.0
numba-special 0.2.0

> import numpy as np
> import numba as nb
> import scipy.special.cython_special

> print(scipy.special.cython_special.eval_legendre)
> print(scipy.special.cython_special.eval_legendre(2.0, 0.546))

I got

>> <cyfunction eval_legendre at 0x0000000012CAF588>
>> -0.05282599999999993

> addr = get_cython_function_address("scipy.special.cython_special", "eval_legendre")
> functype = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double, ctypes.c_double)
> eval_legendre_fn = functype(addr)

I got

File “C:\Users\Enot\PycharmProjects\GPU-DL\venv\lib\site-packages\numba\core\extending.py”, line 407, in get_cython_function_address return _import_cython_function(module_name, function_name)
ValueError: No function ‘eval_legendre’ found in pyx_capi of ‘scipy.special.cython_special’

Good people from https://gitter.im/numba/numba prompt me that:

`eval_legendre is not part of scipy.special.cython_special.pyx_capi (at least not with that name) and therefore numba simply cannot find it. And I have to find some kind of interface to a lower level function with another name.

'__pyx_fuse_0_0eval_legendre': <capsule object "__pyx_t_double_complex (double, __pyx_t_double_complex, int __pyx_skip_dispatch)" at 0x7f575f1c84e0>,
    '__pyx_fuse_0_1eval_legendre': <capsule object "double (double, double, int __pyx_skip_dispatch)" at 0x7f575f1c8510>,
    '__pyx_fuse_1_0eval_legendre': <capsule object "__pyx_t_double_complex (long, __pyx_t_double_complex, int __pyx_skip_dispatch)" at 0x7f575f1c8540>,
    '__pyx_fuse_1_1eval_legendre': <capsule object "double (long, double, int __pyx_skip_dispatch)" at 0x7f575f1c8570>,

I do not need complex numbers so I choose

‘__pyx_fuse_0_1eval_legendre’: <capsule object “double (double, double, int __pyx_skip_dispatch)” at 0x7f575f1c8510>

According to scipy docs: scipy.special.eval_legendre(n, x, out=None) n: array_like degree of the polynomial. xarray_like - points at which to evaluate
It is logical to assume that one of the "double"s - output; other double is input x; and degree of the polynomial is int.
So using an example from the link: support for scipy.special functions : feature request · Issue #3086 · numba/numba · GitHub and link complex argument support for erfc · Issue #1 · person142/numba_special · GitHub

I made smth like:

> import numba as nb
> from numba.extending import intrinsic
> from numba.core import cgutils
> import ctypes
> from numba import types
> from numba.extending import get_cython_function_address
> import scipy.special.cython_special
> from numba import njit
> 
> _PTR = ctypes.POINTER
> _dble = ctypes.c_double
> _ptr_dble = _PTR(_dble)
> 
> @intrinsic
> def val_to_double_ptr(typingctx, data):
>     def impl(context, builder, signature, args):
>         ptr = cgutils.alloca_once_value(builder,args[0])
>         return ptr
>     sig = types.CPointer(types.float64)(types.float64)
>     return sig, impl
> 
> @intrinsic
> def double_ptr_to_val(typingctx, data):
>     def impl(context, builder, signature, args):
>         val = builder.load(args[0])
>         return val
>     sig = types.float64(types.CPointer(types.float64))
>     return sig, impl
> 
> addr = get_cython_function_address("scipy.special.cython_special", "__pyx_fuse_0_1eval_legendre")
> functype = ctypes.CFUNCTYPE(None, _ptr_dble, _dble, _int)
> eval_legendre_float64_fn = functype(addr)
> 
> @njit
> def numba_eval_legendre_float64(n, x):
>     out_p = val_to_double_ptr(0.)
>     eval_legendre_float64_fn(out_p, x, n)
>     out = double_ptr_to_val(out_p)
>     return out
> 
> print("eval legendre: ", numba_eval_legendre_float64(2, 1.567))

This raises no error but return

> >> eval legendre: 0.0

So do nothing.

Any ideas on how to fix this?

Thank You.

2

the following works, but it does not seem to return the correct answer

addr = get_cython_function_address("scipy.special.cython_special", "__pyx_fuse_0_1eval_legendre")
functype = ctypes.CFUNCTYPE(_dble, _dble, _int)
eval_legendre_float64_fn = functype(addr)

@njit
def numba_eval_legendre_float64(n, x):
    out = eval_legendre_float64_fn(x, n)
    return out

I think you should try using numba-scipy, instead of linking to the cython function directly. Numba-scipy should allow you to avoid all the address and functype code, and just write normal python. Just install it with pip or conda and the following should work:

@njit
def numba_eval_legendre_float64(n, x):
    out = eval_legendre(x, n)
    return out
3

Thank You.

This works:

> addr = get_cython_function_address("scipy.special.cython_special", "__pyx_fuse_0_1eval_legendre")
> functype = ctypes.CFUNCTYPE(_dble, _dble, _dble)
> eval_legendre_float64_fn = functype(addr)
> 
> @njit
> def numba_eval_legendre_float64(n, x):
>     return eval_legendre_float64_fn(n, x)
> 
> g = numba_eval_legendre_float64(2.0, 1.567)
> print("eval legendre: ", g)
> print(scipy.special.cython_special.eval_legendre(2.0, 1.567))
> print(scipy.special.eval_legendre(2.0, 1.567))
> print(scipy.special.legendre(2.0)(1.567))

Same result: 3.1832334999999996

4

I install all this. Trying to run example code from numba-special · PyPI
> import numba
> import scipy.special as sc
> import numba_special # The import generates Numba overloads for special
> @numba.njit
> def gamma_plus_1(x): return sc.gamma(x) + 1.0
> gamma_plus_1(5.0)

I got:

> Traceback (most recent call last):
>   File "C:/Users/Enot/PycharmProjects/GPU-DL/cython_w/test_question_cython.py", line 119, in <module>
>     import numba_special  # The import generates Numba overloads for special
>   File "C:\Users\Enot\PycharmProjects\GPU-DL\venv\lib\site-packages\numba_special\__init__.py", line 745, in <module>
>     from . import numba_overloads
>   File "C:\Users\Enot\PycharmProjects\GPU-DL\venv\lib\site-packages\numba_special\numba_overloads.py", line 7, in <module>
>     from . import function_pointers
>   File "numba_special/function_pointers.pyx", line 2, in init numba_special.function_pointers
>     #cython: language_level=3
> TypeError: C function scipy.special.cython_special.__pyx_fuse_1bdtr has wrong signature (expected double (long, long, double, int __pyx_skip_dispatch), got double (double, long, double, int __pyx_skip_dispatch))

so I have to go the other way

6

Hi Racoon,
can you try running this without importing numba-special - I imagine that it is clashing with numba-scipy since they are trying to achieve the same thing. Afaik you do not have to import numba-scipy, it just needs to be installed in order to work because it registers certain extensions in the background.

7

I have to press Flag button at the bottom of the page, under all replies, between button Bookmark and Reply?

8

I try do run example without import numba_special and get Error:

> Traceback (most recent call last):
>   File "C:/Users/Enot/PycharmProjects/GPU-DL/cython_w/test_question_cython.py", line 127, in <module>
>     gamma_plus_1(5.0)
>   File "C:\Users\Enot\PycharmProjects\GPU-DL\venv\lib\site-packages\numba\core\dispatcher.py", line 401, in _compile_for_args
>     error_rewrite(e, 'typing')
>   File "C:\Users\Enot\PycharmProjects\GPU-DL\venv\lib\site-packages\numba\core\dispatcher.py", line 344, in error_rewrite
>     reraise(type(e), e, None)
>   File "C:\Users\Enot\PycharmProjects\GPU-DL\venv\lib\site-packages\numba\core\utils.py", line 80, in reraise
>     raise value.with_traceback(tb)
> numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
> Unknown attribute 'gamma' of type Module(<module 'scipy.special' from 'C:\\Users\\Enot\\PycharmProjects\\GPU-DL\\venv\\lib\\site-packages\\scipy\\special\\__init__.py'>)
> 
> File "test_question_cython.py", line 125:
> def gamma_plus_1(x):
>     return sc.gamma(x) + 1.0
>     ^
> 
> [1] During: typing of get attribute at C:/Users/Enot/PycharmProjects/GPU-DL/cython_w/test_question_cython.py (125)
> 
> File "test_question_cython.py", line 125:
> def gamma_plus_1(x):
>     return sc.gamma(x) + 1.0
9

Nevermind this, I just answered to the wrong thread, that’s why I deleted this answer :slight_smile:

10

mhmmm, I am currently struggling with numba-scipy due to https://github.com/numba/numba-scipy/issues/35
so it’s a little tricky to help out :smiley:
But looking at the git repo, numba-scipy at version 0.1.0 only supported a single function: beta, so I guess that may be related to your problem.

You would probably have to install numba-scipy from the branch in this PR https://github.com/numba/numba-scipy/pull/37 in order to get it running with scipy 1.5.0 (that PR is working on a Patch for the bug I encountered)
Then I expect it should just work - Fingers crossed :smiley:

11

Thank You. I’ll try.

12

Maybe You could help me with the next problem.

Assume that scipy.special.eval_legendre or scipy.special.jv works well with numba or I use trick listed above with “get_cython_function_address”.

I have to perform fast numerical integration of 2D matrix.

> cos_beta1 = 0.99999 
> cos_beta2 = 0.00001
> k = 2* np.pi / 600
> 
> num = random.random()
> ro_p = np.full((101, 101), num)
> 
> z = random.randint(-5, 5)
> z_p = np.full((101, 101), z)
> 
> #### f1, f2, f3 - math expressions that contain legendre polynomial and bessel functions
> 
> @cfunc("float64(float64, float64, float64[:], float64[:])")
> def pi_plus_tau(n, k, ro_p, z_p):
> 
>     @nb.njit()
>     def frst_deriv(x):
>         return math_expression f1(n, x)
>         
>     @nb.njit()
>     def sec_deriv(x):
> 	       return math_expression f2(n, x)
>         
>     @nb.njit()
>     def integrand(x):
>         return f3(k, ro_p, z_p) * frst_deriv(x) * sec_deriv(x)  
> 
>     return integrand
> 
> def intergration_plus(n, angle1, angle2):
>     return (integrate.quad_vec(pi_plus_tau(float(n), k, ro_p, z_p), angle1, angle2, workers=1))[0] 
> 
> for n in range(1, 10):
>     intergration_plus(n, cos_beta1, cos_beta2)

This works well without numba but quite slow.

With @cfunc and @nb.njit this raise an exception:

> numba.core.errors.TypingError: Failed in nopython mode pipeline (step: convert make_function into JIT functions)
> Cannot capture the non-constant value associated with variable 'n' in a function that will escape.
> 
> File "test_question_cython.py", line 215:
> def pi_plus_tau(n, k, ro_p, z_p):
>     <source elided>
> 
>     @nb.njit()

Thank You.

13

have you tried removing the @cfunc from pi_plus_tau? Function factories don’t need to be jitted. What matters is that integrand has @njit but pi_plus_tau can (and should) be pure python.

14

Fine! Without @cfunc above pi_plus_tau code works good. It remains to remake the bessel function numba_jv to work with matrices.
Thanks you!

15

Most of the time program is spent on executing the scipy.integrate.quad_vec function and I cannot speed up the process with the help of numba module. But I want to try to use numba with scipy.integrate.quad function
Do smth. like:

> a = np.empty(shape=(image_size+1, image_size+1), dtype=np.complex)
> for i in range(image_size+1):
>     for j in range(image_size+1):
> 	       a[i, j] = scipy.integrate.quad(some_func, 0, 1)

So to calculate integrals separately.

A simple example to play with:

> import numpy as np
> import scipy
> 
> fun2int = lambda x, a: np.sqrt(x+a)
> intfun = lambda a: scipy.integrate.quad(fun2int, 0, 4, args=(a))[0]
> vec_int = np.vectorize(intfun)
> print("VECTORIZE: ", vec_int(np.arange(0, 1000000).reshape(1000, 1000)))

How I can speed up this small example with numba?

Thank You.

16

Have you read this https://numba.pydata.org/numba-doc/latest/user/cfunc.html?

It has an example of how to speed up quad.

17

Right, sorry. I have already done it.

18

Sorry to disturb you, I want to know how could you find a lower-level function in pyx_capi, I can’t find this file and need to know if there is low-level Bessel functions.

19

hi @rlibd , does scipy have the function you want? Special functions (scipy.special) — SciPy v1.6.0 Reference Guide

if it does, you can check numba-scipy supported list: numba-scipy/signatures.py at master · numba/numba-scipy · GitHub

20

Thanks a lot! I find what I need in the Github list.