diff --git a/code/chirpdetection.py b/code/chirpdetection.py index c4a04e7..63b6ebe 100755 --- a/code/chirpdetection.py +++ b/code/chirpdetection.py @@ -130,7 +130,7 @@ class PlotBuffer: data_oi, self.data.raw_rate, self.t0 - 5, - [np.max(self.frequency) - 200, np.max(self.frequency) + 200] + [np.min(self.frequency) - 200, np.max(self.frequency) + 200] ) for track_id in self.data.ids: @@ -145,14 +145,15 @@ class PlotBuffer: # get tracked frequencies and their times f = self.data.freq[window_idx] - t = self.data.time[ - self.data.idx[self.data.ident == self.track_id]] - tmask = (t >= t0_track) & (t <= (t0_track + dt_track)) + # t = self.data.time[ + # self.data.idx[self.data.ident == self.track_id]] + # tmask = (t >= t0_track) & (t <= (t0_track + dt_track)) + t = self.data.time[self.data.idx[window_idx]] if track_id == self.track_id: - ax0.plot(t[tmask]-self.t0_old, f, lw=lw, + ax0.plot(t-self.t0_old, f, lw=lw, zorder=10, color=ps.gblue1) else: - ax0.plot(t[tmask]-self.t0_old, f, lw=lw, + ax0.plot(t-self.t0_old, f, lw=lw, zorder=10, color=ps.gray, alpha=0.5) ax0.fill_between( @@ -472,7 +473,9 @@ def find_searchband( ) # search window in boolean - search_window_bool = np.ones_like(len(search_window), dtype=bool) + bool_lower = np.ones_like(search_window, dtype=bool) + bool_upper = np.ones_like(search_window, dtype=bool) + search_window_bool = np.ones_like(search_window, dtype=bool) # make seperate arrays from the qartiles q25 = np.asarray([i[0] for i in frequency_percentiles]) @@ -492,11 +495,10 @@ def find_searchband( q25_temp = q25[percentiles_ids == check_track_id] q75_temp = q75[percentiles_ids == check_track_id] - print(q25_temp, q75_temp) - - search_window_bool[ - (search_window > q25_temp) & (search_window < q75_temp) - ] = False + bool_lower[search_window > q25_temp - config.search_res] = False + bool_upper[search_window < q75_temp + config.search_res] = False + search_window_bool[(bool_lower == False) & + (bool_upper == False)] = False # find gaps in search window search_window_indices = np.arange(len(search_window)) @@ -552,7 +554,7 @@ def find_searchband( return config.default_search_freq -def main(datapath: str, plot: str) -> None: +def chirpdetection(datapath: str, plot: str) -> None: assert plot in [ "save", @@ -561,6 +563,7 @@ def main(datapath: str, plot: str) -> None: ], "plot must be 'save', 'show' or 'false'" # load raw file + print('datapath', datapath) data = LoadData(datapath) # load config file @@ -589,8 +592,8 @@ def main(datapath: str, plot: str) -> None: raw_time = np.arange(data.raw.shape[0]) / data.raw_rate # good chirp times for data: 2022-06-02-10_00 - # window_start_index = (3 * 60 * 60 + 6 * 60 + 43.5 + 5) * data.raw_rate - # window_duration_index = 60 * data.raw_rate + window_start_index = (3 * 60 * 60 + 6 * 60 + 43.5 + 5) * data.raw_rate + window_duration_index = 60 * data.raw_rate # t0 = 0 # dt = data.raw.shape[0] @@ -651,14 +654,14 @@ def main(datapath: str, plot: str) -> None: # approximate sampling rate to compute expected durations if there # is data available for this time window for this fish id - track_samplerate = np.mean(1 / np.diff(data.time)) - expected_duration = ( - (window_start_seconds + window_duration_seconds) - - window_start_seconds - ) * track_samplerate +# track_samplerate = np.mean(1 / np.diff(data.time)) +# expected_duration = ( +# (window_start_seconds + window_duration_seconds) +# - window_start_seconds +# ) * track_samplerate # check if tracked data available in this window - if len(current_frequencies) < expected_duration / 2: + if len(current_frequencies) < 3: logger.warning( f"Track {track_id} has no data in window {st}, skipping." ) @@ -918,11 +921,9 @@ def main(datapath: str, plot: str) -> None: multielectrode_chirps.append(singleelectrode_chirps) # only initialize the plotting buffer if chirps are detected - chirp_detected = ( - (el == config.number_electrodes - 1) - & (len(singleelectrode_chirps) > 0) - & (plot in ["show", "save"]) - ) + chirp_detected = (el == (config.number_electrodes - 1) + & (plot in ["show", "save"]) + ) if chirp_detected: @@ -987,11 +988,12 @@ def main(datapath: str, plot: str) -> None: # if chirps are detected and the plot flag is set, plot the # chirps, otheswise try to delete the buffer if it exists - if len(multielectrode_chirps_validated) > 0: + if ((len(multielectrode_chirps_validated) > 0) & (plot in ["show", "save"])): try: buffer.plot_buffer(multielectrode_chirps_validated, plot) + del buffer except NameError: - pass + embed() else: try: del buffer @@ -1049,4 +1051,4 @@ if __name__ == "__main__": datapath = "../data/2022-06-02-10_00/" # datapath = "/home/weygoldt/Data/uni/efishdata/2016-colombia/fishgrid/2016-04-09-22_25/" # datapath = "/home/weygoldt/Data/uni/chirpdetection/GP2023_chirp_detection/data/mount_data/2020-03-13-10_00/" - main(datapath, plot="save") + chirpdetection(datapath, plot="show") diff --git a/code/chirpdetector_conf.yml b/code/chirpdetector_conf.yml index 2f4fc9a..af48384 100755 --- a/code/chirpdetector_conf.yml +++ b/code/chirpdetector_conf.yml @@ -3,7 +3,7 @@ dataroot: "../data/" outputdir: "../output/" # Duration and overlap of the analysis window in seconds -window: 10 +window: 5 overlap: 1 edge: 0.25 diff --git a/code/extract_chirps.py b/code/extract_chirps.py new file mode 100644 index 0000000..2180a7f --- /dev/null +++ b/code/extract_chirps.py @@ -0,0 +1,44 @@ +import os +import numpy as np +from chirpdetection import chirpdetection +from IPython import embed + + +def main(datapaths): + + for path in datapaths: + chirpdetection(path, plot='show') + + +if __name__ == '__main__': + + dataroot = '../data/mount_data/' + + datasets = sorted([name for name in os.listdir(dataroot) if os.path.isdir( + os.path.join(dataroot, name))]) + + valid_datasets = [] + + for dataset in datasets: + + path = os.path.join(dataroot, dataset) + csv_name = '-'.join(dataset.split('-')[:3]) + '.csv' + + if os.path.exists(os.path.join(path, csv_name)) is False: + continue + + if os.path.exists(os.path.join(path, 'ident_v.npy')) is False: + continue + + ident = np.load(os.path.join(path, 'ident_v.npy')) + number_of_fish = len(np.unique(ident[~np.isnan(ident)])) + if number_of_fish != 2: + continue + + valid_datasets.append(dataset) + + datapaths = [os.path.join(dataroot, dataset) + + '/' for dataset in valid_datasets] + embed() + + main(datapaths[3]) diff --git a/code/plot_event_timeline.py b/code/plot_event_timeline.py new file mode 100644 index 0000000..6c984be --- /dev/null +++ b/code/plot_event_timeline.py @@ -0,0 +1,203 @@ +import numpy as np + +import os + +import numpy as np +import matplotlib.pyplot as plt +from thunderfish.powerspectrum import decibel + +from IPython import embed +from pandas import read_csv +from modules.logger import makeLogger +from modules.plotstyle import PlotStyle + +ps = PlotStyle() + +logger = makeLogger(__name__) + + +class Behavior: + """Load behavior data from csv file as class attributes + Attributes + ---------- + behavior: 0: chasing onset, 1: chasing offset, 2: physical contact + behavior_type: + behavioral_category: + comment_start: + comment_stop: + dataframe: pandas dataframe with all the data + duration_s: + media_file: + observation_date: + observation_id: + start_s: start time of the event in seconds + stop_s: stop time of the event in seconds + total_length: + """ + + def __init__(self, folder_path: str) -> None: + + + LED_on_time_BORIS = np.load(os.path.join(folder_path, 'LED_on_time.npy'), allow_pickle=True) + + csv_filename = [f for f in os.listdir(folder_path) if f.endswith('.csv')][0] + logger.info(f'CSV file: {csv_filename}') + self.dataframe = read_csv(os.path.join(folder_path, csv_filename)) + + self.chirps = np.load(os.path.join(folder_path, 'chirps.npy'), allow_pickle=True) + self.chirps_ids = np.load(os.path.join(folder_path, 'chirps_ids.npy'), allow_pickle=True) + + self.ident = np.load(os.path.join(folder_path, 'ident_v.npy'), allow_pickle=True) + self.idx = np.load(os.path.join(folder_path, 'idx_v.npy'), allow_pickle=True) + self.freq = np.load(os.path.join(folder_path, 'fund_v.npy'), allow_pickle=True) + self.time = np.load(os.path.join(folder_path, "times.npy"), allow_pickle=True) + self.spec = np.load(os.path.join(folder_path, "spec.npy"), allow_pickle=True) + + for k, key in enumerate(self.dataframe.keys()): + key = key.lower() + if ' ' in key: + key = key.replace(' ', '_') + if '(' in key: + key = key.replace('(', '') + key = key.replace(')', '') + setattr(self, key, np.array(self.dataframe[self.dataframe.keys()[k]])) + + last_LED_t_BORIS = LED_on_time_BORIS[-1] + real_time_range = self.time[-1] - self.time[0] + factor = 1.034141 + shift = last_LED_t_BORIS - real_time_range * factor + self.start_s = (self.start_s - shift) / factor + self.stop_s = (self.stop_s - shift) / factor + +def correct_chasing_events( + category: np.ndarray, + timestamps: np.ndarray + ) -> tuple[np.ndarray, np.ndarray]: + + onset_ids = np.arange( + len(category))[category == 0] + offset_ids = np.arange( + len(category))[category == 1] + + # Check whether on- or offset is longer and calculate length difference + if len(onset_ids) > len(offset_ids): + len_diff = len(onset_ids) - len(offset_ids) + longer_array = onset_ids + shorter_array = offset_ids + logger.info(f'Onsets are greater than offsets by {len_diff}') + elif len(onset_ids) < len(offset_ids): + len_diff = len(offset_ids) - len(onset_ids) + longer_array = offset_ids + shorter_array = onset_ids + logger.info(f'Offsets are greater than offsets by {len_diff}') + elif len(onset_ids) == len(offset_ids): + logger.info('Chasing events are equal') + return category, timestamps + + + # Correct the wrong chasing events; delete double events + wrong_ids = [] + for i in range(len(longer_array)-(len_diff+1)): + if (shorter_array[i] > longer_array[i]) & (shorter_array[i] < longer_array[i+1]): + pass + else: + wrong_ids.append(longer_array[i]) + longer_array = np.delete(longer_array, i) + + category = np.delete( + category, wrong_ids) + timestamps = np.delete( + timestamps, wrong_ids) + return category, timestamps + + + +def main(datapath: str): + # behabvior is pandas dataframe with all the data + bh = Behavior(datapath) + # chirps are not sorted in time (presumably due to prior groupings) + # get and sort chirps and corresponding fish_ids of the chirps + chirps = bh.chirps[np.argsort(bh.chirps)] + chirps_fish_ids = bh.chirps_ids[np.argsort(bh.chirps)] + category = bh.behavior + timestamps = bh.start_s + # Correct for doubles in chasing on- and offsets to get the right on-/offset pairs + # Get rid of tracking faults (two onsets or two offsets after another) + category, timestamps = correct_chasing_events(category, timestamps) + + # split categories + chasing_onset = (timestamps[category == 0]/ 60) /60 + chasing_offset = (timestamps[category == 1]/ 60) /60 + physical_contact = (timestamps[category == 2] / 60) /60 + + all_fish_ids = np.unique(chirps_fish_ids) + fish1_id = all_fish_ids[0] + fish2_id = all_fish_ids[1] + # Associate chirps to inidividual fish + fish1 = (chirps[chirps_fish_ids == fish1_id] / 60) /60 + fish2 = (chirps[chirps_fish_ids == fish2_id] / 60) /60 + fish1_color = ps.red + fish2_color = ps.orange + + fig, ax = plt.subplots(4, 1, figsize=(10, 5), height_ratios=[0.5, 0.5, 0.5, 6], sharex=True) + # marker size + s = 200 + ax[0].scatter(physical_contact, np.ones(len(physical_contact)), color='firebrick', marker='|', s=s) + ax[1].scatter(chasing_onset, np.ones(len(chasing_onset)), color='green', marker='|', s=s ) + ax[2].scatter(fish1, np.ones(len(fish1))-0.25, color=fish1_color, marker='|', s=s) + ax[2].scatter(fish2, np.zeros(len(fish2))+0.25, color=fish2_color, marker='|', s=s) + + + freq_temp = bh.freq[bh.ident==fish1_id] + time_temp = bh.time[bh.idx[bh.ident==fish1_id]] + ax[3].plot((time_temp/ 60) /60, freq_temp, color=fish1_color) + + freq_temp = bh.freq[bh.ident==fish2_id] + time_temp = bh.time[bh.idx[bh.ident==fish2_id]] + ax[3].plot((time_temp/ 60) /60, freq_temp, color=fish2_color) + + #ax[3].imshow(decibel(bh.spec), extent=[bh.time[0]/60/60, bh.time[-1]/60/60, 0, 2000], aspect='auto', origin='lower') + + # Hide grid lines + ax[0].grid(False) + ax[0].set_frame_on(False) + ax[0].set_xticks([]) + ax[0].set_yticks([]) + ps.hide_ax(ax[0]) + + + ax[1].grid(False) + ax[1].set_frame_on(False) + ax[1].set_xticks([]) + ax[1].set_yticks([]) + ps.hide_ax(ax[1]) + + ax[2].grid(False) + ax[2].set_frame_on(False) + ax[2].set_yticks([]) + ax[2].set_xticks([]) + ps.hide_ax(ax[2]) + + + + ax[3].axvspan(0, 3, 0, 5, facecolor='grey', alpha=0.5) + ax[3].set_xticks(np.arange(0, 6.1, 0.5)) + + labelpad = 40 + ax[0].set_ylabel('Physical contact', rotation=0, labelpad=labelpad) + ax[1].set_ylabel('Chasing events', rotation=0, labelpad=labelpad) + ax[2].set_ylabel('Chirps', rotation=0, labelpad=labelpad) + ax[3].set_ylabel('EODf') + + ax[3].set_xlabel('Time [h]') + + plt.show() + embed() + + # plot chirps + + +if __name__ == '__main__': + # Path to the data + datapath = '../data/mount_data/2020-05-13-10_00/' + main(datapath)