diff --git a/docs/source/eFeatures.rst b/docs/source/eFeatures.rst index 9b9314c..4e72c01 100644 --- a/docs/source/eFeatures.rst +++ b/docs/source/eFeatures.rst @@ -2574,8 +2574,10 @@ See voltage clamp example for more details. # keep trace going from stim_start to voltage max max_idx = np.argmax(voltage_interval) - time_interval = time_interval[:max_idx + 1] - voltage_interval = voltage_interval[:max_idx + 1] + min_idx = np.argmin(voltage_interval[:max_idx + 1]) + t_duration = time_interval[max_idx] - time_interval[min_idx] + time_interval = time_interval[min_idx:max_idx + 1] + voltage_interval = voltage_interval[min_idx:max_idx + 1] # correct time so that it starts from 0 time_interval -= time_interval[0] @@ -2585,11 +2587,14 @@ See voltage clamp example for more details. exp_fit, time_interval, voltage_interval, - p0=(1., voltage_interval[-1], voltage_interval[0] - voltage_interval[-1]), + p0=(1., voltage_interval[-1], np.min(voltage_interval) - voltage_interval[-1]), bounds=((0, -np.inf, -np.inf), (np.inf, np.inf, 0)), # positive tau, negative A1 nan_policy="omit", ) - time_constant = np.array([abs(popt[0])]) + time_constant = np.array([abs(popt[0])]) + # correction for very large values dur to noisy data + if abs(popt[0]) > 10 * t_duration: + time_constant = np.array([t_duration]) deactivation_time_constant ~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/efel/pyfeatures/pyfeatures.py b/efel/pyfeatures/pyfeatures.py index 4e7386b..790b907 100644 --- a/efel/pyfeatures/pyfeatures.py +++ b/efel/pyfeatures/pyfeatures.py @@ -27,11 +27,15 @@ """ from efel.pyfeatures.cppfeature_access import _get_cpp_data, get_cpp_feature from efel.pyfeatures.isi import * +import logging from typing_extensions import deprecated import numpy as np from numpy.fft import * + +logger = logging.getLogger(__name__) + all_pyfeatures = [ 'voltage', 'time', @@ -390,8 +394,13 @@ def activation_time_constant() -> np.ndarray | None: # keep trace going from stim_start to voltage max max_idx = np.argmax(voltage_interval) - time_interval = time_interval[:max_idx + 1] - voltage_interval = voltage_interval[:max_idx + 1] + min_idx = np.argmin(voltage_interval[:max_idx + 1]) + t_duration = time_interval[max_idx] - time_interval[min_idx] + if t_duration == 0: + logger.debug("Activation time constant: interval duration is 0, returning [0]") + return np.array([0.0]) + time_interval = time_interval[min_idx:max_idx + 1] + voltage_interval = voltage_interval[min_idx:max_idx + 1] # correct time so that it starts from 0 time_interval -= time_interval[0] @@ -416,6 +425,15 @@ def activation_time_constant() -> np.ndarray | None: except (ValueError, RuntimeError): return None + if abs(popt[0]) > 10 * t_duration: + logger.debug( + "Activation time constant: unexpected time constant value (%s) " + "for smaller interval duration (%s), returning interval duration", + abs(popt[0]), + t_duration, + ) + return np.array([t_duration]) + return np.array([abs(popt[0])]) diff --git a/tests/test_pyfeatures.py b/tests/test_pyfeatures.py index c13733f..8c22f53 100644 --- a/tests/test_pyfeatures.py +++ b/tests/test_pyfeatures.py @@ -403,11 +403,11 @@ def test_activation_time_constant(): feats = efel.get_feature_values(traces, ["activation_time_constant"]) act_tau = [feat_dict["activation_time_constant"][0] for feat_dict in feats] act_tau_ref = [ - 4.430965e-02, 6.308063e-02, 1.354392e+00, 8.766326e-08, - 1.556848e-06, 1.087861e+02, 3.513249e+01, 1.017114e+01, - 4.616776e+00, 2.836394e+00, 1.972983e+00, 1.590110e+00, - 1.21969111e+00, 1.07930254e+00, 8.55734307e-01, 7.39074539e-01, - 6.58848184e-01, 5.96009883e-01 + 0.3026084023253483, 0.06597094092104987, 0.09179018404730191, + 0.20514886735773108, 0.13922091324594776, 110.55045560417942, + 35.48842150253384, 9.878098476029532, 4.391423684073282, 2.6753411591580787, + 1.8958411673573918, 1.5901099575473472, 1.2196911091114777, 1.0793025478535727, + 0.8557343061582321, 0.7390745386928983, 0.6588481846679133, 0.596009882583066, ] numpy.testing.assert_allclose(act_tau, act_tau_ref, rtol=1e-6, atol=1e-10)