diff --git a/neuroanalysis/event_detection.py b/neuroanalysis/event_detection.py index 9a2e1d4..620d01e 100644 --- a/neuroanalysis/event_detection.py +++ b/neuroanalysis/event_detection.py @@ -95,40 +95,128 @@ def zero_crossing_events(data, min_length=3, min_peak=0.0, min_sum=0.0, noise_th return events +def _deal_unbalanced_initial_off(omit_ends, on_inds, off_inds): + """Deals with situation where there is an "off" crossing from above to below threshold + at the beginning of a trace without there first being an "on" crossing from below to above + threshold. Note that the usage of this function is looking for extreme regions + where a trace is below a negative threshold or above a positive threshold, thus, the + sign of the trace value at *on_inds* and *off_inds* can be positive or negative + """ + if not omit_ends: + on_inds = [0] + on_inds #prepend the edge as on on ind + else: + off_inds = off_inds[1:] #remove the off ind + return on_inds, off_inds + +def _deal_unbalanced_termination_on(omit_ends, on_inds, off_inds, off_to_add): + """Deals with situation where there is an "on" crossing from below to above threshold + toward the end of a trace without an "off" crossing happening thereafter. Note that + the usage of this function is looking for extreme regions + where a trace is below a negative threshold or above a positive threshold, thus, the + sign of the trace value at *on_inds* and *off_inds* can be positive or negative + """ + if not omit_ends: + off_inds = off_inds + [off_to_add] #append the index of the last data point + else: + on_inds = on_inds[:-1] #remove the last on indicie + return on_inds, off_inds + def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_ends=True): """ - Finds regions in a trace that cross a threshold value (as measured by distance from baseline). Returns the index, length, peak, and sum of each event. - Optionally adjusts index to an extrapolated baseline-crossing. + Finds regions in a trace that cross a threshold value (as measured by distance from baseline) and then + recross threshold ('bumps'). If a threshold is crossed at the end of the trace, an event may be excluded + or the beginning/end may be used as the the start/end of the event (depending on the value of *omit_ends*). + + + Parameters + ========== + trace: *Tseries* instance + threshold: float or np.array with dimensions of *trace.data* + Algorithm checks if waveform crosses both positive and negative *threshold* symetrically + around from the y-axis. i.e. if -5. is provided, the algorithm looks for places where + the waveform crosses +/-5. If an array is provided, each index of the *threshold* will + be compared with the data pointwise. + adjust_times: boolean + If True, move the start and end times of the event outward, estimating the zero-crossing point for the event + baseline: float + Value subtracted from the data. + omit_ends: boolean + If true, add the trace endpoint indices to incomplete events, i.e., events that started above threhold at the + beginning of trace, or crossed threshold but did not return below threshold at the end of a trace. If false, + remove the imcomplete events. + + + Returns + ======= + events: numpy structured array. + An event ('bump') is a region of the *trace.data* waveform that crosses above *threshold* and then falls below + threshold again. Each index contains information about an event. Fields as follows: + index: int + index of the initial crossing of the *threshold* + len: int + index length of the event + sum: float + sum of the values in the array between the start and end of the event + peak: float + peak value of event + peak_index: int + index value of the peak + time: float, or np.nan if timing not available + time of the onset of an event + duration: float, or np.nan if timing not available + duration of time of the event + area: float, or np.nan if timing not available + area under the curve of the event + peak_time: float, or np.nan if timing not available + time of peak """ - threshold = abs(threshold) + + data = trace.data data1 = data - baseline - #if (hasattr(data, 'implements') and data.implements('MetaArray')): - ## find all threshold crossings - masks = [(data1 > threshold).astype(np.byte), (data1 < -threshold).astype(np.byte)] + # convert threshold array + if isinstance(threshold, float): + threshold = np.ones(len(data)) * abs(threshold) + + ## find all threshold crossings in both positive and negative directions + ## deal with imcomplete events, and store events + + # -1 (or +1) when crosses from above to below threshold (or visa versa if threshold is negative). Note above threshold refers to value furthest from zero, i.e. it can be positive or negative + masks = [(data1 > threshold).astype(np.byte), (data1 < -threshold).astype(np.byte)] + hits = [] for mask in masks: diff = mask[1:] - mask[:-1] - on_inds = list(np.argwhere(diff==1)[:,0] + 1) - off_inds = list(np.argwhere(diff==-1)[:,0] + 1) - if len(on_inds) == 0 or len(off_inds) == 0: - continue - if off_inds[0] < on_inds[0]: - if omit_ends: - off_inds = off_inds[1:] - if len(off_inds) == 0: - continue - else: - on_inds.insert(0, 0) - if off_inds[-1] < on_inds[-1]: - if omit_ends: - on_inds = on_inds[:-1] - else: - off_inds.append(len(diff)) + # indices where crosses from below to above threshold ('on') + on_inds = list(np.argwhere(diff==1)[:,0] + 1) + # indices where crosses from above to below threshold ('off') + off_inds = list(np.argwhere(diff==-1)[:,0] + 1) + + # deal with cases when there are unmatched on and off indicies + if len(off_inds) > 0: #if there are some off indicies + if len(on_inds) > 0: #and there are also on indicies + if on_inds[0] > off_inds[0]: #check if off happens before on + on_inds, off_inds = _deal_unbalanced_initial_off(omit_ends, on_inds, off_inds) + else: #there are no on indicies + on_inds, off_inds = _deal_unbalanced_initial_off(omit_ends, on_inds, off_inds) + + if len(on_inds) > 0: #if there are some on indicies + if len(off_inds) > 0: #and there are also off indicies + if on_inds[-1] > off_inds[-1]: #check if off happens before on + on_inds, off_inds = _deal_unbalanced_termination_on(omit_ends, on_inds, off_inds, len(data1)) + else: #there are no off indicies + on_inds, off_inds = _deal_unbalanced_termination_on(omit_ends, on_inds, off_inds, len(data1)) + + + # at this point every 'on' should have and 'off' + assert len(on_inds) == len(off_inds) + + # put corresponding on and off indeces in a list for i in range(len(on_inds)): if on_inds[i] == off_inds[i]: + #something wierd happened continue hits.append((on_inds[i], off_inds[i])) @@ -154,7 +242,7 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end ## 2) adjust event times if requested, then recompute parameters for i in range(n_events): ind1, ind2 = hits[i] - ln = ind2-ind1 + ln = ind2 - ind1 ev_data = data1[ind1:ind2] sum = ev_data.sum() if sum > 0: diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py index ed42f59..4a2121e 100644 --- a/neuroanalysis/spike_detection.py +++ b/neuroanalysis/spike_detection.py @@ -50,74 +50,154 @@ def rc_decay(t, tau, Vo): return -(Vo/tau)*np.exp(-t/tau) -def detect_ic_evoked_spikes(trace, pulse_edges, dv2_threshold=40e3, mse_threshold=30., ui=None): - """ +def detect_ic_evoked_spikes(trace, + pulse_edges, + mse_threshold=30., + dv2_event_area=10e-6, + pulse_bounds_move=[.15e-3, 0.02e-3], #0.03e-3 sometimes gets some artifact + double_spike=1.e-3, + ui=None, + artifact_width=.1e-3, #this number seems fairly acurate + dv_thresholds = [125., 15.] #130 is too high + ): + """Return a dict describing an evoked spike in a patch clamp recording. Or Return None + if no spike is detected. + + This function assumes that a square voltage pulse is used to evoke a spike + in a current clamp recording, and that the spike initiation occurs *during* or within a + short region after the stimulation pulse. Note that finding + + + Parameters + ========== + trace: *Trace* instance + The recorded patch clamp data. The recording should be made with a brief pulse intended + to evoke a single spike with short latency. + pulse_edges: (float, float) + The start and end times of the stimulation pulse, relative to the timebase in *trace*. + mse_threshold: float + Value used to determine if there was a spike close to the end of the stimulation pulse. If the + mean square error value of a fit to a RC voltage decay is larger than *mse_threshold* a spike + has happened. + event_area: float + The integral of the 'bump' in dv2 must have at least the following *dv2_event_area*. + pulse_bounds_move: np.array; [float, float] + There are large fluxuations in v, dv, and dv2 around the time of pulse initiation + and termination. *pulse_bounds_move* specifies how much time after the edges of the + stimulation pulse should be considered in the search window. *pulse_bounds_move[0]* is added to + stimulation pulse initiation, *pulse_bounds_move[1]* is added to stimulation pulse termination. + Spikes after the termination of the search window are identified by attempting to fit the dv decay + to the trace that would be seen from standard RC decay (see *mse_threshold*). + double_spike: float + time allowed between two events + ui: + user interface for viewing spike detection + artifact_width: + amount of time that should be not considered for spike metrics after the pulse search window, + i.e. pulse_edges[1] + pulse_bounds_move[1] + artifact_width is the first point where a spike + is searched for after the stimulation pulse + dv_thresholds: [float, float] + *dv_thresholds[0]* (*dv_thresholds[1]*) threshold for potential spike at beginning (end) of pulse window. + + Returns + ======= + spikes: dictionary + contains following information about spikes: + onset_time: float + Onset time of region where searching for spike. Defined as a crossing + of a *threshold* in dv2/dt. + peak_time: float or None + time of spike peak + max_slope_time: float (required output) + time where the voltage is changing most rapidly (i.e. dvdt = 0). Note that + a time during the artifact will never be chosen. The boundry with the largest + slope will be used. Any error due to this estimate will not be + larger than the *artifact_width*. + peak_value: float or None + None or value at *peak_time* + max_slope: float or None + Value at *max_slope_time* """ + if not isinstance(trace, TSeries): + raise TypeError("data must be Trace instance.") + if ui is not None: ui.clear() ui.console.setStack() ui.plt1.plot(trace.time_values, trace.data) assert trace.data.ndim == 1 - pulse_edges = tuple(map(float, pulse_edges)) # make sure pulse_edges is (float, float) - + + pulse_edges = np.array(tuple(map(float, pulse_edges))) # confirms pulse_edges is [float, float] + window_edges = pulse_edges + pulse_bounds_move + + #--------------------------------------------------------- + #----this is were vc and ic code diverge------------------ + #--------------------------------------------------------- + # calculate derivatives within pulse window - diff1 = trace.time_slice(*pulse_edges).diff() - diff2 = diff1.diff() + diff1 = trace.time_slice(window_edges[0], window_edges[1]).diff() + # note that this isnt actually used for anything other than plotting + diff2 = diff1.diff() # mask out pulse artifacts in diff2 before lowpass filtering for edge in pulse_edges: apply_cos_mask(diff2, center=edge + 100e-6, radius=400e-6, power=2) # low pass filter the second derivative diff2 = bessel_filter(diff2, 10e3, order=4, bidir=True) + # low pass filter the second derivative + + diff1 = bessel_filter(diff1, 10e3, order=4, bidir=True) - # look for positive bumps in second derivative - events2 = list(threshold_events(diff2 / dv2_threshold, threshold=1.0, adjust_times=False)) + # create a linear threshold + dv_threshold = dv_thresholds[0] - (dv_thresholds[0] - dv_thresholds[1])/len(diff1.data) * np.arange(len(diff1.data)) + events1 = list(threshold_events(diff1, + threshold=dv_threshold, + adjust_times=False, + omit_ends=False)) if ui is not None: ui.plt2.plot(diff1.time_values, diff1.data) + ui.plt2.plot(diff1.time_values, dv_threshold, pen='y') ui.plt3.plot(diff2.time_values, diff2.data) - ui.plt3.addLine(y=dv2_threshold) - # for each bump in d2, either discard the event or generate spike metrics + # for each bump in dvdt, either discard the event or generate + # spike metrics from v and dvdt spikes = [] - for ev in events2: - total_area = ev['area'] - onset_time = ev['time'] - # ignore events near pulse offset - if abs(onset_time - pulse_edges[1]) < 50e-6: + slope_hit_boundry = False + for ev in events1: + # require 3 indexes are above threshold + if ev['len'] < 3: continue - - # require dv2 bump to be positive, not tiny - if total_area < 10e-6: - continue - - # don't double-count spikes within 1 ms - if len(spikes) > 0 and onset_time < spikes[-1]['onset_time'] + 1e-3: + # require events to be positive + if ev['sum'] < 0: continue - max_slope_window = onset_time, pulse_edges[1]-50e-6 - max_slope_chunk = diff1.time_slice(*max_slope_window) - if len(max_slope_chunk) == 0: + onset_time = ev['time'] #this number is arbitrary + + # don't double-count spikes + if len(spikes) > 0 and onset_time < spikes[-1]['onset_time'] + double_spike: continue - max_slope_idx = np.argmax(max_slope_chunk.data) - max_slope_time = max_slope_chunk.time_at(max_slope_idx) - max_slope_time, is_edge = max_time(diff1.time_slice(onset_time, pulse_edges[1] - 50e-6)) + #max slope must be within the event. + max_slope_time, is_edge = max_time(diff1.time_slice(onset_time, min(onset_time + ev['duration'] + diff1.dt, window_edges[1]))) #window edges will not be relevant if the duration is used max_slope = diff1.value_at(max_slope_time) - # require dv/dt to be above a threshold value - if max_slope <= 30: # mV/ms + + if is_edge ==-1: + #if it is an edge at the beginning of the pulse V is most likely just increasing from pulse initiation continue - if is_edge != 0: - # can't see max slope - max_slope_time = None - max_slope = None - peak_time, is_edge = max_time(trace.time_slice(onset_time, pulse_edges[1] + 2e-3)) - if is_edge != 0 or pulse_edges[1] < peak_time < pulse_edges[1] + 50e-6: - # peak is obscured by pulse edge - peak_time = None + #this should only happen at the end of pulse window since we are looking at dvdt in events + if is_edge == 1: + slope_hit_boundry = True + peak_time = None # if slope is hitting the pulse window boundry don't trust the peak + else: + slope_hit_boundry = False + peak_time, is_edge = max_time(trace.time_slice(onset_time, min(max_slope_time + 1.e-3, window_edges[1]))) + if is_edge != 0 or peak_time > window_edges[1]: + # peak is obscured by pulse edge + peak_time = None spikes.append({ 'onset_time': onset_time, @@ -127,11 +207,21 @@ def detect_ic_evoked_spikes(trace, pulse_edges, dv2_threshold=40e3, mse_threshol 'max_slope': max_slope, }) - # if no spike was found in the pulse region check to see if there is a spike in the pulse termination region - if len(spikes) == 0: - # note that this is using the dvdt with the termination artifact in it to locate where it should start - dv_after_pulse = trace.time_slice(pulse_edges[1] + 100e-6, None).diff() - dv_after_pulse = bessel_filter(dv_after_pulse, 15e3, bidir=True) + + """check for evidence of spike in the decay after the pulse if + 1. there are no previous spikes in the trace. Note this is specifically here so don't get an error from spikes[-1] if it is empty + 2. there are no previous spikes within 1 ms of the boundry + 3. a potential spike was found but a clear max slope was not found because it is hitting window boundry (set where + artifact begins). However, if the max_slope is > 80 it is likely that the spike is too fast to observe in the decay so the values + in the spike dictionary remain. If the max_slope at the boundry is less than 80, the algorigthm checks to see if a spike is + detected in the decay. If the slope is not > 80 it checks to see if the max slope after the artifact is also a boundry, if it is + it takes the largest value of the two, if not it replaces the last spike with new values + """ + if (len(spikes) == 0) or (slope_hit_boundry and (spikes[-1]['max_slope'] < 80.)) or (spikes[-1]['max_slope_time'] < (window_edges[1] - double_spike)): #last spike is more than 1 ms before end + + #TODO: resolve different trace types and filtering + dv_after_pulse = trace.time_slice(window_edges[1] + artifact_width, None).diff() #taking non filtered data + dv_after_pulse = bessel_filter(dv_after_pulse, 15e3, bidir=True) #does this really need a different filtering # create a vector to fit ttofit = dv_after_pulse.time_values # setting time to start at zero, note: +1 because time trace of derivative needs to be one shorter @@ -140,34 +230,84 @@ def detect_ic_evoked_spikes(trace, pulse_edges, dv2_threshold=40e3, mse_threshol # do fit and see if it matches popt, pcov = curve_fit(rc_decay, ttofit, dv_after_pulse.data, maxfev=10000) fit = rc_decay(ttofit, *popt) - if ui is not None: - ui.plt2.plot(dv_after_pulse.time_values, dv_after_pulse.data) - ui.plt2.plot(dv_after_pulse.time_values, fit, pen='b') - diff = dv_after_pulse - fit mse = (diff.data**2).mean() # mean squared error + if ui is not None: + ui.plt2.plot(dv_after_pulse.time_values, dv_after_pulse.data, pen='b') + if mse > mse_threshold: + ui.plt2.plot(dv_after_pulse.time_values, fit, pen='g') + else: + ui.plt2.plot(dv_after_pulse.time_values, fit, pen='r') + + # there is a spike, so identify paramters if mse > mse_threshold: - search_window = 2e-3 - max_slope_time, is_edge = max_time(diff.time_slice(pulse_edges[1], pulse_edges[1] + search_window)) - if is_edge != 0: - max_slope_time = None - peak_time, is_edge = max_time(trace.time_slice(max_slope_time or pulse_edges[1] + 100e-6, pulse_edges[1] + search_window)) - if is_edge != 0: + onset_time = window_edges[1] + artifact_width + max_slope_time, slope_is_edge = max_time(dv_after_pulse.time_slice(onset_time, onset_time + 0.5e-3)) #can't search far out because slope might have dropped and then peak will be found far away. + max_slope = dv_after_pulse.value_at(max_slope_time) + peak_time, peak_is_edge = max_time(trace.time_slice(max_slope_time, max_slope_time + 1e-3)) + peak_value = dv_after_pulse.value_at(peak_time) +# import pdb; pdb.set_trace() + if peak_is_edge != 0: peak_time = None + peak_value = None + + if len(spikes) == 0: #if there are no previous values in the spike array + # append the newly found values + pass + else: #there was a value in spike array + #if max slopes are on the boundry on each side of the artifact + if slope_is_edge and slope_hit_boundry: + # check if the slope before the artifact was larger than after the artifact + if spikes[-1]['max_slope'] > max_slope: + # set values to the values before artifact + onset_time = spikes[-1]['onset_time'] # set the onset to be before the artifact + peak_time = spikes[-1]['peak_time'] + peak_value = spikes[-1]['peak_value'] + max_slope = spikes[-1]['max_slope'] + max_slope_time = spikes[-1]['max_slope_time'] + spikes.pop(-1) # get rid of the last values for spike because they will be updated + + # if there is an obvous minimum after the artifact and there is a previous value in the spike array that hit a boundry + elif slope_hit_boundry and not slope_is_edge: + # get rid of the last recorded spike so it can be replaced + spikes.pop(-1) #so get rid of the last value in spike set off before artifact but not completed because it hit the boundry looking for slope time + + #there was not a spike close to the boundry so no need to get rid of last spike + elif not slope_hit_boundry and slope_is_edge: + pass + elif not slope_hit_boundry and not slope_is_edge: + #there is no spike close to the boundry and there is an obvious min after the artifact + pass + else: + raise Exception('this has not been accounted for') + + # add the found spike to the end spikes.append({ - 'onset_time': None, + 'onset_time': onset_time, 'max_slope_time': max_slope_time, 'peak_time': peak_time, - 'peak_value': None if peak_time is None else trace.value_at(peak_time), - 'max_slope': None if max_slope_time is None else dv_after_pulse.value_at(max_slope_time), + 'peak_value':peak_value, + 'max_slope': max_slope }) + + else: + # there is no spike found from the decay + if slope_hit_boundry: + # if last spike recorded was against the bound delete it + spikes.pop(-1) - for spike in spikes: + + for spike in spikes: + # max_slope_time is how we define spike initiation, so, make sure it is defined. assert 'max_slope_time' in spike return spikes -def detect_vc_evoked_spikes(trace, pulse_edges, ui=None): +def detect_vc_evoked_spikes(trace, + pulse_edges, + pulse_bounds_move=[.1e-3, 0.02e-3], #0.03e-3 sometimes gets some artifact + dv2_threshold=.02, + ui=None): """Return a dict describing an evoked spike in a patch clamp recording, or None if no spike is detected. This function assumes that a square voltage pulse is used to evoke an unclamped spike @@ -181,6 +321,26 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None): to evoke a single spike with short latency. pulse_edges : (float, float) The start and end times of the stimulation pulse, relative to the timebase in *trace*. + pulse_bounds_move: np.array; [float, float] + There are large fluxuations in v, dv, and dv2 around the time of pulse initiation + and termination. *pulse_bounds_move* specifies how much time after the edges of the stimulation pulse should be considered in the search window. *pulse_bounds_move[0]* is added to stimulation pulse initiation, *pulse_bounds_move[1]* is added to stimulation pulse termination. + ui: + user interface for viewing spike detection + + Returns + ======= + spikes: dictionary + contains following information about spikes: + onset_time: float + Onset time of region where searching for a spike (identified in event_detection + module). Defined as a crossing of a *threshold* or *baseline* in dv2/dt. + peak_time: float + time of spike peak + max_slope_time: + time where the voltage is changing most rapidly (i.e. max or min of dvdt) + peak_value: float + None if *peak_time* is None, else *trace.value_at(peak_time)* + max_slope: float """ if not isinstance(trace, TSeries): raise TypeError("data must be TSeries instance.") @@ -191,61 +351,68 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None): ui.plt1.plot(trace.time_values, trace.data) assert trace.ndim == 1 - pulse_edges = tuple(map(float, pulse_edges)) # make sure pulse_edges is (float, float) - - diff1 = trace.time_slice(pulse_edges[0], pulse_edges[1] + 2e-3).diff() + pulse_edges = np.array(tuple(map(float, pulse_edges))) # make sure pulse_edges is (float, float) + window_edges = pulse_edges + pulse_bounds_move + + # calculate derivatives within pulse window + diff1 = trace.time_slice(window_edges[0], window_edges[1]).diff() diff2 = diff1.diff() # crop and filter diff1 - diff1 = diff1.time_slice(pulse_edges[0] + 100e-6, pulse_edges[1]) diff1 = bessel_filter(diff1, cutoff=20e3, order=4, btype='low', bidir=True) # crop and low pass filter the second derivative diff2 = diff2.time_slice(pulse_edges[0] + 150e-6, pulse_edges[1]) - diff2 = bessel_filter(diff2, 20e3, order=4, bidir=True) - # chop off ending transient - diff2 = diff2.time_slice(None, diff2.t_end - 100e-6) + diff2 = bessel_filter(diff2, 10e3, order=4, bidir=True) # look for negative bumps in second derivative - # dv1_threshold = 1e-6 - dv2_threshold = 0.02 - events = list(threshold_events(diff2 / dv2_threshold, threshold=1.0, adjust_times=False, omit_ends=True)) + events = list(threshold_events(diff2 / dv2_threshold, + threshold=1.0, adjust_times=False, omit_ends=True)) + if ui is not None: ui.plt2.plot(diff1.time_values, diff1.data) - # ui.plt2.plot(diff1_hp.time_values, diff1.data) - # ui.plt2.addLine(y=-dv1_threshold) ui.plt3.plot(diff2.time_values, diff2.data) ui.plt3.addLine(y=dv2_threshold) if len(events) == 0: return [] + # for each bump in d2vdt, either discard the event or generate + # spike metrics from v and dvdt spikes = [] for ev in events: + # require 3 indexes are above threshold + if ev['len'] < 3: + continue + # require events to be positive + if np.abs(ev['sum']) < 2.: + continue if ev['sum'] > 0 and ev['peak'] < 5. and ev['time'] < diff2.t0 + 60e-6: # ignore positive bumps very close to the beginning of the trace continue if len(spikes) > 0 and ev['peak_time'] < spikes[-1]['max_slope_time'] + 1e-3: # ignore events that follow too soon after a detected spike continue - if ev['sum'] < 0: - onset_time = ev['peak_time'] - search_time = onset_time + # only accept positive bumps + continue else: - search_time = ev['time'] - 200e-6 - onset_time = None + search_time = ev['time'] - 100e-6 + onset_time = search_time + max_slope_rgn = diff1.time_slice(search_time, search_time + 0.5e-3) - max_slope_time, is_edge = min_time(max_slope_rgn) + max_slope_time, is_edge = min_time(max_slope_rgn) #note this is looking for min because event must be negative in VC max_slope = diff1.value_at(max_slope_time) + if max_slope > 0: # actual slope must be negative at this point # (above we only tested the sign of the high-passed event) continue peak_search_rgn = trace.time_slice(max_slope_time, min(pulse_edges[1], search_time + 1e-3)) + if len(peak_search_rgn) == 0: peak = None peak_time = None @@ -265,6 +432,29 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None): 'peak_value': peak, }) + # remove spikes where the same values were found from two different events + # it is probable that this would no happen in negative bumps in dv2 where + # ignored. + # if len(spikes) > 1: + # slopes=[] + + # # find unique values of slope + # for spike in spikes: + # slopes.append(spike['max_slope_time']) + # uq = np.unique(slopes) + + # # if there are less unique values than spikes there are duplicates + # if len(uq)!=len(spikes): + # #find indicies to remove + # remove_indicies=[] + # for ii in range(1,len(slopes)): + # for jj in range(ii): + # if slopes[ii] == slopes[jj]: + # remove_indicies.append(ii) + # #remove the duplicate spikes + # for ii in remove_indicies: + # del spikes[ii] + for spike in spikes: assert 'max_slope_time' in spike return spikes diff --git a/neuroanalysis/ui/spike_detection.py b/neuroanalysis/ui/spike_detection.py index 28c8cfe..61f277a 100644 --- a/neuroanalysis/ui/spike_detection.py +++ b/neuroanalysis/ui/spike_detection.py @@ -37,11 +37,11 @@ def show_result(self, spikes): continue for spike in spikes: if spike.get('onset_time') is not None: - plt.addLine(x=spike['onset_time']) + plt.addLine(x=spike['onset_time'], pen='g') if spike.get('max_slope_time') is not None: plt.addLine(x=spike['max_slope_time'], pen='b') if spike.get('peak_time') is not None: - plt.addLine(x=spike['peak_time'], pen='g') + plt.addLine(x=spike['peak_time'], pen='r') class SpikeDetectTestUi(UserTestUi):