Merge pull request #78 from CardiacModelling/infer-e-docstring

Added provisional docstring to infer_reversal_potential

Author: MichaelClerx

pcpostprocess/infer_reversal.py

@@ -5,41 +5,55 @@
 import numpy as np
 import numpy.polynomial.polynomial as poly
 
+from .detect_ramp_bounds import detect_ramp_bounds
 
-def infer_reversal_potential(current, times, voltage_segments, voltages,
-                             ax=None, output_path=None, plot=None,
-                             known_Erev=None, figsize=(5, 3)):
-
-    if output_path:
-        dirname = os.path.dirname(output_path)
-        if not os.path.exists(dirname):
-            os.makedirs(dirname)
-
-    if (ax or output_path) and plot is not False:
-        plot = True
-
-    # Find indices of observations during the reversal ramp
-    ramps = [line for line in voltage_segments if line[2] != line[3]]
 
-    # Assume the last ramp is the reversal ramp (convert to ms)
-    tstart, tend = np.array(ramps)[-1, :2]
+def infer_reversal_potential(current, times, voltage_segments, voltages,
+                             output_path=None, known_Erev=None,
+                             figsize=(5, 3)):
+    """
+    Infers a reversal potential in a time series, based on a reversal ramp.
+
+    The data is denoised by fitting a 4-th order polynomial through the ramp
+    data, from which a reversal potential is then detected. If no polynomial
+    can be fit or the resulting zero-crossing is outside of
+    ``min(voltages), max(voltages)``, then ``np.nan`` is returned.
+
+    @param current: The currents that make up a time series with ``times``
+    @param times: The sampled times
+    @param voltage_segments: A list of tuples (tstart, tend, vstart, vend)
+    describing voltage steps or ramps. It is assumed the final ramp is the
+    reversal ramp.
+    @param voltages: The sampled voltages
+    @param output_path: An optional path to store a plot at
+    @param known_Erev: A known reversal potential to include in the plot
+    @param figsize: A size for the plot.
+
+    @return: The inferred reversal potential
+    """
+
+    # Get ramp bounds, assuming final ramp is the reversal ramp
+    tstart, tend = detect_ramp_bounds(times, voltage_segments, -1)
 
     istart = np.argmax(times > tstart)
     iend = np.argmax(times > tend)
 
-    times = times[istart:iend]
     current = current[istart:iend]
     voltages = voltages[istart:iend]
 
+    # Fit a 4-th order polynomial
     try:
         fitted_poly = poly.Polynomial.fit(voltages, current, 4)
     except ValueError as exc:
         logging.warning(str(exc))
         return np.nan
 
+    # Try extracting the polynomial's roots, accepting only ones that are
+    # within the range of sampled voltages (so not using ramp info here!)
     try:
+        vmin, vmax = np.min(voltages), np.max(voltages)
         roots = np.unique([np.real(root) for root in fitted_poly.roots()
-                           if root > np.min(voltages) and root < np.max(voltages)])
+                           if root > vmin and root < vmax])
     except np.linalg.LinAlgError as exc:
         logging.warning(str(exc))
         return np.nan
@@ -50,33 +64,30 @@ def infer_reversal_potential(current, times, voltage_segments, voltages,
 
     if len(roots) == 0:
         return np.nan
+    erev = roots[-1]
 
-    if plot:
-        created_fig = False
-        if ax is None and output_path is not None:
-
-            created_fig = True
-            fig = plt.figure(figsize=figsize)
-            ax = fig.subplots()
+    # Optional plot
+    if output_path is not None:
+        dirname = os.path.dirname(output_path)
+        if not os.path.exists(dirname):
+            os.makedirs(dirname)
 
-        ax.set_xlabel('$V$ (mV)')
+        fig = plt.figure(figsize=figsize)
+        ax = fig.subplots()
+        ax.set_xlabel('$V$ (mV)')  # Assuming mV here
         ax.set_ylabel('$I$ (nA)')
 
         # Now plot current vs voltage
         ax.plot(voltages, current, 'x', markersize=2, color='grey', alpha=.5)
-        ax.axvline(roots[-1], linestyle='--', color='grey', label=r'$E_\mathrm{obs}$')
+        ax.axvline(erev, linestyle='--', color='grey', label=r'$E_\mathrm{obs}$')
         if known_Erev:
             ax.axvline(known_Erev, linestyle='--', color='orange',
                        label="Calculated $E_{Kr}$")
         ax.axhline(0, linestyle='--', color='grey')
         ax.plot(*fitted_poly.linspace())
         ax.legend()
 
-        if output_path is not None:
-            fig = ax.figure
-            fig.savefig(output_path)
-
-        if created_fig:
-            plt.close(fig)
+        fig.savefig(output_path)
+        plt.close(fig)
 
-    return roots[-1]
+    return erev

pcpostprocess/scripts/run_herg_qc.py

@@ -660,24 +660,20 @@ def extract_protocol(readname, savename, time_strs, selected_wells, args):
             after_corrected = after_current[sweep, :] - after_leak
             before_corrected = before_current[sweep, :] - before_leak
 
-            E_rev_before = infer_reversal_potential(before_corrected, times,
-                                                    desc, voltages, plot=True,
-                                                    output_path=os.path.join(reversal_plot_dir,
-                                                                             f"{well}_{savename}_sweep{sweep}_before"),
-                                                    known_Erev=args.Erev)
-
-            E_rev_after = infer_reversal_potential(after_corrected, times,
-                                                   desc, voltages,
-                                                   plot=True,
-                                                   output_path=os.path.join(reversal_plot_dir,
-                                                                            f"{well}_{savename}_sweep{sweep}_after"),
-                                                   known_Erev=args.Erev)
-
-            E_rev = infer_reversal_potential(subtracted_trace, times, desc,
-                                             voltages, plot=True,
-                                             output_path=os.path.join(reversal_plot_dir,
-                                                                      f"{well}_{savename}_sweep{sweep}_subtracted"),
-                                             known_Erev=args.Erev)
+            E_rev_before = infer_reversal_potential(
+                before_corrected, times, desc, voltages,
+                output_path=os.path.join(reversal_plot_dir, f"{well}_{savename}_sweep{sweep}_before"),
+                known_Erev=args.Erev)
+
+            E_rev_after = infer_reversal_potential(
+                after_corrected, times, desc, voltages,
+                output_path=os.path.join(reversal_plot_dir, f"{well}_{savename}_sweep{sweep}_after"),
+                known_Erev=args.Erev)
+
+            E_rev = infer_reversal_potential(
+                subtracted_trace, times, desc, voltages,
+                output_path=os.path.join(reversal_plot_dir, f"{well}_{savename}_sweep{sweep}_subtracted"),
+                known_Erev=args.Erev)
 
             row_dict['R_leftover'] =\
                 np.sqrt(np.sum((after_corrected)**2)/(np.sum(before_corrected**2)))

tests/test_herg_qc.py

@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 import copy
 import os
 import string

tests/test_leak_correct.py

@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 import os
 import unittest
 

tests/test_subtraction_plots.py

@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 import os
 import unittest