Merge pull request #165 from mariarv/master

Impedance ("Z score") feature is added

docs/source/eFeatures.rst

@@ -2000,6 +2000,36 @@ The burst detection can be fine-tuned by changing the setting strict_burst_facto
 
     numpy.array([spikes_per_burst[0] - spikes_per_burst[-1]])
 
+`Python efeature`_ : impedance
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Computes the impedance given a ZAP current input and its voltage response.
+It will return the frequency at which the impedance is maximal, in the range (0, impedance_max_freq] Hz,
+with impedance_max_freq being a setting with 50.0 as a default value.
+
+- **Required features**: current, LibV1:Spikecount, LibV5:voltage_base, LibV5:current_base
+- **Units**: Hz
+- **Pseudocode**: ::
+
+    normalized_voltage = voltage_trace - voltage_base
+    normalized_current = current_trace - current_base
+    if spike_count < 1:  # if there is no spikes in ZAP
+        fft_volt = numpy.fft.fft(normalized_voltage)
+        fft_cur = numpy.fft.fft(normalized_current)
+        if any(fft_cur) == 0:
+            return None
+        # convert dt from ms to s to have freq in Hz
+        freq = numpy.fft.fftfreq(len(normalized_voltage), d=dt / 1000.)
+        Z = fft_volt / fft_cur
+        norm_Z = abs(Z) / max(abs(Z))
+        select_idxs = numpy.swapaxes(numpy.argwhere((freq > 0) & (freq <= impedance_max_freq)), 0, 1)[0]
+        smooth_Z = gaussian_filter1d(norm_Z[select_idxs], 10)
+        ind_max = numpy.argmax(smooth_Z)
+        return freq[ind_max]
+    else:
+        return None
+
+
 
 .. _LibV1: https://github.com/BlueBrain/eFEL/blob/master/efel/cppcore/LibV1.cpp
 .. _LibV2: https://github.com/BlueBrain/eFEL/blob/master/efel/cppcore/LibV2.cpp

efel/api.py

@@ -98,6 +98,7 @@ def reset():
     setDoubleSetting("precision_threshold", 1e-10)
     setDoubleSetting("sahp_start", 5.0)
     setIntSetting("ignore_first_ISI", 1)
+    setDoubleSetting("impedance_max_freq", 50.0)
 
     _initialise()
 

efel/pyfeatures/pyfeatures.py

@@ -28,6 +28,7 @@
 
 import numpy
 import efel.cppcore
+from numpy.fft import *
 
 
 all_pyfeatures = [
@@ -43,7 +44,8 @@
     'spikes_per_burst',
     'spikes_per_burst_diff',
     'spikes_in_burst1_burst2_diff',
-    'spikes_in_burst1_burstlast_diff'
+    'spikes_in_burst1_burstlast_diff',
+    'impedance',
 ]
 
 
@@ -57,6 +59,40 @@ def time():
     return _get_cpp_feature("time")
 
 
+def impedance():
+    from scipy.ndimage.filters import gaussian_filter1d
+
+    dt = _get_cpp_data("interp_step")
+    Z_max_freq = _get_cpp_data("impedance_max_freq")
+    voltage_trace = voltage()
+    holding_voltage = _get_cpp_feature("voltage_base")
+    normalized_voltage = voltage_trace - holding_voltage
+    current_trace = current()
+    if current_trace is not None:
+        holding_current = _get_cpp_feature("current_base")
+        normalized_current = current_trace - holding_current
+        spike_count = _get_cpp_feature("Spikecount")
+        if spike_count < 1:  # if there is no spikes in ZAP
+            fft_volt = numpy.fft.fft(normalized_voltage)
+            fft_cur = numpy.fft.fft(normalized_current)
+            if any(fft_cur) == 0:
+                return None
+            # convert dt from ms to s to have freq in Hz
+            freq = numpy.fft.fftfreq(len(normalized_voltage), d=dt / 1000.)
+            Z = fft_volt / fft_cur
+            norm_Z = abs(Z) / max(abs(Z))
+            select_idxs = numpy.swapaxes(
+                numpy.argwhere((freq > 0) & (freq <= Z_max_freq)), 0, 1
+            )[0]
+            smooth_Z = gaussian_filter1d(norm_Z[select_idxs], 10)
+            ind_max = numpy.argmax(smooth_Z)
+            return freq[ind_max]
+        else:
+            return None
+    else:
+        return None
+
+
 def current():
     """Get current trace"""
     return _get_cpp_feature("current")
@@ -66,7 +102,6 @@ def ISIs():
     """Get all ISIs"""
 
     peak_times = _get_cpp_feature("peak_time")
-
     return numpy.diff(peak_times)
 
 

efel/units/units.json

@@ -143,5 +143,6 @@
     "initburst_sahp": "mV",
     "initburst_sahp_ssse": "mV",
     "initburst_sahp_vb": "mV",
-    "depol_block_bool": "constant"
+    "depol_block_bool": "constant",
+    "impedance": "Hz"
 }

tests/featurenames.json

@@ -115,6 +115,7 @@
     "doublet_ISI",
     "fast_AHP",
     "fast_AHP_change",
+    "impedance",
     "initburst_sahp",
     "initburst_sahp_ssse",
     "initburst_sahp_vb",

tests/test_allfeatures.py

@@ -121,7 +121,7 @@ def get_allfeature_values():
         soma_featurenames.remove(feature_name)
 
     soma_featurenames = [x for x in soma_featurenames
-                         if x not in ['current', 'current_base']]
+                         if x not in ['current', 'current_base', 'impedance']]
 
     import warnings
     with warnings.catch_warnings():

tests/test_pyfeatures.py

@@ -131,6 +131,16 @@
         'v_col': 3,
         'stim_start': 700.0,
         'stim_end': 2700.0},
+    'impedance': {
+        'url': 'file://%s' % os.path.join(
+            os.path.abspath(testdata_dir),
+            'basic',
+            'impedance.txt'),
+        't_col': 1,
+        'v_col': 2,
+        'i_col': 3,
+        'stim_start': 100.0,
+        'stim_end': 51000.0}
 }
 
 
@@ -366,3 +376,12 @@ def interpolate(time, voltage, new_dt):
     assert len(interp_voltage) == len(feature_voltage)
     assert len(feature_voltage) == len(feature_current)
     assert numpy.allclose(interp_current, feature_current, atol=1e-6)
+
+
+def test_impedance():
+    """pyfeatures: Test impedance feature"""
+
+    feature_name = "impedance"
+
+    expected_values = {feature_name: 4.615384615384615}
+    _test_expected_value(feature_name, expected_values)