Bergmann_2022_DataRepository/04_analysis/saccades.py

import numpy as _np
from scipy.signal import find_peaks as _find_peaks
import sliding1d as _sliding

def get_sign(x):
    return x / _np.abs(x)

def sliding_diff(x):
    x   = _np.array(x)
    ret = _np.empty(x.size, dtype=_np.float64)
    ret[1:-1] = (x[2:] - x[:-2])/2
    ret[0]    = x[1] - x[0]
    ret[-1]   = x[-1] - x[-2]
    return ret

def smoothing_diff(x, rad=5, num=3):
    return sliding_diff(_sliding.nanmean(x, rad, num))

def threshold_values(x, threshold):
    ret = _np.array(x, copy=True)
    ret[_np.abs(x) < threshold] = 0
    return ret

def detect(time, left_pupil, right_pupil,
           std_period=None, std_threshold=5, min_distance_seconds=0.2,
           smoothing_radius_seconds=0.05, smoothing_number=3):
    """detect saccade events.
    
    Parameters
    ----------
    
    time        -- the zero-starting trial time array.
    left_pupil  -- the left-pupil position at each time point.
    right_pupil -- the right-pupil position at each time point.
    
    For the keyword arguments, refer to the "Procedures" section below.
    
    Returns
    -------
    
    a numpy.ndarray object whose size equals to the size of the `time` parameter.
    the returned array has all-zero values, except for the positions of the detected saccades.
    the values at the time points of detected saccades represent the average velocity of the two eyes.
    
    Procedures
    ----------
    
    1. smoothing_diff() is used to compute velocity from position.
       - `smoothing_radius_seconds` and `smoothing_number` is used here.
       
    2. the velocity is thresholded by its abosolute value
       with the criterion: `abs(v) > v[std_period].std(ddof=1)*std_threshold`
       - if `std_period` is None, then the initial 1 second of the recording is used.
       
    3. "average velocity" is calculated as the geometric mean 
       of the (thresholded) velocities of the left and the right pupils.
       
    4. `scipy.signal.find_peaks` is used to detect peaks from the average velocity.
       - `min_distance_seconds` is used to set the minimum distance between the peaks.
     
    5. the returned array is formatted. the value of each event is computed based on:
       - amplitude: the amplitude in the average velocity
       - signature: the signature in the left-pupil velocity.
       
    """
    dt = _np.diff(time).mean()
    if std_period is None:
        std_period = time < 1 # the initial 1 second
    smoothing_radius     = int(round(smoothing_radius_seconds / dt))
    min_distance_samples = int(round(min_distance_seconds / dt))
    
    vleft  = smoothing_diff(left_pupil, rad=smoothing_radius, num=smoothing_number) / dt
    vright = smoothing_diff(right_pupil, rad=smoothing_radius, num=smoothing_number) / dt
    vleft_thresholded  = threshold_values(vleft, vleft[std_period].std(ddof=1)*std_threshold)
    vright_thresholded = threshold_values(vright, vright[std_period].std(ddof=1)*std_threshold)
    
    # there are some "negative spikes", possibly due to small unsynchronization of pupil motions
    # we take both negative and positive spikes here
    vamp  = _np.sqrt(_np.abs(vleft_thresholded*vright_thresholded))
    peaks = _find_peaks(vamp, distance=min_distance_samples)[0]
    
    ret   = _np.zeros(left_pupil.size, dtype=_np.float64)
    ret[peaks] = vamp[peaks] * get_sign(vleft_thresholded[peaks])
    return ret

def annotate(events):
    """splits the saccade events array into three traces for the plotting purpose.
    
    Parameter
    ---------
    
    events -- as it is returned from saccades.detect()
    
    Returns
    -------
    
    a `dict` object, which contains the following keys and values:
    
    - none:  no-event period trace, for plotting "baseline" period.
    - left:  leftward saccade event trace
    - right: rightward saccade event trace
    
    For event traces, NaNs are used to represent values at irrelevant time points.
    
    """
    noevents    = _np.zeros(events.size, dtype=_np.float64)
    leftevents  = _np.empty(events.size, dtype=_np.float64); leftevents[:]  = _np.nan
    rightevents = _np.empty(events.size, dtype=_np.float64); rightevents[:] = _np.nan
    for levt in _np.where(events > 0)[0]:
        rng = slice(levt-1, levt+2)
        leftevents[rng] = events[rng]
    for revt in _np.where(events < 0)[0]:
        rng = slice(revt-1, revt+2)
        rightevents[rng] = events[rng]
    return dict(none=noevents, left=leftevents, right=rightevents)