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
217
218 pre_discard = 0
219
220 if onsets is not None:
221 if post is None:
222 raise ValueError, \
223 "Duration post onsets must be provided if onsets are given"
224
225 duration = pre + post
226
227
228 bcm = BoxcarMapper(onsets,
229 boxlength = int(SR * duration),
230 offset = -int(SR * pre))
231 erp_data = bcm(data)
232
233
234 pre_discard = 0
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
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
257 erp_data *= ymult
258
259
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
268 erp_baseline = N.mean(
269 erp_data[:, int((pre_onset-pre_mean)*SR):int(pre_onset*SR)])
270
271
272
273 erp_data = erp_data - erp_baseline
274
275
276
277
278 time_points = N.arange(erp_data.shape[1]) * 1.0 / SR - pre_onset
279
280
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
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
293 erp_mean = N.mean(erp_data, axis=0)
294
295
296 if errtype is None:
297 errtype = []
298 if not isinstance(errtype, list):
299 errtype = [errtype]
300
301 for et in errtype:
302
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
322 pfill = ax.fill(time_points2w, error2w,
323 edgecolor=errcolor, facecolor=errcolor, alpha=0.2,
324 zorder=3)
325
326
327 ax.plot(time_points, erp_mean, lw=2, color=color, zorder=4,
328 *args, **kwargs)
329
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
412 ax.axison = True
413
414 labels = []
415 for erp_def in erps:
416 plot_data = data
417 params = {'ymult' : ymult}
418
419
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
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
442
443 if isinstance(erp_def, dict):
444 erp_def['data'] = plot_data
445
446 props = dict(color='black',
447 linewidth=alinewidth, markeredgewidth=alinewidth,
448 zorder=1, offset=loffset)
449
450 def set_limits():
451 ax.set_xlim( (-pre, post) )
452 if ylim != None:
453 ax.set_ylim(*ylim)
454
455 set_limits()
456 _make_centeredaxis(ax, 0, ai=0, label=xlabel, **props)
457 set_limits()
458 _make_centeredaxis(ax, 0, ai=1, mult=N.sign(ymult), label=ylabel, **props)
459
460 ax.yaxis.set_major_locator(P.NullLocator())
461 ax.xaxis.set_major_locator(P.NullLocator())
462
463
464
465 if legend is not None and N.any(N.array(labels) != ''):
466 P.legend(labels, loc=legend)
467
468 fig.canvas.draw()
469 return fig
470