Merge pull request BlueBrain#275 from AurelienJaquier/postburst_AHP

new feature: postburst_min_values

docs/source/eFeatures.rst

@@ -502,6 +502,36 @@ The first spike is ignored by default. This can be changed by setting ignore_fir
         ) for i in burst_end_indices if i + 1 < len(peak_indices)
     ]
 
+LibV5 : postburst_min_values
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The minimum voltage after the end of a burst.
+
+This implementation does not assume that every spike belongs to a burst.
+
+The first spike is ignored by default. This can be changed by setting ignore_first_ISI to 0.
+
+- **Required features**: peak_indices, burst_end_indices
+- **Units**: mV
+- **Pseudocode**: ::
+
+    interburst_min = [
+        numpy.min(
+            v[peak_indices[i]:peak_indices[i + 1]]
+        ) for i in burst_end_indices if i + 1 < len(peak_indices)
+    ]
+
+    if len(postburst_min) < len(burst_end_indices):
+        if t[burst_end_indices[-1]] < stim_end:
+            end_idx = numpy.where(t >= stim_end)[0][0]
+            postburst_min.append(numpy.min(
+                v[peak_indices[burst_end_indices[-1]]:end_idx]
+            ))
+        else:
+            postburst_min.append(numpy.min(
+                v[peak_indices[burst_end_indices[-1]]:]
+            ))
+
 LibV5 : time_to_interburst_min
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

efel/DependencyV5.txt

@@ -172,4 +172,6 @@ LibV5:ADP_peak_values     #LibV5:ADP_peak_indices #LibV1:interpolate
 LibV5:ADP_peak_amplitude  #LibV5:ADP_peak_values #LibV5:min_AHP_values #LibV1:interpolate
 LibV5:interburst_min_indices     #LibV5:peak_indices #LibV5:burst_end_indices #LibV1:interpolate
 LibV5:interburst_min_values      #LibV5:interburst_min_indices #LibV1: interpolate
+LibV5:postburst_min_indices     #LibV5:peak_indices #LibV5:burst_end_indices #LibV1:interpolate
+LibV5:postburst_min_values      #LibV5:postburst_min_indices #LibV1: interpolate
 LibV5:time_to_interburst_min     #LibV5:interburst_min_indices #LibV1:peak_time #LibV5:burst_end_indices #LibV1:interpolate

efel/cppcore/FillFptrTable.cpp

@@ -285,6 +285,8 @@ int FillFptrTable() {
 
   FptrTableV5["interburst_min_indices"] = &LibV5::interburst_min_indices;
   FptrTableV5["interburst_min_values"] = &LibV5::interburst_min_values;
+  FptrTableV5["postburst_min_indices"] = &LibV5::postburst_min_indices;
+  FptrTableV5["postburst_min_values"] = &LibV5::postburst_min_values;
   FptrTableV5["time_to_interburst_min"] = &LibV5::time_to_interburst_min;
 
   //****************** end of FptrTableV5 *****************************

efel/cppcore/LibV5.cpp

@@ -3863,6 +3863,98 @@ int LibV5::interburst_min_values(mapStr2intVec& IntFeatureData,
 }
 
 
+static int __postburst_min_indices(vector<double>& t,
+                                   vector<double>& v,
+                                   vector<int>& peak_indices,
+                                   vector<int>& burst_end_indices,
+                                   vector<int>& postburst_min_indices,
+                                   vector<double>& postburst_min_values,
+                                   const double stim_end) {
+  unsigned postburst_min_index, stim_end_index, end_index;
+  stim_end_index =
+      distance(t.begin(),
+                find_if(t.begin(), t.end(),
+                        std::bind2nd(std::greater_equal<double>(), stim_end)));
+  end_index = distance(t.begin(), t.end());
+  for (size_t i = 0; i < burst_end_indices.size(); i++) {
+    if (burst_end_indices[i] + 1 < peak_indices.size()){
+      postburst_min_index = min_element(
+        v.begin() + peak_indices[burst_end_indices[i]],
+        v.begin() + peak_indices[burst_end_indices[i] + 1]
+      ) - v.begin();
+    } else if (peak_indices[burst_end_indices[i]] < stim_end_index){
+      postburst_min_index = min_element(
+        v.begin() + peak_indices[burst_end_indices[i]],
+        v.begin() + stim_end_index
+      ) - v.begin();
+    } else {
+      postburst_min_index = min_element(
+        v.begin() + peak_indices[burst_end_indices[i]],
+        v.begin() + end_index
+      ) - v.begin();
+    }
+
+    postburst_min_indices.push_back(postburst_min_index);
+    postburst_min_values.push_back(v[postburst_min_index]);
+  }
+
+  return postburst_min_indices.size();
+}
+
+int LibV5::postburst_min_indices(mapStr2intVec& IntFeatureData,
+                                  mapStr2doubleVec& DoubleFeatureData,
+                                  mapStr2Str& StringData) {
+  int retVal, nSize;
+  retVal = CheckInMap(IntFeatureData, StringData,
+                      "postburst_min_indices", nSize);
+  if (retVal)
+    return nSize;
+
+  double stim_end;
+  vector<int> peak_indices, burst_end_indices, postburst_min_indices;
+  vector<double> postburst_min_values, t, v, stim_end_vec;
+  retVal = getVec(DoubleFeatureData, StringData, "T", t);
+  if (retVal <= 0) return -1;
+  retVal = getVec(DoubleFeatureData, StringData, "V", v);
+  if (retVal <= 0) return -1;
+  retVal = getVec(IntFeatureData, StringData, "peak_indices",
+                  peak_indices);
+  if (retVal < 0) return -1;
+  retVal = getVec(IntFeatureData, StringData, "burst_end_indices",
+                  burst_end_indices);
+  if (retVal < 0) return -1;
+  retVal =
+      getVec(DoubleFeatureData, StringData, "stim_end", stim_end_vec);
+  if (retVal <= 0) {
+    return -1;
+  } else {
+    stim_end = stim_end_vec[0];
+  }
+
+  retVal = __postburst_min_indices(t, v, peak_indices, burst_end_indices,
+                                   postburst_min_indices,
+                                   postburst_min_values,
+                                   stim_end);
+  if (retVal >= 0) {
+    setVec(IntFeatureData, StringData, "postburst_min_indices",
+           postburst_min_indices);
+    setVec(DoubleFeatureData, StringData, "postburst_min_values",
+           postburst_min_values);
+  }
+  return retVal;
+}
+
+int LibV5::postburst_min_values(mapStr2intVec& IntFeatureData,
+                                 mapStr2doubleVec& DoubleFeatureData,
+                                 mapStr2Str& StringData) {
+  int retVal, nSize;
+  retVal =
+      CheckInMap(DoubleFeatureData, StringData, "postburst_min_values", nSize);
+  if (retVal >= 0) return nSize;
+  return -1;
+}
+
+
 int LibV5::time_to_interburst_min(mapStr2intVec& IntFeatureData,
                                   mapStr2doubleVec& DoubleFeatureData,
                                   mapStr2Str& StringData) {

efel/cppcore/LibV5.h

@@ -290,6 +290,12 @@ int interburst_min_indices(mapStr2intVec& IntFeatureData,
 int interburst_min_values(mapStr2intVec& IntFeatureData,
                             mapStr2doubleVec& DoubleFeatureData,
                             mapStr2Str& StringData);
+int postburst_min_indices(mapStr2intVec& IntFeatureData,
+                            mapStr2doubleVec& DoubleFeatureData,
+                            mapStr2Str& StringData);
+int postburst_min_values(mapStr2intVec& IntFeatureData,
+                            mapStr2doubleVec& DoubleFeatureData,
+                            mapStr2Str& StringData);
 int time_to_interburst_min(mapStr2intVec& IntFeatureData,
                             mapStr2doubleVec& DoubleFeatureData,
                             mapStr2Str& StringData);

efel/cppcore/cfeature.cpp

@@ -264,6 +264,8 @@ void cFeature::fillfeaturetypes() {
 
   featuretypes["interburst_min_indices"] = "int";
   featuretypes["interburst_min_values"] = "double";
+  featuretypes["postburst_min_indices"] = "int";
+  featuretypes["postburst_min_values"] = "double";
   featuretypes["time_to_interburst_min"] = "double";
 
   // end of feature types

efel/tests/featurenames.json

@@ -181,6 +181,8 @@
     "ADP_peak_amplitude",
     "interburst_min_indices",
     "interburst_min_values",
+    "postburst_min_indices",
+    "postburst_min_values",
     "time_to_interburst_min",
     "AHP_depth_slow"
 ]

efel/tests/test_basic.py

@@ -3648,6 +3648,91 @@ def test_time_to_interburst_min():
         )
 
 
+def py_postburst_min_values(t, v, peak_indices, burst_end_indices, stim_end):
+    """min voltage between burst and next peak"""
+    if peak_indices is None or burst_end_indices is None:
+        return None
+
+    postburst_min = [
+        numpy.min(
+            v[peak_indices[i]:peak_indices[i + 1]]
+        ) for i in burst_end_indices if i + 1 < len(peak_indices)
+    ]
+    if len(postburst_min) < len(burst_end_indices):
+        if t[burst_end_indices[-1]] < stim_end:
+            end_idx = numpy.where(t >= stim_end)[0][0]
+            postburst_min.append(numpy.min(
+                v[peak_indices[burst_end_indices[-1]]:end_idx]
+            ))
+        else:
+            postburst_min.append(numpy.min(
+                v[peak_indices[burst_end_indices[-1]]:]
+            ))
+
+    return numpy.array(postburst_min)
+
+
+def test_postburst_min_values():
+    """basic: Test interburst_min_values"""
+    urls = [burst1_url, burst2_url, burst3_url]
+    for i, url in enumerate(urls):
+        import efel
+        efel.reset()
+        # use this to have all spikes in burst for burst3_url case
+        efel.setDoubleSetting('strict_burst_factor', 4.0)
+
+        time = efel.io.load_fragment('%s#col=1' % url)
+        voltage = efel.io.load_fragment('%s#col=2' % url)
+
+        interp_time, interp_voltage = interpolate(time, voltage, 0.1)
+
+        trace = {}
+
+        trace['T'] = time
+        trace['V'] = voltage
+        if i in [0, 1]:
+            trace['stim_start'] = [250]
+            trace['stim_end'] = [1600]
+        elif i == 2:
+            trace['stim_start'] = [800]
+            trace['stim_end'] = [2150]
+
+        features = [
+            "peak_indices",
+            "burst_end_indices",
+            "postburst_min_values",
+        ]
+
+        feature_values = efel.getFeatureValues(
+            [trace],
+            features,
+            raise_warnings=False
+        )
+
+        peak_indices = feature_values[0]["peak_indices"]
+        burst_end_indices = feature_values[0]["burst_end_indices"]
+        postburst_min_values = feature_values[0]["postburst_min_values"]
+
+        postburst_min_py = py_postburst_min_values(
+            interp_time,
+            interp_voltage,
+            peak_indices,
+            burst_end_indices,
+            trace["stim_end"][0]
+        )
+
+        # convert to float so that None edge case get converted to nan
+        # and can pass in assert_allclose
+        postburst_min_py = numpy.array(postburst_min_py, dtype=numpy.float64)
+        postburst_min_values = numpy.array(
+            postburst_min_values, dtype=numpy.float64
+        )
+
+        numpy.testing.assert_allclose(
+            postburst_min_py, postburst_min_values
+        )
+
+
 def test_spikes_per_burst_diff():
     """basic: Test spikes_per_burst_diff"""
     import efel

efel/tests/testdata/allfeatures/expectedresults.json

@@ -60679,6 +60679,12 @@
     ],
     "interburst_min_indices": [],
     "interburst_min_values": [],
+    "postburst_min_indices": [
+        26465
+    ],
+    "postburst_min_values": [
+        -41.52919675004766
+    ],
     "time_to_interburst_min": [],
     "AHP_depth_slow": [
         28.81047723194662,