Started python-only features, added initburst_sahp as first example

Makefile

@@ -36,6 +36,9 @@ update_version:
 test: virtualenv install install_test_requirements
 	cd efel/tests; nosetests -s -v -x --with-coverage --cover-xml \
 	   --cover-package efel 
+debugtest: virtualenv install install_test_requirements
+	cd efel/tests; nosetests -a debugtest -s -v -x --with-coverage --cover-xml \
+	   --cover-package efel 
 pypi: test
 	pip install twine --upgrade
 	rm -rf dist

efel/api.py

@@ -33,6 +33,8 @@
 import efel
 import efel.cppcore as cppcore
 
+import efel.pyfeatures as pyfeatures
+
 """
 Disabling cppcore importerror override, it confuses users in case the error
 is caused by something else
@@ -78,6 +80,9 @@ def reset():
     setDoubleSetting("initial_perc", 0.1)
     setDoubleSetting("min_spike_height", 20.0)
     setIntSetting("strict_stiminterval", 0)
+    setDoubleSetting("initburst_freq_threshold", 50)
+    setDoubleSetting("initburst_sahp_start", 5)
+    setDoubleSetting("initburst_sahp_end", 100)
 
     _initialise()
 
@@ -166,6 +171,8 @@ def getFeatureNames():
     feature_names = []
     cppcore.getFeatureNames(feature_names)
 
+    feature_names += pyfeatures.all_pyfeatures
+
     return feature_names
 
 
@@ -331,6 +338,12 @@ def getFeatureValues(
         return map_result
 
 
+def get_py_feature(featureName):
+    """Return python feature"""
+
+    return getattr(pyfeatures, featureName)()
+
+
 def _get_feature_values_serial(trace_featurenames):
     """Single thread of getFeatureValues"""
 
@@ -342,7 +355,7 @@ def _get_feature_values_serial(trace_featurenames):
         try:
             len(trace['stim_start'])
             len(trace['stim_end'])
-        except:
+        except BaseException:
             raise Exception('Unable to determine length of stim_start or '
                             'stim_end, are you sure these are lists ?')
 
@@ -367,19 +380,22 @@ def _get_feature_values_serial(trace_featurenames):
         cppcore.setFeatureDouble(item, [x for x in trace[item]])
 
     for featureName in featureNames:
-        cppcoreFeatureValues = list()
-        exitCode = cppcore.getFeature(featureName, cppcoreFeatureValues)
-
-        if exitCode < 0:
-            if raise_warnings:
-                import warnings
-                warnings.warn(
-                    "Error while calculating feature %s: %s" %
-                    (featureName, cppcore.getgError()),
-                    RuntimeWarning)
-            featureDict[featureName] = None
+        if featureName in pyfeatures.all_pyfeatures:
+            featureDict[featureName] = get_py_feature(featureName)
         else:
-            featureDict[featureName] = numpy.array(cppcoreFeatureValues)
+            cppcoreFeatureValues = list()
+            exitCode = cppcore.getFeature(featureName, cppcoreFeatureValues)
+
+            if exitCode < 0:
+                if raise_warnings:
+                    import warnings
+                    warnings.warn(
+                        "Error while calculating feature %s: %s" %
+                        (featureName, cppcore.getgError()),
+                        RuntimeWarning)
+                featureDict[featureName] = None
+            else:
+                featureDict[featureName] = numpy.array(cppcoreFeatureValues)
 
     return featureDict
 
@@ -427,4 +443,5 @@ def getMeanFeatureValues(traces, featureNames, raise_warnings=True):
 
     return featureDicts
 
+
 reset()

efel/cppcore/cppcore.cpp

@@ -145,6 +145,32 @@ _getfeature(PyObject* self, PyObject* args, const string &type) {
   return Py_BuildValue("i", return_value);
 }
 
+
+static PyObject*
+_getmapdata(PyObject* self, PyObject* args, const string &type) {
+  char* data_name;
+  PyObject* py_values = PyList_New(0);
+
+  if (!PyArg_ParseTuple(args, "s", &data_name)) {
+    return NULL;
+  }
+
+  if (type == "int") {
+    vector<int> return_values;
+    return_values = pFeature->getmapIntData(string(data_name));
+    PyList_from_vectorint(return_values, py_values);
+  } else if (type == "double") {
+    vector<double> return_values;
+    return_values = pFeature->getmapDoubleData(string(data_name));
+    PyList_from_vectordouble(return_values, py_values);
+  } else {
+    PyErr_SetString(PyExc_TypeError, "Unknown data name");
+    return NULL;
+  }
+
+  return py_values;
+}
+
 static PyObject* getfeature(PyObject* self, PyObject* args) {
   const string empty("");
   return _getfeature(self, args, empty);
@@ -169,6 +195,14 @@ static PyObject* getfeatureint(PyObject* self, PyObject* args) {
   const string type("int");
   return _getfeature(self, args, type);
 }
+static PyObject* getmapintdata(PyObject* self, PyObject* args) {
+  const string type("int");
+  return _getmapdata(self, args, type);
+}
+static PyObject* getmapdoubledata(PyObject* self, PyObject* args) {
+  const string type("double");
+  return _getmapdata(self, args, type);
+}
 
 static PyObject* setfeaturedouble(PyObject* self, PyObject* args) {
   char* feature_name;
@@ -190,6 +224,7 @@ static PyObject* getfeaturedouble(PyObject* self, PyObject* args) {
   return _getfeature(self, args, type);
 }
 
+
 static PyObject* getFeatureNames(PyObject* self, PyObject* args) {
   vector<string> feature_names;
   PyObject* py_feature_names;
@@ -253,6 +288,10 @@ static PyMethodDef CppCoreMethods[] = {
       "Get a integer feature."},
     {"getFeatureDouble", getfeaturedouble, METH_VARARGS,
       "Get a double feature."},
+    {"getMapIntData", getmapintdata, METH_VARARGS,
+      "Get a int data."},
+    {"getMapDoubleData", getmapdoubledata, METH_VARARGS,
+      "Get a double data."},
 
     {"setFeatureInt", setfeatureint, METH_VARARGS,
       "Set a integer feature."},

efel/pyfeatures/__init__.py

@@ -0,0 +1,31 @@
+"""Python implementation of features"""
+
+"""
+Copyright (c) 2015, Blue Brain Project/EPFL
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+  * Redistributions of source code must retain the above copyright
+    notice, this list of conditions and the following disclaimer.
+  * Redistributions in binary form must reproduce the above copyright
+    notice, this list of conditions and the following disclaimer in the
+    documentation and/or other materials provided with the distribution.
+  * Neither the name of the copyright holder nor the
+    names of its contributors may be used to endorse or promote products
+    derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+"""
+
+from pyfeatures import *  # NOQA

efel/pyfeatures/pyfeatures.py

@@ -0,0 +1,183 @@
+"""Python implementation of features"""
+
+"""
+Copyright (c) 2015, Blue Brain Project/EPFL
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+  * Redistributions of source code must retain the above copyright
+    notice, this list of conditions and the following disclaimer.
+  * Redistributions in binary form must reproduce the above copyright
+    notice, this list of conditions and the following disclaimer in the
+    documentation and/or other materials provided with the distribution.
+  * Neither the name of the copyright holder nor the
+    names of its contributors may be used to endorse or promote products
+    derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+"""
+
+
+import numpy
+import efel.cppcore
+
+
+all_pyfeatures = [
+    'voltage',
+    'time',
+    'ISIs',
+    'initburst_sahp',
+    'initburst_sahp_vb',
+    'initburst_sahp_ssse']
+
+
+def voltage():
+    """Get voltage trace"""
+    return _get_cpp_feature("voltage")
+
+
+def time():
+    """Get time trace"""
+    return _get_cpp_feature("time")
+
+
+def ISIs():
+    """Get all ISIs"""
+
+    peak_times = _get_cpp_feature("peak_time")
+
+    return numpy.diff(peak_times)
+
+
+def initburst_sahp_vb():
+    """SlowAHP voltage from voltage base after initial burst"""
+
+    # Required cpp features
+    initburst_sahp_value = initburst_sahp()
+    voltage_base = _get_cpp_feature("voltage_base")
+
+    if initburst_sahp_value is None or voltage_base is None or \
+            len(initburst_sahp_value) != 1 or len(voltage_base) != 1:
+        return None
+    else:
+        return numpy.array([initburst_sahp_value[0] - voltage_base[0]])
+
+
+def initburst_sahp_ssse():
+    """SlowAHP voltage from steady_state_voltage_stimend after initial burst"""
+
+    # Required cpp features
+    initburst_sahp_value = initburst_sahp()
+    ssse = _get_cpp_feature("steady_state_voltage_stimend")
+
+    if initburst_sahp_value is None or ssse is None or \
+            len(initburst_sahp_value) != 1 or len(ssse) != 1:
+        return None
+    else:
+        return numpy.array([initburst_sahp_value[0] - ssse[0]])
+
+
+def initburst_sahp():
+    """SlowAHP voltage after initial burst"""
+
+    # Required cpp features
+    voltage = _get_cpp_feature("voltage")
+    time = _get_cpp_feature("time")
+    peak_times = _get_cpp_feature("peak_time")
+
+    # Required python features
+    all_isis = ISIs()
+
+    # Required trace data
+    stim_end = _get_cpp_data("stim_end")
+
+    # Required settings
+    initburst_freq_thresh = _get_cpp_data("initburst_freq_threshold")
+    initburst_sahp_start = _get_cpp_data("initburst_sahp_start")
+    initburst_sahp_end = _get_cpp_data("initburst_sahp_end")
+
+    last_isi = None
+
+    # Loop over ISIs until frequency higher than initburst_freq_threshold
+    for isi_counter, isi in enumerate(all_isis):
+        # Convert to Hz
+        freq = 1000.0 / isi
+        if freq < initburst_freq_thresh:
+            # Threshold reached
+            break
+        else:
+            # Add isi to initburst
+            last_isi = isi_counter
+
+    if last_isi is None:
+        # No initburst found
+        return None
+    else:
+        # Get index of second peak of last ISI
+        last_peak = last_isi + 1
+
+    # Get time of last peak
+    last_peak_time = peak_times[last_peak]
+
+    # Determine start of sahp interval
+    sahp_interval_start = min(
+        last_peak_time +
+        initburst_sahp_start,
+        stim_end)
+
+    # Get next peak, we wont search beyond that
+    next_peak = last_peak + 1
+
+    # Determine end of sahp interval
+    # Add initburst_slow_ahp_max to last peak time
+    # If next peak or stim_end is earlier, use these
+    # If no next peak, use stim end
+    if next_peak < len(peak_times):
+        next_peak_time = peak_times[next_peak]
+
+        sahp_interval_end = min(
+            last_peak_time + initburst_sahp_end, next_peak_time, stim_end)
+    else:
+        sahp_interval_end = min(
+            last_peak_time + initburst_sahp_end, stim_end)
+
+    if sahp_interval_end <= sahp_interval_start:
+        return None
+    else:
+        sahp_interval = voltage[numpy.where(
+            (time <= sahp_interval_end) &
+            (time >= sahp_interval_start))]
+
+        min_volt_index = numpy.argmin(sahp_interval)
+
+        slow_ahp = sahp_interval[min_volt_index]
+
+        return numpy.array([slow_ahp])
+
+
+def _get_cpp_feature(feature_name):
+    """Get cpp feature"""
+    cppcoreFeatureValues = list()
+    exitCode = efel.cppcore.getFeature(feature_name, cppcoreFeatureValues)
+
+    if exitCode < 0:
+        return None
+    else:
+        return numpy.array(cppcoreFeatureValues)
+
+
+def _get_cpp_data(data_name):
+    """Get cpp data value"""
+
+    return efel.cppcore.getMapDoubleData(data_name)[0]

efel/tests/featurenames.json

@@ -152,5 +152,11 @@
     "voltage_base",
     "voltage_deflection",
     "voltage_deflection_begin",
-    "voltage_deflection_vb_ssse"
+    "voltage_deflection_vb_ssse",
+    "voltage",
+    "time",
+    "ISIs",
+    "initburst_sahp",
+    "initburst_sahp_vb",
+    "initburst_sahp_ssse"
 ]