Address nan issue by trimming sweeps

docs/gallery/analysis_examples/lsq_analysis.py

@@ -36,9 +36,7 @@
 )
 long_square_sweeps = data_set.sweep_set(long_square_table.sweep_number)
 
-# Select epoch corresponding to the actual recording from the sweeps
-# and align sweeps so that the experiment would start at the same time
-long_square_sweeps.select_epoch("recording")
+# Align sweeps so that the experiment would start at the same time
 long_square_sweeps.align_to_start_of_epoch("experiment")
 
 # find the start and end time of the stimulus

docs/gallery/analysis_examples/ramp_analysis.py

@@ -37,14 +37,7 @@
 )
 ramp_sweeps = data_set.sweep_set(ramp_table.sweep_number)
 
-# Select epoch corresponding to the actual recording from the sweeps
-# and align sweeps so that the experiment would start at the same time
-ramp_sweeps.select_epoch("recording")
-ramp_sweeps.align_to_start_of_epoch("experiment")
-
-# Select epoch corresponding to the actual recording from the sweeps
-# and align sweeps so that the experiment would start at the same time
-ramp_sweeps.select_epoch("recording")
+# Align sweeps so that the experiment would start at the same time
 ramp_sweeps.align_to_start_of_epoch("experiment")
 
 # find the start and end time of the stimulus

docs/gallery/analysis_examples/short_square_analysis.py

@@ -35,9 +35,7 @@
 )
 short_square_sweeps = data_set.sweep_set(short_square_table.sweep_number)
 
-# Select epoch corresponding to the actual recording from the sweeps
-# and align sweeps so that the experiment would start at the same time
-short_square_sweeps.select_epoch("recording")
+# Align sweeps so that the experiment would start at the same time
 short_square_sweeps.align_to_start_of_epoch("experiment")
 
 # find the start and end time of the stimulus

ipfx/bin/plot_ephys_nwb.py

@@ -1,34 +1,68 @@
 import sys
+import warnings
 import matplotlib.pyplot as plt
 import numpy as np
 from ipfx.dataset.create import create_ephys_data_set
+from ipfx.qc_feature_extractor import sweep_qc_features
+from ipfx.utilities import drop_failed_sweeps
 from ipfx.stimulus import StimulusOntology
 import allensdk.core.json_utilities as ju
-
-
-def plot_data_set(data_set, sweep_table, nwb_file_name):
-
+from typing import (
+    Optional, List, Dict, Tuple, Collection, Sequence, Union
+)
+
+
+def plot_data_set(data_set,
+            clamp_mode: Optional[str] = None,
+            stimuli: Optional[Collection[str]] = None,
+            stimuli_exclude: Optional[Collection[str]] = None,
+            show_amps: Optional[bool] = True,
+            qc_sweeps: Optional[bool] = True,
+            figsize=(15, 7),
+    ):
+    nwb_file_name = str(data_set._data.nwb_file)
+    if qc_sweeps:
+        drop_failed_sweeps(data_set)
+    elif show_amps:
+        data_set.sweep_info = sweep_qc_features(data_set)
+
+    sweep_table = data_set.filtered_sweep_table(clamp_mode=clamp_mode, stimuli=stimuli, stimuli_exclude=stimuli_exclude)
+
+    if len(sweep_table)==0:
+        warnings.warn("No sweeps to plot")
+        return
+
     height_ratios, width_ratios = axes_ratios(sweep_table)
 
     fig, ax = plt.subplots(len(height_ratios), 3,
-                           figsize=(15, 7),
+                           figsize=figsize,
                            gridspec_kw={'height_ratios': height_ratios, 'width_ratios': width_ratios}
                            )
+    if len(height_ratios)==1:
+        # ensure 2d array
+        ax = ax[None, :]
 
     for fig_row, (stimulus_code, sweep_set_table) in enumerate(sweep_table.groupby("stimulus_code")):
-        sweep_numbers = sweep_set_table["sweep_number"].sort_values().values
+        sweep_set_table = sweep_set_table.copy().sort_values("sweep_number", ascending=False)
+        sweep_numbers = sweep_set_table["sweep_number"]
         ss = data_set.sweep_set(sweep_numbers)
+        if qc_sweeps:
+            ss.select_epoch('experiment')
+        annot = sweep_numbers.astype(str)
+        if show_amps:
+            annot += sweep_set_table['stimulus_amplitude'].apply(": {:.3g} pA".format)
+
 
         ax_a = ax[fig_row,0]
         ax_i = ax[fig_row,1]
         ax_v = ax[fig_row,2]
 
-        plot_waveforms(ax_i, ss.i, ss.sampling_rate, ss.sweep_number)
-        plot_waveforms(ax_v, ss.v, ss.sampling_rate, ss.sweep_number)
+        plot_waveforms(ax_i, ss.i, ss.sampling_rate, annot)
+        plot_waveforms(ax_v, ss.v, ss.sampling_rate)
         ax_v.get_shared_x_axes().join(ax_i, ax_v)
 
         clamp_mode = sweep_set_table["clamp_mode"].values[0]
-        ax_a.text(0, 0.0, "%s, %s " % (stimulus_code, clamp_mode))
+        ax_a.text(0, 0.0, "%s \n%s " % (stimulus_code, clamp_mode))
         ax_a.axis('off')
 
     ax[0,0].set_title("Description")
@@ -40,41 +74,48 @@ def plot_data_set(data_set, sweep_table, nwb_file_name):
     fig.suptitle("file: " + nwb_file_name, fontsize=12)
 
     mng = plt.get_current_fig_manager()
-    mng.resize(*mng.window.maxsize())
+    if hasattr(mng, 'window'):
+        mng.resize(*mng.window.maxsize())
     plt.subplots_adjust(left=0.01, right=0.98, bottom=0.02,top=0.92)
 
 
 def axes_ratios(sweep_table):
 
     height_ratios = []
-    width_ratios = [1,3,3]
+    width_ratios = [1,4,4]
 
     for _, sweep_set_table in sweep_table.groupby("stimulus_code"):
         height_ratios.append(len(sweep_set_table.index))
 
     return height_ratios, width_ratios
 
 
-def plot_waveforms(ax, ys, rs, sn):
+def plot_waveforms(ax, ys, rs, annotations=None):
 
     offset = 0
-
-    for y, r, sn in zip(ys, rs, sn):
-
-        offset += get_vertical_offset(y)
+    dy = 0
+    for i, (y, r) in enumerate(zip(ys, rs)):
+        if len(y)==0:
+            continue
+        y -= y[0]
+        dy = max(dy, get_vertical_offset(y))
+        offset += dy
         y += offset
         x = np.arange(0, len(y)) / r
 
         ax.plot(x, y)
-        ax.text(x[0] - 0.05 * (x[-1] - x[0]), y[0], str(sn), fontsize=8)
-        ax.set_xlim(x[0], x[-1])
+        if annotations is not None:
+            ax.text(x[0] - 0.01 * (x[-1] - x[0]), y[0], annotations.iloc[i], fontsize=8, ha='right')
+    # need to set limits to show all sweeps
+    # mode would be best but mean will do fine
+    ax.set_xlim(0, np.max([len(y) for y in ys])/r)
 
     customize_axis(ax)
 
 
 def customize_axis(ax):
 
-    ax.tick_params(axis="x", direction="in", pad=-10,labelsize=8)
+    ax.tick_params(axis="x", direction="in", pad=3, labelsize=8)
     ax.get_yaxis().set_visible(False)
     ax.spines['top'].set_visible(False)
     ax.spines['right'].set_visible(False)
@@ -102,12 +143,8 @@ def main():
     ont = StimulusOntology(ju.read(stimulus_ontology_file))
 
     data_set = create_ephys_data_set(nwb_file=nwb_file)
-    vclamp_sweep_table = data_set.filtered_sweep_table(clamp_mode=data_set.VOLTAGE_CLAMP)
-    plot_data_set(data_set, vclamp_sweep_table, nwb_file)
-
-    data_set = create_ephys_data_set(nwb_file=nwb_file)
-    iclamp_sweep_table = data_set.filtered_sweep_table(clamp_mode=data_set.CURRENT_CLAMP)
-    plot_data_set(data_set, iclamp_sweep_table, nwb_file)
+    plot_data_set(data_set, clamp_mode=data_set.VOLTAGE_CLAMP)
+    plot_data_set(data_set, clamp_mode=data_set.CURRENT_CLAMP)
 
     plt.show()
 

ipfx/data_set_features.py

@@ -160,7 +160,6 @@ def extract_sweep_features(data_set, sweep_table):
         sweep_numbers = sorted(sweep_numbers)
 
         sweep_set = data_set.sweep_set(sweep_numbers)
-        sweep_set.select_epoch("recording")
         sweep_set.align_to_start_of_epoch("experiment")
 
         dp = detection_parameters(stimulus_name).copy()
@@ -238,7 +237,6 @@ def extract_cell_long_square_features(data_set, subthresh_min_amp=None):
             logging.info("Assigned subthreshold minimum amplitude of %f.", subthresh_min_amp)
 
     lsq_sweeps = data_set.sweep_set(long_square_sweep_numbers)
-    lsq_sweeps.select_epoch("recording")
     lsq_sweeps.align_to_start_of_epoch("experiment")
 
     lsq_start, lsq_dur, _, _, _ = stf.get_stim_characteristics(
@@ -278,7 +276,6 @@ def extract_cell_short_square_features(data_set):
         raise er.FeatureError("No short square sweeps available for feature extraction")
 
     ssq_sweeps = data_set.sweep_set(short_square_sweep_numbers)
-    ssq_sweeps.select_epoch("recording")
     ssq_sweeps.align_to_start_of_epoch("experiment")
 
     ssq_start, ssq_dur, _, _, _ = stf.get_stim_characteristics(
@@ -312,7 +309,6 @@ def extract_cell_ramp_features(data_set):
         raise er.FeatureError("No ramp sweeps available for feature extraction")
 
     ramp_sweeps = data_set.sweep_set(ramp_sweep_numbers)
-    ramp_sweeps.select_epoch("recording")
     ramp_sweeps.align_to_start_of_epoch("experiment")
 
     ramp_start, ramp_dur, _, _, _ = stf.get_stim_characteristics(ramp_sweeps.sweeps[0].i, ramp_sweeps.sweeps[0].t)

ipfx/dataset/ephys_data_set.py

@@ -316,11 +316,12 @@ def get_sweep_data(self, sweep_number: int) -> Dict:
 
         response = sweep_data['response']
 
-        nonzero = np.flatnonzero(response)
-        if len(nonzero) == 0:
+        valid_data = (response != 0) & ~np.isnan(response)
+        valid_i = np.flatnonzero(valid_data)
+        if len(valid_i) == 0:
             recording_end_idx = 0
         else:
-            recording_end_idx = nonzero[-1] + 1
+            recording_end_idx = valid_i[-1] + 1
 
         sweep_data["response"] = response[:recording_end_idx]
         sweep_data["stimulus"] = sweep_data["stimulus"][:recording_end_idx]

ipfx/epochs.py

@@ -51,27 +51,6 @@ def get_last_noise_epoch(idx1, hz):
     return idx1-int(NOISE_EPOCH * hz), idx1
 
 
-def get_recording_epoch(response):
-    """
-    Detect response epoch defined as interval from start to the last non-nan value of the response
-
-    Parameters
-    ----------
-    response: float np.array
-
-    Returns
-    -------
-    start,end: int
-        indices of the epoch
-    """
-
-    if len(tsu.flatnotnan(response)) == 0:
-        end_idx = 0
-    else:
-        end_idx = tsu.flatnotnan(response)[-1]
-    return 0, end_idx
-
-
 def get_sweep_epoch(response):
     """
     Defined as interval including entire sweep

ipfx/qc_feature_extractor.py

@@ -287,7 +287,7 @@ def sweep_qc_features(data_set):
             qc_features = current_clamp_sweep_qc_features(sweep, is_ramp)
             sweep_features.update(qc_features)
         else:
-            logging.warning("sweep {}: {}".format(sweep_num, tags))
+            logging.info("sweep {}: {}".format(sweep_num, tags))
 
         sweeps_features.append(sweep_features)
 
@@ -303,8 +303,8 @@ def check_sweep_integrity(sweep, is_ramp):
             tags.append(F"{k} epoch is missing")
 
     if not is_ramp:
-        if sweep.epochs["recording"] and sweep.epochs["experiment"]:
-            if sweep.epochs["recording"][1] < sweep.epochs["experiment"][1]:
+        if sweep.epochs["sweep"] and sweep.epochs["experiment"]:
+            if sweep.epochs["sweep"][1] < sweep.epochs["experiment"][1]:
                 tags.append("Recording stopped before completing the experiment epoch")
 
     if np.isnan(sweep.i).any():
@@ -353,7 +353,7 @@ def current_clamp_sweep_qc_features(sweep, is_ramp):
     # measure mean and rms of Vm at end of recording
     # do not check for ramps, because they do not have enough time to recover
 
-    _, rec_end_idx = ep.get_recording_epoch(voltage)
+    _, rec_end_idx = ep.get_sweep_epoch(voltage)
 
     if not is_ramp:
         idx0, idx1 = ep.get_last_stability_epoch(rec_end_idx, hz)

ipfx/script_utils.py

@@ -179,7 +179,6 @@ def categorize_iclamp_sweeps(data_set, stimuli_names, sweep_qc_option="none", sp
 
 def validate_sweeps(data_set, sweep_numbers, extra_dur=0.2):
     check_sweeps = data_set.sweep_set(sweep_numbers)
-    check_sweeps.select_epoch("recording")
     valid_sweep_stim = []
     start = None
     dur = None
@@ -216,7 +215,6 @@ def preprocess_long_square_sweeps(data_set, sweep_numbers, extra_dur=0.2, subthr
     lsq_sweeps, lsq_start, lsq_end = validate_sweeps(data_set, sweep_numbers, extra_dur=extra_dur)
     if len(lsq_sweeps.sweeps) == 0:
         raise er.FeatureError("No long square sweeps were long enough or did not end early")
-    lsq_sweeps.select_epoch("recording")
 
     lsq_spx, lsq_spfx = dsf.extractors_for_sweeps(
         lsq_sweeps,
@@ -239,7 +237,6 @@ def preprocess_short_square_sweeps(data_set, sweep_numbers, extra_dur=0.2, spike
     ssq_sweeps, ssq_start, ssq_end  = validate_sweeps(data_set, sweep_numbers, extra_dur=extra_dur)
     if len(ssq_sweeps.sweeps) == 0:
         raise er.FeatureError("No short square sweeps were long enough or did not end early")
-    ssq_sweeps.select_epoch("recording")
 
     ssq_spx, ssq_spfx = dsf.extractors_for_sweeps(ssq_sweeps,
                                                   est_window = [ssq_start, ssq_start + 0.001],
@@ -258,7 +255,6 @@ def preprocess_ramp_sweeps(data_set, sweep_numbers):
         raise er.FeatureError("No ramp sweeps available for feature extraction")
 
     ramp_sweeps = data_set.sweep_set(sweep_numbers)
-    ramp_sweeps.select_epoch("recording")
 
     ramp_start, ramp_dur, _, _, _ = stf.get_stim_characteristics(ramp_sweeps.sweeps[0].i, ramp_sweeps.sweeps[0].t)
     ramp_spx, ramp_spfx = dsf.extractors_for_sweeps(ramp_sweeps,

ipfx/stim_features.py

@@ -32,6 +32,7 @@ def get_stim_characteristics(i, t, test_pulse=True):
     else:
         amplitude = float(peak_low)
 
+    # this is very unlikely, but possible for zero stim with mid-sweep NaNs
     if np.isnan(amplitude):
         return None, None, 0.0, None, None
 

ipfx/sweep.py

@@ -62,7 +62,6 @@ def detect_epochs(self):
 
         epoch_detectors = {
             "sweep": ep.get_sweep_epoch(self.response),
-            "recording": ep.get_recording_epoch(self.response),
             "experiment": ep.get_experiment_epoch(self._i, self.sampling_rate,test_pulse),
             "stim": ep.get_stim_epoch(self.stimulus, test_pulse),
         }

tests/test_epochs.py

@@ -116,25 +116,6 @@ def test_get_stim_epoch(i, stim_epoch):
     assert stim_epoch == ep.get_stim_epoch(i)
 
 
-@pytest.mark.parametrize('response,'
-                         'recording_epoch',
-                         [
-                             (
-                                 [0, 0, 1, 1.5, 0, 0, 2, 3, 4, 1, np.nan, np.nan],
-                                 (0, 9)
-                             ),
-
-                             #   zero array
-                             (
-                                 [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-                                 (0, 0)
-                             ),
-
-                         ]
-                         )
-def test_get_recording_epoch(response, recording_epoch):
-    assert recording_epoch == ep.get_recording_epoch(response)
-
 
 @pytest.mark.parametrize('response,'
                          'sweep_epoch',

tests/test_sweep.py

@@ -6,8 +6,8 @@
 def sweep():
 
     i = [0,0,1,1,0,0,0,2,2,2,2,2,0,0,0,0]
-    v = [0,0,1,2,1,0,0,1,2,3,1,np.nan,np.nan,np.nan,np.nan,np.nan]
-    sampling_rate = 2
+    v = [0,0,1,2,1,0,0,1,2,3,1,1,np.nan,np.nan,np.nan,np.nan]
+    sampling_rate = 10
     dt = 1./sampling_rate
     t = np.arange(0,len(v))*dt
 
@@ -19,9 +19,9 @@ def test_select_epoch(sweep):
     i_sweep = sweep.i
     v_sweep = sweep.v
 
-    sweep.select_epoch("recording")
-    assert np.all(sweep.i == [0,0,1,1,0,0,0,2,2,2,2])
-    assert np.all(sweep.v == [0,0,1,2,1,0,0,1,2,3,1])
+    sweep.select_epoch("stim")
+    assert np.all(sweep.i == [2,2,2,2,2])
+    assert np.all(sweep.v == [1,2,3,1,1])
 
     sweep.select_epoch("sweep")
     assert np.all(sweep.i == i_sweep)