Package mvpa :: Package mappers :: Module wavelet

[hide private]
[frames] | no frames]

# Source Code for Module mvpa.mappers.wavelet

      1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
      2  # vi: set ft=python sts=4 ts=4 sw=4 et:
      3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
      4  #
      5  #   See COPYING file distributed along with the PyMVPA package for the
      6  #   copyright and license terms.
      7  #
      8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
      9  """Wavelet mappers"""
     10
     11  from mvpa.base import externals
     12
     13  if externals.exists('pywt', raiseException=True):
     14      # import conditional to be able to import the whole module while building
     15      # the docs even if pywt is not installed
     16      import pywt
     17
     18  import numpy as N
     19
     20  from mvpa.base import warning
     21  from mvpa.mappers.base import Mapper
     22  from mvpa.base.dochelpers import enhancedDocString
     23
     24  if __debug__:
     25      from mvpa.base import debug
     26
     27  # WaveletPacket and WaveletTransformation mappers share lots of common
     28  # functionality at the moment
     29
     30 -class _WaveletMapper(Mapper):
     31      """Generic class for Wavelet mappers (decomposition and packet)
     32      """
     33
     34 -    def __init__(self, dim=1, wavelet='sym4', mode='per', maxlevel=None):
     35          """Initialize _WaveletMapper mapper
     36
     37          :Parameters:
     38            dim : int or tuple of int
     39              dimensions to work across (for now just scalar value, ie 1D
     40              transformation) is supported
     41            wavelet : basestring
     42              one from the families available withing pywt package
     43            mode : basestring
     44              periodization mode
     45            maxlevel : int or None
     46              number of levels to use. If None - automatically selected by pywt
     47          """
     48          Mapper.__init__(self)
     49
     50          self._dim = dim
     51          """Dimension to work along"""
     52
     53          self._maxlevel = maxlevel
     54          """Maximal level of decomposition. None for automatic"""
     55
     56          if not wavelet in pywt.wavelist():
     57              raise ValueError, \
     58                    "Unknown family of wavelets '%s'. Please use one " \
     59                    "available from the list %s" % (wavelet, pywt.wavelist())
     60          self._wavelet = wavelet
     61          """Wavelet family to use"""
     62
     63          if not mode in pywt.MODES.modes:
     64              raise ValueError, \
     65                    "Unknown periodization mode '%s'. Please use one " \
     66                    "available from the list %s" % (mode, pywt.MODES.modes)
     67          self._mode = mode
     68          """Periodization mode"""
     69
     70
     71 -    def forward(self, data):
     72          data = N.asanyarray(data)
     73          self._inshape = data.shape
     74          self._intimepoints = data.shape[self._dim]
     75          res = self._forward(data)
     76          self._outshape = res.shape
     77          return res
     78
     79
     80 -    def reverse(self, data):
     81          data = N.asanyarray(data)
     82          return self._reverse(data)
     83
     84
     85 -    def _forward(self, *args):
     86          raise NotImplementedError
     87
     88
     89 -    def _reverse(self, *args):
     90          raise NotImplementedError
     91
     92
     93 -    def getInSize(self):
     94          """Returns the number of original features."""
     95          return self._inshape[1:]
     96
     97
     98 -    def getOutSize(self):
     99          """Returns the number of wavelet components."""
    100          return self._outshape[1:]
    101
    102
    103 -    def selectOut(self, outIds):
    104          """Choose a subset of components...
    105
    106          just use MaskMapper on top?"""
    107          raise NotImplementedError, "Please use in conjunction with MaskMapper"
    108
    109
    110      __doc__ = enhancedDocString('_WaveletMapper', locals(), Mapper)
    111
    112
    113 -def _getIndexes(shape, dim):
    114      """Generator for coordinate tuples providing slice for all in `dim`
    115
    116      XXX Somewhat sloppy implementation... but works...
    117      """
    118      if len(shape) < dim:
    119          raise ValueError, "Dimension %d is incorrect for a shape %s" % \
    120                (dim, shape)
    121      n = len(shape)
    122      curindexes = [0] * n
    123      curindexes[dim] = Ellipsis#slice(None)       # all elements for dimension dim
    124      while True:
    125          yield tuple(curindexes)
    126          for i in xrange(n):
    127              if i == dim and dim == n-1:
    128                  return                  # we reached it -- thus time to go
    129              if curindexes[i] == shape[i] - 1:
    130                  if i == n-1:
    131                      return
    132                  curindexes[i] = 0
    133              else:
    134                  if i != dim:
    135                      curindexes[i] += 1
    136                      break
    137
    138
    139 -class WaveletPacketMapper(_WaveletMapper):
    140      """Convert signal into an overcomplete representaion using Wavelet packet
    141      """
    142
    143 -    def __init__(self, level=None, **kwargs):
    144          """Initialize WaveletPacketMapper mapper
    145
    146          :Parameters:
    147            level : int or None
    148              What level to decompose at. If 'None' data for all levels
    149              is provided, but due to different sizes, they are placed
    150              in 1D row.
    151          """
    152
    153          _WaveletMapper.__init__(self,**kwargs)
    154
    155          self.__level = level
    156
    157
    158      # XXX too much of duplications between such methods -- it begs
    159      #     refactoring
    160 -    def __forwardSingleLevel(self, data):
    161          if __debug__:
    162              debug('MAP', "Converting signal using DWP (single level)")
    163
    164          wp = None
    165
    166          level = self.__level
    167          wavelet = self._wavelet
    168          mode = self._mode
    169          dim = self._dim
    170
    171          level_paths = None
    172          for indexes in _getIndexes(data.shape, self._dim):
    173              if __debug__:
    174                  debug('MAP_', " %s" % (indexes,), lf=False, cr=True)
    175              WP = pywt.WaveletPacket(
    176                  data[indexes], wavelet=wavelet,
    177                  mode=mode, maxlevel=level)
    178
    179              level_nodes = WP.get_level(level)
    180              if level_paths is None:
    181                  # Needed for reconstruction
    182                  self.__level_paths = N.array([node.path for node in level_nodes])
    183              level_datas = N.array([node.data for node in level_nodes])
    184
    185              if wp is None:
    186                  newdim = data.shape
    187                  newdim = newdim[:dim] + level_datas.shape + newdim[dim+1:]
    188                  if __debug__:
    189                      debug('MAP_', "Initializing storage of size %s for single "
    190                            "level (%d) mapping of data of size %s" % (newdim, level, data.shape))
    191                  wp = N.empty( tuple(newdim) )
    192
    193              wp[indexes] = level_datas
    194
    195          return wp
    196
    197
    198 -    def __forwardMultipleLevels(self, data):
    199          wp = None
    200          levels_length = None                # total length at each level
    201          levels_lengths = None                # list of lengths per each level
    202          for indexes in _getIndexes(data.shape, self._dim):
    203              if __debug__:
    204                  debug('MAP_', " %s" % (indexes,), lf=False, cr=True)
    205              WP = pywt.WaveletPacket(
    206                  data[indexes],
    207                  wavelet=self._wavelet,
    208                  mode=self._mode, maxlevel=self._maxlevel)
    209
    210              if levels_length is None:
    211                  levels_length = [None] * WP.maxlevel
    212                  levels_lengths = [None] * WP.maxlevel
    213
    214              levels_datas = []
    215              for level in xrange(WP.maxlevel):
    216                  level_nodes = WP.get_level(level+1)
    217                  level_datas = [node.data for node in level_nodes]
    218
    219                  level_lengths = [len(x) for x in level_datas]
    220                  level_length = N.sum(level_lengths)
    221
    222                  if levels_lengths[level] is None:
    223                      levels_lengths[level] = level_lengths
    224                  elif levels_lengths[level] != level_lengths:
    225                      raise RuntimeError, \
    226                            "ADs of same level of different samples should have same number of elements." \
    227                            " Got %s, was %s" % (level_lengths, levels_lengths[level])
    228
    229                  if levels_length[level] is None:
    230                      levels_length[level] = level_length
    231                  elif levels_length[level] != level_length:
    232                      raise RuntimeError, \
    233                            "Levels of different samples should have same number of elements." \
    234                            " Got %d, was %d" % (level_length, levels_length[level])
    235
    236                  level_data = N.hstack(level_datas)
    237                  levels_datas.append(level_data)
    238
    239              # assert(len(data) == levels_length)
    240              # assert(len(data) >= Ntimepoints)
    241              if wp is None:
    242                  newdim = list(data.shape)
    243                  newdim[self._dim] = N.sum(levels_length)
    244                  wp = N.empty( tuple(newdim) )
    245              wp[indexes] = N.hstack(levels_datas)
    246
    247          self.levels_lengths, self.levels_length = levels_lengths, levels_length
    248          if __debug__:
    249              debug('MAP_', "")
    250              debug('MAP', "Done convertion into wp. Total size %s" % str(wp.shape))
    251          return wp
    252
    253
    254 -    def _forward(self, data):
    255          if __debug__:
    256              debug('MAP', "Converting signal using DWP")
    257
    258          if self.__level is None:
    259              return self.__forwardMultipleLevels(data)
    260          else:
    261              return self.__forwardSingleLevel(data)
    262
    263      #
    264      # Reverse mapping
    265      #
    266 -    def __reverseSingleLevel(self, wp):
    267
    268          # local bindings
    269          level_paths = self.__level_paths
    270
    271          # define wavelet packet to use
    272          WP = pywt.WaveletPacket(
    273              data=None, wavelet=self._wavelet,
    274              mode=self._mode, maxlevel=self.__level)
    275
    276          # prepare storage
    277          signal_shape = wp.shape[:1] + self.getInSize()
    278          signal = N.zeros(signal_shape)
    279          Ntime_points = self._intimepoints
    280          for indexes in _getIndexes(signal_shape,
    281                                     self._dim):
    282              if __debug__:
    283                  debug('MAP_', " %s" % (indexes,), lf=False, cr=True)
    284
    285              for path, level_data in zip(level_paths, wp[indexes]):
    286                  WP[path] = level_data
    287
    288              signal[indexes] = WP.reconstruct(True)[:Ntime_points]
    289
    290          return signal
    291
    292
    293 -    def _reverse(self, data):
    294          if __debug__:
    295              debug('MAP', "Converting signal back using DWP")
    296
    297          if self.__level is None:
    298              raise NotImplementedError
    299          else:
    300              if not externals.exists('pywt wp reconstruct'):
    301                  raise NotImplementedError, \
    302                        "Reconstruction for a single level for versions of " \
    303                        "pywt < 0.1.7 (revision 103) is not supported"
    304              if not externals.exists('pywt wp reconstruct fixed'):
    305                  warning("Reconstruction using available version of pywt might "
    306                          "result in incorrect data in the tails of the signal")
    307              return self.__reverseSingleLevel(data)
    308
    309
    310
    311
    312
    313 -class WaveletTransformationMapper(_WaveletMapper):
    314      """Convert signal into wavelet representaion
    315      """
    316
    317 -    def _forward(self, data):
    318          """Decompose signal into wavelets's coefficients via dwt
    319          """
    320          if __debug__:
    321              debug('MAP', "Converting signal using DWT")
    322          wd = None
    323          coeff_lengths = None
    324          for indexes in _getIndexes(data.shape, self._dim):
    325              if __debug__:
    326                  debug('MAP_', " %s" % (indexes,), lf=False, cr=True)
    327              coeffs = pywt.wavedec(
    328                  data[indexes],
    329                  wavelet=self._wavelet,
    330                  mode=self._mode,
    331                  level=self._maxlevel)
    332              # Silly Yarik embedds extraction of statistics right in place
    333              #stats = []
    334              #for coeff in coeffs:
    335              #    stats_ = [N.std(coeff),
    336              #              N.sqrt(N.dot(coeff, coeff)),
    337              #              ]# + list(N.histogram(coeff, normed=True)[0]))
    338              #    stats__ = list(coeff) + stats_[:]
    339              #    stats__ += list(N.log(stats_))
    340              #    stats__ += list(N.sqrt(stats_))
    341              #    stats__ += list(N.array(stats_)**2)
    342              #    stats__ += [  N.median(coeff), N.mean(coeff), scipy.stats.kurtosis(coeff) ]
    343              #    stats.append(stats__)
    344              #coeffs = stats
    345              coeff_lengths_ = N.array([len(x) for x in coeffs])
    346              if coeff_lengths is None:
    347                  coeff_lengths = coeff_lengths_
    348              assert((coeff_lengths == coeff_lengths_).all())
    349              if wd is None:
    350                  newdim = list(data.shape)
    351                  newdim[self._dim] = N.sum(coeff_lengths)
    352                  wd = N.empty( tuple(newdim) )
    353              coeff = N.hstack(coeffs)
    354              wd[indexes] = coeff
    355          if __debug__:
    356              debug('MAP_', "")
    357              debug('MAP', "Done DWT. Total size %s" % str(wd.shape))
    358          self.lengths = coeff_lengths
    359          return wd
    360
    361
    362 -    def _reverse(self, wd):
    363          if __debug__:
    364              debug('MAP', "Performing iDWT")
    365          signal = None
    366          wd_offsets = [0] + list(N.cumsum(self.lengths))
    367          Nlevels = len(self.lengths)
    368          Ntime_points = self._intimepoints #len(time_points)
    369          # unfortunately sometimes due to padding iDWT would return longer
    370          # sequences, thus we just limit to the right ones
    371
    372          for indexes in _getIndexes(wd.shape, self._dim):
    373              if __debug__:
    374                  debug('MAP_', " %s" % (indexes,), lf=False, cr=True)
    375              wd_sample = wd[indexes]
    376              wd_coeffs = [wd_sample[wd_offsets[i]:wd_offsets[i+1]] for i in xrange(Nlevels)]
    377              # need to compose original list
    378              time_points = pywt.waverec(
    379                  wd_coeffs, wavelet=self._wavelet, mode=self._mode)
    380              if signal is None:
    381                  newdim = list(wd.shape)
    382                  newdim[self._dim] = Ntime_points
    383                  signal = N.empty(newdim)
    384              signal[indexes] = time_points[:Ntime_points]
    385          if __debug__:
    386              debug('MAP_', "")
    387              debug('MAP', "Done iDWT. Total size %s" % (signal.shape, ))
    388          return signal
    389