From 39c75916be1a5862afc04683d5dad146ffb622d0 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Tue, 31 Oct 2023 13:16:27 -0400 Subject: [PATCH 01/47] refactor + add event-by-event spectrogram --- neuropy/utils/signal_process.py | 40 ++++++++++++++++++++++++++++----- 1 file changed, 35 insertions(+), 5 deletions(-) diff --git a/neuropy/utils/signal_process.py b/neuropy/utils/signal_process.py index 65772411..829f0e8e 100644 --- a/neuropy/utils/signal_process.py +++ b/neuropy/utils/signal_process.py @@ -103,14 +103,14 @@ def get_noisy_spect_bool(self, thresh=5): spect_sum = self.traces.sum(axis=0) return (stats.zscore(spect_sum) >= thresh) | (spect_sum <= 0) - def get_pe_mean_spec( + def get_pe_spec( self, event_times, buffer_sec=(0.5, 0.5), ignore_epochs: core.Epoch = None, print_ignored_frames: bool = True, ): - """Get peri-event mean spectrogram + """Get peri-event spectrogram for every event. Parameters ---------- @@ -122,8 +122,8 @@ def get_pe_mean_spec( Returns ------- - Spectrogram class with traces = mean spectrogram and times ranging from -buffer_sec[0] to buffer_sec[1] - """ + spec: nfreq x nbins x nevents array containing spectrogram for each event + ranging from -buffer_sec[0] to buffer_sec[1]""" event_times = event_times.squeeze() assert event_times.ndim == 1, "event_times must be broadcastable to ndim=1" @@ -171,7 +171,37 @@ def get_pe_mean_spec( ) sxx_list.append(sxx_temp) - sxx_mean = np.nanmean(np.stack(sxx_list, axis=2), axis=2) + return np.stack(sxx_list, axis=2) + + def get_pe_mean_spec( + self, + event_times, + buffer_sec=(0.5, 0.5), + ignore_epochs: core.Epoch = None, + print_ignored_frames: bool = True, + ): + """Get peri-event mean spectrogram + + Parameters + ---------- + event_times: ndarray of floats times of each event in seconds, will be time 0 in mean spectrogram + + buffer_sec: tuple of floats defining amount of time before/after event to grab. + + ignore_epochs: core.Epoch class of epochs to ignore when calculating mean spectrogram + + Returns + ------- + Spectrogram class with traces = mean spectrogram and times ranging from -buffer_sec[0] to buffer_sec[1] + """ + + # Get spectrogram for each event + sxx_by_event = self.get_pe_spec( + event_times, buffer_sec, ignore_epochs, print_ignored_frames + ) + + # Take the mean + sxx_mean = np.nanmean(sxx_by_event, axis=2) return Spectrogram( sxx_mean, From 3288bc9b75c93f42f33cffd139353eea8a2a4e86 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Tue, 31 Oct 2023 13:17:09 -0400 Subject: [PATCH 02/47] enhancement --- neuropy/utils/manipulate_files.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/neuropy/utils/manipulate_files.py b/neuropy/utils/manipulate_files.py index e1e00699..5200fd73 100644 --- a/neuropy/utils/manipulate_files.py +++ b/neuropy/utils/manipulate_files.py @@ -92,6 +92,7 @@ def save_notebook( base_dir: str or Path, save_name: str or None = None, save_prepend: str or None = None, + save_append: str or None = None, ): """Save a jupyter notebook to base directory noted. Make sure you save it beforehand! (Autosave should take care of this, but please check. @@ -101,14 +102,16 @@ def save_notebook( assert isinstance(save_name, (str, type(None))) assert isinstance(save_prepend, (str, type(None))) + assert isinstance(save_append, (str, type(None))) # Assemble/get save_name variables nb_fname = ipynbname.name() if save_name is None else save_name nb_path = Path(base_dir) save_prepend = "" if save_prepend is None else save_prepend + save_append = "" if save_append is None else save_append # Create savename and save file - nb_copy_savename = nb_path / f"{save_prepend}{nb_fname}.ipynb" + nb_copy_savename = nb_path / f"{save_prepend}{nb_fname}{save_append}.ipynb" shutil.copy(str(ipynbname.path()), str(nb_copy_savename)) # Spit it out for confirmation From 3765b1acf3a5596ed5ae99216bef203ea2d36987 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 1 Nov 2023 15:57:47 -0400 Subject: [PATCH 03/47] add Epoch.add_epoch_by_index method --- neuropy/core/epoch.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index e8d4d662..ca6acee8 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -5,6 +5,7 @@ from pathlib import Path import scipy.signal as sg import typing +from copy import deepcopy def _unpack_args(values, fs=1): @@ -105,6 +106,15 @@ def add_epoch_manually(self, start, stop, label="", merge_dt: float or None = 0) else: return self.__add__(Epoch(comb_df)) + def add_epoch_by_index(self, index, start, stop, label=""): + assert np.mod(index, 1) > 0, "index must be a non-integer, e.g. -0.5 or 11.5" + epochs_df = deepcopy(self._epochs) + line = pd.DataFrame( + {"start": start, "stop": stop, "label": label}, index=[index] + ) + epochs_df = pd.concat((epochs_df, line), ignore_index=False) + self._epochs = epochs_df.sort_index().reset_index(drop=True) + def shift(self, dt): epochs = self._epochs.copy() epochs[["start", "stop"]] += dt From 3fa88de368b0b698aea5fb3b2865dd6274aea0a6 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 1 Nov 2023 16:24:47 -0400 Subject: [PATCH 04/47] enh - keep missing start times or those outside time range as nan --- neuropy/utils/signal_process.py | 69 ++++++++++++++++++--------------- 1 file changed, 37 insertions(+), 32 deletions(-) diff --git a/neuropy/utils/signal_process.py b/neuropy/utils/signal_process.py index 829f0e8e..3fdbee85 100644 --- a/neuropy/utils/signal_process.py +++ b/neuropy/utils/signal_process.py @@ -133,42 +133,47 @@ def get_pe_spec( keep_bool = (epoch_start_stops[:, 0] > self.t_start) & ( epoch_start_stops[:, 1] < self.t_stop ) - if np.sum(keep_bool) < len(event_times): - event_times = event_times[keep_bool] - print( - f"Events {np.where(~keep_bool)[0]} are outside of data range and were dropped" - ) + # if np.sum(keep_bool) < len(event_times): + # event_times = event_times[keep_bool] + # print( + # f"Events {np.where(~keep_bool)[0]} are outside of data range and were dropped" + # ) sxx_list = [] ntime_bins = int(np.ceil(np.sum(buffer_sec) * self.sampling_rate)) for t_event in event_times: - start_time = t_event - buffer_sec[0] - stop_time = t_event + buffer_sec[1] - wvlet_temp = self.time_slice(t_start=start_time, t_stop=stop_time) - - # Add/remove one frame at end if array size doesn't match expected number of time bins - if len(wvlet_temp.time) != ntime_bins: - if len(wvlet_temp.time) == (ntime_bins + 1): - stop_time -= 0.5 / self.sampling_rate - elif len(wvlet_temp.time) == (ntime_bins - 1): - stop_time += 0.5 / self.sampling_rate - else: - print("Error - time bins off by more than 1") + if not np.isnan(t_event): + start_time = t_event - buffer_sec[0] + stop_time = t_event + buffer_sec[1] wvlet_temp = self.time_slice(t_start=start_time, t_stop=stop_time) - sxx_temp = wvlet_temp.traces - - # Ignore times if specified - these are likely already NaN in the spectrogram - if ignore_epochs is not None: - time_bins = np.linspace(start_time, stop_time, ntime_bins) - ignore_bool, ignore_times, _ = ignore_epochs.contains(time_bins) - - # Display ignored frames - if np.sum(ignore_bool) > 0: - sxx_temp[:, ignore_bool] = np.nan - if print_ignored_frames: - print( - f"{np.sum(ignore_bool)} frames between {ignore_times.min():.1F} and {ignore_times.max():.1F} ignored (sent to nan)" - ) + + # Add/remove one frame at end if array size doesn't match expected number of time bins + if len(wvlet_temp.time) != ntime_bins: + if len(wvlet_temp.time) == (ntime_bins + 1): + stop_time -= 0.5 / self.sampling_rate + elif len(wvlet_temp.time) == (ntime_bins - 1): + stop_time += 0.5 / self.sampling_rate + else: + print("Error - time bins off by more than 1") + wvlet_temp = self.time_slice(t_start=start_time, t_stop=stop_time) + sxx_temp = wvlet_temp.traces + + # Ignore times if specified - these are likely already NaN in the spectrogram + if ignore_epochs is not None: + time_bins = np.linspace(start_time, stop_time, ntime_bins) + ignore_bool, ignore_times, _ = ignore_epochs.contains(time_bins) + + # Display ignored frames + if np.sum(ignore_bool) > 0: + sxx_temp[:, ignore_bool] = np.nan + if print_ignored_frames: + print( + f"{np.sum(ignore_bool)} frames between {ignore_times.min():.1F} and {ignore_times.max():.1F} ignored (sent to nan)" + ) + + else: # Fill in with nans if event is missing + wvlet_temp = self.time_slice(t_start=1, t_stop=1 + np.sum(buffer_sec)) + sxx_temp = np.ones_like(wvlet_temp.traces) * np.nan sxx_list.append(sxx_temp) return np.stack(sxx_list, axis=2) @@ -180,7 +185,7 @@ def get_pe_mean_spec( ignore_epochs: core.Epoch = None, print_ignored_frames: bool = True, ): - """Get peri-event mean spectrogram + """Get peri-event mean spectrogram. Parameters ---------- From e45b61f5f48121416c77ea034370e8b08ece9962 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 1 Nov 2023 16:42:04 -0400 Subject: [PATCH 05/47] fixed bug from previous commit - not well tested yet --- neuropy/utils/signal_process.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/neuropy/utils/signal_process.py b/neuropy/utils/signal_process.py index 3fdbee85..70d21801 100644 --- a/neuropy/utils/signal_process.py +++ b/neuropy/utils/signal_process.py @@ -128,16 +128,19 @@ def get_pe_spec( event_times = event_times.squeeze() assert event_times.ndim == 1, "event_times must be broadcastable to ndim=1" - # Keep only events within time limits of recording + # Should not be necessary, keep in just in case + # Keep only events within time limits of recording, retaining any nans epoch_start_stops = event_times[:, None] + np.multiply(buffer_sec, (-1, 1)) keep_bool = (epoch_start_stops[:, 0] > self.t_start) & ( epoch_start_stops[:, 1] < self.t_stop ) - # if np.sum(keep_bool) < len(event_times): - # event_times = event_times[keep_bool] - # print( - # f"Events {np.where(~keep_bool)[0]} are outside of data range and were dropped" - # ) + nan_bool = np.where(np.isnan(event_times))[0] + keep_bool = keep_bool | nan_bool + if np.sum(keep_bool) < len(event_times): + event_times = event_times[keep_bool] + print( + f"Events {np.where(~keep_bool)[0]} are outside of data range and were dropped" + ) sxx_list = [] ntime_bins = int(np.ceil(np.sum(buffer_sec) * self.sampling_rate)) From e32f4716ff3591d1dd210cf478a89a0975192ce9 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Thu, 2 Nov 2023 14:05:28 -0400 Subject: [PATCH 06/47] mediocre enhancement --- neuropy/plotting/signals.py | 62 +++++++++++++++++++++---------------- 1 file changed, 36 insertions(+), 26 deletions(-) diff --git a/neuropy/plotting/signals.py b/neuropy/plotting/signals.py index 962ee4e0..79e3b81a 100644 --- a/neuropy/plotting/signals.py +++ b/neuropy/plotting/signals.py @@ -14,6 +14,7 @@ def plot_spectrogram( cmap="jet", sigma=None, std_sxx=None, + widget=True, ): """Generating spectrogram plot for given channel @@ -29,6 +30,8 @@ def plot_spectrogram( overlap between windows, by default 2 ax : [obj], optional if none generates a new figure, by default None + widget: [bool], optional + enable use of sliders for adjust clim and zooming in on a frequency range """ # Figure out if using legacy functionality or updated @@ -72,19 +75,22 @@ def plotspec(n_std, freq_lim): ax.set_ylabel("Frequency (Hz)") # ---- updating plotting values for interaction ------------ - ipywidgets.interact( - plotspec, - n_std=ipywidgets.FloatSlider( - value=6, - min=0.1, - max=30, - step=0.1, - description="Clim :", - ), - freq_lim=ipywidgets.IntRangeSlider( - value=freq_lims, min=0, max=625, step=1, description="Freq. range:" - ), - ) + if widget: + ipywidgets.interact( + plotspec, + n_std=ipywidgets.FloatSlider( + value=6, + min=0.1, + max=30, + step=0.1, + description="Clim :", + ), + freq_lim=ipywidgets.IntRangeSlider( + value=freq_lims, min=0, max=625, step=1, description="Freq. range:" + ), + ) + else: + plotspec(6, freq_lims) else: def plotspec(n_std, freq): @@ -105,19 +111,23 @@ def plotspec(n_std, freq): ax.set_ylabel("Frequency (Hz)") # ---- updating plotting values for interaction ------------ - ipywidgets.interact( - plotspec, - n_std=ipywidgets.FloatSlider( - value=6, - min=0.1, - max=30, - step=0.1, - description="Clim :", - ), - freq=ipywidgets.IntRangeSlider( - value=freq_lims, min=0, max=625, step=1, description="Freq. range:" - ), - ) + if widget: + ipywidgets.interact( + plotspec, + n_std=ipywidgets.FloatSlider( + value=6, + min=0.1, + max=30, + step=0.1, + description="Clim :", + ), + freq=ipywidgets.IntRangeSlider( + value=freq_lims, min=0, max=625, step=1, description="Freq. range:" + ), + ) + else: + plotspec(6, freq_lims) + return ax From 364bd294739f6f03bbbb48c668658c79512b0912 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Thu, 2 Nov 2023 14:06:06 -0400 Subject: [PATCH 07/47] small bugfixes --- neuropy/utils/signal_process.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/neuropy/utils/signal_process.py b/neuropy/utils/signal_process.py index 70d21801..e7d88e71 100644 --- a/neuropy/utils/signal_process.py +++ b/neuropy/utils/signal_process.py @@ -128,13 +128,13 @@ def get_pe_spec( event_times = event_times.squeeze() assert event_times.ndim == 1, "event_times must be broadcastable to ndim=1" - # Should not be necessary, keep in just in case + # Should not be necessary, keep ls - just in case # Keep only events within time limits of recording, retaining any nans epoch_start_stops = event_times[:, None] + np.multiply(buffer_sec, (-1, 1)) keep_bool = (epoch_start_stops[:, 0] > self.t_start) & ( epoch_start_stops[:, 1] < self.t_stop ) - nan_bool = np.where(np.isnan(event_times))[0] + nan_bool = np.isnan(event_times) keep_bool = keep_bool | nan_bool if np.sum(keep_bool) < len(event_times): event_times = event_times[keep_bool] From 95fd4d45c719402da2aafd659689215249495c63 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Fri, 3 Nov 2023 11:38:52 -0400 Subject: [PATCH 08/47] speed up by changing pcolormesh to pcolorfast (note pcolorfast is experimental though) --- neuropy/plotting/signals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neuropy/plotting/signals.py b/neuropy/plotting/signals.py index 79e3b81a..b5652c90 100644 --- a/neuropy/plotting/signals.py +++ b/neuropy/plotting/signals.py @@ -96,7 +96,7 @@ def plotspec(n_std, freq_lim): def plotspec(n_std, freq): """Plots data from Spectrogram class and preserves time and frequency info on axes""" spec_use = spec.time_slice(t_start=time_lims[0], t_stop=time_lims[1]) - ax.pcolormesh( + ax.pcolorfast( spec_use.time, spec_use.freqs, spec_use.traces, From 3597095e01f688f26e8de816b5da898aafef21a1 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Tue, 21 Nov 2023 08:45:19 -0500 Subject: [PATCH 09/47] enh - return t_stop of wav file --- neuropy/io/usvio.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/neuropy/io/usvio.py b/neuropy/io/usvio.py index 06c2aaf2..e6780ca8 100644 --- a/neuropy/io/usvio.py +++ b/neuropy/io/usvio.py @@ -62,8 +62,9 @@ def detect_tone( smooth_window: float = 0.05, thresh: float = 5, tone_length: float = 0.5, - tone_label:str = "start_tone", + tone_label: str = "start_tone", plot_check: bool = True, + return_stop: bool = False, ): """Detect USV start tone - typically a 0.5 sec 500Hz pure tone""" """NRK todo: load in with a downsample. From https://stackoverflow.com/questions/30619740/downsampling-wav-audio-file @@ -117,7 +118,19 @@ def detect_tone( ax.set_xlabel("Time (sec)") ax.set_ylabel("Smoothed power ratio") - return Epoch({"start": tone_cands[:, 0], "stop": tone_cands[:, 1], "label": tone_label}) + if not return_stop: + return Epoch({"start": tone_cands[:, 0], "stop": tone_cands[:, 1], "label": tone_label}) + else: + return Epoch({"start": tone_cands[:, 0], "stop": tone_cands[:, 1], "label": tone_label}), audio_sig.t_stop + + +def load_wav(filename, sr_ds=17750): + """Load in a wav file after downsampling""" + + Audiodata, fs = librosa.load(filename, sr=sr_ds) + audio_sig = Signal(Audiodata, fs) + + return audio_sig if __name__ == "__main__": From 4d9315dfadc216b4ed2c00bf3941072d68ae4c3a Mon Sep 17 00:00:00 2001 From: nkinsky Date: Tue, 21 Nov 2023 09:52:52 -0500 Subject: [PATCH 10/47] enhancement --- neuropy/core/epoch.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index ca6acee8..49d9a7f5 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -430,13 +430,15 @@ def merge_neighbors(self): return Epoch.from_array(epochs_arr[:, 0], epochs_arr[:, 1], labels_arr) - def contains(self, t): + def contains(self, t, return_closest: bool = False): """Check if timepoints lie within epochs, must be non-overlapping epochs Parameters ---------- t : array timepoints in seconds + return_closest: bool + True = return closest epoch before to all points in t even if t is outside epoch Returns ------- @@ -450,11 +452,14 @@ def contains(self, t): bin_loc = np.digitize(t, self.flatten()) indx_bool = bin_loc % 2 == 1 - return ( - indx_bool, - t[indx_bool], - labels[((bin_loc[indx_bool] - 1) / 2).astype("int")], - ) + if not return_closest: + return ( + indx_bool, + t[indx_bool], + labels[((bin_loc[indx_bool] - 1) / 2).astype("int")], + ) + else: + return indx_bool, t, labels[indx_bool] def delete_in_between(self, t1, t2): epochs_df = self.to_dataframe()[["start", "stop", "label"]] From 47a90f20f4b8d921cb3a0df068aa4d8944299511 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Tue, 21 Nov 2023 14:46:11 -0500 Subject: [PATCH 11/47] small enhancement --- neuropy/core/epoch.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index 49d9a7f5..b199cedb 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -459,7 +459,7 @@ def contains(self, t, return_closest: bool = False): labels[((bin_loc[indx_bool] - 1) / 2).astype("int")], ) else: - return indx_bool, t, labels[indx_bool] + return indx_bool, t, labels[bin_loc], bin_loc def delete_in_between(self, t1, t2): epochs_df = self.to_dataframe()[["start", "stop", "label"]] From 9a82b542a4daaa4429743945a7f54af9b4d3b6e4 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 22 Nov 2023 13:58:04 -0500 Subject: [PATCH 12/47] outline future OEsync class --- neuropy/io/openephysio.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/neuropy/io/openephysio.py b/neuropy/io/openephysio.py index 987a9a6a..b8da5750 100644 --- a/neuropy/io/openephysio.py +++ b/neuropy/io/openephysio.py @@ -12,6 +12,21 @@ import matplotlib.pyplot as plt +class OESyncIO: + """Class to synchronize external data sources to Open-Ephys recordings.""" + def __init__(self, basepath) -> None: + pass + + def rough_align_w_datetime(self): + """Perform a rough first-pass alignment of recordings using datetimes""" + pass + + def align_w_TTL(self): + """Align TTLs in OE with external timestamps. + + :return pd.DataFrame OR np.array with OE timestamps matching external TTL events""" + pass + def get_us_start(settings_file: str, from_zone="UTC", to_zone="America/Detroit"): """Get microsecond time second precision start time from Pho/Timestamp plugin""" From 03074a95c372ebcc8a232571a8f4eaf143b6e83c Mon Sep 17 00:00:00 2001 From: Nat Date: Tue, 28 Nov 2023 11:52:23 -0500 Subject: [PATCH 13/47] refactor to rescue missing functions from last merge --- neuropy/core/epoch.py | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index b199cedb..4018fdde 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -607,20 +607,6 @@ def add_epoch_buffer(self, buffer_sec: float or int or tuple or list): self.starts self.stops print(f"Buffer of {buffer_sec} added before/after each epoch") - - -def add_epoch_buffer(epoch_df: pd.DataFrame, buffer_sec: float or int or tuple or list): - """Extend each epoch by buffer_sec before/after start/stop of each epoch""" - if type(buffer_sec) in [int, float]: - buffer_sec = (buffer_sec, buffer_sec) - else: - assert len(buffer_sec) == 2 - - epoch_df["start"] -= buffer_sec[0] - epoch_df["stop"] += buffer_sec[1] - - return epoch_df - @staticmethod def from_peaks(arr: np.ndarray, thresh, length, sep=0, boundary=0, fs=1): hmin, hmax = _unpack_args(thresh) # does not need fs @@ -714,6 +700,17 @@ def get_indices_for_time(self, t: np.array): return time_bool.astype("bool") +def add_epoch_buffer(epoch_df: pd.DataFrame, buffer_sec: float or int or tuple or list): + """Extend each epoch by buffer_sec before/after start/stop of each epoch""" + if type(buffer_sec) in [int, float]: + buffer_sec = (buffer_sec, buffer_sec) + else: + assert len(buffer_sec) == 2 + + epoch_df["start"] -= buffer_sec[0] + epoch_df["stop"] += buffer_sec[1] + + return epoch_df if __name__ == "__main__": art_file = "/data3/Trace_FC/Recording_Rats/Finn2/2023_05_06_habituation1/Finn2_habituation1_denoised.art_epochs.npy" From d1cd37652e79535b4f31028fde9fa4cbeab359b6 Mon Sep 17 00:00:00 2001 From: Nat Date: Tue, 5 Dec 2023 12:59:00 -0500 Subject: [PATCH 14/47] debug code at bottom --- neuropy/utils/signal_process.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/neuropy/utils/signal_process.py b/neuropy/utils/signal_process.py index e7d88e71..dd258a22 100644 --- a/neuropy/utils/signal_process.py +++ b/neuropy/utils/signal_process.py @@ -1436,4 +1436,10 @@ def plot_miniscope_noise( if __name__ == "__main__": - pass + import DataPaths.subjects as subjects + sess = subjects.sd.ratKday1[0] + post = sess.paradigm["post"].flatten() + period = [post[0] + 5 * 3600, post[0] + 6 * 3600] + chan = sess.best_channels.slow_wave + lfp = sess.eegfile.get_signal(chan, *period) + irasa(lfp.traces.squeeze(), lfp.sampling_rate, band=(1, 100), return_fit=True) From 481b023cc7895655f9dfd954c44b18691e2b8129 Mon Sep 17 00:00:00 2001 From: Nat Date: Tue, 5 Dec 2023 13:00:20 -0500 Subject: [PATCH 15/47] nada --- neuropy/analyses/brainstates.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neuropy/analyses/brainstates.py b/neuropy/analyses/brainstates.py index 779d06b1..3ad6e8da 100644 --- a/neuropy/analyses/brainstates.py +++ b/neuropy/analyses/brainstates.py @@ -439,4 +439,4 @@ def plot_thresh(x, fit_params, x_label, low_label, high_label, title): ) print(f"{fp_bokeh_plot} saved") - return epochs + return epochs \ No newline at end of file From 9870159514f216dae92b3511823e6bc7a70a771e Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 8 Dec 2023 14:23:55 -0500 Subject: [PATCH 16/47] bugfix --- neuropy/analyses/decoders.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neuropy/analyses/decoders.py b/neuropy/analyses/decoders.py index 76994e32..227f3705 100644 --- a/neuropy/analyses/decoders.py +++ b/neuropy/analyses/decoders.py @@ -232,7 +232,7 @@ def _estimate(self): """Estimates position within each bin""" tuning_curves = self.ratemap.tuning_curves - bincntr = self.ratemap.x + bincntr = self.ratemap.x_coords() if self.epochs is not None: spkcount, nbins = self.neurons.get_spikes_in_epochs( From 189aa8a1f2df55699d5b28aa0c27120d086c3916 Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 8 Dec 2023 14:24:21 -0500 Subject: [PATCH 17/47] bugfix --- neuropy/core/ratemap.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/neuropy/core/ratemap.py b/neuropy/core/ratemap.py index e1b3b4ed..38e0cf16 100644 --- a/neuropy/core/ratemap.py +++ b/neuropy/core/ratemap.py @@ -138,12 +138,12 @@ def copy(self): @property def x_binsize(self): - return np.diff(self.x_coords)[0] + return np.diff(self.x_coords())[0] @property def y_binsize(self): if self.y is not None: - return np.diff(self.y_coords)[0] + return np.diff(self.y_coords())[0] @property def n_neurons(self): From 9e8ca449183310ef3432193161acf65df0bdaf8e Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 15 Dec 2023 10:32:50 -0500 Subject: [PATCH 18/47] docstring wrapped for readability --- neuropy/analyses/spkepochs.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/neuropy/analyses/spkepochs.py b/neuropy/analyses/spkepochs.py index 5353ad97..65572f30 100644 --- a/neuropy/analyses/spkepochs.py +++ b/neuropy/analyses/spkepochs.py @@ -7,7 +7,8 @@ def detect_off_epochs(mua: core.Mua, ignore_epochs: core.Epoch = None): - """Detects OFF periods using multiunit activity. During these epochs neurons stop almost stop firing. These off periods were reported by Vyazovskiy et al. 2011 in cortex for sleep deprived animals. + """Detects OFF periods using multiunit activity. During these epochs neurons stop almost stop firing. + These off periods were reported by Vyazovskiy et al. 2011 in cortex for sleep deprived animals. Parameters ---------- From 7a2b421db353cd9055ca652fc12df248510b00e9 Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 15 Dec 2023 10:35:18 -0500 Subject: [PATCH 19/47] Finally fixed 'SettingWithCopyWarning in Epoch._validate. My work here is done. --- neuropy/core/epoch.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index 4018fdde..a8092d91 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -5,7 +5,7 @@ from pathlib import Path import scipy.signal as sg import typing -from copy import deepcopy +from copy import deepcopy, copy def _unpack_args(values, fs=1): @@ -47,8 +47,18 @@ def _validate(self, epochs): assert ( pd.Series(["start", "stop", "label"]).isin(epochs.columns).all() ), "epochs should at least have columns/keys with names: start, stop, label" - epochs.loc[:, "label"] = epochs["label"].astype("str") + + ## Make sure labels are formatted correctly as strings. + # Note that the following would be MUCH simpler but throws a "SettingWithCopyWarning" + # so we have to add the convoluted code below to avoid it + # epochs.loc[:, "label"] = epochs.loc[:, "label"].astype("str") + epochs_labels_str = copy(epochs["label"].astype("str")) + epochs = epochs.drop(columns="label", inplace=False) # this also throws a warning if used with inplace=True + epochs.loc[:, "label"] = epochs_labels_str + + # Sort epochs = epochs.sort_values(by=["start"]).reset_index(drop=True) + return epochs.copy() @property From a2179f000ec12484d914219681ad7b265f11c3da Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 15 Dec 2023 10:38:38 -0500 Subject: [PATCH 20/47] enh/bugfix to id gaps in epochs as zeros --- neuropy/core/epoch.py | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index a8092d91..556c541c 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -505,7 +505,7 @@ def delete_in_between(self, t1, t2): epochs_df = pd.concat([epochs_df, flank_start, flank_stop], ignore_index=True) return Epoch(epochs_df) - def proportion_by_label(self, t_start=None, t_stop=None): + def proportion_by_label(self, t_start=None, t_stop=None, ignore_gaps=False): """Get porportion of time for each label type Parameters @@ -529,25 +529,29 @@ def proportion_by_label(self, t_start=None, t_stop=None): ep = self._epochs.copy() ep = ep[(ep.stop > t_start) & (ep.start < t_stop)].reset_index(drop=True) + if not ignore_gaps: + assert ep.shape[0] > 0, "cannot have empty time gaps between epoch labels with ignore_gaps=False" + elif ignore_gaps and (ep.shape[0] > 0): + if ep["start"].iloc[0] < t_start: + ep.at[0, "start"] = t_start - if ep["start"].iloc[0] < t_start: - ep.at[0, "start"] = t_start + if ep["stop"].iloc[-1] > t_stop: + ep.at[ep.index[-1], "stop"] = t_stop - if ep["stop"].iloc[-1] > t_stop: - ep.at[ep.index[-1], "stop"] = t_stop + ep["duration"] = ep.stop - ep.start - ep["duration"] = ep.stop - ep.start + ep_group = ep.groupby("label").sum(numeric_only=True).duration / duration - ep_group = ep.groupby("label").sum().duration / duration + label_proportion = {} + for label in self.get_unique_labels(): + label_proportion[label] = 0.0 - label_proportion = {} - for label in self.get_unique_labels(): - label_proportion[label] = 0.0 + for state in ep_group.index.values: + label_proportion[state] = ep_group[state] - for state in ep_group.index.values: - label_proportion[state] = ep_group[state] - - return label_proportion + return label_proportion + else: + return None def durations_by_label(self): """Return total duration for each unique label From de0ef52b951b2f7c5b21838d2c9cc19ae26bc0de Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 15 Dec 2023 10:39:05 -0500 Subject: [PATCH 21/47] fix typos and docstring --- neuropy/core/epoch.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index 556c541c..31774f7f 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -506,7 +506,7 @@ def delete_in_between(self, t1, t2): return Epoch(epochs_df) def proportion_by_label(self, t_start=None, t_stop=None, ignore_gaps=False): - """Get porportion of time for each label type + """Get proportion of time for each label type Parameters ---------- @@ -514,6 +514,7 @@ def proportion_by_label(self, t_start=None, t_stop=None, ignore_gaps=False): start time in seconds, by default None t_stop : float, optional stop time in seconds, by default None + ignore_gaps: will return None if set and there is no epoch in the time period selected. Returns ------- @@ -559,7 +560,7 @@ def durations_by_label(self): Returns ------- dict - dictionary contating duration of each unique label + dictionary containing duration of each unique label """ labels = self.labels durations = self.durations From c710a47c77bbdaf02eded945cf1884e701d2182f Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 15 Dec 2023 10:39:30 -0500 Subject: [PATCH 22/47] add resample_labeled_epochs method --- neuropy/core/epoch.py | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index 31774f7f..4b7ea977 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -571,6 +571,43 @@ def durations_by_label(self): return label_durations + def resample_labeled_epochs(self, res, t_start=None, t_stop=None, merge_neighbors=True): + """Resample epochs to different size blocks using a winner take all method to assign + a label name. e.g. if the first 100-second epoch is 40% quiet wake, 50% REM, and 10% NREM + it would get labeled as REM. Pretty slow, even slower with merge_neighbors=True + + :param: res: block size in seconds + :param: t_start: start time in seconds, default = start of first epoch + :param: t_stop : stop time in seconds, default = stop of last epoch + :param merge_neighbors: combine adjacent epochs of the same label, default=True""" + + if t_start is None: + t_start = self.starts[0] + elif t_start < self.starts[0]: + t_start = self.starts[0] + print('t_start < start time of first epoch, reassigned to match first epoch start time') + + if t_stop is None: + t_stop = self.stops[-1] + if t_stop > self.stops[-1]: + t_stop = self.stops[-1] + print('t_stop > stop time of first epoch, reassigned to match last epoch stop time') + bins = np.arange(t_start, t_stop + res, res) + start_rs = bins[:-1] + stop_rs = bins[1:] + label_rs = [] + for start, stop in zip(start_rs, stop_rs): + props = self.proportion_by_label(start, stop, ignore_gaps=True) + label_add = list(props.keys())[np.argmax(list(props.values()))] if props is not None else "" + label_rs.append(label_add) + # except AssertionError: # Append nothing if gap found in epochs + # label_rs.append("") + + epoch_rs = Epoch(pd.DataFrame({"start": start_rs, "stop": stop_rs, "label": label_rs})) + epoch_rs = epoch_rs.merge_neighbors() if merge_neighbors else epoch_rs + + return epoch_rs + def count(self, t_start=None, t_stop=None, binsize=300): if t_start is None: t_start = 0 From 2b9adc6c7fded7c4915027996732874690dfaed9 Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 15 Dec 2023 16:46:21 -0500 Subject: [PATCH 23/47] bugfix --- neuropy/plotting/epochs.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/neuropy/plotting/epochs.py b/neuropy/plotting/epochs.py index 5b2b5f8c..a6f12842 100644 --- a/neuropy/plotting/epochs.py +++ b/neuropy/plotting/epochs.py @@ -30,7 +30,6 @@ def plot_epochs( [type] [description] """ - if isinstance(epochs, pd.DataFrame): epochs = Epoch(epochs) @@ -47,9 +46,10 @@ def plot_epochs( elif isinstance(colors, dict): colors = [colors[label] for label in epochs.labels] - if epochs.has_labels: + if epochs.has_labels or (len(epochs.to_dataframe().label) > 1): labels = epochs.labels - unique_labels = np.unique(epochs.labels) + # unique_labels = np.unique(epochs.labels) + unique_labels = epochs.to_dataframe().label.unique() n_labels = len(unique_labels) if labels_order is not None: From 9c9480182610b20f8b3636efd7094de6cdcd61b1 Mon Sep 17 00:00:00 2001 From: Nat Date: Tue, 19 Dec 2023 09:10:49 -0500 Subject: [PATCH 24/47] enhcancement - color by key/value other than 'label' --- neuropy/plotting/epochs.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/neuropy/plotting/epochs.py b/neuropy/plotting/epochs.py index a6f12842..27a84cc1 100644 --- a/neuropy/plotting/epochs.py +++ b/neuropy/plotting/epochs.py @@ -7,7 +7,7 @@ def plot_epochs( - epochs: Epoch, labels_order=None, colors="Set3", alpha=1, collapsed=False, ax=None + epochs: Epoch, labels_order=None, colors="Set3", alpha=1, collapsed=False, colorby="label", ax=None ): """Plots epochs on a given axis, with different style of plotting @@ -21,8 +21,10 @@ def plot_epochs( [description], by default 0.5 ymax : float, optional [description], by default 0.55 - color : str, optional - [description], by default "gray" + color : str or dict, optional + [description], by default "gray", if dict = {"value1": color, "value2": color2} + where value1, value2, ... are values in the column defined by colorby param + colorby: str, column in epochs to map colors to collapsed: Returns @@ -33,7 +35,7 @@ def plot_epochs( if isinstance(epochs, pd.DataFrame): epochs = Epoch(epochs) - assert isinstance(epochs, Epoch), "epochs must be neuropy.Epoch object" + # assert isinstance(epochs, Epoch), "epochs must be neuropy.Epoch object" n_epochs = epochs.n_epochs @@ -44,7 +46,7 @@ def plot_epochs( except: colors = [colors] * n_epochs elif isinstance(colors, dict): - colors = [colors[label] for label in epochs.labels] + colors = [colors[label] for label in epochs.to_dataframe()[colorby]] if epochs.has_labels or (len(epochs.to_dataframe().label) > 1): labels = epochs.labels From 2d8380f61850ecf3bc674f84b8b62137e1b63a48 Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 22 Dec 2023 17:28:06 -0500 Subject: [PATCH 25/47] bugfix --- neuropy/plotting/epochs.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/neuropy/plotting/epochs.py b/neuropy/plotting/epochs.py index 27a84cc1..72dc471c 100644 --- a/neuropy/plotting/epochs.py +++ b/neuropy/plotting/epochs.py @@ -46,7 +46,8 @@ def plot_epochs( except: colors = [colors] * n_epochs elif isinstance(colors, dict): - colors = [colors[label] for label in epochs.to_dataframe()[colorby]] + # Define colors, sending those not in the colors dict to white + colors = [colors[label] if label in colors.keys() else "#ffffff" for label in epochs.to_dataframe()[colorby]] if epochs.has_labels or (len(epochs.to_dataframe().label) > 1): labels = epochs.labels From 02dbe3fa843c7e519c63a04fc0c51560d4b50708 Mon Sep 17 00:00:00 2001 From: Nat Date: Mon, 1 Jan 2024 14:54:04 -0500 Subject: [PATCH 26/47] bugfix --- neuropy/plotting/epochs.py | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/neuropy/plotting/epochs.py b/neuropy/plotting/epochs.py index 72dc471c..50a2237d 100644 --- a/neuropy/plotting/epochs.py +++ b/neuropy/plotting/epochs.py @@ -55,11 +55,22 @@ def plot_epochs( unique_labels = epochs.to_dataframe().label.unique() n_labels = len(unique_labels) + # Update to order labels correctly if labels_order is not None: - assert np.array_equal( - np.sort(labels_order), np.sort(unique_labels) - ), "labels_order does not match with epochs labels" - unique_labels = labels_order + # assert np.array_equal( + # np.sort(labels_order), np.sort(unique_labels) + # ), "labels_order does not match with epochs labels" + + # Make sure all labels are in labels_order + + # This code might be necessary, keep for potential debugging + # if np.array_equal(np.sort(labels_order), np.sort(unique_labels)) or \ + # np.all([label in labels_order for label in unique_labels]): + if np.all([label in labels_order for label in unique_labels]): + unique_labels = labels_order + n_labels = len(unique_labels) + else: + assert False, "labels_order does not match with epochs labels" dh = 1 if collapsed else 1 / n_labels y_min = np.zeros(len(epochs)) @@ -75,11 +86,12 @@ def plot_epochs( epoch.start, epoch.stop, y_min[i], - y_min[i] + dh, + y_min[ i] + dh, facecolor=colors[i], edgecolor=None, alpha=alpha, ) + ax.set_ylim([0, 1]) return ax From a3cc7732534b00c8a6583bc851fc7af9a073ffe8 Mon Sep 17 00:00:00 2001 From: Nat Date: Thu, 4 Jan 2024 11:59:39 -0500 Subject: [PATCH 27/47] bugfix --- neuropy/core/epoch.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index 4b7ea977..e18cf717 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -402,7 +402,7 @@ def merge(self, dt): return Epoch.from_array(epochs_arr[:, 0], epochs_arr[:, 1]) - def merge_neighbors(self): + def merge_neighbors(self, max_epoch_sep=1e-6): """Epochs of same label and common boundary will be merged. For example, [1,2] and [2,3] --> [1,3] @@ -411,7 +411,7 @@ def merge_neighbors(self): core.Epoch epochs after merging neigbours sharing same label and boundary """ - ep_times, ep_stops, ep_labels = (self.starts, self.stops, self.labels) + ep_times, ep_stops, ep_labels = (deepcopy(self.starts), deepcopy(self.stops), deepcopy(self.labels)) ep_durations = self.durations ind_delete = [] @@ -420,7 +420,7 @@ def merge_neighbors(self): for i in range(len(inds) - 1): # if two sequentially adjacent epochs with the same label # overlap or have less than 1 microsecond separation, merge them - if ep_times[inds[i + 1]] - ep_stops[inds[i]] < 1e-6: + if ((inds[i+1] - inds[i]) == 1) & ((ep_times[inds[i + 1]] - ep_stops[inds[i]]) < max_epoch_sep): # stretch the second epoch to cover the range of both epochs ep_times[inds[i + 1]] = min( ep_times[inds[i]], ep_times[inds[i + 1]] From e2553d15210cfc3d817eacceb1cb3a3307dca660 Mon Sep 17 00:00:00 2001 From: Nat Date: Thu, 4 Jan 2024 17:06:30 -0500 Subject: [PATCH 28/47] enh to compute overlap time of two different Epoch objects --- neuropy/core/epoch.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index e18cf717..ea5f6af9 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -764,6 +764,24 @@ def add_epoch_buffer(epoch_df: pd.DataFrame, buffer_sec: float or int or tuple o return epoch_df + +def get_epoch_overlap_duration(epochs1: Epoch, epochs2: Epoch): + """Calculate time of overlapping epochs""" + e_array1 = epochs1.to_dataframe().loc[:, ["start", "stop"]].values + e_array2 = epochs2.to_dataframe().loc[:, ["start", "stop"]].values + overlaps = [] + for e1 in e_array1: + for e2 in e_array2: + overlaps.append(getOverlap(e1, e2)) + + return np.array(overlaps).sum() + + +def getOverlap(a, b): + """From https://stackoverflow.com/questions/2953967/built-in-function-for-computing-overlap-in-python""" + return max(0, min(a[1], b[1]) - max(a[0], b[0])) + + if __name__ == "__main__": art_file = "/data3/Trace_FC/Recording_Rats/Finn2/2023_05_06_habituation1/Finn2_habituation1_denoised.art_epochs.npy" art_epochs = Epoch(epochs=None, file=art_file) From 4784b4d8e76da4d0888f095e8ac572b811fc363e Mon Sep 17 00:00:00 2001 From: Nat Date: Tue, 9 Jan 2024 16:21:07 -0500 Subject: [PATCH 29/47] type fix --- neuropy/core/epoch.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index ea5f6af9..646c0b51 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -409,7 +409,7 @@ def merge_neighbors(self, max_epoch_sep=1e-6): Returns ------- core.Epoch - epochs after merging neigbours sharing same label and boundary + epochs after merging neighbours sharing same label and boundary """ ep_times, ep_stops, ep_labels = (deepcopy(self.starts), deepcopy(self.stops), deepcopy(self.labels)) ep_durations = self.durations From 8b317f520fac0b264892d35e371276cd6bd3ffef Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 26 Jan 2024 09:06:54 -0500 Subject: [PATCH 30/47] allow slicing of epochs by list of labels --- neuropy/core/epoch.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index 646c0b51..d372f73d 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -163,6 +163,9 @@ def __getitem__(self, i): data = self._epochs[self._epochs["label"] == i].copy() # elif all(isinstance(_, str) for _ in i): # data = self._epochs[self._epochs["label"].isin(i)].copy() + elif isinstance(i, list): + assert all(isinstance(_, str) for _ in i), "All entries in epochs slicing list must be str" + data = self._epochs[[label in i for label in self._epochs["label"]]] elif isinstance(i, (int, np.integer)): data = self._epochs.iloc[[i]].copy() else: From 595977306e2c272c5e7abb9f1bd5067c1a38ecc3 Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 26 Jan 2024 09:07:22 -0500 Subject: [PATCH 31/47] add less agressive flatten funciton (only flattens one level) --- neuropy/utils/misc.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/neuropy/utils/misc.py b/neuropy/utils/misc.py index 9a4042d1..4723b831 100644 --- a/neuropy/utils/misc.py +++ b/neuropy/utils/misc.py @@ -1,6 +1,7 @@ import numpy as np import math from collections.abc import Iterable +from itertools import chain def find_nearest(array, value): @@ -34,10 +35,16 @@ def get_interval(period, nwindows): return interval -def flatten(xs): +def flatten(list_in): + """Flatten a ragged list of different sized lists into one continuous list. + Unlike `flatten_all` this only flattens the top level.""" + return list(chain.from_iterable(list_in)) + +def flatten_all(xs): + """Completely flattens an iterable of iterables into one long generator""" for x in xs: if isinstance(x, Iterable) and not isinstance(x, (str, bytes)): - yield from flatten(x) + yield from flatten_all(x) else: yield x From 0a63936e1a23fbcbe68f689d70e9512d5d754586 Mon Sep 17 00:00:00 2001 From: Nat Date: Wed, 20 Mar 2024 12:58:54 -0400 Subject: [PATCH 32/47] add get_speed() and update get_timestamps() to grab from miniscope timestamp files --- neuropy/io/dlcio.py | 74 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 56 insertions(+), 18 deletions(-) diff --git a/neuropy/io/dlcio.py b/neuropy/io/dlcio.py index de664dda..85556b5b 100644 --- a/neuropy/io/dlcio.py +++ b/neuropy/io/dlcio.py @@ -12,6 +12,8 @@ from pickle import dump, load import scipy.ndimage as simage +from neuropy.io.miniscopeio import MiniscopeIO + from .optitrackio import getStartTime, posfromCSV @@ -76,6 +78,9 @@ def __init__(self, base_path, search_str=None, movie_type=".mp4", pix2cm=1): # Grab metadata for later access self.get_metadata() + # Initialize other fields + self.timestamp_type = None + @property def SampleRate(self): try: @@ -114,27 +119,37 @@ def get_metadata(self): with open(meta_file, "rb") as f: self.meta = load(f) - def get_timestamps(self): + def get_timestamps(self, camera_type: str in ['optitrack', 'ms_webcam', 'ms_webcam1', 'ms_webcam2'] = 'optitrack'): """Tries to import timestamps from CSV file from optitrack, if not, infers it from timestamp in filename, - sample rate, and nframes""" - - opti_file = self.tracking_file.parent / ( - self.tracking_file.stem[: self.tracking_file.stem.find("-Camera")] + ".csv" - ) - if opti_file.is_file(): - start_time = getStartTime(opti_file) - _, _, _, t = posfromCSV(opti_file) + sample rate, and nframes + :param camera_type: 'optitrack' looks for optitrack csv file with tracking data, other options look for + UCLA miniscope webcam files""" + + assert camera_type in ['optitrack', 'ms_webcam', 'ms_webcam1', 'ms_webcam2'] + self.timestamp_type = camera_type + if camera_type == 'optitrack': + opti_file = self.tracking_file.parent / ( + self.tracking_file.stem[: self.tracking_file.stem.find("-Camera")] + ".csv" + ) + if opti_file.is_file(): + start_time = getStartTime(opti_file) + _, _, _, t = posfromCSV(opti_file) - self.timestamps = start_time + pd.to_timedelta(t, unit="sec") + self.timestamps = start_time + pd.to_timedelta(t, unit="sec") + else: + print( + "No Optitrack csv file found, inferring start time from file name. SECOND PRECISION IN START TIME!!!" + ) + start_time = deduce_starttime_from_file(self.tracking_file) + time_deltas = pd.to_timedelta( + np.arange(self.nframes) / self.SampleRate, unit="sec" + ) + self.timestamps = start_time + time_deltas else: - print( - "No Optitrack csv file found, inferring start time from file name. SECOND PRECISION IN START TIME!!!" - ) - start_time = deduce_starttime_from_file(self.tracking_file) - time_deltas = pd.to_timedelta( - np.arange(self.nframes) / self.SampleRate, unit="sec" - ) - self.timestamps = start_time + time_deltas + mio = MiniscopeIO(self.base_dir) + webcam_flag = True if camera_type == "ms_webcam" else int(camera_type[-1]) + self.timestamps = mio.load_all_timestamps(webcam=webcam_flag) + def to_cm(self): """Convert pixels to centimeters in pos_data""" @@ -176,6 +191,29 @@ def smooth_pos( SampleRate=self.SampleRate, ) + def get_speed(self, bodypart): + """Get speed of animal""" + assert self.pos_smooth is not None, "You must smooth data first" + assert self.timestamps is not None, "You must get timestamps first" + if self.timestamp_type == "optitrack": + print('get_speed not yet tested for optitrack data, use with caution') + + # Get position of bodypart + data_use = self.pos_smooth[bodypart] + x = data_use["x"] + y = data_use["y"] + total_seconds = (self.timestamps['Timestamps'] - self.timestamps['Timestamps'][0]).dt.total_seconds() + + # Get delta in position and time from frame-to-frame + delta_pos = np.sqrt(np.square(np.diff(x)) + np.square(np.diff(y))) + delta_t = np.diff(total_seconds) + + # Calculate speed + speed = delta_pos / delta_t + + return speed + + def plot1d( self, bodyparts=None, From 60f37d75d94a2b57376282c0769768b0273bc46a Mon Sep 17 00:00:00 2001 From: Nat Date: Wed, 20 Mar 2024 12:59:37 -0400 Subject: [PATCH 33/47] add ability to get WebCam timestamps --- neuropy/io/miniscopeio.py | 29 ++++++++++++++++++++++------- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/neuropy/io/miniscopeio.py b/neuropy/io/miniscopeio.py index af00d049..9e902846 100644 --- a/neuropy/io/miniscopeio.py +++ b/neuropy/io/miniscopeio.py @@ -16,16 +16,18 @@ def __init__(self, basedir) -> None: self.basedir = Path(basedir) self.times_all = None self.orient_all = None + self.webcam_times_all = {} pass def load_all_timestamps( - self, format="UCLA", exclude_str: str = "WebCam", print_included_folders=False + self, format="UCLA", webcam: bool or int = False, exclude_str: str = "WebCam", print_included_folders=False ): """Loads multiple timestamps from multiple videos in the UCLA miniscope software file format (folder = ""hh_mm_ss") :param format: str, either 'UCLA' (default) to use the UCLA miniscope folder format or a regex if you are using a different folder naming convention. + :param webcam see "load_timestamps" function :param exclude_str: exclude any folders containing this string from being loaded in. :param print_included_folders: bool, True = prints folders included in generating timestamps, useful for debugging mismatches in nframes, default=False @@ -61,13 +63,19 @@ def load_all_timestamps( times_list = [] for rec_folder in self.rec_folders: times_temp, _, _, _ = load_timestamps( - rec_folder, corrupted_videos="from_file" + rec_folder, webcam=webcam, corrupted_videos="from_file" ) times_list.append(times_temp) - self.times_all = pd.concat(times_list) + if not webcam: + self.times_all = pd.concat(times_list) + + return self.times_all + else: + self.webcam_number = webcam + self.webcam_times_all = pd.concat(times_list) - return self.times_all + return self.webcam_times_all def load_all_orientation(self, format="UCLA", exclude_str: str = "WebCam"): """Loads head orientation data from multiple videos in the UCLA miniscope @@ -131,11 +139,12 @@ def get_recording_metadata(rec_folder: pathlib.Path): def load_timestamps( - rec_folder, corrupted_videos=None, print_success=False, print_corrupt_success=True + rec_folder, webcam: False or True or int = False, corrupted_videos=None, print_success=False, print_corrupt_success=True ): """Loads in timestamps corresponding to all videos in rec_folder. :param rec_folder: str or path + :param webcam: if True load "My_WebCam", if int load that # Webcam, e.g. 1 = 'My_WebCam1' :param corrupted_videos: (default) list of indices of corrupted files (e.g. [4, 6] would mean the 5th and 7th videos were corrupted and should be skipped - their timestamps will be omitted from the output. 'from_file' will automatically grab any video indices provided as a csv file named 'corrupted_videos.csv' @@ -151,6 +160,10 @@ def load_timestamps( # Grab metadata rec_metadata, vid_metadata, vid_folder = get_recording_metadata(rec_folder) + if webcam: + folder_name = Path("My_WebCam" if isinstance(webcam, bool) else f"My_WebCam{webcam}") + vid_folder = vid_folder.parent / folder_name + # Derive start_time from rec_metadata rec_start = rec_metadata["recordingStartTime"] del rec_start["msecSinceEpoch"] @@ -332,5 +345,7 @@ def euler_from_quaternion(qx, qy, qz, qw): if __name__ == "__main__": - parent_folder = "/data2/Trace_FC/Recording_Rats/Rose/2022_06_20_habituation1" - move_files_to_combined_folder(parent_folder) + dir_use = '/data/Nat/Trace_FC/Recording_Rats/Django/2023_03_08_training' + mio_dj = MiniscopeIO(dir_use) + mio_dj.load_all_timestamps(webcam=True, exclude_str="My_V4_Miniscope") + From ca8cefd722795bf40e36982fdf43169722d516ff Mon Sep 17 00:00:00 2001 From: Nat Date: Wed, 20 Mar 2024 13:00:02 -0400 Subject: [PATCH 34/47] small commenting update --- neuropy/core/session.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neuropy/core/session.py b/neuropy/core/session.py index 111d2610..6bdae240 100644 --- a/neuropy/core/session.py +++ b/neuropy/core/session.py @@ -9,7 +9,7 @@ def __init__(self, basepath): basepath = Path(basepath) self.basepath = basepath xml_files = sorted(basepath.glob("*.xml")) - assert len(xml_files) == 1, "Found more than one .xml file" + assert len(xml_files) == 1, "Found fewer/more than one .xml file" fp = xml_files[0].with_suffix("") self.filePrefix = fp From 06774893263ddaab9a6987998ae622384336bbde Mon Sep 17 00:00:00 2001 From: Nat Date: Wed, 20 Mar 2024 13:01:15 -0400 Subject: [PATCH 35/47] option to return ripple power in detect_ripple_epochs --- neuropy/analyses/oscillations.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/neuropy/analyses/oscillations.py b/neuropy/analyses/oscillations.py index af0829c8..0c06ea8c 100644 --- a/neuropy/analyses/oscillations.py +++ b/neuropy/analyses/oscillations.py @@ -237,6 +237,7 @@ def detect_ripple_epochs( sigma=0.0125, ignore_epochs: Epoch = None, ripple_channel: int or list = None, + return_power: bool = False, ): # TODO chewing artifact frequency (>300 Hz) or emg based rejection of ripple epochs @@ -291,10 +292,18 @@ def detect_ripple_epochs( fs=signal.sampling_rate, sigma=sigma, ignore_times=ignore_times, + return_power=return_power, ) - epochs = epochs.shift(dt=signal.t_start) - epochs.metadata = dict(channels=selected_chans) - return epochs + + if not return_power: + epochs = epochs.shift(dt=signal.t_start) + epochs.metadata = dict(channels=selected_chans) + return epochs + else: + ripple_power = epochs[1] + epochs = epochs[0].shift(dt=signal.t_start) + epochs.metadata = dict(channels=selected_chans) + return epochs, ripple_power def detect_sharpwave_epochs( From 4489e3b27ef7f05d5b69eaee26db3bc40c3acc62 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 20 Mar 2024 13:04:55 -0400 Subject: [PATCH 36/47] add detect_beta_epochs --- neuropy/analyses/oscillations.py | 83 +++++++++++++++++++++++++++++--- 1 file changed, 77 insertions(+), 6 deletions(-) diff --git a/neuropy/analyses/oscillations.py b/neuropy/analyses/oscillations.py index 0c06ea8c..4b756b50 100644 --- a/neuropy/analyses/oscillations.py +++ b/neuropy/analyses/oscillations.py @@ -123,7 +123,7 @@ def _detect_freq_band_epochs( def detect_hpc_delta_wave_epochs( signal: Signal, - freq_band=(0.5, 4), + freq_band=(0.2, 5), min_dur=0.15, max_dur=0.5, ignore_epochs: Epoch = None, @@ -225,6 +225,72 @@ def detect_hpc_delta_wave_epochs( return Epoch(epochs=epochs, metadata=params) +def detect_beta_epochs( + signal: Signal, + probegroup: ProbeGroup = None, + freq_band=(15, 40), + thresh=(0, 0.5), + mindur=0.25, + maxdur=5, + mergedist=0.5, + sigma=0.125, + edge_cutoff=-0.25, + ignore_epochs: Epoch = None, + return_power=False, +): + if probegroup is None: + selected_chan = signal.channel_id + traces = signal.traces + else: + if isinstance(probegroup, np.ndarray): + changrps = np.array(probegroup, dtype="object") + if isinstance(probegroup, ProbeGroup): + changrps = probegroup.get_connected_channels(groupby="shank") + channel_ids = np.concatenate(changrps).astype("int") + + duration = signal.duration + t1, t2 = signal.t_start, signal.t_start + np.min([duration, 3600]) + signal_slice = signal.time_slice(channel_id=channel_ids, t_start=t1, t_stop=t2) + hil_stat = signal_process.hilbert_amplitude_stat( + signal_slice.traces, + freq_band=freq_band, + fs=signal.sampling_rate, + statistic="mean", + ) + selected_chan = channel_ids[np.argmax(hil_stat)].reshape(-1) + traces = signal.time_slice(channel_id=selected_chan).traces.reshape(1, -1) + + print(f"Best channel for beta: {selected_chan}") + if ignore_epochs is not None: + ignore_times = ignore_epochs.as_array() + else: + ignore_times = None + + epochs = _detect_freq_band_epochs( + signals=traces, + freq_band=freq_band, + thresh=thresh, + mindur=mindur, + maxdur=maxdur, + mergedist=mergedist, + fs=signal.sampling_rate, + ignore_times=ignore_times, + sigma=sigma, + edge_cutoff=edge_cutoff, + return_power=return_power, + ) + + if not return_power: + epochs = epochs.shift(dt=signal.t_start) + epochs.metadata = dict(channels=selected_chan) + return epochs + else: + beta_power = epochs[1] + epochs = epochs[0].shift(dt=signal.t_start) + epochs.metadata = dict(channels=selected_chan) + return epochs, beta_power + + def detect_ripple_epochs( signal: Signal, probegroup: ProbeGroup = None, @@ -514,13 +580,16 @@ class Gamma: def get_peak_intervals( self, lfp, - band=(25, 50), + band=(40, 80), lowthresh=0, highthresh=1, minDistance=300, minDuration=125, + return_amplitude=False, + ): - """Returns strong theta lfp. If it has multiple channels, then strong theta periods are calculated from that channel which has highest area under the curve in the theta frequency band. Parameters are applied on z-scored lfp. + """Returns strong theta lfp. If it has multiple channels, then strong theta periods are calculated from that + channel which has highest area under the curve in the theta frequency band. Parameters are applied on z-scored lfp. Parameters ---------- @@ -554,10 +623,12 @@ def get_peak_intervals( minDistance=minDistance, minDuration=minDuration, ) + if not return_amplitude: + return peakevents + else: + return peakevents, gamma_amp - return peakevents - - def csd(self, period, refchan, chans, band=(25, 50), window=1250): + def csd(self, period, refchan, chans, band=(40, 80), window=1250): """Calculating current source density using laplacian method Parameters From 75072769874a0c7056f4c2f7fa0ea9953e7c951d Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 20 Mar 2024 13:05:27 -0400 Subject: [PATCH 37/47] add crucial assertion --- neuropy/core/epoch.py | 1 + 1 file changed, 1 insertion(+) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index d372f73d..7f578c15 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -460,6 +460,7 @@ def contains(self, t, return_closest: bool = False): """ assert self.is_overlapping == False, "Epochs must be non overlapping" + assert isinstance(t, np.ndarray), "t must be a numpy.ndarray" labels = self.labels bin_loc = np.digitize(t, self.flatten()) From 79a3349f2f9c315c2cc9d98760b836d6eb181b03 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 20 Mar 2024 13:11:07 -0400 Subject: [PATCH 38/47] merge --- neuropy/io/miniscopeio.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/neuropy/io/miniscopeio.py b/neuropy/io/miniscopeio.py index 9e902846..3b81ad28 100644 --- a/neuropy/io/miniscopeio.py +++ b/neuropy/io/miniscopeio.py @@ -266,6 +266,7 @@ def load_orientation(rec_folder, corrupted_videos=None): def move_files_to_combined_folder( parent_folder, re_pattern="**/My_V4_Miniscope/*.avi", + prepend_rec_date=False, prepend_rec_time=True, copy_during_prepend=False, ): @@ -295,10 +296,19 @@ def move_files_to_combined_folder( [file.parent for file in movie_files] ) # get unique folder names for folder in folder_names: + + # Prepend time of day prepend_time_from_folder_to_file( folder, ext=re_pattern.split(".")[-1], copy=copy_during_prepend ) + # Prepend recording date + if prepend_rec_date: + prepend_time_from_folder_to_file( + folder, ext=re_pattern.split(".")[-1], copy=False, + time_str="[0-9]{4}_[0-9]{2}_[0-9]{2}$", + ) + # Finally, get updated file names movie_files = sorted(parent_folder.glob(re_pattern)) @@ -345,7 +355,6 @@ def euler_from_quaternion(qx, qy, qz, qw): if __name__ == "__main__": - dir_use = '/data/Nat/Trace_FC/Recording_Rats/Django/2023_03_08_training' - mio_dj = MiniscopeIO(dir_use) - mio_dj.load_all_timestamps(webcam=True, exclude_str="My_V4_Miniscope") + parent_folder = '/data2/Psilocybin/Recording_Rats/Finn2' + move_files_to_combined_folder(parent_folder, prepend_rec_date=True, copy_during_prepend=True) From aa5fb61b7f4814063012c1dd70ab9498cffc46c4 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 20 Mar 2024 13:06:47 -0400 Subject: [PATCH 39/47] bugfix --- neuropy/plotting/signals.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/neuropy/plotting/signals.py b/neuropy/plotting/signals.py index b5652c90..04a73abe 100644 --- a/neuropy/plotting/signals.py +++ b/neuropy/plotting/signals.py @@ -48,6 +48,7 @@ def plot_spectrogram( sxx = gaussian_filter(sxx, sigma=sigma) if std_sxx is None: # Calculate standard deviation if needed for plotting purposes. std_sxx = np.std(sxx) + mean_sxx = np.mean(sxx) # time = np.linspace(time[0], time[1], sxx.shape[1]) # freq = np.linspace(freq[0], freq[1], sxx.shape[0]) @@ -56,13 +57,13 @@ def plot_spectrogram( # ---------- plotting ---------------- if legacy: - + print('test') def plotspec(n_std, freq_lim): """Plots data fine but doesn't preserve time and frequency info on axes""" ax.imshow( sxx, cmap=cmap, - vmax=n_std * std_sxx, + vmax= mean_sxx + n_std * std_sxx, rasterized=True, origin="lower", extent=[time_lims[0], time_lims[-1], freq_lims[0], freq_lims[-1]], @@ -76,6 +77,7 @@ def plotspec(n_std, freq_lim): # ---- updating plotting values for interaction ------------ if widget: + print('test2') ipywidgets.interact( plotspec, n_std=ipywidgets.FloatSlider( @@ -101,7 +103,7 @@ def plotspec(n_std, freq): spec_use.freqs, spec_use.traces, cmap=cmap, - vmax=n_std * std_sxx, + vmax=mean_sxx + n_std * std_sxx, rasterized=True, ) ax.set_ylim(freq) From 379f625fe5db0a9f09ab1922cba5c97a6dec8739 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 20 Mar 2024 13:07:46 -0400 Subject: [PATCH 40/47] enh to prepend date to files when moving --- neuropy/utils/manipulate_files.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/neuropy/utils/manipulate_files.py b/neuropy/utils/manipulate_files.py index 5200fd73..12b0a44c 100644 --- a/neuropy/utils/manipulate_files.py +++ b/neuropy/utils/manipulate_files.py @@ -24,7 +24,7 @@ def get_record_time_from_pathname( ext = "*." + ext dir = Path(folder) files = sorted(dir.glob("**/" + ext)) - folder_dates = [] + folder_times = [] # Iterate through all parent folders and find the one matching your time_str for file in files: @@ -33,18 +33,19 @@ def get_record_time_from_pathname( match_obj = re.match(time_str, fold) if match_obj is not None: dir_date.append(match_obj.group(0)) + break # Except for loop for that folder now that you found its time! try: assert ( len(dir_date) == 1 ), "time_str you entered returns too many or no matches" except AssertionError: pass - folder_dates.extend(dir_date) + folder_times.extend(dir_date) - if len(files) != len(folder_dates): - print("Obtained fewer dates than files, double-check") + if len(files) != len(folder_times): + print("Obtained fewer times than files, double-check") - return files, folder_dates + return files, folder_times def prepend_time_from_folder_to_file( @@ -58,13 +59,13 @@ def prepend_time_from_folder_to_file( ## NRK todo: add in user prompt if copy=True entered!!! # Get names to prepend - files, folder_dates = get_record_time_from_pathname(folder, ext, time_str) + files, folder_times = get_record_time_from_pathname(folder, ext, time_str) assert len(files) == len( - folder_dates + folder_times ), "#files does not match #folders, check time_str and try again" success = 0 - for file, date in zip(files, folder_dates): + for file, date in zip(files, folder_times): new_name = file.with_name(f"{date}_{file.name}") if copy: shutil.copy(str(file), new_name) From 1874897dab4e506882f536a4f379502d08b14219 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 20 Mar 2024 13:09:11 -0400 Subject: [PATCH 41/47] small bugfix --- neuropy/utils/signal_process.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neuropy/utils/signal_process.py b/neuropy/utils/signal_process.py index dd258a22..4204c956 100644 --- a/neuropy/utils/signal_process.py +++ b/neuropy/utils/signal_process.py @@ -186,7 +186,7 @@ def get_pe_mean_spec( event_times, buffer_sec=(0.5, 0.5), ignore_epochs: core.Epoch = None, - print_ignored_frames: bool = True, + print_ignored_frames: bool = False, ): """Get peri-event mean spectrogram. From a7e667756039258e8fb44b9611555f2952c25df6 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 20 Mar 2024 15:00:15 -0400 Subject: [PATCH 42/47] small bugfix --- neuropy/core/epoch.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/neuropy/core/epoch.py b/neuropy/core/epoch.py index 7f578c15..df35b16e 100644 --- a/neuropy/core/epoch.py +++ b/neuropy/core/epoch.py @@ -723,6 +723,8 @@ def from_boolean_array(arr, t=None): core.Epoch epochs where the arr is high """ + if isinstance(t, pd.Series): # grab values only + t = t.values assert np.array_equal(arr, arr.astype(bool)), "Only boolean array accepted" int_arr = arr.astype("int") pad_arr = np.pad(int_arr, 1) @@ -732,6 +734,7 @@ def from_boolean_array(arr, t=None): if t is not None: assert len(t) == len(arr), "time length should be same as input array" + starts, stops = t[starts], t[stops] return Epoch.from_array(starts, stops, "high") From e4eefc94b964d3f56b52627d1673d2de56b0c74d Mon Sep 17 00:00:00 2001 From: nkinsky Date: Wed, 20 Mar 2024 15:00:55 -0400 Subject: [PATCH 43/47] enh - allow filtering by include_str to select only specific folders when loading miniscope timestamps --- neuropy/io/dlcio.py | 20 +++++++++++++++++--- neuropy/io/miniscopeio.py | 20 +++++++++++++++++--- 2 files changed, 34 insertions(+), 6 deletions(-) diff --git a/neuropy/io/dlcio.py b/neuropy/io/dlcio.py index 85556b5b..362472b0 100644 --- a/neuropy/io/dlcio.py +++ b/neuropy/io/dlcio.py @@ -80,6 +80,7 @@ def __init__(self, base_path, search_str=None, movie_type=".mp4", pix2cm=1): # Initialize other fields self.timestamp_type = None + self.speed = None @property def SampleRate(self): @@ -119,10 +120,12 @@ def get_metadata(self): with open(meta_file, "rb") as f: self.meta = load(f) - def get_timestamps(self, camera_type: str in ['optitrack', 'ms_webcam', 'ms_webcam1', 'ms_webcam2'] = 'optitrack'): + def get_timestamps(self, camera_type: str in ['optitrack', 'ms_webcam', 'ms_webcam1', 'ms_webcam2'] = 'optitrack', + exclude_str: str or None = None, include_str: str or None = None): """Tries to import timestamps from CSV file from optitrack, if not, infers it from timestamp in filename, sample rate, and nframes :param camera_type: 'optitrack' looks for optitrack csv file with tracking data, other options look for + :params exclude_str, include_str: see MiniscopeIO.load_all_timestamps - can include or exclude folders UCLA miniscope webcam files""" assert camera_type in ['optitrack', 'ms_webcam', 'ms_webcam1', 'ms_webcam2'] @@ -148,7 +151,8 @@ def get_timestamps(self, camera_type: str in ['optitrack', 'ms_webcam', 'ms_webc else: mio = MiniscopeIO(self.base_dir) webcam_flag = True if camera_type == "ms_webcam" else int(camera_type[-1]) - self.timestamps = mio.load_all_timestamps(webcam=webcam_flag) + self.timestamps = mio.load_all_timestamps(webcam=webcam_flag, exclude_str=exclude_str, + include_str=include_str) def to_cm(self): @@ -191,6 +195,16 @@ def smooth_pos( SampleRate=self.SampleRate, ) + def get_all_speed(self): + """Get speed of all bodyparts""" + df_list = [] + for bodypart in self.bodyparts: + df_list.append(pd.DataFrame({bodypart: self.get_speed(bodypart)})) + + self.speed = pd.concat(df_list, axis=1) + + return self.speed + def get_speed(self, bodypart): """Get speed of animal""" assert self.pos_smooth is not None, "You must smooth data first" @@ -210,7 +224,7 @@ def get_speed(self, bodypart): # Calculate speed speed = delta_pos / delta_t - + speed = np.hstack((speed, np.nan)) # add in Nan at end to match pos data shape return speed diff --git a/neuropy/io/miniscopeio.py b/neuropy/io/miniscopeio.py index 3b81ad28..3103d6ff 100644 --- a/neuropy/io/miniscopeio.py +++ b/neuropy/io/miniscopeio.py @@ -20,7 +20,8 @@ def __init__(self, basedir) -> None: pass def load_all_timestamps( - self, format="UCLA", webcam: bool or int = False, exclude_str: str = "WebCam", print_included_folders=False + self, format="UCLA", webcam: bool or int = False, exclude_str: str = "WebCam", + include_str = None, print_included_folders=False ): """Loads multiple timestamps from multiple videos in the UCLA miniscope software file format (folder = ""hh_mm_ss") @@ -59,6 +60,17 @@ def load_all_timestamps( rec_folder2.append(folder) self.rec_folders = rec_folder2 + # Include only folders with "include_str" in their name + if include_str is not None: + assert exclude_str is None, "cannot combine 'include_str' and 'exclude_str'" + rec_folder2 = [] + for folder in self.rec_folders: + if re.search(include_str, str(folder)) is not None: + if print_included_folders: + print("including folder " + str(folder)) + rec_folder2.append(folder) + self.rec_folders = rec_folder2 + # Loop through and load all timestamps, then concatenate times_list = [] for rec_folder in self.rec_folders: @@ -356,5 +368,7 @@ def euler_from_quaternion(qx, qy, qz, qw): if __name__ == "__main__": - parent_folder = '/data2/Psilocybin/Recording_Rats/Finn2' - move_files_to_combined_folder(parent_folder, prepend_rec_date=True, copy_during_prepend=True) + parent_folder = '/data2/Trace_FC/Recording_Rats/Django/2023_03_08_training' + mio = MiniscopeIO(parent_folder) + mio.load_all_timestamps(webcam=True, print_included_folders=True) + From 6033b4719f29e40d8441717276684278d00deae1 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Thu, 28 Mar 2024 13:38:30 -0400 Subject: [PATCH 44/47] lots of bugfixes for importing multiple DLC files for same session --- neuropy/io/dlcio.py | 169 +++++++++++++++++++++++++++----------------- 1 file changed, 105 insertions(+), 64 deletions(-) diff --git a/neuropy/io/dlcio.py b/neuropy/io/dlcio.py index 362472b0..2f0e8e17 100644 --- a/neuropy/io/dlcio.py +++ b/neuropy/io/dlcio.py @@ -1,6 +1,6 @@ """Class to deal with DeepLabCut data""" -from .movie import tracking_movie, deduce_starttime_from_file +from neuropy.io.movie import tracking_movie, deduce_starttime_from_file import numpy as np import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec @@ -14,7 +14,7 @@ from neuropy.io.miniscopeio import MiniscopeIO -from .optitrackio import getStartTime, posfromCSV +from neuropy.io.optitrackio import getStartTime, posfromCSV class DLC: @@ -52,73 +52,105 @@ def __init__(self, base_path, search_str=None, movie_type=".mp4", pix2cm=1): # Pull out only the specific file you want if specified with a search_str above, otherwise grab first one only! # NRK todo: use walrus operator here to check and make sure output file list len == 1? if search_str is not None: - self.tracking_file = Path( - get_matching_files(tracking_files, match_str=search_str)[0] + self.tracking_files = Path( + get_matching_files(tracking_files, match_str=search_str) ) else: - self.tracking_file = Path(tracking_files[0]) - self.movie_file = self.tracking_file.with_suffix(movie_type) - print("Using tracking file " + str(self.tracking_file)) + self.tracking_files = [Path(file) for file in tracking_files] + self.movie_files = [tracking_file.with_suffix(movie_type) for tracking_file in self.tracking_files] + + self.movies = [] + # self.nframes = self.nframes + # self.SampleRate = [] + self.timestamps = [] + + pos_list, scorernames = [], [] + for idf, file in enumerate(self.tracking_files): + print(f"Using tracking file #{idf+1}: {str(file)}") + # Import position data as-is + pos_data = import_tracking_data(file) + scorername = get_scorername(pos_data) + scorernames.append(scorername) + # Now grab the data for the scorer identified above - easier access later! + pos_data = pos_data[scorername] + pos_list.append(pos_data) + self.pos_data = pd.concat(pos_list, axis=0) self.pix2cm = pix2cm - # First import position data as-is - pos_data = import_tracking_data(self.tracking_file) - # Assume you've only run one scorer and that bodyparts are the same for all files - self.scorername = get_scorername(pos_data) + self.scorernames = scorername self.bodyparts = get_bodyparts(pos_data) - # Now grab the data for the scorer identified above - easier access later! - self.pos_data = pos_data[self.scorername] + # # Now grab the data for the scorer identified above - easier access later! + # self.pos_data = pos_data[self.scorername] # Now convert from pixels to centimeters self.to_cm() # Grab metadata for later access + self.meta = [] self.get_metadata() + self.nframes = self.get_nframes() + self.SampleRate = self.get_SampleRate() # Initialize other fields self.timestamp_type = None self.speed = None - @property - def SampleRate(self): - try: - SampleRate = self.meta["data"]["fps"] - except KeyError: # try to import from movie directly + def get_SampleRate(self): + """Get mean sample rate across all movies""" + SampleRate = [] + for idm, movie_file in enumerate(self.movie_files): + try: + SampleRate.append(self.meta[idm]["data"]["fps"]) + meta_found = True + except KeyError: # try to import from movie directly + meta_found = False # Initialize movie - if self.movie_file.is_file(): - self.movie = tracking_movie(self.movie_file) - SampleRate = self.movie.get_sample_rate() - else: - SampleRate = None + if movie_file.is_file(): + self.movies[idm] = tracking_movie(movie_file) + SampleRate.append(self.movies[idm].get_sample_rate()) + else: + SampleRate.append(None) + + SampleRate = np.mean(SampleRate) + if not meta_found: + print("No metadata found - taking mean sample rate from all raw movies in self.movie_files") + if idm > 0: + print("Multiple videos found - taking mean sample rate from all videos") return SampleRate - @property - def nframes(self): - try: - nframes = self.meta["data"]["nframes"] - except KeyError: # try to import from movie directly - # Initialize movie - if self.movie_file.is_file(): - self.movie = tracking_movie(self.movie_file) - nframes = self.movie.get_nframes() - else: - nframes = None + def get_nframes(self): + + nframes = [] + for idm, movie_file in enumerate(self.movie_files): + try: + nframes.append(self.meta[idm]["data"]["nframes"]) + meta_found = True + except KeyError: # try to import from movie directly + meta_found = False + # Initialize movie + if movie_file.is_file(): + self.movies[idm] = tracking_movie(movie_file) + nframes.append(self.movies[idm].get_nframes()) + else: + nframes.append(None) + + if not meta_found: + print("No metadata found - getting nframes directly from all raw movies in self.movie_files") return nframes def get_metadata(self): """Load in meta-data corresponding to tracking file""" + self.meta = [] + for tracking_file in self.tracking_files: + meta_file = tracking_file.parent / (tracking_file.stem + "_meta.pickle") - meta_file = self.tracking_file.parent / ( - self.tracking_file.stem + "_meta.pickle" - ) - - with open(meta_file, "rb") as f: - self.meta = load(f) + with open(meta_file, "rb") as f: + self.meta.append(load(f)) def get_timestamps(self, camera_type: str in ['optitrack', 'ms_webcam', 'ms_webcam1', 'ms_webcam2'] = 'optitrack', exclude_str: str or None = None, include_str: str or None = None): @@ -130,24 +162,26 @@ def get_timestamps(self, camera_type: str in ['optitrack', 'ms_webcam', 'ms_webc assert camera_type in ['optitrack', 'ms_webcam', 'ms_webcam1', 'ms_webcam2'] self.timestamp_type = camera_type + self.timestamps = [] if camera_type == 'optitrack': - opti_file = self.tracking_file.parent / ( - self.tracking_file.stem[: self.tracking_file.stem.find("-Camera")] + ".csv" - ) - if opti_file.is_file(): - start_time = getStartTime(opti_file) - _, _, _, t = posfromCSV(opti_file) - - self.timestamps = start_time + pd.to_timedelta(t, unit="sec") - else: - print( - "No Optitrack csv file found, inferring start time from file name. SECOND PRECISION IN START TIME!!!" - ) - start_time = deduce_starttime_from_file(self.tracking_file) - time_deltas = pd.to_timedelta( - np.arange(self.nframes) / self.SampleRate, unit="sec" + for idt, tracking_file in enumerate(self.tracking_files): + opti_file = tracking_file.parent / ( + tracking_file.stem[: tracking_file.stem.find("-Camera")] + ".csv" ) - self.timestamps = start_time + time_deltas + if opti_file.is_file(): + start_time = getStartTime(opti_file) + _, _, _, t = posfromCSV(opti_file) + + self.timestamps.append(start_time + pd.to_timedelta(t, unit="sec")) + else: + print(f"No Optitrack csv file found at {tracking_file.stem}.") + print("Inferring start time from file name. SECOND PRECISION IN START TIME!!!") + start_time = deduce_starttime_from_file(tracking_file) + time_deltas = pd.to_timedelta( + np.arange(self.nframes[idt]) / self.SampleRate, unit="sec" + ) + self.timestamps.append(start_time + time_deltas) + self.timestamps = pd.concat(self.timestamps, axis=0) else: mio = MiniscopeIO(self.base_dir) webcam_flag = True if camera_type == "ms_webcam" else int(camera_type[-1]) @@ -216,7 +250,7 @@ def get_speed(self, bodypart): data_use = self.pos_smooth[bodypart] x = data_use["x"] y = data_use["y"] - total_seconds = (self.timestamps['Timestamps'] - self.timestamps['Timestamps'][0]).dt.total_seconds() + total_seconds = (self.timestamps.reset_index()['Timestamps'] - self.timestamps.reset_index()['Timestamps'][0]).dt.total_seconds() # Get delta in position and time from frame-to-frame delta_pos = np.sqrt(np.square(np.diff(x)) + np.square(np.diff(y))) @@ -529,11 +563,18 @@ def plot_2d( if __name__ == "__main__": - dlc = DLC( - "/data2/Trace_FC/Pilot1/Rat700/2021_02_22_habituation", - search_str="Camera 1", - pix2cm=0.13, - ) - dlc.smooth_pos() - dlc.get_freezing_epochs() + # dlc_path = '/data2/Trace_FC/Recording_Rats/Finn2/2023_05_08_training' + # dlc_path = '/data2/Trace_FC/Recording_Rats/Han/2022_08_03_training' + dlc_path = '/data2/Trace_FC/Recording_Rats/Rose/2022_06_22_training' + arena_side_pix = 60 # Keep this + arena_side_cm = 25.4 # Update this after measuring!!! + pix2cm = arena_side_cm / arena_side_pix + + # Read in DLC data + dlc = DLC(dlc_path, pix2cm=pix2cm) + + # Smooth position, get timestamps, and get speed + dlc.get_timestamps('ms_webcam', include_str="/2_training/") + dlc.smooth_pos(bodyparts=["crown_middle", "back_middle"]) + dlc.get_all_speed() pass From 57269de11e34097b70385abf28491ccb5e27ee22 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Thu, 28 Mar 2024 13:38:37 -0400 Subject: [PATCH 45/47] lots of bugfixes for importing multiple DLC files for same session --- neuropy/io/miniscopeio.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/neuropy/io/miniscopeio.py b/neuropy/io/miniscopeio.py index 3103d6ff..6c173e02 100644 --- a/neuropy/io/miniscopeio.py +++ b/neuropy/io/miniscopeio.py @@ -173,7 +173,13 @@ def load_timestamps( rec_metadata, vid_metadata, vid_folder = get_recording_metadata(rec_folder) if webcam: - folder_name = Path("My_WebCam" if isinstance(webcam, bool) else f"My_WebCam{webcam}") + if isinstance(webcam, bool): + folder_name = "My_WebCam" # Path("My_WebCam") + elif webcam == 1: + folder_name = "My_First_Webcam" + elif webcam == 2: + folder_name = "My_Second_Webcam" + # folder_name = Path("My_WebCam" if isinstance(webcam, bool) else f"My_WebCam{webcam}") vid_folder = vid_folder.parent / folder_name # Derive start_time from rec_metadata From 2b089b35fcf50efe09231c2bfd3b57a9804ccfd5 Mon Sep 17 00:00:00 2001 From: nkinsky Date: Thu, 28 Mar 2024 13:39:01 -0400 Subject: [PATCH 46/47] add useful assertion --- neuropy/plotting/signals.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/neuropy/plotting/signals.py b/neuropy/plotting/signals.py index 04a73abe..6eb10112 100644 --- a/neuropy/plotting/signals.py +++ b/neuropy/plotting/signals.py @@ -43,6 +43,9 @@ def plot_spectrogram( legacy = False spec = sxx sxx = spec.traces + assert time_lims[0] >= spec.t_start, "time_lims[0] must be greater than or equal to t_start in input Spectrogram" + assert time_lims[1] <= spec.t_stop, "time_lims[1] must be less than or equal to t_start in input Spectrogram" + if sigma is not None: sxx = gaussian_filter(sxx, sigma=sigma) @@ -57,7 +60,6 @@ def plot_spectrogram( # ---------- plotting ---------------- if legacy: - print('test') def plotspec(n_std, freq_lim): """Plots data fine but doesn't preserve time and frequency info on axes""" ax.imshow( From b91dd1a14e5ca48533314f1ae871adbf7fbfeb6a Mon Sep 17 00:00:00 2001 From: Nat Date: Fri, 5 Apr 2024 16:58:01 -0400 Subject: [PATCH 47/47] commenting --- neuropy/analyses/oscillations.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/neuropy/analyses/oscillations.py b/neuropy/analyses/oscillations.py index 4b756b50..a221c7e4 100644 --- a/neuropy/analyses/oscillations.py +++ b/neuropy/analyses/oscillations.py @@ -28,7 +28,7 @@ def _detect_freq_band_epochs( low and high threshold for detection mindur : float, optional minimum duration of epoch - maxdur : float, optiona + maxdur : float, optional chans : list channels used for epoch detection, if None then chooses best chans """ @@ -41,6 +41,8 @@ def _detect_freq_band_epochs( # Because here one shank is selected per shank, based on visualization: # mean: very conservative in cases where some shanks may not have that strong ripple # max: works well but may have ocassional false positives + + # First, bandpass the signal in the range of interest power = np.zeros(signals.shape[1]) for sig in signals: yf = signal_process.filter_sig.bandpass(sig, lf=lf, hf=hf, fs=fs) @@ -48,9 +50,10 @@ def _detect_freq_band_epochs( # zscsignal[sig_i] = zsc_chan power += np.abs(signal_process.hilbertfast(yf)) - # mean and smooth + # Second, take the mean and smooth the signal with a sigma wide gaussian kernel power = smooth(power / signals.shape[0]) + # Third, exclude any noisy periods due to motion or other artifact # ---------setting noisy periods zero -------- if ignore_times is not None: assert ignore_times.ndim == 2, "ignore_times should be 2 dimensional array" @@ -65,9 +68,12 @@ def _detect_freq_band_epochs( noisy_frames = noisy_frames[noisy_frames < len(power)] power[noisy_frames] = 0 + # Fourth, identify candidate epochs above edge_cutoff threshold # ---- thresholding and detection ------ power = stats.zscore(power) power_thresh = np.where(power >= edge_cutoff, power, 0) + + # Fifth, refine candidate epochs to periods between lowthresh and highthresh peaks, props = sg.find_peaks( power_thresh, height=[lowthresh, highthresh], prominence=0 ) @@ -75,6 +81,7 @@ def _detect_freq_band_epochs( peaks_power = power_thresh[peaks] # ----- merge overlapping epochs ------ + # Last, merge any epochs that overlap into one longer epoch n_epochs = len(starts) ind_delete = [] for i in range(n_epochs - 1):