Fixed AP_amplitude so that it works when there are spike outside of the

stimulus intervals (which are ignored)
BlueBrain · Oct 5, 2015 · 159cdea · 159cdea
1 parent 55fedb8
commit 159cdea
Show file tree

Hide file tree

Showing 8 changed files with 8,815 additions and 9 deletions.
diff --git a/efel/DependencyV5.txt b/efel/DependencyV5.txt
@@ -15,7 +15,7 @@ LibV1:burst_mean_freq	    #LibV1:burst_ISI_indices    #LibV1:peak_time
 LibV1:burst_number	    #LibV1:burst_mean_freq
 LibV1:interburst_voltage	#LibV1:burst_ISI_indices
 LibV1:AP_height             #LibV1:peak_voltage
-LibV1:AP_amplitude          #LibV5:AP_begin_indices #LibV1:peak_voltage
+LibV1:AP_amplitude          #LibV5:AP_begin_indices #LibV1:peak_voltage #LibV1:peak_time
 LibV1:AHP_depth_abs_slow	#LibV1:peak_indices
 LibV1:AHP_slow_time         #LibV1:AHP_depth_abs_slow
 LibV1:time_constant

diff --git a/efel/cppcore/DependencyTree.cpp b/efel/cppcore/DependencyTree.cpp
@@ -120,7 +120,7 @@ int cTree::setFeaturePointers(
       // Find the feature function pointer map for the library
       mapLibItr = mapFptrLib.find(strLib);
       if (mapLibItr == mapFptrLib.end()) {
-        ErrorStr = ErrorStr + string("\nLibrary [") + strLib + "] is missing";
+        ErrorStr = ErrorStr + string("\nLibrary [") + strLib + "] is missing\n";
         return (-1);
       }
 

diff --git a/efel/cppcore/LibV1.cpp b/efel/cppcore/LibV1.cpp
@@ -456,36 +456,73 @@ int LibV1::AP_amplitude(mapStr2intVec& IntFeatureData,
     return nSize;
   } else {
     vector<double> peakvoltage;
+    vector<double> peaktime;
     vector<int> apbeginindices;
     vector<double> v;
     retVal = getDoubleVec(DoubleFeatureData, StringData, "V", v);
     if (retVal <= 0) {
       GErrorStr += "AP_amplitude: Can't find voltage vector V";  
       return -1;
     }
+
+    vector<double> stimstart;                                                      
+    retVal = getDoubleVec(DoubleFeatureData, StringData, "stim_start", stimstart); 
+    if (retVal != 1) {
+        GErrorStr += "AP_amplitude: Error getting stim_start";  
+        return -1;                                                     
+    }
+
+    vector<double> stimend;                                                        
+    retVal = getDoubleVec(DoubleFeatureData, StringData, "stim_end", stimend);
+    if (retVal != 1) {
+        GErrorStr += "AP_amplitude: Error getting stim_end";  
+        return -1;                                                     
+    }
+
     retVal = getDoubleVec(DoubleFeatureData, StringData, "peak_voltage",
                           peakvoltage);
     if (retVal <= 0) {
       GErrorStr += "AP_amplitude: Error calculating peak_voltage";
       return -1;
     }
+
+    retVal = getDoubleVec(DoubleFeatureData, StringData, "peak_time",
+                          peaktime);
+    if (retVal <= 0) {
+      GErrorStr += "AP_amplitude: Error calculating peak_time";
+      return -1;
+    }
+
     retVal = getIntVec(IntFeatureData, StringData, "AP_begin_indices",
                        apbeginindices);
     if (retVal <= 0) {
       GErrorStr += "AP_amplitude: Error calculating AP_begin_indicies";
       return -1;
     }
 
-    if (peakvoltage.size() > apbeginindices.size()) {
+    if (peakvoltage.size() != peaktime.size()) {
+      GErrorStr += 
+          "AP_amplitude: Not the same amount of peak_time and peak_voltage entries";
+      return -1;
+    }
+
+    vector<double> peakvoltage_duringstim;
+    for (unsigned i = 0; i < peaktime.size(); i++) {
+      if (peaktime[i] >= stimstart[0] && peaktime[i] <= stimend[0]) {
+         peakvoltage_duringstim.push_back(peakvoltage[i]); 
+      }
+    }
+
+    if (peakvoltage_duringstim.size() > apbeginindices.size()) {
       GErrorStr += 
-          "AP_amplitude: More peak_voltage entries than AP_begin_indices entries";
+          "AP_amplitude: More peak_voltage entries during the stimulus than AP_begin_indices entries";
       return -1;
     }
 
     vector<double> apamplitude;
-    apamplitude.resize(peakvoltage.size());
+    apamplitude.resize(peakvoltage_duringstim.size());
     for (unsigned i = 0; i < apamplitude.size(); i++) {
-      apamplitude[i] = peakvoltage[i] - v[apbeginindices[i]];
+      apamplitude[i] = peakvoltage_duringstim[i] - v[apbeginindices[i]];
     }
     setDoubleVec(DoubleFeatureData, StringData, "AP_amplitude", apamplitude);
     return apamplitude.size();

diff --git a/efel/cppcore/mapoperations.cpp b/efel/cppcore/mapoperations.cpp
@@ -89,7 +89,7 @@ int getIntVec(mapStr2intVec& IntFeatureData, mapStr2Str& StringData,
   strFeature += params;
   mapStr2intVec::iterator mapstr2IntItr(IntFeatureData.find(strFeature));
   if (mapstr2IntItr == IntFeatureData.end()) {
-    GErrorStr = GErrorStr + "\nFeature [" + strFeature + "] is missing";
+    GErrorStr = GErrorStr + "\nFeature [" + strFeature + "] is missing\n";
     return -1;
   }
   v = mapstr2IntItr->second;
@@ -104,7 +104,7 @@ int getDoubleVec(mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData,
   mapStr2doubleVec::iterator mapstr2DoubleItr(
       DoubleFeatureData.find(strFeature));
   if (mapstr2DoubleItr == DoubleFeatureData.end()) {
-    GErrorStr = GErrorStr + "\nFeature [" + strFeature + "] is missing";
+    GErrorStr = GErrorStr + "\nFeature [" + strFeature + "] is missing\n";
     return -1;
   }
   v = mapstr2DoubleItr->second;

diff --git a/efel/tests/.gitignore b/efel/tests/.gitignore
@@ -1,3 +1,2 @@
 /coverage.xml
 /.coverage
-/log
diff --git a/efel/tests/log/.gitignore b/efel/tests/log/.gitignore
@@ -0,0 +1 @@
+/fllog.txt
diff --git a/efel/tests/test_basic.py b/efel/tests/test_basic.py
@@ -405,6 +405,44 @@ def test_mean_frequency1():
     nt.assert_almost_equal(feature_values[0]['mean_frequency'][0], 15.2858453)
 
 
+def test_ap_amplitude_outside_stim():
+    """basic: Test AP amplitude with spike outside stim"""
+
+    import efel
+    efel.reset()
+    import numpy
+
+    stim_start = 700.0
+    stim_end = 2700.0
+
+    test_data_path = joinp(testdata_dir, 'basic', 'spike_outside_stim.txt')
+    data = numpy.loadtxt(test_data_path)
+
+    time = data[:, 0]
+    voltage = data[:, 1]
+
+    trace = {}
+
+    trace['T'] = time
+    trace['V'] = voltage
+    trace['stim_start'] = [stim_start]
+    trace['stim_end'] = [stim_end]
+
+    features = ['AP_amplitude', 'peak_time']
+
+    feature_values = \
+        efel.getFeatureValues(
+            [trace],
+            features)
+
+    # Make sure AP_amplitude doesn't pick up the two spikes outside of
+    # the stimulus
+    # (which are present in peak_time)
+    nt.assert_equal(
+        len(feature_values[0]['AP_amplitude']) + 2,
+        len(feature_values[0]['peak_time']))
+
+
 def test_ap_amplitude_from_voltagebase1():
     """basic: Test AP_amplitude_from_voltagebase 1"""