From ffd2d8095c8034cdcaec5bf52c1083389eb943a7 Mon Sep 17 00:00:00 2001 From: Vahid Tavanashad Date: Mon, 24 Jun 2024 02:17:06 -0500 Subject: [PATCH 1/4] implement fftfreq --- dpnp/fft/dpnp_iface_fft.py | 202 +++++++++++++++++-- tests/test_fft.py | 19 ++ tests/test_sycl_queue.py | 14 ++ tests/test_usm_type.py | 14 ++ tests/third_party/cupy/fft_tests/test_fft.py | 3 +- 5 files changed, 238 insertions(+), 14 deletions(-) diff --git a/dpnp/fft/dpnp_iface_fft.py b/dpnp/fft/dpnp_iface_fft.py index c8064e0122ea..6bf8ab292ad7 100644 --- a/dpnp/fft/dpnp_iface_fft.py +++ b/dpnp/fft/dpnp_iface_fft.py @@ -184,19 +184,113 @@ def fft2(x, s=None, axes=(-2, -1), norm=None): return call_origin(numpy.fft.fft2, x, s, axes, norm) -def fftfreq(n=None, d=1.0): +def fftfreq(n, d=1.0, device=None, usm_type=None, sycl_queue=None): """ - Compute the one-dimensional discrete Fourier Transform sample frequencies. + Return the Discrete Fourier Transform sample frequencies. + + The returned float array `f` contains the frequency bin centers in cycles + per unit of the sample spacing (with zero at the start). For instance, if + the sample spacing is in seconds, then the frequency unit is cycles/second. + + Given a window length `n` and a sample spacing `d`:: + + f = [0, 1, ..., n/2-1, -n/2, ..., -1] / (d*n) if n is even + f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n) if n is odd For full documentation refer to :obj:`numpy.fft.fftfreq`. - Limitations - ----------- - Parameter `d` is unsupported. + Parameters + ---------- + n : int + Window length. + d : scalar, optional + Sample spacing (inverse of the sampling rate). + Default: ``1``. + device : {None, string, SyclDevice, SyclQueue}, optional + An array API concept of device where the output array is created. + The `device` can be ``None`` (the default), an OneAPI filter selector + string, an instance of :class:`dpctl.SyclDevice` corresponding to + a non-partitioned SYCL device, an instance of :class:`dpctl.SyclQueue`, + or a `Device` object returned by + :obj:`dpnp.dpnp_array.dpnp_array.device` property. + usm_type : {None, "device", "shared", "host"}, optional + The type of SYCL USM allocation for the output array. + Default: ``None``. + sycl_queue : {None, SyclQueue}, optional + A SYCL queue to use for output array allocation and copying. + Default: ``None``. + + Returns + ------- + f : dpnp.ndarray + Array of length `n` containing the sample frequencies. + + See Also + -------- + :obj:`dpnp.fft.rfftfreq` : Return the Discrete Fourier Transform sample + frequencies (for usage with :obj:`dpnp.fft.rfft` and + :obj:`dpnp.fft.irfft`). + + Examples + -------- + >>> import dpnp as np + >>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5]) + >>> fourier = np.fft.fft(signal) + >>> n = signal.size + >>> timestep = 0.1 + >>> freq = np.fft.fftfreq(n, d=timestep) + >>> freq + array([ 0. , 1.25, 2.5 , 3.75, -5. , -3.75, -2.5 , -1.25]) + + Creating the output array on a different device or with a + specified usm_type: + + >>> x = np.fft.fftfreq(n, d=timestep) # default case + >>> x.shape, x.device, x.usm_type + ((8, ), Device(level_zero:gpu:0), 'device') + + >>> y = np.fft.fftfreq(n, d=timestep, device="cpu") + >>> y.shape, y.device, y.usm_type + ((8, ), Device(opencl:cpu:0), 'device') + + >>> z = np.fft.fftfreq(n, d=timestep, usm_type="host") + >>> z.shape, z.device, z.usm_type + ((8, ), Device(level_zero:gpu:0), 'host') """ - return call_origin(numpy.fft.fftfreq, n, d) + if not isinstance(n, int): + raise ValueError("`n` should be an integer") + if not dpnp.isscalar(d): + raise ValueError("`d` should be an scalar") + val = 1.0 / (n * d) + results = dpnp.empty( + n, + dtype=dpnp.intp, + device=device, + usm_type=usm_type, + sycl_queue=sycl_queue, + ) + N = (n - 1) // 2 + 1 + p1 = dpnp.arange( + 0, + N, + dtype=dpnp.intp, + device=device, + usm_type=usm_type, + sycl_queue=sycl_queue, + ) + results[:N] = p1 + p2 = dpnp.arange( + -(n // 2), + 0, + dtype=dpnp.intp, + device=device, + usm_type=usm_type, + sycl_queue=sycl_queue, + ) + results[N:] = p2 + return results * val def fftn(x, s=None, axes=None, norm=None): @@ -849,19 +943,103 @@ def rfft2(x, s=None, axes=(-2, -1), norm=None): return call_origin(numpy.fft.rfft2, x, s, axes, norm) -def rfftfreq(n=None, d=1.0): +def rfftfreq(n, d=1.0, device=None, usm_type=None, sycl_queue=None): """ - Compute the one-dimensional discrete Fourier Transform sample frequencies. + Return the Discrete Fourier Transform sample frequencies + (for usage with :obj:`dpnp.fft.rfft`, :obj:`dpnp.fft.irfft`). + + The returned float array `f` contains the frequency bin centers in cycles + per unit of the sample spacing (with zero at the start). For instance, if + the sample spacing is in seconds, then the frequency unit is cycles/second. + + Given a window length `n` and a sample spacing `d`:: + + f = [0, 1, ..., n/2-1, n/2] / (d*n) if n is even + f = [0, 1, ..., (n-1)/2-1, (n-1)/2] / (d*n) if n is odd + + Unlike :obj:`dpnp.fft.fftfreq` the Nyquist frequency component is + considered to be positive. For full documentation refer to :obj:`numpy.fft.rfftfreq`. - Limitations - ----------- - Parameter `d` is unsupported. + Parameters + ---------- + n : int + Window length. + d : scalar, optional + Sample spacing (inverse of the sampling rate). + Default: ``1``. + device : {None, string, SyclDevice, SyclQueue}, optional + An array API concept of device where the output array is created. + The `device` can be ``None`` (the default), an OneAPI filter selector + string, an instance of :class:`dpctl.SyclDevice` corresponding to + a non-partitioned SYCL device, an instance of :class:`dpctl.SyclQueue`, + or a `Device` object returned by + :obj:`dpnp.dpnp_array.dpnp_array.device` property. + usm_type : {None, "device", "shared", "host"}, optional + The type of SYCL USM allocation for the output array. + Default: ``None``. + sycl_queue : {None, SyclQueue}, optional + A SYCL queue to use for output array allocation and copying. + Default: ``None``. + + Returns + ------- + f : dpnp.ndarray + Array of length `n` containing the sample frequencies. + + See Also + -------- + :obj:`dpnp.fft.fftfreq` : Return the Discrete Fourier Transform sample + frequencies (for usage with :obj:`dpnp.fft.fft` and + :obj:`dpnp.fft.ifft`). + + Examples + -------- + >>> import dpnp as np + >>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5, -3, 4]) + >>> fourier = np.fft.fft(signal) + >>> n = signal.size + >>> sample_rate = 100 + >>> freq = np.fft.fftfreq(n, d=1./sample_rate) + >>> freq + array([ 0., 10., 20., 30., 40., -50., -40., -30., -20., -10.]) + >>> freq = np.fft.rfftfreq(n, d=1./sample_rate) + >>> freq + array([ 0., 10., 20., 30., 40., 50.]) + + Creating the output array on a different device or with a + specified usm_type: + + >>> x = np.fft.rfftfreq(n, d=1./sample_rate) # default case + >>> x.shape, x.device, x.usm_type + ((6, ), Device(level_zero:gpu:0), 'device') + + >>> y =np.fft.rfftfreq(n, d=1./sample_rate, device="cpu") + >>> y.shape, y.device, y.usm_type + ((6, ), Device(opencl:cpu:0), 'device') + + >>> z = np.fft.rfftfreq(n, d=1./sample_rate, usm_type="host") + >>> z.shape, z.device, z.usm_type + ((6, ), Device(level_zero:gpu:0), 'host') """ - return call_origin(numpy.fft.rfftfreq, n, d) + if not isinstance(n, int): + raise ValueError("`n` should be an integer") + if not dpnp.isscalar(d): + raise ValueError("`d` should be an scalar") + val = 1.0 / (n * d) + N = n // 2 + 1 + results = dpnp.arange( + 0, + N, + dtype=dpnp.intp, + device=device, + usm_type=usm_type, + sycl_queue=sycl_queue, + ) + return results * val def rfftn(x, s=None, axes=None, norm=None): diff --git a/tests/test_fft.py b/tests/test_fft.py index b439ef38cce6..79c0e9dc54a0 100644 --- a/tests/test_fft.py +++ b/tests/test_fft.py @@ -1,5 +1,6 @@ import numpy import pytest +from numpy.testing import assert_raises import dpnp @@ -74,3 +75,21 @@ def test_fft_invalid_dtype(func_name): dpnp_func = getattr(dpnp.fft, func_name) with pytest.raises(TypeError): dpnp_func(a) + + +class TestFftfreq: + @pytest.mark.parametrize("func", ["fftfreq", "rfftfreq"]) + @pytest.mark.parametrize("n", [10, 20]) + @pytest.mark.parametrize("d", [0.5, 2]) + def test_fftfreq(self, func, n, d): + expected = getattr(dpnp.fft, func)(n, d) + result = getattr(numpy.fft, func)(n, d) + assert_dtype_allclose(expected, result) + + @pytest.mark.parametrize("func", ["fftfreq", "rfftfreq"]) + def test_error(self, func): + # n should be an integer + assert_raises(ValueError, getattr(dpnp.fft, func), 10.0) + + # d should be an scalar + assert_raises(ValueError, getattr(dpnp.fft, func), 10, (2,)) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 99334cfabfcd..eb3241b188a1 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1179,6 +1179,20 @@ def test_fft_rfft(type, shape, device): assert_sycl_queue_equal(result_queue, expected_queue) +@pytest.mark.parametrize("func", ["fftfreq", "rfftfreq"]) +@pytest.mark.parametrize( + "device", + valid_devices, + ids=[device.filter_string for device in valid_devices], +) +def test_fftfreq(func, device): + result = getattr(dpnp.fft, func)(10, 0.5, device=device) + expected = getattr(numpy.fft, func)(10, 0.5) + + assert_dtype_allclose(result, expected) + assert result.sycl_device == device + + @pytest.mark.parametrize( "data, is_empty", [ diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 4f7314ff2db0..eae1377e4616 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -914,6 +914,20 @@ def test_eigenvalue(func, shape, usm_type): assert a.usm_type == dp_val.usm_type +@pytest.mark.parametrize("func", ["fftfreq", "rfftfreq"]) +@pytest.mark.parametrize("usm_type", list_of_usm_types + [None]) +def test_fftfreq(func, usm_type): + result = getattr(dp.fft, func)(10, 0.5, usm_type=usm_type) + expected = getattr(numpy.fft, func)(10, 0.5) + + if usm_type is None: + # assert against default USM type + usm_type = "device" + + assert_dtype_allclose(result, expected) + assert result.usm_type == usm_type + + @pytest.mark.parametrize( "usm_type_matrix", list_of_usm_types, ids=list_of_usm_types ) diff --git a/tests/third_party/cupy/fft_tests/test_fft.py b/tests/third_party/cupy/fft_tests/test_fft.py index 4822c5bde87f..c6196236ccad 100644 --- a/tests/third_party/cupy/fft_tests/test_fft.py +++ b/tests/third_party/cupy/fft_tests/test_fft.py @@ -282,8 +282,7 @@ def test_ihfft(self, xp, dtype): {"n": 10, "d": 0.5}, {"n": 100, "d": 2}, ) -@pytest.mark.usefixtures("allow_fall_back_on_numpy") -class TestFftfreq(unittest.TestCase): +class TestFftfreq: @testing.for_all_dtypes() @testing.numpy_cupy_allclose( rtol=1e-4, From 8836bfb089c14753102f5f1558df6e7616b53e54 Mon Sep 17 00:00:00 2001 From: Vahid Tavanashad Date: Fri, 12 Jul 2024 23:09:13 -0500 Subject: [PATCH 2/4] fix pre-commit --- tests/test_fft.py | 4 ++-- tests/test_usm_type.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_fft.py b/tests/test_fft.py index abd089bda1d0..f6818a53214d 100644 --- a/tests/test_fft.py +++ b/tests/test_fft.py @@ -355,7 +355,7 @@ def test_fft_invalid_dtype(self, func_name): with pytest.raises(TypeError): dpnp_func(a) - + class TestFftfreq: @pytest.mark.parametrize("func", ["fftfreq", "rfftfreq"]) @pytest.mark.parametrize("n", [10, 20]) @@ -371,4 +371,4 @@ def test_error(self, func): assert_raises(ValueError, getattr(dpnp.fft, func), 10.0) # d should be an scalar - assert_raises(ValueError, getattr(dpnp.fft, func), 10, (2,)) \ No newline at end of file + assert_raises(ValueError, getattr(dpnp.fft, func), 10, (2,)) diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 4513831833bf..0febedbc6360 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -944,7 +944,7 @@ def test_fftfreq(func, usm_type): assert_dtype_allclose(result, expected) assert result.usm_type == usm_type - + @pytest.mark.parametrize("func", ["fft", "ifft"]) @pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) def test_fft(func, usm_type): From 815e266b596277e0200ff440151d178170dff0c7 Mon Sep 17 00:00:00 2001 From: vtavana <120411540+vtavana@users.noreply.github.com> Date: Sat, 13 Jul 2024 09:36:37 -0500 Subject: [PATCH 3/4] Apply suggestions from code review Co-authored-by: Anton <100830759+antonwolfy@users.noreply.github.com> --- dpnp/fft/dpnp_iface_fft.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/dpnp/fft/dpnp_iface_fft.py b/dpnp/fft/dpnp_iface_fft.py index 3bef71e07341..0fbd15d5cf02 100644 --- a/dpnp/fft/dpnp_iface_fft.py +++ b/dpnp/fft/dpnp_iface_fft.py @@ -218,7 +218,7 @@ def fftfreq(n, d=1.0, device=None, usm_type=None, sycl_queue=None): Window length. d : scalar, optional Sample spacing (inverse of the sampling rate). - Default: ``1``. + Default: ``1.0``. device : {None, string, SyclDevice, SyclQueue}, optional An array API concept of device where the output array is created. The `device` can be ``None`` (the default), an OneAPI filter selector @@ -226,6 +226,7 @@ def fftfreq(n, d=1.0, device=None, usm_type=None, sycl_queue=None): a non-partitioned SYCL device, an instance of :class:`dpctl.SyclQueue`, or a `Device` object returned by :obj:`dpnp.dpnp_array.dpnp_array.device` property. + Default: ``None``. usm_type : {None, "device", "shared", "host"}, optional The type of SYCL USM allocation for the output array. Default: ``None``. @@ -260,15 +261,15 @@ def fftfreq(n, d=1.0, device=None, usm_type=None, sycl_queue=None): >>> x = np.fft.fftfreq(n, d=timestep) # default case >>> x.shape, x.device, x.usm_type - ((8, ), Device(level_zero:gpu:0), 'device') + ((8,), Device(level_zero:gpu:0), 'device') >>> y = np.fft.fftfreq(n, d=timestep, device="cpu") >>> y.shape, y.device, y.usm_type - ((8, ), Device(opencl:cpu:0), 'device') + ((8,), Device(opencl:cpu:0), 'device') >>> z = np.fft.fftfreq(n, d=timestep, usm_type="host") >>> z.shape, z.device, z.usm_type - ((8, ), Device(level_zero:gpu:0), 'host') + ((8,), Device(level_zero:gpu:0), 'host') """ @@ -989,7 +990,7 @@ def rfftfreq(n, d=1.0, device=None, usm_type=None, sycl_queue=None): Window length. d : scalar, optional Sample spacing (inverse of the sampling rate). - Default: ``1``. + Default: ``1.0``. device : {None, string, SyclDevice, SyclQueue}, optional An array API concept of device where the output array is created. The `device` can be ``None`` (the default), an OneAPI filter selector @@ -997,6 +998,7 @@ def rfftfreq(n, d=1.0, device=None, usm_type=None, sycl_queue=None): a non-partitioned SYCL device, an instance of :class:`dpctl.SyclQueue`, or a `Device` object returned by :obj:`dpnp.dpnp_array.dpnp_array.device` property. + Default: ``None``. usm_type : {None, "device", "shared", "host"}, optional The type of SYCL USM allocation for the output array. Default: ``None``. @@ -1007,7 +1009,7 @@ def rfftfreq(n, d=1.0, device=None, usm_type=None, sycl_queue=None): Returns ------- f : dpnp.ndarray - Array of length `n` containing the sample frequencies. + Array of length ``n//2 + 1`` containing the sample frequencies. See Also -------- @@ -1027,22 +1029,22 @@ def rfftfreq(n, d=1.0, device=None, usm_type=None, sycl_queue=None): array([ 0., 10., 20., 30., 40., -50., -40., -30., -20., -10.]) >>> freq = np.fft.rfftfreq(n, d=1./sample_rate) >>> freq - array([ 0., 10., 20., 30., 40., 50.]) + array([ 0., 10., 20., 30., 40., 50.]) Creating the output array on a different device or with a specified usm_type: >>> x = np.fft.rfftfreq(n, d=1./sample_rate) # default case >>> x.shape, x.device, x.usm_type - ((6, ), Device(level_zero:gpu:0), 'device') + ((6,), Device(level_zero:gpu:0), 'device') - >>> y =np.fft.rfftfreq(n, d=1./sample_rate, device="cpu") + >>> y = np.fft.rfftfreq(n, d=1./sample_rate, device="cpu") >>> y.shape, y.device, y.usm_type - ((6, ), Device(opencl:cpu:0), 'device') + ((6,), Device(opencl:cpu:0), 'device') >>> z = np.fft.rfftfreq(n, d=1./sample_rate, usm_type="host") >>> z.shape, z.device, z.usm_type - ((6, ), Device(level_zero:gpu:0), 'host') + ((6,), Device(level_zero:gpu:0), 'host') """ From fd320121c6f439ae8ae9dc5ff9846a38e7ba8e4f Mon Sep 17 00:00:00 2001 From: Vahid Tavanashad Date: Sat, 13 Jul 2024 09:51:16 -0500 Subject: [PATCH 4/4] update cupy tests --- tests/third_party/cupy/fft_tests/test_fft.py | 39 ++++++++++++++++---- 1 file changed, 31 insertions(+), 8 deletions(-) diff --git a/tests/third_party/cupy/fft_tests/test_fft.py b/tests/third_party/cupy/fft_tests/test_fft.py index 71cd2b3961a7..7b354cdab278 100644 --- a/tests/third_party/cupy/fft_tests/test_fft.py +++ b/tests/third_party/cupy/fft_tests/test_fft.py @@ -10,6 +10,14 @@ from tests.third_party.cupy import testing +@pytest.fixture +def skip_forward_backward(request): + if request.instance.norm in ("backward", "forward"): + if not (np.lib.NumpyVersion(np.__version__) >= "1.20.0"): + pytest.skip("forward/backward is supported by NumPy 1.20+") + + +@pytest.mark.usefixtures("skip_forward_backward") @testing.parameterize( *testing.product( { @@ -26,12 +34,16 @@ class TestFft: atol=1e-7, accept_error=ValueError, contiguous_check=False, - type_check=False, + type_check=has_support_aspect64(), ) def test_fft(self, xp, dtype): a = testing.shaped_random(self.shape, xp, dtype) out = xp.fft.fft(a, n=self.n, norm=self.norm) + # np.fft.fft always returns np.complex128 + if xp is np and dtype in [np.float16, np.float32, np.complex64]: + out = out.astype(np.complex64) + return out @testing.for_all_dtypes() @@ -40,12 +52,18 @@ def test_fft(self, xp, dtype): atol=1e-7, accept_error=ValueError, contiguous_check=False, - type_check=False, + type_check=has_support_aspect64(), ) + # NumPy 1.17.0 and 1.17.1 raises ZeroDivisonError due to a bug + @testing.with_requires("numpy!=1.17.0") + @testing.with_requires("numpy!=1.17.1") def test_ifft(self, xp, dtype): a = testing.shaped_random(self.shape, xp, dtype) out = xp.fft.ifft(a, n=self.n, norm=self.norm) + if xp is np and dtype in [np.float16, np.float32, np.complex64]: + out = out.astype(np.complex64) + return out @@ -65,7 +83,7 @@ class TestFftOrder: atol=1e-6, accept_error=ValueError, contiguous_check=False, - type_check=False, + type_check=has_support_aspect64(), ) def test_fft(self, xp, dtype): a = testing.shaped_random(self.shape, xp, dtype) @@ -73,6 +91,10 @@ def test_fft(self, xp, dtype): a = xp.asfortranarray(a) out = xp.fft.fft(a, axis=self.axis) + # np.fft.fft always returns np.complex128 + if xp is np and dtype in [np.float16, np.float32, np.complex64]: + out = out.astype(np.complex64) + return out @testing.for_all_dtypes() @@ -81,7 +103,7 @@ def test_fft(self, xp, dtype): atol=1e-7, accept_error=ValueError, contiguous_check=False, - type_check=False, + type_check=has_support_aspect64(), ) def test_ifft(self, xp, dtype): a = testing.shaped_random(self.shape, xp, dtype) @@ -89,6 +111,9 @@ def test_ifft(self, xp, dtype): a = xp.asfortranarray(a) out = xp.fft.ifft(a, axis=self.axis) + if xp is np and dtype in [np.float16, np.float32, np.complex64]: + out = out.astype(np.complex64) + return out @@ -325,24 +350,22 @@ def test_ihfft(self, xp, dtype): {"n": 100, "d": 2}, ) class TestFftfreq: - @testing.for_all_dtypes() @testing.numpy_cupy_allclose( rtol=1e-4, atol=1e-7, type_check=has_support_aspect64(), ) - def test_fftfreq(self, xp, dtype): + def test_fftfreq(self, xp): out = xp.fft.fftfreq(self.n, self.d) return out - @testing.for_all_dtypes() @testing.numpy_cupy_allclose( rtol=1e-4, atol=1e-7, type_check=has_support_aspect64(), ) - def test_rfftfreq(self, xp, dtype): + def test_rfftfreq(self, xp): out = xp.fft.rfftfreq(self.n, self.d) return out