|
62 | 62 | "amax", |
63 | 63 | "amin", |
64 | 64 | "average", |
| 65 | + "convolve", |
65 | 66 | "corrcoef", |
66 | 67 | "correlate", |
67 | 68 | "cov", |
@@ -357,6 +358,146 @@ def average(a, axis=None, weights=None, returned=False, *, keepdims=False): |
357 | 358 | return avg |
358 | 359 |
|
359 | 360 |
|
| 361 | +def _convolve_impl(a, v, mode, method, rdtype): |
| 362 | + l_pad, r_pad = _get_padding(a.size, v.size, mode) |
| 363 | + |
| 364 | + if method == "auto": |
| 365 | + method = _choose_conv_method(a, v, rdtype) |
| 366 | + |
| 367 | + if method == "direct": |
| 368 | + r = _run_native_sliding_dot_product1d(a, v[::-1], l_pad, r_pad, rdtype) |
| 369 | + elif method == "fft": |
| 370 | + r = _convolve_fft(a, v, l_pad, r_pad, rdtype) |
| 371 | + else: |
| 372 | + raise ValueError( |
| 373 | + f"Unknown method: {method}. Supported methods: auto, direct, fft" |
| 374 | + ) |
| 375 | + |
| 376 | + return r |
| 377 | + |
| 378 | + |
| 379 | +def convolve(a, v, mode="full", method="auto"): |
| 380 | + r""" |
| 381 | + Returns the discrete, linear convolution of two one-dimensional sequences. |
| 382 | + The convolution operator is often seen in signal processing, where it |
| 383 | + models the effect of a linear time-invariant system on a signal [1]_. In |
| 384 | + probability theory, the sum of two independent random variables is |
| 385 | + distributed according to the convolution of their individual |
| 386 | + distributions. |
| 387 | +
|
| 388 | + If `v` is longer than `a`, the arrays are swapped before computation. |
| 389 | +
|
| 390 | + For full documentation refer to :obj:`numpy.convolve`. |
| 391 | +
|
| 392 | + Parameters |
| 393 | + ---------- |
| 394 | + a : {dpnp.ndarray, usm_ndarray} |
| 395 | + First input array. |
| 396 | + v : {dpnp.ndarray, usm_ndarray} |
| 397 | + Second input array. |
| 398 | + mode : {'full', 'valid', 'same'}, optional |
| 399 | + - 'full': This returns the convolution |
| 400 | + at each point of overlap, with an output shape of (N+M-1,). At |
| 401 | + the end-points of the convolution, the signals do not overlap |
| 402 | + completely, and boundary effects may be seen. |
| 403 | + - 'same': Mode 'same' returns output of length ``max(M, N)``. Boundary |
| 404 | + effects are still visible. |
| 405 | + - 'valid': Mode 'valid' returns output of length |
| 406 | + ``max(M, N) - min(M, N) + 1``. The convolution product is only given |
| 407 | + for points where the signals overlap completely. Values outside |
| 408 | + the signal boundary have no effect. |
| 409 | +
|
| 410 | + Default: ``'full'``. |
| 411 | + method : {'auto', 'direct', 'fft'}, optional |
| 412 | + - 'direct': The convolution is determined directly from sums. |
| 413 | + - 'fft': The Fourier Transform is used to perform the calculations. |
| 414 | + This method is faster for long sequences but can have accuracy issues. |
| 415 | + - 'auto': Automatically chooses direct or Fourier method based on |
| 416 | + an estimate of which is faster. |
| 417 | +
|
| 418 | + Note: Use of the FFT convolution on input containing NAN or INF |
| 419 | + will lead to the entire output being NAN or INF. |
| 420 | + Use ``method='direct'`` when your input contains NAN or INF values. |
| 421 | +
|
| 422 | + Default: ``'auto'``. |
| 423 | +
|
| 424 | + Returns |
| 425 | + ------- |
| 426 | + out : dpnp.ndarray |
| 427 | + Discrete, linear convolution of `a` and `v`. |
| 428 | +
|
| 429 | + See Also |
| 430 | + -------- |
| 431 | + :obj:`dpnp.correlate` : Cross-correlation of two 1-dimensional sequences. |
| 432 | +
|
| 433 | + Notes |
| 434 | + ----- |
| 435 | + The discrete convolution operation is defined as |
| 436 | +
|
| 437 | + .. math:: (a * v)_n = \sum_{m = -\infty}^{\infty} a_m v_{n - m} |
| 438 | +
|
| 439 | + It can be shown that a convolution :math:`x(t) * y(t)` in time/space |
| 440 | + is equivalent to the multiplication :math:`X(f) Y(f)` in the Fourier |
| 441 | + domain, after appropriate padding (padding is necessary to prevent |
| 442 | + circular convolution). Since multiplication is more efficient (faster) |
| 443 | + than convolution, the function implements two approaches - direct and fft |
| 444 | + which are regulated by the keyword `method`. |
| 445 | +
|
| 446 | + References |
| 447 | + ---------- |
| 448 | + .. [1] Uncyclopedia, "Convolution", |
| 449 | + https://en.wikipedia.org/wiki/Convolution |
| 450 | +
|
| 451 | + Examples |
| 452 | + -------- |
| 453 | + Note how the convolution operator flips the second array |
| 454 | + before "sliding" the two across one another: |
| 455 | +
|
| 456 | + >>> import dpnp as np |
| 457 | + >>> a = np.array([1, 2, 3], dtype=np.float32) |
| 458 | + >>> v = np.array([0, 1, 0.5], dtype=np.float32) |
| 459 | + >>> np.convolve(a, v) |
| 460 | + array([0. , 1. , 2.5, 4. , 1.5], dtype=float32) |
| 461 | +
|
| 462 | + Only return the middle values of the convolution. |
| 463 | + Contains boundary effects, where zeros are taken |
| 464 | + into account: |
| 465 | +
|
| 466 | + >>> np.convolve(a, v, 'same') |
| 467 | + array([1. , 2.5, 4. ], dtype=float32) |
| 468 | +
|
| 469 | + The two arrays are of the same length, so there |
| 470 | + is only one position where they completely overlap: |
| 471 | +
|
| 472 | + >>> np.convolve(a, v, 'valid') |
| 473 | + array([2.5], dtype=float32) |
| 474 | +
|
| 475 | + """ |
| 476 | + |
| 477 | + a, v = dpnp.atleast_1d(a, v) |
| 478 | + |
| 479 | + if a.size == 0 or v.size == 0: |
| 480 | + raise ValueError( |
| 481 | + f"Array arguments cannot be empty. " |
| 482 | + f"Received sizes: a.size={a.size}, v.size={v.size}" |
| 483 | + ) |
| 484 | + if a.ndim > 1 or v.ndim > 1: |
| 485 | + raise ValueError( |
| 486 | + f"Only 1-dimensional arrays are supported. " |
| 487 | + f"Received shapes: a.shape={a.shape}, v.shape={v.shape}" |
| 488 | + ) |
| 489 | + |
| 490 | + device = a.sycl_device |
| 491 | + rdtype = result_type_for_device([a.dtype, v.dtype], device) |
| 492 | + |
| 493 | + if v.size > a.size: |
| 494 | + a, v = v, a |
| 495 | + |
| 496 | + r = _convolve_impl(a, v, mode, method, rdtype) |
| 497 | + |
| 498 | + return dpnp.asarray(r, dtype=rdtype, order="C") |
| 499 | + |
| 500 | + |
360 | 501 | def corrcoef(x, y=None, rowvar=True, *, dtype=None): |
361 | 502 | """ |
362 | 503 | Return Pearson product-moment correlation coefficients. |
@@ -714,17 +855,7 @@ def correlate(a, v, mode="valid", method="auto"): |
714 | 855 | revert = True |
715 | 856 | a, v = v, a |
716 | 857 |
|
717 | | - l_pad, r_pad = _get_padding(a.size, v.size, mode) |
718 | | - |
719 | | - if method == "auto": |
720 | | - method = _choose_conv_method(a, v, rdtype) |
721 | | - |
722 | | - if method == "direct": |
723 | | - r = _run_native_sliding_dot_product1d(a, v, l_pad, r_pad, rdtype) |
724 | | - elif method == "fft": |
725 | | - r = _convolve_fft(a, v[::-1], l_pad, r_pad, rdtype) |
726 | | - else: # pragma: no cover |
727 | | - raise ValueError(f"Unknown method: {method}") |
| 858 | + r = _convolve_impl(a, v[::-1], mode, method, rdtype) |
728 | 859 |
|
729 | 860 | if revert: |
730 | 861 | r = r[::-1] |
|
0 commit comments