1
2
3
4
5
6
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
23
24 import matplotlib.transforms as mlt
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
42 return mlt.offset_copy(ax.transData, x=x, y=y, units='dots')
43 elif 'BlendedAffine2D' in d:
44
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
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
126 args = [ [x[1], x[0]] for x in args ]
127
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:
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
152 lx, ly,
153 label,
154 horizontalalignment=lhalignment,
155 verticalalignment=lvalignment, fontsize=14,
156
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:
217 if post is None:
218 raise ValueError, \
219 "Duration post onsets must be provided if onsets are given"
220
221 duration = pre + post
222
223
224 bcm = BoxcarMapper(onsets,
225 boxlength = int(SR * duration),
226 offset = -int(SR * pre))
227 erp_data = bcm(data)
228
229
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
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
253 erp_data *= ymult
254
255
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
264 erp_baseline = N.mean(
265 erp_data[:, int((pre_onset-pre_mean)*SR):int(pre_onset*SR)])
266
267
268
269 erp_data = erp_data - erp_baseline
270
271
272
273
274 time_points = N.arange(erp_data.shape[1]) * 1.0 / SR - pre_onset
275
276
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
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
289 erp_mean = N.mean(erp_data, axis=0)
290
291
292 if errtype is None:
293 errtype = []
294 if not isinstance(errtype, list):
295 errtype = [errtype]
296
297 for et in errtype:
298
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
318 pfill = ax.fill(time_points2w, error2w,
319 edgecolor=errcolor, facecolor=errcolor, alpha=0.2,
320 zorder=3)
321
322
323 ax.plot(time_points, erp_mean, lw=2, color=color, zorder=4,
324 *args, **kwargs)
325
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
408 ax.axison = True
409
410 labels = []
411 for erp_def in erps:
412 plot_data = data
413 params = {'ymult' : ymult}
414
415
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
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
438
439 if isinstance(erp_def, dict):
440 erp_def['data'] = plot_data
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
460
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