Merge pull request #189 from AurelienJaquier/impedance-plot

Impedance plot

Author: AurelienJaquier

bluepyefe/extract.py

@@ -1042,7 +1042,11 @@ def extract_efeatures(
 
     if plot:
         plot_all_recordings_efeatures(
-            cells, protocols, output_dir=output_directory, mapper=map_function
+            cells,
+            protocols,
+            output_dir=output_directory,
+            mapper=map_function,
+            efel_settings=efel_settings
         )
 
     if extract_per_cell and write_files:

bluepyefe/plotting.py

@@ -342,15 +342,91 @@ def plot_grouped_efeatures(
     _ = _plot_legend(colors, markers, output_dir, show=show)
 
 
+def plot_impedance(cell, output_dir, efel_settings):
+    """Plots the impedance."""
+    from scipy.ndimage.filters import gaussian_filter1d
+
+    dt = 0.1
+    Z_max_freq = 50.0
+    if efel_settings is not None:
+        dt = efel_settings.get("interp_step", dt)
+        Z_max_freq = efel_settings.get("impedance_max_freq", Z_max_freq)
+
+    for protocol_name in cell.get_protocol_names():
+        if "sinespec" in protocol_name.lower():
+            recordings = cell.get_recordings_by_protocol_name(protocol_name)
+            for i, rec in enumerate(recordings):
+                voltage = rec.voltage
+                current = rec.current
+
+                efel_vals = rec.call_efel(
+                    [
+                        "voltage_base",
+                        "steady_state_voltage_stimend",
+                        "current_base",
+                        "steady_state_current_stimend",
+                    ],
+                    efel_settings
+                )
+                if efel_vals[0]["voltage_base"] is not None:
+                    holding_voltage = efel_vals[0]["voltage_base"][0]
+                else:
+                    holding_voltage = efel_vals[0]["steady_state_voltage_stimend"][0]
+                if efel_vals[0]["current_base"] is not None:
+                    holding_current = efel_vals[0]["current_base"][0]
+                else:
+                    holding_current = efel_vals[0]["steady_state_current_stimend"][0]
+
+                normalized_voltage = voltage - holding_voltage
+                normalized_current = current - holding_current
+
+                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)
+
+                filename = "{}_{}_{}_{}_impedance.pdf".format(
+                    cell.name, protocol_name, i, rec.name
+                )
+                dirname = pathlib.Path(output_dir) / cell.name / "impedance"
+
+                fig = plt.figure()
+                ax = fig.add_subplot(1, 1, 1)
+                ax.plot(freq[:len(smooth_Z)], smooth_Z)
+                ax.set_xlabel("Frequency (Hz)")
+                ax.set_ylabel("normalized Z")
+                fig .suptitle(f"Impedance for {rec.name}\nfor cell {cell.name}")
+                _save_fig(dirname, filename)
+
+
+def plot_all_impedances(cells, output_dir, mapper=map, efel_settings=None):
+    """Plot recordings for all cells and all protocols"""
+    if mapper == map:
+        # For the built-in map(), ensure immediate evaluation as it returns a lazy iterator
+        # which won't execute the function until iterated over. Converting to a list forces this iteration.
+        list(mapper(partial(plot_impedance, output_dir=output_dir, efel_settings=efel_settings), cells))
+    else:
+        mapper(partial(plot_impedance, output_dir=output_dir, efel_settings=efel_settings), cells)
+
+
 def plot_all_recordings_efeatures(
-    cells, protocols, output_dir=None, show=False, mapper=map
+    cells, protocols, output_dir=None, show=False, mapper=map, efel_settings=None
 ):
     """Generate plots for all recordings and efeatures both for individual
     cells and across all cells."""
 
     colors, markers = _get_colors_markers_wheels(cells)
 
     plot_all_recordings(cells, output_dir, mapper=mapper)
+    plot_all_impedances(cells, output_dir, mapper=map, efel_settings=efel_settings)
 
     for key_amp in ["amp", "amp_rel"]:
         plot_individual_efeatures(

bluepyefe/recording.py

@@ -25,6 +25,7 @@
 import numpy
 import efel
 import matplotlib.pyplot as plt
+from pathlib import Path
 
 from .tools import to_ms, to_mV, to_nA, set_efel_settings
 
@@ -81,6 +82,17 @@ def __init__(self, config_data, reader_data, protocol_name):
         self.auto_threshold = None
         self.peak_time = None
 
+    @property
+    def name(self):
+        """Proxy that can be used to name the recording."""
+        if self.id:
+            return self.id
+        if "filepath" in self.config_data and self.config_data["filepath"] is not None:
+            return str(Path(self.config_data["filepath"]).stem)
+        if "v_file" in self.config_data:
+            return str(Path(self.config_data["v_file"]).stem)
+        return ""
+
     @property
     def time(self):
         """Alias of the time attribute"""