fix: bugfixes

Author: jnsbck

ephyspy/allen_sdk/ephys_features.py

@@ -203,10 +203,12 @@ def find_upstroke_indexes(v, t, spike_indexes, peak_indexes, filter=10.0, dvdt=N
     if dvdt is None:
         dvdt = calculate_dvdt(v, t, filter)
 
-    upstroke_indexes = [
-        np.argmax(dvdt[spike:peak]) + spike
-        for spike, peak in zip(spike_indexes, peak_indexes)
-    ]
+    upstroke_indexes = []
+    for spike, peak in zip(spike_indexes, peak_indexes):
+        if len(dvdt[spike:peak]) > 0:
+            upstroke_indexes.append(np.argmax(dvdt[spike:peak]) + spike)
+        else:
+            upstroke_indexes.append(np.nan)
 
     return np.array(upstroke_indexes)
 
@@ -985,7 +987,8 @@ def detect_pauses(isis, isi_types, cost_weight=1.0):
             break
         cv = non_pause_isis.std() / non_pause_isis.mean()
         benefit = all_cv - cv
-        cost = np.sum(non_pause_isis.std() / np.abs(non_pause_isis.mean() - pause_isis))
+        div = non_pause_isis.mean() - pause_isis
+        cost = np.sum(non_pause_isis.std() / np.abs(div if div != 0 else float("nan")))
         cost *= cost_weight
         net = benefit - cost
         if net > 0 and net < best_net:

ephyspy/features/base.py

@@ -291,6 +291,10 @@ def lookup_spike_feature(
         ):
             sweep.process_spikes()
 
+        # Allen dataframe does not always have all columns even though the features
+        # in theory could be computed.
+        if is_allen_ft and feature_name not in sweep._spikes_df.columns:
+            return np.array([])
         return sweep.spike_feature(feature_name, include_clipped=True)
 
     def get_value(

ephyspy/features/spike_features.py

@@ -552,9 +552,11 @@ def __init__(self, data=None, **kwargs):
 
     def _compute(self, recompute=False, store_diagnostics=True):
         width = self.lookup_spike_feature("width", recompute=recompute)
-        trough_idxs = self.lookup_spike_feature("trough_index").astype(int)
-        spike_idxs = self.lookup_spike_feature("threshold_index").astype(int)
-        peak_idxs = self.lookup_spike_feature("peak_index").astype(int)
+
+        # values are floats since they can be nan
+        trough_idxs = self.lookup_spike_feature("trough_index")
+        spike_idxs = self.lookup_spike_feature("threshold_index")
+        peak_idxs = self.lookup_spike_feature("peak_index")
 
         if store_diagnostics:
             self._update_diagnostics(
@@ -574,6 +576,9 @@ def _plot(self, ax: Optional[Axes] = None, selected_idxs=None, **kwargs) -> Axes
             trough_idxs, spike_idxs, peak_idxs = unpack(
                 self.diagnostics, ["trough_idx", "spike_idx", "peak_idx"]
             )
+            trough_idxs = trough_idxs.astype(int)
+            spike_idxs = spike_idxs.astype(int)
+            peak_idxs = peak_idxs.astype(int)
 
             t = self.data.t
             v = self.data.v

ephyspy/features/sweep_features.py

@@ -587,6 +587,11 @@ def _compute(self, recompute=False, store_diagnostics=True):
         isi_ff = float("nan")
         if has_spikes(self.data):
             isi = self.lookup_spike_feature("isi", recompute=recompute)[1:]
+            during_stim = where_spike_during_stimulus(self, recompute=recompute)[1:]
+            if np.any(during_stim):
+                isi = isi[during_stim]
+                isi = isi[1:] if isi[0] == 0 else isi
+
             if len(isi) > 1:
                 isi_ff = np.nanvar(isi) / np.nanmean(isi)
 
@@ -619,6 +624,11 @@ def _compute(self, recompute=False, store_diagnostics=True):
         isi_cv = float("nan")
         if has_spikes(self.data):
             isi = self.lookup_spike_feature("isi", recompute=recompute)[1:]
+            during_stim = where_spike_during_stimulus(self, recompute=recompute)[1:]
+            if np.any(during_stim):
+                isi = isi[during_stim]
+                isi = isi[1:] if isi[0] == 0 else isi
+
             if len(isi) > 1:
                 isi_cv = np.nanstd(isi) / np.nanmean(isi)
 
@@ -1513,15 +1523,16 @@ def __init__(self, data=None, compute_at_init=True, **kwargs):
     def _compute(self, recompute=False, store_diagnostics=True):
         isi_adapt = float("nan")
         if has_spikes(self.data):
-            isi = self.lookup_spike_feature("isi", recompute=recompute)
-            during_stim = where_spike_during_stimulus(self, recompute=recompute)
-            isi = isi[during_stim]
-            isi = isi[1:] if isi[0] == 0 else isi
-            if len(isi) > 1:
-                isi_adapt = isi[1] / isi[0]
+            isi = self.lookup_spike_feature("isi", recompute=recompute)[1:]
+            during_stim = where_spike_during_stimulus(self, recompute=recompute)[1:]
+            if np.any(during_stim):
+                isi = isi[during_stim]
+                isi = isi[1:] if isi[0] == 0 else isi
+                if len(isi) > 1:
+                    isi_adapt = isi[1] / isi[0]
 
-            if store_diagnostics:
-                self._update_diagnostics({"isi": isi})
+                if store_diagnostics:
+                    self._update_diagnostics({"isi": isi})
         return isi_adapt
 
     def _plot(self, ax: Optional[Axes] = None, **kwargs) -> Axes:
@@ -1544,16 +1555,17 @@ def __init__(self, data=None, compute_at_init=True, **kwargs):
     def _compute(self, recompute=False, store_diagnostics=True):
         isi_adapt_avg = float("nan")
         if has_spikes(self.data):
-            isi = self.lookup_spike_feature("isi", recompute=recompute)
-            during_stim = where_spike_during_stimulus(self, recompute=recompute)
-            isi = isi[during_stim]
-            isi = isi[1:] if isi[0] == 0 else isi
-            if len(isi) > 2:
-                isi_changes = isi[1:] / isi[:-1]
-                isi_adapt_avg = isi_changes.mean()
+            isi = self.lookup_spike_feature("isi", recompute=recompute)[1:]
+            during_stim = where_spike_during_stimulus(self, recompute=recompute)[1:]
+            if np.any(during_stim):
+                isi = isi[during_stim]
+                isi = isi[1:] if isi[0] == 0 else isi
+                if len(isi) > 2:
+                    isi_changes = isi[1:] / isi[:-1]
+                    isi_adapt_avg = isi_changes.mean()
 
-                if store_diagnostics:
-                    self._update_diagnostics({"isi": isi})
+                    if store_diagnostics:
+                        self._update_diagnostics({"isi": isi})
         return isi_adapt_avg
 
     def _plot(self, ax: Optional[Axes] = None, **kwargs) -> Axes:

ephyspy/features/sweepset_features.py

@@ -99,7 +99,8 @@ def _select(self, fts):
             peaks_to_low = np.all(self.lookup_sweep_feature("ap_peak") < -30)
 
             if not peaks_to_low:
-                idx = np.where(is_depol & has_spikes)[0][0]
+                idx = np.where(is_depol & has_spikes)[0]
+                idx = idx[0] if len(idx) > 0 else np.array([], dtype=int)
             else:
                 idx = np.array([], dtype=int)
         elif stimulus_type(self.data) == "ramp":
@@ -199,8 +200,9 @@ def _select(self, fts):
         num_spikes = self.lookup_sweep_feature("num_ap")
         wildness = self.lookup_sweep_feature("wildness")
         is_non_wild = np.isnan(wildness)
-        idx = pd.Series(num_spikes)[is_non_wild].idxmax()
-        idx = np.array([], dtype=int) if np.isnan(idx) else idx
+        num_non_wild_spikes = pd.Series(num_spikes)[is_non_wild]
+        no_spikes = np.all(np.isnan(num_non_wild_spikes))
+        idx = np.array([], dtype=int) if no_spikes else num_non_wild_spikes.idxmax()
 
         self._update_diagnostics(
             {
@@ -459,31 +461,33 @@ def _compute(self, recompute=False, store_diagnostics=False):
             i = stim_amp[is_depol]
 
             has_spikes = ~np.isnan(f)
-            # sometimes all depolarization traces spike
-            i_sub = (
-                0 if all(has_spikes) else i[~has_spikes][0]
-            )  # last stim < spike threshold
-            i_sup = i[has_spikes][0]  # first stim > spike threshold
+            if np.any(has_spikes):
+                # sometimes all depolarization traces spike
+                i_sub = (
+                    0 if all(has_spikes) else i[~has_spikes][0]
+                )  # last stim < spike threshold
+                i_sup = i[has_spikes][0]  # first stim > spike threshold
+
+                if not np.isnan(dfdi):
+                    rheobase = float(ransac.predict(np.array([[0]]))) / dfdi
+
+                    if rheobase < i_sub or rheobase > i_sup:
+                        rheobase = i_sup
+                else:
+                    rheobase = i_sup
+                rheobase -= dc_offset
 
-            if not np.isnan(dfdi):
-                rheobase = float(ransac.predict(np.array([[0]]))) / dfdi
+                if store_diagnostics:
+                    self._update_diagnostics(
+                        {
+                            "i_sub": i_sub,
+                            "i_sup": i_sup,
+                            "f_sup": f[has_spikes][0],
+                            "dfdi": dfdi,
+                            "dc_offset": dc_offset,
+                        }
+                    )
 
-                if rheobase < i_sub or rheobase > i_sup:
-                    rheobase = i_sup
-            else:
-                rheobase = i_sup
-            rheobase -= dc_offset
-
-            if store_diagnostics:
-                self._update_diagnostics(
-                    {
-                        "i_sub": i_sub,
-                        "i_sup": i_sup,
-                        "f_sup": f[has_spikes][0],
-                        "dfdi": dfdi,
-                        "dc_offset": dc_offset,
-                    }
-                )
         if stimulus_type(self.data) == "ramp":
             has_ap = self.lookup_sweep_feature("num_ap", recompute=recompute) > 0
             if np.any(has_ap):
@@ -685,19 +689,22 @@ def _compute(self, recompute=False, store_diagnostics=False):
         slow_hyperpolarization = float("nan")
         if stimulus_type(self.data) == "long_square":
             has_aps = self.lookup_sweep_feature("num_ap", recompute=recompute) > 0
-            v_baseline = self.lookup_sweep_feature("v_baseline", recompute=recompute)
-            v_baseline = v_baseline[has_aps]
-
-            slow_hyperpolarization = v_baseline.max() - v_baseline.min()
-
-            if store_diagnostics:
-                self._update_diagnostics(
-                    {
-                        "v_baseline": v_baseline,
-                        "v_baseline_max": v_baseline.max(),
-                        "v_baseline_min": v_baseline.min(),
-                    }
+            if np.any(has_aps):
+                v_baseline = self.lookup_sweep_feature(
+                    "v_baseline", recompute=recompute
                 )
+                v_baseline = v_baseline[has_aps]
+
+                slow_hyperpolarization = v_baseline.max() - v_baseline.min()
+
+                if store_diagnostics:
+                    self._update_diagnostics(
+                        {
+                            "v_baseline": v_baseline,
+                            "v_baseline_max": v_baseline.max(),
+                            "v_baseline_min": v_baseline.min(),
+                        }
+                    )
         return slow_hyperpolarization
 
     def _plot(self, ax: Optional[Axes] = None, **kwargs) -> Axes:
@@ -730,26 +737,29 @@ def _compute(self, recompute=False, store_diagnostics=False):
         slow_hyperpolarization_slope = float("nan")
         if stimulus_type(self.data) == "long_square":
             has_aps = self.lookup_sweep_feature("num_ap", recompute=recompute) > 0
-            v_baseline = self.lookup_sweep_feature("v_baseline", recompute=recompute)
-            v_baseline = v_baseline[has_aps]
-
-            v_baseline = v_baseline.reshape(-1, 1)
-            sweep_idx = np.arange(len(v_baseline)).reshape(-1, 1)
-
-            if len(v_baseline) >= 3:
-                ransac.fit(sweep_idx, v_baseline)
-                slope = ransac.coef_[0, 0] * 1000
-                intercept = ransac.intercept_[0]
-                slow_hyperpolarization_slope = slope
-
-            if store_diagnostics:
-                self._update_diagnostics(
-                    {
-                        "v_baseline": v_baseline,
-                        "sweep_idx": sweep_idx,
-                        "v_intercept": intercept,
-                    }
+            if np.any(has_aps):
+                v_baseline = self.lookup_sweep_feature(
+                    "v_baseline", recompute=recompute
                 )
+                v_baseline = v_baseline[has_aps]
+
+                v_baseline = v_baseline.reshape(-1, 1)
+                sweep_idx = np.arange(len(v_baseline)).reshape(-1, 1)
+
+                if len(v_baseline) >= 3:
+                    ransac.fit(sweep_idx, v_baseline)
+                    slope = ransac.coef_[0, 0] * 1000
+                    intercept = ransac.intercept_[0]
+                    slow_hyperpolarization_slope = slope
+
+                    if store_diagnostics:
+                        self._update_diagnostics(
+                            {
+                                "v_baseline": v_baseline,
+                                "sweep_idx": sweep_idx,
+                                "v_intercept": intercept,
+                            }
+                        )
         return slow_hyperpolarization_slope
 
     def _plot(self, ax: Optional[Axes] = None, **kwargs) -> Axes: