Merge pull request #196 from BlueBrain/ap-end

add AP end indices with a variable derivative threshold in LibV5 #88
BlueBrain · Jan 21, 2021 · 83d27ff · 83d27ff
2 parents fee2a0b + 12fab4b
commit 83d27ff
Show file tree

Hide file tree

Showing 8 changed files with 145 additions and 7 deletions.
diff --git a/.gitignore b/.gitignore
@@ -9,3 +9,6 @@
 /.python-version
 /.tox
 /.coverage
+.vscode
+*.so
+coverage.xml
diff --git a/efel/DependencyV5.txt b/efel/DependencyV5.txt
@@ -33,15 +33,15 @@ LibV1:ISI_CV	#LibV1:ISI_values #LibV1:interpolate
 LibV1:Spikecount	#LibV5:peak_indices #LibV1:interpolate 
 LibV5:Spikecount_stimint	#LibV1:peak_time #LibV1:interpolate 
 LibV1:AHP_depth	#LibV5:voltage_base	#LibV5:min_AHP_values #LibV1:interpolate 
-LibV2:AP_rise_indices       #LibV5:peak_indices     #LibV5:AP_begin_indices #LibV1:interpolate 
-LibV2:AP_end_indices        #LibV5:peak_indices #LibV1:interpolate 
-LibV2:AP_fall_indices       #LibV5:peak_indices     #LibV5:AP_begin_indices     #LibV2:AP_end_indices #LibV1:interpolate 
-LibV2:AP_duration	        #LibV5:AP_begin_indices #LibV2:AP_end_indices #LibV1:interpolate 
+LibV2:AP_rise_indices       #LibV5:peak_indices     #LibV5:AP_begin_indices #LibV1:interpolate
+LibV5:AP_end_indices        #LibV5:peak_indices #LibV1:interpolate
+LibV2:AP_fall_indices       #LibV5:peak_indices     #LibV5:AP_begin_indices     #LibV5:AP_end_indices #LibV1:interpolate 
+LibV2:AP_duration	        #LibV5:AP_begin_indices #LibV5:AP_end_indices #LibV1:interpolate 
 LibV2:AP_duration_half_width	#LibV2:AP_rise_indices	#LibV2:AP_fall_indices #LibV1:interpolate 
 LibV2:AP_rise_time	#LibV5:AP_begin_indices	#LibV5:peak_indices #LibV1:interpolate 
-LibV2:AP_fall_time	#LibV5:peak_indices	#LibV2:AP_end_indices #LibV1:interpolate 
+LibV2:AP_fall_time	#LibV5:peak_indices	#LibV5:AP_end_indices #LibV1:interpolate 
 LibV2:AP_rise_rate	#LibV5:AP_begin_indices	#LibV5:peak_indices #LibV1:interpolate 
-LibV2:AP_fall_rate	#LibV5:peak_indices	#LibV2:AP_end_indices #LibV1:interpolate 
+LibV2:AP_fall_rate	#LibV5:peak_indices	#LibV5:AP_end_indices #LibV1:interpolate 
 LibV2:fast_AHP	#LibV5:AP_begin_indices	#LibV5:min_AHP_indices #LibV1:interpolate 
 LibV2:AP_amplitude_change	#LibV1:AP_amplitude #LibV1:interpolate 
 LibV2:AP_duration_change	#LibV2:AP_duration #LibV1:interpolate 
@@ -119,7 +119,7 @@ LibV5:inv_fifth_ISI   #LibV5:all_ISI_values #LibV1:interpolate
 LibV5:inv_last_ISI   #LibV5:all_ISI_values #LibV1:interpolate 
 LibV5:inv_time_to_first_spike   #LibV1:time_to_first_spike #LibV1:interpolate 
 LibV5:spike_half_width 	    #LibV5:min_AHP_indices	#LibV5:peak_indices #LibV1:interpolate 
-LibV5:AP_begin_indices      #LibV5:min_AHP_indices #LibV1:interpolate 
+LibV5:AP_begin_indices      #LibV5:min_AHP_indices #LibV1:interpolate
 LibV5:AHP_depth_abs         #LibV5:min_AHP_values #LibV1:interpolate 
 LibV5:AP_begin_width 	    #LibV5:min_AHP_indices	#LibV5:AP_begin_indices #LibV1:interpolate 
 LibV5:AP_begin_voltage  #LibV5:AP_begin_indices #LibV1:interpolate 

diff --git a/efel/api.py b/efel/api.py
@@ -77,6 +77,9 @@ def reset():
     setIntSetting('max_spike_skip', 2)
     setDoubleSetting('Threshold', _settings.threshold)
     setDoubleSetting('DerivativeThreshold', _settings.derivative_threshold)
+    setDoubleSetting(
+        'DownDerivativeThreshold',
+        _settings.down_derivative_threshold)
     setDoubleSetting('interp_step', 0.1)
     setDoubleSetting('burst_factor', 1.5)
     setDoubleSetting('voltage_base_start_perc', 0.9)

diff --git a/efel/cppcore/FillFptrTable.cpp b/efel/cppcore/FillFptrTable.cpp
@@ -191,6 +191,7 @@ int FillFptrTable() {
   FptrTableV5["AHP_depth_abs"] = &LibV5::AHP_depth_abs;
   FptrTableV5["spike_half_width"] = &LibV5::spike_width1;
   FptrTableV5["AP_begin_indices"] = &LibV5::AP_begin_indices;
+  FptrTableV5["AP_end_indices"] = &LibV5::AP_end_indices;
   FptrTableV5["irregularity_index"] = &LibV5::irregularity_index;
   FptrTableV5["number_initial_spikes"] = &LibV5::number_initial_spikes;
   FptrTableV5["AP1_amp"] = &LibV5::AP1_amp;

diff --git a/efel/cppcore/LibV5.cpp b/efel/cppcore/LibV5.cpp
@@ -751,6 +751,72 @@ int LibV5::AP_begin_indices(mapStr2intVec& IntFeatureData,
   }
   return retVal;
 }
+
+static int __AP_end_indices(const vector<double>& t, const vector<double>& v,
+                            const vector<int>& pi, vector<int>& apei,
+                            double derivativethreshold) {
+
+  vector<double> dvdt(v.size());
+  vector<double> dv;
+  vector<double> dt;
+  getCentralDifferenceDerivative(1., v, dv);
+  getCentralDifferenceDerivative(1., t, dt);
+  transform(dv.begin(), dv.end(), dt.begin(), dvdt.begin(),
+            std::divides<double>());
+
+  apei.resize(pi.size());
+  vector<int> picopy(pi.begin(), pi.end());
+  picopy.push_back(v.size() - 1);
+
+  for (size_t i = 0; i < apei.size(); i++) {
+    // assure that the width of the slope is bigger than 4
+    apei[i] = std::distance(
+        dvdt.begin(),
+        std::find_if(dvdt.begin() + picopy[i] + 1, dvdt.begin() + picopy[i + 1],
+                std::bind2nd(std::greater_equal<double>(), derivativethreshold)));
+  }
+  return apei.size();
+}
+
+
+int LibV5::AP_end_indices(mapStr2intVec& IntFeatureData,
+                          mapStr2doubleVec& DoubleFeatureData,
+                          mapStr2Str& StringData) {
+
+  int retVal;
+  int nSize;
+  retVal = CheckInMap(IntFeatureData, StringData, "AP_end_indices", nSize);
+  if (retVal) {
+    return nSize;
+  }
+
+  vector<double> t;
+  retVal = getVec(DoubleFeatureData, StringData, "T", t);
+  if (retVal < 0) return -1;
+  vector<double> v;
+  retVal = getVec(DoubleFeatureData, StringData, "V", v);
+  if (retVal < 0) return -1;
+  vector<int> pi;
+  retVal = getVec(IntFeatureData, StringData, "peak_indices", pi);
+  if (retVal < 0) return -1;
+
+  // Get DerivativeThreshold
+  vector<double> dTh;
+  retVal = getDoubleParam(DoubleFeatureData, "DownDerivativeThreshold", dTh);
+  if (retVal <= 0) {
+    // derivative at peak end
+    dTh.push_back(-12.0);
+  }
+
+  vector<int> apei;
+  retVal = __AP_end_indices(t, v, pi, apei, dTh[0]);
+  if (retVal >= 0) {
+    setVec(IntFeatureData, StringData, "AP_end_indices", apei);
+  }
+  return retVal;
+}
+
+
 static int __irregularity_index(vector<double>& isiValues,
                                 vector<double>& irregularity_index) {
   double ISISub, iRI;

diff --git a/efel/cppcore/LibV5.h b/efel/cppcore/LibV5.h
@@ -74,6 +74,10 @@ int AP_begin_indices(mapStr2intVec& IntFeatureData,
                      mapStr2doubleVec& DoubleFeatureData,
                      mapStr2Str& StringData);
 
+int AP_end_indices(mapStr2intVec& IntFeatureData,
+                   mapStr2doubleVec& DoubleFeatureData,
+                   mapStr2Str& StringData);
+
 int irregularity_index(mapStr2intVec& IntFeatureData,
                        mapStr2doubleVec& DoubleFeatureData,
                        mapStr2Str& StringData);

diff --git a/efel/settings.py b/efel/settings.py
@@ -38,6 +38,7 @@ class Settings(object):
     def __init__(self):
         self.threshold = -20.0
         self.derivative_threshold = 10.0
+        self.down_derivative_threshold = -12.0
         self.dependencyfile_path = os.path.join(
             _get_script_path(),
             'DependencyV5.txt')
diff --git a/efel/tests/test_basic.py b/efel/tests/test_basic.py
@@ -778,6 +778,66 @@ def test_AP_begin_indices1():
         len(feature_values[0]['AP_duration_half_width']))
 
 
+def test_AP_end_indices():
+    """basic: Test AP end indices."""
+    import efel
+    efel.reset()
+
+    stim_start = 31.2
+    stim_end = 431.2
+
+    test_data_path = os.path.join(
+        testdata_dir,
+        'basic',
+        'AP_begin_indices_95810005.abf.csv')
+    voltage = numpy.loadtxt(test_data_path)
+
+    time = numpy.arange(len(voltage)) * 0.1
+
+    trace = {}
+
+    trace['V'] = voltage
+    trace['T'] = time
+    trace['stim_start'] = [stim_start]
+    trace['stim_end'] = [stim_end]
+
+    features = [
+        'AP_begin_indices',
+        'peak_indices',
+        'AP_end_indices']
+
+    feature_values = \
+        efel.getFeatureValues(
+            [trace],
+            features,
+            raise_warnings=False)
+
+    begin_indices = feature_values[0]["AP_begin_indices"]
+    peak_indices = feature_values[0]["peak_indices"]
+    end_indices = feature_values[0]["AP_end_indices"]
+
+    for begin, peak, end in zip(begin_indices, peak_indices, end_indices):
+        # the voltage value for the end index should be closer than that of
+        # begin index than the peak
+        nt.assert_true(abs(voltage[begin] - voltage[end])
+                       < abs(voltage[peak] - voltage[end]))
+        nt.assert_true(end > begin)
+
+    efel.reset()
+
+    efel.setDoubleSetting("DownDerivativeThreshold", -24)
+    feature_values = \
+        efel.getFeatureValues(
+            [trace],
+            features,
+            raise_warnings=False)
+
+    updated_end_indices = feature_values[0]["AP_end_indices"]
+
+    for end_index, updated_end_index in zip(end_indices, updated_end_indices):
+        nt.assert_true(end_index != updated_end_index)
+
+
 def test_mean_frequency1():
     """basic: Test mean_frequency 1"""