From 43ea0ab9697682e8191eb4e4265e2f8914ee4e6f Mon Sep 17 00:00:00 2001 From: Tanguy Damart <53294536+DrTaDa@users.noreply.github.com> Date: Fri, 3 Feb 2023 16:02:02 +0100 Subject: [PATCH] Add step detection sAHP eCode (#124) --- bluepyefe/ecode/ramp.py | 7 ++- bluepyefe/ecode/sAHP.py | 109 +++++++++++++++++++++++++++++++++++----- 2 files changed, 102 insertions(+), 14 deletions(-) diff --git a/bluepyefe/ecode/ramp.py b/bluepyefe/ecode/ramp.py index f1b07234..cd4a84c8 100644 --- a/bluepyefe/ecode/ramp.py +++ b/bluepyefe/ecode/ramp.py @@ -79,7 +79,12 @@ def interpret(self, t, current, config_data, reader_data): # Smooth the current smooth_current = scipy_signal2d(current, 85) - self.set_timing_ecode(["ton", "toff"], config_data) + self.set_timing_ecode(["ton"], config_data) + + if "toff" in config_data and config_data["toff"] is not None: + self.toff = int(round(config_data["toff"] / self.dt)) + else: + self.toff = numpy.argmax(smooth_current[self.ton:]) + self.ton hypamp_value = base_current(current) self.set_amplitudes_ecode("hypamp", config_data, reader_data, hypamp_value) diff --git a/bluepyefe/ecode/sAHP.py b/bluepyefe/ecode/sAHP.py index 430a5879..a4cbdfe3 100644 --- a/bluepyefe/ecode/sAHP.py +++ b/bluepyefe/ecode/sAHP.py @@ -23,6 +23,7 @@ from ..recording import Recording from .tools import scipy_signal2d +from .tools import base_current logger = logging.getLogger(__name__) @@ -79,16 +80,10 @@ def get_stimulus_parameters(self): } return ecode_params - def interpret(self, t, current, config_data, reader_data): - """Analyse a current array and extract from it the parameters - needed to reconstruct the array""" - self.dt = t[1] + def compute_amp(self, current, config_data, reader_data): - # Smooth the current smooth_current = scipy_signal2d(current, 85) - self.set_timing_ecode(["ton", "tmid", "tmid2", "toff"], config_data) - hypamp_value = numpy.median( numpy.concatenate( (smooth_current[: self.ton], smooth_current[self.toff :]) @@ -99,8 +94,8 @@ def interpret(self, t, current, config_data, reader_data): amp_value = numpy.median( numpy.concatenate( ( - smooth_current[self.ton : self.tmid], - smooth_current[self.tmid2 : self.toff], + smooth_current[self.ton: self.tmid], + smooth_current[self.tmid2: self.toff], ) ) ) - self.hypamp @@ -109,10 +104,98 @@ def interpret(self, t, current, config_data, reader_data): amp2_value = numpy.median(smooth_current[self.tmid : self.tmid2]) - self.hypamp self.set_amplitudes_ecode("amp2", config_data, reader_data, amp2_value) - # Converting back to ms - for name_timing in ["ton", "tmid", "tmid2", "toff"]: - self.index_to_ms(name_timing, t) - self.tend = len(t) * self.dt + def step_detection(self, current, config_data, reader_data): + + # Set the threshold to detect the step + noise_level = numpy.std(numpy.concatenate((self.current[:50], self.current[-50:]))) + step_threshold = numpy.clip(4.5 * noise_level, 1e-5, numpy.max(self.current)) + + # The buffer prevent miss-detection of the step when artifacts are + # present at the very start or very end of the current trace + buffer_detect = 2.0 + idx_buffer = int(buffer_detect / self.dt) + idx_buffer = max(1, idx_buffer) + + buffer_step = 50 + smooth_current = scipy_signal2d(current, 85) + + # Infer the beginning of the long step + self.hypamp = base_current(current) + + if "ton" in config_data and config_data["ton"] is not None: + self.ton = int(round(config_data["ton"] / self.dt)) + else: + tmp_current = numpy.abs(smooth_current[idx_buffer:] - self.hypamp) + self.ton = idx_buffer + numpy.argmax(tmp_current > step_threshold) + + # Infer the end of the long step + tmp_current = numpy.flip( + numpy.abs(smooth_current[self.ton:-idx_buffer] - self.hypamp) + ) + self.toff = ( + (len(current) - numpy.argmax(tmp_current > step_threshold)) - 1 - idx_buffer + ) + + # Get the amplitude of the step current (relative to hypamp) + self.amp = numpy.median( + numpy.concatenate((smooth_current[self.ton:self.ton + 50], + smooth_current[self.toff - 50:self.toff])) + ) - self.hypamp + + # Infer the beginning of the short step + tmp_current = numpy.abs( + smooth_current[self.ton + buffer_step:self.toff - buffer_step] - + self.amp - self.hypamp + ) + self.tmid = self.ton + buffer_step + numpy.argmax(tmp_current > step_threshold) + + # Infer the end of the long step + tmp_current = numpy.flip( + numpy.abs( + smooth_current[self.ton + buffer_step:self.toff - 50] - + self.amp - self.hypamp + ) + ) + self.tmid2 = ( + (self.toff - numpy.argmax(tmp_current > step_threshold)) - 1 - buffer_step + ) + + self.amp2 = numpy.median(smooth_current[self.tmid:self.tmid2]) - self.hypamp + + # Converting back ton and toff to ms + self.ton = self.t[int(round(self.ton))] + self.toff = self.t[int(round(self.toff))] + self.tmid = self.t[int(round(self.tmid))] + self.tmid2 = self.t[int(round(self.tmid2))] + self.tend = len(self.t) * self.dt + + # Check for some common step detection failures when the current + # is constant. + if self.ton > self.toff or self.ton > self.tend or \ + self.toff > self.tend or self.tmid == self.ton \ + or self.tmid2 == self.toff: + + self.ton = 0. + self.toff = self.tend + + logger.warning( + "The automatic step detection failed for the recording " + f"{self.protocol_name} in files {self.files}. You should " + "specify ton and toff by hand in your files_metadata " + "for this file." + ) + + def interpret(self, t, current, config_data, reader_data): + """Analyse a current with a step and extract from it the parameters + needed to reconstruct the array""" + + self.dt = t[1] + + if "ton" in config_data and "tmid" in config_data: + self.set_timing_ecode(["ton", "tmid", "tmid2", "toff"], config_data) + self.compute_amp(current, config_data, reader_data) + else: + self.step_detection(current, config_data, reader_data) def generate(self): """Generate the current array from the parameters of the ecode"""