Here, two figures are created:
LE:
In [1]:
%cd
%cd mmnotebooks
In [2]:
%%capture
import mutils.misc as mi
import libshai.util as ut
nbfile = "AnkleSLIP - find minimal model- Version 3.ipynb"
cell_ids = ['0.1','1','3','3.1','4', '4.2a', '4.2b', '4.2c',] # '3.2', '4.3', '4.3.1']
mi.run_nbcells(nbfile, ['0'])
conf.subject = 3
conf.exclude_IC_from_factors = False
conf.startup_n_full_maps = 200
conf.cslip_forceZeroRef = False
conf.quiet = True
conf.startup_compute_full_maps = True
mi.run_nbcells(nbfile, cell_ids)
In [3]:
# using the same font in mathematical expressions as elsewhere
import matplotlib as mpl
mpl.rcParams['mathtext.default'] = 'regular' # this is the key
import matplotlib.font_manager as fmt
FM = fmt.FontManager()
fnt = FM.findfont('Arial')
fnt_s = FM.findfont('Symbol')
fp = fmt.FontProperties(fname=fnt)
fp.set_size(9.0)
fp_s = fmt.FontProperties(fname=fnt_s)
fp_s.set_size(9.0)
d = ws1.dataset_full # shortcut
if conf.exclude_IC_from_factors:
def condition(lbl):
return '-' in lbl
else:
def condition(lbl):
return True
com_free_dat = vstack([d.s_kin_r[:, dim] for dim in range(d.s_kin_r.shape[1]) if
condition(d.kin_labels[dim])]).T
com_free_labels = [lbl for lbl in d.kin_labels if condition(lbl)]
_, aA, _ = fda.fitData(com_free_dat, ws1.dataset_full.s_param_r, nrep=100)
com_free_dat = vstack([d.s_kin_l[:, dim] for dim in range(d.s_kin_l.shape[1]) if
condition(d.kin_labels[dim])]).T
com_free_labels = [lbl for lbl in d.kin_labels if condition(lbl)]
_, aA_L, _ = fda.fitData(com_free_dat, ws1.dataset_full.s_param_l, nrep=100)
A = fda.meanMat(aA)
A_L = fda.meanMat(aA_L)
#u,s,v = svd(A, full_matrices=False)
# v is the row space projector (actually: v.T v would be the projector)
#v = dot(u,v) # now, each row corresponds to a specific SLIP parameter
v = vstack([A[row, :] / norm(A[row, :]) for row in range(A.shape[0])])
v_L = vstack([A_L[row, :] / norm(A_L[row, :]) for row in range(A_L.shape[0])])
# --- prepare colormap
from matplotlib.colors import LinearSegmentedColormap
cdict = {'red' : ( (0.0, 0.0, 0.0),
(0.5, 1.0, 1.0),
(1.0, 1.0, 1.0) ),
'green' : ( (0.0, 0.0, 0.0),
(0.5, 1.0, 1.0),
(1.0, 0.0, 0.0) ),
'blue' : ( (0.0, 1.0, 1.0),
(0.5, 1.0, 1.0),
(1.0, 0.0, 0.0) ),
}
map_bwr = LinearSegmentedColormap('BWR1', cdict)
# prepare font (for proc R soc)
#import matplotlib.font_manager as fmt
#FM = fmt.FontManager()
#fnt = FM.findfont('Times New Roman')
#fp = fmt.FontProperties(fname=fnt)
#fp.set_size(9.0)
In [4]:
# permute labels and columns
perm = arange(A.shape[1])
def swap(a,b):
tmp = perm[a]
perm[a] = perm[b]
perm[b] = tmp
swap(2, 8)
swap(36, 32)
swap(39, 33)
swap(46,40)
swap(49,57)
swap(53,59)
swap(62,65)
com_labels_reord = [com_free_labels[x] for x in perm]
lbls_r = [(nr, lbl) for nr, lbl in enumerate(com_labels_reord) if 'r_anl' in lbl or '-' not in lbl]
lbls_l = [(nr, lbl) for nr, lbl in enumerate(com_labels_reord) if 'l_anl' in lbl or '-' not in lbl]
# skip figure
if False:
figure()
subplot(1,2,1)
for elem in lbls_r:
text(.5, elem[0], elem[1])
gca().set_ylim(0,100)
subplot(1,2,2)
for elem in lbls_l:
text(.5, elem[0], elem[1])
gca().set_ylim(0,100)
In [5]:
# --- right step
f1 = figure(figsize=(6.83,1.8))
p = pcolor(v[:, perm], cmap=map_bwr)
# magnitude of columns
ax = gca()
ax2 = ax.twinx()
ax2.plot(arange(v.shape[1]) + .5, sqrt(sum(v[:, perm]**2, axis=0)),'v', markersize=3, color='#ffffff')
ax2.set_yticks([])
axes(ax)
gca().set_xlim(0,v.shape[1])
gca().set_xlim(0,v.shape[1])
# labels including CoM
#lbls = [(nr, lbl) for nr, lbl in enumerate(com_free_labels) if 'r_anl' in lbl or '-' not in lbl]
gca().set_xticks([lbl[0]+.5 for lbl in lbls_r])
# right labels:
r_lbls = [r'CoM$_y$', r'r. ankle$_z$', r'r. ankle$_x$', r'r. ankle$_y$', r'CoM$_{vz}$', r'CoM$_{vx}$',
r'r. ankle$_{vz}$', r'r. ankle$_{vx}$', r'r. ankle$_{vy}$',]
#gca().set_xticklabels([lbl[1].split(' ')[0] for lbl in lbls], rotation=45, fontproperties=fp,
# ha='right')
fp.set_size(8.0)
gca().set_xticklabels(r_lbls, fontproperties=fp, rotation=45, ha='right')
fp.set_size(9.0)
gca().set_yticks(arange(5) + .4)
gca().set_yticklabels(['k', r'$\alpha$', r'$L_0$', r'$\beta$', r'$\Delta E$'], fontproperties=fp_s,
)
p.set_clim(-1,1)
#title(r'normalized rows of the map "full state" $\to$ "SLIP parameter" for a right step',
# fontproperties=fp)
title('visualization of the Factor-SLIP feedback matrix for a right step', fontproperties=fp)
ylabel('parameter to predict', fontproperties=fp)
xlabel('state space variable', fontproperties=fp)
subplots_adjust(left=.075, right=.90, bottom=.4,wspace=.01)
# create colorbar manually
ax0 = axes([.925,.4,.025,.5])
cdat = linspace(-1,1,201)
ax0.set_xticks([])
ax0.set_yticks([])
ax = ax0.twinx()
ax.pcolor(array([[0,1],]),cdat, cdat[:, newaxis], cmap=map_bwr)
ax.set_xticks([])
ax.set_yticks([-1,-.5,0,.5,1])
fp.set_size(8.0)
ax.set_yticklabels([-1,-.5,0,.5,1], fontproperties=fp, )
ax.set_xlim(0.2, 0.8)
fp.set_size(9.0)
savefig('img/predictors_subj{}_r.pdf'.format(conf.subject))
savefig('img/predictors_subj{}_r.svg'.format(conf.subject))
# --- left step
f2 = figure(figsize=(6.83,2))
p = pcolor(v_L[:, perm], cmap=map_bwr)
ax = gca()
ax2 = ax.twinx()
ax2.plot(arange(v_L.shape[1]) + .5, sqrt(sum(v_L[:, perm]**2, axis=0)),'v', markersize=3, color='#ffffff')
ax2.set_yticks([])
axes(ax)
gca().set_xlim(0,v_L.shape[1])
# labels including CoM
#lbls_l = [(nr, lbl) for nr, lbl in enumerate(com_free_labels) if 'l_anl' in lbl or '-' not in lbl]
gca().set_xticks([lbl[0]+.5 for lbl in lbls_l])
# right labels:
l_lbls = [r'CoM$_z$', r'l. ankle$_x$', r'l. ankle$_y$', r'l. ankle$_z$', r'CoM$_{vx}$', r'CoM$_{vy}$',
r'l. ankle$_{vx}$', r'l. ankle$_{vy}$', r'l. ankle$_{vz}$',]
#gca().set_xticklabels([lbl[1].split(' ')[0] for lbl in lbls], rotation=45, fontproperties=fp,
# ha='right')
fp.set_size(8.0)
gca().set_xticklabels(l_lbls, fontproperties=fp, rotation=45, ha='right')
fp.set_size(9.0)
gca().set_yticks(arange(5) + .4)
gca().set_yticklabels(['k', r'$\alpha$', r'$L_0$', r'$\beta$', r'$\Delta E$'], fontproperties=fp_s,
)
p.set_clim(-1,1)
#title(r'normalized rows of the map "full state" $\to$ "SLIP parameter" for a left step',
# fontproperties=fp)
title('visualization of the Factor-SLIP feedback matrix for a left step', fontproperties=fp)
ylabel('parameter to predict', fontproperties=fp)
xlabel('state space variable', fontproperties=fp)
subplots_adjust(left=.075, right=.90, bottom=.4,wspace=.01)
# create colorbar manually
ax0 = axes([.925,.4,.025,.5])
cdat = linspace(-1,1,201)
ax0.set_xticks([])
ax0.set_yticks([])
ax = ax0.twinx()
ax.pcolor(array([[0,1],]),cdat, cdat[:, newaxis], cmap=map_bwr)
ax.set_xticks([])
ax.set_yticks([-1,-.5,0,.5,1])
fp.set_size(8.0)
ax.set_yticklabels([-1,-.5,0,.5,1], fontproperties=fp, )
ax.set_xlim(0.2, 0.8)
fp.set_size(9.0)
savefig('img/predictors_subj{}_l.pdf'.format(conf.subject))
savefig('img/predictors_subj{}_l.svg'.format(conf.subject))
print 'figures stored as img/predictors_subj{}_l/r.pdf'.format(conf.subject)
pass
In [6]:
#t = ws1.cslip.fac_slip
#IC = ws1.cslip.fac_slip_refstate_r
#print t(IC) - IC
print ws1.cslip.ankle_slip_refstate_r
print ws1.cslip.aug_slip_refstate_r
print ws1.cslip.fac_slip_refstate_r
In [7]:
mean(ws1.fscore_l, axis=0)
mi.run_nbcells(nbfile, ['3.2'])
In [8]:
# --- compute eigenvalues of simulated and affine models
ws = mio.saveable()
from mutils.io import cacherun
# --- --- model 1: factor SLIP
J1 = []
h = .01
f = ws1.cslip.fac_slip
for elem in range(len(ws1.cslip.fac_slip_refstate_r)):
IC = ws1.cslip.fac_slip_refstate_r.copy()
IC[elem] += h
fsp = f(IC)
IC[elem] -= 2*h
fsn = f(IC)
J1.append((fsp - fsn) / (2.*h))
ws.J1 = vstack(J1).T
ws.J1_doc = "linearized map of controlled SLIP 1: factor SLIP"
# --- --- affine model
_, maps_1, _ = cacherun(fda.fitData,ws1.reddat_r, ws1.reddat_l, nps=1, nrep=200, rcond=1e-7)
_, maps_2, _ = cacherun(fda.fitData,ws1.reddat_l[:-1, :], ws1.reddat_r[1:, :], nps=1, nrep=200, rcond=1e-7, )
ws.maps_1 = [dot(m2, m1) for m1, m2 in zip(maps_1, maps_2)]
ws.maps_1_doc = "maps for model 1: factor SLIP (affine)"
# --- --- model 2: ankle slip
J2 = []
h = .001
f = ws1.cslip.ankle_slip
for elem in range(len(ws1.cslip.ankle_slip_refstate_r)):
IC = ws1.cslip.ankle_slip_refstate_r.copy()
IC[elem] += h
fsp = f(IC)
IC[elem] -= 2*h
fsn = f(IC)
J2.append((fsp - fsn) / (2.*h))
ws.J2 = vstack(J2).T
ws.J2_doc = "linearized map of controlled SLIP 2: ankle SLIP"
# --- --- affine model
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'r_anl_y - com_y',
'r_anl_z - com_z', 'r_anl_x - com_x',]
tmp_dsr = build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'l_anl_y - com_y',
'l_anl_z - com_z', 'l_anl_x - com_x',]
tmp_dsl = build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
# Re-ordering the dimensions is just for completeness. It does not affect the eigenvalues (verified).
rix=reord_idx(tmp_dsr.kin_labels)
lix=reord_idx(tmp_dsl.kin_labels)
idat_ra = fda.dt_movingavg(hstack([tmp_dsr.all_IC_rc, tmp_dsr.all_kin_r[:, rix][:,3:]]),
conf.dt_window, conf.dt_medfilter)
idat_la = fda.dt_movingavg(hstack([tmp_dsl.all_IC_lc, tmp_dsl.all_kin_l[:, lix][:,3:]]),
conf.dt_window, conf.dt_medfilter)
#idat_r[:,3:] -= mean(idat_r[:, 3:], axis=0)
#idat_l[:,3:] -= mean(idat_l[:, 3:], axis=0)
_, maps_1a, _ = cacherun( fda.fitData, idat_ra, idat_la, nps=1, nrep=200, rcond=1e-8)
_, maps_2a, _ = cacherun( fda.fitData, idat_la[:-1, :], idat_ra[1:, :], nps=1, nrep=200, rcond=1e-8,)
ws.maps_2 = [dot(m2, m1) for m1, m2 in zip(maps_1a, maps_2a)]
ws.maps_2_doc = "maps for model 2: ankle SLIP (affine)"
# --- --- model 3: augmented SLIP
J3 = []
h = .001
h_scalings = ones(8) # possibility to adapt different magnitudes of data
f = ws1.cslip.aug_slip
for elem in range(len(ws1.cslip.aug_slip_refstate_r)):
IC = ws1.cslip.aug_slip_refstate_r.copy()
IC[elem] += h * h_scalings[elem]
fsp = f(IC)
IC[elem] -= 2*h* h_scalings[elem]
fsn = f(IC)
J3.append((fsp - fsn) / (2.*h*h_scalings[elem]))
ws.J3 = vstack(J3).T
ws.J3_doc = "linearized map of controlled SLIP 3: augmented SLIP"
idat_r = hstack([ws1.dataset_full.all_IC_r, roll(ws1.dataset_full.all_param_l, 1, axis=0)])
idat_l = hstack([ws1.dataset_full.all_IC_l, ws1.dataset_full.all_param_r])
idat_r = fda.dt_movingavg(idat_r, conf.dt_window, conf.dt_medfilter)
idat_l = fda.dt_movingavg(idat_l, conf.dt_window, conf.dt_medfilter)
_, maps_1, _ = fda.fitData(idat_r, idat_l, nps=1, nrep=200, rcond=1e-7)
_, maps_2, _ = fda.fitData(idat_l[:-1, :], idat_r[1:, :], nps=1, nrep=200, rcond=1e-7, )
ws.maps_3 = [dot(m2, m1) for m1, m2 in zip(maps_1, maps_2)]
ws.maps_3_doc = "maps for model 3: augmented SLIP (affine)"
In [9]:
# --- compute eigenvalues for different numbers of sections: 1,2,5
# --- --- compute eigenvalues of full state with multiple sections (try n=7!)
import sys
import mutils.io as mio
reload(mio)
# --- --- compute multiple-section map n=3
nps = 3
ws1.k.selection = ws1.full_markers
d1d = fda.dt_movingavg(ws1.k.make_1D(nps=nps, phases_list=ws1.dataset_full.all_phases_r),
conf.dt_window, conf.dt_medfilter)[:,2*nps:]
print "3 sections:"
print "=" * nps
pt_maps = []
for rep in range(nps):
erep = rep + 1 if rep < nps - 1 else 0
skip = None if rep < nps - 1 else -1
start = None if rep < nps - 1 else 1
sys.stdout.write('.')
sys.stdout.flush()
_, maps1, _ = mio.cacherun(fda.fitData, d1d[:skip,rep::nps], d1d[start:,erep::nps],
nps=1, rcond=1e-7, nrep=200)
pt_maps.append(maps1)
ws.fmaps_3 = [reduce(dot, ptm) for ptm in zip(*pt_maps[::-1])]
print "done"
# recompute these mappings
print "computing partial mappings with same samples"
sctMaps, fmaps, _ = mio.cacherun(fda.fitData, d1d[:-1,:], d1d[1:,:], sections=[0,1,2],
nps=3, nrep=200, rcond=1e-7)
ws.fmaps_3s = [reduce(dot, ptm) for ptm in zip(*sctMaps[::-1])]
print "done"
# --- --- compute multiple-section map n=5
nps = 5
ws1.k.selection = ws1.full_markers
d1d = fda.dt_movingavg(ws1.k.make_1D(nps=nps, phases_list=ws1.dataset_full.all_phases_r),
conf.dt_window, conf.dt_medfilter)[:,2*nps:]
print "5 sections:"
print "=" * nps
pt_maps = []
for rep in range(nps):
erep = rep + 1 if rep < nps - 1 else 0
skip = None if rep < nps - 1 else -1
start = None if rep < nps - 1 else 1
sys.stdout.write('.')
sys.stdout.flush()
_, maps1, _ = mio.cacherun( fda.fitData, d1d[:skip,rep::nps], d1d[start:,erep::nps],
nps=1, rcond=1e-7, nrep=200)
pt_maps.append(maps1)
ws.fmaps_5 = [reduce(dot, ptm) for ptm in zip(*pt_maps[::-1])]
print "done"
print "computing partial mappings with same samples"
sctMaps, fmaps, _ = mio.cacherun(fda.fitData, d1d[:-1,:], d1d[1:,:], sections=[0,1,2,3,4],
nps=5, nrep=200, rcond=1.1e-7)
ws.fmaps_5s = [reduce(dot, ptm) for ptm in zip(*sctMaps[::-1])]
print "done"
# --- for a single section
print "1 section ...",
_, ws.fmaps_1, _ = mio.cacherun(fda.fitData, d1d[:-1, 0::nps], d1d[1:, 0::nps], nps=1,
rcond=1e-7, nrep=200)
print "done"
In [10]:
vi3 = ut.VisEig(127, False)
vi3s = ut.VisEig(127, False)
vi5 = ut.VisEig(127, False)
vi5s = ut.VisEig(127, False)
for A in ws.fmaps_3s:
vi3s.add(eig(A)[0])
for A in ws.fmaps_3:
vi3.add(eig(A)[0])
for A in ws.fmaps_5:
vi5.add(eig(A)[0])
for A in ws.fmaps_5s:
vi5s.add(eig(A)[0])
In [11]:
figure(figsize=(14,6))
subplot(1,2,1)
v = vi3s.vis1(rlabels=[.5,])
v[0].set_cmap(cm.winter)
v = vi3.vis1(rlabels=[.5,])
v[0].set_cmap(cm.spring)
title('3 sections\nsame sample (red), random sample (blue)')
gca().set_ylim(-.5, .5)
gca().set_xlim(-.5, 1)
subplot(1,2,2)
v = vi5.vis1(rlabels=[.5,])
v[0].set_cmap(cm.winter)
v = vi5s.vis1(rlabels=[.5,])
v[0].set_cmap(cm.spring)
title('5 sections\nsame sample (red), random sample (blue)')
gca().set_ylim(-.5, .5)
gca().set_xlim(-.5, 1)
#savefig('img/ev_sameSamples_s{}.pdf'.format(conf.subject))
Out[11]:
In [12]:
conf
Out[12]:
In [13]:
# --- create 2D histograms of EV distributions (affine models)
# --- --- "base" model
reload(ut)
vi_2s = ut.VisEig(127, False)
for A in ws1.maps_full:
vi_2s.add(eig(A)[0])
# --- --- factor SLIP (affine)
vi1 = ut.VisEig(127, False)
for A in ws.maps_1:
vi1.add(eig(A)[0])
# --- --- ankle SLIP (affine)
vi2 = ut.VisEig(127, False)
for A in ws.maps_2:
vi2.add(eig(A)[0])
# --- --- augmented SLIP (affine)
vi3 = ut.VisEig(127, False)
for A in ws.maps_3:
vi3.add(eig(A)[0])
# --- other number of sections Floquet models
vi_1s = ut.VisEig(127, False)
for A in ws.fmaps_1:
vi_1s.add(eig(A)[0])
vi_3s = ut.VisEig(127, False)
for A in ws.fmaps_3:
vi_3s.add(eig(A)[0])
vi_5s = ut.VisEig(127, False)
for A in ws.fmaps_5:
vi_5s.add(eig(A)[0])
In [14]:
import mutils.plotting as mplt
colors_ = mplt.colorset_distinct
from matplotlib.colors import LinearSegmentedColormap
# continous black-red-black
cdict_1 = {'red' : ( (0.0, 0.9, 0.9),
(1.0, 0.6, 0.6)),
'green' : ((0.0, 0.3, 0.3),
(1.0, 0.15, 0.15),
),
'blue' : ( (0.0, 0.34, 0.34),
(1.0, 0.2, 0.2),
),
}
map_r1 = LinearSegmentedColormap('r1', cdict_1)
cdict_2 = {'red' : ( (0.0, 0.3, 0.3),
(1.0, 0.2, 0.2) ),
'green' : ( (0.0, 0.65, 0.65),
(1.0, 0.43, 0.43) ),
'blue' : ( (0.0, 0.94, 0.94),
(1.0, 0.65, 0.65) ),
}
map_b1 = LinearSegmentedColormap('b1', cdict_2)
figure()
data = arange(50)[newaxis, :]
subplot(2, 1, 1)
pcolor(data, cmap=map_r1)
subplot(2, 1, 2)
pcolor(data, cmap=map_b1)
# compute corresponding hex colors (mean color)
# rgb2hex
import struct
rgbs = [int(elem) for elem in (0.75*255, 0.225*255, 0.27*255)]
color_red = "#{}".format(struct.pack('BBB',*rgbs).encode('hex').upper())
rgbs = [int(elem) for elem in (0.25*255, 0.54*255, 0.795*255)]
color_blue = "#{}".format(struct.pack('BBB',*rgbs).encode('hex').upper())
pass
In [13]:
In [14]:
# --- create the figure
# OLD CODE (only average for simulated model is given) - see step4!
# MAYBE(!) this is required because of imports ... I have to test this.
%config InlineBackend.close_figures = False
import matplotlib.gridspec as gridspec
#Colorbrewer: Someone has already figured out what colors are good to use, so use them.
#'Data-pixel ratio': 'Erase' chart items until the next thing you erase would remove data
# --- define some formats
#ffmt2 = {'family' : 'serif',
# 'serif' : ['Times', 'Times New Roman'],
# 'size' : 9}
#matplotlib.rc('font', **ffmt2)
import matplotlib.font_manager as fmt
FM = fmt.FontManager()
#fnt = FM.findfont('Times New Roman') # ProcRSoc B
fnt = FM.findfont('Arial')
fp = fmt.FontProperties(fname=fnt)
fp.set_size(9.0)
import matplotlib as mpl
mpl.rcParams['mathtext.default'] = 'regular' # this is the key
gs1 = gridspec.GridSpec(1,3, left=.025, right=.7, wspace=0.05, top=.85, bottom=.05)
gs2 = gridspec.GridSpec(2,2, left=.75, right=.975, wspace=0.05, hspace=.05, top=.85, bottom=.05)
figure(figsize=(6.83,1.6)) # width of Proc R Soc 17.5 cm double column (84mm single)
def fmt(ax):
ax.set_xlim(-.55,.85)
ax.set_ylim(-.55,.55)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.grid('on')
def fmt2(ax):
ax.set_xlim(-.55,.85)
ax.set_ylim(-.55,.55)
ax.set_xticks([])
ax.set_yticks([])
def ntxt(ax, txt):
ax.text(.35,-.5, txt, fontproperties=fp, va='bottom',
bbox=dict(facecolor='#ffffff', edgecolor='#ffffff', pad=1.0))
# --- --- plot 1
subplot(gs1[0,0])
v0 = vi_2s.vis1(N=3, rlabels=[.5,], fontproperties=fp)
v0[0].set_cmap(cm.gray)
v = vi1.vis1(N=3, rlabels=[],fontproperties=fp, lw=1)
v[0].set_cmap(map_r1)
for ev in eig(ws.J1)[0]:
plot(ev.real, ev.imag, 'bx', markeredgewidth=1)
fmt(gca())
title('Factor-SLIP', fontproperties=fp)
# --- --- plot 2
subplot(gs1[0,1])
v0 = vi_2s.vis1(N=3, rlabels=[.5,], fontproperties=fp)
v0[0].set_cmap(cm.gray)
v = vi2.vis1(N=5, rlabels=[], fontproperties=fp, lw=2)
v[0].set_cmap(map_r1)
for ev in eig(ws.J2)[0]:
plot(ev.real, ev.imag, 'bx', markeredgewidth=1)
fmt(gca())
title('Ankle-SLIP', fontproperties=fp)
# --- --- plot 3
subplot(gs1[0,2])
v0 = vi_2s.vis1(N=3, rlabels=[.5,], fontproperties=fp)
v0[0].set_cmap(cm.gray)
v = vi3.vis1(N=4, rlabels=[], fontproperties=fp)
v[0].set_cmap(map_r1)
for ev in eig(ws.J3)[0]:
plot(ev.real, ev.imag, 'bx', markeredgewidth=1)
fmt(gca())
title('Augmented-SLIP', fontproperties=fp)
# --- --- plots for different sections
#ax_x = subplot(1,4,4, frameon=False)
#ax_x.set_title('different sections')
#ax_x.set_xticks([])
#ax_x.set_yticks([])
# only if multi maps have been computed
if 'vi_1s' in locals() and not allclose(vi_1s.H, 0):
subplot(gs2[0,0])
v0 = vi_1s.vis1(N=5, rlabels=[ .5], fontproperties=fp)
v0[0].set_cmap(cm.gray)
fmt2(gca())
ntxt(gca(), "N=1")
subplot(gs2[0,1])
v = vi_2s.vis1(N=5, rlabels=[ .5], fontproperties=fp)
v[0].set_cmap(cm.gray)
fmt2(gca())
ntxt(gca(), "N=2")
subplot(gs2[1,0])
v = vi_3s.vis1(N=5, rlabels=[ .5], fontproperties=fp)
v[0].set_cmap(cm.gray)
fmt2(gca())
ntxt(gca(), "N=3")
subplot(gs2[1,1])
v = vi_5s.vis1(N=5, rlabels=[ .5], fontproperties=fp)
v[0].set_cmap(cm.gray)
fmt2(gca())
ntxt(gca(), "N=5")
# --- --- general formatting
#savefig('img/fig_ev_subj{}.pdf'.format(conf.subject))
#savefig('img/fig_ev_subj{}.svg'.format(conf.subject))
pass
In [19]:
savefig('tmp/tmp.pdf')
In [15]:
import libshai.util as ut
In [16]:
# hint: edit conf.cslip_forceZeroRef (True / False)
vi1_fac = ut.VisEig(127, False)
vi2_fac = ut.VisEig(127, False)
# create factor-model
# old code
# cell_ID 4.2a
dt_fstate_r = hstack([fda.dt_movingavg(ws1.dataset.all_IC_r,
conf.dt_window, conf.dt_medfilter), ws1.fscore_r])
dt_fstate_l = hstack([fda.dt_movingavg(ws1.dataset.all_IC_l,
conf.dt_window, conf.dt_medfilter), ws1.fscore_l])
_, m1, _ = fda.fitData(dt_fstate_r, dt_fstate_l, nrep=200, rcond=1e-8)
_, m2, _ = fda.fitData(dt_fstate_l[:-1, :], dt_fstate_r[1:, :], nrep=200, rcond=1e-8)
for A in [dot(ma, mb) for ma, mb in zip(m2, m1)]:
vi1_fac.add(eig(A)[0])
if conf.cslip_forceZeroRef:
dt_fstate_r -= mean(dt_fstate_r, axis=0)
dt_fstate_l -= mean(dt_fstate_l, axis=0)
# compute parameter prediction maps
_, all_ctrl_r, _ = fda.fitData(dt_fstate_r, fda.dt_movingavg(ws1.dataset.all_param_r,
conf.dt_window, conf.dt_medfilter), nps=1,
nrep=200, sections=[0,], rcond=1e-8)
_, all_ctrl_l, _ = fda.fitData(dt_fstate_l, fda.dt_movingavg(ws1.dataset.all_param_l,
conf.dt_window, conf.dt_medfilter), nps=1,
nrep=200, sections=[0,], rcond=1e-8)
# compute factors prediction maps - "propagators" of factor's state
_, all_facprop_r, _ = fda.fitData(dt_fstate_r, ws1.fscore_l, nps=1, nrep=200, sections=[0,], rcond=1e-8)
_, all_facprop_l, _ = fda.fitData(dt_fstate_l[:-1, :], ws1.fscore_r[1:, :], nps=1, nrep=200, sections=[0,],
rcond=1e-8)
# set this to zero to force "zero" as reference for the fscores!
allow_nonzero_ref = 0. if conf.cslip_forceZeroRef else 1.
fac_slip_refstate_r = hstack(
[ws1.cslip.ICp_r, allow_nonzero_ref * mean(ws1.fscore_r, axis=0)]) # the latter should be zero anyway
fac_slip_refstate_l = hstack(
[ws1.cslip.ICp_l, allow_nonzero_ref * mean(ws1.fscore_l, axis=0)]) # the latter should be zero anyway
all_J = []
for cr, cl, pr, pl in zip(all_ctrl_r, all_ctrl_l,
all_facprop_r, all_facprop_l):
facslip = get_auto_sys(ws1.cslip.Pp_r, ws1.cslip.Pp_l,
fac_slip_refstate_r, fac_slip_refstate_l,
cr,
cl,
pr,
pl,
{'m': ws1.cslip.mass, 'g':-9.81})
J1 = []
h = .01
f = facslip
for elem in range(len(ws1.cslip.fac_slip_refstate_r)):
IC = ws1.cslip.fac_slip_refstate_r.copy()
IC[elem] += h
fsp = f(IC)
IC[elem] -= 2*h
fsn = f(IC)
J1.append((fsp - fsn) / (2.*h))
all_J.append(vstack(J1).T)
for A in all_J:
vi2_fac.add(eig(A)[0])
print "eigenvalues for factor-SLIP computed"
In [17]:
# hint: edit conf.cslip_forceZeroRef (True / False)
vi1_ank = ut.VisEig(127, False)
vi2_ank = ut.VisEig(127, False)
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'r_anl_y - com_y', 'r_anl_z - com_z', 'r_anl_x - com_x']
dataset_r = build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'l_anl_y - com_y', 'l_anl_z - com_z', 'l_anl_x - com_x']
dataset_l= build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
ws1.k_doc = "KinData object for s:" + str(ws1.k.subject) + " t:" + str(ws1.k.ttype) + " with left ankle"
states_r = dataset_r.all_kin_r[:, reord_idx(dataset_r.kin_labels)]
states_l = dataset_l.all_kin_l[:, reord_idx(dataset_l.kin_labels)]
# --- --- additional data
d = ws1.dataset_full # shortcut
# use interpolated apex states
states_r[:,:3] = d.all_IC_r
states_l[:,:3] = d.all_IC_l
dt_fstate_r = fda.dt_movingavg(states_r, conf.dt_window, conf.dt_medfilter)
dt_fstate_l = fda.dt_movingavg(states_l, conf.dt_window, conf.dt_medfilter)
#dt_fstate_l = hstack([fda.dt_movingavg(ws1.dataset.all_IC_l,
# conf.dt_window, conf.dt_medfilter), ws1.fscore_l])
_, m1, _ = fda.fitData(dt_fstate_r, dt_fstate_l, nrep=200, rcond=1e-8)
_, m2, _ = fda.fitData(dt_fstate_l[:-1, :], dt_fstate_r[1:, :], nrep=200, rcond=1e-8)
for A in [dot(ma, mb) for ma, mb in zip(m2, m1)]:
vi1_ank.add(eig(A)[0])
if conf.cslip_forceZeroRef:
dt_fstate_r -= mean(dt_fstate_r, axis=0)
dt_fstate_l -= mean(dt_fstate_l, axis=0)
# compute parameter prediction maps
_, all_ctrl_r, _ = fda.fitData(dt_fstate_r, fda.dt_movingavg(ws1.dataset.all_param_r,
conf.dt_window, conf.dt_medfilter), nps=1,
nrep=200, sections=[0,], rcond=1e-8)
_, all_ctrl_l, _ = fda.fitData(dt_fstate_l, fda.dt_movingavg(ws1.dataset.all_param_l,
conf.dt_window, conf.dt_medfilter), nps=1,
nrep=200, sections=[0,], rcond=1e-8)
# compute factors prediction maps - "propagators" of factor's state
_, all_facprop_r, _ = fda.fitData(dt_fstate_r, dt_fstate_l[:, 3:], nps=1, nrep=200, sections=[0,], rcond=1e-8)
_, all_facprop_l, _ = fda.fitData(dt_fstate_l[:-1, :], dt_fstate_r[1:, 3:], nps=1, nrep=200, sections=[0,],
rcond=1e-8)
# set this to zero to force "zero" as reference for the fscores!
allow_nonzero_ref = 0. if conf.cslip_forceZeroRef else 1.
fac_slip_refstate_r = hstack(
[ws1.cslip.ICp_r, allow_nonzero_ref * mean(dt_fstate_r[:, 3:], axis=0)]) # the latter should be zero anyway
fac_slip_refstate_l = hstack(
[ws1.cslip.ICp_l, allow_nonzero_ref * mean(dt_fstate_l[:, 3:], axis=0)]) # the latter should be zero anyway
all_J = []
for cr, cl, pr, pl in zip(all_ctrl_r, all_ctrl_l,
all_facprop_r, all_facprop_l):
facslip = get_auto_sys(ws1.cslip.Pp_r, ws1.cslip.Pp_l,
fac_slip_refstate_r, fac_slip_refstate_l,
cr,
cl,
pr,
pl,
{'m': ws1.cslip.mass, 'g':-9.81})
J1 = []
h = .01
f = facslip
for elem in range(len(fac_slip_refstate_r)):
IC = fac_slip_refstate_r.copy()
IC[elem] += h
fsp = f(IC)
IC[elem] -= 2*h
fsn = f(IC)
J1.append((fsp - fsn) / (2.*h))
all_J.append(vstack(J1).T)
for A in all_J:
vi2_ank.add(eig(A)[0])
print "eigenvalues for ankle-SLIP computed"
In [18]:
# hint: edit conf.cslip_forceZeroRef (True / False)
vi1_aug = ut.VisEig(127, False)
vi2_aug = ut.VisEig(127, False)
states_r = hstack([ws1.dataset_full.all_IC_r, roll(ws1.dataset_full.s_param_l, 1, axis=0)])
states_l = hstack([ws1.dataset_full.all_IC_l, ws1.dataset_full.s_param_r])
dt_fstate_r = fda.dt_movingavg(states_r, conf.dt_window, conf.dt_medfilter)
dt_fstate_l = fda.dt_movingavg(states_l, conf.dt_window, conf.dt_medfilter)
#dt_fstate_l = hstack([fda.dt_movingavg(ws1.dataset.all_IC_l,
# conf.dt_window, conf.dt_medfilter), ws1.fscore_l])
_, m1, _ = fda.fitData(dt_fstate_r, dt_fstate_l, nrep=200, rcond=1e-8)
_, m2, _ = fda.fitData(dt_fstate_l[:-1, :], dt_fstate_r[1:, :], nrep=200, rcond=1e-8)
for A in [dot(ma, mb) for ma, mb in zip(m2, m1)]:
vi1_aug.add(eig(A)[0])
if conf.cslip_forceZeroRef:
dt_fstate_r -= mean(dt_fstate_r, axis=0)
dt_fstate_l -= mean(dt_fstate_l, axis=0)
# compute parameter prediction maps
_, all_ctrl_r, _ = fda.fitData(dt_fstate_r, fda.dt_movingavg(ws1.dataset.all_param_r,
conf.dt_window, conf.dt_medfilter), nps=1,
nrep=200, sections=[0,], rcond=1e-8)
_, all_ctrl_l, _ = fda.fitData(dt_fstate_l, fda.dt_movingavg(ws1.dataset.all_param_l,
conf.dt_window, conf.dt_medfilter), nps=1,
nrep=200, sections=[0,], rcond=1e-8)
# compute factors prediction maps - "propagators" of factor's state
_, all_facprop_r, _ = fda.fitData(dt_fstate_r, dt_fstate_l[:, 3:], nps=1, nrep=200, sections=[0,], rcond=1e-8)
_, all_facprop_l, _ = fda.fitData(dt_fstate_l[:-1, :], dt_fstate_r[1:, 3:], nps=1, nrep=200, sections=[0,],
rcond=1e-8)
# set this to zero to force "zero" as reference for the fscores!
allow_nonzero_ref = 0. if conf.cslip_forceZeroRef else 1.
fac_slip_refstate_r = hstack(
[ws1.cslip.ICp_r, allow_nonzero_ref * mean(dt_fstate_r[:, 3:], axis=0)]) # the latter should be zero anyway
fac_slip_refstate_l = hstack(
[ws1.cslip.ICp_l, allow_nonzero_ref * mean(dt_fstate_l[:, 3:], axis=0)]) # the latter should be zero anyway
all_J = []
for cr, cl, pr, pl in zip(all_ctrl_r, all_ctrl_l,
all_facprop_r, all_facprop_l):
facslip = get_auto_sys(ws1.cslip.Pp_r, ws1.cslip.Pp_l,
fac_slip_refstate_r, fac_slip_refstate_l,
cr,
cl,
pr,
pl,
{'m': ws1.cslip.mass, 'g':-9.81})
J1 = []
h = .005
f = facslip
for elem in range(len(fac_slip_refstate_r)):
IC = fac_slip_refstate_r.copy()
IC[elem] += h
fsp = f(IC)
IC[elem] -= 2*h
fsn = f(IC)
J1.append((fsp - fsn) / (2.*h))
all_J.append(vstack(J1).T)
for A in all_J:
vi2_aug.add(eig(A)[0])
print "eigenvalues for augmented SLIP computed"
In [19]:
# --- create the figure
%config InlineBackend.close_figures = False
import matplotlib.gridspec as gridspec
thin_black = True
thin_lw = .65
#Colorbrewer: Someone has already figured out what colors are good to use, so use them.
#'Data-pixel ratio': 'Erase' chart items until the next thing you erase would remove data
# --- define some formats
#ffmt2 = {'family' : 'serif',
# 'serif' : ['Times', 'Times New Roman'],
# 'size' : 9}
#matplotlib.rc('font', **ffmt2)
import matplotlib.font_manager as fmt_m
FM = fmt_m.FontManager()
#fnt = FM.findfont('Times New Roman') # ProcRSoc B
fnt = FM.findfont('Arial')
fp = fmt_m.FontProperties(fname=fnt)
fp.set_size(9.0)
import matplotlib as mpl
mpl.rcParams['mathtext.default'] = 'regular' # this is the key
gs1 = gridspec.GridSpec(1,3, left=.05, right=.7125, wspace=0.05, top=.75, bottom=.125)
gs2 = gridspec.GridSpec(2,2, left=.75, right=.975, wspace=0.05, hspace=.05, top=.75, bottom=.125)
#figure(figsize=(6.83,1.6)) # width of Proc R Soc 17.5 cm double column (84mm single)
fig = figure(figsize=(18./2.56, 2)) # Ncomms format
def fmt(ax):
ax.set_xlim(-.55,.85)
ax.set_ylim(-.55,.55)
ax.set_xticks([0,])
ax.set_yticks([0,])
ax.set_xticklabels([])
ax.set_yticklabels([])
for elem in ax.get_children():
if isinstance(elem, matplotlib.collections.LineCollection):
# adjust linewidhts of colored lines only?
if norm(elem.get_edgecolor()[:3] - mean(elem.get_edgecolor()[:3])) > 0.1:
#if thin_black:
elem.set_linewidths(thin_lw)
ax.grid('on')
def fmt2(ax):
ax.set_xlim(-.55,.85)
ax.set_ylim(-.55,.55)
ax.set_xticks([])
ax.set_yticks([])
for elem in ax.get_children():
if isinstance(elem, matplotlib.collections.LineCollection):
# adjust linewidhts of colored lines only?
if thin_black:
elem.set_linewidths(thin_lw)
def ntxt(ax, txt):
ax.text(.35,-.5, txt, fontproperties=fp, va='bottom',
bbox=dict(facecolor='#ffffff', edgecolor='#ffffff', pad=1.0))
# --- --- plot 1
subplot(gs1[0,0])
v0 = vi_2s.vis1(N=3, rlabels=[.5,], fontproperties=fp)
v0[0].set_cmap(cm.gray)
v = vi1_fac.vis1(N=5, rlabels=[],fontproperties=fp,)
v[0].set_cmap(map_r1)
v = vi2_fac.vis1(N=3, rlabels=[ .5], fontproperties=fp)
v[0].set_cmap(map_b1)
fmt(gca())
gca().set_ylabel('imaginary axis', fontproperties=fp)
title('Factor-SLIP', fontproperties=fp)
# --- --- plot 2
subplot(gs1[0,1])
v0 = vi_2s.vis1(N=3, rlabels=[.5,], fontproperties=fp)
v0[0].set_cmap(cm.gray)
v = vi1_ank.vis1(N=5, rlabels=[], fontproperties=fp, lw=2)
v[0].set_cmap(map_r1)
v = vi2_ank.vis1(N=3, rlabels=[], fontproperties=fp, lw=2)
v[0].set_cmap(map_b1)
#for ev in eig(ws.J2)[0]:
# plot(ev.real, ev.imag, 'bx', markeredgewidth=1)
fmt(gca())
gca().set_xlabel('real axis', fontproperties=fp)
title('Ankle-SLIP', fontproperties=fp)
# --- --- plot 3
subplot(gs1[0,2])
v0 = vi_2s.vis1(N=3, rlabels=[.5,], fontproperties=fp)
v0[0].set_cmap(cm.gray)
v = vi1_aug.vis1(N=5, rlabels=[], fontproperties=fp, lw=2)
v[0].set_cmap(map_r1)
v = vi2_aug.vis1(N=3, rlabels=[], fontproperties=fp, lw=2)
v[0].set_cmap(map_b1)
fmt(gca())
title('Augmented-SLIP', fontproperties=fp)
# --- --- plots for different sections
#ax_x = subplot(1,4,4, frameon=False)
#ax_x.set_title('different sections')
#ax_x.set_xticks([])
#ax_x.set_yticks([])
# only if multi maps have been computed
if 'vi_1s' in locals() and not allclose(vi_1s.H, 0):
ax1 = subplot(gs2[0,0])
v0 = vi_1s.vis1(N=5, rlabels=[ .5], fontproperties=fp,)
v0[0].set_cmap(cm.gray)
tt = ax1.set_title('different number of sections', fontproperties=fp)
tt.set_position((1, 1))
#ax1.get_xaxis().set_label_coords(1, 0.9) # (x,y) from 0 to 1
fmt2(gca())
ntxt(gca(), "N=1")
subplot(gs2[0,1])
v = vi_2s.vis1(N=5, rlabels=[ .5], fontproperties=fp, )
v[0].set_cmap(cm.gray)
fmt2(gca())
ntxt(gca(), "N=2")
subplot(gs2[1,0])
v = vi_3s.vis1(N=5, rlabels=[ .5], fontproperties=fp)
v[0].set_cmap(cm.gray)
fmt2(gca())
ntxt(gca(), "N=3")
subplot(gs2[1,1])
v = vi_5s.vis1(N=5, rlabels=[ .5], fontproperties=fp)
v[0].set_cmap(cm.gray)
fmt2(gca())
ntxt(gca(), "N=5")
# --- legend axis
axl = axes([.05, .9125, .9, .075], frameon=False)
axl.set_xticks([])
axl.set_yticks([])
# red lines
axl.plot([0.1,0.3],[0, 0], '-', color=color_red)
# blue lines
axl.plot([2.1,2.3],[0, 0], '-', color=color_blue)
# black lines
axl.plot([4.1,4.3],[0, 0], 'k-')
text(0.35, 0, 'reduced linear model', va='center', fontproperties=fp)
text(2.35, 0, 'simulated SLIP', va='center', fontproperties=fp)
text(4.35, 0, 'full kinematic state', va='center', fontproperties=fp)
axl.set_xlim((-1.5,7))
# --- --- general formatting
savefig('img/fig_ev_subj{}.pdf'.format(conf.subject))
#savefig('img/fig_ev_subj{}.svg'.format(conf.subject))
pass
In [70]:
In [71]: