added wavedecn and waverecn to multilevel.py

pywt/__init__.py

@@ -19,10 +19,10 @@
 from .thresholding import *
 from .wavelet_packets import *
 
-
-__all__ = [s for s in dir() if not s.startswith('_')]
-
 from pywt.version import version as __version__
 
 from numpy.testing import Tester
+
+__all__ = [s for s in dir() if not s.startswith('_')]
+
 test = Tester().test

pywt/multidim.py

@@ -9,14 +9,15 @@
 
 from __future__ import division, print_function, absolute_import
 
-__all__ = ['dwt2', 'idwt2', 'swt2', 'dwtn', 'idwtn']
-
 from itertools import cycle, product, repeat, islice
 
 import numpy as np
+import warnings
 
 from ._pywt import Wavelet, MODES
-from ._pywt import dwt, idwt, swt, downcoef, upcoef
+from ._pywt import (dwt, idwt, swt, downcoef, upcoef)
+
+__all__ = ['dwt2', 'idwt2', 'swt2', 'dwtn', 'idwtn']
 
 
 def dwt2(data, wavelet, mode='sym'):
@@ -198,6 +199,7 @@ def idwt2(coeffs, wavelet, mode='sym'):
 def dwtn(data, wavelet, mode='sym'):
     """
     Single-level n-dimensional Discrete Wavelet Transform.
+    This version does the looping in C if the array is <= 4D.
 
     Parameters
     ----------
@@ -226,26 +228,24 @@ def dwtn(data, wavelet, mode='sym'):
 
     """
     data = np.asarray(data)
-    dim = len(data.shape)
+    dim = data.ndim
     if dim < 1:
         raise ValueError("Input data must be at least 1D")
-    coeffs = [('', data)]
 
-    def _downcoef(data, wavelet, mode, type):
-        """Adapts pywt.downcoef call for apply_along_axis"""
-        return downcoef(type, data, wavelet, mode, level=1)
+    # have to upgrade integer types to float
+    if data.dtype != np.float32:
+        data = np.asarray(data, dtype=np.float64)
 
+    coeffs = [('', data)]
     for axis in range(dim):
         new_coeffs = []
         for subband, x in coeffs:
-            new_coeffs.extend([
-                (subband + 'a', np.apply_along_axis(_downcoef, axis, x,
-                                                    wavelet, mode, 'a')),
-                (subband + 'd', np.apply_along_axis(_downcoef, axis, x,
-                                                    wavelet, mode, 'd'))])
-
+            a = downcoef('a', x, wavelet=wavelet, mode=mode, level=1,
+                         axis=axis)
+            d = downcoef('d', x, wavelet=wavelet, mode=mode, level=1,
+                         axis=axis)
+            new_coeffs.extend([(subband + 'a', a), (subband + 'd', d)])
         coeffs = new_coeffs
-
     return dict(coeffs)
 
 
@@ -286,11 +286,12 @@ def idwtn(coeffs, wavelet, mode='sym', take=None):
     dims = max(len(key) for key in coeffs.keys())
 
     try:
-        coeff_shapes = ( v.shape for k, v in coeffs.items()
-                         if v is not None and len(k) == dims )
+        coeff_shapes = (v.shape for k, v in coeffs.items()
+                        if v is not None and len(k) == dims)
         coeff_shape = next(coeff_shapes)
     except StopIteration:
-        raise ValueError("`coeffs` must contain at least one non-null wavelet band")
+        raise ValueError("`coeffs` must contain at least one non-null wavelet "
+                         "band")
     if any(s != coeff_shape for s in coeff_shapes):
         raise ValueError("`coeffs` must all be of equal size (or None)")
 
@@ -307,24 +308,20 @@ def idwtn(coeffs, wavelet, mode='sym', take=None):
         else:
             takes = [2*s - wavelet.rec_len + 2 for s in reversed(coeff_shape)]
 
-    def _upcoef(coeffs, wavelet, take, type):
-        """Adapts pywt.upcoef call for apply_along_axis"""
-        return upcoef(type, coeffs, wavelet, level=1, take=take)
-
     for axis, take in zip(reversed(range(dims)), takes):
         new_coeffs = {}
-        new_keys = [''.join(coeff) for coeff in product('ad', repeat=axis)]
-
+        new_keys = \
+            [''.join(coeff) for coeff in product('ad', repeat=axis)]
         for key in new_keys:
             L = coeffs.get(key + 'a')
             H = coeffs.get(key + 'd')
-
+            # add new axes up to 4D for compatibility with the cython code
             if L is not None:
-                L = np.apply_along_axis(_upcoef, axis, L, wavelet, take, 'a')
-
+                L = upcoef('a', L, wavelet=wavelet, level=1, take=take,
+                           axis=axis)
             if H is not None:
-                H = np.apply_along_axis(_upcoef, axis, H, wavelet, take, 'd')
-
+                H = upcoef('d', H, wavelet=wavelet, level=1, take=take,
+                           axis=axis)
             if H is None and L is None:
                 new_coeffs[key] = None
             elif H is None:
@@ -333,7 +330,6 @@ def _upcoef(coeffs, wavelet, take, type):
                 new_coeffs[key] = H
             else:
                 new_coeffs[key] = L + H
-
         coeffs = new_coeffs
 
     return coeffs['']
@@ -414,7 +410,8 @@ def swt2(data, wavelet, level, start_level=0):
 
         # build result structure: (approx, (horizontal, vertical, diagonal))
         approx = np.transpose(LL)
-        ret.append((approx, (np.transpose(LH), np.transpose(HL), np.transpose(HH))))
+        ret.append((approx,
+                    (np.transpose(LH), np.transpose(HL), np.transpose(HH))))
 
         # for next iteration
         data = approx

pywt/multilevel.py

@@ -10,13 +10,14 @@
 
 from __future__ import division, print_function, absolute_import
 
-__all__ = ['wavedec', 'waverec', 'wavedec2', 'waverec2']
-
 import numpy as np
 
 from ._pywt import Wavelet
 from ._pywt import dwt, idwt, dwt_max_level
-from .multidim import dwt2, idwt2
+from .multidim import dwt2, idwt2, dwtn, idwtn
+
+__all__ = ['wavedec', 'waverec', 'wavedec2', 'waverec2', 'wavedecn',
+           'waverecn']
 
 
 def wavedec(data, wavelet, mode='sym', level=None):
@@ -39,9 +40,9 @@ def wavedec(data, wavelet, mode='sym', level=None):
     -------
     [cA_n, cD_n, cD_n-1, ..., cD2, cD1] : list
         Ordered list of coefficients arrays
-        where `n` denotes the level of decomposition. The first element (`cA_n`) of
-        the result is approximation coefficients array and the following elements
-        (`cD_n` - `cD_1`) are details coefficients arrays.
+        where `n` denotes the level of decomposition. The first element
+        (`cA_n`) of the result is approximation coefficients array and the
+        following elements (`cD_n` - `cD_1`) are details coefficients arrays.
 
     Examples
     --------
@@ -220,3 +221,154 @@ def waverec2(coeffs, wavelet, mode='sym'):
         a = idwt2((a, d), wavelet, mode)
 
     return a
+
+
+def wavedecn(data, wavelet, mode='sym', level=None):
+    """
+    Multilevel nD Discrete Wavelet Transform.
+
+    Parameters
+    ----------
+    data : ndarray
+        nD input data
+    wavelet : Wavelet object or name string
+        Wavelet to use
+    mode : str, optional
+        Signal extension mode, see MODES (default: 'sym')
+    level : int, optional
+        Decomposition level. If level is None (default) then it will be
+        calculated using `dwt_max_level` function.
+
+    Returns
+    -------
+    [cAn, {diff_coeffs_leveln}, ... {diff_coeffs_level1}] : list
+        Coefficients list
+
+    Examples
+    --------
+    >>> from pywt import multilevel
+    >>> coeffs = multilevel.wavedecn(np.ones((4,4,4)), 'db1')
+    >>> # Levels:
+    >>> len(coeffs)-1
+    3
+    >>> multilevel.waverecn(coeffs, 'db1')
+    array([[[ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.]],
+
+       [[ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.]],
+
+       [[ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.]],
+
+       [[ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.]]])
+
+    """
+
+    data = np.asarray(data, np.float64)
+
+    if len(data.shape) < 1:
+        raise ValueError("Expected at least 1D input data.")
+
+    if not isinstance(wavelet, Wavelet):
+        wavelet = Wavelet(wavelet)
+
+    if level is None:
+        size = min(data.shape)
+        level = dwt_max_level(size, wavelet.dec_len)
+    elif level < 0:
+        raise ValueError(
+            "Level value of %d is too low . Minimum level is 0." % level)
+
+    coeffs_list = []
+
+    a = data
+    for i in range(level):
+        coeffs = dwtn(a, wavelet, mode)
+        a = coeffs.pop('a'*data.ndim)
+        coeffs_list.append(coeffs)
+
+    coeffs_list.append(a)
+    coeffs_list.reverse()
+
+    return coeffs_list
+
+
+def waverecn(coeffs, wavelet, mode='sym'):
+    """
+    Multilevel nD Inverse Discrete Wavelet Transform.
+
+    coeffs : array_like
+        Coefficients list [cAn, {diff_coeffs_leveln}, ... {diff_coeffs_level1}]
+    wavelet : Wavelet object or name string
+        Wavelet to use
+    mode : str, optional
+        Signal extension mode, see MODES (default: 'sym')
+
+    Returns
+    -------
+    nD array of reconstructed data.
+
+    Examples
+    --------
+    >>> from pywt import multilevel
+    >>> coeffs = multilevel.wavedecn(np.ones((4,4,4)), 'db1')
+    >>> # Levels:
+    >>> len(coeffs)-1
+    2
+    >>> multilevel.waverecn(coeffs, 'db1')
+    array([[[ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.]],
+
+       [[ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.]],
+
+       [[ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.]],
+
+       [[ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.],
+        [ 1.,  1.,  1.,  1.]]])
+    """
+
+    if not isinstance(coeffs, (list, tuple)):
+        raise ValueError("Expected sequence of coefficient arrays.")
+
+    if len(coeffs) < 2:
+        raise ValueError(
+            "Coefficient list too short (minimum 2 arrays required).")
+
+    a, ds = coeffs[0], coeffs[1:]
+
+    dims = max(len(key) for key in ds[0].keys())
+
+    for idx, d in enumerate(ds):
+
+        # determine the shape at the next level for idwtn
+        if idx < (len(ds)-1):
+            next_shape = [
+                v.shape for k, v in ds[idx+1].items() if v is not None]
+            take = next_shape[0]
+        else:
+            take = None
+
+        d['a'*dims] = a
+        a = idwtn(d, wavelet, mode, take=take)
+
+    return a

pywt/setup.py

@@ -2,7 +2,7 @@
 from __future__ import division, print_function, absolute_import
 
 
-def configuration(parent_package='',top_path=None):
+def configuration(parent_package='', top_path=None):
     from numpy.distutils.misc_util import Configuration
     import numpy as np