Package mvpa :: Package datasets :: Module miscfx

[hide private]
[frames] | no frames]

# Source Code for Module mvpa.datasets.miscfx

      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.
     10
     11  All the functions defined in this module must accept dataset as the
     12  first argument since they are bound to Dataset class in the trailer.
     13  """
     14
     15  __docformat__ = 'restructuredtext'
     16
     17  from operator import isSequenceType
     18
     19  import numpy as N
     20
     21  from mvpa.datasets.base import Dataset, datasetmethod
     22  from mvpa.base.dochelpers import table2string
     23  from mvpa.misc.support import getBreakPoints
     24
     25  from mvpa.base import externals, warning
     26
     27  if __debug__:
     28      from mvpa.base import debug
     29
     30  if externals.exists('scipy'):
     31      from mvpa.datasets.miscfx_sp import detrend
     32
     33
     34  @datasetmethod
     35 -def zscore(dataset, mean=None, std=None,
     36             perchunk=True, baselinelabels=None,
     37             pervoxel=True, targetdtype='float64'):
     38      """Z-Score the samples of a `Dataset` (in-place).
     39
     40      `mean` and `std` can be used to pass custom values to the z-scoring.
     41      Both may be scalars or arrays.
     42
     43      All computations are done *in place*. Data upcasting is done
     44      automatically if necessary into `targetdtype`
     45
     46      If `baselinelabels` provided, and `mean` or `std` aren't provided, it would
     47      compute the corresponding measure based only on labels in `baselinelabels`
     48
     49      If `perchunk` is True samples within the same chunk are z-scored independent
     50      of samples from other chunks, e.i. mean and standard deviation are
     51      calculated individually.
     52      """
     53
     54      if __debug__ and perchunk \
     55        and N.array(dataset.samplesperchunk.values()).min() <= 2:
     56          warning("Z-scoring chunk-wise and one chunk with less than three "
     57                  "samples will set features in these samples to either zero "
     58                  "(with 1 sample in a chunk) "
     59                  "or -1/+1 (with 2 samples in a chunk).")
     60
     61      # cast to floating point datatype if necessary
     62      if str(dataset.samples.dtype).startswith('uint') \
     63         or str(dataset.samples.dtype).startswith('int'):
     64          dataset.setSamplesDType(targetdtype)
     65
     66      def doit(samples, mean, std, statsamples=None):
     67          """Internal method."""
     68
     69          if statsamples is None:
     70              # if nothing provided  -- mean/std on all samples
     71              statsamples = samples
     72
     73          if pervoxel:
     74              axisarg = {'axis':0}
     75          else:
     76              axisarg = {}
     77
     78          # calculate mean if necessary
     79          if mean is None:
     80              mean = statsamples.mean(**axisarg)
     81
     82          # de-mean
     83          samples -= mean
     84
     85          # calculate std-deviation if necessary
     86          # XXX YOH: would that be actually what we want?
     87          #          may be we want actually estimate of deviation from the mean,
     88          #          which per se might be not statsamples.mean (see above)?
     89          #          if logic to be changed -- adjust ZScoreMapper as well
     90          if std is None:
     91              std = statsamples.std(**axisarg)
     92
     93          # do the z-scoring
     94          if pervoxel:
     95              # Assure std being an array
     96              if N.isscalar(std):
     97                  std = N.ones(samples.shape[1])
     98              else:
     99                  # and so we don't perform list comparison to 0
    100                  std = N.asanyarray(std)
    101              samples[:, std != 0] /= std[std != 0]
    102          else:
    103              samples /= std
    104
    105          return samples
    106
    107      if baselinelabels is None:
    108          statids = None
    109      else:
    110          statids = set(dataset.idsbylabels(baselinelabels))
    111
    112      # for the sake of speed yoh didn't simply create a list
    113      # [True]*dataset.nsamples to provide easy selection of everything
    114      if perchunk:
    115          for c in dataset.uniquechunks:
    116              slicer = N.where(dataset.chunks == c)[0]
    117              if not statids is None:
    118                  statslicer = list(statids.intersection(set(slicer)))
    119                  dataset.samples[slicer] = doit(dataset.samples[slicer],
    120                                                 mean, std,
    121                                                 dataset.samples[statslicer])
    122              else:
    123                  slicedsamples = dataset.samples[slicer]
    124                  dataset.samples[slicer] = doit(slicedsamples,
    125                                                 mean, std,
    126                                                 slicedsamples)
    127      elif statids is None:
    128          doit(dataset.samples, mean, std, dataset.samples)
    129      else:
    130          doit(dataset.samples, mean, std, dataset.samples[list(statids)])
    131
    132
    133  @datasetmethod
    134 -def aggregateFeatures(dataset, fx=N.mean):
    135      """Apply a function to each row of the samples matrix of a dataset.
    136
    137      The functor given as `fx` has to honour an `axis` keyword argument in the
    138      way that NumPy used it (e.g. NumPy.mean, var).
    139
    140      :Returns:
    141         a new `Dataset` object with the aggregated feature(s).
    142      """
    143      agg = fx(dataset.samples, axis=1)
    144
    145      return Dataset(samples=N.array(agg, ndmin=2).T,
    146                     labels=dataset.labels,
    147                     chunks=dataset.chunks)
    148
    149
    150  @datasetmethod
    151 -def removeInvariantFeatures(dataset):
    152      """Returns a new dataset with all invariant features removed.
    153      """
    154      return dataset.selectFeatures(dataset.samples.std(axis=0).nonzero()[0])
    155
    156
    157  @datasetmethod
    158 -def coarsenChunks(source, nchunks=4):
    159      """Change chunking of the dataset
    160
    161      Group chunks into groups to match desired number of chunks. Makes
    162      sense if originally there were no strong groupping into chunks or
    163      each sample was independent, thus belonged to its own chunk
    164
    165      :Parameters:
    166        source : Dataset or list of chunk ids
    167          dataset or list of chunk ids to operate on. If Dataset, then its chunks
    168          get modified
    169        nchunks : int
    170          desired number of chunks
    171      """
    172
    173      if isinstance(source, Dataset):
    174          chunks = source.chunks
    175      else:
    176          chunks = source
    177      chunks_unique = N.unique(chunks)
    178      nchunks_orig = len(chunks_unique)
    179
    180      if nchunks_orig < nchunks:
    181          raise ValueError, \
    182                "Original number of chunks is %d. Cannot coarse them " \
    183                "to get %d chunks" % (nchunks_orig, nchunks)
    184
    185      # figure out number of samples per each chunk
    186      counts = dict(zip(chunks_unique, [ 0 ] * len(chunks_unique)))
    187      for c in chunks:
    188          counts[c] += 1
    189
    190      # now we need to group chunks to get more or less equalized number
    191      # of samples per chunk. No sophistication is done -- just
    192      # consecutively group to get close to desired number of samples
    193      # per chunk
    194      avg_chunk_size = N.sum(counts.values())*1.0/nchunks
    195      chunks_groups = []
    196      cur_chunk = []
    197      nchunks = 0
    198      cur_chunk_nsamples = 0
    199      samples_counted = 0
    200      for i, c in enumerate(chunks_unique):
    201          cc = counts[c]
    202
    203          cur_chunk += [c]
    204          cur_chunk_nsamples += cc
    205
    206          # time to get a new chunk?
    207          if (samples_counted + cur_chunk_nsamples
    208              >= (nchunks+1)*avg_chunk_size) or i==nchunks_orig-1:
    209              chunks_groups.append(cur_chunk)
    210              samples_counted += cur_chunk_nsamples
    211              cur_chunk_nsamples = 0
    212              cur_chunk = []
    213              nchunks += 1
    214
    215      if len(chunks_groups) != nchunks:
    216          warning("Apparently logic in coarseChunks is wrong. "
    217                  "It was desired to get %d chunks, got %d"
    218                  % (nchunks, len(chunks_groups)))
    219
    220      # remap using groups
    221      # create dictionary
    222      chunks_map = {}
    223      for i, group in enumerate(chunks_groups):
    224          for c in group:
    225              chunks_map[c] = i
    226
    227      chunks_new = [chunks_map[x] for x in chunks]
    228
    229      if __debug__:
    230          debug("DS_", "Using dictionary %s to remap old chunks %s into new %s"
    231                % (chunks_map, chunks, chunks_new))
    232
    233      if isinstance(source, Dataset):
    234          if __debug__:
    235              debug("DS", "Coarsing %d chunks into %d chunks for %s"
    236                    %(nchunks_orig, len(chunks_new), source))
    237          source.chunks = chunks_new
    238          return
    239      else:
    240          return chunks_new
    241
    242
    243  @datasetmethod
    244 -def getSamplesPerChunkLabel(dataset):
    245      """Returns an array with the number of samples per label in each chunk.
    246
    247      Array shape is (chunks x labels).
    248
    249      :Parameters:
    250        dataset: Dataset
    251          Source dataset.
    252      """
    253      ul = dataset.uniquelabels
    254      uc = dataset.uniquechunks
    255
    256      count = N.zeros((len(uc), len(ul)), dtype='uint')
    257
    258      for cc, c in enumerate(uc):
    259          for lc, l in enumerate(ul):
    260              count[cc, lc] = N.sum(N.logical_and(dataset.labels == l,
    261                                                  dataset.chunks == c))
    262
    263      return count
    264
    265
    266 -class SequenceStats(dict):
    267      """Simple helper to provide representation of sequence statistics
    268
    269      Matlab analog:
    270      http://cfn.upenn.edu/aguirre/code/matlablib/mseq/mtest.m
    271
    272      WARNING: Experimental -- API might change without warning!
    273      Current implementation is ugly!
    274      """
    275
    276 -    def __init__(self, seq, order=2):#, chunks=None, perchunk=False):
    277          """Initialize SequenceStats
    278
    279          :Parameters:
    280            seq : list or ndarray
    281              Actual sequence of labels
    282
    283          :Keywords:
    284            order : int
    285              Maximal order of counter-balancing check. For perfect
    286              counterbalancing all matrices should be identical
    287          """
    288          """
    289            chunks : None or list or ndarray
    290              Chunks to use if `perchunk`=True
    291            perchunk .... TODO
    292            """
    293          dict.__init__(self)
    294          self.order = order
    295          self._seq = seq
    296          self.stats = None
    297          self._str_stats = None
    298          self.__compute()
    299
    300
    301 -    def __repr__(self):
    302          """Representation of SequenceStats
    303          """
    304          return "SequenceStats(%s, order=%d)" % (repr(self._seq), self.order)
    305
    306 -    def __str__(self):
    307          return self._str_stats
    308
    309 -    def __compute(self):
    310          """Compute stats and string representation
    311          """
    312          # Do actual computation
    313          order = self.order
    314          seq = list(self._seq)               # assure list
    315          nsamples = len(seq)                 # # of samples/labels
    316          ulabels = sorted(list(set(seq)))    # unique labels
    317          nlabels = len(ulabels)              # # of labels
    318
    319          # mapping for labels
    320          labels_map = dict([(l, i) for i,l in enumerate(ulabels)])
    321
    322          # map sequence first
    323          seqm = [labels_map[i] for i in seq]
    324          nperlabel = N.bincount(seqm)
    325
    326          res = dict(ulabels=ulabels)
    327          # Estimate counter-balance
    328          cbcounts = N.zeros((order, nlabels, nlabels), dtype=int)
    329          for cb in xrange(order):
    330              for i,j in zip(seqm[:-(cb+1)], seqm[cb+1:]):
    331                  cbcounts[cb, i, j] += 1
    332          res['cbcounts'] = cbcounts
    333
    334          """
    335          Lets compute relative counter-balancing
    336          Ideally, nperlabel[i]/nlabels should precede each label
    337          """
    338          # Autocorrelation
    339          corr = []
    340          # for all possible shifts:
    341          for shift in xrange(1, nsamples):
    342              shifted = seqm[shift:] + seqm[:shift]
    343              # ??? User pearsonsr with p may be?
    344              corr += [N.corrcoef(seqm, shifted)[0, 1]]
    345              # ??? report high (anti)correlations?
    346          res['corrcoef'] = corr = N.array(corr)
    347          res['sumabscorr'] = sumabscorr = N.sum(N.abs(corr))
    348          self.update(res)
    349
    350          # Assign textual summary
    351          # XXX move into a helper function and do on demand
    352          t = [ [""] * (1 + self.order*(nlabels+1)) for i in xrange(nlabels+1) ]
    353          t[0][0] = "Labels/Order"
    354          for i, l  in enumerate(ulabels):
    355              t[i+1][0] = '%s:' % l
    356          for cb in xrange(order):
    357              t[0][1+cb*(nlabels+1)] = "O%d" % (cb+1)
    358              for i  in xrange(nlabels+1):
    359                  t[i][(cb+1)*(nlabels+1)] = " | "
    360              m = cbcounts[cb]
    361              # ??? there should be better way to get indexes
    362              ind = N.where(~N.isnan(m))
    363              for i, j in zip(*ind):
    364                  t[1+i][1+cb*(nlabels+1)+j] = '%d' % m[i, j]
    365
    366          sout = "Original sequence had %d entries from set %s\n" \
    367                 % (len(seq), ulabels) + \
    368                 "Counter-balance table for orders up to %d:\n" % order \
    369                 + table2string(t)
    370          sout += "Correlations: min=%.2g max=%.2g mean=%.2g sum(abs)=%.2g" \
    371                  % (min(corr), max(corr), N.mean(corr), sumabscorr)
    372          self._str_stats = sout
    373
    374
    375 -    def plot(self):
    376          """Plot correlation coefficients
    377          """
    378          externals.exists('pylab', raiseException=True)
    379          import pylab as P
    380          P.plot(self['corrcoef'])
    381          P.title('Auto-correlation of the sequence')
    382          P.xlabel('Offset')
    383          P.ylabel('Correlation Coefficient')
    384          P.show()
    385