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
43
44 return mlt.offset_copy(ax.transData, ax, x=x, y=y, units='dots')
45 elif 'BlendedAffine2D' in d:
46
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
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
128 args = [ [x[1], x[0]] for x in args ]
129
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:
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
154 lx, ly,
155 label,
156 horizontalalignment=lhalignment,
157 verticalalignment=lvalignment, fontsize=14,
158
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
219 pre_discard = 0
220
221 if onsets is not None:
222 if post is None:
223 raise ValueError, \
224 "Duration post onsets must be provided if onsets are given"
225
226 duration = pre + post
227
228
229 bcm = BoxcarMapper(onsets,
230 boxlength = int(SR * duration),
231 offset = -int(SR * pre))
232 erp_data = bcm(data)
233
234
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 """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
465
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