diff --git a/do_make_pancake.py b/do_make_pancake.py new file mode 100644 index 0000000..04f8ba4 --- /dev/null +++ b/do_make_pancake.py @@ -0,0 +1,152 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +from IPython import embed +import helper_functions as hf +import os +import datetime +import time +import pandas as pd + +if __name__ == '__main__': + + ################################################################################################################### + # load data + ################################################################################################################### + # load all the data of one day + start = time.time() + + filename = sorted(os.listdir('../../../data/'))[6] + + ident = np.load('../../../data/'+filename+'/all_ident_v.npy', allow_pickle=True) + power = np.load('../../../data/'+filename+'/all_sign_v.npy', allow_pickle=True) + freq = np.load('../../../data/'+filename+'/all_fund_v.npy', allow_pickle=True) + timeidx = np.load('../../../data/'+filename+'/all_idx_v.npy', allow_pickle=True) + realtime = np.load('../../../data/'+filename+'/all_times.npy', allow_pickle=True) + temp = pd.read_csv('../../../data/' + filename + '/temperatures.csv', sep=';') + temp_l = np.array(temp.values.tolist()) + + aifl = np.load('../data/'+filename+'/aifl2.npy', allow_pickle=True) + end = time.time() + print(end - start) + ################################################################################################################### + # parameter and variables + ################################################################################################################### + # lists + fish_in_aifl = list(np.unique(np.where(~np.isnan(aifl[:, :, 0]))[1])) + + # variables and empty lists + sampling_rate = 1 / np.diff(realtime)[0] # in sec + # ipp + power_means = [] + all_ipp = [] + all_xtickses = [] + all_run_mean = [] + all_run_std = [] + all_Ctime_gesamt = [] + all_max_ch = [] + all_hms = [] + thresholds = np.full(3, np.nan) + + ################################################################################################################### + # analysis + ################################################################################################################### + + for fish_number in fish_in_aifl: + + ############################################################################################################### + # power cube: 1 dim: different channel for different traces + # 2 dim: time dimension + # 3 dim: again channel but in this direction the power of the traces of one channel are written into + power_cube = np.full([16, len(realtime), 16], np.nan) + + for channel in range(16): + fish_IDs = aifl[channel, fish_number, ~np.isnan(aifl[channel, fish_number])] + for len_idx1 in range(len(fish_IDs)): + ID = fish_IDs[len_idx1] + t = timeidx[channel][ident[channel] == ID] + p = power[channel][ident[channel] == ID] + power_cube[channel, t] = p + + ############################################################################################################### + # interpolated power pancake = ipp, heatmap of the trajectory of the fish over the time (2 dimension) and + # channels (1. dimension) + power_pancake = np.nanmax(power_cube, axis=0) + Ctime = realtime[~np.isnan(power_pancake[:, 0])] + all_Ctime = realtime[int(np.where(realtime == Ctime[0])[0]):int(np.where(realtime == Ctime[-1])[0] + 1)] + Cpancake = power_pancake[~np.isnan(power_pancake[:, 0])] + + try: + ipp = np.array(list(map(lambda x: np.interp(all_Ctime, Ctime, x), Cpancake.T))).T + except: + continue + all_ipp.append(ipp) + all_Ctime_gesamt.append(all_Ctime) + + ############################################################################################################### + # trajectories (channel) of the fish over time by using the maximum power on each time point + max_ch = np.argmax(ipp, axis=1) + all_max_ch.append(max_ch) + + ############################################################################################################### + # power means of each fish + power_means.append(np.mean(ipp[range(len(max_ch)), max_ch])) + + ############################################################################################################### + # all x ticks of the fish in datetime format + datetime_v = [] + # hour minutes seconds vector + hms = [] + dt = datetime.datetime.strptime(filename[-5:], '%H_%M') + for i in list(np.arange(int(np.where(realtime == Ctime[0])[0]), int(np.where(realtime == Ctime[-1])[0] + 1))): + current_time = dt + datetime.timedelta(seconds=realtime[i]) # converts the time point in datetime points + datetime_v.append(current_time) + hms.append(current_time.hour * 60 * 60 + current_time.minute * 60 + current_time.second + + float('0.' + str(current_time.microsecond))) # converts the time points into a second format + + x_tickses = mdates.date2num(datetime_v) + all_xtickses.append(x_tickses) + all_hms.append(np.array(hms)) + + ############################################################################################################### + # running mean and std + # run_mean, run_std = hf.running_3binsizes(ipp, sampling_rate) + # TODO: not absolute max, instead average of the three strongest channel + # bin_size = int(np.floor((Ctime[-1]-Ctime[0])/50)) + # + # max_ch_mean = np.full([len(max_ch)], np.nan) + # for i in range(int(len(max_ch)-np.floor(bin_size))): + # max_ch_mean[int(i + np.floor(bin_size / 2))] = np.mean(max_ch[i:bin_size + i]) + + # all_run_mean.append(run_mean) + # all_run_std.append(run_std) + + ############################################################################################################### + # threshold for the activity -- OUT DATED!!!!!! + # rs = np.array(all_run_std) + # + # for idx in range(3): + # rs1 = np.hstack(rs[:, idx]) + # rs1 = rs1[~np.isnan(rs1)] + # rs1 = rs1[rs1>0] + # print(np.percentile(rs1, [5, 95, 50])) + # thresholds[idx] = np.percentile(rs1, 50) + ################################################################################################################### + # save + ################################################################################################################### + if filename not in os.listdir('../data/'): + os.mkdir('../data/'+filename) + + np.save('../data/' + filename + '/all_Ctime_v.npy', all_Ctime_gesamt) + np.save('../data/' + filename + '/all_max_ch.npy', all_max_ch) + # np.save('../data/' + filename + '/all_run_mean.npy', all_run_mean) + # np.save('../data/' + filename + '/all_run_std.npy', all_run_std) + # np.save('../data/' + filename + '/thresholds.npy', thresholds) + np.save('../data/' + filename + '/all_xtickses.npy', all_xtickses) + np.save('../data/' + filename + '/all_ipp.npy', all_ipp) + np.save('../data/' + filename + '/power_means.npy', power_means) + np.save('../data/' + filename + '/all_hms.npy', all_hms) + + + embed() + quit() diff --git a/do_match_eod_waveform.py b/do_match_eod_waveform.py new file mode 100644 index 0000000..791f0a1 --- /dev/null +++ b/do_match_eod_waveform.py @@ -0,0 +1,143 @@ +import sys +import os +import numpy as np +from tqdm import tqdm +import matplotlib.pyplot as plt +from thunderfish.dataloader import open_data +from thunderfish.eodanalysis import eod_waveform +from IPython import embed +import matplotlib.gridspec as gridspec + + +def plot_mean_eod(EOD, fish_number, ax): + + ax.plot(EOD.T[0], EOD.T[1], color='firebrick', lw=2) + ax.fill_between(EOD.T[0], EOD.T[1] + EOD.T[2], + EOD.T[1] - EOD.T[2], + color='firebrick', alpha=0.7) + + +def calc_mean_eod(t0, f, data, dt = 10): + channel_list = np.arange(data.channels) + samplerate = data.samplerate + + start_i = t0 * samplerate + end_i = t0 * samplerate + dt * samplerate + 1 + t = np.arange(0, dt, 1 / f) + + mean_EODs = [] + for c in channel_list: + mean_eod, eod_times = eod_waveform(data[start_i:end_i, c], samplerate, t) + # TODO: unfilter function + mean_EODs.append(mean_eod) + + max_size = list(map(lambda x: np.max(x.T[1]) - np.min(x.T[1]), mean_EODs)) + EOD = mean_EODs[np.argmax(max_size)] + + return EOD + + +def main(folder, filename): + # folder = path_to_files + data = open_data(os.path.join(folder, 'traces-grid1.raw'), -1, 60.0, 10.0) + + power_means = np.load('../data/' + filename + '/power_means.npy', allow_pickle=True) + all_q10 = np.load('../data/' + filename + '/fish_freq_q10.npy', allow_pickle=True) + all_t = np.load('../data/' + filename + '/eod_times_new_new.npy', allow_pickle=True) + all_f = np.load('../data/' + filename + '/eod_freq_new_new.npy', allow_pickle=True) + + # all_fish = np.load('../data/' + filename + '/eod_fishnumber.npy', allow_pickle=True) + + name = [] + for fish_number in range(len(power_means)): + # print('fish:' + fish_number) + if 400 <= all_q10[fish_number, 2] <= 600: + + # print(len(all_t[fish_number])) + fig = plt.figure(constrained_layout=True) + gs_pic = np.sqrt(len(all_t[fish_number])) + gs = gridspec.GridSpec(int(np.ceil(gs_pic)), int(np.ceil(gs_pic)), figure=fig) + # max_std = [] + # EODs = [] + c = 0 + for t, f in zip(all_t[fish_number], all_f[fish_number]): + EOD = calc_mean_eod(t, f, data) + + ax = fig.add_subplot(gs[c]) + ax.plot(EOD.T[0], EOD.T[1], color='firebrick', lw=2) + ax.fill_between(EOD.T[0], EOD.T[1] + EOD.T[2], + EOD.T[1] - EOD.T[2], + color='firebrick', alpha=0.7) + ax.set_xticks([]) + ax.set_yticks([]) + fig.add_subplot(ax) + c+=1 + fig.suptitle(all_q10[fish_number, 2]) + plt.show() + + in_str = input('Which fish? 1 = Eigen, 2 = Apt') + if in_str == '1': + name.append('Eigenmannia') + elif in_str == '2': + name.append('Apteronotus') + else: + name.append('unknown') + elif all_q10[fish_number, 2] < 500: + name.append('Eigenmannia') + elif all_q10[fish_number, 2] > 600: + name.append('Apteronotus') + + return name + + +def plot_eod_for_material_and_methods(folder, filename): + data = open_data(os.path.join(folder, 'traces-grid1.raw'), -1, 60.0, 10.0) + + power_means = np.load('../data/' + filename + '/power_means.npy', allow_pickle=True) + all_q10 = np.load('../data/' + filename + '/fish_freq_q10.npy', allow_pickle=True) + all_t = np.load('../data/' + filename + '/eod_times_new_new.npy', allow_pickle=True) + all_f = np.load('../data/' + filename + '/eod_freq_new_new.npy', allow_pickle=True) + + # all_fish = np.load('../data/' + filename + '/eod_fishnumber.npy', allow_pickle=True) + + name = [] + for fish_number in range(len(power_means)): + # print('fish:' + fish_number) + if 400 <= all_q10[fish_number, 2] <= 600: + + # print(len(all_t[fish_number])) + fig = plt.figure(constrained_layout=True) + gs_pic = np.sqrt(len(all_t[fish_number])) + gs = gridspec.GridSpec(int(np.ceil(gs_pic)), int(np.ceil(gs_pic)), figure=fig) + # max_std = [] + # EODs = [] + c = 0 + for t, f in zip(all_t[fish_number], all_f[fish_number]): + EOD = calc_mean_eod(t, f, data) + + ax = fig.add_subplot(gs[c]) + ax.plot(EOD.T[0], EOD.T[1], color='firebrick', lw=2) + ax.fill_between(EOD.T[0], EOD.T[1] + EOD.T[2], + EOD.T[1] - EOD.T[2], + color='firebrick', alpha=0.7) + ax.set_xticks([]) + ax.set_yticks([]) + fig.add_subplot(ax) + c += 1 + fig.suptitle(all_q10[fish_number, 2]) + plt.show() + + +if __name__ == '__main__': + + + for index, filename_idx in enumerate([0, 2, 3, 5]): + filename = sorted(os.listdir('../../../data/mount_data/sanmartin/softgrid_1x16/'))[filename_idx] + folder = '../../../data/mount_data/sanmartin/softgrid_1x16/'+filename + print('new file: ' + filename) + name = main(folder, filename) + + if filename not in os.listdir('../data/'): + os.mkdir('../data/' + filename) + + np.save('../data/' + filename + '/fish_species.npy', name) diff --git a/eventdetection.py b/eventdetection.py new file mode 100644 index 0000000..ff2e36f --- /dev/null +++ b/eventdetection.py @@ -0,0 +1,1438 @@ +""" +Detect and handle peaks and troughs as well as threshold crossings in data arrays. + +## Peak detection + +- `detect_peaks()`: peak and trough detection with a relative threshold. +- `peak_width()`: compute width of each peak. +- `peak_size_width()`: compute for each peak its size and width. + +## Threshold crossings + +- `threshold_crossings()`: detect crossings of an absolute threshold. +- `threshold_crossing_times()`: compute times of threshold crossings by linear interpolation. + +## Event manipulation + +- `trim()`: make the list of peaks and troughs the same length. +- `trim_to_peak()`: ensure that the peak is first. +- `trim_closest()`: ensure that peaks minus troughs is smallest. + +- `merge_events()`: Merge events if they are closer than a minimum distance. +- `remove_events()`: Remove events that are too short or too long. +- `widen_events()`: Enlarge events on both sides without overlap. + +## Threshold estimation + +- `std_threshold()`: estimate detection threshold based on the standard deviation. +- `median_std_threshold()`: estimate detection threshold based on the median standard deviation of data snippets. +- `hist_threshold()`: esimate detection threshold based on a histogram of the data. +- `minmax_threshold()`: estimate detection threshold based on maximum minus minimum value. +- `percentile_threshold()`: estimate detection threshold based on interpercentile range. + +## Snippets + +- `snippets()`: cut out data snippets around a list of indices. + +## Peak detection with dynamic threshold: + +- `detect_dynamic_peaks()`: peak and trough detection with a dynamically adapted threshold. +- `accept_peak_size_threshold()`: adapt the dection threshold to the size of the detected peaks. +""" + +import sys +import numpy as np + +try: + from numba import jit, int64 + index_type = int64 +except ImportError: + def jit(*args, **kwargs): + def decorator_jit(func): + return func + return decorator_jit + index_type = np.int + + +def detect_peaks(data, threshold): + """ + Detect peaks and troughs using a relative threshold according to + Bryan S. Todd and David C. Andrews (1999): The identification of peaks in physiological signals. + Computers and Biomedical Research 32, 322-335. + + Parameters + ---------- + data: array + An 1-D array of input data where peaks are detected. + threshold: float or array + A positive number or array of numbers setting the detection threshold, + i.e. the minimum distance between peaks and troughs. + In case of an array make sure that the threshold does not change faster + than the expected intervals between peaks and troughs. + + Returns + ------- + peak_array: array of ints + A list of indices of detected peaks. + trough_array: array of ints + A list of indices of detected troughs. + + Raises + ------ + ValueError: If `threshold <= 0`. + IndexError: If `data` and `threshold` arrays differ in length. + """ + + if np.isscalar(threshold): + if threshold <= 0: + raise ValueError('input argument threshold must be positive!') + return detect_peaks_flat(data, threshold) + else: + if len(data) != len(threshold): + raise IndexError('input arrays data and threshold must have same length!') + if np.min(threshold) <= 0: + raise ValueError('input argument threshold must be positive!') + return detect_peaks_array(data, threshold) + + +@jit(nopython=True) +def detect_peaks_flat(data, thresh): + """ + Detect peaks and troughs using a fixed, relative threshold. + + Parameters + ---------- + data: array + An 1-D array of input data where peaks are detected. + threshold: float + A positive number setting the detection threshold, + i.e. the minimum distance between peaks and troughs. + + Returns + ------- + peak_array: array of ints + A list of indices of detected peaks. + trough_array: array of ints + A list of indices of detected troughs. + + """ + + # initialize: + direction, min_inx, max_inx, pi, ti = 0, 0, 0, 0, 0 + min_value, max_value = data[0], data[0] + peaks_list = np.zeros(len(data)//2,dtype=index_type) + troughs_list = np.zeros(len(data)//2,dtype=index_type) + + # loop through the data: + for index, value in enumerate(data): + + peaks_list, troughs_list, pi, ti, direction, min_inx, max_inx, min_value, max_value = analyse_for_peaks(index, value, thresh, peaks_list, + troughs_list, pi, ti, direction, min_inx, max_inx, min_value, max_value) + + return peaks_list[:pi], troughs_list[:ti] + + +@jit(nopython=True) +def detect_peaks_array(data, threshold): + """ + Detect peaks and troughs using a variable relative threshold. + + Parameters + ---------- + data: array + An 1-D array of input data where peaks are detected. + threshold: array + A positive array of numbers setting the detection threshold, + i.e. the minimum distance between peaks and troughs. + + Returns + ------- + peak_array: array of ints + A list of indices of detected peaks. + trough_array: array of ints + A list of indices of detected troughs. + + """ + + # initialize: + direction, min_inx, max_inx, pi, ti = 0, 0, 0, 0, 0 + min_value, max_value = data[0], data[0] + peaks_list = np.zeros(len(data)//2,dtype=index_type) + troughs_list = np.zeros(len(data)//2,dtype=index_type) + + + # loop through the data: + for index, value in enumerate(data): + + thresh = threshold[index] + peaks_list, troughs_list, pi, ti, direction, min_inx, max_inx, min_value, max_value = analyse_for_peaks(index, value, thresh, peaks_list, + troughs_list, pi, ti, direction, min_inx, max_inx, min_value, max_value) + + return peaks_list[:pi], troughs_list[:ti] + + +@jit(nopython=True) +def analyse_for_peaks(index, value, thresh, peaks_list, troughs_list, pi, ti, direction, min_inx, max_inx, min_value, max_value): + """ + Detect if the current datapoint is a peak or threshold and add it to the existing peak or trough list. Method is taken from: + Bryan S. Todd and David C. Andrews (1999): The identification of peaks in physiological signals. + Computers and Biomedical Research 32, 322-335. + + Parameters + ---------- + index: int + Index of current datapoint. + value: float + Value of current datapoint. + thresh: float + Peak/trough threshold for current datapoint. + peaks_list: numpy array + Numpy array with peak indices prior to current datapoint. + troughs_list: numpy array + Numpy array with trough indices prior to current datapoint. + pi: int + Index of next peak in peaks_list. + ti: int + Index of next trough in troughs_list. + direction: int + Direction of data derivative, e.g. +1 for up and -1 for down. + min_inx: int + Index of the minimum value in the current direction. + max_inx: int + Index of the maximum value in the current direction. + min_value: float + Minimum value in the current direction. + max_value: float + Maximum value in the current direction. + + Returns + ------- + peaks_list: numpy array + Numpy array with peak indices including current datapoint. + troughs_list: numpy array + Numpy array with trough indices including current datapoint. + pi: int + Index of next peak in peaks_list. + ti: int + Index of next trough in troughs_list. + direction: int + Direction of data derivative, e.g. +1 for up and -1 for down. + min_inx: int + Index of the minimum value in the current direction. + max_inx: int + Index of the maximum value in the current direction. + min_value: float + Minimum value in the current direction. + max_value: float + Maximum value in the current direction. + """ + + # rising? + if direction > 0: + if value > max_value: + # update maximum element: + max_inx = index + max_value = value + # otherwise, if the new value is falling below + # the maximum value minus the threshold: + # the maximum is a peak! + elif value <= max_value - thresh: + peaks_list[pi] = max_inx + pi = pi + 1 + # change direction: + direction = -1 + # store minimum element: + min_inx = index + min_value = value + + # falling? + elif direction < 0: + if value < min_value: + # update minimum element: + min_inx = index + min_value = value + # otherwise, if the new value is rising above + # the minimum value plus the threshold: + # the minimum is a trough! + elif value >= min_value + thresh: + troughs_list[ti] = min_inx + ti = ti + 1 + # change direction: + direction = +1 + # store maximum element: + max_inx = index + max_value = value + + # don't know direction yet: + else: + if value <= max_value - thresh: + direction = -1 # falling + elif value >= min_value + thresh: + direction = 1 # rising + + if value > max_value: + # update maximum element: + max_inx = index + max_value = value + elif value < min_value: + # update minimum element: + min_inx = index + min_value = value + + return peaks_list, troughs_list, pi, ti ,direction, min_inx, max_inx,min_value, max_value + + +def detect_peaks_fast(data, threshold): + """Experimental. Try to make algorithm faster. + + Yeah, this is more than three times as fast as `detect_peaks()` with the for loops! + """ + + peaks_list = [] + troughs_list = [] + + # initialize: + max_value = np.maximum.accumulate(data) + min_value = np.minimum.accumulate(data) + falling_idx = np.where(data <= max_value - threshold)[0] + raising_idx = np.where(data >= min_value + threshold)[0] + direction = 1 + if len(falling_idx) > 0 and (len(raising_idx) == 0 or falling_idx[0] < raising_idx[0]): + direction = -1 + index = falling_idx[0] + else: + index = raising_idx[0] + # find peaks and troughs: + while index < len(data): + if direction > 0: + max_value = np.maximum.accumulate(data[index:]) + idx = np.argmax(data[index:] <= max_value - threshold) + if data[index+idx] <= max_value[idx] - threshold: + indices = np.where(max_value[:idx] != max_value[idx])[0] + if len(indices) > 0: + index += indices[-1] + 1 + peaks_list.append(index) + direction = -1 + continue + else: + min_value = np.minimum.accumulate(data[index:]) + idx = np.argmax(data[index:] >= min_value + threshold) + if data[index+idx] >= min_value[idx] + threshold: + indices = np.where(min_value[:idx] != min_value[idx])[0] + if len(indices) > 0: + index += indices[-1] + 1 + troughs_list.append(index) + direction = +1 + continue + break + + return np.asarray(peaks_list), np.asarray(troughs_list) + + +def peak_width(time, data, peak_indices, trough_indices, + peak_frac=0.5, base='max'): + """ + Width of each peak. + + Peak width is computed from interpolated threshold crossings at + `peak_frac` hieght of each peak. + + Parameters + ---------- + time: array + Time, must not be `None`. + data: array + The data with the peaks. + peak_indices: array + Indices of the peaks. + trough_indices: array + Indices of corresponding troughs. + peak_frac: float + Fraction of peak height where its width is measured. + base: string + Height and width of peak is measured relative to + - 'left': trough to the left + - 'right': trough to the right + - 'min': the minimum of the two troughs to the left and to the right + - 'max': the maximum of the two troughs to the left and to the right + - 'mean': mean of the throughs to the left and to the rigth + - 'closest': trough that is closest to peak + + Returns + ------- + widths: array + Width at `peak_frac` height of each peak. + + Raises + ------ + ValueError: + If an invalid value is passed to `base`. + """ + def left_base(data, left_inx, right_inx, peak_inx): + return data[left_inx] + def right_base(data, left_inx, right_inx, peak_inx): + return data[right_inx] + def min_base(data, left_inx, right_inx, peak_inx): + return min(data[left_inx], data[right_inx]) + def max_base(data, left_inx, right_inx, peak_inx): + return max(data[left_inx], data[right_inx]) + def mean_base(data, left_inx, right_inx, peak_inx): + return np.mean((data[left_inx], data[right_inx])) + def closest_base(data, left_inx, right_inx, peak_inx): + return data[left_inx] if peak_inx-left_inx <= right_inx-peak_inx else data[right_inx] + + widths = np.zeros(len(peak_indices)) + if len(peak_indices) == 0: + return widths + # we need a trough before and after each peak: + peak_inx = np.asarray(peak_indices, dtype=int) + trough_inx = np.asarray(trough_indices, dtype=int) + if len(trough_inx) == 0 or peak_inx[0] < trough_inx[0]: + trough_inx = np.hstack((0, trough_inx)) + if peak_inx[-1] > trough_inx[-1]: + trough_inx = np.hstack((trough_inx, len(data)-1)) + # base for size of peaks: + base_func = closest_base + if base == 'left': + base_func = left_base + elif base == 'right': + base_func = right_base + elif base == 'min': + base_func = min_base + elif base == 'max': + base_func = max_base + elif base == 'mean': + base_func = mean_base + elif base == 'closest': + base_func = closest_base + else: + raise ValueError('Invalid value for base (%s)' % base) + # width of peaks: + for j in range(len(peak_inx)): + li = trough_inx[j] + ri = trough_inx[j+1] + baseval = base_func(data, li, ri, peak_inx[j]) + thresh = baseval*(1.0-peak_frac) + data[peak_inx[j]]*peak_frac + inx = li + np.argmax(data[li:ri] > thresh) + if inx > 0: + ti0 = np.interp(thresh, data[inx-1:inx+1], time[inx-1:inx+1]) + else: + ti0 = time[0] + inx = ri - np.argmax(data[ri:li:-1] > thresh) + if inx+1 < len(data): + ti1 = np.interp(thresh, data[inx+1:inx-1:-1], time[inx+1:inx-1:-1]) + else: + ti1 = time[-1] + widths[j] = ti1 - ti0 + return widths + + +def peak_size_width(time, data, peak_indices, trough_indices, + peak_frac=0.75, base='closest'): + """ + Compute for each peak its size and width. + + Parameters + ---------- + time: array + Time, must not be `None`. + data: array + The data with the peaks. + peak_indices: array + Indices of the peaks. + trough_indices: array + Indices of the troughs. + peak_frac: float + Fraction of peak height where its width is measured. + base: string + Height and width of peak is measured relative to + - 'left': trough to the left + - 'right': trough to the right + - 'min': the minimum of the two troughs to the left and to the right + - 'max': the maximum of the two troughs to the left and to the right + - 'mean': mean of the throughs to the left and to the rigth + - 'closest': trough that is closest to peak + + Returns + ------- + peaks: 2-D array + First dimension is the peak index. Second dimension is + time, height (value of data at the peak), + size (peak height minus height of closest trough), + width (at `peak_frac` size), 0.0 (count) of the peak. See `peak_width()`. + + Raises + ------ + ValueError: + If an invalid value is passed to `base`. + """ + def left_base(data, left_inx, right_inx, peak_inx): + return data[left_inx] + def right_base(data, left_inx, right_inx, peak_inx): + return data[right_inx] + def min_base(data, left_inx, right_inx, peak_inx): + return min(data[left_inx], data[right_inx]) + def max_base(data, left_inx, right_inx, peak_inx): + return max(data[left_inx], data[right_inx]) + def mean_base(data, left_inx, right_inx, peak_inx): + return np.mean((data[left_inx], data[right_inx])) + def closest_base(data, left_inx, right_inx, peak_inx): + return data[left_inx] if peak_inx-left_inx <= right_inx-peak_inx else data[right_inx] + + peaks = np.zeros((len(peak_indices), 5)) + if len(peak_indices) == 0: + return peaks + # time point of peaks: + peaks[:, 0] = time[peak_indices] + # height of peaks: + peaks[:, 1] = data[peak_indices] + # we need a trough before and after each peak: + peak_inx = np.asarray(peak_indices, dtype=int) + trough_inx = np.asarray(trough_indices, dtype=int) + if len(trough_inx) == 0 or peak_inx[0] < trough_inx[0]: + trough_inx = np.hstack((0, trough_inx)) + + if peak_inx[-1] > trough_inx[-1]: + trough_inx = np.hstack((trough_inx, len(data)-1)) + # base for size of peaks: + base_func = closest_base + if base == 'left': + base_func = left_base + elif base == 'right': + base_func = right_base + elif base == 'min': + base_func = min_base + elif base == 'max': + base_func = max_base + elif base == 'mean': + base_func = mean_base + elif base == 'closest': + base_func = closest_base + + else: + raise ValueError('Invalid value for base (%s)' % base) + # size and width of peaks: + for j, pi in enumerate(peak_inx): + li = trough_inx[j] + ri = trough_inx[j+1] + baseval = base_func(data, li, ri, pi) + thresh = baseval*(1.0-peak_frac) + data[pi]*peak_frac + inx = li + np.argmax(data[li:ri] > thresh) + if inx > 0: + ti0 = np.interp(thresh, data[inx-1:inx+1], time[inx-1:inx+1]) + else: + ti0 = time[0] + inx = ri - np.argmax(data[ri:li:-1] > thresh) + if inx+1 < len(data): + ti1 = np.interp(thresh, data[inx+1:inx-1:-1], time[inx+1:inx-1:-1]) + else: + ti1 = time[-1] + if np.any(np.isfinite((data[pi], baseval))): + peaks[j, 2] = data[pi] - baseval + peaks[j, 3] = ti1 - ti0 + return peaks + + +def threshold_crossings(data, threshold): + """ + Detect crossings of a threshold with positive and negative slope. + + Parameters + ---------- + data: array + An 1-D array of input data where threshold crossings are detected. + threshold: float or array + A number or array of numbers setting the threshold + that needs to be crossed. + + Returns + ------- + up_indices: array of ints + A list of indices where the threshold is crossed with positive slope. + down_indices: array of ints + A list of indices where the threshold is crossed with negative slope. + + Raises + ------ + IndexError: If `data` and `threshold` arrays differ in length. + """ + + if np.isscalar(threshold): + up_indices = np.nonzero((data[1:]>threshold) & (data[:-1]<=threshold))[0] + down_indices = np.nonzero((data[1:]<=threshold) & (data[:-1]>threshold))[0] + else: + if len(data) != len(threshold): + raise IndexError('input arrays data and threshold must have same length!') + up_indices = np.nonzero((data[1:]>threshold[1:]) & (data[:-1]<=threshold[:-1]))[0] + down_indices = np.nonzero((data[1:]<=threshold[1:]) & (data[:-1]>threshold[:-1]))[0] + return up_indices, down_indices + + +def threshold_crossing_times(time, data, threshold, up_indices, down_indices): + """ + Compute times of threshold crossings by linear interpolation. + + Parameters + ---------- + time: array + Time, must not be `None`. + data: array + The data. + up_indices: array of ints + A list of indices where the threshold is crossed with positive slope. + down_indices: array of ints + A list of indices where the threshold is crossed with negative slope. + + Returns + ------- + up_times: array of floats + Interpolated times where the threshold is crossed with positive slope. + down_times: array of floats + Interpolated times where the threshold is crossed with negative slope. + """ + up_times = np.zeros(len(up_indices)) + for k, inx in enumerate(up_indices): + up_times[k] = np.interp(threshold, data[inx:inx+2], time[inx:inx+2]) + down_times = np.zeros(len(down_indices)) + for k, inx in enumerate(down_indices): + down_times[k] = np.interp(-threshold, -data[inx:inx+2], time[inx:inx+2]) + return up_times, down_times + + +def trim(peaks, troughs): + """ + Trims the peaks and troughs arrays such that they have the same length. + + Parameters + ---------- + peaks: array + List of peak indices or times. + troughs: array + List of trough indices or times. + + Returns + ------- + peaks: array + List of peak indices or times. + troughs: array + List of trough indices or times. + """ + # common len: + n = min(len(peaks), len(troughs)) + # align arrays: + return peaks[:n], troughs[:n] + + +def trim_to_peak(peaks, troughs): + """ + Trims the peaks and troughs arrays such that they have the same length + and the first peak comes first. + + Parameters + ---------- + peaks: array + List of peak indices or times. + troughs: array + List of trough indices or times. + + Returns + ------- + peaks: array + List of peak indices or times. + troughs: array + List of trough indices or times. + """ + # start index for troughs: + tidx = 0 + if len(peaks) > 0 and len(troughs) > 0 and troughs[0] < peaks[0]: + tidx = 1 + # common len: + n = min(len(peaks), len(troughs[tidx:])) + # align arrays: + return peaks[:n], troughs[tidx:tidx + n] + + +def trim_closest(peaks, troughs): + """ + Trims the peaks and troughs arrays such that they have the same length + and that peaks-troughs is on average as small as possible. + + Parameters + ---------- + peaks: array + List of peak indices or times. + troughs: array + List of trough indices or times. + + Returns + ------- + peaks: array + List of peak indices or times. + troughs: array + List of trough indices or times. + """ + pidx = 0 + tidx = 0 + nn = min(len(peaks), len(troughs)) + if nn == 0: + return np.array([]), np.array([]) + dist = np.abs(np.mean(peaks[:nn] - troughs[:nn])) + if len(peaks) == 0 or len(troughs) == 0: + nn = 0 + else: + if peaks[0] < troughs[0]: + nnp = min(len(peaks[1:]), len(troughs)) + distp = np.abs(np.mean(peaks[1:nnp] - troughs[:nnp - 1])) + if distp < dist: + pidx = 1 + nn = nnp + else: + nnt = min(len(peaks), len(troughs[1:])) + distt = np.abs(np.mean(peaks[:nnt - 1] - troughs[1:nnt])) + if distt < dist: + tidx = 1 + nn = nnt + # align arrays: + return peaks[pidx:pidx + nn], troughs[tidx:tidx + nn] + + +def merge_events(onsets, offsets, min_distance): + """Merge events if they are closer than a minimum distance. + + If the beginning of an event (onset, peak, or positive threshold crossing, + is too close to the end of the previous event (offset, trough, or negative + threshold crossing) the two events are merged into a single one that begins + with the first one and ends with the second one. + + Parameters + ---------- + onsets: 1-D array + The onsets (peaks, or positive threshold crossings) of the events + as indices or times. + offsets: 1-D array + The offsets (troughs, or negative threshold crossings) of the events + as indices or times. + min_distance: int or float + The minimum distance between events. If the beginning of an event is separated + from the end of the previous event by less than this distance then the two events + are merged into one. If the event onsets and offsets are given in indices than + min_distance is also in indices. + + Returns + ------- + merged_onsets: 1-D array + The onsets (peaks, or positive threshold crossings) of the merged events + as indices or times according to onsets. + merged_offsets: 1-D array + The offsets (troughs, or negative threshold crossings) of the merged events + as indices or times according to offsets. + """ + onsets, offsets = trim_to_peak(onsets, offsets) + if len(onsets) == 0 or len(offsets) == 0: + return np.array([]), np.array([]) + else: + diff = onsets[1:] - offsets[:-1] + indices = diff > min_distance + merged_onsets = onsets[np.hstack([True, indices])] + merged_offsets = offsets[np.hstack([indices, True])] + return merged_onsets, merged_offsets + + +def remove_events(onsets, offsets, min_duration, max_duration=None): + """Remove events that are too short or too long. + + If the length of an event, i.e. `offset` (offset, trough, or negative + threshold crossing) minus `onset` (onset, peak, or positive threshold crossing), + is shorter than `min_duration` or longer than `max_duration`, then this event is + removed. + + Parameters + ---------- + onsets: 1-D array + The onsets (peaks, or positive threshold crossings) of the events + as indices or times. + offsets: 1-D array + The offsets (troughs, or negative threshold crossings) of the events + as indices or times. + min_duration: int, float, or None + The minimum duration of events. If the event offset minus the event onset + is less than `min_duration`, then the event is removed from the lists. + If the event onsets and offsets are given in indices than + `min_duration` is also in indices. If `None` then this test is skipped. + max_duration: int, float, or None + The maximum duration of events. If the event offset minus the event onset + is larger than `max_duration`, then the event is removed from the lists. + If the event onsets and offsets are given in indices than + `max_duration` is also in indices. If `None` then this test is skipped. + + Returns + ------- + onsets: 1-D array + The onsets (peaks, or positive threshold crossings) of the events + with too short and too long events removed as indices or times according to onsets. + offsets: 1-D array + The offsets (troughs, or negative threshold crossings) of the events + with too short and too long events removed as indices or times according to offsets. + """ + onsets, offsets = trim_to_peak(onsets, offsets) + if len(onsets) == 0 or len(offsets) == 0: + return np.array([]), np.array([]) + elif min_duration is not None or max_duration is not None: + diff = offsets - onsets + if min_duration is not None and max_duration is not None: + indices = (diff > min_duration) & (diff < max_duration) + elif min_duration is not None: + indices = diff > min_duration + else: + indices = diff < max_duration + onsets = onsets[indices] + offsets = offsets[indices] + return onsets, offsets + + +def widen_events(onsets, offsets, max_time, duration): + """Enlarge events on both sides without overlap. + + Subtracts `duration` from the `onsets` and adds `duration` to the offsets. + If two succeeding events are separated by less than two times the `duration`, + then the offset of the previous event and the onset of the following event are + set at the center between the two events. + + Parameters + ---------- + onsets: 1-D array + The onsets (peaks, or positive threshold crossings) of the events + as indices or times. + offsets: 1-D array + The offsets (troughs, or negative threshold crossings) of the events + as indices or times. + max_time: int or float + The maximum value for the end of the last event. + If the event onsets and offsets are given in indices than + max_time is the maximum possible index, i.e. the len of the + data array on which the events where detected. + duration: int or float + The number of indices or the time by which the events should be enlarged. + If the event onsets and offsets are given in indices than + duration is also in indices. + + Returns + ------- + onsets: 1-D array + The onsets (peaks, or positive threshold crossings) of the enlarged events. + offsets: 1-D array + The offsets (troughs, or negative threshold crossings) of the enlarged events. + """ + new_onsets = [] + new_offsets = [] + if len(onsets) > 0: + on_idx = onsets[0] + new_onsets.append( on_idx - duration if on_idx >= duration else 0 ) + for off_idx, on_idx in zip(offsets[:-1], onsets[1:]): + if on_idx - off_idx < 2*duration: + mid_idx = (on_idx + off_idx)//2 + new_offsets.append(mid_idx) + new_onsets.append(mid_idx) + else: + new_offsets.append(off_idx + duration) + new_onsets.append(on_idx - duration) + if len(offsets) > 0: + off_idx = offsets[-1] + new_offsets.append(off_idx + duration if off_idx + duration < max_time else max_time) + return new_onsets, new_offsets + + +def std_threshold(data, samplerate=None, win_size=None, thresh_fac=5.): + """Estimates a threshold for peak detection based on the standard deviation of the data. + + The threshold is computed as the standard deviation of the data + multiplied with `thresh_fac`. + + In case of Gaussian distributed data, setting `thresh_fac=2.0` (two standard deviations) + captures 68% of the data, `thresh_fac=4.0` captures 95%, and `thresh_fac=6.0` 99.7%. + + If `samplerate` and `win_size` is given, then the threshold is computed for + each half-overlapping window of duration `win_size` separately. + In this case the returned threshold is an array of the same size as data. + Without a `samplerate` and `win_size` a single threshold value determined from + the whole data array is returned. + + Parameters + ---------- + data: 1-D array + The data to be analyzed. + samplerate: float or None + Sampling rate of the data in Hz. + win_size: float or None + Size of window in which a threshold value is computed. + thresh_fac: float + Factor by which the standard deviation is multiplied to set the threshold. + + Returns + ------- + threshold: float or 1-D array + The computed threshold. + """ + + if samplerate and win_size: + threshold = np.zeros(len(data)) + win_size_indices = int(win_size * samplerate) + + for inx0 in range(0, len(data), win_size_indices//2): + inx1 = inx0 + win_size_indices + std = np.std(data[inx0:inx1], ddof=1) + threshold[inx0:inx1] = std * thresh_fac + return threshold + else: + return np.std(data, ddof=1) * thresh_fac + + +@jit(nopython=True) +def median_std_threshold(data, samplerate, win_size=0.0005, n_snippets=1000, thresh_fac=6.0): + """Estimate a threshold for peak detection based on the median standard deviation of data snippets. + + On `n_snippets` snippets of `win_size` duration the standard + deviation of the data is estimated. The returned threshold is the + median of these standard deviations multiplied by `thresh_fac`. + + Parameters + ---------- + data: 1-D array of float + The data to be analysed. + samplerate: int or float + Sampling rate of the data + win_size: float + Duration of windows on which standarad deviations are computed in seconds. + n_snippets: int + Number of snippets on which the standard deviations are estimated. + thresh_fac: float + Factor by which the median standard deviation is multiplied to set the threshold. + + Returns + ------- + threshold: float + The computed threshold. + """ + win_size_indices = int(win_size * samplerate) + if win_size_indices < 10: + win_size_indices = 10 + step = len(data)//n_snippets + if step < win_size_indices//2: + step = win_size_indices//2 + stds = np.array([np.std(data[i:i+win_size_indices]) + for i in range(0, len(data)-win_size_indices, step)]) + return np.median(stds)*thresh_fac + + +def hist_threshold(data, samplerate=None, win_size=None, thresh_fac=5., + nbins=100, hist_height=1.0/np.sqrt(np.e)): + """Estimate a threshold for peak detection based on a histogram of the data. + + The standard deviation of the data is estimated from half the + width of the histogram of the data at `hist_height` relative height. + This estimates the data's standard deviation by ignoring tails of the distribution. + + However, you need enough data to robustly estimate the histogram. + + If `samplerate` and `win_size` is given, then the threshold is computed for + each half-overlapping window of duration `win_size` separately. + In this case the returned threshold is an array of the same size as data. + Without a samplerate and win_size a single threshold value determined from + the whole data array is returned. + + Parameters + ---------- + data: 1-D array + The data to be analyzed. + samplerate: float or None + Sampling rate of the data in Hz. + win_size: float or None + Size of window in which a threshold value is computed in sec. + thresh_fac: float + Factor by which the width of the histogram is multiplied to set the threshold. + nbins: int or list of floats + Number of bins or the bins for computing the histogram. + hist_height: float + Height between 0 and 1 at which the width of the histogram is computed. + + Returns + ------- + threshold: float or 1-D array + The computed threshold. + center: float or 1-D array + The center (mean) of the width of the histogram. + """ + + if samplerate and win_size: + threshold = np.zeros(len(data)) + centers = np.zeros(len(data)) + win_size_indices = int(win_size * samplerate) + + for inx0 in range(0, len(data), win_size_indices//2): + inx1 = inx0 + win_size_indices + std, center = hist_threshold(data[inx0:inx1], samplerate=None, win_size=None, + thresh_fac=thresh_fac, nbins=nbins, + hist_height=hist_height) + threshold[inx0:inx1] = std + centers[inx0:inx1] = center + return threshold, centers + else: + maxd = np.max(data) + mind = np.min(data) + contrast = np.abs((maxd - mind)/(maxd + mind)) + if contrast > 1e-8: + hist, bins = np.histogram(data, nbins, density=False) + inx = hist > np.max(hist) * hist_height + lower = bins[0:-1][inx][0] + upper = bins[1:][inx][-1] # needs to return the next bin + center = 0.5 * (lower + upper) + std = 0.5 * (upper - lower) + else: + std = np.std(data) + center = np.mean(data) + return std * thresh_fac, center + + +def minmax_threshold(data, samplerate=None, win_size=None, thresh_fac=0.8): + """Estimate a threshold for peak detection based on minimum and maximum values of the data. + + The threshold is computed as the difference between maximum and + minimum value of the data multiplied with `thresh_fac`. + + If `samplerate` and `win_size` is given, then the threshold is computed for + each half-overlapping window of duration `win_size` separately. + In this case the returned threshold is an array of the same size as data. + Without a samplerate and win_size a single threshold value determined from + the whole data array is returned. + + Parameters + ---------- + data: 1-D array + The data to be analyzed. + samplerate: float or None + Sampling rate of the data in Hz. + win_size: float or None + Size of window in which a threshold value is computed. + thresh_fac: float + Factor by which the difference between minimum and maximum data value + is multiplied to set the threshold. + + Returns + ------- + threshold: float or 1-D array + The computed threshold. + """ + if samplerate and win_size: + threshold = np.zeros(len(data)) + win_size_indices = int(win_size * samplerate) + + for inx0 in range(0, len(data), win_size_indices//2): + inx1 = inx0 + win_size_indices + window_min = np.min(data[inx0:inx1]) + window_max = np.max(data[inx0:inx1]) + threshold[inx0:inx1] = (window_max - window_min) * thresh_fac + return threshold + + else: + return (np.max(data) - np.min(data)) * thresh_fac + + +def percentile_threshold(data, samplerate=None, win_size=None, thresh_fac=1.0, percentile=1.0): + """Estimate a threshold for peak detection based on an inter-percentile range of the data. + + The threshold is computed as the range between the percentile and + 100.0-percentile percentiles of the data multiplied with + thresh_fac. + + For very small values of `percentile` the estimated threshold + approaches the one returned by `minmax_threshold()` (for same values + of `thresh_fac`). For `percentile=16.0` and Gaussian distributed data, + the returned theshold is twice the one returned by `std_threshold()` + or `hist_threshold()`, i.e. twice the standard deviation. + + If you have knowledge about how many data points are in the tails of + the distribution, then this method is preferred over + `hist_threhsold()`. For example, if you expect peaks that you want + to detect using `detect_peaks()` at an average rate of 10Hz and + these peaks are about 1ms wide, then you have a 1ms peak per 100ms + period, i.e. the peaks make up 1% of the distribution. So you should + set `percentile=1.0` or lower. For much lower percentile values, you + might choose `thresh_fac` slightly smaller than one to capture also + smaller peaks. Setting `percentile` slightly higher might not change + the estimated threshold too much, since you start incorporating the + noise floor with a large density, but you may want to set + `thresh_fac` larger than one to reduce false detections. + + If `samplerate` and `win_size` is given, then the threshold is computed for + each half-overlapping window of duration `win_size` separately. + In this case the returned threshold is an array of the same size as data. + Without a samplerate and win_size a single threshold value determined from + the whole data array is returned. + + Parameters + ---------- + data: 1-D array + The data to be analyzed. + samplerate: float or None + Sampling rate of the data in Hz. + win_size: float or None + Size of window in which a threshold value is computed. + percentile: float + The interpercentile range is computed at percentile and 100.0-percentile. + If zero, compute maximum minus minimum data value as the interpercentile range. + thresh_fac: float + Factor by which the inter-percentile range of the data is multiplied to set the threshold. + + Returns + ------- + threshold: float or 1-D array + The computed threshold. + """ + if percentile < 1e-8: + return minmax_threshold(data, samplerate=samplerate, win_size=win_size, + thresh_fac=thresh_fac) + if samplerate and win_size: + threshold = np.zeros(len(data)) + win_size_indices = int(win_size * samplerate) + for inx0 in range(0, len(data), win_size_indices//2): + inx1 = inx0 + win_size_indices + threshold[inx0:inx1] = np.squeeze(np.abs(np.diff( + np.percentile(data[inx0:inx1], [100.0 - percentile, percentile])))) * thresh_fac + return threshold + else: + return np.squeeze(np.abs(np.diff( + np.percentile(data, [100.0 - percentile, percentile])))) * thresh_fac + + +def snippets(data, indices, start=-10, stop=10): + """ + Cut out data arround each position given in indices. + + Parameters + ---------- + data: 1-D array + Data array from which snippets are extracted. + indices: list of int + Indices around which snippets are cut out. + start: int + Each snippet starts at index + start. + stop: int + Each snippet ends at index + stop. + + Returns + ------- + snippet_data: 2-D array + The snippets: first index number of snippet, second index time. + """ + idxs = indices[(indices>=-start) & (indices= len(data): + idx = len(data) - 2 + threshold += (min_thresh - threshold) * (time[idx + 1] - time[idx]) / tau + + # rising? + if direction > 0: + # if the new value is bigger than the old maximum: set it as new maximum: + if value > max_value: + max_inx = index # maximum element + max_value = value + + # otherwise, if the new value is falling below the maximum value minus the threshold: + # the maximum is a peak! + elif max_value >= value + threshold: + # check and update peak with the check_peak_func function: + if check_peak_func: + r, th = check_peak_func(time, data, max_inx, index, + min_inx, threshold, + min_thresh=min_thresh, tau=tau, **kwargs) + if r is not None: + # this really is a peak: + peaks_list.append(r) + if th is not None: + threshold = th + if threshold < min_thresh: + threshold = min_thresh + else: + # this is a peak: + if time is None: + peaks_list.append(max_inx) + else: + peaks_list.append(time[max_inx]) + + # change direction: + min_inx = index # minimum element + min_value = value + direction = -1 + + # falling? + elif direction < 0: + if value < min_value: + min_inx = index # minimum element + min_value = value + + elif value >= min_value + threshold: + # there was a trough: + + # check and update trough with the check_trough function: + if check_trough_func: + r, th = check_trough_func(time, data, min_inx, index, + max_inx, threshold, + min_thresh=min_thresh, tau=tau, **kwargs) + if r is not None: + # this really is a trough: + troughs_list.append(r) + if th is not None: + threshold = th + if threshold < min_thresh: + threshold = min_thresh + else: + # this is a trough: + if time is None: + troughs_list.append(min_inx) + else: + troughs_list.append(time[min_inx]) + + # change direction: + max_inx = index # maximum element + max_value = value + direction = 1 + + # don't know direction yet: + else: + if max_value >= value + threshold: + direction = -1 # falling + elif value >= min_value + threshold: + direction = 1 # rising + + if max_value < value: + max_inx = index # maximum element + max_value = value + + elif value < min_value: + min_inx = index # minimum element + min_value = value + + return np.asarray(peaks_list), np.asarray(troughs_list) + + +def accept_peak_size_threshold(time, data, event_inx, index, min_inx, threshold, + min_thresh, tau, thresh_ampl_fac=0.75, thresh_weight=0.02): + """Accept each detected peak/trough and return its index or time. + Adjust the threshold to the size of the detected peak. + To be passed to the `detect_dynamic_peaks()` function. + + Parameters + ---------- + time: array + Time values, can be `None`. + data: array + The data in wich peaks and troughs are detected. + event_inx: int + Index of the current peak/trough. + index: int + Current index. + min_inx: int + Index of the previous trough/peak. + threshold: float + Threshold value. + min_thresh: float + The minimum value the threshold is allowed to assume.. + tau: float + The time constant of the the decay of the threshold value + given in indices (`time` is `None`) or time units (`time` is not `None`). + thresh_ampl_fac: float + The new threshold is `thresh_ampl_fac` times the size of the current peak. + thresh_weight: float + New threshold is weighted against current threshold with `thresh_weight`. + + Returns + ------- + index: int + Index of the peak/trough if `time` is `None`. + time: int + Time of the peak/trough if `time` is not `None`. + threshold: float + The new threshold to be used. + """ + size = data[event_inx] - data[min_inx] + threshold += thresh_weight * (thresh_ampl_fac * size - threshold) + if time is None: + return event_inx, threshold + else: + return time[event_inx], threshold + + +if __name__ == "__main__": + import matplotlib.pyplot as plt + + print("Checking eventetection module ...") + print('') + # generate data: + dt = 0.001 + time = np.arange(0.0, 10.0, dt) + f = 2.0 + data = (0.5 * np.sin(2.0 * np.pi * f * time) + 0.5) ** 4.0 + data += -0.1 * time * (time - 10.0) + data += 0.1 * np.random.randn(len(data)) + + print("generated waveform with %d peaks" % int(np.round(time[-1] * f))) + plt.plot(time, data) + + print('') + print('check detect_peaks(data, 1.0)...') + peaks, troughs = detect_peaks(data, 1.0) + print(peaks) + print(troughs) + p, t = detect_peaks_fast(data, 1.0) + print(p) + print(t) + # print peaks: + print('detected %d peaks with period %g that differs from the real frequency by %g' % ( + len(peaks), np.mean(np.diff(peaks)), f - 1.0 / np.mean(np.diff(peaks)) / np.mean(np.diff(time)))) + # print troughs: + print('detected %d troughs with period %g that differs from the real frequency by %g' % ( + len(troughs), np.mean(np.diff(troughs)), f - 1.0 / np.mean(np.diff(troughs)) / np.mean(np.diff(time)))) + + # plot peaks and troughs: + plt.plot(time[peaks], data[peaks], '.r', ms=20) + plt.plot(time[troughs], data[troughs], '.g', ms=20) + + # detect threshold crossings: + onsets, offsets = threshold_crossings(data, 3.0) + onsets, offsets = merge_events(onsets, offsets, int(0.5/f/dt)) + plt.plot(time, 3.0*np.ones(len(time)), 'k') + plt.plot(time[onsets], data[onsets], '.c', ms=20) + plt.plot(time[offsets], data[offsets], '.b', ms=20) + + plt.ylim(-0.5, 4.0) + plt.show() + + # check the faster algorithm: + import timeit + def wrapper(func, *args, **kwargs): + def wrapped(): + return func(*args, **kwargs) + return wrapped + wrapped = wrapper(detect_peaks, data, 1.0) + t1 = timeit.timeit(wrapped, number=200) + print(t1) + wrapped = wrapper(detect_peaks_fast, data, 1.0) + t2 = timeit.timeit(wrapped, number=200) + print(t2) + print('new algorithm takes %.0f%% of old algorithm' % (100.0*t2/t1)) diff --git a/helper_functions.py b/helper_functions.py new file mode 100644 index 0000000..cb87737 --- /dev/null +++ b/helper_functions.py @@ -0,0 +1,545 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +import matplotlib.gridspec as gridspec + +from IPython import embed + + +# fill aifl with the IDs +def fill_aifl(id0, id1, aifl, Ch, Ch_connect, matrix_fish_counter, time, frequency, identity, faifl): + """ his function checks where to adds the two identities to the aifl + + 3 cases: + 1: both ids do not exist in aifl - new fish and both ids are added to that fish number + 2: both ids already exist in aifl - check whether the ids exist in the same fish number -- old fish + - if not the same number the traces of both fish numbers are plotted and than + by input of the person either but together as on fish number or added to the faifl + 3: one of the ids exist in aifl - add the other id to the same fish number + + Parameters + ---------- + id0: int + ID of the first trace + id1: int + ID of the second trace + aifl: 3-D array + all identity fish list; 1_D: channel; 2_D: fish number; 3_D: fish identities + Ch: int + first channel + Ch_connect: int + channel of id1 and to which we connect + matrix_fish_counter: int + counter where the new fish should be registered + time: 2-D array + loaded time array; 1_D: channel; 2_D: arbitrary time points [int] + frequency: 2-D array + loaded frequency array; 1_D: channel; 2_D: frequencies [Hz] + identity: 2-D array + loaded ident array; 1_D: channel; 2_D: identities [int] + faifl: 2-D array + false all identity fish list; 1_D: number of how many false enteries in the aifl exist + 2_D: channels and fish number that are falsely connected; + 4 columns: Ch, fishN0, Ch_connect, fishN1 + + Returns + ------- + matrix_fish_counter: int + updated counter + aifl: 3-D array + updated aifl + faifl 2-D array + updated faifl + """ + + # fish number + fishN0 = np.where(aifl[Ch][:][:] == id0) + fishN1 = np.where(aifl[Ch_connect][:][:] == id1) + + # when none of the IDs existes -- new fish + if aifl[Ch][fishN0].size == 0 and aifl[Ch_connect][fishN1].size == 0: + aifl[Ch][matrix_fish_counter][0] = id0 + aifl[Ch_connect][matrix_fish_counter][0] = id1 + matrix_fish_counter += 1 + print('new fish - ', matrix_fish_counter-1) + + # what happens when both IDs already exist in the channels + elif aifl[Ch][fishN0].size != 0 and aifl[Ch_connect][fishN1].size != 0: + try: + # both IDs are already together as one fish + if fishN0[0] == fishN1[0]: + print('old fish') + # the IDs are in two different fish, we have to check is the fish should be merged + # plotting to identify by eye + # then input if the traces should be merged: + # yes/no/false --- if the merge before was already false + else: + print('confused fish', fishN0[0][0], fishN1[0][0]) + + a1 = aifl[:, fishN0[0], :] + a2 = aifl[:, fishN1[0], :] + + for i in range(a1.shape[0]): + if a2[i, 0, ~np.isnan(a2[i, 0, :])].size != 0: + append_counter = 1 + for j in range(a2[i, 0, ~np.isnan(a2[i, 0, :])].size): + print(a2[i][0][j]) + nan_pos = np.where(a1[i, 0, ~np.isnan(a1[i, 0, :])])[0] + if nan_pos.size != 0: + aifl[i, fishN0[0], nan_pos[-1] + append_counter] = a2[i, 0, j] + append_counter += 1 + else: + aifl[i, fishN0[0], 0] = a2[i, 0, j] + aifl[:, fishN1[0], :] = np.nan + + # plot_confused_fish(time, frequency, identity, aifl, id0, id1, Ch, Ch_connect) + # go_signal = input('Are the traces matching? [y/n/f]') + # + # # if the merge before was already false and the fish_number is already in the faifl + # if np.any(faifl[:, [1, 3]] == fishN0[0][0]) or np.any(faifl[:, [1, 3]] == fishN1[0][0]): + # print('Traces are not matching') + # faifl = np.append(faifl, [[Ch, fishN0[0][0], Ch_connect, fishN1[0][0]]], axis=0) + # # go_signal = yes: merge the two fish to one + # elif go_signal == 'y': + # a1 = aifl[:, fishN0[0], :] + # a2 = aifl[:, fishN1[0], :] + # + # for i in range(a1.shape[0]): + # if a2[i, 0, ~np.isnan(a2[i, 0, :])].size != 0: + # append_counter = 1 + # for j in range(a2[i, 0, ~np.isnan(a2[i, 0, :])].size): + # print(a2[i][0][j]) + # nan_pos = np.where(a1[i, 0, ~np.isnan(a1[i, 0, :])])[0] + # if nan_pos.size != 0: + # aifl[i, fishN0[0], nan_pos[-1] + append_counter] = a2[i, 0, j] + # append_counter += 1 + # else: + # aifl[i, fishN0[0], 0] = a2[i, 0, j] + # aifl[:, fishN1[0], :] = np.nan + # + # # go_signal = false: do not merge and put the fish_number into the faifl + # elif go_signal == 'f': + # faifl = np.append(faifl, [[Ch, fishN0[0][0], Ch_connect, fishN1[0][0]]], axis=0) + # + # # go_signal = everything else/no: no merge + # else: + # print('no merge') + except: + embed() + quit() + + # if one of the fish does exist but the other one not: + # if fish0 exists assign fish1 the same fish_number + elif aifl[Ch][fishN0].size != 0: + aifl = add_id_to_aifl(aifl, Ch_connect, id1, fishN0) + + # if fish1 exists assign fish0 the same fish_number + elif aifl[Ch_connect][fishN1].size != 0: + aifl = add_id_to_aifl(aifl, Ch, id0, fishN1) + + return matrix_fish_counter, aifl, faifl + + +def add_id_to_aifl(aifl, Ch, ID, fishN): + """ adds the ID to the fishN in aifl + + Parameters + ---------- + aifl: 3-D array + all identity fish list; 1_D: channel; 2_D: fish number; 3_D: fish identities + Ch: int + Channel + ID: int + the fishID which is added to aifl + fishN: int + the fish number to which we add the ID + + Returns + ------- + aifl: 3-D array + all identity fish list; 1_D: channel; 2_D: fish number; 3_D: fish identities + with the new ID at the fishN + """ + + nan_pos = np.where(aifl[Ch][fishN[0]][~np.isnan(aifl[Ch][fishN[0]])])[0] + if nan_pos.size != 0: + aifl[Ch][fishN[0][0]][nan_pos[-1] + 1] = ID + else: + aifl[Ch][fishN[0][0]][0] = ID + return aifl + + +def plot_confused_fish(time, frequency, identity, aifl, id0, id1, Ch, Ch_next): + """ plots the two traces which should be connected but already exists in two different fish numbers in the aifl + it plots all the traces of this channel and the next of both fish numbers + + plot both traces in question: fish0 -- red; fish1 -- blue + plot the existing traces of the fish_numbers in question + plot traces of fish0 in channel and channel+1 in grey and black + plot traces of fish1 in channel and channel+1 in green and yellow + + Parameters + ---------- + time: 2-D array + loaded time array; 1_D: channel; 2_D: arbitrary time points [int] + frequency: 2-D array + loaded frequency array; 1_D: channel; 2_D: frequencies [Hz] + identity: 2-D array + loaded ident array; 1_D: channel; 2_D: identities [int] + aifl: 3-D array + all identity fish list; 1_D: channel; 2_D: fish number; 3_D: fish identities + id0: int + identity of the first trace + id1: int + identity of the second trace + Ch: int + current channel + Ch_next: int + next channel to which we connect + + Returns + ------- + nothing + + """ + # parameter + t0 = time[Ch][identity[Ch] == id0] + f0 = frequency[Ch][identity[Ch] == id0] + t1 = time[Ch_next][identity[Ch_next] == id1] + f1 = frequency[Ch_next][identity[Ch_next] == id1] + fishN0 = np.where(aifl[Ch][:][:] == id0) + fishN1 = np.where(aifl[Ch_next][:][:] == id1) + + # ToDo: ich hab das == True entfernt, wenns nicht mehr läuft zurück ändern + # ~np.isnan(aifl[connect_channel][where_id1[0]]) == True + for plotidx0 in range(len(np.where(~np.isnan(aifl[Ch][fishN0[0]]))[1])): + plt.plot(time[Ch][identity[Ch] == aifl[Ch][fishN0[0][0]][plotidx0]], + frequency[Ch][identity[Ch] == aifl[Ch][fishN0[0][0]][plotidx0]], + 'grey', LineWidth=6, alpha=0.5) + + for plotidx00 in range(len(np.where(~np.isnan(aifl[Ch_next][fishN0[0]]))[1])): + plt.plot(time[Ch_next][identity[Ch_next] == aifl[Ch_next][fishN0[0][0]][plotidx00]], + frequency[Ch_next][ + identity[Ch_next] == aifl[Ch_next][fishN0[0][0]][plotidx00]], + 'black', LineWidth=5, alpha=0.5) + + for plotidx1 in range(len(np.where(~np.isnan(aifl[Ch_next][fishN1[0]]))[1])): + plt.plot(time[Ch_next][identity[Ch_next] == aifl[Ch_next][fishN1[0][0]][plotidx1]], + frequency[Ch_next][ + identity[Ch_next] == aifl[Ch_next][fishN1[0][0]][plotidx1]], + 'green', LineWidth=4, alpha=0.5) + + for plotidx11 in range(len(np.where(~np.isnan(aifl[Ch][fishN1[0]]))[1])): + plt.plot(time[Ch][identity[Ch] == aifl[Ch][fishN1[0][0]][plotidx11]], + frequency[Ch][identity[Ch] == aifl[Ch][fishN1[0][0]][plotidx11]], + 'yellow', LineWidth=3, alpha=0.5) + + plt.plot(t0, f0, 'red', LineWidth=2) + plt.plot(t1, f1, 'blue', LineWidth=1) + plt.title('ch {}: fish {}; ch {}: fish {}'.format(Ch, fishN0[0][0], Ch_next, fishN1[0][0])) + plt.show() + + +def plot_together(time, frequency, identity, aifl, fishN, farbe): + """ plots all the frequency-time traces of all identities of one fish + + Parameters + ---------- + time: 2-D array + loaded time array; 1_D: channel; 2_D: arbitrary time points [int] + frequency: 2-D array + loaded frequency array; 1_D: channel; 2_D: frequencies [Hz] + identity: 2-D array + loaded ident array; 1_D: channel; 2_D: identities [int] + aifl: 3-D array + all identity fish list; 1_D: channel; 2_D: fish number; 3_D: fish identities + fishN: int + number of the fish for aifl + farbe: str + hex color for the plot + + Returns + ------- + nothing + """ + + fig = plt.figure(figsize=(16, 14)) + fig.suptitle('fish {}'.format(fishN)) + for channel_idx in range(16): + fish1 = aifl[channel_idx, fishN, ~np.isnan(aifl[channel_idx, fishN])] + r1 = len(fish1) + print(fish1) + for len_idx1 in range(r1): + plt.plot(time[channel_idx][identity[channel_idx] == fish1[len_idx1]], + frequency[channel_idx][identity[channel_idx] == fish1[len_idx1]], + Linewidth=1, label=str(channel_idx) + ',' + str(fish1[len_idx1]), color=farbe) + # plt.subplots_adjust(bottom=0.2) + ax = plt.gca() + legend_without_duplicate_labels(fig, ax) + plt.show() + + return ax + + +def plot_all_channels(time, frequency, identity, aifl, fishN1, fishN2=None): + """ plots all the traces of each channel in a different subfigures + + Parameters + ---------- + time: 2-D array + loaded time array; 1_D: channel; 2_D: arbitrary time points [int] + frequency: 2-D array + loaded frequency array; 1_D: channel; 2_D: frequencies [Hz] + identity: 2-D array + loaded ident array; 1_D: channel; 2_D: identities [int] + aifl: 3-D array + all identity fish list; 1_D: channel; 2_D: fish number; 3_D: fish identities + fishN1: int + fish number of the first fish + fishN2: int; optional + fish number of the second fish + + Returns + ------- + nothing + """ + + xcounter = 0 + ycounter = 0 + + fig1, axs1 = plt.subplots(nrows=4, ncols=4, sharex=True, sharey=True, figsize=(16, 14)) + fig1.suptitle('fish {}; fish {}'.format(fishN1, fishN2)) + for channel_idx in range(16): + fish1 = aifl[channel_idx, fishN1, ~np.isnan(aifl[channel_idx, fishN1])] + r1 = len(fish1) + print(fish1) + for len_idx1 in range(r1): + axs1[ycounter, xcounter].plot(time[channel_idx][identity[channel_idx] == fish1[len_idx1]], + frequency[channel_idx][identity[channel_idx] == fish1[len_idx1]], + 'gray', Linewidth=3) + # if fishN2 is given + if fishN2 is not None: + fish2 = aifl[channel_idx, fishN2, ~np.isnan(aifl[channel_idx, fishN2])] + r2 = len(fish2) + for len_idx2 in range(r2): + axs1[ycounter, xcounter].plot(time[channel_idx][identity[channel_idx] == fish2[len_idx2]], + frequency[channel_idx][identity[channel_idx] == fish2[len_idx2]], + 'blue', Linewidth=1) + if xcounter == 3: + xcounter = 0 + ycounter = ycounter + 1 + else: + xcounter = xcounter + 1 + plt.show() + + +def legend_without_duplicate_labels(fig, ax): + """ creats a legend without duplicated labels + + Parameters + ---------- + fig: figure handle + ax: ax handle + + """ + handles, labels = ax.get_legend_handles_labels() + unique = [(h, l) for lwdl_idx, (h, l) in enumerate(zip(handles, labels)) if l not in labels[:lwdl_idx]] + ax.legend(*zip(*unique), bbox_to_anchor=(1, 0.4), loc="center right", bbox_transform=fig.transFigure, ncol=1) + + +def running_3binsizes(ipp, sampling_rate): + """ + + Parameters + ---------- + ipp: 2-D array + interpolated power pancake; power of the signal over time for each channel + 1_D: channel; 2_D: power of the signal over time + sampling_rate: int + sampling rate of the recordings + + Returns + ------- + run_mean: 2-D array + running mean with different bin sizes, + first dimension: mean, second: bin size (from the running function) + run_std: 2-D array + running std with different bin sizes, + first dimension: std, second: bin size (from the running function, same bin size as mean) + + """ + + # lists + running_mean = [] + running_std = [] + + # max power in which channel of ipp + max_ch = np.argmax(ipp, axis=1) + + # bin sizes + bin1 = int(np.floor(15 * 60 / sampling_rate)) + bin2 = int(np.floor(30 * 60 / sampling_rate)) + bin3 = int(np.floor(60 * 60 / sampling_rate)) + + # steps with steps 1/2 of bin size + steps1 = np.arange(0, int(len(max_ch) - bin1), 7.5 * 60) + steps2 = np.arange(0, int(len(max_ch) - bin2), 15 * 60) + steps3 = np.arange(0, int(len(max_ch) - bin2), 30 * 60) + + # make running mean, std and sem + for bins, step in zip([bin1, bin2, bin3], [steps1, steps2, steps3]): + bin_mean = np.full([len(max_ch)], np.nan) + bin_std = np.full([len(max_ch)], np.nan) + + # for i in range(int(len(max_ch) - bins)): + for i in step: + i = int(i) + bin_mean[int(i + np.floor(bins / 2))] = np.mean(max_ch[i:bins + i]) + bin_std[int(i + np.floor(bins / 2))] = np.std(max_ch[i:bins + i]) + + running_mean.append(bin_mean) + running_std.append(bin_std) + + return running_mean, running_std + + +def running_variablebinsize(ipp, sampling_rate): + ''' calculates the running mean and running std of the interpolated power pancake + with two different bin sizes and a step size half of the bin size + + Parameters + ---------- + ipp: 2-D array + interpolated power pancake; power of the signal over time for each channel + 1_D: channel; 2_D: power of the signal over time + sampling_rate: int + sampling rate of the recordings + + Returns + ------- + running_mean: 2-D array + running mean over the time axis with two different bin sizes and step sizes, which are also + dependent on the length of the time trace; + 1_D: different bin size [2]; 2_D: MEAN, (with steps half of the bin size) + running_std: 2-D array + running std over the time axis with two different bin sizes and step sizes, which are also + dependent on the length of the time trace; + 1_D: different bin size [2]; 2_D: STD, (with steps half of the bin size) + + ''' + + # lists + running_mean = [] + running_std = [] + + # max power in which channel of ipp + max_ch = np.argmax(ipp, axis=1) + + # calculate bin sizes + if len(ipp) <= 1800 / sampling_rate: # all time traces < 30 min + bin1 = int(np.floor(60 / sampling_rate)) + bin2 = int(np.floor(2 * 60 / sampling_rate)) + steps1 = np.arange(0, int(len(max_ch) - bin1), 0.5 * 60) + steps2 = np.arange(0, int(len(max_ch) - bin2), 1 * 60) + + elif len(ipp) <= 14400 / sampling_rate: # all time traces < 4 h + bin1 = int(np.floor(10 * 60 / sampling_rate)) + bin2 = int(np.floor(20 * 60 / sampling_rate)) + steps1 = np.arange(0, int(len(max_ch) - bin1), 5 * 60) + steps2 = np.arange(0, int(len(max_ch) - bin2), 10 * 60) + + elif len(ipp) > 14400 / sampling_rate: # all time traces > 4 h + bin1 = int(np.floor(60 * 60 / sampling_rate)) + bin2 = int(np.floor(180 * 60 / sampling_rate)) + steps1 = np.arange(0, int(len(max_ch) - bin1), 30 * 60) + steps2 = np.arange(0, int(len(max_ch) - bin2), 90 * 60) + + # make running mean, std and sem + for bins, step in zip([bin1, bin2], [steps1, steps2]): + bin_mean = np.full([len(max_ch)], np.nan) + bin_std = np.full([len(max_ch)], np.nan) + + # for i in range(int(len(max_ch) - bins)): + for i in step: + i = int(i) + bin_mean[int(i + np.floor(bins / 2))] = np.mean(max_ch[i:bins + i]) + bin_std[int(i + np.floor(bins / 2))] = np.std(max_ch[i:bins + i]) + + running_mean.append(bin_mean) + running_std.append(bin_std) + running_mean = np.array(running_mean) + running_std = np.array(running_std) + + return running_mean, running_std + + +def plot_running(x, run_mean, run_std, threshold, fish, farbe, fs, fs_ticks, lw): + """ + + Parameters + ---------- + x: 1-D array + x axis points in date time format + run_mean: 2-D array + running mean with different bin sizes, + first dimension: mean, second: bin size (from the running function) + run_std: 2-D array + running std with different bin sizes, + first dimension: std, second: bin size (from the running function, same bin size as mean) + threshold: 1-D array + threshold for the activity phases dependent on the run mean, + length of array same as how many bin sizes where calculated + fish: int + title for the figure + farbe: 1-D list + list with all available colors from the matplotlib.colors in hex-color code + fs: int; fontsize + fs_ticks: int; fontsize of ticks + lw: int; line width + + Returns + ------- + + """ + + fig = plt.figure(figsize=(16, 14)) + spec = gridspec.GridSpec(ncols=2, nrows=2, figure=fig) + ax1 = fig.add_subplot(spec[0, 0]) + ax2 = fig.add_subplot(spec[0, 1]) + ax3 = fig.add_subplot(spec[1, 0]) + ax4 = fig.add_subplot(spec[1, 1]) + + plt.suptitle('fish {}'.format(fish), fontsize=fs + 2) + ax_cntr = 0 + + for ax in [ax1, ax2, ax3]: + non_nan = np.arange(len(run_mean[ax_cntr]))[~np.isnan(run_mean[ax_cntr])] + + ax.plot(x[fish][non_nan], run_mean[ax_cntr][non_nan], '.', color=farbe[ax_cntr + 4]) + ax.fill_between(x[fish][non_nan], run_mean[ax_cntr][non_nan] + run_std[ax_cntr][non_nan], + run_mean[ax_cntr][non_nan] - run_std[ax_cntr][non_nan], + facecolor=farbe[ax_cntr + 4], alpha=0.5) + + ax4.plot(x[fish][non_nan], run_std[ax_cntr][non_nan], color=farbe[ax_cntr + 4], + label=threshold[ax_cntr]) + + ax4.plot([x[fish][non_nan][0], x[fish][non_nan][-1]], + [threshold[ax_cntr], threshold[ax_cntr]], + color=farbe[ax_cntr + 4]) + ax_cntr += 1 + + ax.set_ylim([0, 15]) + ax.invert_yaxis() + + for ax in [ax1, ax2, ax3, ax4]: + ax.tick_params(width=lw - 2) + ax.tick_params(axis='both', which='major', labelsize=fs_ticks) + ax.xaxis_date() + date_format = mdates.DateFormatter('%H:%M') + ax.xaxis.set_major_formatter(date_format) + fig.autofmt_xdate() + + ax1.set_xlabel('time [h]', fontsize=fs) + ax1.set_ylabel('channel', fontsize=fs) + + plt.legend() diff --git a/match_oofl.py b/match_oofl.py new file mode 100644 index 0000000..6a03c9f --- /dev/null +++ b/match_oofl.py @@ -0,0 +1,130 @@ +import numpy as np +import matplotlib.pyplot as plt +from IPython import embed +import helper_functions as hf + +if __name__ == '__main__': + + ################################################################################################################### + # load data + ################################################################################################################### + # load all the data of one day + ident = np.load('../../../data/2019-10-07-18_28/all_ident_v.npy', allow_pickle=True) + # power = np.load('../../../data/2019-10-07-18_28/all_sign_v.npy', allow_pickle=True) + freq = np.load('../../../data/2019-10-07-18_28/all_fund_v.npy', allow_pickle=True) + timeidx = np.load('../../../data/2019-10-07-18_28/all_idx_v.npy', allow_pickle=True) + + aifl = np.load('../data/aifl.npy', allow_pickle=True) + oofl = np.load('../data/oofl2.npy', allow_pickle=True) + faifl = np.load('../data/faifl.npy', allow_pickle=True) + faifl = np.delete(faifl, [0], axis=0) + + uni_ident_list = [] + for index in range(len(ident)): + uni_ident = np.unique(ident[index][~np.isnan(ident[index])]) + uni_ident_list.append(uni_ident) + for k in range(len(oofl)-len(oofl[np.where(oofl[:, 0] == 15)[0][0]:, :])): + channel = int(oofl[k, 0]) + + # find for each identity the first and last time stamp in the given channel and the next + # channel 0 + t_Ch0 = np.empty((1, 2)) + id0 = oofl[k,1] + t = timeidx[channel][ident[channel] == id0] + t_Ch0[0][0] = t[0] + t_Ch0[0][1] = t[-1] + + # channel 1 + t_Ch1 = np.empty((len(uni_ident_list[channel + 1]), 2)) + for index in range(len(uni_ident_list[channel + 1])): + id1 = uni_ident_list[channel + 1][index] + t1 = timeidx[channel + 1][ident[channel + 1] == id1] + t_Ch1[index, 0] = t1[0] + t_Ch1[index, 1] = t1[-1] + + # parameters + jitter_time = 5 + tolerance_time = 120 + + false_count = 0 + true_count = 0 + for i in range(len(t_Ch0)): + for j in range(len(t_Ch1)): + t0 = timeidx[channel][ident[channel] == id0] + f0 = freq[channel][ident[channel] == id0] + + id1 = uni_ident_list[channel + 1][j] + t1 = timeidx[channel + 1][ident[channel + 1] == id1] + f1 = freq[channel + 1][ident[channel + 1] == id1] + + t00 = t_Ch0[i][0] - jitter_time + t0n = t_Ch0[i][1] + jitter_time + + t10 = t_Ch1[j][0] - jitter_time + t1n = t_Ch1[j][1] + jitter_time + + window_times = sorted(np.array([t00, t0n, t10, t1n]))[1:3] + + sort_it = False + case = np.nan + if (t00 <= t10) and (t00 <= t1n) and (t0n >= t1n): + sort_it = True + case = 1 + elif (t10 <= t00) and (t10 <= t0n) and (t1n >= t0n): + sort_it = True + case = 2 + elif (t00 <= t10) and (t0n >= t10) and (t0n <= t1n): + sort_it = True + case = 3 + elif (t10 <= t00) and (t1n >= t00) and (t1n <= t0n): + sort_it = True + case = 4 + else: + false_count += 1 + pass + + # if t_Ch0[i][0] - jitter_time <= t_Ch1[j][0] <= t_Ch0[i][1] + jitter_time: + if sort_it: + + true_count+=1 + f0_box_min = np.median(f0[np.isclose(t0, window_times[0], atol=tolerance_time)]) + f0_box_max = np.median(f0[np.isclose(t0, window_times[1], atol=tolerance_time)]) + f1_box_min = np.median(f1[np.isclose(t1, window_times[0], atol=tolerance_time)]) + f1_box_max = np.median(f1[np.isclose(t1, window_times[1], atol=tolerance_time)]) + if f0_box_min.size > 0 and f1_box_min.size > 0: + fdiff0 = abs(f0_box_min - f1_box_min) + if fdiff0 <= 10: + plt.plot(timeidx[channel][ident[channel] == id0], + freq[channel][ident[channel] == id0], label=id0) + plt.plot(timeidx[channel + 1][ident[channel + 1] == id1], + freq[channel + 1][ident[channel + 1] == id1], label=id1) + plt.legend() + plt.show() + embed() + #hf.fill_aifl(id0, id1, aifl, channel, 4, timeidx, freq, ident, faifl) + else: + continue + elif f0_box_max.size > 0 and f1_box_max.size > 0: + fdiff1 = abs(f0_box_max - f1_box_max) + if fdiff1 <= 10: + plt.plot(timeidx[channel][ident[channel] == id0], + freq[channel][ident[channel] == id0], label=id0) + plt.plot(timeidx[channel + 1][ident[channel + 1] == id1], + freq[channel + 1][ident[channel + 1] == id1], label=id1) + plt.legend() + plt.show() + embed() + #hf.fill_aifl(id0,id1,aifl,channel,4,timeidx,freq, ident,faifl) + else: + continue + + oofl = [] + counter = 0 + for idx in range(len(uni_ident_list)): + for j in range(len(uni_ident_list[idx])): + if not np.any(aifl[idx] == uni_ident_list[idx][j]): + oofl.append([idx, uni_ident_list[idx][j]]) + counter = counter + 1 + + embed() + quit() \ No newline at end of file