Coverage for src/thunderfish/pulseplots.py: 0%
524 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 20:05 +0000
1"""
2Plot and save key steps in pulses.py for visualizing the alorithm.
3"""
5import numpy as np
6import matplotlib.pyplot as plt
8from scipy import stats
9from matplotlib import gridspec, ticker
10try:
11 from matplotlib.colors import colorConverter as cc
12except ImportError:
13 import matplotlib.colors as cc
14try:
15 from matplotlib.colors import to_hex
16except ImportError:
17 from matplotlib.colors import rgb2hex as to_hex
18from matplotlib.patches import ConnectionPatch, Rectangle
19from matplotlib.lines import Line2D
21import warnings
22def warn(*args, **kwargs):
23 """
24 Ignore all warnings.
25 """
26 pass
27warnings.warn=warn
30# plotting parameters and colors:
31cmap = plt.get_cmap("Dark2")
32c_g = cmap(0)
33c_o = cmap(1)
34c_grey = cmap(7)
35cmap_pts = [cmap(2), cmap(3)]
38def darker(color, saturation):
39 """ Make a color darker.
41 From bendalab/plottools package.
43 Parameters
44 ----------
45 color: dict or matplotlib color spec
46 A matplotlib color (hex string, name color string, rgb tuple)
47 or a dictionary with an 'color' or 'facecolor' key.
48 saturation: float
49 The smaller the saturation, the darker the returned color.
50 A saturation of 0 returns black.
51 A saturation of 1 leaves the color untouched.
52 A saturation of 2 returns white.
54 Returns
55 -------
56 color: string or dictionary
57 The darker color as a hexadecimal RGB string (e.g. '#rrggbb').
58 If `color` is a dictionary, a copy of the dictionary is returned
59 with the value of 'color' or 'facecolor' set to the darker color.
60 """
61 try:
62 c = color['color']
63 cd = dict(**color)
64 cd['color'] = darker(c, saturation)
65 return cd
66 except (KeyError, TypeError):
67 try:
68 c = color['facecolor']
69 cd = dict(**color)
70 cd['facecolor'] = darker(c, saturation)
71 return cd
72 except (KeyError, TypeError):
73 if saturation > 2:
74 sauration = 2
75 if saturation > 1:
76 return lighter(color, 2.0-saturation)
77 if saturation < 0:
78 saturation = 0
79 r, g, b = cc.to_rgb(color)
80 rd = r*saturation
81 gd = g*saturation
82 bd = b*saturation
83 return to_hex((rd, gd, bd)).upper()
86def lighter(color, lightness):
87 """Make a color lighter
89 From bendalab/plottools package.
91 Parameters
92 ----------
93 color: dict or matplotlib color spec
94 A matplotlib color (hex string, name color string, rgb tuple)
95 or a dictionary with an 'color' or 'facecolor' key.
96 lightness: float
97 The smaller the lightness, the lighter the returned color.
98 A lightness of 0 returns white.
99 A lightness of 1 leaves the color untouched.
100 A lightness of 2 returns black.
102 Returns
103 -------
104 color: string or dict
105 The lighter color as a hexadecimal RGB string (e.g. '#rrggbb').
106 If `color` is a dictionary, a copy of the dictionary is returned
107 with the value of 'color' or 'facecolor' set to the lighter color.
108 """
109 try:
110 c = color['color']
111 cd = dict(**color)
112 cd['color'] = lighter(c, lightness)
113 return cd
114 except (KeyError, TypeError):
115 try:
116 c = color['facecolor']
117 cd = dict(**color)
118 cd['facecolor'] = lighter(c, lightness)
119 return cd
120 except (KeyError, TypeError):
121 if lightness > 2:
122 lightness = 2
123 if lightness > 1:
124 return darker(color, 2.0-lightness)
125 if lightness < 0:
126 lightness = 0
127 r, g, b = cc.to_rgb(color)
128 rl = r + (1.0-lightness)*(1.0 - r)
129 gl = g + (1.0-lightness)*(1.0 - g)
130 bl = b + (1.0-lightness)*(1.0 - b)
131 return to_hex((rl, gl, bl)).upper()
134def xscalebar(ax, x, y, width, wunit=None, wformat=None, ha='left', va='bottom',
135 lw=None, color=None, capsize=None, clw=None, **kwargs):
136 """Horizontal scale bar with label.
138 From bendalab/plottools package.
140 Parameters
141 ----------
142 ax: matplotlib axes
143 Axes where to draw the scale bar.
144 x: float
145 x-coordinate where to draw the scale bar in relative units of the axes.
146 y: float
147 y-coordinate where to draw the scale bar in relative units of the axes.
148 width: float
149 Length of the scale bar in units of the data's x-values.
150 wunit: string or None
151 Optional unit of the data's x-values.
152 wformat: string or None
153 Optional format string for formatting the label of the scale bar
154 or simply a string used for labeling the scale bar.
155 ha: 'left', 'right', or 'center'
156 Scale bar aligned left, right, or centered to (x, y)
157 va: 'top' or 'bottom'
158 Label of the scale bar either above or below the scale bar.
159 lw: int, float, None
160 Line width of the scale bar.
161 color: matplotlib color
162 Color of the scalebar.
163 capsize: float or None
164 If larger then zero draw cap lines at the ends of the bar.
165 The length of the lines is given in points (same unit as linewidth).
166 clw: int, float, None
167 Line width of the cap lines.
168 kwargs: key-word arguments
169 Passed on to `ax.text()` used to print the scale bar label.
170 """
171 ax.autoscale(False)
172 # ax dimensions:
173 pixelx = np.abs(np.diff(ax.get_window_extent().get_points()[:,0]))[0]
174 pixely = np.abs(np.diff(ax.get_window_extent().get_points()[:,1]))[0]
175 xmin, xmax = ax.get_xlim()
176 ymin, ymax = ax.get_ylim()
177 unitx = xmax - xmin
178 unity = ymax - ymin
179 dxu = np.abs(unitx)/pixelx
180 dyu = np.abs(unity)/pixely
181 # transform x, y from relative units to axis units:
182 x = xmin + x*unitx
183 y = ymin + y*unity
184 # bar length:
185 if wformat is None:
186 wformat = '%.0f'
187 if width < 1.0:
188 wformat = '%.1f'
189 try:
190 ls = wformat % width
191 width = float(ls)
192 except TypeError:
193 ls = wformat
194 # bar:
195 if ha == 'left':
196 x0 = x
197 x1 = x+width
198 elif ha == 'right':
199 x0 = x-width
200 x1 = x
201 else:
202 x0 = x-0.5*width
203 x1 = x+0.5*width
204 # line width:
205 if lw is None:
206 lw = 2
207 # color:
208 if color is None:
209 color = 'k'
210 # scalebar:
211 lh = ax.plot([x0, x1], [y, y], '-', color=color, lw=lw,
212 solid_capstyle='butt', clip_on=False)
213 # get y position of line in figure pixel coordinates:
214 ly = np.array(lh[0].get_window_extent(ax.get_figure().canvas.get_renderer()))[0,1]
215 # caps:
216 if capsize is None:
217 capsize = 0
218 if clw is None:
219 clw = 0.5
220 if capsize > 0.0:
221 dy = capsize*dyu
222 ax.plot([x0, x0], [y-dy, y+dy], '-', color=color, lw=clw,
223 solid_capstyle='butt', clip_on=False)
224 ax.plot([x1, x1], [y-dy, y+dy], '-', color=color, lw=clw,
225 solid_capstyle='butt', clip_on=False)
226 # label:
227 if wunit:
228 ls += u'\u2009%s' % wunit
229 if va == 'top':
230 th = ax.text(0.5*(x0+x1), y, ls, clip_on=False,
231 ha='center', va='bottom', **kwargs)
232 # get y coordinate of text bottom in figure pixel coordinates:
233 ty = np.array(th.get_window_extent(ax.get_figure().canvas.get_renderer()))[0,1]
234 dty = ly+0.5*lw + 2.0 - ty
235 else:
236 th = ax.text(0.5*(x0+x1), y, ls, clip_on=False,
237 ha='center', va='top', **kwargs)
238 # get y coordinate of text bottom in figure pixel coordinates:
239 ty = np.array(th.get_window_extent(ax.get_figure().canvas.get_renderer()))[1,1]
240 dty = ly-0.5*lw - 2.0 - ty
241 th.set_position((0.5*(x0+x1), y+dyu*dty))
242 return x0, x1, y
245def yscalebar(ax, x, y, height, hunit=None, hformat=None, ha='left', va='bottom',
246 lw=None, color=None, capsize=None, clw=None, **kwargs):
248 """Vertical scale bar with label.
250 From bendalab/plottools package.
252 Parameters
253 ----------
254 ax: matplotlib axes
255 Axes where to draw the scale bar.
256 x: float
257 x-coordinate where to draw the scale bar in relative units of the axes.
258 y: float
259 y-coordinate where to draw the scale bar in relative units of the axes.
260 height: float
261 Length of the scale bar in units of the data's y-values.
262 hunit: string
263 Unit of the data's y-values.
264 hformat: string or None
265 Optional format string for formatting the label of the scale bar
266 or simply a string used for labeling the scale bar.
267 ha: 'left' or 'right'
268 Label of the scale bar either to the left or to the right
269 of the scale bar.
270 va: 'top', 'bottom', or 'center'
271 Scale bar aligned above, below, or centered on (x, y).
272 lw: int, float, None
273 Line width of the scale bar.
274 color: matplotlib color
275 Color of the scalebar.
276 capsize: float or None
277 If larger then zero draw cap lines at the ends of the bar.
278 The length of the lines is given in points (same unit as linewidth).
279 clw: int, float
280 Line width of the cap lines.
281 kwargs: key-word arguments
282 Passed on to `ax.text()` used to print the scale bar label.
283 """
285 ax.autoscale(False)
286 # ax dimensions:
287 pixelx = np.abs(np.diff(ax.get_window_extent().get_points()[:,0]))[0]
288 pixely = np.abs(np.diff(ax.get_window_extent().get_points()[:,1]))[0]
289 xmin, xmax = ax.get_xlim()
290 ymin, ymax = ax.get_ylim()
291 unitx = xmax - xmin
292 unity = ymax - ymin
293 dxu = np.abs(unitx)/pixelx
294 dyu = np.abs(unity)/pixely
295 # transform x, y from relative units to axis units:
296 x = xmin + x*unitx
297 y = ymin + y*unity
298 # bar length:
299 if hformat is None:
300 hformat = '%.0f'
301 if height < 1.0:
302 hformat = '%.1f'
303 try:
304 ls = hformat % height
305 width = float(ls)
306 except TypeError:
307 ls = hformat
308 # bar:
309 if va == 'bottom':
310 y0 = y
311 y1 = y+height
312 elif va == 'top':
313 y0 = y-height
314 y1 = y
315 else:
316 y0 = y-0.5*height
317 y1 = y+0.5*height
318 # line width:
319 if lw is None:
320 lw = 2
321 # color:
322 if color is None:
323 color = 'k'
324 # scalebar:
325 lh = ax.plot([x, x], [y0, y1], '-', color=color, lw=lw,
326 solid_capstyle='butt', clip_on=False)
327 # get x position of line in figure pixel coordinates:
328 lx = np.array(lh[0].get_window_extent(ax.get_figure().canvas.get_renderer()))[0,0]
329 # caps:
330 if capsize is None:
331 capsize = 0
332 if clw is None:
333 clw = 0.5
334 if capsize > 0.0:
335 dx = capsize*dxu
336 ax.plot([x-dx, x+dx], [y0, y0], '-', color=color, lw=clw, solid_capstyle='butt',
337 clip_on=False)
338 ax.plot([x-dx, x+dx], [y1, y1], '-', color=color, lw=clw, solid_capstyle='butt',
339 clip_on=False)
340 # label:
341 if hunit:
342 ls += u'\u2009%s' % hunit
343 if ha == 'right':
344 th = ax.text(x, 0.5*(y0+y1), ls, clip_on=False, rotation=90.0,
345 ha='left', va='center', **kwargs)
346 # get x coordinate of text bottom in figure pixel coordinates:
347 tx = np.array(th.get_window_extent(ax.get_figure().canvas.get_renderer()))[0,0]
348 dtx = lx+0.5*lw + 2.0 - tx
349 else:
350 th = ax.text(x, 0.5*(y0+y1), ls, clip_on=False, rotation=90.0,
351 ha='right', va='center', **kwargs)
352 # get x coordinate of text bottom in figure pixel coordinates:
353 tx = np.array(th.get_window_extent(ax.get_figure().canvas.get_renderer()))[1,0]
354 dtx = lx-0.5*lw - 1.0 - tx
355 th.set_position((x+dxu*dtx, 0.5*(y0+y1)))
356 return x, y0, y1
359def arrowed_spines(ax, ms=10):
360 """ Spine with arrow on the y-axis of a plot.
362 Parameters
363 ----------
364 ax : matplotlib figure axis
365 Axis on which the arrow should be plot.
366 """
367 xmin, xmax = ax.get_xlim()
368 ymin, ymax = ax.get_ylim()
369 ax.scatter([xmin], [ymax], s=ms, marker='^', clip_on=False, color='k')
370 ax.set_xlim(xmin, xmax)
371 ax.set_ylim(ymin, ymax)
374def loghist(ax, x, bmin, bmax, n, c, orientation='vertical', label=''):
375 """ Plot histogram with logarithmic scale.
377 Parameters
378 ----------
379 ax : matplotlib axis
380 Axis to plot the histogram on.
381 x : numpy array
382 Input data for histogram.
383 bmin : float
384 Minimum value for the histogram bins.
385 bmax : float
386 Maximum value for the histogram bins.
387 n : int
388 Number of bins.
389 c : matplotlib color
390 Color of histogram.
391 orientation : string (optional)
392 Histogram orientation.
393 Defaults to 'vertical'.
394 label : string (optional)
395 Label for x.
396 Defaults to '' (no label).
398 Returns
399 -------
400 n : array
401 The values of the histogram bins.
402 bins : array
403 The edges of the bins.
404 patches : BarContainer
405 Container of individual artists used to create the histogram.
406 """
407 return ax.hist(x, bins=np.exp(np.linspace(np.log(bmin), np.log(bmax), n)),
408 color=c, orientation=orientation, label=label)
411def plot_all(data, rate, eod_p_times, eod_tr_times, mean_eods):
412 """Quick way to view the output of extract_pulsefish in a single plot.
414 Parameters
415 ----------
416 data: array
417 Recording data.
418 eod_p_times: array of ints
419 EOD peak indices.
420 eod_tr_times: array of ints
421 EOD trough indices.
422 rate: float
423 Sampling rate.
424 mean_eods: list of numpy arrays
425 Mean EODs of each pulsefish found in the recording.
427 Returns
428 -------
429 fig: matplotlib.Figure
430 The figure with the plots.
431 """
432 if len(eod_p_times) > 0:
433 fig = plt.figure(figsize=(10, 5), layout='tight')
434 gs = gridspec.GridSpec(2, len(eod_p_times))
435 ax = fig.add_subplot(gs[0,:])
436 ax.plot(np.arange(len(data))/rate, data, c='k', alpha=0.3)
438 for i, (pt, tt) in enumerate(zip(eod_p_times, eod_tr_times)):
439 ax.plot(pt, data[(pt*rate).astype('int')], 'o', label=i+1, ms=10, c=cmap(i))
440 ax.plot(tt, data[(tt*rate).astype('int')], 'o', label=i+1, ms=10, c=cmap(i))
442 ax.set_xlabel('time [s]')
443 ax.set_ylabel('amplitude [V]')
445 for i, m in enumerate(mean_eods):
446 ax = fig.add_subplot(gs[1,i])
447 ax.plot(1000*m[0], 1000*m[1], c='k')
449 ax.fill_between(1000*m[0], 1000*(m[1]-m[2]), 1000*(m[1]+m[2]), color=cmap(i))
450 ax.set_xlabel('time [ms]')
451 ax.set_ylabel('amplitude [mV]')
452 else:
453 fig, ax = plt.subplots(figsize=(10, 5), layout='tight')
454 ax.plot(np.arange(len(data))/rate, data, c='k', alpha=0.3)
455 return fig
458def plot_clustering(rate, eod_widths, eod_hights, eod_shapes, disc_masks, merge_masks):
459 """Plot all clustering steps.
461 Plot clustering steps on width, height and shape. Then plot the remaining EODs after
462 the EOD assessment step and the EODs after the merge step.
464 Parameters
465 ----------
466 rate : float
467 Sampling rate of EOD snippets.
468 eod_widths : list of three 1D numpy arrays
469 The first list entry gives the unique labels of all width clusters as a list of ints.
470 The second list entry gives the width values for each EOD in samples as a
471 1D numpy array of ints.
472 The third list entry gives the width labels for each EOD as a 1D numpy array of ints.
473 eod_hights : nested lists (2 layers) of three 1D numpy arrays
474 The first list entry gives the unique labels of all height clusters as a list of ints
475 for each width cluster.
476 The second list entry gives the height values for each EOD as a 1D numpy array
477 of floats for each width cluster.
478 The third list entry gives the height labels for each EOD as a 1D numpy array
479 of ints for each width cluster.
480 eod_shapes : nested lists (3 layers) of three 1D numpy arrays
481 The first list entry gives the raw EOD snippets as a 2D numpy array for each
482 height cluster in a width cluster.
483 The second list entry gives the snippet PCA values for each EOD as a 2D numpy array
484 of floats for each height cluster in a width cluster.
485 The third list entry gives the shape labels for each EOD as a 1D numpy array of ints
486 for each height cluster in a width cluster.
487 disc_masks : Nested lists (two layers) of 1D numpy arrays
488 The masks of EODs that are discarded by the discarding step of the algorithm.
489 The masks are 1D boolean arrays where
490 instances that are set to True are discarded by the algorithm. Discarding masks
491 are saved in nested lists that represent the width and height clusters.
492 merge_masks : Nested lists (two layers) of 2D numpy arrays
493 The masks of EODs that are discarded by the merging step of the algorithm.
494 The masks are 2D boolean arrays where
495 for each sample point `i` either `merge_mask[i,0]` or `merge_mask[i,1]` is set to True.
496 Here, merge_mask[:,0] represents the
497 peak-centered clusters and `merge_mask[:,1]` represents the trough-centered clusters.
498 Merge masks are saved in nested lists
499 that represent the width and height clusters.
501 Returns
502 -------
503 fig: matplotlib.Figure
504 The figure with the plots.
505 """
506 # create figure + transparant figure.
507 fig = plt.figure(figsize=(12, 7))
508 transFigure = fig.transFigure.inverted()
510 # set up the figure layout
511 outer = gridspec.GridSpec(1, 5, width_ratios=[1, 1, 2, 1, 2], left=0.05, right=0.95)
513 # set titles for each clustering step
514 titles = ['1. Widths', '2. Heights', '3. Shape', '4. Pulse EODs', '5. Merge']
515 for i, title in enumerate(titles):
516 title_ax = gridspec.GridSpecFromSubplotSpec(1, 1, subplot_spec = outer[i])
517 ax = fig.add_subplot(title_ax[0])
518 ax.text(0, 110, title, ha='center', va='bottom', clip_on=False)
519 ax.set_xlim(-100, 100)
520 ax.set_ylim(-100, 100)
521 ax.axis('off')
523 # compute sizes for each axis
524 w_size = 1
525 h_size = len(eod_hights[1])
527 shape_size = np.sum([len(sl) for sl in eod_shapes[0]])
529 # count required axes sized for the last two plot columns.
530 disc_size = 0
531 merge_size= 0
532 for shapelabel, dmasks, mmasks in zip(eod_shapes[2], disc_masks, merge_masks):
533 for sl, dm, mm in zip(shapelabel, dmasks, mmasks):
534 uld1 = np.unique((sl[0]+1)*np.invert(dm[0]))
535 uld2 = np.unique((sl[1]+1)*np.invert(dm[1]))
536 disc_size = disc_size+len(uld1[uld1>0])+len(uld2[uld2>0])
538 uld1 = np.unique((sl[0]+1)*mm[0])
539 uld2 = np.unique((sl[1]+1)*mm[1])
540 merge_size = merge_size+len(uld1[uld1>0])+len(uld2[uld2>0])
542 # set counters to keep track of the plot axes
543 disc_block = 0
544 merge_block = 0
545 shape_count = 0
547 # create all axes
548 width_hist_ax = gridspec.GridSpecFromSubplotSpec(w_size, 1, subplot_spec = outer[0])
549 hight_hist_ax = gridspec.GridSpecFromSubplotSpec(h_size, 1, subplot_spec = outer[1])
550 shape_ax = gridspec.GridSpecFromSubplotSpec(shape_size, 1, subplot_spec = outer[2])
551 shape_windows = [gridspec.GridSpecFromSubplotSpec(2, 2, hspace=0.0, wspace=0.0,
552 subplot_spec=shape_ax[i])
553 for i in range(shape_size)]
555 EOD_delete_ax = gridspec.GridSpecFromSubplotSpec(disc_size, 1, subplot_spec=outer[3])
556 EOD_merge_ax = gridspec.GridSpecFromSubplotSpec(merge_size, 1, subplot_spec=outer[4])
558 # plot width labels histogram
559 ax1 = fig.add_subplot(width_hist_ax[0])
560 # set axes features.
561 ax1.set_xscale('log')
562 ax1.spines['top'].set_visible(False)
563 ax1.spines['right'].set_visible(False)
564 ax1.spines['bottom'].set_visible(False)
565 ax1.axes.xaxis.set_visible(False)
566 ax1.set_yticklabels([])
568 # indices for plot colors (dark to light)
569 colidxsw = -np.linspace(-1.25, -0.5, h_size)
571 for i, (wl, colw, uhl, eod_h, eod_h_labs, w_snip, w_feat, w_lab, w_dm, w_mm) in enumerate(zip(eod_widths[0], colidxsw, eod_hights[0], eod_hights[1], eod_hights[2], eod_shapes[0], eod_shapes[1], eod_shapes[2], disc_masks, merge_masks)):
573 # plot width hist
574 hw, _, _ = ax1.hist(eod_widths[1][eod_widths[2]==wl],
575 bins=np.linspace(np.min(eod_widths[1]), np.max(eod_widths[1]), 100),
576 color=lighter(c_o, colw), orientation='horizontal')
578 # set arrow when the last hist is plot so the size of the axes are known.
579 if i == h_size-1:
580 arrowed_spines(ax1, ms=20)
582 # determine total size of the hight histograms now.
583 my, b = np.histogram(eod_h, bins=np.exp(np.linspace(np.min(np.log(eod_h)),
584 np.max(np.log(eod_h)), 100)))
585 maxy = np.max(my)
587 # set axes features for hight hist.
588 ax2 = fig.add_subplot(hight_hist_ax[h_size-i-1])
589 ax2.set_xscale('log')
590 ax2.spines['top'].set_visible(False)
591 ax2.spines['right'].set_visible(False)
592 ax2.spines['bottom'].set_visible(False)
593 ax2.set_xlim(0.9, maxy)
594 ax2.axes.xaxis.set_visible(False)
595 ax2.set_yscale('log')
596 ax2.yaxis.set_major_formatter(ticker.NullFormatter())
597 ax2.yaxis.set_minor_formatter(ticker.NullFormatter())
599 # define colors for plots
600 colidxsh = -np.linspace(-1.25, -0.5, len(uhl))
602 for n, (hl, hcol, snippets, features, labels, dmasks, mmasks) in enumerate(zip(uhl, colidxsh, w_snip, w_feat, w_lab, w_dm, w_mm)):
604 hh, _, _ = loghist(ax2, eod_h[eod_h_labs==hl], np.min(eod_h), np.max(eod_h), 100,
605 lighter(c_g, hcol), orientation='horizontal')
607 # set arrow spines only on last plot
608 if n == len(uhl)-1:
609 arrowed_spines(ax2, ms=10)
611 # plot line from the width histogram to the height histogram.
612 if n == 0:
613 coord1 = transFigure.transform(ax1.transData.transform([np.median(hw[hw!=0]),
614 np.median(eod_widths[1][eod_widths[2]==wl])]))
615 coord2 = transFigure.transform(ax2.transData.transform([0.9, np.mean(eod_h)]))
616 line = Line2D((coord1[0], coord2[0]), (coord1[1], coord2[1]),
617 transform=fig.transFigure, color='grey', linewidth=0.5)
618 fig.lines.append(line)
620 # compute sizes of the eod_discarding and merge steps
621 s1 = np.unique((labels[0]+1)*(~dmasks[0]))
622 s2 = np.unique((labels[1]+1)*(~dmasks[1]))
623 disc_block = disc_block + len(s1[s1>0]) + len(s2[s2>0])
625 s1 = np.unique((labels[0]+1)*(mmasks[0]))
626 s2 = np.unique((labels[1]+1)*(mmasks[1]))
627 merge_block = merge_block + len(s1[s1>0]) + len(s2[s2>0])
629 axs = []
630 disc_count = 0
631 merge_count = 0
633 # now plot the clusters for peak and trough centerings
634 for pt, cmap_pt in zip([0, 1], cmap_pts):
636 ax3 = fig.add_subplot(shape_windows[shape_size-1-shape_count][pt,0])
637 ax4 = fig.add_subplot(shape_windows[shape_size-1-shape_count][pt,1])
639 # remove axes
640 ax3.axes.xaxis.set_visible(False)
641 ax4.axes.yaxis.set_visible(False)
642 ax3.axes.yaxis.set_visible(False)
643 ax4.axes.xaxis.set_visible(False)
645 # set color indices
646 colidxss = -np.linspace(-1.25, -0.5, len(np.unique(labels[pt][labels[pt]>=0])))
647 j=0
648 for c in np.unique(labels[pt]):
650 if c<0:
651 # plot noise features + snippets
652 ax3.plot(features[pt][labels[pt]==c,0], features[pt][labels[pt]==c,1],
653 '.', color='lightgrey', label='-1', rasterized=True)
654 ax4.plot(snippets[pt][labels[pt]==c].T, linewidth=0.1,
655 color='lightgrey', label='-1', rasterized=True)
656 else:
657 # plot cluster features and snippets
658 ax3.plot(features[pt][labels[pt]==c,0], features[pt][labels[pt]==c,1],
659 '.', color=lighter(cmap_pt, colidxss[j]), label=c,
660 rasterized=True)
661 ax4.plot(snippets[pt][labels[pt]==c].T, linewidth=0.1,
662 color=lighter(cmap_pt, colidxss[j]), label=c, rasterized=True)
664 # check if the current cluster is an EOD, if yes, plot it.
665 if np.sum(dmasks[pt][labels[pt]==c]) == 0:
667 ax = fig.add_subplot(EOD_delete_ax[disc_size-disc_block+disc_count])
668 ax.axis('off')
670 # plot mean EOD snippet
671 ax.plot(np.mean(snippets[pt][labels[pt]==c], axis=0),
672 color=lighter(cmap_pt, colidxss[j]))
673 disc_count = disc_count + 1
675 # match colors and draw line..
676 coord1 = transFigure.transform(ax4.transData.transform([ax4.get_xlim()[1],
677 ax4.get_ylim()[0] + 0.5*(ax4.get_ylim()[1]-ax4.get_ylim()[0])]))
678 coord2 = transFigure.transform(ax.transData.transform([ax.get_xlim()[0],ax.get_ylim()[0] + 0.5*(ax.get_ylim()[1]-ax.get_ylim()[0])]))
679 line = Line2D((coord1[0], coord2[0]), (coord1[1], coord2[1]),
680 transform=fig.transFigure, color='grey',
681 linewidth=0.5)
682 fig.lines.append(line)
683 axs.append(ax)
685 # check if the current EOD survives the merge step
686 # if so, plot it.
687 if np.sum(mmasks[pt, labels[pt]==c])>0:
689 ax = fig.add_subplot(EOD_merge_ax[merge_size-merge_block+merge_count])
690 ax.axis('off')
692 ax.plot(np.mean(snippets[pt][labels[pt]==c], axis=0),
693 color=lighter(cmap_pt, colidxss[j]))
694 merge_count = merge_count + 1
696 j=j+1
698 if pt==0:
699 # draw line from hight cluster to EOD shape clusters.
700 coord1 = transFigure.transform(ax2.transData.transform([np.median(hh[hh!=0]),
701 np.median(eod_h[eod_h_labs==hl])]))
702 coord2 = transFigure.transform(ax3.transData.transform([ax3.get_xlim()[0],
703 ax3.get_ylim()[0]]))
704 line = Line2D((coord1[0], coord2[0]), (coord1[1], coord2[1]),
705 transform=fig.transFigure, color='grey', linewidth=0.5)
706 fig.lines.append(line)
708 shape_count = shape_count + 1
710 if len(axs)>0:
711 # plot lines that indicate the merged clusters.
712 coord1 = transFigure.transform(axs[0].transData.transform([axs[0].get_xlim()[1]+0.1*(axs[0].get_xlim()[1]-axs[0].get_xlim()[0]),
713 axs[0].get_ylim()[1]-0.25*(axs[0].get_ylim()[1]-axs[0].get_ylim()[0])]))
714 coord2 = transFigure.transform(axs[-1].transData.transform([axs[-1].get_xlim()[1]+0.1*(axs[-1].get_xlim()[1]-axs[-1].get_xlim()[0]),
715 axs[-1].get_ylim()[0]+0.25*(axs[-1].get_ylim()[1]-axs[-1].get_ylim()[0])]))
716 line = Line2D((coord1[0], coord2[0]), (coord1[1], coord2[1]),
717 transform=fig.transFigure, color='grey', linewidth=1)
718 fig.lines.append(line)
719 return fig
722def plot_bgm(x, means, variances, weights, use_log, labels, labels_am, xlabel):
723 """Plot a BGM clustering step either on EOD width or height.
725 Parameters
726 ----------
727 x: 1D array of float
728 BGM input values.
729 means: list of float
730 BGM Gaussian means
731 variances: list of float
732 BGM Gaussian variances.
733 weights: list of float
734 BGM Gaussian weights.
735 use_log: bool
736 True if the z-scored logarithm of the data was used as BGM input.
737 labels: 1D array of int
738 Labels defined by BGM model (before merging based on merge factor).
739 labels_am: 1D array of int
740 Labels defined by BGM model (after merging based on merge factor).
741 xlabel: str
742 Label for plot (defines the units of the BGM data).
744 Returns
745 -------
746 fig: matplotlib.Figure
747 The figure with the plots.
748 """
749 if xlabel.lower().startswith('width'):
750 ccol = c_o
751 elif xlabel.lower().startswith('height'):
752 ccol = c_g
753 else:
754 ccol = 'b'
756 # get the transform that was used as BGM input
757 if use_log:
758 x_transform = stats.zscore(np.log(x))
759 xplot = np.exp(np.linspace(np.log(np.min(x)), np.log(np.max(x)), 1000))
760 else:
761 x_transform = stats.zscore(x)
762 xplot = np.linspace(np.min(x), np.max(x), 1000)
764 # compute the x values and gaussians
765 x2 = np.linspace(np.min(x_transform), np.max(x_transform), 1000)
766 gaussians = []
767 gmax = 0
768 for i, (w, m, std) in enumerate(zip(weights, means, variances)):
769 gaus = np.sqrt(w*stats.norm.pdf(x2, m, np.sqrt(std)))
770 gaussians.append(gaus)
771 gmax = max(np.max(gaus), gmax)
773 # compute classes defined by gaussian intersections
774 classes = np.argmax(np.vstack(gaussians), axis=0)
776 # find the minimum of any gaussian that is within its class
777 gmin = 100
778 for i, c in enumerate(np.unique(classes)):
779 gmin=min(gmin, np.min(gaussians[c][classes==c]))
781 # set up the figure
782 fig, ax1 = plt.subplots(figsize=(8, 4.8), layout='tight')
783 fig_ysize = 4
784 ax2 = ax1.twinx()
785 ax1.spines['top'].set_visible(False)
786 ax2.spines['top'].set_visible(False)
787 ax1.set_ylabel('#')
788 ax2.set_ylabel('Likelihood')
789 ax2.set_yscale('log')
790 ax1.set_yscale('log')
791 if use_log:
792 ax1.set_xscale('log')
793 ax1.set_xlabel(xlabel)
795 # define colors for plotting gaussians
796 colidxs = -np.linspace(-1.25, -0.5, len(np.unique(classes)))
798 # plot the gaussians
799 for i, c in enumerate(np.unique(classes)):
800 ax2.plot(xplot, gaussians[c], c=lighter(c_grey, colidxs[i]),
801 linewidth=2, label=r'$N(\mu_%i, \sigma_%i)$'%(c, c))
803 # plot intersection lines
804 ax2.vlines(xplot[1:][np.diff(classes)!=0], 0, gmax/gmin,
805 color='k', linewidth=2, linestyle='--')
806 ax2.set_ylim(gmin, np.max(np.vstack(gaussians))*1.1)
808 # plot data distributions and classes
809 colidxs = -np.linspace(-1.25, -0.5, len(np.unique(labels)))
810 for i, l in enumerate(np.unique(labels)):
811 if use_log:
812 h, binn, _ = loghist(ax1, x[labels==l], np.min(x), np.max(x), 100,
813 lighter(ccol, colidxs[i]), label=r'$x_%i$'%l)
814 else:
815 h, binn, _ = ax1.hist(x[labels==l],
816 bins=np.linspace(np.min(x), np.max(x), 100),
817 color=lighter(ccol, colidxs[i]),
818 label=r'$x_%i$'%l)
820 # annotate merged clusters
821 for l in np.unique(labels_am):
822 maps = np.unique(labels[labels_am==l])
823 if len(maps) > 1:
824 x1 = x[labels==maps[0]]
825 x2 = x[labels==maps[1]]
827 ax2.plot([np.median(x1), np.median(x2)], [1.2*gmax, 1.2*gmax], c='k', clip_on=False)
828 ax2.plot([np.median(x1), np.median(x1)], [1.1*gmax, 1.2*gmax], c='k', clip_on=False)
829 ax2.plot([np.median(x2), np.median(x2)], [1.1*gmax, 1.2*gmax], c='k', clip_on=False)
830 ax2.annotate(r'$\frac{|{\tilde{x}_%i-\tilde{x}_%i}|}{max(\tilde{x}_%i, \tilde{x}_%i)} < \epsilon$' % (maps[0], maps[1], maps[0], maps[1]), [np.median(x1)*1.1, gmax*1.2], xytext=(10, 10), textcoords='offset points', fontsize=12, annotation_clip=False, ha='center')
832 # add legends and plot.
833 ax2.legend(loc='lower left', frameon=False, bbox_to_anchor=(-0.05, 1.3),
834 ncol=len(np.unique(classes)))
835 ax1.legend(loc='upper left', frameon=False, bbox_to_anchor=(-0.05, 1.3),
836 ncol=len(np.unique(labels)))
837 return fig
840def plot_feature_extraction(raw_snippets, normalized_snippets, features, labels, dt, pt):
841 """Plot clustering step on EOD shape.
843 Parameters
844 ----------
845 raw_snippets : 2D numpy array
846 Raw EOD snippets.
847 normalized_snippets : 2D numpy array
848 Normalized EOD snippets.
849 features : 2D numpy array
850 PCA values for each normalized EOD snippet.
851 labels : 1D numpy array of ints
852 Cluster labels.
853 dt : float
854 Sample interval of snippets.
855 pt : int
856 Set to 0 for peak-centered EODs and set to 1 for trough-centered EODs.
858 Returns
859 -------
860 fig: matplotlib.Figure
861 The figure with the plots.
862 """
863 ccol = cmap_pts[pt]
865 # set up the figure layout
866 fig = plt.figure(figsize=(((2+0.2)*4.8), 4.8))
867 outer = gridspec.GridSpec(1, 2, wspace=0.2, hspace=0)
869 x = np.arange(-dt*1000*raw_snippets.shape[1]/2, dt*1000*raw_snippets.shape[1]/2, dt*1000)
871 snip_ax = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec = outer[0], hspace=0.35)
872 pc_ax = gridspec.GridSpecFromSubplotSpec(features.shape[1]-1, features.shape[1]-1,
873 subplot_spec = outer[1], hspace=0, wspace=0)
875 # 3 plots: raw snippets, normalized, pcs.
876 ax_raw_snip = fig.add_subplot(snip_ax[0])
877 ax_normalized_snip = fig.add_subplot(snip_ax[1])
879 colidxs = -np.linspace(-1.25, -0.5, len(np.unique(labels[labels>=0])))
880 j=0
882 for c in np.unique(labels):
883 if c<0:
884 color='lightgrey'
885 else:
886 color = lighter(ccol, colidxs[j])
887 j=j+1
889 ax_raw_snip.plot(x, raw_snippets[labels==c].T, color=color, label='-1',
890 rasterized=True, alpha=0.25)
891 ax_normalized_snip.plot(x, normalized_snippets[labels==c].T, color=color, alpha=0.25)
892 ax_raw_snip.spines['top'].set_visible(False)
893 ax_raw_snip.spines['right'].set_visible(False)
894 ax_raw_snip.get_xaxis().set_ticklabels([])
895 ax_raw_snip.set_title('Raw snippets')
896 ax_raw_snip.set_ylabel('Amplitude [a.u.]')
897 ax_normalized_snip.spines['top'].set_visible(False)
898 ax_normalized_snip.spines['right'].set_visible(False)
899 ax_normalized_snip.set_title('Normalized snippets')
900 ax_normalized_snip.set_ylabel('Amplitude [a.u.]')
901 ax_normalized_snip.set_xlabel('Time [ms]')
903 ax_raw_snip.axis('off')
904 ax_normalized_snip.axis('off')
906 ax_overlay = fig.add_subplot(pc_ax[:,:])
907 ax_overlay.set_title('Features')
908 ax_overlay.axis('off')
910 for n in range(features.shape[1]):
911 for m in range(n):
912 ax = fig.add_subplot(pc_ax[n-1,m])
913 ax.scatter(features[labels==c,m], features[labels==c,n], marker='.',
914 color=color, alpha=0.25)
915 ax.set_xlim(np.min(features), np.max(features))
916 ax.set_ylim(np.min(features), np.max(features))
917 ax.get_xaxis().set_ticklabels([])
918 ax.get_yaxis().set_ticklabels([])
919 ax.get_xaxis().set_ticks([])
920 ax.get_yaxis().set_ticks([])
922 if m==0:
923 ax.set_ylabel('PC %i'%(n+1))
925 if n==features.shape[1]-1:
926 ax.set_xlabel('PC %i'%(m+1))
928 ax = fig.add_subplot(pc_ax[0,features.shape[1]-2])
929 ax.set_xlim(np.min(features), np.max(features))
930 ax.set_ylim(np.min(features), np.max(features))
932 size = max(1, int(np.ceil(-np.log10(np.max(features)-np.min(features)))))
933 wbar = np.floor((np.max(features)-np.min(features))*10**size)/10**size
935 # should be smaller than the actual thing! so like x% of it?
936 xscalebar(ax, 0, 0, wbar, wformat='%%.%if'%size)
937 yscalebar(ax, 0, 0, wbar, hformat='%%.%if'%size)
938 ax.axis('off')
939 return fig
942def plot_moving_fish(ws, dts, clusterss, ts, fishcounts, T, ignore_stepss):
943 """Plot moving fish detection step.
945 Parameters
946 ----------
947 ws : list of floats
948 Median width for each width cluster that the moving fish algorithm is computed on
949 (in seconds).
950 dts : list of floats
951 Sliding window size (in seconds) for each width cluster.
952 clusterss : list of 1D numpy int arrays
953 Cluster labels for each EOD cluster in a width cluster.
954 ts : list of 1D numpy float arrays
955 EOD emission times for each EOD in a width cluster.
956 fishcounts : list of lists
957 Sliding window timepoints and fishcounts for each width cluster.
958 T : float
959 Lenght of analyzed recording in seconds.
960 ignore_stepss : list of 1D int arrays
961 Mask for fishcounts that were ignored (ignored if True) in the moving_fish analysis.
963 Returns
964 -------
965 fig: matplotlib.Figure
966 The figure with the plots.
967 """
968 fig = plt.figure()
970 # create gridspec
971 outer = gridspec.GridSpec(len(ws), 1)
973 for i, (w, dt, clusters, t, fishcount, ignore_steps) in enumerate(zip(ws, dts, clusterss, ts, fishcounts, ignore_stepss)):
975 gs = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec = outer[i])
977 # axis for clusters
978 ax1 = fig.add_subplot(gs[0])
979 # axis for fishcount
980 ax2 = fig.add_subplot(gs[1])
982 # plot clusters as eventplot
983 for cnum, c in enumerate(np.unique(clusters[clusters>=0])):
984 ax1.eventplot(t[clusters==c], lineoffsets=cnum, linelengths=0.5, color=cmap(i))
985 cnum = cnum + 1
987 # Plot the sliding window
988 rect=Rectangle((0, -0.5), dt, cnum, linewidth=1, linestyle='--', edgecolor='k',
989 facecolor='none', clip_on=False)
990 ax1.add_patch(rect)
991 ax1.arrow(dt+0.1, -0.5, 0.5, 0, head_width=0.1, head_length=0.1, facecolor='k',
992 edgecolor='k')
994 # plot parameters
995 ax1.set_title(r'$\tilde{w}_%i = %.3f ms$'%(i, 1000*w))
996 ax1.set_ylabel('cluster #')
997 ax1.set_yticks(range(0, cnum))
998 ax1.set_xlabel('time')
999 ax1.set_xlim(0, T)
1000 ax1.axes.xaxis.set_visible(False)
1001 ax1.spines['bottom'].set_visible(False)
1002 ax1.spines['top'].set_visible(False)
1003 ax1.spines['right'].set_visible(False)
1004 ax1.spines['left'].set_visible(False)
1006 # plot for fishcount
1007 x = fishcount[0]
1008 y = fishcount[1]
1010 ax2 = fig.add_subplot(gs[1])
1011 ax2.spines['top'].set_visible(False)
1012 ax2.spines['right'].set_visible(False)
1013 ax2.spines['bottom'].set_visible(False)
1014 ax2.axes.xaxis.set_visible(False)
1016 yplot = np.array(y, dtype=float)
1017 ax2.plot(x+dt/2, yplot, linestyle='-', marker='.', c=cmap(i), alpha=0.25)
1018 yplot[ignore_steps.astype(bool)] = np.nan
1019 ax2.plot(x+dt/2, yplot, linestyle='-', marker='.', c=cmap(i))
1020 ax2.set_ylabel('Fish count')
1021 ax2.set_yticks(range(int(np.min(y)), 1+int(np.max(y))))
1022 ax2.set_xlim(0, T)
1024 if i < len(ws)-1:
1025 ax2.axes.xaxis.set_visible(False)
1026 else:
1027 ax2.axes.xaxis.set_visible(False)
1028 xscalebar(ax2, 1, 0, 1, wunit='s', ha='right')
1030 con = ConnectionPatch([0, -0.5], [dt/2, y[0]], "data", "data",
1031 axesA=ax1, axesB=ax2, color='k')
1032 ax2.add_artist(con)
1033 con = ConnectionPatch([dt, -0.5], [dt/2, y[0]], "data", "data",
1034 axesA=ax1, axesB=ax2, color='k')
1035 ax2.add_artist(con)
1037 ax2.set_xlim(0, T)
1038 return fig
1041if __name__ == '__main__':
1042 import sys
1043 from pathlib import Path
1044 from thunderlab.dataloader import load_data
1045 from thunderlab.configfile import ConfigFile
1046 from .bestwindow import analysis_window
1047 from .pulses import extract_pulsefish
1049 if len(sys.argv) < 2:
1050 print('pulseplots recording.wav')
1051 exit()
1052 data_file = Path(sys.argv[1])
1053 cfg = ConfigFile()
1054 raw_data, rate, unit, amax = load_data(data_file)
1055 data, idx0, idx1, clipped, min_clip, max_clip = \
1056 analysis_window(raw_data[:, 0], rate, amax, 'best', cfg)
1057 frate = 0.5e6 # TODO: make parameter
1058 mean_eods, eod_times, eod_peaktimes, pdict = \
1059 extract_pulsefish(data, rate, frate,
1060 verbose=0, plot_level=0,
1061 return_data=['all_eod_times',
1062 'peak_detection',
1063 'all_cluster_steps',
1064 'BGM_width',
1065 'BGM_height',
1066 'snippet_clusters',
1067 'eod_deletion',
1068 'masks',
1069 'moving_fish'])
1070 bgm = pdict['BGM_width']
1071 plot_bgm(bgm['x'], bgm['means'], bgm['variances'],
1072 bgm['weights'], bgm['use_log'], bgm['labels'],
1073 bgm['labels_bgm'], 'width [ms]')
1074 for wi in range(100):
1075 k = f'BGM_height_{wi}'
1076 if not k in pdict:
1077 break
1078 bgm = pdict[k]
1079 plot_bgm(bgm['x'], bgm['means'], bgm['variances'],
1080 bgm['weights'], bgm['use_log'], bgm['labels'],
1081 bgm['labels_bgm'], f'height (width cluster {wi}) [{unit}]')
1082 for hi in range(100):
1083 k = f'snippet_clusters_{wi}_{hi}_peak'
1084 if not k in pdict:
1085 break
1086 snip = pdict[k]
1087 plot_feature_extraction(snip['raw_snippets'], snip['snippets'],
1088 snip['features'], snip['clusters'], 1/pdict['i_rate'], 0)
1089 k = f'snippet_clusters_{wi}_{hi}_trough'
1090 snip = pdict[k]
1091 plot_feature_extraction(snip['raw_snippets'], snip['snippets'],
1092 snip['features'], snip['clusters'], 1/pdict['i_rate'], 1)
1093 plot_clustering(rate, pdict['EOD_widths'], pdict['EOD_heights'], pdict['EOD_shapes'],
1094 pdict['discarding_masks'], pdict['merge_masks'])
1095 plot_all(data, rate, pdict['eod_peaktimes'], pdict['eod_troughtimes'], mean_eods)
1096 move = pdict['moving_fish']
1097 plot_moving_fish(move['w'], move['dt'], move['clusters'], move['t'],
1098 move['fishcount'], move['T'], move['ignore_steps'])
1099 plt.show()