Coverage for src/thunderfish/pulses.py: 5%
790 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 20:05 +0000
1"""
2Extract and cluster EOD waverforms of pulse-type electric fish.
4## Main function
6- `extract_pulsefish()`: checks for pulse-type fish based on EOD amplitude and shape.
8"""
10import numpy as np
12try:
13 import matplotlib.pyplot as plt
14except ImportError:
15 pass
17from pathlib import Path
18from scipy import stats
19from scipy.interpolate import interp1d
20from scipy.signal import savgol_filter
21from sklearn.decomposition import PCA
22from sklearn.cluster import DBSCAN
23from sklearn.mixture import BayesianGaussianMixture
24from sklearn.metrics import pairwise_distances
26from thunderlab.eventdetection import detect_peaks, median_std_threshold
27from thunderlab.fourier import fourier_coeffs
29import warnings
30def warn(*args, **kwargs):
31 """
32 Ignore all warnings.
33 """
34 pass
35warnings.warn = warn
37try:
38 from numba import jit
39except ImportError:
40 def jit(*args, **kwargs):
41 def decorator_jit(func):
42 return func
43 return decorator_jit
46# upgrade numpy functions for backwards compatibility:
48if not hasattr(np, 'isin'):
49 np.isin = np.in1d
52def unique_counts(ar):
53 """ Find the unique elements of an array and their counts, ignoring shape.
55 The code is condensed from numpy version 1.17.0.
57 Parameters
58 ----------
59 ar : numpy array
60 Input array
62 Returns
63 -------
64 unique_vaulues : numpy array
65 Unique values in array ar.
66 unique_counts : numpy array
67 Number of instances for each unique value in ar.
68 """
69 try:
70 return np.unique(ar, return_counts=True)
71 except TypeError:
72 ar = np.asanyarray(ar).flatten()
73 ar.sort()
74 mask = np.empty(ar.shape, dtype=bool_)
75 mask[:1] = True
76 mask[1:] = ar[1:] != ar[:-1]
77 idx = np.concatenate(np.nonzero(mask) + ([mask.size],))
78 return ar[mask], np.diff(idx)
81###########################################################################
84def extract_pulsefish(data, rate, frate=0.5e6, width_factor_shape=3,
85 width_factor_wave=8, max_width=0.03, align='fourier',
86 verbose=0, plot_level=0, return_data=[]):
87 """Extract and cluster pulse-type fish EODs from data.
89 Takes recording data containing an unknown number of pulsefish and
90 extracts the mean EOD and EOD timepoints for each fish present in
91 the recording.
93 Parameters
94 ----------
95 data: 1-D array of float
96 The data to be analysed.
97 rate: float
98 Sampling rate of the data in Hertz.
99 frate: float
100 Sampling rate used for the returned waveform estimates.
101 width_factor_shape : float
102 Width multiplier used for EOD shape analysis.
103 EOD snippets are extracted based on width between the
104 peak and trough multiplied by the width factor.
105 width_factor_wave : float
106 Width multiplier used for wavefish detection.
107 max_width : float
108 Maximum width of returned EOD mean.
109 align: 'max', 'savgol', 'fourier'
110 Method used to align peak snippets. See extract_snippets() for details.
111 verbose : int
112 Verbosity level.
113 plot_level : int
114 Similar to verbosity levels, but with plots.
115 Only set to > 0 for debugging purposes.
116 return_data : list of str
117 Specify data that should be logged and returned in a dictionary. Each clustering
118 step has a specific keyword that results in adding different variables to the log dictionary.
119 Optional keys for return_data and the resulting additional key-value pairs to the log dictionary are:
121 - 'all_eod_times':
122 - 'all_times': list of two lists of float.
123 All peak (`all_times[0]`) and trough times (`all_times[1]`) extracted
124 by the peak detection algorithm. Times are given in seconds.
125 - 'eod_peaktimes': list of 1D arrays.
126 The timepoints in seconds of each unique extracted EOD cluster,
127 where each 1D array encodes one cluster.
128 - 'eod_troughtimes': list of 1D arrays.
129 The timepoints in seconds of each unique extracted EOD cluster,
130 where each 1D array encodes one cluster.
132 - 'peak_detection':
133 - "data": 1D numpy array of float.
134 Quadratically interpolated data which was used for peak detection.
135 - "peaks_1": 1D numpy array of int.
136 Peak indices on interpolated data after first peak detection step.
137 - "troughs_1": 1D numpy array of int.
138 Peak indices on interpolated data after first peak detection step.
139 - "peaks_2": 1D numpy array of int.
140 Peak indices on interpolated data after second peak detection step.
141 - "troughs_2": 1D numpy array of int.
142 Peak indices on interpolated data after second peak detection step.
143 - "peaks_3": 1D numpy array of int.
144 Peak indices on interpolated data after third peak detection step.
145 - "troughs_3": 1D numpy array of int.
146 Peak indices on interpolated data after third peak detection step.
147 - "peaks_4": 1D numpy array of int.
148 Peak indices on interpolated data after fourth peak detection step.
149 - "troughs_4": 1D numpy array of int.
150 Peak indices on interpolated data after fourth peak detection step.
152 - 'all_cluster_steps':
153 - 'rate': float.
154 Sampling rate of interpolated data.
155 - 'EOD_widths': list of three 1D numpy arrays.
156 The first list entry gives the unique labels of all width clusters
157 as a list of int.
158 The second list entry gives the width values for each EOD in samples
159 as a 1D numpy array of int.
160 The third list entry gives the width labels for each EOD
161 as a 1D numpy array of int.
162 - 'EOD_heights': nested lists (2 layers) of three 1D numpy arrays.
163 The first list entry gives the unique labels of all height clusters
164 as a list of int for each width cluster.
165 The second list entry gives the height values for each EOD
166 as a 1D numpy array of floats for each width cluster.
167 The third list entry gives the height labels for each EOD
168 as a 1D numpy array of int for each width cluster.
169 - 'EOD_shapes': nested lists (3 layers) of three 1D numpy arrays
170 The first list entry gives the raw EOD snippets as a 2D numpy array
171 for each height cluster in a width cluster.
172 The second list entry gives the snippet PCA values for each EOD
173 as a 2D numpy array of floats for each height cluster in a width cluster.
174 The third list entry gives the shape labels for each EOD as a 1D numpy array
175 of int for each height cluster in a width cluster.
176 - 'discarding_masks': Nested lists (two layers) of 1D numpy arrays.
177 The masks of EODs that are discarded by the discarding step of the algorithm.
178 The masks are 1D boolean arrays where instances that are set to True are
179 discarded by the algorithm. Discarding masks are saved in nested lists
180 that represent the width and height clusters.
181 - 'merge_masks': Nested lists (two layers) of 2D numpy arrays.
182 The masks of EODs that are discarded by the merging step of the algorithm.
183 The masks are 2D boolean arrays where for each sample point `i` either
184 `merge_mask[i,0]` or `merge_mask[i,1]` is set to True. Here, merge_mask[:,0]
185 represents the peak-centered clusters and `merge_mask[:,1]` represents the
186 trough-centered clusters. Merge masks are saved in nested lists that
187 represent the width and height clusters.
189 - 'BGM_width':
190 - 'BGM_width': dictionary
191 - 'x': 1D numpy array of float.
192 BGM input values (in this case the EOD widths),
193 - 'use_log': boolean.
194 True if the z-scored logarithm of the data was used as BGM input.
195 - 'means': 1D array of float
196 Means of the Gaussian fits.
197 - 'variances': 1D array of float
198 Variances of the Gaussian fits.
199 - 'weights': 1D array of float
200 Weights of the Gaussian fits.
201 - 'labels_bgm': 1D array of int
202 Labels resulting from BGM clustering.
203 - 'labels_split': 1D array of int
204 Labels after splitting at overlapping Gaussians.
205 - 'labels': 1D array of int
206 Final labels (after BGM clustering, splitting, and merging).
208 - 'BGM_height':
209 This key adds a new dictionary for each width cluster.
210 - 'BGM_height_*n*' : dictionary, where *n* defines the width cluster as an int.
211 - 'x': 1D numpy array of float.
212 BGM input values (in this case the EOD heights),
213 - 'use_log': boolean.
214 True if the z-scored logarithm of the data was used as BGM input.
215 - 'means': 1D array of float
216 Means of the Gaussian fits.
217 - 'variances': 1D array of float
218 Variances of the Gaussian fits.
219 - 'weights': 1D array of float
220 Weights of the Gaussian fits.
221 - 'labels_bgm': 1D array of int
222 Labels resulting from BGM clustering.
223 - 'labels_split': 1D array of int
224 Labels after splitting at overlapping Gaussians.
225 - 'labels': 1D array of int
226 Final labels (after BGM clustering, splitting, and merging).
228 - 'snippet_clusters':
229 This key adds a new dictionary for each height cluster.
230 - 'snippet_clusters*_n_m_p*' : dictionary, where *n* defines the width cluster
231 (int), *m* defines the height cluster (int) and *p* defines shape clustering
232 on peak or trough centered EOD snippets (str: 'peak' or 'trough').
233 - 'raw_snippets': 2D numpy array (nsamples, nfeatures).
234 Raw EOD snippets.
235 - 'snippets': 2D numpy array.
236 Normalized EOD snippets.
237 - 'features': 2D numpy array.(nsamples, nfeatures)
238 PCA values for each normalized EOD snippet.
239 - 'clusters': 1D numpy array of int.
240 Cluster labels.
241 - 'rate': float.
242 Sampling rate of snippets.
244 - 'artefacts':
245 This key adds two dictionaries for each width cluster,
246 one for peak-centered and one for trough-centered shape clusters:
247 - 'artefacts_peaks_*width*'
248 - 'artefacts_troughs_*width*'
249 where *width* (int) is the label of the width cluster.
250 Both of them have the same entries
251 where *cluster* (int) is the unique shape cluster label:
252 - 'mean_eod_*cluster*': 1D array of float
253 Mean EOD waveform of cluster.
254 - 'fft_power_*cluster*': 1D array of float
255 Power spectrum of mean EOD waveform.
256 - 'fft_freqs_*cluster*': 1D array of float
257 Frequencies of power spectrum.
258 - 'low_frequency': float
259 Frequency defining low frequencies.
260 - 'low_frequency_ratio_*cluster*': float
261 Ratio of low frequency power to total power.
262 - 'mask_*cluster*': float
264 - 'eod_deletion':
265 This key adds two dictionaries for each (peak centered) shape cluster,
266 where *cluster* (int) is the unique shape cluster label.
267 - 'mask_*cluster*' : list of four booleans.
268 The mask for each cluster discarding step.
269 The first instance represents the artefact masks, where artefacts
270 are set to True.
271 The second instance represents the unreliable cluster masks,
272 where unreliable clusters are set to True.
273 The third instance represents the wavefish masks, where wavefish
274 are set to True.
275 The fourth instance represents the sidepeak masks, where sidepeaks
276 are set to True.
277 - 'vals_*cluster*' : list of lists.
278 All variables that are used for each cluster deletion step.
279 The first instance is a list of two 1D numpy arrays: the mean EOD and
280 the FFT of that mean EOD.
281 The second instance is a 1D numpy array with all EOD width to ISI ratios.
282 The third instance is a list with three entries:
283 The first entry is a 1D numpy array zoomed out version of the mean EOD.
284 The second entry is a list of two 1D numpy arrays that define the peak
285 and trough indices of the zoomed out mean EOD.
286 The third entry contains a list of two values that represent the
287 peak-trough pair in the zoomed out mean EOD with the largest height
288 difference.
289 - 'rate' : float.
290 EOD snippet sampling rate.
292 - 'masks':
293 - 'masks' : 2D numpy array (4,N).
294 Each row contains masks for each EOD detected by the EOD peakdetection step.
295 The first row defines the artefact masks, the second row defines the
296 unreliable EOD masks,
297 the third row defines the wavefish masks and the fourth row defines
298 the sidepeak masks.
300 - 'moving_fish':
301 - 'moving_fish': dictionary.
302 - 'w' : list of float.
303 Median width for each width cluster that the moving fish algorithm is
304 computed on (in seconds).
305 - 'T' : list of float.
306 Lenght of analyzed recording for each width cluster (in seconds).
307 - 'dt' : list of float.
308 Sliding window size (in seconds) for each width cluster.
309 - 'clusters' : list of 1D numpy int arrays.
310 Cluster labels for each EOD cluster in a width cluster.
311 - 't' : list of 1D numpy float arrays.
312 EOD emission times for each EOD in a width cluster.
313 - 'fishcount' : list of lists.
314 Sliding window timepoints and fishcounts for each width cluster.
315 - 'ignore_steps' : list of 1D int arrays.
316 Mask for fishcounts that were ignored (ignored if True) in the
317 moving_fish analysis.
319 Returns
320 -------
321 mean_eods: list of 2D array of float
322 The average EOD for each detected fish. First column is time in seconds,
323 second column the mean eod, third column the standard deviation.
324 eod_times: list of 1D arrays
325 For each detected fish the times of EOD peaks or troughs in seconds.
326 Use these timepoints for EOD averaging.
327 eod_peaktimes: list of 1D arrays
328 For each detected fish the times of EOD peaks in seconds.
329 log_dict: dictionary
330 Dictionary with logged variables, where variables to log are specified
331 by `return_data`.
333 """
334 if verbose > 0:
335 print('')
336 if verbose > 1:
337 print(70*'#')
338 print('##### extract_pulsefish', 46*'#')
340 mean_eods, eod_times, eod_peaktimes = [], [], []
341 log_dict = {}
343 # interpolate:
344 i_rate = frate
345 try:
346 f = interp1d(np.arange(len(data))/rate, data, kind='quadratic')
347 i_data = f(np.arange(0.0, (len(data)-1)/rate, 1.0/i_rate))
348 except MemoryError:
349 i_rate = rate
350 i_data = data
351 log_dict['data'] = i_data # TODO: could be removed
352 log_dict['rate'] = i_rate # TODO: could be removed
353 log_dict['i_data'] = i_data
354 log_dict['i_rate'] = i_rate
356 # standard deviation of data in small snippets:
357 win_size = int(0.002*rate) # 2ms windows
358 threshold = median_std_threshold(data, win_size) # TODO make this a parameter
360 # extract peaks:
361 if 'peak_detection' in return_data:
362 x_peak, x_trough, eod_heights, eod_widths, pd_log_dict = \
363 detect_pulses(i_data, i_rate, threshold,
364 width_fac=width_factor_shape,
365 verbose=verbose, return_data=True)
366 log_dict.update(pd_log_dict)
367 else:
368 x_peak, x_trough, eod_heights, eod_widths = \
369 detect_pulses(i_data, i_rate, threshold,
370 width_fac=width_factor_shape,
371 verbose=verbose, return_data=False)
373 if len(x_peak) > 0:
374 # cluster:
375 min_samples = 5 # TODO make parameter
376 min_samples_frac = 0.05
377 n_gaus_width = 3
378 merge_thresh_width = 0.5
379 n_gaus_height = 10
380 merge_thresh_height = 0.1
381 n_pca = 5
382 shape_eps = 0.3
383 low_freq = 5000
384 min_power_ratio = 0.5
385 pulse_thresh = 0.2
386 clusters, x_merge, c_log_dict = \
387 cluster(i_data, i_rate, x_peak, x_trough, eod_heights, eod_widths,
388 width_factor_shape, width_factor_wave,
389 min_samples=min_samples, min_samples_frac=min_samples_frac,
390 n_gaus_width=n_gaus_width,
391 merge_thresh_width=merge_thresh_width,
392 n_gaus_height=n_gaus_height,
393 merge_thresh_height=merge_thresh_height,
394 n_pca=n_pca, shape_eps=shape_eps,
395 low_freq=low_freq, min_power_ratio=min_power_ratio,
396 pulse_thresh=pulse_thresh, align=align,
397 verbose=verbose, plot_level=plot_level-1,
398 return_data=return_data)
400 # extract mean eods and times:
401 max_samples = int(max_width*i_rate)
402 mean_eods, eod_times, eod_peaktimes, eod_troughtimes, cluster_labels = \
403 extract_mean_waveforms(i_data, i_rate, x_merge, x_peak, x_trough,
404 eod_widths, clusters, max_samples, align)
406 """
407 # TODO: leave that to eodanalysis
408 # determine clipped clusters (save them, but ignore in other steps):
409 clusters, clipped_eods, clipped_times, clipped_peaktimes, clipped_troughtimes = \
410 find_clipped_clusters(clusters, mean_eods, eod_times,
411 eod_peaktimes, eod_troughtimes,
412 cluster_labels, clip_threshold=clip_thresh,
413 verbose=verbose)
414 """
416 # delete moving fish:
417 clusters, mf_log_dict = \
418 delete_moving_fish(clusters, x_merge/i_rate, len(data)/rate,
419 eod_heights, eod_widths/i_rate, i_rate,
420 verbose=verbose,
421 return_data='moving_fish' in return_data)
423 if 'moving_fish' in return_data:
424 log_dict['moving_fish'] = mf_log_dict
426 clusters = remove_sparse_detections(clusters, eod_widths, i_rate,
427 len(data)/rate, verbose=verbose)
429 # extract mean eods
430 mean_eods, eod_times, eod_peaktimes, eod_troughtimes, cluster_labels = \
431 extract_mean_waveforms(i_data, i_rate, x_merge, x_peak, x_trough,
432 eod_widths, clusters, max_samples, align)
434 """
435 mean_eods.extend(clipped_eods)
436 eod_times.extend(clipped_times)
437 eod_peaktimes.extend(clipped_peaktimes)
438 eod_troughtimes.extend(clipped_troughtimes)
439 """
441 if 'all_eod_times' in return_data:
442 log_dict['all_times'] = [x_peak/i_rate, x_trough/i_rate]
443 log_dict['eod_peaktimes'] = eod_peaktimes
444 log_dict['eod_troughtimes'] = eod_troughtimes
446 log_dict.update(c_log_dict)
448 if verbose > 0:
449 print('')
451 return mean_eods, eod_times, eod_peaktimes, log_dict
454def detect_pulses(data, rate, thresh, min_rel_slope_diff=0.25,
455 min_width=0.00005, max_width=0.01, width_fac=5.0,
456 verbose=0, return_data=False):
457 """Detect pulses in data.
459 Parameters
460 ----------
461 data: 1-D array of float
462 The data to be analysed.
463 rate: float
464 Sampling rate of the data.
465 thresh: float
466 Threshold for peak and trough detection via `detect_peaks()`.
467 Must be a positive number that sets the minimum difference
468 between a peak and a trough.
469 min_rel_slope_diff: float
470 Minimum required difference between left and right slope (between
471 peak and troughs) relative to mean slope for deciding which trough
472 to take besed on slope difference.
473 min_width: float
474 Minimum width (peak-trough distance) of pulses in seconds.
475 max_width: float
476 Maximum width (peak-trough distance) of pulses in seconds.
477 width_fac: float
478 Pulses extend plus or minus `width_fac` times their width
479 (distance between peak and assigned trough).
480 Only pulses are returned that can fully be analysed with this width.
481 verbose : int
482 Verbosity level.
483 return_data : bool
484 If `True` data of this function is logged and returned (see
485 extract_pulsefish()).
487 Returns
488 -------
489 peak_indices: array of int
490 Indices of EOD peaks in data.
491 trough_indices: array of int
492 Indices of EOD troughs in data. There is one x_trough for each x_peak.
493 heights: array of float
494 EOD heights for each x_peak.
495 widths: array of int
496 EOD widths for each x_peak (in samples).
497 peak_detection_result : dictionary
498 Key value pairs of logged data.
499 This is only returned if `return_data` is `True`.
501 """
502 peak_detection_result = {}
504 # detect peaks and troughs in the data:
505 peak_indices, trough_indices = detect_peaks(data, thresh)
506 if verbose > 0:
507 print('Peaks/troughs detected in data: %5d %5d'
508 % (len(peak_indices), len(trough_indices)))
509 if return_data:
510 peak_detection_result.update(peaks_1=np.array(peak_indices),
511 troughs_1=np.array(trough_indices))
512 if len(peak_indices) < 2 or \
513 len(trough_indices) < 2 or \
514 len(peak_indices) > len(data)/20:
515 # TODO: if too many peaks increase threshold!
516 if verbose > 0:
517 print('No or too many peaks/troughs detected in data.')
518 if return_data:
519 return np.array([], dtype=int), np.array([], dtype=int), \
520 np.array([]), np.array([], dtype=int), peak_detection_result
521 else:
522 return np.array([], dtype=int), np.array([], dtype=int), \
523 np.array([]), np.array([], dtype=int)
525 # assign troughs to peaks:
526 peak_indices, trough_indices, heights, widths, slopes = \
527 assign_side_peaks(data, peak_indices, trough_indices, min_rel_slope_diff)
528 if verbose > 1:
529 print('Number of peaks after assigning side-peaks: %5d'
530 % (len(peak_indices)))
531 if return_data:
532 peak_detection_result.update(peaks_2=np.array(peak_indices),
533 troughs_2=np.array(trough_indices))
535 # check widths:
536 keep = ((widths>min_width*rate) & (widths<max_width*rate))
537 peak_indices = peak_indices[keep]
538 trough_indices = trough_indices[keep]
539 heights = heights[keep]
540 widths = widths[keep]
541 slopes = slopes[keep]
542 if verbose > 1:
543 print('Number of peaks after checking pulse width: %5d'
544 % (len(peak_indices)))
545 if return_data:
546 peak_detection_result.update(peaks_3=np.array(peak_indices),
547 troughs_3=np.array(trough_indices))
549 # discard connected peaks:
550 same = np.nonzero(trough_indices[:-1] == trough_indices[1:])[0]
551 keep = np.ones(len(trough_indices), dtype=bool)
552 for i in same:
553 # same troughs at trough_indices[i] and trough_indices[i+1]:
554 s = slopes[i:i+2]
555 rel_slopes = np.abs(np.diff(s))[0]/np.mean(s)
556 if rel_slopes > min_rel_slope_diff:
557 keep[i+(s[1]<s[0])] = False
558 else:
559 keep[i+(heights[i+1]<heights[i])] = False
560 peak_indices = peak_indices[keep]
561 trough_indices = trough_indices[keep]
562 heights = heights[keep]
563 widths = widths[keep]
564 if verbose > 1:
565 print('Number of peaks after merging pulses: %5d'
566 % (len(peak_indices)))
567 if return_data:
568 peak_detection_result.update(peaks_4=np.array(peak_indices),
569 troughs_4=np.array(trough_indices))
570 if len(peak_indices) == 0:
571 if verbose > 0:
572 print('No peaks remain as pulse candidates.')
573 if return_data:
574 return np.array([], dtype=int), np.array([], dtype=int), \
575 np.array([]), np.array([], dtype=int), peak_detection_result
576 else:
577 return np.array([], dtype=int), np.array([], dtype=int), \
578 np.array([]), np.array([], dtype=int)
580 # only take those where the maximum cutwidth does not cause issues -
581 # if the width_fac times the width + x is more than length.
582 keep = ((peak_indices - widths > 0) &
583 (peak_indices + widths < len(data)) &
584 (trough_indices - widths > 0) &
585 (trough_indices + widths < len(data)))
587 if verbose > 0:
588 print('Remaining peaks after EOD extraction: %5d'
589 % (np.sum(keep)))
590 print('')
592 if return_data:
593 return peak_indices[keep], trough_indices[keep], \
594 heights[keep], widths[keep], peak_detection_result
595 else:
596 return peak_indices[keep], trough_indices[keep], \
597 heights[keep], widths[keep]
600@jit(nopython=True)
601def assign_side_peaks(data, peak_indices, trough_indices,
602 min_rel_slope_diff=0.25):
603 """Assign to each peak the trough resulting in a pulse with the steepest slope or largest height.
605 The slope between a peak and a trough is computed as the height
606 difference divided by the distance between peak and trough. If the
607 slopes between the left and the right trough differ by less than
608 `min_rel_slope_diff`, then just the heigths between and the two
609 troughs relative to the peak are compared.
611 Was `def detect_eod_peaks(data, main_indices, side_indices,
612 max_width=20, min_width=2, verbose=0)` before.
614 Parameters
615 ----------
616 data: array of float
617 Data in which the events were detected.
618 peak_indices: array of int
619 Indices of the detected peaks in the data time series.
620 trough_indices: array of int
621 Indices of the detected troughs in the data time series.
622 min_rel_slope_diff: float
623 Minimum required difference of left and right slope relative
624 to mean slope.
626 Returns
627 -------
628 peak_indices: array of int
629 Peak indices. Same as input `peak_indices` but potentially shorter
630 by one or two elements.
631 trough_indices: array of int
632 Corresponding trough indices of trough to the left or right
633 of the peaks.
634 heights: array of float
635 Peak heights (distance between peak and corresponding trough amplitude)
636 widths: array of int
637 Peak widths (distance between peak and corresponding trough indices)
638 slopes: array of float
639 Peak slope (height divided by width)
640 """
641 # is a main or side peak first?
642 peak_first = int(peak_indices[0] < trough_indices[0])
643 # is a main or side peak last?
644 peak_last = int(peak_indices[-1] > trough_indices[-1])
645 # ensure all peaks to have side peaks (troughs) at both sides,
646 # i.e. troughs at same index and next index are before and after peak:
647 peak_indices = peak_indices[peak_first:len(peak_indices)-peak_last]
648 y = data[peak_indices]
650 # indices of troughs on the left and right side of main peaks:
651 l_indices = np.arange(len(peak_indices))
652 r_indices = l_indices + 1
654 # indices, distance to peak, height, and slope of left troughs:
655 l_side_indices = trough_indices[l_indices]
656 l_distance = np.abs(peak_indices - l_side_indices)
657 l_height = np.abs(y - data[l_side_indices])
658 l_slope = np.abs(l_height/l_distance)
660 # indices, distance to peak, height, and slope of right troughs:
661 r_side_indices = trough_indices[r_indices]
662 r_distance = np.abs(r_side_indices - peak_indices)
663 r_height = np.abs(y - data[r_side_indices])
664 r_slope = np.abs(r_height/r_distance)
666 # which trough to assign to the peak?
667 # - either the one with the steepest slope, or
668 # - when slopes are similar on both sides
669 # (within `min_rel_slope_diff` difference),
670 # the trough with the maximum height difference to the peak:
671 rel_slopes = np.abs(l_slope-r_slope)/(0.5*(l_slope + r_slope))
672 take_slopes = rel_slopes > min_rel_slope_diff
673 take_left = l_height > r_height
674 take_left[take_slopes] = l_slope[take_slopes] > r_slope[take_slopes]
676 # assign troughs, heights, widths, and slopes:
677 trough_indices = np.where(take_left,
678 trough_indices[:-1], trough_indices[1:])
679 heights = np.where(take_left, l_height, r_height)
680 widths = np.where(take_left, l_distance, r_distance)
681 slopes = np.where(take_left, l_slope, r_slope)
683 return peak_indices, trough_indices, heights, widths, slopes
686def cluster(data, rate, eod_xp, eod_xt, eod_heights, eod_widths,
687 width_factor_shape, width_factor_wave,
688 min_samples=5, min_samples_frac=0.05,
689 n_gaus_width=3, merge_thresh_width=0.6,
690 n_gaus_height=10, merge_thresh_height=0.1,
691 n_pca=5, shape_eps=0.05,
692 low_freq=10000, min_power_ratio=0.5,
693 pulse_thresh=0.2, align='fourier',
694 verbose=0, plot_level=0, return_data=[]):
695 """Cluster EODs.
697 First cluster on EOD widths using a Bayesian Gaussian
698 Mixture (BGM) model, then cluster on EOD heights using a
699 BGM model. Lastly, cluster on EOD waveform with DBSCAN.
700 Clustering on EOD waveform is performed twice, once on
701 peak-centered EODs and once on trough-centered EODs.
702 Non-pulsetype EOD clusters are deleted, and clusters are
703 merged afterwards.
705 Parameters
706 ----------
707 data: array of float
708 Data in which to detect pulse EODs.
709 rate : float
710 Sampling rate of `data`.
711 eod_xp : list of int
712 Location of EOD peaks in indices.
713 eod_xt: list of int
714 Locations of EOD troughs in indices.
715 eod_heights: list of float
716 EOD heights.
717 eod_widths: list of int
718 EOD widths in samples.
719 width_factor_shape : float
720 Multiplier for snippet extraction width. This factor is
721 multiplied with the width between the peak and through of a
722 single EOD.
723 width_factor_wave : float
724 Multiplier for wavefish extraction width.
725 min_samples: int
726 Minimum number of samples required for a valid cluster.
727 min_samples_frac: float
728 Cluster with less samples than this fraction of all the samples are removed.
729 n_gaus_width : int
730 Number of gaussians to use for the clustering based on EOD width.
731 merge_thresh_width : float
732 Threshold for merging clusters that are similar in width.
733 n_gaus_height : int
734 Number of gaussians to use for the clustering based on EOD height.
735 merge_thresh_height : float
736 Threshold for merging clusters that are similar in height.
737 n_pca: int
738 Number of PCs to use for PCA.
739 shape_eps : float
740 Epsilon to use for DBSCAN clustering of EOD shapes.
741 low_freq: float
742 Frequency defining low frequencies in Hertz.
743 min_power_ratio: float
744 Minimum ratio of low-frequency power to total power
745 for mean EOD snippets not being an artefact.
746 pulse_thresh: float
747 Maximum ratio of pulse width relative to inter-pulse interval.
748 align: 'max', 'savgol', 'fourier'
749 Method used to align peak snippets. See extract_snippets() for details.
750 verbose : int
751 Verbosity level.
752 plot_level : int
753 Similar to verbosity levels, but with plots.
754 Only set to > 0 for debugging purposes.
755 return_data : list of str
756 Keys that specify data to be logged. Keys that can be used to log data
757 in this function are: 'all_cluster_steps', 'BGM_width', 'BGM_height',
758 'snippet_clusters', 'artefacts', 'eod_deletion'
759 (see extract_pulsefish()).
761 Returns
762 -------
763 labels : list of int
764 EOD cluster labels based on height and EOD waveform.
765 x_merge : list of int
766 Locations of EODs in clusters.
767 saved_data : dictionary
768 Key value pairs of logged data. Data to be logged is specified
769 by return_data.
771 """
772 saved_data = {}
774 if plot_level > 0 or 'all_cluster_steps' in return_data:
775 all_heightlabels = []
776 all_shapelabels = []
777 all_snippets = []
778 all_features = []
779 all_heights = []
780 all_unique_heightlabels = []
782 all_p_labels = -1*np.ones(len(eod_xp), dtype=int)
783 all_t_labels = -1*np.ones(len(eod_xp), dtype=int)
784 artefact_masks_p = np.ones(len(eod_xp), dtype=bool)
785 artefact_masks_t = np.ones(len(eod_xp), dtype=bool)
787 x_merge = -1*np.ones(len(eod_xp), dtype=int)
789 # keep track of the labels so that no labels are overwritten:
790 max_label_p = 0
791 max_label_t = 0
793 if verbose > 0:
794 print('clusters generated based on EOD width:')
796 # first cluster on width:
797 width_labels, bgm_log_dict = BGM(1000*eod_widths/rate,
798 min_samples=min_samples,
799 min_samples_frac=min_samples_frac,
800 merge_thresh=merge_thresh_width,
801 n_gaus=n_gaus_width, use_log=False,
802 verbose=verbose - 1,
803 return_data='BGM_width' in return_data or plot_level > 0)
805 if 'BGM_width' in return_data:
806 saved_data['BGM_width'] = bgm_log_dict
808 if verbose > 0:
809 # report width clusters:
810 for l in np.unique(width_labels):
811 print(f' {l:2d}: num={len(width_labels[width_labels == l]):5d}, width={np.mean(1000*eod_widths[width_labels == l]/rate):6.3f} +- {np.std(1000*eod_widths[width_labels == l]/rate):6.3f}ms')
813 if plot_level > 0:
814 fig, axs = plt.subplots(3, 1, sharex=True)
815 plot_bgm_cluster(axs, 'width', 'ms', **bgm_log_dict)
816 if plot_level > 1:
817 plt.show()
819 # loop over width clusters:
820 unique_width_labels = np.unique(width_labels[width_labels != -1])
821 for wi, width_label in enumerate(unique_width_labels):
822 # select only features in one width cluster at a time:
823 w_eod_widths = eod_widths[width_labels == width_label]
824 w_eod_heights = eod_heights[width_labels == width_label]
825 w_eod_xp = eod_xp[width_labels == width_label]
826 w_eod_xt = eod_xt[width_labels == width_label]
827 width = int(width_factor_shape*np.median(w_eod_widths))
828 if width > w_eod_xp[0]:
829 width = w_eod_xp[0]
830 if width > w_eod_xt[0]:
831 width = w_eod_xt[0]
832 if width > len(data) - w_eod_xp[-1]:
833 width = len(data) - w_eod_xp[-1]
834 if width > len(data) - w_eod_xt[-1]:
835 width = len(data) - w_eod_xt[-1]
837 wp_labels = -1*np.ones(len(w_eod_xp), dtype=int)
838 wt_labels = -1*np.ones(len(w_eod_xp), dtype=int)
839 wartefact_mask = np.ones(len(w_eod_xp), dtype=int)
841 # determine height labels:
842 raw_p_snippets, p_snippets, p_features, p_bg_ratio = \
843 extract_snippet_features(data, w_eod_xp, w_eod_widths,
844 w_eod_heights, width, n_pca, align)
845 raw_t_snippets, t_snippets, t_features, t_bg_ratio = \
846 extract_snippet_features(data, w_eod_xt, w_eod_widths,
847 w_eod_heights, width, n_pca, align)
849 # TODO: rather keep the height threshold independent of slopes!
850 #median_bg_ratios = np.median(np.min(np.vstack([p_bg_ratio, t_bg_ratio]),
851 # axis=0))
852 #merge_thresh_height = min(merge_thresh_height, median_bg_ratios)
854 if verbose > 0:
855 print(f' clusters generated based on EOD height in width cluster {width_label}:')
857 height_labels, bgm_log_dict = \
858 BGM(w_eod_heights, min_samples=min_samples, min_samples_frac=min_samples_frac,
859 merge_thresh=merge_thresh_height, n_gaus=n_gaus_height,
860 use_log=True, verbose=verbose - 1,
861 return_data='BGM_height' in return_data or plot_level > 0)
863 if 'BGM_height' in return_data:
864 saved_data[f'BGM_height_{wi}'] = bgm_log_dict
866 if verbose > 0:
867 # report height clusters:
868 for l in np.unique(height_labels):
869 print(f' {l:2d}: num={len(height_labels[height_labels == l]):5d}, height={np.mean(w_eod_heights[height_labels == l]):7.4g} +- {np.std(w_eod_heights[height_labels == l]):7.4g}')
871 if plot_level > 0:
872 prop_cycle = plt.rcParams['axes.prop_cycle']
873 colors = prop_cycle.by_key()['color']
874 if np.min(height_labels) == -1:
875 colors = ['black'] + colors
876 fig = plt.figure(layout='constrained')
877 gs = fig.add_gridspec(6, 2)
878 axs = [fig.add_subplot(gs[k*2:(k + 1)*2, 0]) for k in range(3)]
879 axp = fig.add_subplot(gs[0:3, 1])
880 axt = fig.add_subplot(gs[3:6, 1])
881 plot_bgm_cluster(axs, 'height', 'a.u.', **bgm_log_dict)
882 axp.set_title(f'Width cluster {width_label}: raw peak snippets')
883 for l, c in zip(np.unique(height_labels), colors):
884 n = np.sum(height_labels == l)
885 ll = [f'{l} ({n})'] + [None]*(n - 1)
886 axp.plot(raw_p_snippets[height_labels == l, :].T,
887 label=ll, color=c, lw=0.2)
888 axp.set_xlabel('index')
889 axp.legend(title='height label')
890 axt.set_title(f'Width cluster {width_label}: raw trough snippets')
891 for l, c in zip(np.unique(height_labels), colors):
892 n = np.sum(height_labels == l)
893 ll = [f'{l} ({n})'] + [None]*(n - 1)
894 axt.plot(raw_t_snippets[height_labels == l, :].T,
895 label=ll, color=c, lw=0.2)
896 axt.set_xlabel('index')
897 axt.legend(title='height label')
898 plt.show()
900 unique_height_labels = np.unique(height_labels)
901 unique_height_labels = unique_height_labels[unique_height_labels != -1]
903 if plot_level > 0 or 'all_cluster_steps' in return_data:
904 all_heightlabels.append(height_labels)
905 all_heights.append(w_eod_heights)
906 all_unique_heightlabels.append(unique_height_labels)
907 shape_labels = []
908 cfeatures = []
909 csnippets = []
911 for hi, height_label in enumerate(unique_height_labels):
913 h_eod_widths = w_eod_widths[height_labels == height_label]
914 h_eod_heights = w_eod_heights[height_labels == height_label]
915 h_eod_xp = w_eod_xp[height_labels == height_label]
916 h_eod_xt = w_eod_xt[height_labels == height_label]
918 p_feats = p_features[height_labels == height_label]
919 t_feats = t_features[height_labels == height_label]
921 p_labels = cluster_on_shape(p_feats, shape_eps, min_samples,
922 min_samples_frac)
923 t_labels = cluster_on_shape(t_feats, shape_eps, min_samples,
924 min_samples_frac)
926 if verbose > 0:
927 print(f' clusters generated based on EOD shape in width cluster {width_label}, height cluster {height_label}: {len(np.unique(p_labels[p_labels != -1]))} {len(np.unique(t_labels[t_labels != -1]))}')
929 if plot_level > 0:
930 p_snips = p_snippets[height_labels == height_label]
931 t_snips = t_snippets[height_labels == height_label]
932 prop_cycle = plt.rcParams['axes.prop_cycle']
933 color_cycle = prop_cycle.by_key()['color']
934 fig, axs = plt.subplots(2, 2, layout='constrained', sharex='col', sharey='col')
935 axs[0, 0].set_title(f'Width cluster {width_label}, height cluster {height_label}: peak snippets')
936 if np.min(p_labels) == -1:
937 colors = ['black'] + color_cycle
938 else:
939 colors = color_cycle
940 for l, c in zip(np.unique(p_labels), colors):
941 n = np.sum(p_labels == l)
942 ll = [f'{l} ({n})'] + [None]*(n - 1)
943 axs[0, 0].plot(p_feats[p_labels == l, 0], p_feats[p_labels == l, 1], 'o', label=ll, color=c)
944 axs[0, 1].plot(p_snips[p_labels == l, :].T, label=ll, color=c, lw=0.2)
945 axs[0, 1].legend(title='shape label')
946 if np.min(t_labels) == -1:
947 colors = ['black'] + color_cycle
948 else:
949 colors = color_cycle
950 axs[1, 0].set_title(f'Width cluster {width_label}, height cluster {height_label}: trough snippets')
951 for l, c in zip(np.unique(t_labels), colors):
952 n = np.sum(t_labels == l)
953 ll = [f'{l} ({n})'] + [None]*(n - 1)
954 axs[1, 0].plot(t_feats[t_labels == l, 0], t_feats[t_labels == l, 1], 'o', label=ll, color=c)
955 axs[1, 1].plot(t_snips[t_labels == l, :].T, label=ll, color=c, lw=0.2)
956 axs[1, 1].legend(title='shape label')
957 plt.show()
959 if 'snippet_clusters' in return_data:
960 saved_data[f'snippet_clusters_{width_label}_{height_label}_peak'] = dict(
961 raw_snippets=raw_p_snippets[height_labels == height_label],
962 snippets=p_snippets[height_labels == height_label],
963 features=p_features[height_labels == height_label],
964 clusters=p_labels,
965 rate=rate)
966 saved_data[f'snippet_clusters_{width_label}_{height_label}_trough'] = dict(
967 raw_snippets=raw_t_snippets[height_labels == height_label],
968 snippets=t_snippets[height_labels == height_label],
969 features=t_features[height_labels == height_label],
970 clusters=t_labels,
971 rate=rate)
973 if 'all_cluster_steps' in return_data:
974 shape_labels.append([p_labels, t_labels])
975 cfeatures.append([p_features[height_labels == height_label],
976 t_features[height_labels == height_label]])
977 csnippets.append([p_snippets[height_labels == height_label],
978 t_snippets[height_labels == height_label]])
980 p_labels[p_labels == -1] = -max_label_p - 1
981 wp_labels[height_labels == height_label] = p_labels + max_label_p
982 max_label_p = max(np.max(wp_labels), np.max(all_p_labels)) + 1
984 t_labels[t_labels == -1] = -max_label_t - 1
985 wt_labels[height_labels == height_label] = t_labels + max_label_t
986 max_label_t = max(np.max(wt_labels), np.max(all_t_labels)) + 1
988 if verbose > 0:
989 unique_clusters = np.unique(width_label[width_label != -1])
990 print(f' clusters generated based on EOD shape in width cluster {width_label}: {len(unique_clusters):2d}')
991 for l in unique_clusters:
992 print(f' {l:2d}: num={len(wp_labels[wp_labels == l]):5d}')
994 if 'all_cluster_steps' in return_data:
995 all_shapelabels.append(shape_labels)
996 all_snippets.append(csnippets)
997 all_features.append(cfeatures)
999 # remove artefacts here, based on the mean snippets ffts.
1000 artefact_masks_p[width_labels == width_label], sdict = \
1001 remove_artefacts(p_snippets, rate, wp_labels,
1002 threshold=min_power_ratio,
1003 low_freq=low_freq,
1004 verbose=verbose-1,
1005 return_data='artefacts' in return_data)
1006 if 'artefacts' in return_data:
1007 saved_data[f'artefacts_peaks_{width_label}'] = sdict
1008 artefact_masks_t[width_labels == width_label], _ = \
1009 remove_artefacts(t_snippets, rate, wt_labels,
1010 low_freq=low_freq,
1011 threshold=min_power_ratio,
1012 verbose=verbose-1,
1013 return_data='artefacts' in return_data)
1014 if 'artefacts' in return_data:
1015 saved_data[f'artefacts_troughs_{width_label}'] = sdict
1017 # update maxlab so that no clusters are overwritten
1018 all_p_labels[width_labels == width_label] = wp_labels
1019 all_t_labels[width_labels == width_label] = wt_labels
1021 if verbose > 1:
1022 print()
1024 # remove all non-reliable clusters
1025 unreliable_fish_mask_p, saved_data = \
1026 delete_unreliable_fish(all_p_labels, eod_xp, eod_widths,
1027 pulse_thresh=pulse_thresh,
1028 verbose=verbose - 1, sdict=saved_data)
1029 unreliable_fish_mask_t, _ = \
1030 delete_unreliable_fish(all_t_labels, eod_xt, eod_widths,
1031 pulse_thresh=pulse_thresh,
1032 verbose=verbose - 1)
1034 wave_mask_p, sidepeak_mask_p, saved_data = \
1035 delete_wavefish_and_sidepeaks(data, all_p_labels, eod_xp, eod_widths,
1036 width_factor_wave, align=align,
1037 verbose=verbose-1, sdict=saved_data)
1038 wave_mask_t, sidepeak_mask_t, _ = \
1039 delete_wavefish_and_sidepeaks(data, all_t_labels, eod_xt, eod_widths,
1040 width_factor_wave, align=align,
1041 verbose=verbose-1)
1043 og_clusters = [np.copy(all_p_labels), np.copy(all_t_labels)]
1044 og_labels = np.copy(all_p_labels + all_t_labels)
1046 # go through all clusters and masks??
1047 all_p_labels[(artefact_masks_p | unreliable_fish_mask_p | wave_mask_p | sidepeak_mask_p)] = -1
1048 all_t_labels[(artefact_masks_t | unreliable_fish_mask_t | wave_mask_t | sidepeak_mask_t)] = -1
1050 # merge here:
1051 all_clusters, x_merge, mask = merge_clusters(np.copy(all_p_labels),
1052 np.copy(all_t_labels),
1053 eod_xp, eod_xt,
1054 verbose=verbose - 1)
1056 if 'all_cluster_steps' in return_data or plot_level > 0:
1057 all_dmasks = []
1058 all_mmasks = []
1060 discarding_masks = \
1061 np.vstack(((artefact_masks_p | unreliable_fish_mask_p | wave_mask_p | sidepeak_mask_p),
1062 (artefact_masks_t | unreliable_fish_mask_t | wave_mask_t | sidepeak_mask_t)))
1063 merge_mask = mask
1065 # save the masks in the same formats as the snippets
1066 for wi, (width_label, w_shape_label, heightlabels, unique_height_labels) in enumerate(zip(unique_width_labels, all_shapelabels, all_heightlabels, all_unique_heightlabels)):
1067 w_dmasks = discarding_masks[:,width_labels == width_label]
1068 w_mmasks = merge_mask[:,width_labels == width_label]
1070 wd_2 = []
1071 wm_2 = []
1073 for hi, (height_label, h_shape_label) in enumerate(zip(unique_height_labels, w_shape_label)):
1075 h_dmasks = w_dmasks[:,heightlabels == height_label]
1076 h_mmasks = w_mmasks[:,heightlabels == height_label]
1078 wd_2.append(h_dmasks)
1079 wm_2.append(h_mmasks)
1081 all_dmasks.append(wd_2)
1082 all_mmasks.append(wm_2)
1084 if 'all_cluster_steps' in return_data:
1085 saved_data.update(rate=rate,
1086 EOD_widths=[unique_width_labels, eod_widths, width_labels],
1087 EOD_heights=[all_unique_heightlabels, all_heights, all_heightlabels],
1088 EOD_shapes=[all_snippets, all_features, all_shapelabels],
1089 discarding_masks=all_dmasks,
1090 merge_masks=all_mmasks)
1092 if 'masks' in return_data:
1093 saved_data.update(masks=np.vstack(((artefact_masks_p & artefact_masks_t),
1094 (unreliable_fish_mask_p & unreliable_fish_mask_t),
1095 (wave_mask_p & wave_mask_t),
1096 (sidepeak_mask_p & sidepeak_mask_t),
1097 (all_p_labels+all_t_labels))))
1099 if verbose > 0:
1100 unique_clusters = np.unique(all_clusters[all_clusters != -1])
1101 print(f'clusters generated based on height, width and shape: {len(unique_clusters)}')
1102 for l in unique_clusters:
1103 print(f' {l:2d}: num={len(all_clusters[all_clusters == l]):5d}')
1105 return all_clusters, x_merge, saved_data
1108def BGM(x, min_samples=5, min_samples_frac=0.05,
1109 merge_thresh=0.1, n_gaus=5, max_iter=200, n_init=5,
1110 use_log=False, verbose=0, return_data=False):
1111 """Use a Bayesian Gaussian Mixture Model to cluster one-dimensional data.
1113 The data are clustered on their z-scores or on the z-scores of the
1114 log-transformed data if `use_log`is true.
1115 Broad gaussian fits that cover one or more other gaussian fits are
1116 split by their intersections with the other gaussians.
1117 Clusters that are closer than `merge_thresh` on the original scale
1118 of the data (not the z-scores) are merged.
1120 Parameters
1121 ----------
1122 x: 1D numpy array
1123 Features to compute clustering on.
1124 min_samples: int
1125 Minimum number of samples required for a valid cluster.
1126 min_samples_frac: float
1127 Cluster with less samples than this fraction of all the samples are removed.
1128 merge_thresh: float
1129 Ratio for merging nearby gaussians.
1130 n_gaus: int
1131 Maximum number of gaussians to fit on data.
1132 max_iter: int
1133 Maximum number of iterations for gaussian fit.
1134 n_init: int
1135 Number of initializations for the gaussian fit.
1136 use_log: boolean
1137 Set to True to compute the gaussian fit on the logarithm of x.
1138 Can improve clustering on features with nonlinear relationships
1139 such as peak height.
1140 verbose: int
1141 Verbosity level.
1142 return_data: bool
1143 True if additional data should be returned in bgm_dict.
1145 Returns
1146 -------
1147 labels : 1D numpy array
1148 Cluster labels for each sample in x.
1149 bgm_dict : dictionary
1150 Key value pairs of logged data if `return_data` is True,
1151 otherwise empty dictionary.
1153 """
1155 bgm_dict = {}
1157 if len(np.unique(x)) <= n_gaus:
1158 return np.zeros(len(x), dtype=int), bgm_dict
1160 if use_log:
1161 z = stats.zscore(np.log(x)).reshape(-1, 1)
1162 else:
1163 z = stats.zscore(x).reshape(-1, 1)
1164 BGM_model = BayesianGaussianMixture(n_components=n_gaus,
1165 max_iter=max_iter,
1166 n_init=n_init)
1167 labels = BGM_model.fit_predict(z)
1169 if not BGM_model.converged_ and verbose > 0:
1170 print(' !!! Bayesian Gaussian mixture did not converge !!!')
1172 labels_bgm = np.copy(labels)
1174 # separate gaussian clusters that can be split by other clusters:
1175 sidx_x = np.argsort(x)
1176 splits = x[sidx_x][1:][np.diff(labels[sidx_x]) != 0]
1177 labels[:] = 0
1178 for i, split in enumerate(splits):
1179 labels[x >= split] = i + 1
1181 labels_split = np.copy(labels)
1183 # merge gaussian clusters that are closer than merge_thresh:
1184 labels = merge_gaussians(x, labels, min_samples, min_samples_frac,
1185 merge_thresh, verbose=verbose - 1)
1187 # sort model attributes by model.means_:
1188 sidx = np.argsort(BGM_model.means_[:, 0])
1189 means = BGM_model.means_[sidx, 0]
1190 variances = BGM_model.covariances_[sidx, 0, 0]
1191 weights = BGM_model.weights_[sidx]
1193 if return_data:
1194 bgm_dict = dict(x=x,
1195 use_log=use_log,
1196 weights=weights,
1197 means=means,
1198 variances=variances,
1199 labels_bgm=labels_bgm,
1200 labels_split=labels_split,
1201 labels=labels)
1203 return labels, bgm_dict
1206def merge_gaussians(x, labels, min_samples=5, min_samples_frac=0.05,
1207 merge_thresh=0.1, verbose=0):
1208 """ Merge all clusters which have medians which are near.
1210 Only works in 1D.
1212 First, cluster with less than `min_samples` samples are removed -
1213 they get their label set to -1.
1214 Then, clusters where the difference between their medians relative to
1215 the larger of the two medians are smaller than `merge_thresh` are
1216 merged - all their members get the smaller of the two labels assignd.
1218 Parameters
1219 ----------
1220 x : 1D array of int or float
1221 Features used for clustering.
1222 labels : 1D array of int
1223 Labels for each sample in `x`.
1224 min_samples: int
1225 Minimum number of samples required for a valid cluster.
1226 min_samples_frac: float
1227 Cluster with less samples than this fraction of all the samples are removed.
1228 merge_thresh : float
1229 Similarity threshold to merge clusters. The difference between
1230 median values relative to the larger median needs to be smaler
1231 than this threshold for merging.
1232 verbose: int
1233 Verbosity level.
1235 Returns
1236 -------
1237 labels : 1D array of int
1238 Merged labels for each sample in x.
1239 """
1240 # remove small clusters:
1241 u_labels, u_counts = unique_counts(labels[labels != -1])
1242 for l, c in zip(u_labels, u_counts):
1243 if c < min_samples:
1244 labels[labels == l] = -1
1245 if verbose > 0:
1246 print(f' removed cluster {l:2d}: number of samples {c:2d} smaller than {min_samples:2d}')
1247 elif c/len(labels) < min_samples_frac:
1248 labels[labels == l] = -1
1249 if verbose > 0:
1250 print(f' removed cluster {l:2d}: fraction of samples {100*c/len(labels):4.1f}% smaller than {100*min_samples_frac:4.1f}%')
1251 u_labels = np.unique(labels[labels != -1])
1252 if len(u_labels) == 0:
1253 return labels
1254 while len(u_labels) >= 2:
1255 # medians for each label:
1256 x_medians = np.array([np.median(x[labels == l]) for l in u_labels])
1257 sidx = np.argsort(x_medians)
1258 # fill a dict with label mappings:
1259 mapping = {}
1260 for k in range(len(sidx[:-1])):
1261 i0 = sidx[k]
1262 i1 = sidx[k + 1]
1263 label_1 = u_labels[i0]
1264 label_2 = u_labels[i1]
1265 median_1 = x_medians[i0]
1266 median_2 = x_medians[i1]
1267 rel_diff = np.abs(median_2 - median_1)/max(median_1, median_2)
1268 if rel_diff < merge_thresh:
1269 mapping[label_2] = label_1
1270 if verbose > 0:
1271 print(f' merge cluster {label_1:2d} (median={median_1:6.4g}) with cluster {label_2:2d} (median={median_2:6.4g}) because relative difference {rel_diff:4.2f} is smaller than {merge_thresh:4.2f}')
1272 elif verbose > 0:
1273 print(f' keep cluster {label_1:2d} (median={median_1:6.4g}) and cluster {label_2:2d} (median={median_2:6.4g}) because relative difference {rel_diff:4.2f} is larger than {merge_thresh:4.2f}')
1274 # apply mapping:
1275 if len(mapping) > 0:
1276 for key in mapping:
1277 labels[labels == key] = mapping[key]
1278 u_labels = np.unique(labels[labels != -1])
1279 else:
1280 break
1281 return labels
1284def plot_bgm_cluster(axs, xlabel, xunit, x, use_log,
1285 weights, means, variances,
1286 labels_bgm, labels_split, labels):
1287 """Plot histogram of x with BGM clusters.
1288 """
1289 titles = ['BGM', 'split','merge']
1290 all_labels = [labels_bgm, labels_split, labels]
1291 bins = np.linspace(np.min(x), np.max(x), 100)
1292 xx = np.linspace(np.min(x), np.max(x), 500)
1293 for k in range(len(all_labels)):
1294 ax = axs[k]
1295 labs = all_labels[k]
1296 ax.set_title(titles[k])
1297 for l in np.unique(labs[labs != -1]):
1298 xl = x[labs == l]
1299 ax.hist(xl, bins, label=f'{l} ({np.sum(labs == l)})')
1300 if k == len(all_labels) - 1:
1301 ax.set_xlabel(f'{xlabel} [{xunit}]')
1302 ax.set_ylabel('counts')
1303 ax.set_ylim(bottom=0.3)
1304 ax.set_yscale('log')
1305 ax.legend(title=f'{xlabel} labels')
1306 ax = axs[0].twinx()
1307 if use_log:
1308 means_logx = means*np.std(np.log(x)) + np.mean(np.log(x))
1309 stds_logx = np.sqrt(variances)*np.std(np.log(x))
1310 for m, s, w in zip(means_logx, stds_logx, weights):
1311 gg = np.exp(-0.5*((np.log(xx) - m)/s)**2)
1312 ax.plot(xx, w*gg, 'k')
1313 else:
1314 means_x = means*np.std(x) + np.mean(x)
1315 stds_x = np.sqrt(variances)*np.std(x)
1316 for m, s, w in zip(means_x, stds_x, weights):
1317 ax.plot(xx, w*np.exp(-0.5*((xx - m)/s)**2), 'k')
1318 ax.set_ylim(bottom=0)
1321def extract_snippets(data, eod_idx, eod_widths, left, right,
1322 align='fourier'):
1323 """Extract and align snippets from recording data.
1325 Parameters
1326 ----------
1327 data: 1D numpy array of float
1328 Recording data.
1329 eod_idx: 1D array of int
1330 Locations of EODs in data as indices.
1331 eod_widths: 1D array of int
1332 EOD widths (distance between peak and trough) in samples.
1333 left: int
1334 Width to cut out to the left in samples.
1335 right: int
1336 Width to cut out to the right in samples.
1337 align: 'max', 'savgol', 'fourier'
1338 Method used to align peak snippets:
1339 'max': keep detected peak position.
1340 'savgol': use maximum of snippets smoothed with third order
1341 Savitzky-Golay filter.
1342 'fourier': use phase of Fourier coefficient.
1344 Returns
1345 -------
1346 snippets : 2D array of float
1347 Extracted EOD snippets (rows are snippets, columns are time).
1349 Raises
1350 ------
1351 ValueError:
1352 Invalid value for `align`.
1354 """
1355 # extract snippets with corresponding width:
1356 xleft = left + min(10, left//3)
1357 xright = right + min(10, right//3)
1358 snippets = np.zeros((len(eod_idx), xleft + xright))
1359 # first snippets:
1360 k0 = 0
1361 for idx in eod_idx:
1362 if idx >= xleft or len(data) - idx < xright:
1363 break
1364 l = min(idx, xleft)
1365 snippets[k0, xleft - l:] = data[idx - l:idx + xright]
1366 k0 += 1
1367 # last snippets:
1368 k1 = len(eod_idx)
1369 for idx in eod_idx[::-1]:
1370 if k1 < k0 or len(data) - idx >= xright:
1371 break
1372 k1 -= 1
1373 r = min(len(data) - idx, xright)
1374 snippets[k1, :xleft + r] = data[idx - xleft:idx + r]
1375 # middle snippets:
1376 for k, idx in enumerate(eod_idx[k0:k1]):
1377 snippets[k0 + k, :] = data[idx - xleft:idx + xright]
1379 # align snippets on peak:
1380 if align.lower() == 'savgol':
1381 win = 5
1382 filtered = savgol_filter(snippets, win, 2, axis=1)
1383 for k in range(len(snippets)):
1384 peak = np.abs(snippets[k,xleft - win:xleft + win])
1385 ishift = win - np.argmax(peak)
1386 snippets[k] = np.roll(snippets[k], ishift)
1387 elif align.lower() == 'fourier':
1388 # align snippets on phase of first Fourier coefficient:
1389 # (aligning on maximum of higher harmonics is much less robust)
1390 n = snippets.shape[1]
1391 dist = int(np.median(eod_widths))
1392 freq = 1/(2*dist)
1393 idist = min(dist + 1, int(1.2*dist))
1394 ileft = xleft
1395 i0 = ileft - idist
1396 i1 = ileft + idist
1397 if i0 < 0:
1398 ileft += i0
1399 i0 = 0
1400 if i1 > n:
1401 i1 = n
1402 if 2*dist > n:
1403 freq = 1.1/n
1404 coefs = np.zeros(len(snippets), dtype=complex)
1405 for k in range(len(snippets)):
1406 snippet = snippets[k, i0:i1]
1407 m = len(snippet)
1408 coef = fourier_coeffs(snippet, np.arange(m) - ileft, freq, 1)[1]
1409 coefs[k] = coef/np.abs(coef) # only keep the phase
1410 coefs *= np.conjugate(np.mean(coefs)) # deviations from mean phase
1411 for k in range(len(snippets)):
1412 tshift = np.angle(coefs[k])/(2*np.pi*freq)
1413 ishift = int(np.round(tshift))
1414 snippets[k] = np.roll(snippets[k], ishift)
1415 elif align.lower() != 'max':
1416 # max: align snippets on maximum - which they are already!
1417 raise ValueError(f'ERROR in thunderfish.pulses.extract_snippets(): invalid value "{align}" for align parameter - should be "fourier".')
1418 snippets = snippets[:, xleft - left:xleft + right]
1419 return snippets
1422def extract_snippet_features(data, eod_idx, eod_widths, eod_heights,
1423 width, n_pca=5, align='fourier'):
1424 """Extract, align, normalize, snippets from recording data, normalize them, and perform PCA.
1426 Parameters
1427 ----------
1428 data: 1D numpy array of float
1429 Recording data.
1430 eod_idx: 1D array of int
1431 Locations of EODs as indices.
1432 eod_widths: 1D array of int
1433 EOD widths (distance between peak and trough) in samples.
1434 eod_heights: 1D array of float
1435 EOD heights.
1436 width: int
1437 Width to cut out to each side in samples.
1438 n_pca: int
1439 Number of PCs to use for PCA.
1440 align: 'max', 'savgol', 'fourier'
1441 Method used to align peak snippets. See extract_snippets() for details.
1443 Returns
1444 -------
1445 raw_snippets : 2D numpy array (N, EOD_width)
1446 Raw extracted EOD snippets.
1447 snippets : 2D numpy array (N, EOD_width)
1448 Normalized EOD snippets
1449 features : 2D numpy array (N, n_pca)
1450 PC values of EOD snippets
1451 bg_ratio : 1D numpy array (N)
1452 Ratio of the background activity slopes compared to EOD height.
1454 """
1455 raw_snippets = extract_snippets(data, eod_idx, eod_widths, width, width,
1456 align)
1458 # subtract background slope:
1459 snippets, bg_ratio = subtract_slope(np.copy(raw_snippets), eod_heights)
1461 # scale so that their euclidian norm equals one:
1462 snippets = (snippets.T/np.linalg.norm(snippets, axis=1)).T
1464 # compute PCA features for clustering on waveform:
1465 pca = PCA(n_pca)
1466 features = pca.fit_transform(snippets)
1468 return raw_snippets, snippets, features, bg_ratio
1471def cluster_on_shape(features, epsilon=0.05,
1472 min_samples=5, min_samples_frac=0.05):
1473 """Separate EODs by their shape using DBSCAN.
1475 Parameters
1476 ----------
1477 features : 2D numpy array of float (N, n_pc)
1478 PCA features of each EOD in a recording.
1479 epsilon : float
1480 Epsilon to use for DBSCAN clustering.
1481 min_samples: int
1482 Minimum number of samples required for a valid cluster.
1483 min_samples_frac: float
1484 Cluster with less samples than this fraction of all the samples are removed.
1486 Returns
1487 -------
1488 labels : 1D array of int
1489 Merged labels for each sample in x.
1490 """
1491 # minimum samples for core points:
1492 min_smpl = int(len(features)*min_samples_frac)
1493 if min_smpl > min_samples:
1494 min_samples = min_smpl
1496 return DBSCAN(eps=epsilon, min_samples=min_samples).fit(features).labels_
1499def subtract_slope(snippets, heights):
1500 """ Subtract underlying slope from all EOD snippets.
1502 Parameters
1503 ----------
1504 snippets: 2-D numpy array
1505 All EODs in a recorded stacked as snippets.
1506 Shape = (number of EODs, EOD width)
1507 heights: 1D numpy array
1508 EOD heights.
1510 Returns
1511 -------
1512 snippets: 2-D numpy array
1513 EOD snippets with underlying slope subtracted.
1514 bg_ratio : 1-D numpy array
1515 EOD height/background activity height.
1516 """
1518 left_y = snippets[:, 0]
1519 right_y = snippets[:, -1]
1521 try:
1522 slopes = np.linspace(left_y, right_y, snippets.shape[1])
1523 except ValueError:
1524 delta = (right_y - left_y)/snippets.shape[1]
1525 slopes = np.arange(0, snippets.shape[1], dtype=snippets.dtype).reshape((-1,) + (1,) * np.ndim(delta))*delta + left_y
1527 return snippets - slopes.T, np.abs(left_y-right_y)/heights
1530def remove_artefacts(snippets, rate, clusters, low_freq=10000,
1531 threshold=0.5, verbose=0, return_data=False):
1532 """ Create a mask for EOD clusters that result from artefacts, based on power in low frequency spectrum.
1534 Parameters
1535 ----------
1536 snippets: 2D array of float
1537 EOD snippets (each row is an EOD snippet).
1538 rate: float
1539 Sampling rate of snippets.
1540 clusters: array of int
1541 EOD cluster labels.
1542 low_freq: float
1543 Frequency defining low frequencies in Hertz.
1544 threshold: float
1545 Minimum value for sum of low frequency power relative to
1546 sum over full power spectrum that separates artefact from
1547 clean pulsefish clusters.
1548 verbose: int
1549 Verbosity level.
1550 return_data: bool
1551 If True log some data into `adict`.
1553 Returns
1554 -------
1555 mask: numpy array of bool
1556 Set to True for every EOD which is an artefact.
1557 adict: dictionary
1558 Key-value pairs of logged data.
1559 """
1560 adict = {}
1562 mask = np.zeros(clusters.shape, dtype=bool)
1564 if low_freq > rate/2:
1565 return mask, adict
1567 for cluster in np.unique(clusters):
1568 if cluster == -1:
1569 continue
1570 cluster_snippets = snippets[clusters == cluster]
1571 mean_eod = np.mean(cluster_snippets, axis=0)
1572 mean_eod -= np.mean(mean_eod)
1573 mean_eod_power = np.abs(np.fft.rfft(mean_eod))**2
1574 freqs = np.fft.rfftfreq(len(mean_eod), 1/rate)
1575 low_power = np.sum(mean_eod_power[freqs < low_freq])
1576 total_power = np.sum(mean_eod_power)
1577 low_freq_ratio = low_power/total_power
1578 if low_freq_ratio < threshold:
1579 mask[clusters == cluster] = True
1580 if verbose > 0:
1581 print(f' delete cluster {cluster}: low frequency ratio of {low_freq_ratio:4.2f} for frequencies below {low_freq:4.0f}Hz is SMALLER than {threshold:4.2f}')
1582 elif verbose > 1:
1583 print(f' keep cluster {cluster}: low frequency ratio of {low_freq_ratio:4.2f} for frequencies below {low_freq:4.0f}Hz is larger than {threshold:4.2f}')
1585 if return_data:
1586 adict[f'mean_eod_{cluster}'] = mean_eod
1587 adict[f'fft_power_{cluster}'] = mean_eod_power
1588 adict[f'fft_freqs_{cluster}'] = freqs
1589 adict[f'low_frequency_{cluster}'] = low_freq
1590 adict[f'low_frequency_ratio_{cluster}'] = low_freq_ratio
1591 adict[f'mask_{cluster}'] = np.any(mask[clusters == cluster])
1593 return mask, adict
1596def delete_unreliable_fish(clusters, eod_x, eod_widths, pulse_thresh=0.2,
1597 verbose=0, sdict={}):
1598 """ Create a mask for EOD clusters that are either mixed with noise or other fish, or wavefish.
1600 This is the case when the ration between the EOD width and the ISI is too large.
1602 Parameters
1603 ----------
1604 clusters : array of int
1605 Cluster labels.
1606 eod_x : array of int or float
1607 EOD times in samples or seconds.
1608 eod_widths : array of int or float
1609 EOD widths in samples or seconds.
1610 pulse_thresh: float
1611 Maximum ratio of pulse width relative to inter-pulse interval.
1612 verbose : int
1613 Verbosity level.
1614 sdict : dictionary
1615 Dictionary that is used to log data. This is only used if a dictionary
1616 was created by remove_artefacts().
1617 For logging data in noise and wavefish discarding steps,
1618 see remove_artefacts().
1620 Returns
1621 -------
1622 mask : numpy array of booleans
1623 Set to True for every unreliable EOD.
1624 sdict : dictionary
1625 Key value pairs of logged data. Data is only logged if a dictionary
1626 was instantiated by remove_artefacts().
1627 """
1628 mask = np.zeros(clusters.shape, dtype=bool)
1629 for cluster in np.unique(clusters[clusters >= 0]):
1630 # we need more than two samples for diff:
1631 if len(eod_x[clusters == cluster]) < 2:
1632 mask[clusters == cluster] = True
1633 if verbose > 0:
1634 print(f'delete unreliable cluster {cluster}: number of EOD times {len(eod_x[clusters == cluster])} less than 2')
1635 continue
1636 median_width = np.median(eod_widths[clusters == cluster])
1637 ipis = np.diff(eod_x[clusters == cluster])
1638 score = np.max(median_width/ipis)
1639 if score > pulse_thresh:
1640 if verbose > 0:
1641 print(f'delete unreliable cluster {cluster}:maximum width over ipi ratio {score:4.2f} greater than {pulse_thresh:4.2f}')
1642 mask[clusters == cluster] = True
1643 if 'vals_%d' % cluster in sdict:
1644 sdict['vals_%d' % cluster].append(np.median(eod_widths[clusters == cluster])/np.diff(eod_x[clusters == cluster]))
1645 sdict['mask_%d' % cluster].append(any(mask[clusters == cluster]))
1646 return mask, sdict
1649def delete_wavefish_and_sidepeaks(data, clusters, eod_x, eod_widths,
1650 width_fac, max_slope_deviation=0.5,
1651 max_phases=4, align='fourier',
1652 verbose=0, sdict={}):
1653 """ Create a mask for EODs that are likely from wavefish, or sidepeaks of bigger EODs.
1655 Parameters
1656 ----------
1657 data : list of float
1658 Raw recording data.
1659 clusters : list of int
1660 Cluster labels.
1661 eod_x : list of int
1662 Indices of EOD times.
1663 eod_widths : list of int
1664 EOD widths in samples.
1665 width_fac : float
1666 Multiplier for EOD analysis width.
1668 max_slope_deviation: float
1669 Maximum deviation of position of maximum slope in snippets from
1670 center position in multiples of mean width of EOD.
1671 max_phases : int
1672 Maximum number of phases for any EOD.
1673 If the mean EOD has more phases than this, it is not a pulse EOD.
1674 align: 'max', 'savgol', 'fourier'
1675 Method used to align peak snippets. See extract_snippets() for details.
1676 verbose : int
1677 Verbosity level.
1678 sdict : dictionary
1679 Dictionary that is used to log data. This is only used if a dictionary
1680 was created by remove_artefacts().
1681 For logging data in noise and wavefish discarding steps, see remove_artefacts().
1683 Returns
1684 -------
1685 mask_wave: numpy array of booleans
1686 Set to True for every EOD which is a wavefish EOD.
1687 mask_sidepeak: numpy array of booleans
1688 Set to True for every snippet which is centered around a sidepeak of an EOD.
1689 sdict : dictionary
1690 Key value pairs of logged data. Data is only logged if a dictionary
1691 was instantiated by remove_artefacts().
1692 """
1693 mask_wave = np.zeros(clusters.shape, dtype=bool)
1694 mask_sidepeak = np.zeros(clusters.shape, dtype=bool)
1695 sdict = {}
1697 for i, cluster in enumerate(np.unique(clusters[clusters >= 0])):
1698 median_width = np.median(eod_widths[clusters == cluster])
1699 cutwidth = int(median_width*width_fac)
1700 snippets = extract_snippets(data, eod_x[clusters == cluster],
1701 median_width, cutwidth, cutwidth, align)
1703 # extract information on main peaks and troughs:
1704 mean_eod = np.mean(snippets, axis=0)
1705 mean_eod = mean_eod - np.mean(mean_eod)
1707 # detect peaks and troughs on data:
1708 thresh = np.std(mean_eod)
1709 mean_eod[0] += 2*thresh
1710 mean_eod[-1] -= 2*thresh
1711 pk, tr = detect_peaks(mean_eod, thresh)
1712 pk = pk[(pk > 0) & (pk < len(mean_eod))]
1713 tr = tr[(tr > 0) & (tr < len(mean_eod))]
1715 # not enough peaks or troughs:
1716 if len(pk) <= 1 or len(tr) <= 1:
1717 return mask_wave, mask_sidepeak, sdict
1719 idxs = np.sort(np.concatenate((pk, tr)))
1720 slopes = np.abs(np.diff(mean_eod[idxs]))
1721 max_slope = np.argmax(slopes)
1722 # check for side peaks:
1723 centered = np.min(np.abs(idxs[max_slope:max_slope + 2] - len(mean_eod)//2))
1724 if centered > max_slope_deviation*median_width: # TODO: check, factor was probably 0.16
1725 if verbose > 0:
1726 print(f'delete cluster {cluster}, which is a sidepeak')
1727 mask_sidepeak[clusters == cluster] = True
1729 # compute all height differences of peaks and troughs within snippets.
1730 # if they are all similar, it is probably noise or a wavefish.
1731 h_diffs = np.diff(mean_eod[idxs]) # TODO: almost the same as slopes!
1732 w_diffs = np.diff(idxs)
1734 #if np.abs(np.diff(idxs[max_slope:max_slope + 2])) < np.mean(eod_widths[clusters == cluster])*0.5 or len(pk) + len(tr)>max_phases or np.min(w_diffs)>2*cutwidth/width_fac: #or len(h_diffs[np.abs(h_diffs)>0.5*(np.max(mean_eod)-np.min(mean_eod))])>max_phases:
1735 if w_diffs[max_slope] < np.mean(eod_widths[clusters == cluster])/2 or len(idxs) > max_phases or np.min(w_diffs) > 2*cutwidth/width_fac: #or len(h_diffs[np.abs(h_diffs)>0.5*(np.max(mean_eod)-np.min(mean_eod))])>max_phases: # Jan
1736 if verbose > 0:
1737 print(f'delete cluster {cluster}, which is a wavefish')
1738 mask_wave[clusters == cluster] = True
1740 if f'vals_{cluster}' in sdict:
1741 sdict[f'vals_{cluster}'].append([mean_eod, [pk, tr],
1742 idxs[max_slope:max_slope + 2]])
1743 sdict[f'mask_{cluster}'].append(any(mask_wave[clusters == cluster]))
1744 sdict[f'mask_{cluster}'].append(any(mask_sidepeak[clusters == cluster]))
1746 return mask_wave, mask_sidepeak, sdict
1749def merge_clusters(clusters_1, clusters_2, x_1, x_2, verbose=0):
1750 """ Merge clusters resulting from two clustering methods.
1752 This method only works if clustering is performed on the same EODs
1753 with the same ordering, where there is a one to one mapping from
1754 clusters_1 to clusters_2.
1756 Parameters
1757 ----------
1758 clusters_1: list of int
1759 EOD cluster labels for cluster method 1.
1760 clusters_2: list of int
1761 EOD cluster labels for cluster method 2.
1762 x_1: list of int
1763 Indices of EODs for cluster method 1 (clusters_1).
1764 x_2: list of int
1765 Indices of EODs for cluster method 2 (clusters_2).
1766 verbose : int
1767 Verbosity level.
1769 Returns
1770 -------
1771 clusters : list of int
1772 Merged clusters.
1773 x_merged : list of int
1774 Merged cluster indices.
1775 mask : 2d numpy array of int (N, 2)
1776 Mask for clusters that are selected from clusters_1 (mask[:,0]) and
1777 from clusters_2 (mask[:,1]).
1778 """
1779 if verbose > 0:
1780 print('\nmerge cluster:')
1782 # these arrays become 1 for each EOD that is chosen from that array
1783 c1_keep = np.zeros(len(clusters_1), dtype=int)
1784 c2_keep = np.zeros(len(clusters_2), dtype=int)
1786 # add n to one of the cluster lists to avoid overlap
1787 ovl = np.max(clusters_1) + 1
1788 clusters_2[clusters_2 != -1] = clusters_2[clusters_2 !=- 1] + ovl
1790 remove_clusters = [[]]
1791 keep_clusters = []
1792 og_clusters = [np.copy(clusters_1), np.copy(clusters_2)]
1794 # loop until done
1795 while True:
1797 # compute unique clusters and cluster sizes
1798 # of cluster that have not been iterated over:
1799 c1_labels, c1_size = unique_counts(clusters_1[(clusters_1 != -1) & (c1_keep == 0)])
1800 c2_labels, c2_size = unique_counts(clusters_2[(clusters_2 != -1) & (c2_keep == 0)])
1802 # if all clusters are done, break from loop:
1803 if len(c1_size) == 0 and len(c2_size) == 0:
1804 break
1806 # if the biggest cluster is in c_p, keep this one and discard all clusters
1807 # on the same indices in c_t:
1808 elif np.argmax([np.max(np.append(c1_size, 0)), np.max(np.append(c2_size, 0))]) == 0:
1810 # remove all the mappings from the other indices
1811 cluster_mappings, _ = unique_counts(clusters_2[clusters_1 == c1_labels[np.argmax(c1_size)]])
1813 clusters_2[np.isin(clusters_2, cluster_mappings)] = -1
1815 c1_keep[clusters_1 == c1_labels[np.argmax(c1_size)]] = 1
1817 remove_clusters.append(cluster_mappings)
1818 keep_clusters.append(c1_labels[np.argmax(c1_size)])
1820 if verbose > 0:
1821 print(f' group 1: keep cluster {c1_labels[np.argmax(c1_size)]:2d}, group 2: delete clusters {cluster_mappings[cluster_mappings != -1] - ovl}')
1823 # if the biggest cluster is in c_t, keep this one and discard all mappings in c_p
1824 elif np.argmax([np.max(np.append(c1_size, 0)), np.max(np.append(c2_size, 0))]) == 1:
1826 # remove all the mappings from the other indices
1827 cluster_mappings, _ = unique_counts(clusters_1[clusters_2 == c2_labels[np.argmax(c2_size)]])
1829 clusters_1[np.isin(clusters_1, cluster_mappings)] = -1
1831 c2_keep[clusters_2 == c2_labels[np.argmax(c2_size)]] = 1
1833 remove_clusters.append(cluster_mappings)
1834 keep_clusters.append(c2_labels[np.argmax(c2_size)])
1836 if verbose > 0:
1837 print(f' group 2: keep cluster {c2_labels[np.argmax(c2_size)] - ovl:2d}, group 1: delete clusters {cluster_mappings[cluster_mappings!=-1]}')
1839 # combine results
1840 clusters = (clusters_1 + 1)*c1_keep + (clusters_2 + 1)*c2_keep - 1
1841 x_merged = (x_1)*c1_keep + (x_2)*c2_keep
1843 return clusters.astype(int), x_merged, np.vstack([c1_keep, c2_keep])
1846def extract_mean_waveforms(data, rate, pos_inx, peak_inx, trough_inx,
1847 widths, labels, max_width, align='fourier'):
1848 """ Extract time points and mean waveforms for each cluster.
1850 Parameters
1851 ----------
1852 data: 1-D array of float
1853 Timeseries with the detected pulses.
1854 rate: float
1855 Sampling rate of time series.
1856 pos_inx: list of int
1857 Locations of pulses in samples.
1858 peak_inx : list of int
1859 Locations of peaks in samples.
1860 trough_inx : list of int
1861 Locations of troughs in samples.
1862 widths: list of int
1863 Widths in samples.
1864 labels: list of int
1865 Cluster labels.
1866 max_width : int
1867 Maximum width of returned mean waveform in samples.
1868 align: 'max', 'savgol', 'fourier'
1869 Method used to align peak snippets. See extract_snippets() for details.
1871 Returns
1872 -------
1873 mean_waveforms: list of 2D array of float
1874 The average waveform for each detected cluster.
1875 First column is time in seconds, second column mean waveform,
1876 third column the corresponding standard error of the mean.
1877 pos_times: list of 1D array of float
1878 For each cluster the time points of waveforms in seconds.
1879 peak__times: list of 1D array of float
1880 For each cluster the time points of waveform peaks in seconds.
1881 trough_times: list of 1D array of float
1882 For each detected fish the time points of waveform troughs in seconds.
1883 ulabels: list of int
1884 Label for each cluster.
1885 """
1886 waveforms = []
1887 pos_times = []
1888 peak_times = []
1889 trough_times = []
1890 heights = []
1891 ulabels = []
1893 for l in np.unique(labels):
1894 if l == -1:
1895 continue
1896 w = widths[labels == l]
1897 median_ipi = int(np.median(np.diff(pos_inx[labels == l])))
1898 left = median_ipi//3
1899 right = median_ipi - left
1900 if left + right > max_width:
1901 left = max_width//3
1902 right = max_width - left
1903 snippets = extract_snippets(data, pos_inx[labels == l], w,
1904 left, right, align)
1906 mean = np.mean(snippets, axis=0)
1907 sem = np.std(snippets, axis=0)/np.sqrt(len(snippets))
1908 height = np.max(mean) - np.min(mean)
1909 time = np.arange(len(mean))/rate - left/rate
1910 waveform = np.column_stack([time, mean, sem])
1912 waveforms.append(waveform)
1913 pos_times.append(pos_inx[labels == l]/rate)
1914 heights.append(height)
1915 peak_times.append(peak_inx[labels == l]/rate)
1916 trough_times.append(trough_inx[labels == l]/rate)
1917 ulabels.append(l)
1919 sidx = list(reversed(np.argsort(heights)))
1920 return [waveforms[i] for i in sidx], [pos_times[i] for i in sidx], [peak_times[i] for i in sidx], [trough_times[i] for i in sidx], [ulabels[i] for i in sidx]
1923"""
1924TODO: this is done in eodanalysis
1925def find_clipped_clusters(clusters, mean_eods, eod_times,
1926 eod_peaktimes, eod_troughtimes,
1927 cluster_labels, clip_threshold=0.9, verbose=0):
1928 "" Detect EODs that are clipped and set all clusterlabels of these clipped EODs to -1.
1930 Also return the mean EODs and timepoints of these clipped EODs.
1932 Parameters
1933 ----------
1934 clusters: array of int
1935 Cluster labels for each EOD in a recording.
1936 mean_eods: list of numpy arrays
1937 Mean EOD waveform for each cluster.
1938 eod_times: list of numpy arrays
1939 EOD timepoints for each EOD cluster.
1940 eod_peaktimes
1941 EOD peaktimes for each EOD cluster.
1942 eod_troughtimes
1943 EOD troughtimes for each EOD cluster.
1944 cluster_labels: numpy array
1945 Unique EOD clusterlabels.
1946 clip_threshold: float
1947 Threshold for detecting clipped EODs.
1949 verbose: int
1950 Verbosity level.
1952 Returns
1953 -------
1954 clusters : array of int
1955 Cluster labels for each EOD in the recording, where clipped EODs have been set to -1.
1956 clipped_eods : list of numpy arrays
1957 Mean EOD waveforms for each clipped EOD cluster.
1958 clipped_times : list of numpy arrays
1959 EOD timepoints for each clipped EOD cluster.
1960 clipped_peaktimes : list of numpy arrays
1961 EOD peaktimes for each clipped EOD cluster.
1962 clipped_troughtimes : list of numpy arrays
1963 EOD troughtimes for each clipped EOD cluster.
1964 ""
1965 clipped_eods = []
1966 clipped_times = []
1967 clipped_peaktimes = []
1968 clipped_troughtimes = []
1969 clipped_labels = []
1971 for mean_eod, eod_time, eod_peaktime, eod_troughtime,label in zip(mean_eods, eod_times, eod_peaktimes, eod_troughtimes, cluster_labels):
1972 print(mean_eod.shape)
1974 if (np.count_nonzero(mean_eod[:, 1] >= clip_threshold) > 0 or (np.count_nonzero(mean_eod[:, 1] <= -clip_threshold) > 0):
1975 clipped_eods.append(mean_eod)
1976 clipped_times.append(eod_time)
1977 clipped_peaktimes.append(eod_peaktime)
1978 clipped_troughtimes.append(eod_troughtime)
1979 clipped_labels.append(label)
1980 if verbose > 0:
1981 print('clipped pulsefish')
1983 clusters[np.isin(clusters, clipped_labels)] = -1
1985 return clusters, clipped_eods, clipped_times, clipped_peaktimes, clipped_troughtimes
1986"""
1988def delete_moving_fish(clusters, eod_t, T, eod_heights, eod_widths,
1989 rate, min_dt=0.25, stepsize=0.05,
1990 sliding_window_factor=2000, verbose=0, return_data=[]):
1991 """
1992 Use a sliding window to detect the minimum number of fish detected simultaneously,
1993 then delete all other EOD clusters.
1995 Do this only for EODs within the same width clusters, as a
1996 moving fish will preserve its EOD width.
1998 Parameters
1999 ----------
2000 clusters: list of int
2001 EOD cluster labels.
2002 eod_t: list of float
2003 Timepoints of the EODs (in seconds).
2004 T: float
2005 Length of recording (in seconds).
2006 eod_heights: list of float
2007 EOD amplitudes.
2008 eod_widths: list of float
2009 EOD widths (in seconds).
2010 rate: float
2011 Recording data sampling rate.
2013 min_dt : float
2014 Minimum sliding window size (in seconds).
2015 stepsize : float
2016 Sliding window stepsize (in seconds).
2017 sliding_window_factor : float
2018 Multiplier for sliding window width,
2019 where the sliding window width = median(EOD_width)*sliding_window_factor.
2020 verbose : int
2021 Verbosity level.
2022 return_data : bool
2023 If True fill mf_dict with data.
2025 Returns
2026 -------
2027 clusters : list of int
2028 Cluster labels, where deleted clusters have been set to -1.
2029 mf_dict : dictionary
2030 Key value pairs of logged data. Data to be logged is specified by return_data.
2031 """
2032 mf_dict = {}
2034 if len(np.unique(clusters[clusters != -1])) == 0:
2035 return clusters, {}
2037 all_keep_clusters = []
2038 width_classes = merge_gaussians(eod_widths, np.copy(clusters), 0.75)
2040 ev_num = 0
2041 for iw, w in enumerate(np.unique(width_classes[clusters >= 0])):
2042 # initialize variables
2043 min_clusters = 10000
2044 average_height = 0
2045 sparse_clusters = 10000
2046 keep_clusters = []
2048 dt = max(min_dt, np.median(eod_widths[width_classes == w])*sliding_window_factor)
2049 window_start = 0
2050 window_end = dt
2052 wclusters = clusters[width_classes == w]
2053 weod_t = eod_t[width_classes == w]
2054 weod_heights = eod_heights[width_classes == w]
2055 weod_widths = eod_widths[width_classes == w]
2057 if verbose > 1:
2058 print(f'moving_fish: sliding window dt = {dt:.3f}s')
2060 x = np.arange(0, T - dt + stepsize, stepsize)
2061 y = np.ones(len(x), dtype=int)
2063 if len(x) == 0:
2064 keep_clusters.append(np.unique(wclusters))
2065 else:
2066 # make W dependent on width??
2067 ignore_steps = np.zeros(len(x), dtype=int)
2069 for i, t in enumerate(x):
2070 current_clusters = wclusters[(weod_t>=t)&(weod_t<t+dt)&(wclusters!=-1)]
2071 if len(np.unique(current_clusters)) == 0:
2072 ignore_steps[i-int(dt/stepsize):i+int(dt/stepsize)] = 1
2073 if verbose > 0:
2074 print(f'moving fish: no pulsefish in recording at T={t:.2f}:{t + dt:.2f}')
2076 running_sum = np.ones(len(x), dtype=int)
2077 ulabs = np.unique(wclusters[wclusters>=0])
2079 # sliding window
2080 for j, (t, ignore_step) in enumerate(zip(x, ignore_steps)):
2081 current_clusters = wclusters[(weod_t >= t) & (weod_t < t+dt) & (wclusters != -1)]
2082 current_widths = weod_widths[(weod_t >= t) & (weod_t < t+dt) & (wclusters != -1)]
2084 unique_clusters = np.unique(current_clusters)
2085 y[j] = len(unique_clusters)
2087 if (len(unique_clusters) <= min_clusters) and \
2088 (ignore_step == 0) and \
2089 (len(unique_clusters != 1)):
2091 current_labels = np.isin(wclusters, unique_clusters)
2092 current_height = np.mean(weod_heights[current_labels])
2094 # compute nr of clusters that are too sparse:
2095 clusters_after_deletion = np.unique(remove_sparse_detections(np.copy(clusters[np.isin(clusters, unique_clusters)]), rate*eod_widths[np.isin(clusters, unique_clusters)], rate, T))
2096 current_sparse_clusters = len(unique_clusters) - len(clusters_after_deletion[clusters_after_deletion != -1])
2098 if current_sparse_clusters <= sparse_clusters and \
2099 ((current_sparse_clusters < sparse_clusters) or
2100 (current_height > average_height) or
2101 (len(unique_clusters) < min_clusters)):
2102 keep_clusters = [unique_clusters]
2103 min_clusters = len(unique_clusters)
2104 average_height = current_height
2105 window_end = t+dt
2106 sparse_clusters = current_sparse_clusters
2108 all_keep_clusters.extend(keep_clusters)
2110 if return_data:
2111 if 'w' in mf_dict:
2112 mf_dict['w'].append(np.median(eod_widths[width_classes == w]))
2113 mf_dict['T'] = T
2114 mf_dict['dt'].append(dt)
2115 mf_dict['clusters'].append(wclusters)
2116 mf_dict['t'].append(weod_t)
2117 mf_dict['fishcount'].append([x+0.5*(x[1]-x[0]), y])
2118 mf_dict['ignore_steps'].append(ignore_steps)
2119 else:
2120 mf_dict['w'] = [np.median(eod_widths[width_classes == w])]
2121 mf_dict['T'] = T
2122 mf_dict['dt'] = [dt]
2123 mf_dict['clusters'] = [wclusters]
2124 mf_dict['t'] = [weod_t]
2125 mf_dict['fishcount'] = [[x+0.5*(x[1]-x[0]), y]]
2126 mf_dict['ignore_steps'] = [ignore_steps]
2128 # delete all clusters that are not selected:
2129 clusters[np.invert(np.isin(clusters, np.concatenate(all_keep_clusters)))] = -1
2130 if verbose > 0:
2131 print(f'moving fish: estimated number of pulsefish: {len(np.unique(clusters[clusters != -1]))}')
2133 return clusters, mf_dict
2136def remove_sparse_detections(clusters, eod_widths, rate, T,
2137 min_density=0.0005, verbose=0):
2138 """ Remove all EOD clusters that are too sparse
2140 Parameters
2141 ----------
2142 clusters : list of int
2143 Cluster labels.
2144 eod_widths : list of int
2145 Cluster widths in samples.
2146 rate : float
2147 Sampling rate.
2148 T : float
2149 Lenght of recording in seconds.
2150 min_density : float
2151 Minimum density for realistic EOD detections.
2152 verbose : int
2153 Verbosity level.
2155 Returns
2156 -------
2157 clusters : list of int
2158 Cluster labels, where sparse clusters have been set to -1.
2159 """
2160 for c in np.unique(clusters):
2161 if c == -1:
2162 continue
2163 n = len(clusters[clusters == c])
2164 w = np.median(eod_widths[clusters == c])/rate
2165 if n*w < T*min_density:
2166 if verbose > 0:
2167 print('removed cluster {c:2d}: median coverage of {n*w:.2f}s is less than {100*min_density:.1f}% of data')
2168 clusters[clusters == c] = -1
2169 return clusters
2172def add_extract_pulsefish_config(cfg, frate=0.5e6, width_factor_shape=3,
2173 width_factor_wave=8, max_width=0.03,
2174 align='fourier'):
2175 """Add all parameters needed for extract_pulsefish() to a configuration.
2177 Parameters
2178 ----------
2179 cfg: ConfigFile
2180 The configuration.
2182 See `extract_pulsefish()` for details on the remaining arguments.
2183 """
2184 cfg.add_section('Pulse-type fish extraction:')
2185 cfg.add('upSamplePulse', frate, 'Hz', 'Sampling rate to which EOD snippets are upsampled.')
2186 cfg.add('pulseShapeFactor', width_factor_shape, '', 'Size of snippets for shape clustering as a multiple of peak-trough distance.')
2187 cfg.add('pulseWaveFactor', width_factor_wave, '', 'Size of snippets for removing wavefish as a multiple of peak-trough distance.')
2188 cfg.add('pulseMaximumWidth', max_width, 's', 'Maximum width of final snippets for further analysis and display.')
2189 cfg.add('alignPulses', align, '', 'Method used to align EOD snippets: "fourier", "savgol", or "max".')
2192def extract_pulsefish_args(cfg):
2193 """Retrieve parameters for `extract_pulsefish()` from configuration.
2195 The return value can then be passed as key-word arguments to this
2196 function.
2198 Parameters
2199 ----------
2200 cfg: ConfigFile
2201 The configuration.
2203 Returns
2204 -------
2205 a: dict
2206 Dictionary with names of arguments of the `extract_pulsefish()` function
2207 and their values as supplied by `cfg`.
2208 """
2209 a = cfg.map(frate='upSamplePulse',
2210 width_factor_shape='pulseShapeFactor',
2211 width_factor_wave='pulseWaveFactor',
2212 max_width='pulseMaximumWidth',
2213 align='alignPulses')
2214 return a