diff --git a/pcpostprocess/infer_reversal.py b/pcpostprocess/infer_reversal.py index 102c51cc..fd5b892b 100644 --- a/pcpostprocess/infer_reversal.py +++ b/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,21 +64,22 @@ 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}$") @@ -72,11 +87,7 @@ def infer_reversal_potential(current, times, voltage_segments, voltages, 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 diff --git a/pcpostprocess/scripts/run_herg_qc.py b/pcpostprocess/scripts/run_herg_qc.py index 4020d4fb..5ef4a52c 100644 --- a/pcpostprocess/scripts/run_herg_qc.py +++ b/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))) diff --git a/tests/test_herg_qc.py b/tests/test_herg_qc.py old mode 100644 new mode 100755 index 9a060179..ebf64ef0 --- a/tests/test_herg_qc.py +++ b/tests/test_herg_qc.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python3 import copy import os import string diff --git a/tests/test_leak_correct.py b/tests/test_leak_correct.py old mode 100644 new mode 100755 index 2c984219..7b213a83 --- a/tests/test_leak_correct.py +++ b/tests/test_leak_correct.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python3 import os import unittest diff --git a/tests/test_subtraction_plots.py b/tests/test_subtraction_plots.py old mode 100644 new mode 100755 index 06521129..85bf47b6 --- a/tests/test_subtraction_plots.py +++ b/tests/test_subtraction_plots.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python3 import os import unittest