Working with wheel data

This example will give you information about how the rotary encoder records wheel movements, how to find out the units and spatial resolution of the standard wheel data, and how to load the wheel data. There are also examples of how to load wheel movements and rection times from ALF datasets and DataJoint tables, as well as how to calculate these from scratch.

[1]:
%matplotlib notebook

import re
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

from brainbox.io.one import load_wheel_reaction_times
import brainbox.behavior.wheel as wh
from ibllib.io.extractors.ephys_fpga import extract_wheel_moves
from ibllib.io.extractors.training_wheel import extract_first_movement_times
# from ibllib.misc.exp_ref import eid2ref
from oneibl.one import ONE

one = ONE()
sns.set_style('whitegrid')
Connected to https://alyx.internationalbrainlab.org as mayo
[2]:
# NB: This function will soon be available from ibllib.misc.exp_ref
def eid2ref(eid):
    """
    Get human-readable session ref from path
    :param eid: The experiment uuid to find reference for
    :return: dict containing 'subject', 'date' and 'sequence'
    """
    path_str = str(one.path_from_eid(eid))
    pattern = r'(?P<subject>[\w-]+)([\\/])(?P<date>\d{4}-\d{2}-\d{2})(\2)(?P<sequence>\d{3})'
    match = re.search(pattern, path_str)
    return match.groupdict()
[3]:
eid = 'eafbdb1a-8776-4390-b210-76b7509e31d0'
eid2ref(eid)
Out[3]:
{'subject': 'CSHL055', 'date': '2020-02-10', 'sequence': '001'}

Device

The standard resolution of the rotary encoder is 1024 ‘ticks’ per revolulation. With quadrature (X4) encoding we measure four fronts giving from the device: low -> high and high -> low, from channels A and B. Therefore the number of measured ‘ticks’ per revolution becomes 4096.

X4 encoding is used in ephys sessions, while in training sessions the encoding is X1, i.e. four times fewer ticks per revolution.

Quadrature encoding schematic

For more information on the rotary encoder see these links:

Units

The wheel module contains some default values that are useful for determining the units of the raw data, and for interconverting your ALF position units.

[4]:
device_info = ('The wheel diameter is {} cm and the number of ticks is {} per revolution'
               .format(wh.WHEEL_DIAMETER, wh.ENC_RES))
print(device_info)
The wheel diameter is 6.2 cm and the number of ticks is 4096 per revolution

Loading the wheel data

The wheel traces can be accessed via ONE. There are two ALF files: wheel.position and wheel.timestamps. The timestamps are, as usual, in seconds (same as the trials object) and the positions are in radians. The positions are not evenly sampled, meaning that when the wheel doesn’t move, no position is registered. For reference,

More information on the ALF dataset types can be found here.

NB: There are some old wheel.velocity ALF files that contain the volocity measured as the diff between neighbouring samples. Later on we describe a better way to calculate the velocity.

[5]:
wheel = one.load_object(eid, 'wheel')

print('wheel.position: \n', wheel.position)
print('wheel.timestamps: \n', wheel.timestamps)
wheel.position:
 [-3.86563159e-01 -3.92699082e-01 -3.98835005e-01 ... -5.65799610e+02
 -5.65805746e+02 -5.65811882e+02]
wheel.timestamps:
 [6.79604098e-01 6.86161141e-01 6.92652183e-01 ... 3.73136665e+03
 3.73138577e+03 3.73141459e+03]

Loading movement data

There is also an ALF dataset type called ‘wheelMoves’ with two attributes:
1. wheelMoves.intervals - An N-by-2 arrays where N is the number of movements detected. The first column contains the movement onset time in absolute seconds; the second contains the offset. 2. wheelMoves.peakAmplitude - The absolute maximum amplitude of each detected wheel movement, relative to onset position. This can be used to determine whether the movement was particularly large and therefore whether it was a flinch vs a determined movement.

If the dataset doesn’t exist you can also extract the wheel moves with a single function. Below we attempt to load the wheelMoves ALF and upon failing, extract it ourselves.

[6]:
try:
    # Warning: Some older sessions may not have a wheelMoves dataset
    wheel_moves = one.load_object(eid, 'wheelMoves')
    assert wheel_moves, 'object not found'
except AssertionError:
    wheel_moves = extract_wheel_moves(wheel.timestamps, wheel.position)

The movements algorithm

The wheel movement onsets and offsets are calculated using the wheel.movements function. The output of this function is saved in the ‘wheelMoves’ ALF.

Refer to the wheel.movements function docstring for details of how the movements are detected. In addition to the data found in the wheelMoves object, the function outputs an array of peak velocity times. Also the function has a make_plots flag which will produce plots of the wheel position and velocity with the movement onsets and offsets highlighted (see below).

The movements algorithm requires the positions and timestamps to be evenly sampled so they should be interpolated first, which can be done with the wheel.interpolate_position function. The default sampling frequency is 1000Hz:

[7]:
pos, t = wh.interpolate_position(wheel.timestamps, wheel.position)
sec = 5  # Number of seconds to plot
plt.figure()

# Plot the interpolated data points
mask = t < (t[0] + sec)
plt.plot(t[mask], pos[mask], '.', markeredgecolor='lightgrey', markersize=1)

# Plot the original data
mask = wheel.timestamps < (wheel.timestamps[0] + sec)
plt.plot(wheel.timestamps[mask], wheel.position[mask], 'r+', markersize=6)

# Labels etc.
plt.xlabel('time / sec')
plt.ylabel('position / rad')
plt.box(on=None)

Once interpolated, the movements can be extracted from the position trace. NB: The position thresholds are dependant on the wheel position units. The defaults values are for the raw input of a X4 1024 encoder, so we will convert them to radians:

[8]:
# Convert the pos threshold defaults from samples to correct unit
thresholds_cm = wh.samples_to_cm(np.array([8, 1.5]), resolution=wh.ENC_RES)
thresholds = wh.cm_to_rad(thresholds_cm)
[9]:
# Detect wheel movements for the first 5 seconds
mask = t < (t[0] + sec)

onsets, offsets, peak_amp, peak_vel_times = wh.movements(
    t[mask], pos[mask], pos_thresh=thresholds[0], pos_thresh_onset=thresholds[0], make_plots=True)

For scale, the stimulus must be moved 35 visual degrees to reach threshold. The wheel gain is 4°/mm (NB: the gain is double for the first session or so, see Appendix 2 of the behavior paper)

[10]:
threshold_deg = 35 # visual degrees
gain = 4  # deg / mm
threshold_rad = wh.cm_to_rad(1e-1) * (threshold_deg / gain)  # rad

print('The wheel must be turned ~%.1f rad to move the stimulus to threshold' % threshold_rad)
The wheel must be turned ~0.3 rad to move the stimulus to threshold

Calculating velocity

Wheel velocity can be calculated using the velocity_smoothed function, which returns the velocity and acceleration convolved with a Gaussian window. As with the movements function, the input is expected to be evenly sampled, therefore you should interpolate the wheel data before calling this function. The default window size of 3ms is reasonable, and interpolating at a frequency of 1000 (the default) is sufficiently high.

[11]:
# pos was the output of interpolate_position using the default frequency of 1000Hz
Fs = 1000
pos, t = wh.interpolate_position(wheel.timestamps, wheel.position, freq=Fs)
vel, acc = wh.velocity_smoothed(pos, Fs)

Last move onset

The movements algorithm is the recommended way of detecting movement onsets because it is quicker and more accurate, however there is another function that will return the last movement onset before a particular event. This is useful for quickly finding the movement that reached threshold for a given trial. This function finds the first sample after the velocity has been zero for at least 50ms. Because it uses velocity, the smoothed derivative of position, it is less accurate. NB: The more accurate approach is to find all moves for which the onset occured before feedback time and the offset occured afterwards.

[12]:
trial_data = one.load_object(eid, 'trials')
idx = 23 # trial index
ts = wh.last_movement_onset(t, vel, trial_data['feedback_times'][idx]);

mask = np.logical_and(trial_data['goCue_times'][idx] < t, t < trial_data['feedback_times'][idx])
plt.figure();
plt.plot(t[mask], pos[mask]);
plt.axvline(x=ts);

Calculating reaction times

Reaction times based on wheel movements can be calculated with the load_wheel_reaction_times function which is located in the behavior module of brainbox.

Reaction times are defined as the time between the go cue (onset tone) and the onset of the first substantial wheel movement. A movement is considered sufficiently large if its peak amplitude is at least 1/3rd of the distance to threshold (~0.1 radians).

Negative times mean the onset of the movement occurred before the go cue. Nans may occur if there was no detected movement withing the period, or when the goCue_times or feedback_times are nan.

The function loads the trials object and if firstMovement_times is not present it loads the wheel moves and extracts these times.

[13]:
# Load the reaction times
# brainbox.io.one.load_wheel_reaction_times
rt = load_wheel_reaction_times(eid)
Connected to https://alyx.internationalbrainlab.org as mayo
[14]:
trial_data = one.load_object(eid, 'trials')

# # Replace nans with zeros
trial_data.contrastRight[np.isnan(trial_data.contrastRight)] = 0
trial_data.contrastLeft[np.isnan(trial_data.contrastLeft)] = 0

contrast = trial_data.contrastRight - trial_data.contrastLeft
mean_rt = [np.nanmean(rt[contrast == c]) for c in set(contrast)]

# RT may be nan if there were no detected movements, or if the goCue or stimOn times were nan
xdata = np.unique(contrast)
plt.figure(figsize=(4, 3))  # Some sort of strange behaviour in this cell's output
plt.plot(xdata, mean_rt);

plt.xlabel('contrast')
plt.ylabel('mean rt / s')
plt.ylim(bottom=0);

Finding reaction time and ‘determined’ movements

Below is an example of how you might filter trials by those responses that are unambiguous. The function extract_first_movement_times is used for calculating the trial reaction times but also returns whether the first significant movement was ‘final’, i.e. was the one that reached threshold. For details of how ‘first movement’ is defined, see the function docstring.

[15]:
firstMove_times, is_final_movement, ids = extract_first_movement_times(wheel_moves, trial_data)
2020-10-29 14:54:32.974 INFO     [training_wheel.py:337] minimum quiescent period assumed to be 200ms

Direction changes

Below is an example of how to plot the times that the wheel changed direction. Changing the smoothing window when calculating the velocity may improve the detected changes, depending on what your goal is.

[16]:
n = 569
on, off = wheel_moves['intervals'][n,]
mask = np.logical_and(t > on, t < off)
sng = np.sign(vel[mask])
idx, = np.where(np.diff(sng) != 0)

plt.figure()
plt.plot(t[mask], pos[mask], 'k')
for i in idx:
    plt.axvline(x=t[mask][i], color='k', linestyle=':')

plt.title('Movement #%s' % n)
plt.xlabel('time / s')
plt.ylabel('position / rad');

The function direction_changes does the same as above, returning a list of times and indices of each movement’s direction changes.

[17]:
n = 403  # trial number
start, end = trial_data['intervals'][n,]  # trial intervals
intervals = wheel_moves['intervals']  # movement onsets and offsets

# Find direction changes for a given trial
mask = np.logical_and(intervals[:,0] > start, intervals[:,0] < end)
change_times, idx, = wh.direction_changes(t, vel, intervals[mask])

plt.figure()
mask = np.logical_and(t > start, t < end)  # trial intervals mask
plt.plot(t[mask], pos[mask], 'k')  # plot wheel trace for trial
for i in np.concatenate(change_times):
    plt.axvline(x=i, color='k', linestyle=':')

plt.title('Trial #%s' % n)
plt.xlabel('time / s')
plt.ylabel('position / rad');

Splitting wheel trace by trial

To plot a selection of ‘determined’ movements, we can split the traces using the traces_by_trial function.

NB: This using the within_ranges function which is generic and can be used to detect which points are within a range. This is useful for returning masks for slicing (must be cast to bool), to label points within ranges or to dectect whether any points belong to more than one range.

[18]:
n_trials = 3  # Number of trials to plot
# Randomly select the trials to plot
trial_ids = np.random.randint(trial_data['choice'].size, size=n_trials)
fig, axs = plt.subplots(1, n_trials, figsize=(8.5,2.5))
plt.tight_layout()

# Plot go cue and response times
goCues = trial_data['goCue_times'][trial_ids]
responses = trial_data['response_times'][trial_ids]

# Plot traces between trial intervals
starts = trial_data['intervals'][trial_ids, 0]
ends = trial_data['intervals'][trial_ids, 1]
# Cut up the wheel vectors
traces = wh.traces_by_trial(t, pos, start=starts, end=ends)
zipped = zip(traces, axs, goCues, responses, trial_ids)

for (trace, ax, go, resp, n) in zipped:
    ax.plot(trace[0], trace[1], 'k-')
    ax.axvline(x=go, color='g', label='go cue', linestyle=':')
    ax.axvline(x=resp, color='r', label='threshold', linestyle=':')
    ax.set_title('Trial #%s' % n)

    # Turn off tick labels
    ax.set_yticklabels([])
    ax.set_xticklabels([])

# Add labels to first
axs[0].set_xlabel('time / sec')
axs[0].set_ylabel('position / rad')
plt.legend();
plt.tight_layout()

There is also a shared DJ table that has been populated using the same functions as above. Note that trials.firstMovement_times is strictly calculated using the goCue_times and feedback_times, whereas the DJ tables use the minimum of stimOn and goCue, feedback and response. This means that if for example the goCue wasn’t detected for a given trial, the firstMovement time will still be present in the table.

[19]:
import datajoint as dj
from ibl_pipeline import acquisition, behavior
from uuid import UUID

dj_wheel = dj.create_virtual_module('wheel_moves', 'group_shared_wheel')
Connecting mayofaulkner@datajoint.internationalbrainlab.org:3306
Connected to https://alyx.internationalbrainlab.org as mayo

WheelMoveSet

This table contains some movement information at the session level.

[20]:
dj_wheel.WheelMoveSet.describe();
# Wheel movements occurring within a session
-> acquisition.Session
---
n_movements          : int                          # total number of movements within the session
total_displacement   : float                        # total displacement of the wheel during session in radians
total_distance       : float                        # total movement of the wheel in radians
n_direction_changes  : int                          # total number of direction changes within a session

WheelMoveSet.Move

This table contains each detected movement for a session.

[21]:
dj_wheel.WheelMoveSet.Move.describe();
-> dj_wheel.WheelMoveSet
move_id              : int                          # movement id
---
movement_onset       : float                        # time of movement onset in seconds from session start
movement_offset      : float                        # time of movement offset in seconds from session start
max_velocity         : float                        # time of movement's peak velocity
movement_amplitude   : float                        # the absolute peak amplitude relative to onset position

MovementTimes

This table contains movement information at the level of the trial. Trials may not have an associated movement if the goCue, stimOn, response and feedback times were missing, or if no movements were detected for that trial. movement_onset is the only field that is in absolute seconds.

[22]:
dj_wheel.MovementTimes.describe();
# Trial movements table
-> behavior.TrialSet.Trial
---
-> dj_wheel.WheelMoveSet.Move
reaction_time        : float                        # time in seconds from go cue to first sufficiently large movement onset of the trial
final_movement       : tinyint                      # indicates whether movement onset was the one that reached threshold
movement_time        : float                        # time in seconds from first movement onset to feedback time
response_time        : float                        # time in seconds from go cue to feedback time
movement_onset       : float                        # time in seconds when first movement onset occurred

[23]:
# Load the trials and trial movements for the same session we looked at above
session = acquisition.Session & {'session_uuid': UUID(eid)}
(dj_wheel.MovementTimes * behavior.TrialSet.Trial) & session
Out[23]:

subject_uuid

session_start_time

start time

trial_id

trial identification number

move_id

movement id

reaction_time

time in seconds from go cue to first sufficiently large movement onset of the trial

final_movement

indicates whether movement onset was the one that reached threshold

movement_time

time in seconds from first movement onset to feedback time

response_time

time in seconds from go cue to feedback time

movement_onset

time in seconds when first movement onset occurred

trial_start_time

beginning of quiescent period time (seconds)

trial_end_time

end of iti (seconds)

trial_response_time

trial_response_choice

which choice was made in choiceworld

trial_stim_on_time

Time of stimulus in choiceworld (seconds)

trial_stim_contrast_left

contrast of the stimulus on the left

trial_stim_contrast_right

contrast of the stimulus on the right

trial_go_cue_time

trial_go_cue_trigger_time

trial_feedback_time

trial_feedback_type

whether feedback is positive or negative in choiceworld (-1 for negative, +1 for positive)

trial_rep_num

the repetition number of the trial, i.e. how many trials have been repeated on this side (counting from 1)

trial_stim_prob_left

probability of the stimulus being present on left

trial_reward_volume

trial_iti_duration

trial_included

whether the trial should be included
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 1 107 7.958 0 4.491 12.449 308.928 0.0 316.019402 313.419 CW nan 0.0 0.25 300.97 300.97 313.4201 -1 None 0.5 0.0 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 2 111 0.2082 1 0.538 0.7462 318.761 316.8516 320.868602 319.299 CW 318.55280000000005 0.125 0.0 318.554 318.553 319.29900000000004 1 None 0.5 1.5 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 3 113 12.2207 1 0.268 12.4887 333.99 321.2422 335.818502 334.258 CW 321.76930000000004 0.0 0.0 321.77 321.769 334.2585 1 None 0.5 1.5 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 4 115 0.1908 1 0.12 0.3108 336.96 336.1835 338.652502 337.08 CCW 336.76919999999996 0.0 0.0625 336.77 336.769 337.08009999999996 1 None 0.5 1.5 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 5 117 0.1445 1 0.1109 0.2554 339.946 339.0159 341.618802 340.057 CCW 339.8015 0.0 0.0625 339.802 339.802 340.0569 1 None 0.5 1.5 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 6 119 0.1744 1 0.1639 0.3383 343.176 341.9859 344.902002 343.34 CW 343.0016 0.25 0.0 343.003 343.002 343.3399 1 None 0.5 1.5 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 7 122 0.127 1 0.097 0.224 346.145 345.3023 347.801402 346.242 CCW 346.01800000000003 0.0 0.25 346.019 346.018 346.242 1 None 0.5 1.5 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 8 123 0.129 1 0.126 0.255 348.781 348.1644 350.468702 348.907 CCW 348.652 0.0 0.0 348.653 348.652 348.9071 1 None 0.5 1.5 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 9 124 0.0845 1 0.108 0.1925 352.136 350.8114 353.817802 352.244 CCW 352.0515 0.0 1.0 352.052 352.052 352.2445 1 None 0.5 1.5 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 10 131 0.0091 1 0.259 0.2681 368.094 367.4547 369.917602 368.353 CW 368.0849 0.0 0.0 368.086 368.085 368.35360000000003 1 None 0.5 1.5 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 11 133 0.1853 1 0.079 0.2643 370.97 370.2891 373.618002 371.049 CCW 370.78470000000004 0.0625 0.0 370.786 370.785 371.0502 -1 None 0.5 0.0 nan None
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 12 134 0.0759 1 0.115 0.1909 375.077 373.9837 377.751502 375.192 CW 375.0011 0.0 0.0 375.002 375.001 375.1928 -1 None 0.5 0.0 nan None

...

Total: 735

[24]:
# Load reaction times and use NaN for trials with no detected movements
session = acquisition.Session & {'session_uuid': UUID(eid)}
(behavior.TrialSet.Trial & session).aggr(dj_wheel.MovementTimes, rt='reaction_time', keep_all_rows=True)
Out[24]:

subject_uuid

session_start_time

start time

trial_id

trial identification number

rt

time in seconds from go cue to first sufficiently large movement onset of the trial
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 1 7.958
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 2 0.2082
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 3 12.2207
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 4 0.1908
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 5 0.1445
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 6 0.1744
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 7 0.127
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 8 0.129
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 9 0.0845
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 10 0.0091
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 11 0.1853
3614ee37-c892-4eb4-ab83-f6c1f3804977 2020-02-10 09:50:44 12 0.0759

...

Total: 739