In [24]:
%reset -f
%pylab inline
import mutils.io as mio
import mutils.FDatAn as fda
import mutils.misc as mi
from mutils.io import build_dataset
import libshai.util as ut
import sys
conf = mio.saveable()
conf.subject = 3
conf.ttype=1
conf.dt_window = 30
conf.dt_medfilter = False
conf.startup_n_full_maps = 100
# cell_ID 1
# load kinematic data
ws1 = mio.saveable()
ws1.k = mio.KinData(data_dir = 'data/2011-mmcl_mat/')
ws1.k_doc = "generic kinematic data access object (empty selection)"
sys.stdout.write('loading data')
sys.stdout.flush()
ws1.k.load(conf.subject, conf.ttype)
ws1.SlipData = [mi.Struct(mio.mload('data/2011-mmcl_mat/SLIP/new/params3D_s%it%ir%i.dict' %
(conf.subject, conf.ttype, rep)))
for rep in ws1.k.reps]
sys.stdout.write(' done\n')
sys.stdout.flush()
# huge data selection, but not every dimension of every marker
ws1.full_markers = [ 'com_x', 'com_y', 'com_z',
'r_mtv_x - com_x',
'r_anl_x - com_x', 'r_kne_x - com_x',
'r_sia_x - com_x',
'l_mtv_x - com_x',
'l_anl_x - com_x', 'l_kne_x - com_x', 'l_sia_x - com_x',
'sacr_x - com_x',
'cvii_x - com_x',
'r_anl_y - com_y', 'r_kne_y - com_y', 'r_trc_y - com_y',
'r_sia_y - com_y',
'sacr_y - com_y', 'l_anl_y - com_y',
'l_kne_y - com_y', 'l_trc_y - com_y',
'l_sia_y - com_y',
'cvii_y - com_y',
'r_hea_x - com_x', 'l_hea_x - com_x', 'r_hea_y - com_y', 'l_hea_y - com_y',
'r_hea_z - com_z', 'l_hea_z - com_z',
'l_wrl_x - com_x', 'r_wrl_x - com_x',
'r_wrl_z - com_z', 'l_wrl_z - com_z',
'r_elb_y - com_y','l_elb_y - com_y',
'r_elb_x - com_x', 'l_elb_x - com_x',
'r_acr_y - com_y', 'l_acr_y - com_y',
'r_acr_z - com_z', 'l_acr_z - com_z',
'r_mtv_z - com_z',
'r_anl_z - com_z',
'r_trc_z - com_z',
'l_mtv_z - com_z',
'l_anl_z - com_z', 'l_trc_z - com_z',
'sacr_z - com_z',
]
ws1.notrunk_markers = [ 'com_x', 'com_y', 'com_z',
'r_mtv_x - com_x',
'r_anl_x - com_x', 'r_kne_x - com_x',
'r_sia_x - com_x',
'l_mtv_x - com_x',
'l_anl_x - com_x', 'l_kne_x - com_x', 'l_sia_x - com_x',
'sacr_x - com_x',
#'cvii_x - com_x',
'r_anl_y - com_y', 'r_kne_y - com_y', 'r_trc_y - com_y',
'r_sia_y - com_y',
'sacr_y - com_y', 'l_anl_y - com_y',
'l_kne_y - com_y', 'l_trc_y - com_y',
'l_sia_y - com_y',
#'cvii_y - com_y',
#'r_hea_x - com_x', 'l_hea_x - com_x', 'r_hea_y - com_y', 'l_hea_y - com_y',
#'r_hea_z - com_z', 'l_hea_z - com_z',
#'l_wrl_x - com_x', 'r_wrl_x - com_x',
#'r_wrl_z - com_z', 'l_wrl_z - com_z',
#'r_elb_y - com_y','l_elb_y - com_y',
#'r_elb_x - com_x', 'l_elb_x - com_x',
#'r_acr_y - com_y', 'l_acr_y - com_y',
#'r_acr_z - com_z', 'l_acr_z - com_z',
'r_mtv_z - com_z',
'r_anl_z - com_z',
'r_trc_z - com_z',
'l_mtv_z - com_z',
'l_anl_z - com_z', 'l_trc_z - com_z',
'sacr_z - com_z',
]
sys.stdout.write('building dataset "full"')
sys.stdout.flush()
ws1.k.selection = ws1.full_markers
ws1.dataset_full = build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
sys.stdout.write(' done\n')
sys.stdout.flush()
sys.stdout.write('computing maps for "full"')
sys.stdout.flush()
_, maps_1, _ = fda.fitData(ws1.dataset_full.s_kin_r, ws1.dataset_full.s_kin_l, nps=1,
nrep=conf.startup_n_full_maps, rcond=1e-7)
_, maps_2, _ = fda.fitData(ws1.dataset_full.s_kin_l[:-1, :], ws1.dataset_full.s_kin_r[1:, :],
nps=1, nrep=conf.startup_n_full_maps, rcond=1e-7, )
ws1.maps_full = [dot(m2, m1) for m1, m2 in zip(maps_1, maps_2)]
sys.stdout.write(' done\n')
sys.stdout.flush()
sys.stdout.write('building dataset "no trunk"')
sys.stdout.flush()
ws1.k.selection = ws1.notrunk_markers
ws1.dataset_full = build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
sys.stdout.write(' done\n')
sys.stdout.flush()
sys.stdout.write('computing maps for "no trunk"')
sys.stdout.flush()
_, maps_1, _ = fda.fitData(ws1.dataset_full.s_kin_r, ws1.dataset_full.s_kin_l, nps=1,
nrep=conf.startup_n_full_maps, rcond=1e-7)
_, maps_2, _ = fda.fitData(ws1.dataset_full.s_kin_l[:-1, :], ws1.dataset_full.s_kin_r[1:, :],
nps=1, nrep=conf.startup_n_full_maps, rcond=1e-7, )
ws1.maps_notrunk = [dot(m2, m1) for m1, m2 in zip(maps_1, maps_2)]
sys.stdout.write(' done\n')
sys.stdout.flush()
In [25]:
vif = ut.VisEig(127, False)
for A in ws1.maps_full:
vif.add(eig(A)[0])
vin = ut.VisEig(127, False)
for A in ws1.maps_notrunk:
vin.add(eig(A)[0])
figure(figsize=(6,6))
v = vif.vis1(rlabels=[.5])
v[0].set_cmap(cm.gray)
v = vin.vis1(rlabels=[.5])
v[0].set_cmap(cm.spring)
gca().set_xlim(-.5,.75)
gca().set_ylim(-.5,.5)
gca().set_xticklabels([])
gca().set_yticklabels([])
grid('on')
title('eigenvalues subject {} with and without trunk markers'.format(conf.subject))
subplots_adjust(bottom=.05,top=.9,left=.05,right=.95)
savefig('img/ev_notrunk_s{}.pdf'.format(conf.subject))
savefig('img/ev_notrunk_s{}.svg'.format(conf.subject))
In [13]:
#--------- visualize largest eigenvectors
m = fda.meanMat(ws1.maps_full)
ev, ew = eig(m)
order = argsort(abs(ev))[::-1]
print abs(ev[order[:3]])
figure(figsize=(24,3))
pcolor(abs(ew[:, order[:3]]).T)
gca().set_xticks(arange(len(ws1.dataset_full.kin_labels)) + .5)
gca().set_xticklabels(ws1.dataset_full.kin_labels, rotation=90)
gca().set_yticks(arange(3) + .5)
gca().set_yticklabels(arange(3) + 1)
gca().set_ylabel('ordinal of eigenvector')
title("magnitude of eigenvectors to largest eigenvalues" +
"\nmag. of eigenvalues: " + ' '.join(['{0:.3g}'.format(x,) for x in abs(ev[order[:3]])]))
pass
In [3]:
figure()
subplot(2,1,1)
plot(ws1.dataset_full.all_IC_rc[:,0],'r.',ws1.dataset_full.all_IC_lc[:,0],'g.')
title('CoM apex height r:right, g:left', fontsize=9)
gca().set_xticklabels([])
grid('on')
subplot(2,1,2)
plot(vstack(ws1.dataset_full.TR)[:,0],'r.',vstack(ws1.dataset_full.TL)[:,0],'g.')
grid('on')
title('step duration r:right, g:left', fontsize=9)
xlabel('stride #', fontsize=9)
subplots_adjust(hspace=.5, bottom=.15,top=.95,left=.1,right=.95)
savefig('img/trends_subj{}.pdf'.format(conf.subject))
In [8]:
import mutils.misc as mi
dt1 = lambda x: mi.dt_movingavg(x, 10, True)
dt2 = lambda x: mi.dt_movingavg(x, 10, False)
datr = vstack(ws1.dataset_full.TR)[:,0]
datl = vstack(ws1.dataset_full.TL)[:,0]
figure(figsize=(12,3))
plot(datr - mean(datr),'b.')
plot(dt1(datr),'r.')
plot(dt2(datr),'g.')
xlim(600,800)
Out[8]:
In [ ]: