Skip to content
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

Add numpy windowing functions support (np.bartlett, np.hamming, np.blackman, np.hanning, np.kaiser) #4076

Merged
merged 11 commits into from
May 24, 2019
5 changes: 5 additions & 0 deletions docs/source/reference/numpysupported.rst
Original file line number Diff line number Diff line change
Expand Up @@ -307,7 +307,9 @@ The following top-level functions are supported:
* :func:`numpy.atleast_1d`
* :func:`numpy.atleast_2d`
* :func:`numpy.atleast_3d`
* :func:`numpy.bartlett`
* :func:`numpy.bincount` (only the 2 first arguments)
* :func:`numpy.blackman`
* :func:`numpy.column_stack`
* :func:`numpy.concatenate`
* :func:`numpy.convolve` (only the 2 first arguments)
Expand All @@ -332,9 +334,12 @@ The following top-level functions are supported:
* :func:`numpy.frombuffer` (only the 2 first arguments)
* :func:`numpy.full` (only the 3 first arguments)
* :func:`numpy.full_like` (only the 3 first arguments)
* :func:`numpy.hamming`
* :func:`numpy.hanning`
* :func:`numpy.histogram` (only the 3 first arguments)
* :func:`numpy.hstack`
* :func:`numpy.identity`
* :func:`numpy.kaiser`
* :func:`numpy.interp` (only the 3 first arguments; requires NumPy >= 1.10)
* :func:`numpy.linspace` (only the 3-argument form)
* :class:`numpy.ndenumerate`
Expand Down
173 changes: 172 additions & 1 deletion numba/targets/arraymath.py
Original file line number Diff line number Diff line change
Expand Up @@ -2857,7 +2857,7 @@ def np_delete(arr, obj):

if isinstance(obj, (types.Array, types.Sequence, types.SliceType)):
if isinstance(obj, (types.SliceType)):
handler = np_delete_handler_isslice
handler = np_delete_handler_isslice
else:
if not isinstance(obj.dtype, types.Integer):
raise TypingError('obj should be of Integer dtype')
Expand Down Expand Up @@ -3540,3 +3540,174 @@ def np_extract_impl(condition, arr):
return np.array(out)

return np_extract_impl

stuartarchibald marked this conversation as resolved.
Show resolved Hide resolved
#----------------------------------------------------------------------------
# Windowing functions
# - translated from the numpy implementations found in:
# https://github.com/numpy/numpy/blob/v1.16.1/numpy/lib/function_base.py#L2543-L3233
# at commit: f1c4c758e1c24881560dd8ab1e64ae750

@register_jitable
def np_bartlett_impl(M):
n = np.arange(M)
return np.where(np.less_equal(n, (M - 1) / 2.0), 2.0 * n / (M - 1),
2.0 - 2.0 * n / (M - 1))


@register_jitable
def np_blackman_impl(M):
n = np.arange(M)
return (0.42 - 0.5 * np.cos(2.0 * np.pi * n / (M - 1)) +
0.08 * np.cos(4.0* np.pi * n / (M - 1)))


@register_jitable
def np_hamming_impl(M):
n = np.arange(M)
return 0.54 - 0.46 * np.cos(2.0 * np.pi * n / (M - 1))


@register_jitable
def np_hanning_impl(M):
n = np.arange(M)
return 0.5 - 0.5 * np.cos(2.0 * np.pi * n / (M - 1))


def window_generator(func):
def window_overload(M):
if not isinstance(M, types.Integer):
raise TypingError('M must be an integer')

def window_impl(M):

if M < 1:
return np.array((), dtype=np.float_)
if M == 1:
return np.ones(1, dtype=np.float_)
return func(M)

return window_impl
return window_overload

overload(np.bartlett)(window_generator(np_bartlett_impl))
overload(np.blackman)(window_generator(np_blackman_impl))
overload(np.hamming)(window_generator(np_hamming_impl))
overload(np.hanning)(window_generator(np_hanning_impl))


_i0A = np.array([
-4.41534164647933937950E-18,
3.33079451882223809783E-17,
-2.43127984654795469359E-16,
1.71539128555513303061E-15,
-1.16853328779934516808E-14,
7.67618549860493561688E-14,
-4.85644678311192946090E-13,
2.95505266312963983461E-12,
-1.72682629144155570723E-11,
9.67580903537323691224E-11,
-5.18979560163526290666E-10,
2.65982372468238665035E-9,
-1.30002500998624804212E-8,
6.04699502254191894932E-8,
-2.67079385394061173391E-7,
1.11738753912010371815E-6,
-4.41673835845875056359E-6,
1.64484480707288970893E-5,
-5.75419501008210370398E-5,
1.88502885095841655729E-4,
-5.76375574538582365885E-4,
1.63947561694133579842E-3,
-4.32430999505057594430E-3,
1.05464603945949983183E-2,
-2.37374148058994688156E-2,
4.93052842396707084878E-2,
-9.49010970480476444210E-2,
1.71620901522208775349E-1,
-3.04682672343198398683E-1,
6.76795274409476084995E-1
])

_i0B = np.array([
-7.23318048787475395456E-18,
-4.83050448594418207126E-18,
4.46562142029675999901E-17,
3.46122286769746109310E-17,
-2.82762398051658348494E-16,
-3.42548561967721913462E-16,
1.77256013305652638360E-15,
3.81168066935262242075E-15,
-9.55484669882830764870E-15,
-4.15056934728722208663E-14,
1.54008621752140982691E-14,
3.85277838274214270114E-13,
7.18012445138366623367E-13,
-1.79417853150680611778E-12,
-1.32158118404477131188E-11,
-3.14991652796324136454E-11,
1.18891471078464383424E-11,
4.94060238822496958910E-10,
3.39623202570838634515E-9,
2.26666899049817806459E-8,
2.04891858946906374183E-7,
2.89137052083475648297E-6,
6.88975834691682398426E-5,
3.36911647825569408990E-3,
8.04490411014108831608E-1
])


@register_jitable
def _chbevl(x, vals):
b0 = vals[0]
b1 = 0.0

for i in range(1, len(vals)):
b2 = b1
b1 = b0
b0 = x * b1 - b2 + vals[i]

return 0.5 * (b0 - b2)


@register_jitable
stuartarchibald marked this conversation as resolved.
Show resolved Hide resolved
def _i0(x):
if x < 0:
x = -x
if x <= 8.0:
y = (0.5 * x) - 2.0
return np.exp(x) * _chbevl(y, _i0A)

return np.exp(x) * _chbevl(32.0 / x - 2.0, _i0B) / np.sqrt(x)


@register_jitable
def _i0n(n, alpha, beta):
y = np.empty_like(n, dtype=np.float_)
t = _i0(np.float_(beta))
for i in range(len(y)):
y[i] = _i0(beta * np.sqrt(1 - ((n[i] - alpha) / alpha)**2.0)) / t

return y


@overload(np.kaiser)
def np_kaiser(M, beta):
stuartarchibald marked this conversation as resolved.
Show resolved Hide resolved
if not isinstance(M, types.Integer):
raise TypingError('M must be an integer')

if not isinstance(beta, (types.Integer, types.Float)):
raise TypingError('beta must be an integer or float')

def np_kaiser_impl(M, beta):
if M < 1:
return np.array((), dtype=np.float_)
if M == 1:
return np.ones(1, dtype=np.float_)

n = np.arange(0, M)
alpha = (M - 1) / 2.0

return _i0n(n, alpha, beta)

return np_kaiser_impl
64 changes: 64 additions & 0 deletions numba/tests/test_np_functions.py
Original file line number Diff line number Diff line change
Expand Up @@ -193,6 +193,26 @@ def array_repeat(a, repeats):
return np.asarray(a).repeat(repeats)


def np_bartlett(M):
return np.bartlett(M)


def np_blackman(M):
return np.blackman(M)


def np_hamming(M):
return np.hamming(M)


def np_hanning(M):
return np.hanning(M)


def np_kaiser(M, beta):
return np.kaiser(M, beta)


class TestNPFunctions(MemoryLeakMixin, TestCase):
"""
Tests for various Numpy functions.
Expand Down Expand Up @@ -2748,6 +2768,50 @@ def test_repeat_exception(self):
with self.assertRaises(TypingError):
nbfunc(np.ones(1), rep)

def test_windowing(self):
def check_window(func):
np_pyfunc = func
np_nbfunc = njit(func)

for M in [0, 1, 5, 12]:
expected = np_pyfunc(M)
got = np_nbfunc(M)
self.assertPreciseEqual(expected, got)

for M in ['a', 1.1, 1j]:
with self.assertRaises(TypingError) as raises:
np_nbfunc(1.1)
self.assertIn("M must be an integer", str(raises.exception))

check_window(np_bartlett)
check_window(np_blackman)
check_window(np_hamming)
check_window(np_hanning)

# Test np.kaiser separately
np_pyfunc = np_kaiser
np_nbfunc = njit(np_kaiser)

for M in [0, 1, 5, 12]:
for beta in [0.0, 5.0, 14.0]:
expected = np_pyfunc(M, beta)
got = np_nbfunc(M, beta)

if IS_32BITS:
self.assertPreciseEqual(expected, got, prec='double', ulps=2)
else:
self.assertPreciseEqual(expected, got, prec='exact')

for M in ['a', 1.1, 1j]:
with self.assertRaises(TypingError) as raises:
np_nbfunc(M, 1.0)
self.assertIn("M must be an integer", str(raises.exception))

for beta in ['a', 1j]:
with self.assertRaises(TypingError) as raises:
np_nbfunc(5, beta)
self.assertIn("beta must be an integer or float", str(raises.exception))


class TestNPMachineParameters(TestCase):
# tests np.finfo, np.iinfo, np.MachAr
Expand Down