55import numpy as np
66import numpy .polynomial .polynomial as poly
77
8+ from .detect_ramp_bounds import detect_ramp_bounds
89
9- def infer_reversal_potential (current , times , voltage_segments , voltages ,
10- ax = None , output_path = None , plot = None ,
11- known_Erev = None , figsize = (5 , 3 )):
12-
13- if output_path :
14- dirname = os .path .dirname (output_path )
15- if not os .path .exists (dirname ):
16- os .makedirs (dirname )
17-
18- if (ax or output_path ) and plot is not False :
19- plot = True
20-
21- # Find indices of observations during the reversal ramp
22- ramps = [line for line in voltage_segments if line [2 ] != line [3 ]]
2310
24- # Assume the last ramp is the reversal ramp (convert to ms)
25- tstart , tend = np .array (ramps )[- 1 , :2 ]
11+ def infer_reversal_potential (current , times , voltage_segments , voltages ,
12+ output_path = None , known_Erev = None ,
13+ figsize = (5 , 3 )):
14+ """
15+ Infers a reversal potential in a time series, based on a reversal ramp.
16+
17+ The data is denoised by fitting a 4-th order polynomial through the ramp
18+ data, from which a reversal potential is then detected. If no polynomial
19+ can be fit or the resulting zero-crossing is outside of
20+ ``min(voltages), max(voltages)``, then ``np.nan`` is returned.
21+
22+ @param current: The currents that make up a time series with ``times``
23+ @param times: The sampled times
24+ @param voltage_segments: A list of tuples (tstart, tend, vstart, vend)
25+ describing voltage steps or ramps. It is assumed the final ramp is the
26+ reversal ramp.
27+ @param voltages: The sampled voltages
28+ @param output_path: An optional path to store a plot at
29+ @param known_Erev: A known reversal potential to include in the plot
30+ @param figsize: A size for the plot.
31+
32+ @return: The inferred reversal potential
33+ """
34+
35+ # Get ramp bounds, assuming final ramp is the reversal ramp
36+ tstart , tend = detect_ramp_bounds (times , voltage_segments , - 1 )
2637
2738 istart = np .argmax (times > tstart )
2839 iend = np .argmax (times > tend )
2940
30- times = times [istart :iend ]
3141 current = current [istart :iend ]
3242 voltages = voltages [istart :iend ]
3343
44+ # Fit a 4-th order polynomial
3445 try :
3546 fitted_poly = poly .Polynomial .fit (voltages , current , 4 )
3647 except ValueError as exc :
3748 logging .warning (str (exc ))
3849 return np .nan
3950
51+ # Try extracting the polynomial's roots, accepting only ones that are
52+ # within the range of sampled voltages (so not using ramp info here!)
4053 try :
54+ vmin , vmax = np .min (voltages ), np .max (voltages )
4155 roots = np .unique ([np .real (root ) for root in fitted_poly .roots ()
42- if root > np . min ( voltages ) and root < np . max ( voltages ) ])
56+ if root > vmin and root < vmax ])
4357 except np .linalg .LinAlgError as exc :
4458 logging .warning (str (exc ))
4559 return np .nan
@@ -50,33 +64,30 @@ def infer_reversal_potential(current, times, voltage_segments, voltages,
5064
5165 if len (roots ) == 0 :
5266 return np .nan
67+ erev = roots [- 1 ]
5368
54- if plot :
55- created_fig = False
56- if ax is None and output_path is not None :
57-
58- created_fig = True
59- fig = plt .figure (figsize = figsize )
60- ax = fig .subplots ()
69+ # Optional plot
70+ if output_path is not None :
71+ dirname = os .path .dirname (output_path )
72+ if not os .path .exists (dirname ):
73+ os .makedirs (dirname )
6174
62- ax .set_xlabel ('$V$ (mV)' )
75+ fig = plt .figure (figsize = figsize )
76+ ax = fig .subplots ()
77+ ax .set_xlabel ('$V$ (mV)' ) # Assuming mV here
6378 ax .set_ylabel ('$I$ (nA)' )
6479
6580 # Now plot current vs voltage
6681 ax .plot (voltages , current , 'x' , markersize = 2 , color = 'grey' , alpha = .5 )
67- ax .axvline (roots [ - 1 ] , linestyle = '--' , color = 'grey' , label = r'$E_\mathrm{obs}$' )
82+ ax .axvline (erev , linestyle = '--' , color = 'grey' , label = r'$E_\mathrm{obs}$' )
6883 if known_Erev :
6984 ax .axvline (known_Erev , linestyle = '--' , color = 'orange' ,
7085 label = "Calculated $E_{Kr}$" )
7186 ax .axhline (0 , linestyle = '--' , color = 'grey' )
7287 ax .plot (* fitted_poly .linspace ())
7388 ax .legend ()
7489
75- if output_path is not None :
76- fig = ax .figure
77- fig .savefig (output_path )
78-
79- if created_fig :
80- plt .close (fig )
90+ fig .savefig (output_path )
91+ plt .close (fig )
8192
82- return roots [ - 1 ]
93+ return erev
0 commit comments