Merge pull request #284 from jpgill86/amplitude-discriminator-type

Add amplitude discriminator parameter for extrema type (peak or trough), and epoch annotation for SpikeTrain

docs/metadata.rst

@@ -484,16 +484,24 @@ provides bandpass filtering.
 Amplitude Discriminators
 ------------------------
 
-Spikes with peaks that fall within amplitude windows given by
+Spikes with peaks (or troughs) that fall within amplitude windows given by
 ``amplitude_discriminators`` can be automatically detected by *neurotic* on the
-basis of amplitude alone. Note that amplitude discriminators are only applied
-if fast loading is off (``lazy=False``).
+basis of amplitude. Note that amplitude discriminators are only applied if fast
+loading is off (``lazy=False``).
 
 Detected spikes are indicated on the signals with markers, and spike trains are
 displayed in a raster plot. Optionally, a color may be specified for an
 amplitude discriminator using a single letter color code (e.g., ``'b'`` for
 blue or ``'k'`` for black) or a hexadecimal color code (e.g., ``'1b9e77'``).
 
+The algorithm can detect either peaks or troughs in the signal. When both the
+lower and upper bounds for amplitude windows are positive, the default behavior
+is to detect peaks. When both are negative, the default is to detect troughs.
+These defaults can be overridden using `type: trough` or `type: peak`,
+respectively. This is useful when, for example, detecting subthreshold
+excitatory postsynaptic potentials. If the signs of the bounds differ, explicit
+specification of the type is required.
+
 In addition to restricting spike detection for a given unit to an amplitude
 window, detection can also be limited in time to overlap with epochs with a
 given label.
@@ -522,17 +530,25 @@ for each amplitude discriminator.
               epoch: Unit 2 activity
               color: 'e6ab02'
 
-Here two units are detected on the same channel with different amplitude
-windows. Any peaks between 50 and 150 microvolts on the "Extracellular" channel
-will be tagged as a spike belonging to "Unit 1". The discriminator for "Unit 2"
-provides the optional ``epoch`` parameter. This restricts detection of "Unit 2"
-to spikes within the amplitude window that occur at the same time as epochs
-labeled "Unit 2 activity". These epochs can be created by the epoch encoder
-(reload required to rerun spike detection at launch-time), specified in the
-read-only ``annotations_file``, or even be contained in the ``data_file`` if
-the format supports epochs.
-
-Amplitude windows are permitted to be negative.
+            - name: Unit 3
+              channel: Intracellular
+              units: mV
+              amplitude: [-10, 60]
+              type: peak
+
+Here two units are detected on the "Extracellular" channel with different
+amplitude windows, and a third unit is detected on the "Intracellular" channel.
+On the "Extracellular" channel, any peaks between 50 and 150 microvolts will be
+tagged as a spike belonging to "Unit 1". The discriminator for "Unit 2" detects
+smaller peaks, between 20 and 50 microvolts, and it provides the optional
+``epoch`` parameter. This restricts detection of "Unit 2" to spikes within the
+amplitude window that occur at the same time as epochs labeled "Unit 2
+activity". These epochs can be created by the epoch encoder (reload required to
+rerun spike detection at launch-time), specified in the read-only
+``annotations_file``, or even be contained in the ``data_file`` if the format
+supports epochs. Finally, peaks between -10 and +60 millivolts will be detected
+on the "Intracellular" channel; because the signs of these bounds differ, the
+type (peak or trough) must be explicitly given.
 
 .. _config-metadata-tridesclous:
 

neurotic/datasets/data.py

@@ -518,13 +518,16 @@ def _create_neo_spike_trains_from_dataframe(dataframe, metadata, t_start, t_stop
             st = neo.SpikeTrain(
                 name = str(spike_label),
                 file_origin = _abs_path(metadata, 'tridesclous_file'),
-                channels = channels, # custom annotation
-                amplitude = None,    # custom annotation
                 times = t_start + sampling_period * df['index'].values,
                 t_start = t_start,
                 t_stop = t_stop,
             )
 
+            st.annotate(
+                channels=channels,
+                amplitude=None,
+            )
+
             spiketrain_list.append(st)
 
     return spiketrain_list
@@ -601,31 +604,46 @@ def _detect_spikes(sig, discriminator, epochs):
 
     min_threshold = min(discriminator['amplitude'])
     max_threshold = max(discriminator['amplitude'])
-    if min_threshold >= 0 and max_threshold > 0:
+    spike_type = discriminator.get('type', None)
+    if spike_type == 'peak':
         sign = 'above'
-    elif min_threshold < 0 and max_threshold <= 0:
+    elif spike_type == 'trough':
         sign = 'below'
+    elif spike_type is None:
+        # infer type from thresholds
+        if min_threshold >= 0 and max_threshold > 0:
+            spike_type = 'peak'
+            sign = 'above'
+        elif min_threshold < 0 and max_threshold <= 0:
+            spike_type = 'trough'
+            sign = 'below'
+        else:
+            raise ValueError('automatic spike type inference for amplitude discriminator is possible only with two nonnegative thresholds (type=peak) or two nonpositive thresholds (type=trough); otherwise, type must be given explicitly: {}'.format(discriminator))
     else:
-        raise ValueError('amplitude discriminator must have two nonnegative thresholds or two nonpositive thresholds: {}'.format(discriminator))
+        raise ValueError('amplitude discriminator type must be "peak", "trough", or unspecified: {}'.format(discriminator))
 
     spikes_crossing_min = _elephant_tools.peak_detection(sig, pq.Quantity(min_threshold, discriminator['units']), sign, 'raw')
     spikes_crossing_max = _elephant_tools.peak_detection(sig, pq.Quantity(max_threshold, discriminator['units']), sign, 'raw')
-    if sign == 'above':
+    if spike_type == 'peak':
         spikes_between_min_and_max = np.setdiff1d(spikes_crossing_min, spikes_crossing_max)
-    elif sign == 'below':
+    elif spike_type == 'trough':
         spikes_between_min_and_max = np.setdiff1d(spikes_crossing_max, spikes_crossing_min)
     else:
-        raise ValueError('sign should be "above" or "below": {}'.format(sign))
+        raise ValueError('type should be "peak" or "trough": {}'.format(spike_type))
 
     st = neo.SpikeTrain(
         name = discriminator['name'],
-        channels = [discriminator['channel']],  # custom annotation
-        amplitude = pq.Quantity(discriminator['amplitude'], discriminator['units']), # custom annotation
         times = spikes_between_min_and_max * pq.s,
         t_start = sig.t_start,
         t_stop  = sig.t_stop,
     )
 
+    st.annotate(
+        channels=[discriminator['channel']],
+        amplitude=pq.Quantity(discriminator['amplitude'], discriminator['units']),
+        type=spike_type,
+    )
+
     if 'epoch' in discriminator:
 
         time_masks = []
@@ -649,6 +667,8 @@ def _detect_spikes(sig, discriminator, epochs):
         # windows
         st = st[np.any(time_masks, axis=0)]
 
+        st.annotate(epoch=discriminator['epoch'])
+
     return st
 
 def _run_burst_detectors(metadata, blk):