Package mvpa :: Package misc :: Package plot :: Module erp

[hide private]
[frames] | no frames]

# Source Code for Module mvpa.misc.plot.erp

      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  """Basic ERP (here ERP = Event Related Plot ;-)) plotting
     10
     11  Can be used for plotting not only ERP but any event-locked data
     12  """
     13
     14  import pylab as P
     15  import numpy as N
     16  import matplotlib as mpl
     17
     18  from mvpa.base import warning
     19  from mvpa.mappers.boxcar import BoxcarMapper
     20
     21  #
     22  # Few helper functions
     23  #
     24  import matplotlib.transforms as mlt
     25 -def _offset(ax, x, y):
     26      """Provide offset in pixels
     27
     28      :Parameters:
     29        x : int
     30          Offset in pixels for x
     31        y : int
     32          Offset in pixels for y
     33
     34      Idea borrowed from
     35       http://www.scipy.org/Cookbook/Matplotlib/Transformations
     36      but then heavily extended to be compatible with many
     37      reincarnations of matplotlib
     38      """
     39      d = dir(mlt)
     40      if 'offset_copy' in d:
     41          # tested with python-matplotlib 0.98.3-5
     42          # ??? if pukes, might need to replace 2nd parameter from
     43          #     ax to ax.get_figure()
     44          return mlt.offset_copy(ax.transData, ax, x=x, y=y, units='dots')
     45      elif 'BlendedAffine2D' in d:
     46          # some newer versions of matplotlib
     47          return ax.transData + \
     48                 mlt.Affine2D().translate(x,y)
     49      elif 'blend_xy_sep_transform' in d:
     50          trans = mlt.blend_xy_sep_transform(ax.transData, ax.transData)
     51          # Now we set the offset in pixels
     52          trans.set_offset((x, y), mlt.identity_transform())
     53          return trans
     54      else:
     55          raise RuntimeError, \
     56                "Lacking needed functions in matplotlib.transform " \
     57                "for _offset. Please upgrade"
     58
     59
     60 -def _make_centeredaxis(ax, loc, offset=5, ai=0, mult=1.0,
     61                         format='%4g', label=None, **props):
     62      """Plot an axis which is centered at loc (e.g. 0)
     63
     64      :Parameters:
     65       ax
     66         Axes from the figure
     67       loc
     68         Value to center at
     69       offset
     70         Relative offset (in pixels) for the labels
     71       ai : int
     72         Axis index: 0 for x, 1 for y
     73       mult
     74         Multiplier for the axis labels. ERPs for instance need to be
     75         inverted, thus labels too manually here since there is no easy
     76         way in matplotlib to invert an axis
     77       label : basestring or None
     78         If not -- put a label outside of the axis
     79       **props
     80         Given to underlying plotting functions
     81
     82      Idea borrowed from
     83        http://www.mail-archive.com/matplotlib-users@lists.sourceforge.net \
     84        /msg05669.html
     85      It sustained heavy refactoring/extension
     86      """
     87      xmin, xmax = ax.get_xlim()
     88      ymin, ymax = ax.get_ylim()
     89
     90      xlocs = [l for l in ax.xaxis.get_ticklocs()
     91              if l>=xmin and l<=xmax]
     92      ylocs = [l for l in ax.yaxis.get_ticklocs()
     93              if l>=ymin and l<=ymax]
     94
     95      if ai == 0:
     96          hlocs = ylocs
     97          locs = xlocs
     98          vrange = [xmin, xmax]
     99          tdir = mpl.lines.TICKDOWN
    100          halignment = 'center'
    101          valignment = 'top'
    102          lhalignment = 'left'
    103          lvalignment = 'center'
    104          lx, ly = xmax, 0
    105          ticklength = ax.xaxis.get_ticklines()[0]._markersize
    106      elif ai == 1:
    107          hlocs = xlocs
    108          locs = ylocs
    109          vrange = [ymin, ymax]
    110          tdir = mpl.lines.TICKLEFT
    111          halignment = 'right'
    112          valignment = 'center'
    113          lhalignment = 'center'
    114          lvalignment = 'bottom'
    115          lx, ly = 0, ymax
    116          ticklength = ax.yaxis.get_ticklines()[0]._markersize
    117      else:
    118          raise ValueError, "Illegal ai=%s" % ai
    119
    120      args = [ (locs, [loc]*len(locs)),
    121               (vrange, [loc, loc]),
    122               [locs, (loc,)*len(locs)]
    123               ]
    124
    125      offset_abs = offset + ticklength
    126      if ai == 1:
    127          # invert
    128          args = [ [x[1], x[0]] for x in args ]
    129          # shift the tick labels labels
    130          trans = _offset(ax, -offset_abs, 0)
    131          transl = _offset(ax, 0, offset)
    132      else:
    133          trans = _offset(ax, 0, -offset_abs)
    134          transl = _offset(ax, offset, 0)
    135
    136      tickline, = ax.plot(linestyle='', marker=tdir, *args[0], **props)
    137      axline, = ax.plot(*args[1], **props)
    138
    139      tickline.set_clip_on(False)
    140      axline.set_clip_on(False)
    141
    142
    143      for i, l in enumerate(locs):
    144          if l == 0:                    # no origin label
    145              continue
    146          coor = [args[2][0][i], args[2][1][i], format % (mult * l)]
    147          ax.text(horizontalalignment=halignment,
    148                  verticalalignment=valignment, transform=trans, *coor)
    149
    150
    151      if label is not None:
    152          ax.text(
    153              #max(args[2][0]), max(args[2][1]),
    154              lx, ly,
    155              label,
    156              horizontalalignment=lhalignment,
    157              verticalalignment=lvalignment, fontsize=14,
    158              # fontweight='bold',
    159              transform=transl)
    160
    161
    162 -def plotERP(data, SR=500, onsets=None,
    163              pre=0.2, pre_onset=None, post=None, pre_mean=None,
    164              color='r', errcolor=None, errtype=None, ax=P,
    165              ymult=1.0, *args, **kwargs):
    166      """Plot single ERP on existing canvas
    167
    168      :Parameters:
    169        data: 1D or 2D ndarray
    170          The data array can either be 1D (samples over time) or 2D
    171          (trials x samples). In the first case a boxcar mapper is used to
    172          extract the respective trial timecourses given a list of trial onsets.
    173          In the latter case, each row of the data array is taken as the EEG
    174          signal timecourse of a particular trial.
    175        onsets: list(int)
    176          List of onsets (in samples not in seconds).
    177        SR: int
    178          Sampling rate (1/s) of the signal.
    179        pre: float
    180          Duration (in seconds) to be plotted prior to onset.
    181        pre_onset : float or None
    182          If data is already in epochs (2D) then pre_onset provides information
    183          on how many seconds pre-stimulus were used to generate them. If None,
    184          then pre_onset = pre
    185        post: float
    186          Duration (in seconds) to be plotted after the onset.
    187        pre_mean: float
    188          Duration (in seconds) at the beginning of the window which is used
    189          for deriving the mean of the signal. If None, pre_mean = pre
    190        errtype: None | 'ste' | 'std' | 'ci95' | list of previous three
    191          Type of error value to be computed per datapoint.
    192            'ste': standard error of the mean
    193            'std': standard deviation
    194            'ci95': 95% confidence interval (1.96 * ste)
    195            None: no error margin is plotted (default)
    196          Optionally, multiple error types can be specified in a list. In that
    197          case all of them will be plotted.
    198        color: matplotlib color code
    199          Color to be used for plotting the mean signal timecourse.
    200        errcolor: matplotlib color code
    201          Color to be used for plotting the error margin. If None, use main color
    202          but with weak alpha level
    203        ax:
    204          Target where to draw.
    205        ymult: float
    206          Multiplier for the values. E.g. if negative-up ERP plot is needed:
    207          provide ymult=-1.0
    208        *args, **kwargs
    209          Additional arguments to plot().
    210
    211      :Returns:
    212        array
    213          Mean ERP timeseries.
    214      """
    215      if pre_mean is None:
    216          pre_mean = pre
    217
    218      # set default
    219      pre_discard = 0
    220
    221      if onsets is not None: # if we need to extract ERPs
    222          if post is None:
    223              raise ValueError, \
    224                    "Duration post onsets must be provided if onsets are given"
    225          # trial timecourse duration
    226          duration = pre + post
    227
    228          # We are working with a full timeline
    229          bcm = BoxcarMapper(onsets,
    230                             boxlength = int(SR * duration),
    231                             offset = -int(SR * pre))
    232          erp_data = bcm(data)
    233
    234          # override values since we are using Boxcar
    235          pre_onset = pre
    236      else:
    237          if pre_onset is None:
    238              pre_onset = pre
    239
    240          if pre_onset < pre:
    241              warning("Pre-stimulus interval to plot %g is smaller than provided "
    242                      "pre-stimulus captured interval %g, thus plot interval was "
    243                      "adjusted" % (pre, pre_onset))
    244              pre = pre_onset
    245
    246          if post is None:
    247              # figure out post
    248              duration = float(data.shape[1]) / SR - pre_discard
    249              post = duration - pre
    250          else:
    251              duration = pre + post
    252
    253          erp_data = data
    254          pre_discard = pre_onset - pre
    255
    256      # Scale the data appropriately
    257      erp_data *= ymult
    258
    259      # validity check -- we should have 2D matrix (trials x samples)
    260      if len(erp_data.shape) != 2:
    261          raise RuntimeError, \
    262                "plotERP() supports either 1D data with onsets, or 2D data " \
    263                "(trials x sample_points). Shape of the data at the point " \
    264                "is %s" % erp_data.shape
    265
    266      if not (pre_mean == 0 or pre_mean is None):
    267          # mean of pre-onset signal accross trials
    268          erp_baseline = N.mean(
    269              erp_data[:, int((pre_onset-pre_mean)*SR):int(pre_onset*SR)])
    270          # center data on pre-onset mean
    271          # NOTE: make sure that we make a copy of the data to don't
    272          #       alter the original. Better be safe than sorry
    273          erp_data = erp_data - erp_baseline
    274
    275      # generate timepoints and error ranges to plot filled error area
    276      # top ->
    277      # bottom <-
    278      time_points = N.arange(erp_data.shape[1]) * 1.0 / SR - pre_onset
    279
    280      # if pre != pre_onset
    281      if pre_discard > 0:
    282          npoints = int(pre_discard * SR)
    283          time_points = time_points[npoints:]
    284          erp_data = erp_data[:, npoints:]
    285
    286      # select only time points of interest (if post is provided)
    287      if post is not None:
    288          npoints = int(duration * SR)
    289          time_points = time_points[:npoints]
    290          erp_data = erp_data[:, :npoints]
    291
    292      # compute mean signal timecourse accross trials
    293      erp_mean = N.mean(erp_data, axis=0)
    294
    295      # give sane default
    296      if errtype is None:
    297          errtype = []
    298      if not isinstance(errtype, list):
    299          errtype = [errtype]
    300
    301      for et in errtype:
    302          # compute error per datapoint
    303          if et in ['ste', 'ci95']:
    304              erp_stderr = erp_data.std(axis=0) / N.sqrt(len(erp_data))
    305              if et == 'ci95':
    306                  erp_stderr *= 1.96
    307          elif et == 'std':
    308              erp_stderr = erp_data.std(axis=0)
    309          else:
    310              raise ValueError, "Unknown error type '%s'" % errtype
    311
    312          time_points2w = N.hstack((time_points, time_points[::-1]))
    313
    314          error_top = erp_mean + erp_stderr
    315          error_bottom = erp_mean - erp_stderr
    316          error2w = N.hstack((error_top, error_bottom[::-1]))
    317
    318          if errcolor is None:
    319              errcolor = color
    320
    321          # plot error margin
    322          pfill = ax.fill(time_points2w, error2w,
    323                          edgecolor=errcolor, facecolor=errcolor, alpha=0.2,
    324                          zorder=3)
    325
    326      # plot mean signal timecourse
    327      ax.plot(time_points, erp_mean, lw=2, color=color, zorder=4,
    328              *args, **kwargs)
    329  #    ax.xaxis.set_major_locator(P.MaxNLocator(4))
    330      return erp_mean
    331
    332
    333 -def plotERPs(erps, data=None, ax=None, pre=0.2, post=None,
    334               pre_onset=None,
    335               xlabel='time (s)', ylabel='$\mu V$',
    336               ylim=None, ymult=1.0, legend=None,
    337               xlformat='%4g', ylformat='%4g',
    338               loffset=10, alinewidth=2,
    339               **kwargs):
    340      """Plot multiple ERPs on a new figure.
    341
    342      :Parameters:
    343        erps : list of tuples
    344          List of definitions of ERPs. Each tuple should consist of
    345          (label, color, onsets) or a dictionary which defines,
    346          label, color, onsets, data. Data provided in dictionary overrides
    347          'common' data provided in the next argument ``data``
    348        data
    349          Data for ERPs to be derived from 1D (samples)
    350        ax
    351          Where to draw (e.g. subplot instance). If None, new figure is
    352          created
    353        pre : float
    354          Duration (seconds) to be plotted prior to onset
    355        pre_onset : float or None
    356          If data is already in epochs (2D) then pre_onset provides information
    357          on how many seconds pre-stimulus were used to generate them. If None,
    358          then pre_onset = pre
    359        post : float or None
    360          Duration (seconds) to be plotted after the onset. If any data is
    361          provided with onsets, it can't be None. If None -- plots all time
    362          points after onsets
    363        ymult : float
    364          Multiplier for the values. E.g. if negative-up ERP plot is needed:
    365          provide ymult=-1.0
    366        xlformat : basestring
    367          Format of the x ticks
    368        ylformat : basestring
    369          Format of the y ticks
    370        legend : basestring or None
    371          If not None, legend will be plotted with position argument
    372          provided in this argument
    373        loffset : int
    374          Offset in voxels for axes and tick labels. Different
    375          matplotlib frontends might have different opinions, thus
    376          offset value might need to be tuned specifically per frontend
    377        alinewidth : int
    378          Axis and ticks line width
    379        **kwargs
    380          Additional arguments provided to plotERP()
    381
    382
    383      :Examples:
    384          kwargs  = {'SR' : eeg.SR, 'pre_mean' : 0.2}
    385          fig = plotERPs((('60db', 'b', eeg.erp_onsets['60db']),
    386                           ('80db', 'r', eeg.erp_onsets['80db'])),
    387                          data[:, eeg.sensor_mapping['Cz']],
    388                          ax=fig.add_subplot(1,1,1,frame_on=False), pre=0.2,
    389                          post=0.6, **kwargs)
    390
    391          or
    392          fig = plotERPs((('60db', 'b', eeg.erp_onsets['60db']),
    393                            {'color': 'r',
    394                             'onsets': eeg.erp_onsets['80db'],
    395                             'data' : data[:, eeg.sensor_mapping['Cz']]}
    396                           ),
    397                          data[:, eeg.sensor_mapping['Cz']],
    398                          ax=fig.add_subplot(1,1,1,frame_on=False), pre=0.2,
    399                          post=0.6, **kwargs)
    400
    401      :Returns: current fig handler
    402      """
    403
    404      if ax is None:
    405          fig = P.figure(facecolor='white')
    406          fig.clf()
    407          ax = fig.add_subplot(111, frame_on=False)
    408      else:
    409          fig = P.gcf()
    410
    411      # We don't want original axis being on
    412      ax.axison = True
    413
    414      labels = []
    415      for erp_def in erps:
    416          plot_data = data
    417          params = {'ymult' : ymult}
    418
    419          # provide custom parameters per ERP
    420          if isinstance(erp_def, tuple) and len(erp_def) == 3:
    421              params.update(
    422                  {'label': erp_def[0],
    423                   'color': erp_def[1],
    424                   'onsets': erp_def[2]})
    425          elif isinstance(erp_def, dict):
    426              plot_data = erp_def.pop('data', None)
    427              params.update(erp_def)
    428
    429          labels.append(params.get('label', ''))
    430
    431          # absorb common parameters
    432          params.update(kwargs)
    433
    434          if plot_data is None:
    435              raise ValueError, "Channel %s got no data provided" \
    436                    % params.get('label', 'UNKNOWN')
    437
    438
    439          plotERP(plot_data, pre=pre, pre_onset=pre_onset, post=post, ax=ax,
    440                  **params)
    441          #             plot_kwargs={'label':label})
    442
    443          if isinstance(erp_def, dict):
    444              erp_def['data'] = plot_data # return it back
    445
    446      props = dict(color='black',
    447                   linewidth=alinewidth, markeredgewidth=alinewidth,
    448                   zorder=1, offset=loffset)
    449
    450      def set_limits():
    451          """Helper to set x and y limits"""
    452          ax.set_xlim( (-pre, post) )
    453          if ylim != None:
    454              ax.set_ylim(*ylim)
    455
    456      set_limits()
    457      _make_centeredaxis(ax, 0, ai=0, label=xlabel, **props)
    458      set_limits()
    459      _make_centeredaxis(ax, 0, ai=1, mult=N.sign(ymult), label=ylabel, **props)
    460
    461      ax.yaxis.set_major_locator(P.NullLocator())
    462      ax.xaxis.set_major_locator(P.NullLocator())
    463
    464      # legend obscures plotting a bit... seems to be plotting
    465      # everything twice. Thus disabled by default
    466      if legend is not None and N.any(N.array(labels) != ''):
    467          P.legend(labels, loc=legend)
    468
    469      fig.canvas.draw()
    470      return fig
    471