Merge branch 'master' into remove-workaraound-dpctl-2030

Author: vtavana

CHANGELOG.md

@@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Added implementation of `dpnp.hanning` [#2358](https://github.com/IntelPython/dpnp/pull/2358)
 * Added implementation of `dpnp.blackman` [#2363](https://github.com/IntelPython/dpnp/pull/2363)
 * Added implementation of `dpnp.bartlett` [#2366](https://github.com/IntelPython/dpnp/pull/2366)
+* Added implementation of `dpnp.convolve` [#2205](https://github.com/IntelPython/dpnp/pull/2205)
 * Added implementation of `dpnp.kaiser` [#2387](https://github.com/IntelPython/dpnp/pull/2387)
 
 ### Changed

dpnp/dpnp_iface_manipulation.py

@@ -958,6 +958,9 @@ def atleast_1d(*arys):
     dpnp.check_supported_arrays_type(*arys)
     for ary in arys:
         if ary.ndim == 0:
+            # 0-d arrays cannot be empty
+            # 0-d arrays always have a size of 1, so
+            # reshape(1) is guaranteed to succeed
             result = ary.reshape(1)
         else:
             result = ary

dpnp/dpnp_iface_mathematical.py

@@ -91,7 +91,6 @@
     "clip",
     "conj",
     "conjugate",
-    "convolve",
     "copysign",
     "cross",
     "cumprod",
@@ -792,24 +791,6 @@ def clip(a, /, min=None, max=None, *, out=None, order="K", **kwargs):
 
 conj = conjugate
 
-
-def convolve(a, v, mode="full"):
-    """
-    Returns the discrete, linear convolution of two one-dimensional sequences.
-
-    For full documentation refer to :obj:`numpy.convolve`.
-
-    Examples
-    --------
-    >>> ca = dpnp.convolve([1, 2, 3], [0, 1, 0.5])
-    >>> print(ca)
-    [0. , 1. , 2.5, 4. , 1.5]
-
-    """
-
-    return call_origin(numpy.convolve, a=a, v=v, mode=mode)
-
-
 _COPYSIGN_DOCSTRING = """
 Composes a floating-point value with the magnitude of `x1_i` and the sign of
 `x2_i` for each element of input arrays `x1` and `x2`.

dpnp/dpnp_iface_statistics.py

@@ -62,6 +62,7 @@
     "amax",
     "amin",
     "average",
+    "convolve",
     "corrcoef",
     "correlate",
     "cov",
@@ -357,6 +358,146 @@ def average(a, axis=None, weights=None, returned=False, *, keepdims=False):
     return avg
 
 
+def _convolve_impl(a, v, mode, method, rdtype):
+    l_pad, r_pad = _get_padding(a.size, v.size, mode)
+
+    if method == "auto":
+        method = _choose_conv_method(a, v, rdtype)
+
+    if method == "direct":
+        r = _run_native_sliding_dot_product1d(a, v[::-1], l_pad, r_pad, rdtype)
+    elif method == "fft":
+        r = _convolve_fft(a, v, l_pad, r_pad, rdtype)
+    else:
+        raise ValueError(
+            f"Unknown method: {method}. Supported methods: auto, direct, fft"
+        )
+
+    return r
+
+
+def convolve(a, v, mode="full", method="auto"):
+    r"""
+    Returns the discrete, linear convolution of two one-dimensional sequences.
+    The convolution operator is often seen in signal processing, where it
+    models the effect of a linear time-invariant system on a signal [1]_. In
+    probability theory, the sum of two independent random variables is
+    distributed according to the convolution of their individual
+    distributions.
+
+    If `v` is longer than `a`, the arrays are swapped before computation.
+
+    For full documentation refer to :obj:`numpy.convolve`.
+
+    Parameters
+    ----------
+    a : {dpnp.ndarray, usm_ndarray}
+        First input array.
+    v : {dpnp.ndarray, usm_ndarray}
+        Second input array.
+    mode : {'full', 'valid', 'same'}, optional
+        - 'full': This returns the convolution
+          at each point of overlap, with an output shape of (N+M-1,). At
+          the end-points of the convolution, the signals do not overlap
+          completely, and boundary effects may be seen.
+        - 'same': Mode 'same' returns output of length ``max(M, N)``. Boundary
+          effects are still visible.
+        - 'valid': Mode 'valid' returns output of length
+          ``max(M, N) - min(M, N) + 1``. The convolution product is only given
+          for points where the signals overlap completely. Values outside
+          the signal boundary have no effect.
+
+          Default: ``'full'``.
+    method : {'auto', 'direct', 'fft'}, optional
+        - 'direct': The convolution is determined directly from sums.
+        - 'fft': The Fourier Transform is used to perform the calculations.
+         This method is faster for long sequences but can have accuracy issues.
+        - 'auto': Automatically chooses direct or Fourier method based on
+         an estimate of which is faster.
+
+        Note: Use of the FFT convolution on input containing NAN or INF
+        will lead to the entire output being NAN or INF.
+        Use ``method='direct'`` when your input contains NAN or INF values.
+
+        Default: ``'auto'``.
+
+    Returns
+    -------
+    out : dpnp.ndarray
+        Discrete, linear convolution of `a` and `v`.
+
+    See Also
+    --------
+    :obj:`dpnp.correlate` : Cross-correlation of two 1-dimensional sequences.
+
+    Notes
+    -----
+    The discrete convolution operation is defined as
+
+    .. math:: (a * v)_n = \sum_{m = -\infty}^{\infty} a_m v_{n - m}
+
+    It can be shown that a convolution :math:`x(t) * y(t)` in time/space
+    is equivalent to the multiplication :math:`X(f) Y(f)` in the Fourier
+    domain, after appropriate padding (padding is necessary to prevent
+    circular convolution). Since multiplication is more efficient (faster)
+    than convolution, the function implements two approaches - direct and fft
+    which are regulated by the keyword `method`.
+
+    References
+    ----------
+    .. [1] Uncyclopedia, "Convolution",
+        https://en.wikipedia.org/wiki/Convolution
+
+    Examples
+    --------
+    Note how the convolution operator flips the second array
+    before "sliding" the two across one another:
+
+    >>> import dpnp as np
+    >>> a = np.array([1, 2, 3], dtype=np.float32)
+    >>> v = np.array([0, 1, 0.5], dtype=np.float32)
+    >>> np.convolve(a, v)
+    array([0. , 1. , 2.5, 4. , 1.5], dtype=float32)
+
+    Only return the middle values of the convolution.
+    Contains boundary effects, where zeros are taken
+    into account:
+
+    >>> np.convolve(a, v, 'same')
+    array([1. , 2.5, 4. ], dtype=float32)
+
+    The two arrays are of the same length, so there
+    is only one position where they completely overlap:
+
+    >>> np.convolve(a, v, 'valid')
+    array([2.5], dtype=float32)
+
+    """
+
+    a, v = dpnp.atleast_1d(a, v)
+
+    if a.size == 0 or v.size == 0:
+        raise ValueError(
+            f"Array arguments cannot be empty. "
+            f"Received sizes: a.size={a.size}, v.size={v.size}"
+        )
+    if a.ndim > 1 or v.ndim > 1:
+        raise ValueError(
+            f"Only 1-dimensional arrays are supported. "
+            f"Received shapes: a.shape={a.shape}, v.shape={v.shape}"
+        )
+
+    device = a.sycl_device
+    rdtype = result_type_for_device([a.dtype, v.dtype], device)
+
+    if v.size > a.size:
+        a, v = v, a
+
+    r = _convolve_impl(a, v, mode, method, rdtype)
+
+    return dpnp.asarray(r, dtype=rdtype, order="C")
+
+
 def corrcoef(x, y=None, rowvar=True, *, dtype=None):
     """
     Return Pearson product-moment correlation coefficients.
@@ -714,17 +855,7 @@ def correlate(a, v, mode="valid", method="auto"):
         revert = True
         a, v = v, a
 
-    l_pad, r_pad = _get_padding(a.size, v.size, mode)
-
-    if method == "auto":
-        method = _choose_conv_method(a, v, rdtype)
-
-    if method == "direct":
-        r = _run_native_sliding_dot_product1d(a, v, l_pad, r_pad, rdtype)
-    elif method == "fft":
-        r = _convolve_fft(a, v[::-1], l_pad, r_pad, rdtype)
-    else:  # pragma: no cover
-        raise ValueError(f"Unknown method: {method}")
+    r = _convolve_impl(a, v[::-1], mode, method, rdtype)
 
     if revert:
         r = r[::-1]

dpnp/tests/helper.py

@@ -185,6 +185,19 @@ def generate_random_numpy_array(
     return a
 
 
+def factor_to_tol(dtype, factor):
+    """
+    Calculate the tolerance for comparing floating point and complex arrays.
+    The tolerance is based on the maximum resolution of the input dtype multiplied by the factor.
+    """
+
+    tol = 0
+    if numpy.issubdtype(dtype, numpy.inexact):
+        tol = numpy.finfo(dtype).resolution
+
+    return factor * tol
+
+
 def get_abs_array(data, dtype=None):
     if numpy.issubdtype(dtype, numpy.unsignedinteger):
         data = numpy.abs(data)
@@ -312,6 +325,36 @@ def get_integer_dtypes(all_int_types=False, no_unsigned=False):
     return dtypes
 
 
+def get_integer_float_dtypes(
+    all_int_types=False,
+    no_unsigned=False,
+    no_float16=True,
+    device=None,
+    xfail_dtypes=None,
+    exclude=None,
+):
+    """
+    Build a list of integer and float types supported by DPNP.
+    """
+    dtypes = get_integer_dtypes(
+        all_int_types=all_int_types, no_unsigned=no_unsigned
+    )
+    dtypes += get_float_dtypes(no_float16=no_float16, device=device)
+
+    def mark_xfail(dtype):
+        if xfail_dtypes is not None and dtype in xfail_dtypes:
+            return pytest.param(dtype, marks=pytest.mark.xfail)
+        return dtype
+
+    def not_excluded(dtype):
+        if exclude is None:
+            return True
+        return dtype not in exclude
+
+    dtypes = [mark_xfail(dtype) for dtype in dtypes if not_excluded(dtype)]
+    return dtypes
+
+
 def has_support_aspect16(device=None):
     """
     Return True if the device supports 16-bit precision floating point operations,

dpnp/tests/test_binary_ufuncs.py

@@ -20,6 +20,7 @@
     get_float_complex_dtypes,
     get_float_dtypes,
     get_integer_dtypes,
+    get_integer_float_dtypes,
     has_support_aspect16,
     numpy_version,
 )
@@ -141,7 +142,7 @@ def test_invalid_out(self, xp, out):
 @pytest.mark.parametrize("func", ["fmax", "fmin", "maximum", "minimum"])
 class TestBoundFuncs:
     @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True))
-    def test_out(self, func, dtype):
+    def test_basic(self, func, dtype):
         a = generate_random_numpy_array(10, dtype)
         b = generate_random_numpy_array(10, dtype)
         expected = getattr(numpy, func)(a, b)
@@ -278,7 +279,7 @@ def test_invalid_out(self, xp, out):
 
 @pytest.mark.parametrize("func", ["floor_divide", "remainder"])
 class TestFloorDivideRemainder:
-    ALL_DTYPES = get_all_dtypes(no_none=True, no_bool=True, no_complex=True)
+    ALL_DTYPES = get_integer_float_dtypes()
 
     def do_inplace_op(self, base, other, func):
         if func == "floor_divide":

dpnp/tests/test_histogram.py

@@ -19,15 +19,14 @@
     get_float_complex_dtypes,
     get_float_dtypes,
     get_integer_dtypes,
+    get_integer_float_dtypes,
     has_support_aspect64,
     numpy_version,
 )
 
 
 class TestDigitize:
-    @pytest.mark.parametrize(
-        "dtype", get_all_dtypes(no_bool=True, no_complex=True)
-    )
+    @pytest.mark.parametrize("dtype", get_integer_float_dtypes())
     @pytest.mark.parametrize("right", [True, False])
     @pytest.mark.parametrize(
         "x, bins",
@@ -73,12 +72,8 @@ def test_digitize_inf(self, dtype, right):
         expected = numpy.digitize(x, bins, right=right)
         assert_dtype_allclose(result, expected)
 
-    @pytest.mark.parametrize(
-        "dtype_x", get_all_dtypes(no_bool=True, no_complex=True)
-    )
-    @pytest.mark.parametrize(
-        "dtype_bins", get_all_dtypes(no_bool=True, no_complex=True)
-    )
+    @pytest.mark.parametrize("dtype_x", get_integer_float_dtypes())
+    @pytest.mark.parametrize("dtype_bins", get_integer_float_dtypes())
     @pytest.mark.parametrize("right", [True, False])
     def test_digitize_diff_types(self, dtype_x, dtype_bins, right):
         x = numpy.array([1, 2, 3, 4, 5], dtype=dtype_x)
@@ -90,9 +85,7 @@ def test_digitize_diff_types(self, dtype_x, dtype_bins, right):
         expected = numpy.digitize(x, bins, right=right)
         assert_dtype_allclose(result, expected)
 
-    @pytest.mark.parametrize(
-        "dtype", get_all_dtypes(no_bool=True, no_complex=True)
-    )
+    @pytest.mark.parametrize("dtype", get_integer_float_dtypes())
     @pytest.mark.parametrize(
         "x, bins",
         [

dpnp/tests/test_linalg.py

@@ -22,6 +22,7 @@
     get_all_dtypes,
     get_complex_dtypes,
     get_float_complex_dtypes,
+    get_integer_float_dtypes,
     has_support_aspect64,
     is_cpu_device,
     is_cuda_device,
@@ -1409,9 +1410,7 @@ def test_einsum_tensor(self):
         result = dpnp.einsum("ijij->", tensor_dp)
         assert_dtype_allclose(result, expected)
 
-    @pytest.mark.parametrize(
-        "dtype", get_all_dtypes(no_bool=True, no_complex=True, no_none=True)
-    )
+    @pytest.mark.parametrize("dtype", get_integer_float_dtypes())
     def test_different_paths(self, dtype):
         # Simple test, designed to exercise most specialized code paths,
         # note the +0.5 for floats.  This makes sure we use a float value