Merge pull request #478 from grlee77/pad

ENH: update the pad utility function and make it part of the public API

Author: rgommers

doc/source/conf.py

@@ -36,6 +36,7 @@
 extensions = ['sphinx.ext.doctest', 'sphinx.ext.autodoc', 'sphinx.ext.todo',
               'sphinx.ext.extlinks', 'sphinx.ext.mathjax',
               'sphinx.ext.autosummary', 'numpydoc',
+              'sphinx.ext.intersphinx',
               'matplotlib.sphinxext.plot_directive']
 
 # Add any paths that contain templates here, relative to this directory.
@@ -224,3 +225,9 @@
 plot_formats = [('png', 96), 'pdf']
 plot_html_show_formats = False
 plot_html_show_source_link = False
+
+# -- Options for intersphinx extension ---------------------------------------
+
+# Intersphinx to get Numpy and other targets
+intersphinx_mapping = {
+    'numpy': ('https://docs.scipy.org/doc/numpy/', None)}

doc/source/pyplots/plot_boundary_modes.py

@@ -5,7 +5,7 @@
 
 In practice, which signal extension mode is beneficial will depend on the
 signal characteristics.  For this particular signal, some modes such as
-"periodic",  "antisymmetric" and "zeros" result in large discontinuities that
+"periodic",  "antisymmetric" and "zero" result in large discontinuities that
 would lead to large amplitude boundary coefficients in the detail coefficients
 of a discrete wavelet transform.
 """
@@ -28,5 +28,5 @@
 boundary_mode_subplot(x, 'periodization', axes[5], symw=False)
 boundary_mode_subplot(x, 'smooth', axes[6], symw=False)
 boundary_mode_subplot(x, 'constant', axes[7], symw=False)
-boundary_mode_subplot(x, 'zeros', axes[8], symw=False)
+boundary_mode_subplot(x, 'zero', axes[8], symw=False)
 plt.show()

doc/source/ref/signal-extension-modes.rst

@@ -136,3 +136,15 @@ periodization      per           N/A
 antisymmetric      asym, asymh   N/A
 antireflect        asymw         reflect, reflect_type='odd'
 ================== ============= ===========================
+
+Padding using PyWavelets Signal Extension Modes - ``pad``
+---------------------------------------------------------
+
+.. autofunction:: pad
+
+Pywavelets provides a function, :func:`pad`, that operate like
+:func:`numpy.pad`, but supporting the PyWavelets signal extension modes
+discussed above. For efficiency, the DWT routines in PyWavelets do not
+expclitly create padded signals using this function. It can be used to manually
+prepad signals to reduce boundary effects in functions such as :func:`cwt` and
+:func:`swt` that do not currently support all of these signal extension modes.

pywt/_doc_utils.py

@@ -4,8 +4,10 @@
 import numpy as np
 from matplotlib import pyplot as plt
 
+from ._dwt import pad
+
 __all__ = ['wavedec_keys', 'wavedec2_keys', 'draw_2d_wp_basis',
-           'draw_2d_fswavedecn_basis', 'pad', 'boundary_mode_subplot']
+           'draw_2d_fswavedecn_basis', 'boundary_mode_subplot']
 
 
 def wavedec_keys(level):
@@ -149,63 +151,6 @@ def draw_2d_fswavedecn_basis(shape, levels, fmt='k', plot_kwargs={}, ax=None,
     return fig, ax
 
 
-def pad(x, pad_widths, mode):
-    """Extend a 1D signal using a given boundary mode.
-
-    Like numpy.pad but supports all PyWavelets boundary modes.
-    """
-    if np.isscalar(pad_widths):
-        pad_widths = (pad_widths, pad_widths)
-
-    if x.ndim > 1:
-        raise ValueError("This padding function is only for 1D signals.")
-
-    if mode in ['symmetric', 'reflect']:
-        xp = np.pad(x, pad_widths, mode=mode)
-    elif mode in ['periodic', 'periodization']:
-        if mode == 'periodization' and x.size % 2 == 1:
-            raise ValueError("periodization expects an even length signal.")
-        xp = np.pad(x, pad_widths, mode='wrap')
-    elif mode == 'zeros':
-        xp = np.pad(x, pad_widths, mode='constant', constant_values=0)
-    elif mode == 'constant':
-        xp = np.pad(x, pad_widths, mode='edge')
-    elif mode == 'smooth':
-        xp = np.pad(x, pad_widths, mode='linear_ramp',
-                    end_values=(x[0] + pad_widths[0] * (x[0] - x[1]),
-                                x[-1] + pad_widths[1] * (x[-1] - x[-2])))
-    elif mode == 'antisymmetric':
-        # implement by flipping portions symmetric padding
-        npad_l, npad_r = pad_widths
-        xp = np.pad(x, pad_widths, mode='symmetric')
-        r_edge = npad_l + x.size - 1
-        l_edge = npad_l
-        # width of each reflected segment
-        seg_width = x.size
-        # flip reflected segments on the right of the original signal
-        n = 1
-        while r_edge <= xp.size:
-            segment_slice = slice(r_edge + 1,
-                                  min(r_edge + 1 + seg_width, xp.size))
-            if n % 2:
-                xp[segment_slice] *= -1
-            r_edge += seg_width
-            n += 1
-
-        # flip reflected segments on the left of the original signal
-        n = 1
-        while l_edge >= 0:
-            segment_slice = slice(max(0, l_edge - seg_width), l_edge)
-            if n % 2:
-                xp[segment_slice] *= -1
-            l_edge -= seg_width
-            n += 1
-    elif mode == 'antireflect':
-        npad_l, npad_r = pad_widths
-        xp = np.pad(x, pad_widths, mode='reflect', reflect_type='odd')
-    return xp
-
-
 def boundary_mode_subplot(x, mode, ax, symw=True):
     """Plot an illustration of the boundary mode in a subplot axis."""
 
@@ -236,7 +181,7 @@ def boundary_mode_subplot(x, mode, ax, symw=True):
         left -= 0.5
         step = len(x)
         rng = range(-2, 4)
-    if mode in ['smooth', 'constant', 'zeros']:
+    if mode in ['smooth', 'constant', 'zero']:
         rng = range(0, 2)
     for rep in rng:
         ax.plot((left + rep * step) * o2, [xp.min() - .5, xp.max() + .5], 'k-')

pywt/_dwt.py

@@ -12,7 +12,7 @@
 
 
 __all__ = ["dwt", "idwt", "downcoef", "upcoef", "dwt_max_level",
-           "dwt_coeff_len"]
+           "dwt_coeff_len", "pad"]
 
 
 def dwt_max_level(data_len, filter_len):
@@ -135,7 +135,6 @@ def dwt(data, wavelet, mode='symmetric', axis=-1):
         Axis over which to compute the DWT. If not given, the
         last axis is used.
 
-
     Returns
     -------
     (cA, cD) : tuple
@@ -211,7 +210,6 @@ def idwt(cA, cD, wavelet, mode='symmetric', axis=-1):
         Axis over which to compute the inverse DWT. If not given, the
         last axis is used.
 
-
     Returns
     -------
     rec: array_like
@@ -401,3 +399,119 @@ def upcoef(part, coeffs, wavelet, level=1, take=0):
     if part not in 'ad':
         raise ValueError("Argument 1 must be 'a' or 'd', not '%s'." % part)
     return np.asarray(_upcoef(part == 'a', coeffs, wavelet, level, take))
+
+
+def pad(x, pad_widths, mode):
+    """Extend a 1D signal using a given boundary mode.
+
+    This function operates like :func:`numpy.pad` but supports all signal
+    extension modes that can be used by PyWavelets discrete wavelet transforms.
+
+    Parameters
+    ----------
+    x : ndarray
+        The array to pad
+    pad_widths : {sequence, array_like, int}
+        Number of values padded to the edges of each axis.
+        ``((before_1, after_1), … (before_N, after_N))`` unique pad widths for
+        each axis. ``((before, after),)`` yields same before and after pad for
+        each axis. ``(pad,)`` or int is a shortcut for
+        ``before = after = pad width`` for all axes.
+    mode : str, optional
+        Signal extension mode, see :ref:`Modes <ref-modes>`.
+
+    Returns
+    -------
+    pad : ndarray
+        Padded array of rank equal to array with shape increased according to
+        ``pad_widths``.
+
+    Notes
+    -----
+    The performance of padding in dimensions > 1 may be substantially slower
+    for modes ``'smooth'`` and ``'antisymmetric'`` as these modes are not
+    supported efficiently by the underlying :func:`numpy.pad` function.
+
+    Note that the behavior of the ``'constant'`` mode here follows the
+    PyWavelets convention which is different from NumPy (it is equivalent to
+    ``mode='edge'`` in :func:`numpy.pad`).
+    """
+    x = np.asanyarray(x)
+
+    # process pad_widths exactly as in numpy.pad
+    pad_widths = np.array(pad_widths)
+    pad_widths = np.round(pad_widths).astype(np.intp, copy=False)
+    if pad_widths.min() < 0:
+        raise ValueError("pad_widths must be > 0")
+    pad_widths = np.broadcast_to(pad_widths, (x.ndim, 2)).tolist()
+
+    if mode in ['symmetric', 'reflect']:
+        xp = np.pad(x, pad_widths, mode=mode)
+    elif mode in ['periodic', 'periodization']:
+        if mode == 'periodization':
+            # Promote odd-sized dimensions to even length by duplicating the
+            # last value.
+            edge_pad_widths = [(0, x.shape[ax] % 2)
+                               for ax in range(x.ndim)]
+            x = np.pad(x, edge_pad_widths, mode='edge')
+        xp = np.pad(x, pad_widths, mode='wrap')
+    elif mode == 'zero':
+        xp = np.pad(x, pad_widths, mode='constant', constant_values=0)
+    elif mode == 'constant':
+        xp = np.pad(x, pad_widths, mode='edge')
+    elif mode == 'smooth':
+        def pad_smooth(vector, pad_width, iaxis, kwargs):
+            # smooth extension to left
+            left = vector[pad_width[0]]
+            slope_left = (left - vector[pad_width[0] + 1])
+            vector[:pad_width[0]] = \
+                left + np.arange(pad_width[0], 0, -1) * slope_left
+
+            # smooth extension to right
+            right = vector[-pad_width[1] - 1]
+            slope_right = (right - vector[-pad_width[1] - 2])
+            vector[-pad_width[1]:] = \
+                right + np.arange(1, pad_width[1] + 1) * slope_right
+            return vector
+        xp = np.pad(x, pad_widths, pad_smooth)
+    elif mode == 'antisymmetric':
+        def pad_antisymmetric(vector, pad_width, iaxis, kwargs):
+            # smooth extension to left
+            # implement by flipping portions symmetric padding
+            npad_l, npad_r = pad_width
+            vsize_nonpad = vector.size - npad_l - npad_r
+            # Note: must modify vector in-place
+            vector[:] = np.pad(vector[pad_width[0]:-pad_width[-1]],
+                               pad_width, mode='symmetric')
+            vp = vector
+            r_edge = npad_l + vsize_nonpad - 1
+            l_edge = npad_l
+            # width of each reflected segment
+            seg_width = vsize_nonpad
+            # flip reflected segments on the right of the original signal
+            n = 1
+            while r_edge <= vp.size:
+                segment_slice = slice(r_edge + 1,
+                                      min(r_edge + 1 + seg_width, vp.size))
+                if n % 2:
+                    vp[segment_slice] *= -1
+                r_edge += seg_width
+                n += 1
+
+            # flip reflected segments on the left of the original signal
+            n = 1
+            while l_edge >= 0:
+                segment_slice = slice(max(0, l_edge - seg_width), l_edge)
+                if n % 2:
+                    vp[segment_slice] *= -1
+                l_edge -= seg_width
+                n += 1
+            return vector
+        xp = np.pad(x, pad_widths, pad_antisymmetric)
+    elif mode == 'antireflect':
+        xp = np.pad(x, pad_widths, mode='reflect', reflect_type='odd')
+    else:
+        raise ValueError(
+            ("unsupported mode: {}. The supported modes are {}").format(
+                mode, Modes.modes))
+    return xp

pywt/_utils.py

@@ -2,6 +2,7 @@
 #                    <https://github.com/PyWavelets/pywt>
 # See COPYING for license details.
 import inspect
+import numpy as np
 import sys
 from collections.abc import Iterable
 
@@ -17,7 +18,7 @@
 
 
 def _as_wavelet(wavelet):
-    """Convert wavelet name to a Wavelet object"""
+    """Convert wavelet name to a Wavelet object."""
     if not isinstance(wavelet, (ContinuousWavelet, Wavelet)):
         wavelet = DiscreteContinuousWavelet(wavelet)
     if isinstance(wavelet, ContinuousWavelet):

pywt/tests/test_dwt_idwt.py

@@ -2,8 +2,8 @@
 from __future__ import division, print_function, absolute_import
 
 import numpy as np
-from numpy.testing import assert_allclose, assert_, assert_raises
-
+from numpy.testing import (assert_allclose, assert_, assert_raises,
+                           assert_array_equal)
 import pywt
 
 # Check that float32, float64, complex64, complex128 are preserved.
@@ -228,8 +228,72 @@ def test_error_on_continuous_wavelet():
 def test_dwt_zero_size_axes():
     # raise on empty input array
     assert_raises(ValueError, pywt.dwt, [], 'db2')
-  
+
     # >1D case uses a different code path so check there as well
     x = np.ones((1, 4))[0:0, :]  # 2D with a size zero axis
     assert_raises(ValueError, pywt.dwt, x, 'db2', axis=0)
 
+
+def test_pad_1d():
+    x = [1, 2, 3]
+    assert_array_equal(pywt.pad(x, (4, 6), 'periodization'),
+                       [1, 2, 3, 3, 1, 2, 3, 3, 1, 2, 3, 3, 1, 2])
+    assert_array_equal(pywt.pad(x, (4, 6), 'periodic'),
+                       [3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3])
+    assert_array_equal(pywt.pad(x, (4, 6), 'constant'),
+                       [1, 1, 1, 1, 1, 2, 3, 3, 3, 3, 3, 3, 3])
+    assert_array_equal(pywt.pad(x, (4, 6), 'zero'),
+                       [0, 0, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0])
+    assert_array_equal(pywt.pad(x, (4, 6), 'smooth'),
+                       [-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+    assert_array_equal(pywt.pad(x, (4, 6), 'symmetric'),
+                       [3, 3, 2, 1, 1, 2, 3, 3, 2, 1, 1, 2, 3])
+    assert_array_equal(pywt.pad(x, (4, 6), 'antisymmetric'),
+                       [3, -3, -2, -1, 1, 2, 3, -3, -2, -1, 1, 2, 3])
+    assert_array_equal(pywt.pad(x, (4, 6), 'reflect'),
+                       [1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 1])
+    assert_array_equal(pywt.pad(x, (4, 6), 'antireflect'),
+                       [-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+
+    # equivalence of various pad_width formats
+    assert_array_equal(pywt.pad(x, 4, 'periodic'),
+                       pywt.pad(x, (4, 4), 'periodic'))
+
+    assert_array_equal(pywt.pad(x, (4, ), 'periodic'),
+                       pywt.pad(x, (4, 4), 'periodic'))
+
+    assert_array_equal(pywt.pad(x, [(4, 4)], 'periodic'),
+                       pywt.pad(x, (4, 4), 'periodic'))
+
+
+def test_pad_errors():
+    # negative pad width
+    x = [1, 2, 3]
+    assert_raises(ValueError, pywt.pad, x, -2, 'periodic')
+
+    # wrong length pad width
+    assert_raises(ValueError, pywt.pad, x, (1, 1, 1), 'periodic')
+
+    # invalid mode name
+    assert_raises(ValueError, pywt.pad, x, 2, 'bad_mode')
+
+
+def test_pad_nd():
+    for ndim in [2, 3]:
+        x = np.arange(4**ndim).reshape((4, ) * ndim)
+        if ndim == 2:
+            pad_widths = [(2, 1), (2, 3)]
+        else:
+            pad_widths = [(2, 1), ] * ndim
+        for mode in pywt.Modes.modes:
+            xp = pywt.pad(x, pad_widths, mode)
+
+            # expected result is the same as applying along axes separably
+            xp_expected = x.copy()
+            for ax in range(ndim):
+                xp_expected = np.apply_along_axis(pywt.pad,
+                                                  ax,
+                                                  xp_expected,
+                                                  pad_widths=[pad_widths[ax]],
+                                                  mode=mode)
+            assert_array_equal(xp, xp_expected)