<jupyter_start><jupyter_text>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.![Quadrature encoding schematic](http://www.ni.com/cms/images/devzone/tut/ipuuzhqc3503.jpg)For more information on the rotary encoder see these links:  [National Instruments guide to rotary encoders](http://www.ni.com/tutorial/7109/en/)  [Datasheet for the Kübler](https://www.kuebler.com/pdf?2400-2420_en.pdf)  UnitsThe wheel module contains some default values that are useful for determining the units of the raw data, and for interconverting your ALF position units.<jupyter_code>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)<jupyter_output>The wheel diameter is 6.2 cm and the number of ticks is 4096 per revolution<jupyter_text>Loading the wheel dataThe 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](https://docs.google.com/document/d/1OqIqqakPakHXRAwceYLwFY9gOrm8_P62XIfCTnHwstg/editheading=h.hnjqyfnroyya).  **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.<jupyter_code>wheel = one.load_object(eid, 'wheel')

print('wheel.position: \n', wheel.position)
print('wheel.timestamps: \n', wheel.timestamps)<jupyter_output>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]<jupyter_text>Loading movement dataThe wheel movement onsets and offsets can be calculated using the `wheel.movements` function.  The output of this function will soon be extracted as the 'wheelMoves' ALF (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.<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>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:<jupyter_code># 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)
# 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)<jupyter_output><empty_output><jupyter_text>Although currently not extracted, there will soon be an ALF dataset type called 'wheelMoves'.  These will have 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.The peak velocity times will not be saved into an ALF file.You can extract the wheel moves with a single function.  Below we attempt to load the wheelMoves ALF and upon failing, extract it ourselves.<jupyter_code>try:
    # Warning: This object is currently not extracted (as of 2020-06-01)
    wheel_moves = one.load_object(eid, 'wheelMoves')
    assert wheel_moves, 'object not found'
except AssertionError:
    wheel_moves = extract_wheel_moves(wheel.timestamps, wheel.position)<jupyter_output>2020-06-19 16:48:52.874 WARNING  [one.py:267] eafbdb1a-8776-4390-b210-76b7509e31d0 does not contain any wheelMoves object datasets
2020-06-19 16:48:52.898 INFO     [training_wheel.py:266] Wheel in rad units using X4 encoding<jupyter_text>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](https://docs.google.com/document/d/1RA6wgbWfxD2kGlpNxt0n3HVcW4TEIx8e-YO7k_W1pHs/edit))<jupyter_code>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)<jupyter_output>The wheel must be turned ~0.3 rad to move the stimulus to threshold<jupyter_text>Calculating velocityWheel 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.<jupyter_code># 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)<jupyter_output><empty_output><jupyter_text>Last move onsetThe `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.<jupyter_code>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);<jupyter_output><empty_output><jupyter_text>Calculating reaction timesReaction 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.<jupyter_code># Load the reaction times
# brainbox.io.one.load_wheel_reaction_times
rt = load_wheel_reaction_times(eid)
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);<jupyter_output><empty_output><jupyter_text>Finding 'determined' movementsBelow 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.<jupyter_code>firstMove_times, is_final_movement, ids = FirstMovementTimes.extract_first_movement_times(wheel_moves, trial_data)<jupyter_output>2020-06-19 16:50:31.975 INFO     [training_trials.py:706] minimum quiescent period assumed to be 200ms<jupyter_text>Direction changesBelow 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 you're goal is.<jupyter_code>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');<jupyter_output><empty_output><jupyter_text>The function `direction_changes` does the same as above, returning a list of times and indices of each movement's direction changes.<jupyter_code>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');<jupyter_output><empty_output><jupyter_text>Splitting wheel trace by trialTo 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.<jupyter_code>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();<jupyter_output><empty_output><jupyter_text>DataJoint SchemaThere is also a share DJ table that has been populated using the same functions as above.<jupyter_code>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')<jupyter_output>Connecting <EMAIL>:3306
Connected to https://alyx.internationalbrainlab.org as miles<jupyter_text>WheelMoveSetThis table contains some movement information at the session level.<jupyter_code>dj_wheel.WheelMoveSet.describe();<jupyter_output># 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<jupyter_text>WheelMoveSet.MoveThis table contains each detected movement for a session.<jupyter_code>dj_wheel.WheelMoveSet.Move.describe();<jupyter_output>-> 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<jupyter_text>MovementTimesThis table contains movement information at the level of the trial.  Trials may not have an associated movement if the goCue or stimOn times were missing, or if no movements were detected for that trial.  `movement_onset` is the only field that is in absolute seconds.<jupyter_code>dj_wheel.MovementTimes.describe();
# 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
# 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)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0881807340665165, 'fraction_numerical': 0.023651782165098716, 'mean_word_length': 5.021140805743917, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 69, 'https://': 4, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '218365', 'license': 'MIT', 'n_tokens_mistral': 4899, 'n_tokens_neox': 4497, 'n_words': 1641, 'path': 'brainbox/examples/wheel_moves.ipynb', 'repo_name': 'k1o0/ibllib', 'size': '15134'}