Merge pull request #143 from nkinsky/main

Merge Nat fork

neuropy/analyses/artifact.py

@@ -13,7 +13,6 @@ def detect_artifact_epochs(
     edge_cutoff=2,
     merge=5,
     filt: list or np.ndarray = None,
-    data_use: str in ["raw_only", "filt_only", "both"] = "raw_only",
 ):
     """
     calculating artifact periods using z-score measure
@@ -33,7 +32,7 @@ def detect_artifact_epochs(
     filt : list, optional
         lower and upper limits with which to filter signal, e.g. 3, 3000] ->
         bandpass between and 3000 Hz while [45, None] -> high-pass above 45.
-    data_use: str, optional
+    data_use: str, optional (Not currently implemented)
         'raw_only' (default): z-score raw only
         'filt_only': z-score filtered data only
         'both': z-score both raw and filtered data
@@ -119,4 +118,4 @@ def detect_artifact_epochs(
 
 
 if __name__ == "__main__":
-    print('test')
+    print("test")

neuropy/core/epoch.py

@@ -31,6 +31,7 @@ def __init__(
                 file is not None
             ), "Must specify file to load if no epochs dataframe entered"
             epochs = np.load(file, allow_pickle=True).item()["epochs"]
+            self.metadata = np.load(file, allow_pickle=True).item()["metadata"]
 
         self._epochs = self._validate(epochs)
 
@@ -91,7 +92,13 @@ def __add__(self, epochs):
         return Epoch(epochs=df_new)
 
     def add_epoch_manually(self, start, stop, label="", merge_dt: float or None = 0):
-        comb_df = pd.DataFrame({"start": [start], "stop": [stop], "label": label})
+        comb_df = pd.DataFrame(
+            {
+                "start": np.array(start).reshape(-1),
+                "stop": np.array(stop).reshape(-1),
+                "label": label,
+            }
+        )
 
         if merge_dt is not None:
             return self.__add__(Epoch(comb_df)).merge(merge_dt)
@@ -691,3 +698,10 @@ def get_indices_for_time(self, t: np.array):
             time_bool[np.where((t >= e[0]) & (t <= e[1]))[0]] = 1
 
         return time_bool.astype("bool")
+
+
+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)
+    epochs_to_add = np.array([[1291, 1291.2], [2734, 2734.5], [1622, 1623.4]])
+    art_epochs.add_epoch_manually(epochs_to_add[:, 0], epochs_to_add[:, 1])

neuropy/core/session.py

@@ -0,0 +1,35 @@
+from pathlib import Path
+
+from neuropy.io.neuroscopeio import NeuroscopeIO
+from neuropy.io.binarysignalio import BinarysignalIO
+
+
+class ProcessData:
+    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"
+
+        fp = xml_files[0].with_suffix("")
+        self.filePrefix = fp
+
+        self.recinfo = NeuroscopeIO(xml_files[0])
+        eegfiles = sorted(basepath.glob("*.eeg"))
+        assert len(eegfiles) == 1, "Fewer/more than one .eeg file detected"
+        self.eegfile = BinarysignalIO(
+            eegfiles[0],
+            n_channels=self.recinfo.n_channels,
+            sampling_rate=self.recinfo.eeg_sampling_rate,
+        )
+        try:
+            self.datfile = BinarysignalIO(
+                eegfiles[0].with_suffix(".dat"),
+                n_channels=self.recinfo.n_channels,
+                sampling_rate=self.recinfo.dat_sampling_rate,
+            )
+        except FileNotFoundError:
+            print("No dat file found, not loading")
+
+    def __repr__(self) -> str:
+        return f"{self.__class__.__name__}({self.recinfo.source_file.name})"

neuropy/core/signal.py

@@ -76,7 +76,13 @@ def time_slice(self, channel_id=None, t_start=None, t_stop=None):
             channel_indx = [list(self.channel_id).index(_) for _ in channel_id]
             traces = self.traces[channel_indx, frame_start:frame_stop]
 
-        return Signal(traces, self.sampling_rate, t_start, channel_id)
+        return Signal(
+            traces,
+            self.sampling_rate,
+            t_start,
+            channel_id,
+            source_file=self.source_file,
+        )
 
     def rescale(self, factor=0.95 * 1e-3):
         """scales signal, use it for converting raw signal to volts, but can consume too much memory and time when used on large memmap arrays

neuropy/core/sync.py

@@ -5,37 +5,6 @@
 from neuropy.io.neuroscopeio import NeuroscopeIO
 
 
-class ProcessData:
-    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"
-
-        fp = xml_files[0].with_suffix("")
-        self.filePrefix = fp
-
-        self.recinfo = NeuroscopeIO(xml_files[0])
-        eegfiles = sorted(basepath.glob("*.eeg"))
-        assert len(eegfiles) == 1, "Fewer/more than one .eeg file detected"
-        self.eegfile = BinarysignalIO(
-            eegfiles[0],
-            n_channels=self.recinfo.n_channels,
-            sampling_rate=self.recinfo.eeg_sampling_rate,
-        )
-        try:
-            self.datfile = BinarysignalIO(
-                eegfiles[0].with_suffix(".dat"),
-                n_channels=self.recinfo.n_channels,
-                sampling_rate=self.recinfo.dat_sampling_rate,
-            )
-        except FileNotFoundError:
-            print("No dat file found, not loading")
-
-    def __repr__(self) -> str:
-        return f"{self.__class__.__name__}({self.recinfo.source_file.name})"
-
-
 class OESync(DataWriter):
     """Class to synchronize different systems with ephys data collected with OpenEphys. Most useful when multiple dat
     files have been concatenated into one, but also useful for individual recordings."""

neuropy/io/neuroscopeio.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 from pathlib import Path
 import xml.etree.ElementTree as Etree
 from .. import core
@@ -114,7 +115,14 @@ def write_neurons(self, neurons: core.Neurons, suffix_num: int = 1):
     def write_epochs(self, epochs: core.Epoch, ext="epc"):
         with self.source_file.with_suffix(f".evt.{ext}").open("w") as a:
             for event in epochs.to_dataframe().itertuples():
-                a.write(f"{event.start*1000} start\n{event.stop*1000} stop\n")
+                # First attempt to fix bug where Neuropy exported .evt files get broken after manual
+                # adjustment in NeuroScope - does not seem to work
+                event_start, event_stop = event.start * 1000, event.stop * 1000
+                if np.mod(event_start, 1) == 0:
+                    event_start += 0.2
+                if np.mod(event_stop, 1) == 0:
+                    event_stop += 0.2
+                a.write(f"{event_start}\tstart\n{event_stop}\tstop\n")
 
     def write_position(self, position: core.Position):
         """Writes core.Position object to neuroscope compatible format
@@ -143,3 +151,25 @@ def to_dict(self):
             "dat_sampling_rate": self.dat_sampling_rate,
             "eeg_sampling_rate": self.eeg_sampling_rate,
         }
+
+    def event_to_epochs(self, evt_file, label=""):
+        """Read in an event file and convert to an epochs object"""
+        with open(evt_file, "r") as f:
+            Lines = f.readlines()
+
+        # Neuropy output saves file without tab separators
+        if Lines[0].find("\t") > -1:
+            split_str = "\t"
+        else:  # if you savne in Neuroscope the event file now has tab separators
+            split_str = " "
+
+        starts, stops = [], []
+        for line in Lines:
+            if line.find("start") > -1:
+                starts.append(float(line.split(f"{split_str}start")[0]) / 1000)
+            elif line.find("stop") > -1:
+                stops.append(float(line.split(f"{split_str}stop")[0]) / 1000)
+
+        return core.Epoch(
+            pd.DataFrame({"start": starts, "stop": stops, "label": label})
+        )

neuropy/io/openephysio.py

@@ -28,7 +28,12 @@ def get_us_start(settings_file: str, from_zone="UTC", to_zone="America/Detroit")
     return dt_start_utc.astimezone(to_zone)
 
 
-def get_dat_timestamps(basepath: str or Path, sync: bool = False, start_end_only=False):
+def get_dat_timestamps(
+    basepath: str or Path,
+    sync: bool = False,
+    start_end_only=False,
+    local_time="America/Detroit",
+):
     """
     Gets timestamps for each frame in your dat file(s) in a given directory.
 
@@ -56,8 +61,8 @@ def get_dat_timestamps(basepath: str or Path, sync: bool = False, start_end_only
         except KeyError:
             try:
                 experiment_meta = XML2Dict(set_folder / set_file)  # Get meta data
-                start_time = pd.Timestamp(
-                    experiment_meta["INFO"]["DATE"]
+                start_time = pd.Timestamp(experiment_meta["INFO"]["DATE"]).tz_localize(
+                    local_time
                 )  # get start time from meta-data
             except FileNotFoundError:
                 print(
@@ -70,7 +75,9 @@ def get_dat_timestamps(basepath: str or Path, sync: bool = False, start_end_only
                     "[0-9]{4}-[0-9]{2}-[0-9]{2}_[0-9]{2}-[0-9]{2}-[0-9]{2}",
                     str(set_folder),
                 )
-                start_time = pd.to_datetime(m.group(0), format="%Y-%m-%d_%H-%M-%S")
+                start_time = pd.to_datetime(
+                    m.group(0), format="%Y-%m-%d_%H-%M-%S"
+                ).tz_localize(local_time)
 
         SR, sync_frame = parse_sync_file(
             file.parents[3] / "recording1/sync_messages.txt"
@@ -367,10 +374,21 @@ def recording_events_to_combined_time(
         event_time_comb = np.array(event_time_comb)
 
     else:
+        good_bool = [start == stop for start, stop in zip(nrec_start, nrec_stop)]
+        good_events = np.where(good_bool)[0]
+        bad_events = np.where(~np.array(good_bool))[0]
+
         print(
-            f"Recording start and end numbers do not all match. starts = {nrec_start}, ends = {nrec_stop}."
+            f"Event(s) # {bad_events + 1} occurs in between recordings and has(have) been left out"
+        )
+        # print(
+        #     f"Recording start and end numbers do not all match. starts = {nrec_start}, ends = {nrec_stop}."
+        # )
+        # event_time_comb = np.nan
+
+        event_time_comb = recording_events_to_combined_time(
+            event_df.iloc[good_events], sync_df, time_out, event_ts_key, sync_ts_key
         )
-        event_time_comb = np.nan
 
     return event_time_comb
 
@@ -716,15 +734,17 @@ def GetRecChs(File):
 
 if __name__ == "__main__":
     import tracefc.io.traceio as traceio
-    basepath= Path("/Users/nkinsky/Documents/UM/Working/Trace_FC/Recording_Rats/Finn/2022_01_21_recall1/")
+
+    basepath = Path("/data3/Trace_FC/Recording_Rats/Finn2/2023_05_06_habituation1")
     ttl_df = load_all_ttl_events(basepath, sanity_check_channel=1, zero_timestamps=True)
-    cs_starts, cs_ends, cs_df = traceio.load_trace_events(basepath, session_type="tone_recall",
-                                                          event_type="CS+", return_df=True)
+    cs_starts, cs_ends, cs_df = traceio.load_trace_events(
+        basepath, session_type="tone_recall", event_type="CS+", return_df=True
+    )
     sync_df = create_sync_df(basepath)
     ttl_lag_use = ttl_lag = pd.Timedelta(0.8, unit="seconds")
-    cs_oe_start_df = traceio.trace_ttl_to_openephys(cs_starts,
-                                                    ttl_df[ttl_df['channel_states'].abs() == 2],
-                                                    ttl_lag=ttl_lag_use)
+    cs_oe_start_df = traceio.trace_ttl_to_openephys(
+        cs_starts, ttl_df[ttl_df["channel_states"].abs() == 2], ttl_lag=ttl_lag_use
+    )
     # Convert to times in combined eeg file
     cs_starts_combined = recording_events_to_combined_time(cs_oe_start_df, sync_df)
 

neuropy/utils/signal_process.py

@@ -104,7 +104,11 @@ def get_noisy_spect_bool(self, thresh=5):
         return (stats.zscore(spect_sum) >= thresh) | (spect_sum <= 0)
 
     def get_pe_mean_spec(
-        self, event_times, buffer_sec=(0.5, 0.5), ignore_epochs: core.Epoch = None
+        self,
+        event_times,
+        buffer_sec=(0.5, 0.5),
+        ignore_epochs: core.Epoch = None,
+        print_ignored_frames: bool = True,
     ):
         """Get peri-event mean spectrogram
 
@@ -161,9 +165,10 @@ def get_pe_mean_spec(
                 # Display ignored frames
                 if np.sum(ignore_bool) > 0:
                     sxx_temp[:, ignore_bool] = np.nan
-                    print(
-                        f"{np.sum(ignore_bool)} frames between {ignore_times.min():.1F} and {ignore_times.max():.1F} ignored (sent to 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)"
+                        )
             sxx_list.append(sxx_temp)
 
         sxx_mean = np.nanmean(np.stack(sxx_list, axis=2), axis=2)