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