Bug fixes: AP_end_indices, AP_fall_indices, depolarized_base (#374)

* Fix index out of bounds error in __depolarized_base

* prevent division by zero and add bound check in depolarized_base

* fix vpeak construction to handle pi[i] > apei[i] and prevent out-of-bounds in AP_fall_indices

* fix bug in __AP_end_indices to avoid invalid min_element range

* update depolarized_base test

efel/cppcore/LibV2.cpp

@@ -76,6 +76,9 @@ static int __AP_fall_indices(const vector<double>& v, const vector<int>& apbi,
                              vector<int>& apfi) {
   apfi.resize(std::min(apbi.size(), pi.size()));
   for (size_t i = 0; i < apfi.size(); i++) {
+    if (pi[i] >= v.size() || apbi[i] >= v.size() || apei[i] >= v.size() || pi[i] > apei[i]) {
+        continue;
+    }
     double halfheight = (v[pi[i]] + v[apbi[i]]) / 2.;
     vector<double> vpeak(&v[pi[i]], &v[apei[i]]);
     transform(vpeak.begin(), vpeak.end(), vpeak.begin(),

efel/cppcore/LibV3.cpp

@@ -36,26 +36,25 @@ static int __depolarized_base(const vector<double>& t, const vector<double>& v,
   int i, n, k, startIndex, endIndex, nPt;
   double baseValue;
   // to make sure it access minimum index of both length
-  if (apendi.size() < apbi.size())
-    n = apendi.size();
-  else
-    n = apbi.size();
-
-  if (apendi.size() == apbi.size()) n = apendi.size() - 1;
+  n = std::min(apendi.size(), apbi.size());
 
   if (n > 2) {
     dep_base.clear();
-    for (i = 0; i < n; i++) {
+    for (i = 0; i < n - 1; i++) {
       nPt = 0;
       baseValue = 0;
       startIndex = apendi[i];
       endIndex = apbi[i + 1];
       for (k = startIndex; k < endIndex; k++) {
-        baseValue += v[k];
-        nPt++;
+        if (k >= 0 && k < v.size()) {
+          baseValue += v[k];
+          ++nPt;
+        }
+      }
+      if (nPt > 0) {
+        baseValue /= nPt;
+        dep_base.push_back(baseValue);
       }
-      baseValue = baseValue / nPt;
-      dep_base.push_back(baseValue);
     }
     return dep_base.size();
   }

efel/cppcore/LibV5.cpp

@@ -402,15 +402,22 @@ static int __AP_end_indices(const vector<double>& t, const vector<double>& v,
   peak_indices.push_back(v.size() - 1);
 
   for (size_t i = 0; i < peak_indices.size() - 1; i++) {
-    max_slope =
-        distance(dvdt.begin(), std::min_element(dvdt.begin() + peak_indices[i] + 1,
-                                                dvdt.begin() + peak_indices[i + 1]));
+    size_t start_index = peak_indices[i] + 1;
+    size_t end_index = peak_indices[i + 1];
+
+    if (start_index >= end_index || start_index >= dvdt.size() || end_index >= dvdt.size()) {
+      continue;
+    }
+
+    auto min_element_it = std::min_element(dvdt.begin() + start_index, dvdt.begin() + end_index);
+    auto max_slope = std::distance(dvdt.begin(), min_element_it);
     // assure that the width of the slope is bigger than 4
-    apei.push_back(distance(dvdt.begin(), find_if(dvdt.begin() + max_slope,
-                                             dvdt.begin() + peak_indices[i + 1],
-                                             [derivativethreshold](double x) {
-                                               return x >= derivativethreshold;
-                                             })));
+    auto threshold_it = std::find_if(dvdt.begin() + max_slope, dvdt.begin() + end_index,
+                                    [derivativethreshold](double x) { return x >= derivativethreshold; });
+
+    if (threshold_it != dvdt.begin() + end_index) {
+      apei.push_back(std::distance(dvdt.begin(), threshold_it));
+    }
   }
   return apei.size();
 }
@@ -1343,7 +1350,7 @@ double __decay_time_constant_after_stim(const vector<double>& times,
 
   if (decayTimes.size() < 1 || decayValues.size() < 1) {
     throw FeatureComputationError("No data points to calculate decay_time_constant_after_stim");
-  } 
+  }
   linear_fit_result fit;
   fit = slope_straight_line_fit(decayTimes, decayValues);
 

tests/test_basic.py

@@ -47,6 +47,7 @@
 
 testdata_dir = Path(__file__).parent / 'testdata'
 meanfrequency1_url = testdata_dir / 'basic' / 'mean_frequency_1.txt'
+meanfrequency2_url = testdata_dir / 'basic' / 'mean_frequency_2.txt'
 ahptest1_url = testdata_dir / 'basic' / 'ahptest_1.txt'
 tau20_0_url = testdata_dir / 'basic' / 'tau20.0.csv'
 spikeoutsidestim_url = testdata_dir / 'basic' / 'spike_outside_stim.txt'
@@ -72,6 +73,10 @@ def load_data(data_name, interp=False, interp_dt=0.1):
         stim_start = 500.0
         stim_end = 900.0
         time, voltage = load_ascii_input(meanfrequency1_url)
+    elif data_name == 'mean_frequency2':
+        stim_start = 500.0
+        stim_end = 900.0
+        time, voltage = load_ascii_input(meanfrequency2_url)
     elif data_name == 'tau20.0':
         stim_start = 100.0
         stim_end = 1000.0
@@ -2627,38 +2632,49 @@ def test_time_constant():
     numpy.testing.assert_allclose(time_cst, py_tau, rtol=1e-3)
 
 
-def test_depolarized_base():
-    """basic: Test depolarized base"""
+def assert_depolarized_base(trace_name, interp=True):
+    """
+    Calculates depolarized base using given trace data and checks for correctness.
+
+    :param trace_name: Name of the trace data to load.
+    :param interp: Whether to interpolate the data.
+    :return: None
+    """
 
     import efel
     efel.reset()
 
-    trace, time, voltage, stim_start, stim_end = load_data(
-        'mean_frequency1', interp=True)
+    trace, time, voltage, stim_start, stim_end = load_data(trace_name, interp=interp)
 
     features = ["depolarized_base", "AP_begin_time", "AP_duration"]
-
-    feature_values = \
-        get_feature_values(
-            [trace],
-            features, raise_warnings=False)
+    feature_values = get_feature_values([trace], features, raise_warnings=False)
 
     depolarized_base = feature_values[0]['depolarized_base']
     AP_begin_times = feature_values[0]['AP_begin_time']
     AP_durations = feature_values[0]['AP_duration']
 
     py_dep_base = []
-    for i, (AP_begin, AP_dur) in enumerate(
-        zip(AP_begin_times[:-1], AP_durations[:-1])
-    ):
+    for i, (AP_begin, AP_dur) in enumerate(zip(AP_begin_times[:-1], AP_durations[:-1])):
         dep_start_time = AP_begin + AP_dur
         dep_end_time = AP_begin_times[i + 1]
+        if dep_end_time <= dep_start_time:
+            continue
         start_idx = numpy.argwhere(time > dep_start_time)[0][0] - 1
         end_idx = numpy.argwhere(time > dep_end_time)[0][0] - 1
 
         py_dep_base.append(numpy.mean(voltage[start_idx:end_idx]))
 
-    numpy.testing.assert_allclose(depolarized_base, py_dep_base)
+    numpy.testing.assert_allclose(depolarized_base, py_dep_base, rtol=1e-3)
+
+
+def test_depolarized_base():
+    """Test depolarized base with standard data."""
+    assert_depolarized_base('mean_frequency1', interp=True)
+
+
+def test_depolarized_base_outlier():
+    """Test depolarized base with outlier data."""
+    assert_depolarized_base('mean_frequency2', interp=False)
 
 
 def test_AP_duration():