Rebase on master

Author: AlexanderKalistratov

dpnp/dpnp_iface_mathematical.py

@@ -90,7 +90,6 @@
     "clip",
     "conj",
     "conjugate",
-    "convolve",
     "copysign",
     "cross",
     "cumprod",
@@ -791,24 +790,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,138 @@ 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 1-D array.
+    v : {dpnp.ndarray, usm_ndarray}
+        Second 1-D array. The length of `v` must be less than or equal to
+        the length of `a`.
+    mode : {'full', 'valid', 'same'}, optional
+        'full':
+          By default, mode is '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.
+    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 : 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)
+    """
+
+    dpnp.check_supported_arrays_type(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}"
+        )
+
+    if a.ndim == 0:
+        a = dpnp.reshape(a, (1,))
+    if v.ndim == 0:
+        v = dpnp.reshape(v, (1,))
+
+    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 +847,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/test_mathematical.py

@@ -110,35 +110,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_almost_equal(dpnp.convolve(d, k)[2:-2], dpnp.full(98, 3))
-
-    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, mode="full")
-        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)

dpnp/tests/test_statistics.py

@@ -127,6 +127,164 @@ def test_avg_error(self):
             dpnp.average(a, axis=0, weights=w)
 
 
+class TestConvolve:
+    def setup_method(self):
+        numpy.random.seed(0)
+
+    @pytest.mark.parametrize(
+        "a, v", [([1], [1, 2, 3]), ([1, 2, 3], [1]), ([1, 2, 3], [1, 2])]
+    )
+    @pytest.mark.parametrize("mode", [None, "full", "valid", "same"])
+    @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True))
+    @pytest.mark.parametrize("method", [None, "auto", "direct", "fft"])
+    def test_convolve(self, a, v, mode, dtype, method):
+        an = numpy.array(a, dtype=dtype)
+        vn = numpy.array(v, dtype=dtype)
+        ad = dpnp.array(an)
+        vd = dpnp.array(vn)
+
+        dpnp_kwargs = {}
+        numpy_kwargs = {}
+        if mode is not None:
+            dpnp_kwargs["mode"] = mode
+            numpy_kwargs["mode"] = mode
+        if method is not None:
+            dpnp_kwargs["method"] = method
+
+        expected = numpy.convolve(an, vn, **numpy_kwargs)
+        result = dpnp.convolve(ad, vd, **dpnp_kwargs)
+
+        assert_dtype_allclose(result, expected)
+
+    @pytest.mark.parametrize("a_size", [1, 100, 10000])
+    @pytest.mark.parametrize("v_size", [1, 100, 10000])
+    @pytest.mark.parametrize("mode", ["full", "valid", "same"])
+    @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True))
+    @pytest.mark.parametrize("method", ["auto", "direct", "fft"])
+    def test_convolve_random(self, a_size, v_size, mode, dtype, method):
+        if dtype == dpnp.bool:
+            an = numpy.random.rand(a_size) > 0.9
+            vn = numpy.random.rand(v_size) > 0.9
+        else:
+            an = (100 * numpy.random.rand(a_size)).astype(dtype)
+            vn = (100 * numpy.random.rand(v_size)).astype(dtype)
+
+            if dpnp.issubdtype(dtype, dpnp.complexfloating):
+                an = an + 1j * (100 * numpy.random.rand(a_size)).astype(dtype)
+                vn = vn + 1j * (100 * numpy.random.rand(v_size)).astype(dtype)
+
+        ad = dpnp.array(an)
+        vd = dpnp.array(vn)
+
+        dpnp_kwargs = {}
+        numpy_kwargs = {}
+        if mode is not None:
+            dpnp_kwargs["mode"] = mode
+            numpy_kwargs["mode"] = mode
+        if method is not None:
+            dpnp_kwargs["method"] = method
+
+        result = dpnp.convolve(ad, vd, **dpnp_kwargs)
+        expected = numpy.convolve(an, vn, **numpy_kwargs)
+
+        rdtype = result.dtype
+        if dpnp.issubdtype(rdtype, dpnp.integer):
+            rdtype = dpnp.default_float_type(ad.device)
+
+        if method != "fft" and (
+            dpnp.issubdtype(dtype, dpnp.integer) or dtype == dpnp.bool
+        ):
+            # For 'direct' and 'auto' methods, we expect exact results for integer types
+            assert_array_equal(result, expected)
+        else:
+            result = result.astype(rdtype)
+            if method == "direct":
+                expected = numpy.convolve(an, vn, **numpy_kwargs)
+                # For 'direct' method we can use standard validation
+                assert_dtype_allclose(result, expected, factor=30)
+            else:
+                rtol = 1e-3
+                atol = 1e-10
+
+                if rdtype == dpnp.float64 or rdtype == dpnp.complex128:
+                    rtol = 1e-6
+                    atol = 1e-12
+                elif rdtype == dpnp.bool:
+                    result = result.astype(dpnp.int32)
+                    rdtype = result.dtype
+
+                expected = expected.astype(rdtype)
+
+                diff = numpy.abs(result.asnumpy() - expected)
+                invalid = diff > atol + rtol * numpy.abs(expected)
+
+                # When using the 'fft' method, we might encounter outliers.
+                # This usually happens when the resulting array contains values close to zero.
+                # For these outliers, the relative error can be significant.
+                # We can tolerate a few such outliers.
+                max_outliers = 8 if expected.size > 1 else 0
+                if invalid.sum() > max_outliers:
+                    assert_dtype_allclose(result, expected, factor=1000)
+
+    def test_convolve_mode_error(self):
+        a = dpnp.arange(5)
+        v = dpnp.arange(3)
+
+        # invalid mode
+        with pytest.raises(ValueError):
+            dpnp.convolve(a, v, mode="unknown")
+
+    @pytest.mark.parametrize("a, v", [([], [1]), ([1], []), ([], [])])
+    def test_convolve_empty(self, a, v):
+        a = dpnp.asarray(a)
+        v = dpnp.asarray(v)
+
+        with pytest.raises(ValueError):
+            dpnp.convolve(a, v)
+
+    @pytest.mark.parametrize(
+        "a, v",
+        [
+            ([[1, 2], [2, 3]], [1]),
+            ([1], [[1, 2], [2, 3]]),
+            ([[1, 2], [2, 3]], [[1, 2], [2, 3]]),
+        ],
+    )
+    def test_convolve_shape_error(self, a, v):
+        a = dpnp.asarray(a)
+        v = dpnp.asarray(v)
+
+        with pytest.raises(ValueError):
+            dpnp.convolve(a, v)
+
+    @pytest.mark.parametrize("size", [2, 10**1, 10**2, 10**3, 10**4, 10**5])
+    def test_convolve_different_sizes(self, size):
+        a = numpy.random.rand(size).astype(numpy.float32)
+        v = numpy.random.rand(size // 2).astype(numpy.float32)
+
+        ad = dpnp.array(a)
+        vd = dpnp.array(v)
+
+        expected = numpy.convolve(a, v)
+        result = dpnp.convolve(ad, vd, method="direct")
+
+        assert_dtype_allclose(result, expected, factor=20)
+
+    def test_convolve_another_sycl_queue(self):
+        a = dpnp.arange(5, sycl_queue=dpctl.SyclQueue())
+        v = dpnp.arange(3, sycl_queue=dpctl.SyclQueue())
+
+        with pytest.raises(ValueError):
+            dpnp.convolve(a, v)
+
+    def test_convolve_unkown_method(self):
+        a = dpnp.arange(5)
+        v = dpnp.arange(3)
+
+        with pytest.raises(ValueError):
+            dpnp.convolve(a, v, method="unknown")
+
+
 class TestCorrcoef:
     @pytest.mark.usefixtures(
         "suppress_divide_invalid_numpy_warnings",