Implementation of convolve

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] Wikipedia, "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)

dpnp/tests/test_mathematical.py

@@ -104,35 +104,6 @@ def test_conj_out(self, dtype):
         assert_dtype_allclose(result, expected)
 
 
-@pytest.mark.usefixtures("allow_fall_back_on_numpy")
-class TestConvolve:
-    def test_object(self):
-        d = [1.0] * 100
-        k = [1.0] * 3
-        assert_array_equal(dpnp.convolve(d, k)[2:-2], dpnp.full(98, 3.0))
-
-    def test_no_overwrite(self):
-        d = dpnp.ones(100)
-        k = dpnp.ones(3)
-        dpnp.convolve(d, k)
-        assert_array_equal(d, dpnp.ones(100))
-        assert_array_equal(k, dpnp.ones(3))
-
-    def test_mode(self):
-        d = dpnp.ones(100)
-        k = dpnp.ones(3)
-        default_mode = dpnp.convolve(d, k)
-        full_mode = dpnp.convolve(d, k, mode="full")
-        assert_array_equal(full_mode, default_mode)
-        # integer mode
-        with assert_raises(ValueError):
-            dpnp.convolve(d, k, mode=-1)
-        assert_array_equal(dpnp.convolve(d, k, mode=2), full_mode)
-        # illegal arguments
-        with assert_raises(TypeError):
-            dpnp.convolve(d, k, mode=None)
-
-
 class TestClip:
     @pytest.mark.parametrize(
         "dtype", get_all_dtypes(no_bool=True, no_none=True, no_complex=True)