Merge branch 'master' into eventtriggeredchirps

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[
+            # t = self.data.time[
-                self.data.idx[self.data.ident == self.track_id]]
+            #     self.data.idx[self.data.ident == self.track_id]]
-            tmask = (t >= t0_track) & (t <= (t0_track + dt_track))
+            # 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)
+            bool_lower[search_window > q25_temp - config.search_res] = False
-
+            bool_upper[search_window < q75_temp + config.search_res] = False
-            search_window_bool[
+            search_window_bool[(bool_lower == False) &
-                (search_window > q25_temp) & (search_window < q75_temp)
+                               (bool_upper == False)] = 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_start_index = (3 * 60 * 60 + 6 * 60 + 43.5 + 5) * data.raw_rate
-    # window_duration_index = 60 * 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))
+#             track_samplerate = np.mean(1 / np.diff(data.time))
-            expected_duration = (
+#             expected_duration = (
-                (window_start_seconds + window_duration_seconds)
+#                 (window_start_seconds + window_duration_seconds)
-                - window_start_seconds
+#                 - window_start_seconds
-            ) * track_samplerate
+#             ) * 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,9 +921,7 @@ def main(datapath: str, plot: str) -> None:
                 multielectrode_chirps.append(singleelectrode_chirps)
 
                 # only initialize the plotting buffer if chirps are detected
-                chirp_detected = (
+                chirp_detected = (el == (config.number_electrodes - 1)
-                    (el == config.number_electrodes - 1)
-                    & (len(singleelectrode_chirps) > 0)
                                   & (plot in ["show", "save"])
                                   )
 
@@ -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")

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
 

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])

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)