Package mvpa :: Package datasets :: Module miscfx_sp

[hide private]
[frames] | no frames]

# Source Code for Module mvpa.datasets.miscfx_sp

      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  """Misc function performing operations on datasets which are based on scipy
     10  """
     11
     12  __docformat__ = 'restructuredtext'
     13
     14  from mvpa.base import externals
     15
     16  import numpy as N
     17  np = N                      # just for easy merging of changes into >=0.5
     18
     19  from operator import isSequenceType
     20
     21  from mvpa.datasets.base import Dataset, datasetmethod
     22  from mvpa.misc.support import getBreakPoints
     23
     24  if externals.exists('scipy', raiseException=True):
     25      from scipy import signal
     26      from scipy.linalg import lstsq
     27      from scipy.special import legendre
     28
     29 -    def legendre_(n, x):
     30          # Scipy 0.8.0 (and possibly later) has regression of reporting
     31          # 'inf's for negative boundary. Lets guard against it for now
     32          leg = legendre(n)
     33          r = leg(x)
     34          infs = np.isinf(r)
     35          if np.any(infs):
     36              r[infs] = leg(x[infs] + 1e-10) # offset to try to overcome problems
     37          return r
     38
     39  @datasetmethod
     40 -def detrend(dataset, perchunk=False, model='linear',
     41              polyord=None, opt_reg=None):
     42      """
     43      Given a dataset, detrend the data inplace either entirely
     44      or per each chunk
     45
     46      :Parameters:
     47        dataset : Dataset
     48          dataset to operate on
     49        perchunk : bool
     50          either to operate on whole dataset at once or on each chunk
     51          separately
     52        model
     53          Type of detrending model to run.  If 'linear' or 'constant',
     54          scipy.signal.detrend is used to perform a linear or demeaning
     55          detrend. Polynomial detrending is activated when 'regress' is
     56          used or when polyord or opt_reg are specified.
     57        polyord : int or list
     58          Order of the Legendre polynomial to remove from the data.  This
     59          will remove every polynomial up to and including the provided
     60          value.  For example, 3 will remove 0th, 1st, 2nd, and 3rd order
     61          polynomials from the data.  N.B.: The 0th polynomial is the
     62          baseline shift, the 1st is the linear trend.
     63          If you specify a single int and perchunk is True, then this value
     64          is used for each chunk.  You can also specify a different polyord
     65          value for each chunk by providing a list or ndarray of polyord
     66          values the length of the number of chunks.
     67        opt_reg : ndarray
     68          Optional ndarray of additional information to regress out from the
     69          dataset.  One example would be to regress out motion parameters.
     70          As with the data, time is on the first axis.
     71
     72      """
     73      if polyord is not None or opt_reg is not None:
     74          model='regress'
     75
     76      if model in ['linear', 'constant']:
     77          # perform scipy detrend
     78          bp = 0                              # no break points by default
     79
     80          if perchunk:
     81              try:
     82                  bp = getBreakPoints(dataset.chunks)
     83              except ValueError, e:
     84                  raise ValueError, \
     85                        "Failed to assess break points between chunks. Often " \
     86                        "that is due to discontinuities within a chunk, which " \
     87                        "ruins idea of per-chunk detrending. Original " \
     88                        "exception was: %s" % str(e)
     89
     90          dataset.samples[:] = signal.detrend(dataset.samples, axis=0,
     91                                           type=model, bp=bp)
     92      elif model in ['regress']:
     93          # perform regression-based detrend
     94          return __detrend_regress(dataset, perchunk=perchunk,
     95                                   polyord=polyord, opt_reg=opt_reg)
     96      else:
     97          # raise exception because not found
     98          raise ValueError('Specified model type (%s) is unknown.'
     99                           % (model))
    100
    101
    102
    103 -def __detrend_regress(dataset, perchunk=True, polyord=None, opt_reg=None):
    104      """
    105      Given a dataset, perform a detrend inplace, regressing out polynomial
    106      terms as well as optional regressors, such as motion parameters.
    107
    108      :Parameters:
    109        dataset : Dataset
    110          Dataset to operate on
    111        perchunk : bool
    112          Either to operate on whole dataset at once or on each chunk
    113          separately.  If perchunk is True, all the samples within a
    114          chunk should be contiguous and the chunks should be sorted in
    115          order from low to high.
    116        polyord : int
    117          Order of the Legendre polynomial to remove from the data.  This
    118          will remove every polynomial up to and including the provided
    119          value.  For example, 3 will remove 0th, 1st, 2nd, and 3rd order
    120          polynomials from the data.  N.B.: The 0th polynomial is the
    121          baseline shift, the 1st is the linear trend.
    122          If you specify a single int and perchunk is True, then this value
    123          is used for each chunk.  You can also specify a different polyord
    124          value for each chunk by providing a list or ndarray of polyord
    125          values the length of the number of chunks.
    126        opt_reg : ndarray
    127          Optional ndarray of additional information to regress out from the
    128          dataset.  One example would be to regress out motion parameters.
    129          As with the data, time is on the first axis.
    130      """
    131
    132      # Process the polyord to be a list with length of the number of chunks
    133      if not polyord is None:
    134          if not isSequenceType(polyord):
    135              # repeat to be proper length
    136              polyord = [polyord]*len(dataset.uniquechunks)
    137          elif perchunk and len(polyord) != len(dataset.uniquechunks):
    138              raise ValueError("If you specify a sequence of polyord values " + \
    139                                   "they sequence length must match the " + \
    140                                   "number of unique chunks in the dataset.")
    141
    142      # loop over chunks if necessary
    143      if perchunk:
    144          # get the unique chunks
    145          uchunks = dataset.uniquechunks
    146
    147          # loop over each chunk
    148          reg = []
    149          for n, chunk in enumerate(uchunks):
    150              # get the indices for that chunk
    151              cinds = dataset.chunks == chunk
    152
    153              # see if add in polyord values
    154              if not polyord is None:
    155                  # create the timespan
    156                  x = N.linspace(-1, 1, cinds.sum())
    157                  # create each polyord with the value for that chunk
    158                  for n in range(polyord[n] + 1):
    159                      newreg = N.zeros((dataset.nsamples, 1))
    160                      newreg[cinds, 0] = legendre_(n, x)
    161                      reg.append(newreg)
    162      else:
    163          # take out mean over entire dataset
    164          reg = []
    165          # see if add in polyord values
    166          if not polyord is None:
    167              # create the timespan
    168              x = N.linspace(-1, 1, dataset.nsamples)
    169              for n in range(polyord[0] + 1):
    170                  reg.append(legendre_(n, x)[:, N.newaxis])
    171
    172      # see if add in optional regs
    173      if not opt_reg is None:
    174          # add in the optional regressors, too
    175          reg.append(opt_reg)
    176
    177      # combine the regs
    178      if len(reg) > 0:
    179          if len(reg) > 1:
    180              regs = N.hstack(reg)
    181          else:
    182              regs = reg[0]
    183      else:
    184          # no regs to remove
    185          raise ValueError("You must specify at least one " + \
    186                               "regressor to regress out.")
    187
    188      # perform the regression
    189      res = lstsq(regs, dataset.samples)
    190
    191      # remove all but the residuals
    192      yhat = N.dot(regs, res[0])
    193      dataset.samples -= yhat
    194
    195      # return the results
    196      return res, regs
    197