In [561]:
import numpy as np
import pandas as pd
import load_raw_data as lrd
from transforms3d import euler
from utils import nearest
import matplotlib.pyplot as plt
import seaborn as sns;sns.set(color_codes=True)
%matplotlib notebook
In [496]:
def parse_dates(timestamp):
"""
parse_dates(timestampstring) takes a timestamp string formatted as Year-month-dayThour:min:sec.decimals+Timezone
and converts it into a datetime.datetime format, ignoring the timezone and the last decimal, keeping a microsecond
precision.
"""
return pd.datetime.strptime(timestamp[:26], '%Y-%m-%dT%H:%M:%S.%f')
In [498]:
filepath = './Data/Pilot/Jesse_FirstPilot/head.csv'
hd = lrd.load_head_data(filepath)
In [499]:
hd.head()
Out[499]:
In [500]:
# determine the frequency of the signal
# first, get the minimum difference between timestamps:
np.diff(hd.index.values).min()
# Then, the frequency is 1 over that period
freq = 1e9 / np.diff(hd.index.values).min().astype(int)
In [501]:
np.diff(hd.index.values).min()
Out[501]:
In [502]:
ed = lrd.load_event_data('./Data/Pilot/Jesse_FirstPilot/events.csv')
ed = ed.drop(['Value1','Value2'],axis=1)
In [503]:
ed.head()
Out[503]:
In [504]:
trialstart_times = ed[ed['Name']=='DecisionStart'].index
In [505]:
ser = hd.loc[nearest(hd.index,trialstart_times[0])]
In [506]:
def series2mat4(head_pos):
return np.array([[head_pos.loc['Value.M11'],head_pos.loc['Value.M21'],head_pos.loc['Value.M31'],head_pos.loc['Value.M41']],
[head_pos.loc['Value.M12'],head_pos.loc['Value.M22'],head_pos.loc['Value.M32'],head_pos.loc['Value.M42']],
[head_pos.loc['Value.M13'],head_pos.loc['Value.M23'],head_pos.loc['Value.M33'],head_pos.loc['Value.M43']],
[head_pos.loc['Value.M14'],head_pos.loc['Value.M24'],head_pos.loc['Value.M34'],head_pos.loc['Value.M44']]])
In [531]:
head_position = series2mat4(ser)
In [532]:
## IMPORTANT: Watch out for gimbal lock.
euler_angles = euler.mat2euler(head_position)
np.degrees(euler_angles)
Out[532]:
In [533]:
euler_angles = euler.mat2euler(head_position,'syzx')
In [534]:
np.degrees(euler_angles)
Out[534]:
Let's now find the end point location for the first trial
In [535]:
hitneutral_times = ed[ed['Name']=='Neutral'].index
In [536]:
ser = hd.loc[nearest(hd.index,hitneutral_times[0])]
In [537]:
head_position = series2mat4(ser)
In [538]:
euler_angles = euler.mat2euler(head_position)
np.degrees(euler_angles)
Out[538]:
In [539]:
# now get all of them in between
In [540]:
trial1_trajectory = hd.loc[trialstart_times[0]:hitneutral_times[0]]
In [541]:
list_of_matrices = [series2mat4(hd.iloc[x]) for x in range(trial1_trajectory.shape[0])]
In [542]:
angles = np.array([np.degrees(euler.mat2euler(i,'syxz')) for i in list_of_matrices])
In [543]:
plt.plot(angles[:,1])
plt.title('Rotation trajectory around Y axis in first trial')
plt.show()
In [544]:
targetleft_times = ed[ed['Name']=='TargetLeft'].index
targetright_times = ed[ed['Name']=='TargetRight'].index
In [648]:
trial_numbers = np.argsort(targetleft_times.append(targetright_times))
In [671]:
trial_numbers
Out[671]:
In [649]:
# get the indides (iloc in dataframe) of the end of each trial left and right
end_trial_indices_left = [ed.index.get_loc(trial)+1 for trial in targetleft_times]
end_trial_indices_right = [ed.index.get_loc(trial)+1 for trial in targetright_times]
In [650]:
# and now get the corresponding timestamps
end_trial_times_left = ed.iloc[end_trial_indices_left].index
end_trial_times_right = ed.iloc[end_trial_indices_right].index
In [ ]:
In [694]:
# let's do this differently. All at once, and then determine left and right after
start_trial_times = targetleft_times.append(targetright_times).sort_values()
end_trial_times = end_trial_times_left.append(end_trial_times_right).sort_values()
In [695]:
# here, extract the list of left-right
target_sides = ed[ed.Name.str.get(0).isin(['T'])].reset_index()
In [696]:
trajectories = []
counter = 0
# Left trials
for i, (start, end) in enumerate(zip(start_trial_times,end_trial_times)):
trial_trajectory = hd.loc[start:end]
trial_trajectory = trial_trajectory.resample('0.01S').pad()
trial_trajectory.loc[:,'Trial number'] = i
trial_trajectory.loc[:,'Target side'] = target_sides.iloc[i]['Name']
trial_trajectory['Trial time'] = trial_trajectory.index - trial_trajectory.index[0]
trajectories.append(trial_trajectory)
trajectories_df = pd.concat(trajectories).sort_index()
In [852]:
# convert to matrices and then to angles
list_of_matrices = [series2mat4(trajectories_df.iloc[x]) for x in range(trajectories_df.shape[0])]
angles = np.array([np.degrees(euler.mat2euler(mat,'syzx')) for mat in list_of_matrices])
In [853]:
angles_df = pd.DataFrame(angles,index=trajectories_df.index,columns =['Y rotation','X rotation','Z rotation'])
In [854]:
trajectories_df = trajectories_df.join(angles_df)
In [855]:
trial_starts = trajectories_df[trajectories_df['Trial time']==trajectories_df.iloc[1]['Trial time']]
In [701]:
zeropoint = trial_starts['Y rotation'].mean()
trajectories_df['Y angle'] = trajectories_df['Y rotation'] - zeropoint
In [ ]:
In [797]:
fig = plt.figure()
ax = sns.tsplot(data=trajectories_df, time="Trial time", value='Y angle', unit='Trial number',condition='Target side')
plt.title('Rotation Trajectory')
plt.xlabel('Time (seconds)')
#plt.savefig('./Figures/rotation_trajectory.png')
plt.show()
In [720]:
# TODO: fix trial numbering system so this works
RT=[]
for i in trajectories_df['Trial number'].unique():
idx = trajectories_df['Trial number'] == i
RT.append(trajectories_df[idx]['Trial time'].max())
In [749]:
trials = pd.DataFrame(index=trajectories_df['Trial number'].unique(),
columns=['RT'],
data=np.array(np.array(RT)))
trials.index.name = 'Trial'
In [751]:
# add the target side info to this dataframe
trials['Target side'] = target_sides['Name']
In [778]:
trials['Reaction time (ms)']= trials['RT'].apply(lambda x: x.microseconds/1000)
In [780]:
sns.distplot(trials['Reaction time (ms)'],rug=True)
Out[780]:
In [790]:
# plot left and right separately
sns.distplot(trials.loc[trials['Target side']=='TargetRight','Reaction time (ms)'],
kde_kws={'label':'TargetRight'})
sns.distplot(trials.loc[trials['Target side']=='TargetLeft','Reaction time (ms)'],
kde_kws={'label':'TargetLeft'})
plt.title('Reaction time histograms with kernel density estimates')
# in order to fit a normal distribution instead: >>> from scipy.stats import norm, then fit=norm as argument
Out[790]:
In [795]:
from scipy.stats import norm
# plot left and right separately
sns.distplot(trials.loc[trials['Target side']=='TargetRight','Reaction time (ms)'],
kde=False,
fit=norm,
fit_kws={'color':'b','label':'TargetRight'},
label='TargetRight')
sns.distplot(trials.loc[trials['Target side']=='TargetLeft','Reaction time (ms)'],
kde=False,
fit=norm,
fit_kws={'color':'g','label':'TargetLeft'})
plt.title('Reaction time histograms with Gaussian fit')
Out[795]:
In [809]:
startpoints=[]
endpoints=[]
for i in trajectories_df['Trial number'].unique():
idx = trajectories_df['Trial number'] == i
startpoints.append(trajectories_df[idx].iloc[1]['Y angle'] )
endpoints.append(trajectories_df[idx].iloc[-1]['Y angle'] )
In [810]:
trials['Starting points'] = startpoints
In [811]:
trials['Movement endpoints'] = endpoints
In [813]:
trials.head()
Out[813]:
In [822]:
sns.distplot(trials.loc[trials['Target side']=='TargetRight','Movement endpoints'],
kde_kws={'label':'TargetRight'})
sns.distplot(trials.loc[trials['Target side']=='TargetLeft','Movement endpoints'],
kde_kws={'label':'TargetLeft'})
plt.xlim([-80,80])
plt.xlabel('Movement endpoint (degrees separation from midline)')
Out[822]:
In [839]:
trial_results = ed[(ed['Name']=='Neutral') | (ed['Name']=='Missed') | (ed['Name']=='Hit') | (ed['Name']=='Penalty')]
trials['Outcome'] = np.array(trial_results['Name'])
In [841]:
trials.head()
Out[841]:
In [864]:
mean_end_right = trials.loc[trials['Target side']=='TargetRight','Movement endpoints'].mean()
mean_end_left = trials.loc[trials['Target side']=='TargetLeft','Movement endpoints'].mean()
In [ ]:
In [ ]: