From ba4377c9adb7e3b9d65fda61153edd5ae3477e55 Mon Sep 17 00:00:00 2001
From: Corinne <corinnet@alleninstitute.org>
Date: Sun, 4 Aug 2019 22:22:08 -0700
Subject: [PATCH 01/10] mostly commenting and code restructure

---
 neuroanalysis/event_detection.py | 68 +++++++++++++++++++++-----
 neuroanalysis/spike_detection.py | 83 +++++++++++++++++++++++++++-----
 2 files changed, 125 insertions(+), 26 deletions(-)

diff --git a/neuroanalysis/event_detection.py b/neuroanalysis/event_detection.py
index cd061e6..9cdf4d6 100644
--- a/neuroanalysis/event_detection.py
+++ b/neuroanalysis/event_detection.py
@@ -98,24 +98,66 @@ def zero_crossing_events(data, min_length=3, min_peak=0.0, min_sum=0.0, noise_th
 
 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*). 
+    Optionally adjusts start and end index of an event to an extrapolated baseline-crossing and calculates 
+    values. 
+
+    Parameters                                                                                                                                     
+    ==========
+    trace: *Trace* instance                            
+    threshold: float
+        algorithm checks if waveform crosses both positive and negative thresholds.  
+        i.e. if -5. is provided, the algorithm looks for places where the waveform crosses +/-5 
+    adjust_times: boolean
+        if True, move the start and end times of the event outward, estimating the zero-crossing point for the event
+    
+    Returns
+    =======
+    events: numpy structured array.  
+        An event is a region of the *Trace.data* waveform that crosses above *threshold* and then falls below threshold again
+        Referred to as a 'bump'.  There are additional criteria (not listed here) for a bump to be considered an event.    
+        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
+           length 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)]
+    ## find all positive and negative threshold crossings of baseline adjusted data
+    masks = [(data1 > threshold).astype(np.byte), (data1 < -threshold).astype(np.byte)] # 1 where data is [above threshold, below negative threshold]
     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)
+        diff = mask[1:] - mask[:-1] # -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
+        
+        on_inds = list(np.argwhere(diff==1)[:,0] + 1) #indices where crosses from below to above threshold
+        off_inds = list(np.argwhere(diff==-1)[:,0] + 1) #indices where crosses from above to below threshold
+#        import pdb; pdb.set_trace()
         if len(on_inds) == 0 or len(off_inds) == 0:
             continue
-        if off_inds[0] < on_inds[0]:
+
+        ## if there are unequal number of crossing from one direction, either remove the ends or add the appropriate initial or end index (which will be the beginning or end of the waveform)
+        if off_inds[0] < on_inds[0]:  
             if omit_ends:
                 off_inds = off_inds[1:]
                 if len(off_inds) == 0:
@@ -128,7 +170,7 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
             else:
                 off_inds.append(len(diff))
         for i in range(len(on_inds)):
-            if on_inds[i] == off_inds[i]:
+            if on_inds[i] == off_inds[i]: #if an index is "crossing in both directions" get rid of it
                 continue
             hits.append((on_inds[i], off_inds[i]))
     
@@ -168,12 +210,12 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
         if adjust_times:  ## Move start and end times outward, estimating the zero-crossing point for the event
         
             ## adjust ind1 first
-            mind = np.argmax(ev_data)
-            pdiff = abs(peak - ev_data[0])
-            if pdiff == 0:
+            mind = np.argmax(ev_data)  # max of whole trace
+            pdiff = abs(peak - ev_data[0])  # find the how high the peak is from the front of event 
+            if pdiff == 0:  
                 adj1 = 0
             else:
-                adj1 = int(threshold * mind / pdiff)
+                adj1 = int(threshold * mind / pdiff)  # (max value of whole trace)* 1/(hight of peak from first data point)
                 adj1 = min(ln, adj1)
             ind1 -= adj1
             
diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py
index a8ad480..141b3be 100644
--- a/neuroanalysis/spike_detection.py
+++ b/neuroanalysis/spike_detection.py
@@ -49,9 +49,52 @@ 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, 
+                            dv2_threshold=40.e3, 
+                            mse_threshold=30., 
+                            dv2_event_area=10e-6,
+                            pulse_term_bound=50.e-6,
+                            double_spike=1.e-3,
+                            ui=None):
+    """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.  Currently, if a spike is detected in the region 
+    shortly after the stimulation pulse it is identified as a spike but the initiation index/time
+    is not resolved.
+
+    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*.
+    dv2_threshold: float
+        Value for the second derivative of the voltage must cross to be considered as a possible spike.  
+    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_term_boundry: float
+        There are large fluxuations in v, dv, and dv2 around the time of pulse initiation
+        and termination.  *pulse_term_boundry* specifies how much time before the termination 
+        of the stimulation pulse should not be considered in the pulse window.  Spikes after this
+        point 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 between 
+    ui: 
+        user interface for viewing spike detection
     """
+    if not isinstance(trace, Trace):
+        raise TypeError("data must be Trace instance.")
+
     if ui is not None:
         ui.clear()
         ui.console.setStack()
@@ -60,6 +103,10 @@ def detect_ic_evoked_spikes(trace, pulse_edges, dv2_threshold=40e3, mse_threshol
     assert trace.data.ndim == 1
     pulse_edges = tuple(map(float, pulse_edges))  # make sure pulse_edges is (float, float)
     
+    #---------------------------------------------------------
+    #----this is were vc and ic code diverge------------------
+    #---------------------------------------------------------
+
     # calculate derivatives within pulse window
     diff1 = trace.time_slice(*pulse_edges).diff()
     diff2 = diff1.diff()
@@ -79,32 +126,33 @@ def detect_ic_evoked_spikes(trace, pulse_edges, dv2_threshold=40e3, mse_threshol
         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 d2vdt, 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:
+        if abs(onset_time - pulse_edges[1]) < pulse_term_bound:
             continue
 
         # require dv2 bump to be positive, not tiny
-        if total_area < 10e-6:
+        if total_area < dv2_event_area:
             continue
         
-        # don't double-count spikes within 1 ms
-        if len(spikes) > 0 and onset_time < spikes[-1]['onset_time'] + 1e-3:
+        # don't double-count spikes 
+        if len(spikes) > 0 and onset_time < spikes[-1]['onset_time'] + double_spike:
             continue
 
-        max_slope_window = onset_time, pulse_edges[1]-50e-6
+        max_slope_window = onset_time, pulse_edges[1] - pulse_term_bound
         max_slope_chunk = diff1.time_slice(*max_slope_window)
         if len(max_slope_chunk) == 0:
             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_time, is_edge = max_time(diff1.time_slice(onset_time, pulse_edges[1] - pulse_term_bound))
         max_slope = diff1.value_at(max_slope_time)
         # require dv/dt to be above a threshold value
         if max_slope <= 30:  # mV/ms
@@ -114,7 +162,7 @@ def detect_ic_evoked_spikes(trace, pulse_edges, dv2_threshold=40e3, mse_threshol
             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:
+        if is_edge != 0 or pulse_edges[1] < peak_time < pulse_edges[1] + pulse_term_bound:
             # peak is obscured by pulse edge
             peak_time = None
         
@@ -192,6 +240,11 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None):
     assert trace.ndim == 1
     pulse_edges = tuple(map(float, pulse_edges))  # make sure pulse_edges is (float, float)
 
+    #---------------------------------------------------------
+    #----this is were vc and ic code diverge------------------
+    #---------------------------------------------------------
+
+    # calculate derivatives within pulse window
     diff1 = trace.time_slice(pulse_edges[0], pulse_edges[1] + 2e-3).diff()
     diff2 = diff1.diff()
 
@@ -207,19 +260,23 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None):
 
     # look for negative bumps in second derivative
     # dv1_threshold = 1e-6
-    dv2_threshold = 0.02
+    dv2_threshold = 0.02 
     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)
+        # ui.plt3.plot(diff2.time_values, diff2.data)
+        # ui.plt3.addLine(y=dv2_threshold)
+        ui.plt3.plot(diff2.time_values, diff2.data/dv2_threshold)
+        ui.plt3.addLine(y=1)
 
     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:
         if ev['sum'] > 0 and ev['peak'] < 5. and ev['time'] < diff2.t0 + 60e-6:

From 9486e41e1e1891043d281f75552268608a71a7eb Mon Sep 17 00:00:00 2001
From: Corinne <corinnet@alleninstitute.org>
Date: Thu, 8 Aug 2019 13:20:24 -0700
Subject: [PATCH 02/10] bug fix in ui.spike_detection, also lots of pdb
 comments that will need to be removed

---
 neuroanalysis/event_detection.py    | 77 +++++++++++++++++++++++++----
 neuroanalysis/spike_detection.py    | 67 +++++++++++++++++++++----
 neuroanalysis/ui/spike_detection.py | 12 +++--
 3 files changed, 132 insertions(+), 24 deletions(-)

diff --git a/neuroanalysis/event_detection.py b/neuroanalysis/event_detection.py
index 9cdf4d6..fc0e165 100644
--- a/neuroanalysis/event_detection.py
+++ b/neuroanalysis/event_detection.py
@@ -96,7 +96,7 @@ def zero_crossing_events(data, min_length=3, min_peak=0.0, min_sum=0.0, noise_th
     return events
 
 
-def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_ends=True):
+def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_ends=True, debug=False):
     """
     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
@@ -117,12 +117,12 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
     =======
     events: numpy structured array.  
         An event is a region of the *Trace.data* waveform that crosses above *threshold* and then falls below threshold again
-        Referred to as a 'bump'.  There are additional criteria (not listed here) for a bump to be considered an event.    
+        Sometimes referred to as a 'bump'.  There are additional criteria (not listed here) for a bump to be considered an event.    
         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
+        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
@@ -132,7 +132,7 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
         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
-           length of time of the event
+           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
@@ -145,14 +145,35 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
     #if (hasattr(data, 'implements') and data.implements('MetaArray')):
     
     ## find all positive and negative threshold crossings of baseline adjusted data
+
+    if debug:
+        import pdb; pdb.set_trace()
+        # FYI: can't use matplot lib before debugger is on
+        # type *continue to see plot
+        import matplotlib.pyplot as mpl
+
+
     masks = [(data1 > threshold).astype(np.byte), (data1 < -threshold).astype(np.byte)] # 1 where data is [above threshold, below negative threshold]
     hits = []
     for mask in masks:
         diff = mask[1:] - mask[:-1] # -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
         
-        on_inds = list(np.argwhere(diff==1)[:,0] + 1) #indices where crosses from below to above threshold
-        off_inds = list(np.argwhere(diff==-1)[:,0] + 1) #indices where crosses from above to below threshold
-#        import pdb; pdb.set_trace()
+        # TODO: It might be a good idea to make offindexes one less so that region above threshold looks symmetrical
+        # find start and end inicies (above threshold) where waveform is above threshold
+        on_inds = list(np.argwhere(diff==1)[:,0] + 1) #where crosses from below to above threshold Note taking index after this happens
+        off_inds = list(np.argwhere(diff==-1)[:,0]) #where crosses from above to below threshold. Note taking index before this happens
+
+        if debug:
+            mpl.figure()
+            mpl.plot(data1, '.-') 
+            for on in on_inds:
+                mpl.axvline(x=on, color='g', linestyle='--')
+            for off in off_inds:
+                mpl.axvline(x=off, color='r', linestyle='--')
+            mpl.axhline(y=threshold, color='y', linestyle='--')
+            mpl.axhline(y=-threshold, color='y', linestyle='--')
+            mpl.show(block=False)
+
         if len(on_inds) == 0 or len(off_inds) == 0:
             continue
 
@@ -169,13 +190,34 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
                 on_inds = on_inds[:-1]
             else:
                 off_inds.append(len(diff))
+        
+        # Add both events above +threshold and those below -threshold to a list
         for i in range(len(on_inds)):
-            if on_inds[i] == off_inds[i]: #if an index is "crossing in both directions" get rid of it
+            
+            # remove any point where an on off is seperated by less than 2 indicies
+            if (on_inds[i] + 1) == off_inds[i]: 
+                continue
+            if (on_inds[i] - 1) == off_inds[i]: 
+                continue
+            if on_inds[i] == off_inds[i]: 
                 continue
+            
             hits.append((on_inds[i], off_inds[i]))
     
     ## sort hits  ## NOTE: this can be sped up since we already know how to interleave the events..
     hits.sort(key=lambda a: a[0])
+    if debug is True:
+        # FYI: can't use matplot lib before debugger is on
+        # type *continue to see plot
+        mpl.figure()
+        mpl.title('first round before adjustment')
+        mpl.plot(data1, '.-') 
+        for hit in hits:
+            mpl.axvline(x=hit[0], color='g', linestyle='--')
+            mpl.axvline(x=hit[1], color='r', linestyle='--')
+        mpl.axhline(y=threshold, color='y', linestyle='--')
+        mpl.axhline(y=-threshold, color='y', linestyle='--')
+        mpl.show(block=False)
     
     n_events = len(hits)
     events = np.empty(n_events, dtype=[
@@ -204,7 +246,7 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
         else:
             peak_ind = np.argmin(ev_data)
         peak = ev_data[peak_ind]
-        peak_ind += ind1
+        peak_ind += ind1 # adjust peak_ind from local event data to entire waveform
             
         #print "event %f: %d" % (xvals[ind1], ind1) 
         if adjust_times:  ## Move start and end times outward, estimating the zero-crossing point for the event
@@ -275,6 +317,18 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
     
         ## remove masked events
         events = events[mask]
+        
+        if debug:
+            mpl.figure()
+            mpl.title('adjusted')
+            mpl.plot(data1, '.-') 
+            for on in events['index']:
+                mpl.axvline(x=on, color='g', linestyle='--')
+            # for off in off_inds:
+            #     mpl.axvline(x=off, color='r', linestyle='--')
+            mpl.axhline(y=threshold, color='y', linestyle='--')
+            mpl.axhline(y=-threshold, color='y', linestyle='--')
+            mpl.show(block=False)
 
     # add in timing information if available:
     if trace.has_timing:
@@ -291,6 +345,9 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
         ev['area'] = np.nan
         ev['peak_time'] = np.nan
 
+
+    if debug:
+        pdb.set_trace()
     return events
 
 
diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py
index 141b3be..177f2a4 100644
--- a/neuroanalysis/spike_detection.py
+++ b/neuroanalysis/spike_detection.py
@@ -91,6 +91,22 @@ def detect_ic_evoked_spikes(trace,
         time between 
     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 spike. 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. dvdt = 0) 
+        peak_value: float
+            None if peak_time is None else trace.value_at(peak_time),
+        max_slope: float
+
     """
     if not isinstance(trace, Trace):
         raise TypeError("data must be Trace instance.")
@@ -101,7 +117,7 @@ def detect_ic_evoked_spikes(trace,
         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 = tuple(map(float, pulse_edges))  # confirms pulse_edges is (float, float)
     
     #---------------------------------------------------------
     #----this is were vc and ic code diverge------------------
@@ -228,6 +244,21 @@ 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*.
+    
+    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, Trace):
         raise TypeError("data must be Trace instance.")
@@ -260,25 +291,36 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None):
 
     # 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))
+    dv2_threshold = 0.02
+
+    #=========================================================================
+    debug = False
+    #=========================================================================
+
+    events = list(threshold_events(diff2 / dv2_threshold, 
+        threshold=1.0, adjust_times=False, omit_ends=True, debug=debug))
+
+
 
     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)
-        ui.plt3.plot(diff2.time_values, diff2.data/dv2_threshold)
-        ui.plt3.addLine(y=1)
+        ui.plt3.plot(diff2.time_values, diff2.data)
+        ui.plt3.addLine(y=dv2_threshold)
+        # ui.plt3.plot(diff2.time_values, diff2.data/dv2_threshold)
+        # ui.plt3.addLine(y=1)
 
     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:
+        if np.abs(ev['sum']) <10:
+            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
@@ -286,22 +328,25 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None):
             # ignore events that follow too soon after a detected spike
             continue
 
+        #TODO: What is this doing?
         if ev['sum'] < 0:
             onset_time = ev['peak_time']
             search_time = onset_time
         else:
             search_time = ev['time'] - 200e-6
-            onset_time = None
+            onset_time = ev['peak_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
@@ -312,6 +357,10 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None):
                 peak_time = None
             else:
                 peak = trace.time_at(peak_time)
+        #============================================================================
+        if debug:
+            import pdb; pdb.set_trace()
+        #============================================================================
 
         spikes.append({
             'onset_time': onset_time,
diff --git a/neuroanalysis/ui/spike_detection.py b/neuroanalysis/ui/spike_detection.py
index b3cea99..78a4743 100644
--- a/neuroanalysis/ui/spike_detection.py
+++ b/neuroanalysis/ui/spike_detection.py
@@ -33,13 +33,15 @@ def clear(self):
 
     def show_result(self, spikes):
         for plt in [self.plt1, self.plt2, self.plt3]:
+            if not 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(object):
@@ -47,8 +49,7 @@ def __init__(self):
         pg.mkQApp()
 
         self.widget = pg.QtGui.QSplitter(pg.QtCore.Qt.Horizontal)
-        self.widget.resize(1600, 1000)
-
+        self.widget.resize(2000, 1000)
         self.display1 = SpikeDetectUI('expected result')
         self.display2 = SpikeDetectUI('current result')
         self.widget.addWidget(self.display1.widget)
@@ -70,7 +71,7 @@ def __init__(self):
         self.fail_btn.clicked.connect(self.fail_clicked)
 
         self.last_btn_clicked = None
-        self.widget.setSizes([400, 400, 800])
+        self.widget.setSizes([800, 800, 400])
 
     def pass_clicked(self):
         self.last_btn_clicked = 'pass'
@@ -81,6 +82,7 @@ def fail_clicked(self):
     def user_passfail(self):
         self.widget.show()
         while True:
+            print('in user_passfail')
             pg.QtGui.QApplication.processEvents()
             last_btn_clicked = self.last_btn_clicked
             self.last_btn_clicked = None

From 6a060d7f3dcf9a2af9711baf46e1efc1dfa963bb Mon Sep 17 00:00:00 2001
From: Corinne <corinnet@alleninstitute.org>
Date: Thu, 8 Aug 2019 14:44:50 -0700
Subject: [PATCH 03/10] merge save

---
 neuroanalysis/ui/spike_detection.py | 120 +---------------------------
 1 file changed, 1 insertion(+), 119 deletions(-)

diff --git a/neuroanalysis/ui/spike_detection.py b/neuroanalysis/ui/spike_detection.py
index 5483804..61f277a 100644
--- a/neuroanalysis/ui/spike_detection.py
+++ b/neuroanalysis/ui/spike_detection.py
@@ -33,7 +33,7 @@ def clear(self):
 
     def show_result(self, spikes):
         for plt in [self.plt1, self.plt2, self.plt3]:
-            if not spikes:
+            if spikes is None:
                 continue
             for spike in spikes:
                 if spike.get('onset_time') is not None:
@@ -48,124 +48,6 @@ class SpikeDetectTestUi(UserTestUi):
     """UI for manually pass/failing spike detection unit tests.
     """
     def __init__(self):
-<<<<<<< HEAD
-        pg.mkQApp()
-
-        self.widget = pg.QtGui.QSplitter(pg.QtCore.Qt.Horizontal)
-        self.widget.resize(2000, 1000)
-        self.display1 = SpikeDetectUI('expected result')
-        self.display2 = SpikeDetectUI('current result')
-        self.widget.addWidget(self.display1.widget)
-        self.widget.addWidget(self.display2.widget)
-
-        self.ctrl = pg.QtGui.QWidget()
-        self.widget.addWidget(self.ctrl)
-        self.ctrl_layout = pg.QtGui.QVBoxLayout()
-        self.ctrl.setLayout(self.ctrl_layout)
-        self.diff_widget = pg.DiffTreeWidget()
-        self.ctrl_layout.addWidget(self.diff_widget)
-
-        self.pass_btn = pg.QtGui.QPushButton('pass')
-        self.fail_btn = pg.QtGui.QPushButton('fail')
-        self.ctrl_layout.addWidget(self.pass_btn)
-        self.ctrl_layout.addWidget(self.fail_btn)
-
-        self.pass_btn.clicked.connect(self.pass_clicked)
-        self.fail_btn.clicked.connect(self.fail_clicked)
-
-        self.last_btn_clicked = None
-        self.widget.setSizes([800, 800, 400])
-
-    def pass_clicked(self):
-        self.last_btn_clicked = 'pass'
-
-    def fail_clicked(self):
-        self.last_btn_clicked = 'fail'
-
-    def user_passfail(self):
-        self.widget.show()
-        while True:
-            print('in user_passfail')
-            pg.QtGui.QApplication.processEvents()
-            last_btn_clicked = self.last_btn_clicked
-            self.last_btn_clicked = None
-
-            if last_btn_clicked == 'fail':
-                raise Exception("User rejected test result.")
-            elif last_btn_clicked == 'pass':
-                break
-            
-            time.sleep(0.03)
-
-    def show_results(self, expected, current):
-        self.display1.show_result(expected)
-        self.display2.show_result(current)
-        self.diff_widget.setData(expected, current)
-
-    def clear(self):
-        self.display1.clear()
-        self.display2.clear()
-        self.diff_widget.setData(None, None)
-||||||| merged common ancestors
-        pg.mkQApp()
-
-        self.widget = pg.QtGui.QSplitter(pg.QtCore.Qt.Horizontal)
-        self.widget.resize(1600, 1000)
-
-        self.display1 = SpikeDetectUI('expected result')
-        self.display2 = SpikeDetectUI('current result')
-        self.widget.addWidget(self.display1.widget)
-        self.widget.addWidget(self.display2.widget)
-
-        self.ctrl = pg.QtGui.QWidget()
-        self.widget.addWidget(self.ctrl)
-        self.ctrl_layout = pg.QtGui.QVBoxLayout()
-        self.ctrl.setLayout(self.ctrl_layout)
-        self.diff_widget = pg.DiffTreeWidget()
-        self.ctrl_layout.addWidget(self.diff_widget)
-
-        self.pass_btn = pg.QtGui.QPushButton('pass')
-        self.fail_btn = pg.QtGui.QPushButton('fail')
-        self.ctrl_layout.addWidget(self.pass_btn)
-        self.ctrl_layout.addWidget(self.fail_btn)
-
-        self.pass_btn.clicked.connect(self.pass_clicked)
-        self.fail_btn.clicked.connect(self.fail_clicked)
-
-        self.last_btn_clicked = None
-        self.widget.setSizes([400, 400, 800])
-
-    def pass_clicked(self):
-        self.last_btn_clicked = 'pass'
-
-    def fail_clicked(self):
-        self.last_btn_clicked = 'fail'
-
-    def user_passfail(self):
-        self.widget.show()
-        while True:
-            pg.QtGui.QApplication.processEvents()
-            last_btn_clicked = self.last_btn_clicked
-            self.last_btn_clicked = None
-
-            if last_btn_clicked == 'fail':
-                raise Exception("User rejected test result.")
-            elif last_btn_clicked == 'pass':
-                break
-            
-            time.sleep(0.03)
-
-    def show_results(self, expected, current):
-        self.display1.show_result(expected)
-        self.display2.show_result(current)
-        self.diff_widget.setData(expected, current)
-
-    def clear(self):
-        self.display1.clear()
-        self.display2.clear()
-        self.diff_widget.setData(None, None)
-=======
         expected_display = SpikeDetectUI('expected result')
         current_display = SpikeDetectUI('current result')
         UserTestUi.__init__(self, expected_display, current_display)
->>>>>>> master

From 86f44f934fc37e5ab33cac0e6f697ea3b42120c0 Mon Sep 17 00:00:00 2001
From: Corinne <corinnet@alleninstitute.org>
Date: Fri, 9 Aug 2019 12:28:25 -0700
Subject: [PATCH 04/10] fix Trace to Tseries

---
 neuroanalysis/spike_detection.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py
index fbd9336..4900acb 100644
--- a/neuroanalysis/spike_detection.py
+++ b/neuroanalysis/spike_detection.py
@@ -109,7 +109,7 @@ def detect_ic_evoked_spikes(trace,
         max_slope: float
 
     """
-    if not isinstance(trace, Trace):
+    if not isinstance(trace, TSeries):
         raise TypeError("data must be Trace instance.")
 
     if ui is not None:

From 8f4c0e201b2d0581b0c84f7b6bd28a4fb5357ec4 Mon Sep 17 00:00:00 2001
From: Corinne <corinnet@alleninstitute.org>
Date: Thu, 15 Aug 2019 20:49:48 -0700
Subject: [PATCH 05/10] intermittent checkin

---
 neuroanalysis/event_detection.py |  69 +++++++--
 neuroanalysis/spike_detection.py | 238 +++++++++++++++++++++++--------
 2 files changed, 235 insertions(+), 72 deletions(-)

diff --git a/neuroanalysis/event_detection.py b/neuroanalysis/event_detection.py
index a7d827f..3c634f2 100644
--- a/neuroanalysis/event_detection.py
+++ b/neuroanalysis/event_detection.py
@@ -95,6 +95,26 @@ 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):
+    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, debug=False):
     """
@@ -174,22 +194,41 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
             mpl.axhline(y=-threshold, color='y', linestyle='--')
             mpl.show(block=False)
 
-        if len(on_inds) == 0 or len(off_inds) == 0:
-            continue
 
-        ## if there are unequal number of crossing from one direction, either remove the ends or add the appropriate initial or end index (which will be the beginning or end of the waveform)
-        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))
+
+        # sometimes an event happens at the beginning of the pulse window and the trace hasn't 
+        # been able to drop below threshold because it hasn't recovered from the artifact.
+        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))
+
+        # this is insufficient because it ignores when there are just off or just on
+        # not sure why this is here if haven't decided whether or not to omit ends
+        # if len(on_inds) == 0 or len(off_inds) == 0:
+        #     continue    
+        # ## if there are unequal number of crossing from one direction, either remove the ends or add the appropriate initial or end index (which will be the beginning or end of the waveform)
+        # 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))
         
         # Add both events above +threshold and those below -threshold to a list
         for i in range(len(on_inds)):
diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py
index 4900acb..a8598ee 100644
--- a/neuroanalysis/spike_detection.py
+++ b/neuroanalysis/spike_detection.py
@@ -55,9 +55,12 @@ def detect_ic_evoked_spikes(trace,
                             dv2_threshold=40.e3, 
                             mse_threshold=30., 
                             dv2_event_area=10e-6,
-                            pulse_term_bound=50.e-6,
+                            pulse_bounds_move=[.2e-3, 0.03e-3],  #0.03e-3 is very close to edge of artifact
                             double_spike=1.e-3,
-                            ui=None):
+                            ui=None,
+                            artifact_width=.1e-3,
+                            dv_threshold = 85. 
+                            ):
     """Return a dict describing an evoked spike in a patch clamp recording. Or Return None
     if no spike is detected.
 
@@ -82,16 +85,20 @@ def detect_ic_evoked_spikes(trace,
         has happened.
     event_area: float
         The integral of the 'bump' in dv2 must have at least the following *dv2_event_area*. 
-    pulse_term_boundry: float
+    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_term_boundry* specifies how much time before the termination 
-        of the stimulation pulse should not be considered in the pulse window.  Spikes after this
-        point are identified by attempting to fit the dv decay to the trace that would be seen 
-        from standard RC decay (see *mse_threshold*). 
+        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 between 
     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]
 
     Returns 
     =======
@@ -118,70 +125,108 @@ def detect_ic_evoked_spikes(trace,
         ui.plt1.plot(trace.time_values, trace.data)
 
     assert trace.data.ndim == 1
-    pulse_edges = tuple(map(float, pulse_edges))  # confirms 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()
+    diff1 = trace.time_slice(window_edges[0], window_edges[1]).diff()
     diff2 = diff1.diff()
 
-    # mask out pulse artifacts in diff2 before lowpass filtering
+    #---------------------------------------------------------
+    # -----not sure looking at dv/dt is the best thing----------
+    # # 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)
+        apply_cos_mask(diff2, center=edge + 200e-6, radius=400e-6, power=2)
 
     # low pass filter the second derivative
     diff2 = bessel_filter(diff2, 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))
+    # # # look for positive bumps in second derivative
+    events2 = list(threshold_events(diff2 / dv2_threshold, 
+        threshold=1.0, adjust_times=False, debug=False))
+    #---------------------------------------------------------
+
+    # scale the threshold by the height of the largest pulse
+    #TODO: not sure that this is relevant any more not that conditions at the 
+    #beginning and end of a trace are addressed.
+    # mip, is_edge = max_time(diff1.time_slice(pulse_edges[0]+.1e-3, pulse_edges[1]))
+    # max_dv_in_pulse = diff1.value_at(mip)
+    # dv_threshold = max_dv_in_pulse/3.
+
+    events1 = list(threshold_events(diff1, 
+                                    threshold=dv_threshold, 
+                                    adjust_times=False, 
+                                    omit_ends=False,
+                                    debug=False))
 
     if ui is not None:
         ui.plt2.plot(diff1.time_values, diff1.data)
+        ui.plt2.addLine(y=dv_threshold)
         ui.plt3.plot(diff2.time_values, diff2.data)
-        ui.plt3.addLine(y=dv2_threshold)
     
     # for each bump in d2vdt, either discard the event or generate 
     # spike metrics from v and dvdt  
     spikes = []
-    for ev in events2:
+
+    slope_hit_boundry = False
+    for ev in events1:
+        if ev['sum'] < 0:
+            continue
+
+
+        #TODO: what if it is wider than it is tall
         total_area = ev['area']
-        onset_time = ev['time']
+        onset_time = ev['time'] #this number is arbitrary 
 
         # ignore events near pulse offset
-        if abs(onset_time - pulse_edges[1]) < pulse_term_bound:
-            continue
+        # if abs(onset_time - pulse_edges[1]) < pulse_term_bound:
+        #     continue
 
         # require dv2 bump to be positive, not tiny
-        if total_area < dv2_event_area:
-            continue
+        # if total_area < dv2_event_area:
+        #     continue
         
         # don't double-count spikes 
         if len(spikes) > 0 and onset_time < spikes[-1]['onset_time'] + double_spike:
             continue
 
-        max_slope_window = onset_time, pulse_edges[1] - pulse_term_bound
+
+        max_slope_window = onset_time, window_edges[1] #TODO: don't like this should be an amount of time after onset not the end of the pulse search window
         max_slope_chunk = diff1.time_slice(*max_slope_window)
         if len(max_slope_chunk) == 0:
             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] - pulse_term_bound))
+        #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
+
+        #TODO: currently this is repetitive
         if max_slope <= 30:  # mV/ms
             continue
+        #this should only happen at the end of pulse window since we are looking at dvdt in events
         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] + pulse_term_bound:
-            # peak is obscured by pulse edge
-            peak_time = None
+            slope_hit_boundry = True
+            # max_slope_time = None
+            # max_slope = None
+            peak_time = None # slope cant be found either can peak, will look below
+        else: 
+            slope_hit_boundry = False
+            #TODO: the multiplicitive factor by ev['duration'] does not seem very principled
+            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,
@@ -191,10 +236,14 @@ def detect_ic_evoked_spikes(trace,
             '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()
+    # 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. there a potential spike was found but it appears to have straddled the end of the pulse
+    # TODO: not quite sure it is flushed out yet
+    if (len(spikes) == 0) or slope_hit_boundry or (spikes[-1]['max_slope_time'] < (window_edges[1] - double_spike)): #last spike is more than 1 ms before end
+
+        dv_after_pulse = trace.time_slice(window_edges[1] + artifact_width, None).diff() #note this is removing an area around the temination artifact
         dv_after_pulse = bessel_filter(dv_after_pulse, 15e3, bidir=True)
 
         # create a vector to fit
@@ -204,29 +253,77 @@ def detect_ic_evoked_spikes(trace,
         # 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:
+            
+            #TODO not sure these pulse edges make sense: should they be window edges
             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 = pulse_edges[1] + artifact_width
+            max_slope_time, slope_is_edge = max_time(diff.time_slice(onset_time, pulse_edges[1] + search_window))
+            max_slope = dv_after_pulse.value_at(max_slope_time)
+            peak_time, peak_is_edge = max_time(trace.time_slice(max_slope_time or window_edges[1] + artifact_width, window_edges[1] + search_window))
+            peak_value = dv_after_pulse.value_at(peak_time)
+            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:  
+                        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
+                })
+        
+        # if there is no spike found from the decay 
+        else:
+            # but if the last spike recorded was against the boundry
+            if slope_hit_boundry:  
+                # 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
 
@@ -288,7 +385,7 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None):
     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 = diff2.time_slice(None, diff2.t_end)
 
     # look for negative bumps in second derivative
     # dv1_threshold = 1e-6
@@ -320,7 +417,7 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None):
 
     spikes = []
     for ev in events:
-        if np.abs(ev['sum']) <10:
+        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
@@ -330,12 +427,20 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None):
             continue
 
         #TODO: What is this doing?
+        # if ev['sum'] < 0:
+        #     onset_time = ev['peak_time']
+        #     search_time = onset_time
+        # else:
+        #     search_time = ev['time'] - 200e-6
+        #     onset_time = ev['peak_time']  
+
         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 = ev['peak_time']  
+            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) #note this is looking for min because event must be negative in VC
@@ -358,10 +463,6 @@ def detect_vc_evoked_spikes(trace, pulse_edges, ui=None):
                 peak_time = None
             else:
                 peak = trace.time_at(peak_time)
-        #============================================================================
-        if debug:
-            import pdb; pdb.set_trace()
-        #============================================================================
 
         spikes.append({
             'onset_time': onset_time,
@@ -371,6 +472,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

From ff3579ce872cbc2cb0b49cfa9f25f23a4700e53a Mon Sep 17 00:00:00 2001
From: Corinne <corinnet@alleninstitute.org>
Date: Fri, 16 Aug 2019 08:45:06 -0700
Subject: [PATCH 06/10] working fairly well but a few more changes needed

---
 neuroanalysis/event_detection.py | 23 +++++++++++++++--------
 neuroanalysis/spike_detection.py | 32 +++++++++++++++++++++++---------
 2 files changed, 38 insertions(+), 17 deletions(-)

diff --git a/neuroanalysis/event_detection.py b/neuroanalysis/event_detection.py
index 3c634f2..d6d4299 100644
--- a/neuroanalysis/event_detection.py
+++ b/neuroanalysis/event_detection.py
@@ -127,9 +127,11 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
     Parameters                                                                                                                                     
     ==========
     trace: *Trace* instance                            
-    threshold: float
+    threshold: float or np.array with dimensions of *trace.data*
         algorithm checks if waveform crosses both positive and negative thresholds.  
-        i.e. if -5. is provided, the algorithm looks for places where the waveform crosses +/-5 
+        i.e. if -5. is provided, the algorithm looks for places where the waveform crosses +/-5.
+        If an array is provided the *threshold* is dynamic
+
     adjust_times: boolean
         if True, move the start and end times of the event outward, estimating the zero-crossing point for the event
     
@@ -158,10 +160,15 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
         peak_time: float, or np.nan if timing not available
            time of peak
     """
-    
-    threshold = abs(threshold)
+
+
     data = trace.data
     data1 = data - baseline
+    # convert threshold array
+    if isinstance(threshold, float):
+        threshold = np.ones(len(data)) * abs(threshold)
+
+
     #if (hasattr(data, 'implements') and data.implements('MetaArray')):
     
     ## find all positive and negative threshold crossings of baseline adjusted data
@@ -190,8 +197,8 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
                 mpl.axvline(x=on, color='g', linestyle='--')
             for off in off_inds:
                 mpl.axvline(x=off, color='r', linestyle='--')
-            mpl.axhline(y=threshold, color='y', linestyle='--')
-            mpl.axhline(y=-threshold, color='y', linestyle='--')
+            mpl.plot(threshold, color='y', linestyle='--')
+            mpl.plot(-threshold, color='y', linestyle='--')
             mpl.show(block=False)
 
 
@@ -254,8 +261,8 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
         for hit in hits:
             mpl.axvline(x=hit[0], color='g', linestyle='--')
             mpl.axvline(x=hit[1], color='r', linestyle='--')
-        mpl.axhline(y=threshold, color='y', linestyle='--')
-        mpl.axhline(y=-threshold, color='y', linestyle='--')
+        mpl.plot(threshold, color='y', linestyle='--')
+        mpl.plot(-threshold, color='y', linestyle='--')
         mpl.show(block=False)
     
     n_events = len(hits)
diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py
index a8598ee..1fd04e6 100644
--- a/neuroanalysis/spike_detection.py
+++ b/neuroanalysis/spike_detection.py
@@ -55,11 +55,11 @@ def detect_ic_evoked_spikes(trace,
                             dv2_threshold=40.e3, 
                             mse_threshold=30., 
                             dv2_event_area=10e-6,
-                            pulse_bounds_move=[.2e-3, 0.03e-3],  #0.03e-3 is very close to edge of artifact
+                            pulse_bounds_move=[.15e-3, 0.02e-3],  #0.03e-3 is very close to edge of artifact
                             double_spike=1.e-3,
                             ui=None,
                             artifact_width=.1e-3,
-                            dv_threshold = 85. 
+                            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.
@@ -99,6 +99,8 @@ def detect_ic_evoked_spikes(trace,
     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]
+    dv_thresholds: [float, float]
+        *dv_thresholds[0]* (*dv_thresholds[1]*) threshold for potential spike at beginning (end) of pulse window.
 
     Returns 
     =======
@@ -143,7 +145,9 @@ def detect_ic_evoked_spikes(trace,
     for edge in pulse_edges:
         apply_cos_mask(diff2, center=edge + 200e-6, radius=400e-6, power=2)
 
+
     # low pass filter the second derivative
+    diff1 = bessel_filter(diff1, 10e3, order=4, bidir=True)
     diff2 = bessel_filter(diff2, 10e3, order=4, bidir=True)
 
     # # # look for positive bumps in second derivative
@@ -158,6 +162,10 @@ def detect_ic_evoked_spikes(trace,
     # max_dv_in_pulse = diff1.value_at(mip)
     # dv_threshold = max_dv_in_pulse/3.
 
+    thresh_start=125
+    thresh_end=15
+
+    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, 
@@ -166,7 +174,7 @@ def detect_ic_evoked_spikes(trace,
 
     if ui is not None:
         ui.plt2.plot(diff1.time_values, diff1.data)
-        ui.plt2.addLine(y=dv_threshold)
+        ui.plt2.plot(diff1.time_values, dv_threshold, pen='y')
         ui.plt3.plot(diff2.time_values, diff2.data)
     
     # for each bump in d2vdt, either discard the event or generate 
@@ -210,12 +218,15 @@ def detect_ic_evoked_spikes(trace,
         # require dv/dt to be above a threshold value
 
         #TODO: currently this is repetitive
-        if max_slope <= 30:  # mV/ms
+        # if max_slope <= 30:  # mV/ms
+        #     continue
+        if is_edge ==-1:
+            #if it is an edge at the beginning V is most likely just increasing from pulse initiation  
             continue
         #this should only happen at the end of pulse window since we are looking at dvdt in events
-        if is_edge != 0:
-            # can't see max slope
+        if is_edge == 1:
             slope_hit_boundry = True
+            # can't see max slope
             # max_slope_time = None
             # max_slope = None
             peak_time = None # slope cant be found either can peak, will look below
@@ -236,12 +247,14 @@ def detect_ic_evoked_spikes(trace,
             'max_slope': max_slope,
         })
 
+
     # 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. there a potential spike was found but it appears to have straddled the end of the pulse
+    #3. there 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 it is too fast to see it in decay so use current values
     # TODO: not quite sure it is flushed out yet
-    if (len(spikes) == 0) or slope_hit_boundry or (spikes[-1]['max_slope_time'] < (window_edges[1] - double_spike)): #last spike is more than 1 ms before end
+    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
 
         dv_after_pulse = trace.time_slice(window_edges[1] + artifact_width, None).diff() #note this is removing an area around the temination artifact
         dv_after_pulse = bessel_filter(dv_after_pulse, 15e3, bidir=True)
@@ -284,12 +297,13 @@ def detect_ic_evoked_spikes(trace,
                 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
+                    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: 

From 3fb09733dcd7e2fb975582319e832db4099e0991 Mon Sep 17 00:00:00 2001
From: Corinne Teeter <corinnet@alleninstitute.org>
Date: Fri, 16 Aug 2019 15:02:39 -0700
Subject: [PATCH 07/10] cleaning code but untested

---
 neuroanalysis/spike_detection.py | 116 ++++++++++---------------------
 1 file changed, 36 insertions(+), 80 deletions(-)

diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py
index 1fd04e6..ebacff6 100644
--- a/neuroanalysis/spike_detection.py
+++ b/neuroanalysis/spike_detection.py
@@ -55,10 +55,10 @@ def detect_ic_evoked_spikes(trace,
                             dv2_threshold=40.e3, 
                             mse_threshold=30., 
                             dv2_event_area=10e-6,
-                            pulse_bounds_move=[.15e-3, 0.02e-3],  #0.03e-3 is very close to edge of artifact
+                            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,
+                            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
@@ -66,20 +66,17 @@ def detect_ic_evoked_spikes(trace,
 
     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.  Currently, if a spike is detected in the region 
-    shortly after the stimulation pulse it is identified as a spike but the initiation index/time
-    is not resolved.
+    short region after the stimulation pulse.  Note that finding 
+
 
     Parameters
     ==========
-    trace: Trace instance
+    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*.
-    dv2_threshold: float
-        Value for the second derivative of the voltage must cross to be considered as a possible spike.  
-    mse_threshold: float
+     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.
@@ -137,34 +134,12 @@ def detect_ic_evoked_spikes(trace,
 
     # calculate derivatives within pulse window
     diff1 = trace.time_slice(window_edges[0], window_edges[1]).diff()
-    diff2 = diff1.diff()
-
-    #---------------------------------------------------------
-    # -----not sure looking at dv/dt is the best thing----------
-    # # mask out pulse artifacts in diff2 before lowpass filtering
-    for edge in pulse_edges:
-        apply_cos_mask(diff2, center=edge + 200e-6, radius=400e-6, power=2)
-
+    diff2 = diff1.diff() # note that this isnt actually used for anything other than plotting
 
     # low pass filter the second derivative
     diff1 = bessel_filter(diff1, 10e3, order=4, bidir=True)
-    diff2 = bessel_filter(diff2, 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, debug=False))
-    #---------------------------------------------------------
-
-    # scale the threshold by the height of the largest pulse
-    #TODO: not sure that this is relevant any more not that conditions at the 
-    #beginning and end of a trace are addressed.
-    # mip, is_edge = max_time(diff1.time_slice(pulse_edges[0]+.1e-3, pulse_edges[1]))
-    # max_dv_in_pulse = diff1.value_at(mip)
-    # dv_threshold = max_dv_in_pulse/3.
-
-    thresh_start=125
-    thresh_end=15
 
+    # 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, 
@@ -177,7 +152,7 @@ def detect_ic_evoked_spikes(trace,
         ui.plt2.plot(diff1.time_values, dv_threshold, pen='y')
         ui.plt3.plot(diff2.time_values, diff2.data)
     
-    # for each bump in d2vdt, either discard the event or generate 
+    # for each bump in dvdt, either discard the event or generate 
     # spike metrics from v and dvdt  
     spikes = []
 
@@ -186,55 +161,33 @@ def detect_ic_evoked_spikes(trace,
         if ev['sum'] < 0:
             continue
 
-
-        #TODO: what if it is wider than it is tall
-        total_area = ev['area']
         onset_time = ev['time'] #this number is arbitrary 
 
-        # ignore events near pulse offset
-        # if abs(onset_time - pulse_edges[1]) < pulse_term_bound:
-        #     continue
-
-        # require dv2 bump to be positive, not tiny
-        # if total_area < dv2_event_area:
-        #     continue
-        
         # don't double-count spikes 
         if len(spikes) > 0 and onset_time < spikes[-1]['onset_time'] + double_spike:
             continue
 
-
-        max_slope_window = onset_time, window_edges[1] #TODO: don't like this should be an amount of time after onset not the end of the pulse search window
-        max_slope_chunk = diff1.time_slice(*max_slope_window)
-        if len(max_slope_chunk) == 0:
-            continue
-        max_slope_idx = np.argmax(max_slope_chunk.data)
-        max_slope_time = max_slope_chunk.time_at(max_slope_idx)
+#        max_slope_window = onset_time, window_edges[1] 
+#        max_slope_chunk = diff1.time_slice(*max_slope_window)
+#        if len(max_slope_chunk) == 0:
+#            continue
+#        max_slope_idx = np.argmax(max_slope_chunk.data)
+#        max_slope_time = max_slope_chunk.time_at(max_slope_idx)
 
         #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
-
-        #TODO: currently this is repetitive
-        # if max_slope <= 30:  # mV/ms
-        #     continue
         if is_edge ==-1:
-            #if it is an edge at the beginning V is most likely just increasing from pulse initiation  
+            #if it is an edge at the beginning of the pulse  V is most likely just increasing from pulse initiation  
             continue
         #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
-            # can't see max slope
-            # max_slope_time = None
-            # max_slope = None
-            peak_time = None # slope cant be found either can peak, will look below
+            peak_time = None # if slope is hitting the pulse window boundry don't trust the peak
         else: 
             slope_hit_boundry = False
-            #TODO: the multiplicitive factor by ev['duration'] does not seem very principled
             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
@@ -248,16 +201,20 @@ def detect_ic_evoked_spikes(trace,
         })
 
 
-    # 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. there 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 it is too fast to see it in decay so use current values
-    # TODO: not quite sure it is flushed out yet
-    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
+    """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
 
-        dv_after_pulse = trace.time_slice(window_edges[1] + artifact_width, None).diff() #note this is removing an area around the temination artifact
-        dv_after_pulse = bessel_filter(dv_after_pulse, 15e3, bidir=True)
+        #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) #and then filtering it again
 
         # 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
@@ -279,11 +236,11 @@ def detect_ic_evoked_spikes(trace,
         if mse > mse_threshold:
             
             #TODO not sure these pulse edges make sense: should they be window edges
-            search_window = 2e-3
-            onset_time = pulse_edges[1] + artifact_width
-            max_slope_time, slope_is_edge = max_time(diff.time_slice(onset_time, pulse_edges[1] + search_window))
+            search_window = 2e-3 #this is arbitrary but should be long enough to detect and max slopes and peak from a spike
+            onset_time = window_edges[1] + artifact_width
+            max_slope_time, slope_is_edge = max_time(dv_after_pulse.time_slice(onset_time, onset_time + search_window))
             max_slope = dv_after_pulse.value_at(max_slope_time)
-            peak_time, peak_is_edge = max_time(trace.time_slice(max_slope_time or window_edges[1] + artifact_width, window_edges[1] + search_window))
+            peak_time, peak_is_edge = max_time(dv_after_pulse.time_slice(max_slope_time or window_edges[1] + artifact_width, window_edges[1] + search_window))
             peak_value = dv_after_pulse.value_at(peak_time)
             if peak_is_edge != 0:
                 peak_time = None
@@ -328,11 +285,10 @@ def detect_ic_evoked_spikes(trace,
                 'max_slope': max_slope
                 })
         
-        # if there is no spike found from the decay 
         else:
-            # but if the last spike recorded was against the boundry
+            # there is no spike found from the decay 
             if slope_hit_boundry:  
-                # delete it
+                # if last spike recorded was against the bound delete it
                 spikes.pop(-1)
 
 

From d04c747ba26c2a355c39a4c9772af67f635d4a8f Mon Sep 17 00:00:00 2001
From: Corinne <corinnet@alleninstitute.org>
Date: Fri, 16 Aug 2019 23:30:53 -0700
Subject: [PATCH 08/10] works has reduced unknowns to unclear examples

---
 neuroanalysis/spike_detection.py | 54 +++++++++++++++++---------------
 1 file changed, 28 insertions(+), 26 deletions(-)

diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py
index ebacff6..6393f15 100644
--- a/neuroanalysis/spike_detection.py
+++ b/neuroanalysis/spike_detection.py
@@ -52,7 +52,6 @@ def rc_decay(t, tau, Vo):
 
 def detect_ic_evoked_spikes(trace, 
                             pulse_edges, 
-                            dv2_threshold=40.e3, 
                             mse_threshold=30., 
                             dv2_event_area=10e-6,
                             pulse_bounds_move=[.15e-3, 0.02e-3],  #0.03e-3 sometimes gets some artifact
@@ -134,9 +133,17 @@ def detect_ic_evoked_spikes(trace,
 
     # calculate derivatives within pulse window
     diff1 = trace.time_slice(window_edges[0], window_edges[1]).diff()
-    diff2 = diff1.diff() # note that this isnt actually used for anything other than plotting
 
+    # 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)
 
     # create a linear threshold
@@ -167,13 +174,6 @@ def detect_ic_evoked_spikes(trace,
         if len(spikes) > 0 and onset_time < spikes[-1]['onset_time'] + double_spike:
             continue
 
-#        max_slope_window = onset_time, window_edges[1] 
-#        max_slope_chunk = diff1.time_slice(*max_slope_window)
-#        if len(max_slope_chunk) == 0:
-#            continue
-#        max_slope_idx = np.argmax(max_slope_chunk.data)
-#        max_slope_time = max_slope_chunk.time_at(max_slope_idx)
-
         #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)
@@ -210,11 +210,11 @@ def detect_ic_evoked_spikes(trace,
        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
+    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) #and then filtering it again
+        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
@@ -234,14 +234,12 @@ def detect_ic_evoked_spikes(trace,
 
         # there is a spike, so identify paramters
         if mse > mse_threshold:
-            
-            #TODO not sure these pulse edges make sense: should they be window edges
-            search_window = 2e-3 #this is arbitrary but should be long enough to detect and max slopes and peak from a spike
             onset_time = window_edges[1] + artifact_width
-            max_slope_time, slope_is_edge = max_time(dv_after_pulse.time_slice(onset_time, onset_time + search_window))
+            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(dv_after_pulse.time_slice(max_slope_time or window_edges[1] + artifact_width, window_edges[1] + search_window))
+            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
@@ -298,7 +296,10 @@ def detect_ic_evoked_spikes(trace,
     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
+                            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
@@ -312,6 +313,8 @@ 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*.
+    ui:                                                                                                                            
+        user interface for viewing spike detection 
     
     Returns 
     =======
@@ -337,23 +340,22 @@ 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)
-
-    #---------------------------------------------------------
-    #----this is were vc and ic code diverge------------------
-    #---------------------------------------------------------
-
+    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(pulse_edges[0], pulse_edges[1] + 2e-3).diff()
+    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)
+    # 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)
+    diff2 = bessel_filter(diff2, 10e3, order=4, bidir=True)
     # chop off ending transient
     diff2 = diff2.time_slice(None, diff2.t_end)
 

From 3d945b801f2e7501000afc3aa1d5378d58ca302c Mon Sep 17 00:00:00 2001
From: Corinne <corinnet@alleninstitute.org>
Date: Sat, 17 Aug 2019 11:38:14 -0700
Subject: [PATCH 09/10] current version

---
 neuroanalysis/spike_detection.py | 60 +++++++++++---------------------
 1 file changed, 21 insertions(+), 39 deletions(-)

diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py
index 6393f15..ae7e2ae 100644
--- a/neuroanalysis/spike_detection.py
+++ b/neuroanalysis/spike_detection.py
@@ -75,7 +75,7 @@ def detect_ic_evoked_spikes(trace,
         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
+    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.
@@ -89,12 +89,13 @@ def detect_ic_evoked_spikes(trace,
         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 between 
+        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]
+        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.
 
@@ -104,15 +105,18 @@ def detect_ic_evoked_spikes(trace,
         contains following information about spikes:
         onset_time: float
             Onset time of region where searching for spike. Defined as a crossing 
-            of a *threshold* or *baseline* in dv2/dt.
-        peak_time: float
+            of a *threshold* in dv2/dt.
+        peak_time: float or None
             time of spike peak
-        max_slope_time:
-            time where the voltage is changing most rapidly (i.e. dvdt = 0) 
-        peak_value: float
-            None if peak_time is None else trace.value_at(peak_time),
-        max_slope: float
-
+        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.")
@@ -281,7 +285,7 @@ def detect_ic_evoked_spikes(trace,
                 'peak_time': peak_time,
                 'peak_value':peak_value,
                 'max_slope': max_slope
-                })
+            })
         
         else:
             # there is no spike found from the decay 
@@ -299,6 +303,7 @@ def detect_ic_evoked_spikes(trace,
 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.
 
@@ -313,6 +318,9 @@ def detect_vc_evoked_spikes(trace,
         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 
     
@@ -352,41 +360,24 @@ def detect_vc_evoked_spikes(trace,
 
     # crop and low pass filter the second derivative
     diff2 = diff2.time_slice(pulse_edges[0] + 150e-6, pulse_edges[1])
-    # 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)
     diff2 = bessel_filter(diff2, 10e3, order=4, bidir=True)
-    # chop off ending transient
-    diff2 = diff2.time_slice(None, diff2.t_end)
 
     # look for negative bumps in second derivative
-    # dv1_threshold = 1e-6
-    dv2_threshold = 0.02
-
-    #=========================================================================
-    debug = False
-    #=========================================================================
-
     events = list(threshold_events(diff2 / dv2_threshold, 
-        threshold=1.0, adjust_times=False, omit_ends=True, debug=debug))
+                                   threshold=1.0, adjust_times=False, omit_ends=True, debug=False))
 
 
 
     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)
-        # ui.plt3.plot(diff2.time_values, diff2.data/dv2_threshold)
-        # ui.plt3.addLine(y=1)
 
     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:
         if np.abs(ev['sum']) < 2.:
@@ -397,15 +388,6 @@ def detect_vc_evoked_spikes(trace,
         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
-
-        #TODO: What is this doing?
-        # if ev['sum'] < 0:
-        #     onset_time = ev['peak_time']
-        #     search_time = onset_time
-        # else:
-        #     search_time = ev['time'] - 200e-6
-        #     onset_time = ev['peak_time']  
-
         if ev['sum'] < 0:
             # only accept positive bumps
             continue

From a9724048ba4ca291946c658da6a6af1c69bf8689 Mon Sep 17 00:00:00 2001
From: Corinne <corinnet@alleninstitute.org>
Date: Wed, 21 Aug 2019 18:28:25 -0700
Subject: [PATCH 10/10] moving event size out of event detection and fixing
 index out of range bug

---
 neuroanalysis/event_detection.py | 171 +++++++++++--------------------
 neuroanalysis/spike_detection.py |  14 ++-
 2 files changed, 67 insertions(+), 118 deletions(-)

diff --git a/neuroanalysis/event_detection.py b/neuroanalysis/event_detection.py
index d6d4299..620d01e 100644
--- a/neuroanalysis/event_detection.py
+++ b/neuroanalysis/event_detection.py
@@ -95,7 +95,7 @@ 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):
+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
@@ -108,7 +108,13 @@ def deal_unbalanced_initial_off(omit_ends, on_inds, off_inds):
         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):
+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:
@@ -116,35 +122,40 @@ def deal_unbalanced_termination_on(omit_ends, on_inds, off_inds, off_to_add):
     return on_inds, off_inds
 
 
-def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_ends=True, debug=False):
+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) and then 
-    recross threshold (bumps).  If a threshold is crossed at the end of the trace, an event may be excluded
+    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*). 
-    Optionally adjusts start and end index of an event to an extrapolated baseline-crossing and calculates 
-    values. 
+
 
     Parameters                                                                                                                                     
     ==========
-    trace: *Trace* instance                            
+    trace: *Tseries* instance                            
     threshold: float or np.array with dimensions of *trace.data*
-        algorithm checks if waveform crosses both positive and negative thresholds.  
-        i.e. if -5. is provided, the algorithm looks for places where the waveform crosses +/-5.
-        If an array is provided the *threshold* is dynamic
-
+        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
+        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 is a region of the *Trace.data* waveform that crosses above *threshold* and then falls below threshold again
-        Sometimes referred to as a 'bump'.  There are additional criteria (not listed here) for a bump to be considered an event.    
-        Each index contains information about an event.  Fields as follows:
+        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
+            index length of the event 
         sum: float
            sum of the values in the array between the start and end of the event
         peak: float
@@ -164,106 +175,53 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
 
     data = trace.data
     data1 = data - baseline
+    
     # convert threshold array
     if isinstance(threshold, float):
         threshold = np.ones(len(data)) * abs(threshold)
-
-
-    #if (hasattr(data, 'implements') and data.implements('MetaArray')):
     
-    ## find all positive and negative threshold crossings of baseline adjusted data
-
-    if debug:
-        import pdb; pdb.set_trace()
-        # FYI: can't use matplot lib before debugger is on
-        # type *continue to see plot
-        import matplotlib.pyplot as mpl
-
-
-    masks = [(data1 > threshold).astype(np.byte), (data1 < -threshold).astype(np.byte)] # 1 where data is [above threshold, below negative 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] # -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
+        diff = mask[1:] - mask[:-1]
+        # 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) 
         
-        # TODO: It might be a good idea to make offindexes one less so that region above threshold looks symmetrical
-        # find start and end inicies (above threshold) where waveform is above threshold
-        on_inds = list(np.argwhere(diff==1)[:,0] + 1) #where crosses from below to above threshold Note taking index after this happens
-        off_inds = list(np.argwhere(diff==-1)[:,0]) #where crosses from above to below threshold. Note taking index before this happens
-
-        if debug:
-            mpl.figure()
-            mpl.plot(data1, '.-') 
-            for on in on_inds:
-                mpl.axvline(x=on, color='g', linestyle='--')
-            for off in off_inds:
-                mpl.axvline(x=off, color='r', linestyle='--')
-            mpl.plot(threshold, color='y', linestyle='--')
-            mpl.plot(-threshold, color='y', linestyle='--')
-            mpl.show(block=False)
-
-
-
-        # sometimes an event happens at the beginning of the pulse window and the trace hasn't 
-        # been able to drop below threshold because it hasn't recovered from the artifact.
+        # 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)
+                    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)
+                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))
+                    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))
+                on_inds, off_inds = _deal_unbalanced_termination_on(omit_ends, on_inds, off_inds, len(data1))
 
-        # this is insufficient because it ignores when there are just off or just on
-        # not sure why this is here if haven't decided whether or not to omit ends
-        # if len(on_inds) == 0 or len(off_inds) == 0:
-        #     continue    
-        # ## if there are unequal number of crossing from one direction, either remove the ends or add the appropriate initial or end index (which will be the beginning or end of the waveform)
-        # 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))
-        
-        # Add both events above +threshold and those below -threshold to a list
+
+        # 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)):
-            
-            # remove any point where an on off is seperated by less than 2 indicies
-            if (on_inds[i] + 1) == off_inds[i]: 
-                continue
-            if (on_inds[i] - 1) == off_inds[i]: 
+            if on_inds[i] == off_inds[i]:
+                #something wierd happened
                 continue
-            if on_inds[i] == off_inds[i]: 
-                continue
-            
             hits.append((on_inds[i], off_inds[i]))
     
     ## sort hits  ## NOTE: this can be sped up since we already know how to interleave the events..
     hits.sort(key=lambda a: a[0])
-    if debug is True:
-        # FYI: can't use matplot lib before debugger is on
-        # type *continue to see plot
-        mpl.figure()
-        mpl.title('first round before adjustment')
-        mpl.plot(data1, '.-') 
-        for hit in hits:
-            mpl.axvline(x=hit[0], color='g', linestyle='--')
-            mpl.axvline(x=hit[1], color='r', linestyle='--')
-        mpl.plot(threshold, color='y', linestyle='--')
-        mpl.plot(-threshold, color='y', linestyle='--')
-        mpl.show(block=False)
     
     n_events = len(hits)
     events = np.empty(n_events, dtype=[
@@ -284,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:
@@ -292,18 +250,18 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
         else:
             peak_ind = np.argmin(ev_data)
         peak = ev_data[peak_ind]
-        peak_ind += ind1 # adjust peak_ind from local event data to entire waveform
+        peak_ind += ind1
             
         #print "event %f: %d" % (xvals[ind1], ind1) 
         if adjust_times:  ## Move start and end times outward, estimating the zero-crossing point for the event
         
             ## adjust ind1 first
-            mind = np.argmax(ev_data)  # max of whole trace
-            pdiff = abs(peak - ev_data[0])  # find the how high the peak is from the front of event 
-            if pdiff == 0:  
+            mind = np.argmax(ev_data)
+            pdiff = abs(peak - ev_data[0])
+            if pdiff == 0:
                 adj1 = 0
             else:
-                adj1 = int(threshold * mind / pdiff)  # (max value of whole trace)* 1/(hight of peak from first data point)
+                adj1 = int(threshold * mind / pdiff)
                 adj1 = min(ln, adj1)
             ind1 -= adj1
             
@@ -363,18 +321,6 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
     
         ## remove masked events
         events = events[mask]
-        
-        if debug:
-            mpl.figure()
-            mpl.title('adjusted')
-            mpl.plot(data1, '.-') 
-            for on in events['index']:
-                mpl.axvline(x=on, color='g', linestyle='--')
-            # for off in off_inds:
-            #     mpl.axvline(x=off, color='r', linestyle='--')
-            mpl.axhline(y=threshold, color='y', linestyle='--')
-            mpl.axhline(y=-threshold, color='y', linestyle='--')
-            mpl.show(block=False)
 
     # add in timing information if available:
     if trace.has_timing:
@@ -391,9 +337,6 @@ def threshold_events(trace, threshold, adjust_times=True, baseline=0.0, omit_end
         ev['area'] = np.nan
         ev['peak_time'] = np.nan
 
-
-    if debug:
-        pdb.set_trace()
     return events
 
 
diff --git a/neuroanalysis/spike_detection.py b/neuroanalysis/spike_detection.py
index ae7e2ae..4a2121e 100644
--- a/neuroanalysis/spike_detection.py
+++ b/neuroanalysis/spike_detection.py
@@ -155,8 +155,7 @@ def detect_ic_evoked_spikes(trace,
     events1 = list(threshold_events(diff1, 
                                     threshold=dv_threshold, 
                                     adjust_times=False, 
-                                    omit_ends=False,
-                                    debug=False))
+                                    omit_ends=False))
 
     if ui is not None:
         ui.plt2.plot(diff1.time_values, diff1.data)
@@ -169,6 +168,10 @@ def detect_ic_evoked_spikes(trace,
 
     slope_hit_boundry = False
     for ev in events1:
+        # require 3 indexes are above threshold
+        if ev['len'] < 3:
+            continue
+        # require events to be positive
         if ev['sum'] < 0:
             continue
 
@@ -364,8 +367,7 @@ def detect_vc_evoked_spikes(trace,
 
     # look for negative bumps in second derivative
     events = list(threshold_events(diff2 / dv2_threshold, 
-                                   threshold=1.0, adjust_times=False, omit_ends=True, debug=False))
-
+                    threshold=1.0, adjust_times=False, omit_ends=True))
 
 
     if ui is not None:
@@ -380,6 +382,10 @@ def detect_vc_evoked_spikes(trace,
     # 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: