last edits:
This notebook continues scientific investigations from the original "AnkleSLIP" notebook.
It has been splitted to regain the general overview.
This script tries to identify a minimalistic SLIP extension (e.g., Ankle-Y-position only).
Step 0: configure notebook and define functions
Step 1: import data
Step 2: sanity checks (incl. visualizations)
Step 3: find minimal predictors
Step 3.1: compute factors and scores
Step 3.2: Compute eigenvalues of full state Floqute model and "factors-"Floquet model
Step 3.3: Test different explicit models
Step 4: create controlled SLIP
step 4.1: find periodic solution
step 4.2a: create autonomous factors-slip system
step 4.2b: create autonomous ankle-slip system
step 4.2c: create autonomous model (augmented-SLIP)
step 4.2d: create autonomous model (fullstate-SLIP)
step 4.3: compare eigenvalues of kinematics and models
Step 5: continuous CoM prediction (Floquet models)
In [1]:
%cd
%cd mmnotebooks
In [1]:
# run this to connect an ipython qtconsole to the notebook
%connect_info
In [38]:
# this cell_ID can be used to access this cell from other notebooks.
# cell_ID 0
# system libs
import os
import sys
import re
from copy import deepcopy
# shai's lib
import libshai.util as ut
# my libs
import mutils.io as mio
import mutils.misc as mi
import mutils.statistics as st
import mutils.FDatAn as fda
from mutils.io import build_dataset
import models.sliputil as su
import models.slip as sl
from models.sliputil import getControlMaps
from models.sliputil import get_auto_sys
conf = mio.saveable()
# define and import data
conf.subject = 2 #available subjects: 1,2,3,7. Other subjects have comparatively bad data quality
conf.ttype = 1 # 1: free running, 2: metronome running (data incomplete!)
conf.n_factors = 5 # 1-5 factors to predict SLIP parameters
conf.n_factors_doc = "how many (optimal) factors to select of the full kinematic state"
conf.exclude_IC_from_factors = False # exclude the initial conditions from kinematic state
# detrending options
conf.dt_window = 30
conf.dt_medfilter = False # use median filter instead of mean
# select ankle-SLIP instead of "factors without CoM state" ?
conf.select_ankle_SLIP = True
conf.cslip_forceZeroRef = True # reference values for controlled SLIP maps must be zero or not
conf.cslip_forceZeroRef_doc = "reference values for controlled SLIP maps must be zero or not"
# to compute periodic SLIP orbit: average over parameters or initial conditions?
conf.po_average_over_IC = True
conf.po_average_over_IC_doc = "average over IC's and T,ymin (alt: parameters) for reference SLIP"
conf.normalize_m = True # convert k -> k/m, dE -> dE/m, m -> 1.0
# startup
conf.startup_compute_full_maps = True # compute return maps on data loading?
conf.startup_n_full_maps = 20 # how many full-stride maps to compute for averaging
conf.startup_compute_PCA = False # compute PCA on startup?
conf.quiet = False # suppress a lot of output
conf.display()
In [2]:
# cell_ID 0.1
def test_model(kdata_r, kdata_l, kdata_full, side="r", out_of_sample=True):
"""
This function performs the prediction tests and creates the
corrsponding figures displaying the results.
:args:
in_dat (nxd array): data to predict. *NOTE* it is assumed that the first 3 columns
represent the CoM state
side (str): "r" or "l", look at right or left step. *NOTE* currently ignored.
out_of_sample: whether to use different parts of the data for regression and prediction
:returns:
fig (figure handle): a handle of the figure which will be created
"""
fig = figure(figsize=(17,7))
# identify indices of CoM in kdata_r - data
cidxr = []
for req_lbl in ['com_z', 'v_com_x', 'v_com_y']:
for idx, lbl in enumerate(kdata_r.kin_labels):
if lbl.lower() == req_lbl.lower():
cidxr.append(idx)
break
cidxr = array(cidxr)
print "indices of CoM in right submodel data:", cidxr
# identify indices of CoM in kdata_l - data
cidxl = []
for req_lbl in ['com_z', 'v_com_x', 'v_com_y']:
for idx, lbl in enumerate(kdata_l.kin_labels):
if lbl.lower() == req_lbl.lower():
cidxl.append(idx)
break
cidxl = array(cidxl)
print "indices of CoM in left submodel data:", cidxl
# make prediction for SLIP parameters
ax1 = fig.add_subplot(1,4,1)
pred1 = st.predTest(kdata_r.s_kin_r, kdata_r.s_param_r, rcond=1e-7,
out_of_sample=out_of_sample)
pred2 = st.predTest(kdata_full.s_kin_r, kdata_r.s_param_r, rcond=1e-7,
out_of_sample=out_of_sample)
pred3 = st.predTest(hstack([kdata_r.s_param_l[:-1,:], kdata_r.s_kin_r[1:, cidxr]]),
kdata_r.s_param_r[1:, ], rcond=1e-7, out_of_sample=out_of_sample) # there is
# no "data before first right apex"
b11 = ax1.boxplot(pred1, positions = arange(pred1.shape[1]) * .6 + .25)
mi.recolor(b11,'r')
b12 = ax1.boxplot(pred2, positions = arange(pred2.shape[1]) * .6 + .25)
mi.recolor(b12,'k')
b13 = ax1.boxplot(pred3, positions = arange(pred3.shape[1]) * .6 + .25)
mi.recolor(b13,'g')
ax1.set_title('prediction for SLIP parameters\nred: reduced model\nblack: ' +
'full state information\ngreen: [CoM; last SLIP params]')
ax1.set_xticks(arange(pred1.shape[1]) * .6 + .25)
ax1.set_xticklabels(['$k$', '$\\alpha$', '$L_0$', '$\\beta$', '$\\Delta E$'])
ax1.set_yticks(.1 * arange(12))
ax1.plot(ax1.get_xlim(), [1, 1], 'k--')
ax1.set_ylim(0, 1.15)
# make "self consistency" prediction -> CoM state prediction is part of this
ax2 = fig.add_subplot(1,4,2)
pred1 = st.predTest(kdata_r.s_kin_r, kdata_l.s_kin_l, rcond=1e-7,
out_of_sample=out_of_sample)
pred2 = st.predTest(kdata_full.s_kin_r, kdata_l.s_kin_l, rcond=1e-7,
out_of_sample=out_of_sample)
pred3 = st.predTest(hstack([kdata_r.s_param_l[:-1,:], kdata_r.s_kin_r[1:, cidxr]]),
kdata_l.s_kin_l[1:, cidxl], rcond=1e-7, out_of_sample=out_of_sample) # there is
# no "data before first right apex"
b21 = ax2.boxplot(pred1, positions = arange(pred1.shape[1]) * .6 + .25)
mi.recolor(b21,'r')
b23 = ax2.boxplot(pred3, positions = cidxl * .6 + .25)
mi.recolor(b23,'g')
b22 = ax2.boxplot(pred2, positions = arange(pred2.shape[1]) * .6 + .25)
mi.recolor(b22,'k')
ax2.set_title('state prediction after 1 step \n(self-consistency)\nred: reduced model\nblack: '+
'full state information\ngreen: [CoM; last SLIP params]')
ax2.set_xticks(arange(pred2.shape[1]) * .6 + .25)
ax2.set_xticklabels(kdata_l.kin_labels, rotation='vertical')
ax2.set_yticks(.1 * arange(12))
ax2.plot(ax2.get_xlim(), [1, 1], 'k--')
ax2.set_ylim(0, 1.15)
# make "self consistency" prediction -> CoM state prediction is part of this
ax3 = fig.add_subplot(1,4,3)
pred1 = st.predTest(kdata_r.s_kin_r[:-1, :], kdata_r.s_kin_r[1:, :], rcond=1e-7,
out_of_sample=out_of_sample)
pred2 = st.predTest(kdata_full.s_kin_r[:-1, :], kdata_r.s_kin_r[1:, :], rcond=1e-7,
out_of_sample=out_of_sample)
pred3 = st.predTest(hstack([kdata_r.s_param_l[:-2,:], kdata_r.s_kin_r[1:-1, cidxr]]),
kdata_l.s_kin_r[2:, cidxr], rcond=1e-7, out_of_sample=out_of_sample) # there is no
# "data before first right apex"
b31 = ax3.boxplot(pred1, positions = arange(pred1.shape[1]) * .6 + .2)
mi.recolor(b31,'r')
b33 = ax3.boxplot(pred3, positions = cidxr * .6 + .25)
mi.recolor(b33,'g')
b32 = ax3.boxplot(pred2, positions = arange(pred2.shape[1]) * .6 + .225)
mi.recolor(b32,'k')
ax3.set_title('state prediction after 1 stride\n (self-consistency test)\nred: ' +
'reduced model\nblack: full state information')
ax3.set_xticks(arange(pred2.shape[1]) * .6 + .25)
ax3.set_xticklabels(kdata_r.kin_labels, rotation='vertical')
ax3.set_yticks(.1 * arange(12))
ax3.plot(ax3.get_xlim(), [1, 1], 'k--')
ax3.set_ylim(0, 1.15)
# also: create two-stage prediction model: R->L->R instead of R->->R (directly)
# first: manually do two-stage prediction
# manually do the bootstrap
# define shortcuts
dat_r = kdata_r.s_kin_r
dat_l = kdata_l.s_kin_l
# 2 stages, reduced model
bs_vred = []
for rep in range(100):
idcs = randint(0, dat_r.shape[0] - 1, dat_r.shape[0] - 1)
if out_of_sample:
oidx = fda.otheridx(idcs, dat_r.shape[0] - 1)
else:
oidx = idcs
# regression R->L
A = dot(dat_l[idcs, :].T, pinv(dat_r[idcs, :].T, rcond=1e-7))
# regression L->R
B = dot(dat_r[idcs + 1, :].T, pinv(dat_l[idcs, :].T, rcond=1e-7))
pred = None # debug. dot(A, dat_r[oidx, :].T).T
pred_tot = reduce(dot, [B, A, dat_r[oidx, :].T]).T
vred = var(dat_r[oidx + 1, :] - pred_tot, axis=0) / var(dat_r[oidx + 1, :], axis=0)
bs_vred.append(vred)
bs_mdl = vstack(bs_vred)
# 2 stages, full model, select only which coordinates to display
# identify which components (indices) are taken in the small models:
idx_mdl = []
for idx in arange(kdata_r.s_kin_r.shape[1]):
for idxf in arange(kdata_full.s_kin_r.shape[1]):
if allclose(kdata_full.s_kin_r[:, idxf], kdata_r.s_kin_r[:, idx], atol=1e-7):
idx_mdl.append(idxf)
break
submdl_ident = len(idx_mdl) == kdata_r.s_kin_r.shape[1] # all corresponding dimensions found
if submdl_ident:
idx_mdl = array(idx_mdl)
print "corresponding dimensions of reduced and full model identified"
print "indices:", idx_mdl
dat_r = kdata_full.s_kin_r
dat_l = kdata_full.s_kin_l
bs_vred = []
for rep in range(100):
idcs = randint(0, dat_r.shape[0] - 1, dat_r.shape[0] - 1)
if out_of_sample:
oidx = fda.otheridx(idcs, dat_r.shape[0] - 1)
else:
oidx = idcs
#oidx = fda.otheridx(idcs, dat_r.shape[0] - 1)
# regression R->L
A = dot(dat_l[idcs, :].T, pinv(dat_r[idcs, :].T, rcond=1e-7))
# regression L->R
B = dot(dat_r[idcs + 1, :].T, pinv(dat_l[idcs, :].T, rcond=1e-7))
pred = None # dot(A, dat_r[oidx, :].T).T
pred_tot = reduce(dot, [B, A, dat_r[oidx, :].T]).T
vred = var(dat_r[oidx + 1, :] - pred_tot, axis=0) / var(dat_r[oidx + 1, :], axis=0)
bs_vred.append(vred)
bs_full = vstack(bs_vred)
# 2 stages, augmented SLIP model ( [CoM; last SLIP params] )
dat_r = hstack([ kdata_r.s_kin_r[1:, cidxr], kdata_r.s_param_l[:-1, :] ])
dat_l = hstack([ kdata_l.s_kin_l[1:, cidxl], kdata_l.s_param_r[1:, :] ])
bs_vred = []
for rep in range(100):
idcs = randint(0, dat_r.shape[0] - 1, dat_r.shape[0] - 1)
if out_of_sample:
oidx = fda.otheridx(idcs, dat_r.shape[0] - 1)
else:
oidx = idcs
# regression R->L
A = dot(dat_l[idcs, :].T, pinv(dat_r[idcs, :].T, rcond=1e-7))
# regression L->R
B = dot(dat_r[idcs + 1, :].T, pinv(dat_l[idcs, :].T, rcond=1e-7))
pred = None # debug. remove this line if no error occurs dot(A, dat_r[oidx, :].T).T
pred_tot = reduce(dot, [B, A, dat_r[oidx, :].T]).T
vred = var(dat_r[oidx + 1, :] - pred_tot, axis=0) / var(dat_r[oidx + 1, :], axis=0)
bs_vred.append(vred)
bs_augslip = vstack(bs_vred)[:, :3]
# display results
ax4 = fig.add_subplot(1,4,4)
b41 = ax4.boxplot(bs_mdl, positions = arange(bs_mdl.shape[1]) * .6 + .25)
mi.recolor(b41,'r')
b43 = ax4.boxplot(bs_augslip, positions = cidxr * .6 + .25)
mi.recolor(b43,'g')
if submdl_ident:
b42 = ax4.boxplot(bs_full[:, idx_mdl], positions = arange(len(idx_mdl)) * .6 + .275)
mi.recolor(b42,'k')
ax4.set_title('state prediction after 1 stride using\n two steps for prediction\nred: ' +
'reduced model\nblack: full state information')
ax4.set_xticks(arange(bs_mdl.shape[1]) * .6 + .25)
ax4.set_xticklabels(kdata_r.kin_labels, rotation='vertical')
ax4.set_yticks(.1 * arange(12))
ax4.plot(ax4.get_xlim(), [1, 1], 'k--')
ax4.set_ylim(0, 1.15)
return fig
def reord_idx(labels):
"""
returns indices for re-ordering dataset data such that first three columns are
CoM state (SLIP state)
:args:
labels [list]: a list of labels that describe the ordering of the columns. it is scanned for
com_z (-> idx0) , v_com_x (-> idx2), v_com_y (idx1)
:returns:
idcs (array): an array of the size of labels that starts with [idx0, idx1, idx2] and
contains all numbers from 0 to len(labels)-1.
This can be used for indexing as in the example.
:example:
re_ordered_data = original_data[:, reord_idx(kin_labels)]
"""
idx0 = [idx for idx, val in enumerate(labels) if val.lower() == 'com_z'][0]
idx1 = [idx for idx, val in enumerate(labels) if val.lower() == 'v_com_y'][0]
idx2 = [idx for idx, val in enumerate(labels) if val.lower() == 'v_com_x'][0]
ridx = list(set(range(len(labels))) - set([idx0, idx1, idx2]))
idcs = hstack([idx0, idx1, idx2, ridx])
return idcs
def getPhase(state, system, nstride=10):
"""
computes the phase of the given SLIP system at the state "state".
:args:
state [n-by-1 array]: the initial state of the system
system [autonomous SLIP]: the dynamical system (SLIP model), obtained e.g. by
models.slip.get_auto_sys
*NOTE* the system must be stable!
nstride [int]: number of strides to be simulated before convergence is assumed
:returns:
(phi0, T) [float, float]: the phase correspoding to the initial condition and the average
cycle duration
"""
# run #nstride strides
sim_t_d = []
sim_s_d = []
fs_d = array(state).copy()
for rep in range(nstride):
fs_d, (sim_t, sim_states) = system(fs_d)
oldtime = 0
if len(sim_t_d) > 0:
oldtime = sim_t_d[-1][-1]
sim_t_d.append(sim_t[:-1] + oldtime)
sim_s_d.append(sim_states[:-1,:])
sim_s_d.append(sim_states[-1, :][newaxis, :])
sim_t_d.append(sim_t[-1] + oldtime)
sim_s_d = vstack(sim_s_d)
sim_t_d = hstack(sim_t_d)
# assume that solution is now periodic -> simulate "reference" stride
fs_f, (sim_t_ref, sim_states_ref) = system(fs_d)
# now, compute phase (after simulation is completed)
# at the end, phase is 0 (or 2pi)
# total amount of phase elaped is: (time / sim_t_ref) * 2pi
phi_total = sim_t_d[-1] / sim_t_ref[-1] * (2*pi)
# final phase is: nstride * 2pi
final_phase = nstride * 2*pi
initial_phase = final_phase - phi_total
# disturbed phase is no:
return initial_phase, sim_t_ref[-1]
def get_evScore(rmap, data, dims=3):
"""
This function returns the score of the dataset on the largest eigenvalues of the return map.
:args:
rmap (d-by-d arrary): the return maps from which the eigenvalues will be computed
data (n-by-d array): data from which the map was computed
dims (int, optionally): number of largest ev's to consider
:returns:
sc (dims-by-n array): scores of the data projected to the highest eigenvalues
"""
ev, ew = eig(rmap)
idcs = argsort(abs(ev))[::-1][:dims]
return dot(ew[:, idcs].T, data.T).T
In [1]:
!pwd
In [3]:
# cell_ID 1
# load kinematic data
ws1 = mio.saveable()
#ws1.k = mio.KinData(data_dir = '/home/moritz/data/2011-mmcl_mat/')
ws1.k = mio.KinData(data_dir = 'data/2011-mmcl_mat/')
ws1.k_doc = "generic kinematic data access object (empty selection)"
ws1.k.load(conf.subject, conf.ttype)
ws1.k_doc = "KinData object: s:" + str(conf.subject) + ", t:" + str(conf.ttype) + ", selection: ankles"
ws1.k.selection = ['r_anl_y - com_y', 'l_anl_y - com_y', 'com_z']
if conf.normalize_m:
ws1.SlipData = [mio.normalize_m(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]
else:
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]
ws1.SlipData_doc = "SlipData for s:" + str(conf.subject) + ", t:" + str(conf.ttype)
if not conf.quiet:
ws1.display()
# 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_wrl_x - com_x',
'r_elb_x - com_x', 'l_elb_x - com_x',
'l_wrl_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',
'r_elb_y - com_y',
'r_acr_y - com_y', 'l_acr_y - com_y',
'l_elb_y - com_y', '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',
'r_wrl_z - com_z',
'r_acr_z - com_z', 'l_acr_z - com_z',
'l_wrl_z - com_z']
ws1.full_markers_doc = 'representative selection of markers for quasi-full state space information'
if not conf.quiet:
print "len of 'ws1.full_markers':", len(ws1.full_markers), " (all: 84)"
ws1.k.selection = ws1.full_markers #all_markers
ws1.dataset_full = build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
ws1.dataset_full_doc = ("reference 'dataset' for s:" + str(ws1.k.subject) + " t:" + str(ws1.k.ttype) +
" with 'full' selection")
# "normalized" dataset: velocities scaled by small factor
ws1.dataset_full.n_kin_r = ws1.dataset_full.all_kin_r.copy()
ws1.dataset_full.n_kin_l = ws1.dataset_full.all_kin_l.copy()
ws1.dataset_full.n_kin_r[:, len(ws1.k.selection)-2:] /= 11.
ws1.dataset_full.n_kin_l[:, len(ws1.k.selection)-2:] /= 11.
ws1.dataset_full.n_kin_r -= mean(ws1.dataset_full.n_kin_r, axis=0)
ws1.dataset_full.n_kin_l -= mean(ws1.dataset_full.n_kin_l, axis=0)
ws1.dataset_full.n_kin_r = fda.dt_movingavg(ws1.dataset_full.n_kin_r, conf.dt_window, conf.dt_medfilter)
ws1.dataset_full.n_kin_l = fda.dt_movingavg(ws1.dataset_full.n_kin_l, conf.dt_window, conf.dt_medfilter)
ws1.dataset_full.n_kin_r_doc = "detrended 'all_kin_r'; velocity has been scaled by factor 1/11"
ws1.dataset_full.n_kin_l_doc = "detrended 'all_kin_l'; velocity has been scaled by factor 1/11"
if conf.startup_compute_full_maps:
_, 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)]
ws1.maps_full_doc = "return maps for 'full' selection model (2-section based)"
# --- --- extended consistency check
if not conf.quiet:
print "\n running extended SLIP - EXP consistency check"
print '=============================================='
mass = median(ws1.dataset_full.masses)
scales_l = std(ws1.dataset_full.all_IC_l, axis=0)
scales_r = std(ws1.dataset_full.all_IC_r, axis=0)
print "right steps:"
i_ic = ws1.dataset_full.all_IC_r[30::200, :]
i_fs = ws1.dataset_full.all_IC_l[30::200, :]
i_p = ws1.dataset_full.all_param_r[30::200, :]
for IC, P, fs in zip(i_ic, i_p, i_fs):
sim_fs = su.finalState(IC, P, {'m' : mass, 'g': -9.81})
print " difference (in stds): {:+.3f}, {:+.3f}, {:+.3f}".format( *(sim_fs - fs) / scales_l )
print "left steps:"
i_ic = ws1.dataset_full.all_IC_l[30::200, :]
i_fs = ws1.dataset_full.all_IC_r[30+1::200, :]
i_p = ws1.dataset_full.all_param_l[30::200, :]
for IC, P, fs in zip(i_ic, i_p, i_fs):
sim_fs = su.finalState(IC, P, {'m' : mass, 'g': -9.81})
print " difference (in stds): {:+.3f}, {:+.3f}, {:+.3f}".format( *(sim_fs - fs) / scales_r)
# --- --- END consistency check
if conf.startup_compute_full_maps and not conf.quiet:
vi = ut.VisEig(127)
for A in ws1.maps_full:
vi.add(eig(A)[0])
figure(figsize=(6,6))
vi.vis1()
title('EV distribution of "full" marker model')
gca().set_xlim(-1.1, 1.1)
gca().set_ylim(-1.1, 1.1)
gca().set_ylabel('imaginary part')
gca().set_xlabel('real part')
if conf.startup_compute_PCA and not conf.quiet:
d1d = ws1.k.make_1D(nps=50,)[:, 100:]
d1m = mean(d1d, axis=0)
d1d -= d1m
u,s,v = svd(d1d, full_matrices=False)
figure(figsize=(16,9))
for pc in range(5):
subplot(5,1,pc+1)
if pc==0:
title('scores on the principal components of the motion')
plot(u.T[pc, :] / std(u.T[pc, :]),'k.')
grid('on')
ylabel('pc #%i' % (pc+1))
pass # suppress output function calls
In [ ]:
#--------- visualize largest eigenvectors
m = fda.meanMat(ws1.maps_full)
ev, ew = eig(m)
order = argsort(abs(ev))[::-1]
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 [6]:
p0 = [x['phi2'] for x in ws1.k.raw_dat]
# find phases "closest" to apex
pr = ws1.dataset_full.all_phases_r
pl = ws1.dataset_full.all_phases_l
all_idx_r = []
all_idx_l = []
for ap, pa in zip(p0, pr):
idx_r = []
for p in pa:
idx_r.append(argmin(abs(p - ap)))
print "rep. done"
all_idx_r.append(hstack(idx_r))
for ap, pa in zip(p0, pl):
idx_l = []
for p in pa:
idx_l.append(argmin(abs(p - ap)))
print "rep. done"
all_idx_l.append(hstack(idx_l))
dpr = [diff(ap.squeeze())[idcx]*250 for ap, idcx in zip(p0, all_idx_r)] # 250: sampling rate
dpl = [diff(ap.squeeze())[idcx]*250 for ap, idcx in zip(p0, all_idx_l)]
figure()
plot(hstack(dpr), 'b.', label='right')
plot(hstack(dpl), 'g.', label='left')
title('subject {}: changes in phase velocity (at apex) over time'.format(conf.subject))
ylabel('phase velocity at apex [rad s$^{-1}$]')
xlabel('stride')
ylim(7,11)
legend(loc='best')
# plot phase at apex
figure(),
plot(mod(hstack(pr), 2*pi), 'b.', label='right')
plot(mod(hstack(pl), 2*pi) - pi, 'g.', label='left ($-\pi$)')
title('subject {}: changes in average phase at apex over time'.format(conf.subject))
ylabel('phase at apex [rad]')
xlabel('stride')
legend(loc='best')
pass
In [ ]:
phi1 = vstack(ws1.dataset.all_phi_r)
phi2 = vstack(ws1.dataset.all_phi_l)
In [7]:
if True: # skip cell? True=Execute
if conf.normalize_m:
SlipData = [mio.normalize_m(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]
else:
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]
tr = []
tl = []
for s in SlipData:
ridx = find(mod(s.phases, 2*pi) < pi)
lidx = find(mod(s.phases, 2*pi) > pi)
tr.append(array(s.T_exp)[ridx])
tl.append(array(s.T_exp)[lidx])
tr = hstack(tr)
tl = hstack(tl)
figure(figsize=(8,6))
plot(tr, 'b.', label='right')
plot(tl, 'g.', label='left')
xlabel('stride number')
ylabel('step time [s]')
legend(loc='best')
title('subject {}: changes in step time (frequency) over time'.format(conf.subject))
savefig('tmp/tstep_s{}.png'.format(conf.subject), res=300)
pass
In [ ]:
# select visualization dataset
nps=100 # frames per stride
if not 'k' in locals() or not type(k) == mio.KinData:
k = mio.KinData(data_dir = 'data/2011-mmcl_mat/')
markers = [
['r_mtv', 'r_hee', 'r_anl', 'r_kne', 'r_trc', 'r_sia', 'r_sip', 'sacr'],
['l_mtv', 'l_hee', 'l_anl', 'l_kne', 'l_trc', 'l_sia', 'l_sip', 'sacr'],
['sacr', 'cvii', 'r_hea', 'l_hea','cvii'],
['r_wrl', 'r_elb', 'r_acr', 'l_acr', 'l_elb', 'l_wrl'],
]
#
ymarkers = []
zmarkers = []
for elem in markers:
ymarkers.append([x + '_y - com_y' for x in elem])
zmarkers.append([x + '_z ' for x in elem])
all_markers = [ 'com_x', 'com_y', 'com_z']
for y in ymarkers:
all_markers.extend(y)
for z in zmarkers:
all_markers.extend(z)
k.selection = all_markers
k.load(conf.subject, conf.ttype)
if not 'ws1' in locals():
ws1 = mio.saveable()
ws1.k = mio.KinData(data_dir = 'data/2011-mmcl_mat/')
ws1.k_doc = "generic kinematic data access object (empty selection)"
ws1.k.load(conf.subject, conf.ttype)
if conf.normalize_m:
ws1.SlipData = [mio.normalize_m(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]
else:
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]
dataset = build_dataset(k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
d1d_vis = k.make_1D(nps=nps, phases_list=dataset.all_phases_r )
print "test dataset loaded for subj", conf.subject
In [ ]:
if not 'dataset' in locals():
d1d = k.make_1D(nps=50,)[:, 100:]
print "warning: different phases and number of strides used!"
else:
d1d = k.make_1D(nps=50, phases_list=dataset.all_phases_r)[:, 100:]
#d1d = ws1.k.make_1D(nps=50,)[:, 100:]
d1m = mean(d1d, axis=0)
d1d -= d1m
u,s,v = svd(d1d, full_matrices=False)
figure(figsize=(16,9))
for pc in range(5):
subplot(5,1,pc+1)
plot(u.T[pc, :] / std(u.T[pc, :]),'k.')
grid('on')
ylabel('pc #%i' % (pc+1))
#savefig('tmp/pc-s%i.png' % conf.subject)
In [ ]:
# detrended principal components (!= pc's of detrended data)
pcdt = fda.dt_movingavg(u, conf.dt_window, conf.dt_medfilter)
figure(figsize=(12,8))
thresh = 2.5 #3.5 # times std
ba = []
for pc in range(25):
badidx = find(abs(pcdt.T[pc, :]) / std(pcdt.T[pc, :]) > thresh)
ba.append(badidx)
ba.append(badidx - 1)
goodidx = fda.otheridx( hstack(ba), pcdt.shape[0])
for pc in range(5):
subplot(5,1,pc+1)
plot(u.T[pc, :] / std(u.T[pc, :]),'k.')
grid('on')
ylabel('pc #%i' % (pc+1))
plot(goodidx, pcdt.T[pc, goodidx] / std(u.T[pc, :]),'r.')
if pc == 0:
title('"good" indices marked in red')
In [ ]:
# make prediction test for all markers
d1d_dt = fda.dt_movingavg(d1d, conf.dt_window, conf.dt_medfilter)
res1 = st.predTest(d1d_dt[:-1, ::50], d1d_dt[1:, ::50], rcond=1e-3)
res2 = st.predTest(d1d_dt[goodidx[:-1], ::50], d1d_dt[goodidx[:-1]+1, ::50], rcond=1e-3, out_of_sample=True)
figure(figsize=(24,12))
b1 = boxplot(vstack(res1))
b2 = boxplot(vstack(res2))
mi.recolor(b1, 'k')
mi.recolor(b2, 'r')
gca().set_xticks(range(112))
labels = k.selection[2:] + ['v_' + elem for elem in k.selection]
gca().set_xticklabels(labels, rotation=90)
gca().set_ylim(0.5,1.5)
title('predictability for non-extreme data (red) vs all data (black)')
show()
print "done; shapes:", d1d.shape[0], d1d_dt.shape[0]
In [ ]:
# preparations for the visualization
do_visualize = True
mean_gc = mean(d1d_vis, axis=0)
len_sct = [len(x) for x in markers]
tlen = sum(len_sct)
In [ ]:
# visualize selected strides
if do_visualize:
s0=528 # first stride
se=543 # last stride
fig = figure('anim', figsize=(16,9))
fig.clf()
ax = fig.add_axes([0,0,1,1])
aframe=0 #counter
# loop over these:
frame=0
step=s0
for step in range(s0, se+1):
for frame in range(nps):
ax.cla()
cnt = 3 # skip first three elements - they are CoM
for sct in len_sct:
x = []
y = []
xm = []
ym = []
sgm = []
for mk in range(sct):
idxx0 = (mk + cnt)*nps + frame
idxy0 = (mk + cnt + tlen)*nps + frame
x.append(d1d_vis[step, idxx0])
y.append(d1d_vis[step, idxy0])
xm.append(mean_gc[idxx0])
ym.append(mean_gc[idxy0])
# plot "treadmill"
fill([-1.5, 1.5, 1.5, -1.5], [0, 0, -.1, -.1], color='#b3b3df')
#print sgm
plot(xm, ym, 'o-', color='#30df70', linewidth=2)
plot(x, y, 'ko-', alpha=.9)
plot(0, mean_gc[nps*2 + frame], 'go', markersize=10)
plot(0, d1d_vis[step, nps*2 + frame], 'ko', markersize=8, alpha=.75)
cnt += sct
text(-1.5, 1.5, ('subject: %i\nstride: %i\nframe: %i/%i' % (conf.subject,
step, frame+1, nps,)))
text(-1.2, 1.35, 'treadmill top view')
pos = [-1.5, 1.05] # position of "treadmill top view" plot
scale = .5 # scale of "treadmill top view"
fill(array([0, 1.4, 1.4, 0])*scale + pos[0], array([.54, .54, 0, 0])*scale
+ pos[1], 'w')
plot(pos[0] + scale*d1d_vis[step, frame + nps],
pos[1] + scale*d1d_vis[step, frame], 'ro')
ax.set_xlim(-1.6*1.3, 1.6*1.3)
ax.set_ylim(-.1*1.3, 1.7*1.3)
grid('on')
#axis('equal')
fig.savefig('tmp/frame_%04i.png' % aframe)
aframe +=1
# conversion into video on linux systems
!avconv -y -r 15 -f image2 -i tmp/frame_%04d.png -vb 600k tmp/animstrides.ogg
!rm tmp/frame*.png
print "done! tmp/animstrides.ogg written"
In [ ]:
# visualize two "average" gait cycles
mean_gc1 = mean(d1d_vis[:300, :], axis=0)
mean_gc2 = mean(d1d_vis[-300:, :], axis=0)
if do_visualize:
fig = figure('anim', figsize=(16,9))
fig.clf()
ax = fig.add_axes([0,0,1,1])
aframe=0 #counter
# loop over these:
frame=0
if True: # just for indentation
for frame in range(nps):
ax.cla()
cnt = 3 # skip first three elements - they are CoM
for sct in len_sct:
x = []
y = []
xm = []
ym = []
sgm = []
for mk in range(sct):
idxx0 = (mk + cnt)*nps + frame
idxy0 = (mk + cnt + tlen)*nps + frame
xm.append(mean_gc2[idxx0])
ym.append(mean_gc2[idxy0])
x.append(mean_gc1[idxx0])
y.append(mean_gc1[idxy0])
# plot "treadmill"
fill([-1.5, 1.5, 1.5, -1.5], [0, 0, -.1, -.1], color='#b3b3df')
#print sgm
plot(xm, ym, 'o-', color='#30df70', linewidth=2)
plot(x, y, 'ko-', alpha=.9)
cnt += sct
text(-1.5, 1.5, ('subject: %i\nframe: %i/%i' % (conf.subject, frame+1, nps,)))
ax.set_xlim(-1.6*1.3, 1.6*1.3)
ax.set_ylim(-.1*1.3, 1.7*1.3)
grid('on')
fig.savefig('tmp/frame2avg_%04i.png' % aframe)
aframe +=1
# conversion into video on linux systems
!avconv -y -r 15 -f image2 -i tmp/frame2avg_%04d.png -vb 600k tmp/avgs.ogg
!rm tmp/frame2avg_*.png
print "done! tmp/avgs.ogg written"
In [ ]:
# this is only interesting for subject #3, plotting pc2 vs pc3
if conf.subject==3:
figure()
plot(u.T[2, :], u.T[3, :], '.')
In [ ]:
# visualize selected PC
if do_visualize:
pc = 4 # select which principal component of the motion to visualize (0-based)
# d1m: mean gait cycle
mov_pc = d1m + s[pc]*v[pc, :]
fig = figure('animpc', figsize=(16,9))
fig.clf()
ax = fig.add_axes([0,0,1,1])
aframe=0 # counter
# loop over these:
frame=0
if True: # dummy, for indentation
for frame in range(50):
ax.cla()
cnt = 1 # skip first element - this is vertical CoM position
for sct in len_sct:
x = []
y = []
xm = []
ym = []
sgm = []
for mk in range(sct):
idxx0 = (mk + cnt)*50 + frame
idxy0 = (mk + cnt + tlen)*50 + frame
x.append(mov_pc[idxx0])
y.append(mov_pc[idxy0])
xm.append(d1m[idxx0])
ym.append(d1m[idxy0])
plot(xm, ym, 'ko-', linewidth=2)
plot(x, y, 'o-', alpha=.9, color='#df4070')
cnt += sct
text(-1.5, 1.5, ('subject: %i\n pc: %i\n (nps: 50)' % (conf.subject, pc) ))
ax.set_xlim(-1.6*1.5, 1.6*1.5)
ax.set_ylim(-.1*1.5, 1.7*1.5)
grid('on')
#axis('equal')
fig.savefig('tmp/framepc_%04i.png' % aframe)
aframe +=1
# conversion into video on linux systems
!avconv -y -r 15 -f image2 -i tmp/framepc_%04d.png -vb 600k tmp/pc.ogg
!rm tmp/framepc_*.png
print "done! tmp/pc.ogg written"
In [ ]:
d1d_dt = fda.dt_movingavg(d1d, conf.dt_window, conf.dt_medfilter)
res = st.predTest(d1d_dt[:-1, ::50], d1d_dt[1:, ::50], rcond=1e-3)
figure()
boxplot(vstack(res))
gca().set_xticks(range(112))
labels = k.selection[2:] + ['v_' + elem for elem in k.selection]
gca().set_xticklabels(labels, rotation=90)
print "done"
In [ ]:
# compute return maps for same data
_, maps1, idcs = fda.fitData(d1d_dt[:-1, ::50], d1d_dt[:-1, 25::50], nps=1, rcond=1e-5, nrep=100)
_, maps2, idcs = fda.fitData(d1d_dt[:-1, 25::50], d1d_dt[1:, ::50], nps=1, rcond=1e-5, nrep=100)
maps = [dot(x,y) for x,y in zip(maps2, maps1)]
vi = ut.VisEig(127, False)
for A in maps:
vi.add(eig(A)[0])
fig = figure('sect. 2.3 ev of data', figsize=(9,9))
fig.clf()
vi.vis1()
title('eigenvalues of the 1st return map')
for ev in eig(fda.meanMat(maps))[0]:
plot(ev.real, ev.imag, 'ro')
gca().set_xlim(-1,1)
gca().set_ylim(-1,1)
In [ ]:
# plot histogram of angle distribution between random vectors in R^n (here: n=112)
# technically, the average scalar product can be computed analytically, when we assume
# that the variance is spread equally randomly over all axes
if False:
angs = []
for rep in range(1000):
a,b = randn(112), randn(112)
a /= norm(a)
b /= norm(b)
angs.append(abs(dot(a,b)))
figure()
hist(angs)
xlabel('cosine (a,b)')
ylabel('frequency')
title('angle between two random vectors in $R^{112}$')
In [ ]:
# plot magnitude of projection of principal components at sct #0 on first eigenvectors
ev, ew = eig(fda.meanMat(maps))
order_idx = argsort(abs(ev))[::-1]
ew_s = ew[:, order_idx] # sorted eigenvectors, according to EV magnitude
allval = []
for n in range(5):
pc_direction = v[n, ::50].copy()
pc_direction /= norm(pc_direction)
vals = []
for evnr in range(20):
vals.append(abs(dot(pc_direction, ew_s[:, evnr])))
allval.append(array(vals))
figure()
allval = vstack(allval).T
for n in range(allval.shape[1]):
plot(allval[:, n], 'd', label='PC #%i' % (n+1))
title('magnitude of dot product of first 5 PCs on "largest" 20 EVs')
gca().set_xticks(arange(21))
xlabel('ordinal of eigenvector (ordered by EV magnitude)')
ylabel(r'$||\vec{PC} \cdot \vec{EV}||$')
legend()
In [4]:
## prepare data
# cell_ID 3
# *NOTE* the first three labels have to be "com_x", "com_y", "com_z". If you edit this, make sure to edit
# the cell where data are further processed...
ws1.k.selection = ws1.full_markers
ws1.k_doc = "KinData object: s:" + str(ws1.k.subject) + ", t:" + str(ws1.k.ttype) + ", selection: 'full'"
ws1.dataset = ws1.dataset_full
ws1.dataset_doc = ws1.dataset_full_doc
if not conf.quiet:
# another consistency check
stepnr = 698 # select any stride between 0 and ~1800 (subject dependend)
mass = mean([d.mass for d in ws1.SlipData])
# su.finalState actually performs a one-step integration of SLIP
print "consistency check:"
print "------------------"
print "simres [right step]:", su.finalState(ws1.dataset.all_IC_r[stepnr, :],
ws1.dataset.all_param_r[stepnr, :],
{'m' : mean(ws1.dataset.masses), 'g' : -9.81})
print "subsequent left apex:", ws1.dataset.all_IC_l[stepnr, :]
print "==========="
print "simres [left step]:", su.finalState(ws1.dataset.all_IC_l[stepnr, :],
ws1.dataset.all_param_l[stepnr, :],
{'m' : mean(ws1.dataset.masses), 'g' : -9.81})
print "subsequent right apex:", ws1.dataset.all_IC_r[stepnr + 1, :]
step 3: find minimal predictors
to content
In [5]:
# cell_ID 3.1
# Here, the main predictors ("directions") are computed. This is computationally quite fast
if conf.exclude_IC_from_factors:
# indices that do *not* represent CoM (only relative positions, indicated by "minus" sign)
sel_idx = array([nr for nr, label in enumerate(ws1.dataset.kin_labels) if '-' in label])
ws1.facs_r_doc = ("directions of 'factors' for a right step excluding CoM, s:" +
str(ws1.k.subject) + ", t:" + str(ws1.k.ttype))
ws1.facs_l_doc = ("directions of 'factors' for a left step excluding CoM, s:" +
str(ws1.k.subject) + ", t:" + str(ws1.k.ttype))
else:
sel_idx = arange(len(ws1.dataset.kin_labels))
ws1.facs_r_doc = ("directions of 'factors' for a right step including CoM, s:" +
str(ws1.k.subject) + ", t:" + str(ws1.k.ttype))
ws1.facs_l_doc = ("directions of 'factors' for a left step including CoM, s:" +
str(ws1.k.subject) + ", t:" + str(ws1.k.ttype))
ws1.facs_labels = [ws1.dataset.kin_labels[x] for x in sel_idx]
ws1.facs_labels_doc = "kinematic labels for the components of the 'factors'"
ws1.facs_r = st.find_factors(ws1.dataset.s_kin_r[:, sel_idx].T, ws1.dataset.s_param_r.T, k=conf.n_factors)
ws1.facs_l = st.find_factors(ws1.dataset.s_kin_l[:, sel_idx].T, ws1.dataset.s_param_l.T, k=conf.n_factors)
ws1.fscore_r = dot(ws1.facs_r.T, ws1.dataset.s_kin_r[:, sel_idx].T).T
ws1.fscore_r_doc = "scores of data when projected on the factors (right)"
ws1.fscore_l = dot(ws1.facs_l.T, ws1.dataset.s_kin_l[:, sel_idx].T).T
ws1.fscore_l_doc = "scores of data when projected on the factors (left)"
ws1.raw_factor_space_r = ws1.dataset.s_kin_r[:, sel_idx]
ws1.raw_factor_space_l = ws1.dataset.s_kin_l[:, sel_idx]
ws1.raw_factor_space_l_doc = "part of state space used for computing the 'factors'"
ws1.raw_factor_space_l_doc = "part of state space used for computing the 'factors'"
# ---------------------visualize
if not conf.quiet:
figure(figsize=(20,3))
pcolor(ws1.facs_r.T)
gca().set_xticks(arange(len(ws1.facs_labels)) + .5)
gca().set_xticklabels(ws1.facs_labels, rotation=90)
grid('on')
title('magnitude of the factors')
gca().set_yticks(arange(ws1.facs_r.shape[1]) + .5)
gca().set_yticklabels(arange(ws1.facs_r.shape[1]) + 1)
gca().set_ylabel('factor #')
pass
(moved to separate notebook)
step 3: find minimal predictors
to content
approach:
problem: The map A * B should be identical with the (single) map of a full state model
In [6]:
factor_mdl = 1 # 1: factors, CoM; 2: factors, CoM, phase at apex 3: only factors
# cell_ID 3.2
if factor_mdl == 1: # reduced model: IC, factors
ws1.reddat_r = hstack([ws1.fscore_r, fda.dt_movingavg(ws1.dataset.all_IC_r,
conf.dt_window, conf.dt_medfilter) /
std(fda.dt_movingavg(ws1.dataset.all_IC_r, conf.dt_window,
conf.dt_medfilter), axis=0)])
ws1.reddat_l = hstack([ws1.fscore_l, fda.dt_movingavg(ws1.dataset.all_IC_l,
conf.dt_window, conf.dt_medfilter) /
std(fda.dt_movingavg(ws1.dataset.all_IC_l, conf.dt_window,
conf.dt_medfilter), axis=0)])
elif factor_mdl ==2: # as above + phases at apex mod 2pi
pr = mod(hstack(ws1.dataset_full.all_phases_r), 2*pi)
pl = mod(hstack(ws1.dataset_full.all_phases_l), 2*pi)
pr = fda.dt_movingavg(pr[:, newaxis], conf.dt_window, conf.dt_medfilter)
pl = fda.dt_movingavg(pl[:, newaxis], conf.dt_window, conf.dt_medfilter)
ws1.reddat_r = hstack([ws1.fscore_r, fda.dt_movingavg(ws1.dataset.all_IC_r,
conf.dt_window, conf.dt_medfilter) /
std(fda.dt_movingavg(ws1.dataset.all_IC_r,
conf.dt_window, conf.dt_medfilter), axis=0),
pr])
ws1.reddat_l = hstack([ws1.fscore_l, fda.dt_movingavg(ws1.dataset.all_IC_l,
conf.dt_window, conf.dt_medfilter) /
std(fda.dt_movingavg(ws1.dataset.all_IC_l,
conf.dt_window, conf.dt_medfilter), axis=0),
pl])
elif factor_mdl == 3: # factors only
ws1.reddat_r = ws1.fscore_r
ws1.reddat_l = ws1.fscore_l
else:
raise ValueError("Invalid factor-model selected!")
if not conf.quiet:
print "states computed"
In [7]:
ws1.reddat_r_doc = "reduced data model: [factors; CoM (z, vx, vy; normalized)].T (right)"
ws1.reddat_l_doc = "reduced data model: [factors; CoM (z, vx, vy; normalized)].T (right)"
_, maps_1, _ = fda.fitData(ws1.reddat_r, ws1.reddat_l, nps=1, nrep=200, rcond=1e-7)
_, maps_2, _ = fda.fitData(ws1.reddat_l[:-1, :], ws1.reddat_r[1:, :], nps=1, nrep=200,
rcond=1e-7, )
maps_red = [dot(m2, m1) for m1, m2 in zip(maps_1, maps_2)]
# visualize
vi_full = ut.VisEig(127, False)
for A in ws1.maps_full:
vi_full.add(eig(A)[0])
vi_red = ut.VisEig(127, False)
for A in maps_red:
vi_red.add(eig(A)[0])
figure(figsize=(16,8))
subplot(1,2,1)
vi_full.vis1()
gca().set_xlim(-1,1)
gca().set_ylim(-1,1)
xlabel('real part')
ylabel('imaginary part')
title('subj. {}: eigenvalues of full dim. Floquet model'.format(conf.subject))
subplot(1,2,2)
vi_red.vis1()
gca().set_xlim(-1,1)
gca().set_ylim(-1,1)
xlabel('real part')
ylabel('imaginary part')
title('subj. {}: eigenvalues of CoM + {} factors'.format(conf.subject, conf.n_factors))
pass # suppress output of last command
step 3: find minimal predictors
step 3.3.1: model 1 (x,y,vx,vy) both ankles
step 3.3.2: model 2 : swing leg ankle
step 3.3.3: model 3 : ankle + tilt indicator
step 3.3.4: model 4 : CoM + factors
to content
things to predict:
NOTE There are two different approaches to test the 1-stride autonomy:
NOTE The coordinate system here is:
x - lateral
y - anterior
z - vertical
step 3: find minimal predictors
step 3.7: test different explicit models
to content
In [ ]:
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'r_anl_y - com_y',
'r_anl_z - com_z', 'r_anl_x - com_x',] #'r_sia_y - l_sia_y'] # add hip rotation indicator?
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',] #'r_sia_y - l_sia_y'] # add hip rotation indicator?
tmp_dsl = 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 a 'ankles' selection")
d1 = fda.dt_movingavg(tmp_dsr.all_kin_r, conf.dt_window, conf.dt_medfilter)
d2 = fda.dt_movingavg(tmp_dsl.all_kin_l, conf.dt_window, conf.dt_medfilter)
_, m1, _ = fda.fitData(d1,d2,1, nrep=200, rcond=1e-6)
_, m2, _ = fda.fitData(d1,d2,1, nrep=200, rcond=1e-6)
maps = [reduce(dot, [mb, ma]) for ma, mb in zip(m1, m2)]
vi = ut.VisEig(127, False)
for A in maps:
vi.add(eig(A)[0])
figure(figsize=(7,7))
vi.vis1()
gca().set_xlim(-1,1)
gca().set_ylim(-1,1)
xlabel('real part')
ylabel('imaginary part')
title('subj. {}: eigenvalues of (unilat.) ankle SLIP'.format(conf.subject))
pass # suppress output
In [ ]:
# here: use only information from leg that will touchdown next
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'r_anl_y - com_y',
'r_anl_z - com_z', 'r_anl_x - com_x',] #'r_sia_y - l_sia_y'] # add hip rotation indicator?
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',] #'r_sia_y - l_sia_y'] # add hip rotation indicator?
tmp_dsl = 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 a 'ankles' selection")
fig = test_model(tmp_dsr, tmp_dsl, ws1.dataset, out_of_sample=True)
In [ ]:
# here: use only a/p information from leg that will touchdown next
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'r_anl_y - com_y',
] #'r_sia_y - l_sia_y'] # add hip rotation indicator?
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',
] #'r_sia_y - l_sia_y'] # add hip rotation indicator?
tmp_dsl = 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 a 'ankles' selection")
fig = test_model(tmp_dsr, tmp_dsl, ws1.dataset, out_of_sample=True)
The "tilt indicator" is actually only a set of two additional markers: neck and sacral marker
step 3: find minimal predictors
step 3.7: test different explicit models
to content
In [ ]:
# here: use only information from leg that will touchdown next
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'r_anl_y - com_y',
'r_anl_z - com_z', 'r_anl_x - com_x', 'cvii_y - com_y', 'sacr_y - com_y']
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', 'cvii_y - com_y', 'sacr_y - com_y']
tmp_dsl = 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 a 'ankles' selection")
#fig = test_model(tmp_dsr, tmp_dsl, ws1.dataset_full)
fig = test_model(tmp_dsr, tmp_dsl, ws1.dataset, out_of_sample=True)
In this current version, the prediciton test is performed manually
step 3: find minimal predictors
step 3.7: test different explicit models
to content
In [11]:
out_of_sample = True # predict out of sample?
r1a = st.predTest(ws1.reddat_r, ws1.reddat_l, out_of_sample=out_of_sample)
r1b = st.predTest(ws1.dataset.s_kin_r, ws1.reddat_l, out_of_sample=out_of_sample)
r2a = st.predTest(ws1.dataset.s_kin_r, ws1.dataset.s_param_r, out_of_sample=out_of_sample)
r2b = st.predTest(ws1.reddat_r, ws1.dataset.s_param_r, out_of_sample=out_of_sample)
r2c = st.predTest(hstack([ws1.dataset.all_IC_r[1:,:], ws1.dataset.s_param_l[:-1,:]]),
ws1.dataset.s_param_r[1:,:], out_of_sample=out_of_sample)
figure(figsize=(18,6))
subplot(1,2,1)
b1 = boxplot(vstack(r1a))
mi.recolor(b1,'r')
b2 = boxplot(vstack(r1b))
mi.recolor(b2,'k')
gca().set_xticks(arange(conf.n_factors + 3) + 1)
gca().set_xticklabels(['factor {}'.format(x) for x in arange(conf.n_factors) + 1] + ['com_z', 'v_com_x', 'v_com_y'],
rotation=90)
gca().set_ylabel('rel. rem. Variance')
title('\n'.join(['subject {}: model self-prediction (consistency)'.format(conf.subject),
'black: full state', 'red: CoM {} factors'.format(conf.n_factors)]))
ylim(0,1)
subplot(1,2,2)
b1 = boxplot(vstack(r2a))
mi.recolor(b1,'k')
b2 = boxplot(vstack(r2b))
mi.recolor(b2,'r')
b2 = boxplot(vstack(r2c))
mi.recolor(b2,'g')
gca().set_xticks(arange(5)+1)
gca().set_xticklabels([r'k', r'$\alpha$', r'$L_0$', r'$\beta$', r'$\Delta E$'])
gca().set_ylabel('rel. rem. Variance')
title('\n'.join(['subject {}: SLIP parameter prediction'.format(conf.subject), 'green: augmented SLIP',
'black: full state', 'red: CoM {} factors'.format(conf.n_factors)]))
ylim(0,1)
pass #suppress output
step 4.1: find periodic solution
step 4.2a: create autonomous factors-slip system
step 4.2b: create autonomous ankle-slip system
step 4.2c: create autonomous model (augmented-SLIP)
step 4.2d: create autonomous model (fullstate-SLIP)
Step 4.3: compare eigenvalues of kinematics and models
Step 3 (for factors and scores on factors)
The closed loop model reads as follows:
Step 1 ("right"):
state = [IC; Factors] where IC = initial CoM state at apex
params = P1 * state where P is a regression from data
new_state = [ode(IC, params); Q1 * [IC; Factors]] where Q is a regressed map from data
Step 2 ("left"):
state = [IC; Factors] where IC = initial CoM state at apex
params = P2 * state where P is a regression from data
new_state = [ode(IC, params); Q2 * [IC; Factors]] where Q is a regressed map from data
step 4: create controlled SLIP
to content
Approach 1: find periodic solution that corresponds to the average parameters, and verify that it's close to the average initial conditions.
Approach 2: find periodic solution that corresponds to the average initial conditions (and time and ymin), and verify that the corresponding parameters are close to the average parameters
to select approach 1 or 2, go to the "init" block at the top of this notebook
In [8]:
# cell_ID 4
ws1.cslip = mio.saveable()
ws1.cslip_doc = "'subdirectory' for the controlled SLIP that corresponds to the Floquet data"
ws1.cslip.mean_pr = mean(ws1.dataset_full.all_param_r, axis=0)
ws1.cslip.mean_pr_doc = "average parameter for right leg"
ws1.cslip.mean_pl = mean(ws1.dataset_full.all_param_l, axis=0)
ws1.cslip.mean_pl_doc = "average parameter for left leg"
ws1.cslip.mean_ICr = mean(ws1.dataset_full.all_IC_r, axis=0)
ws1.cslip.mean_ICr_doc = "mean CoM apex state right"
ws1.cslip.mean_ICl = mean(ws1.dataset_full.all_IC_l, axis=0)
ws1.cslip.mean_ICl_doc = "mean CoM apex state left"
ws1.cslip.mass = mean(ws1.dataset_full.masses)
ws1.cslip.mass_doc = "mass of SLIP (mean mass of subject over all datasets)"
if not conf.po_average_over_IC: # average parameters, compute corresponding periodic solution
raise NotImplementedError, ("ERROR: periodic solution for average parameters " +
"not tested for new SLIP implementation!")
ws1.cslip.mean_pl[4] = -ws1.cslip.mean_pr[4] # set total energy change to exactly zero.
# Note: typically, the difference is low, roughly ~0.01 J
ws1.cslip.ICp_r, ws1.cslip.ICp_l = su.getPeriodicOrbit_p(ws1.cslip.mean_pr, ws1.cslip.mean_pl,
aug_dict={'m': ws1.cslip.mass, 'g' : -9.81}, startState=ws1.cslip.mean_ICr)
ws1.cslip.ICp_r_doc = "initial conditions (r) for periodic solution (mean-parameter periodic)"
ws1.cslip.ICp_l_doc = "initial conditions (l) for periodic solution (mean-parameter periodic)"
ws1.cslip.Pp_r = ws1.cslip.mean_pr
ws1.cslip.Pp_l = ws1.cslip.mean_pl
ws1.cslip.Pp_r_doc = "parameters (r) for periodic solution (mean-parameter periodic)"
ws1.cslip.Pp_l_doc = "parameters (l) for periodic solution (mean-parameter periodic)"
else: # average initial conditions, compute corresponding periodic solution
#[ws1.cslip.ICp_r, Pp_r, dE_r], [ws1.cslip.ICp_l, Pp_l, dE_l] = su.getPeriodicOrbit(ws1.dataset_full.all_IC_r,
# vstack(ws1.dataset_full.TR), vstack(ws1.dataset_full.yminR),
# ws1.dataset_full.all_IC_l, vstack(ws1.dataset_full.TL),
# vstack(ws1.dataset_full.yminL), baseParams ={'m': ws1.cslip.mass, 'g' : -9.81},
# startParams=mean(vstack(ws1.dataset_full.all_param_r), axis=0)[:4])
m = ws1.cslip.mass
ICr = mean(ws1.dataset_full.all_IC_r, axis=0)
FSr = mean(ws1.dataset_full.all_IC_l, axis=0)
yminr = mean(vstack(ws1.dataset_full.yminR))
Tr = mean(vstack(ws1.dataset_full.TR))
ICl = mean(ws1.dataset_full.all_IC_l, axis=0)
FSl = mean(ws1.dataset_full.all_IC_r, axis=0)
yminl = mean(vstack(ws1.dataset_full.yminL))
Tl = mean(vstack(ws1.dataset_full.TL))
mpr = mean(ws1.dataset_full.all_param_r, axis=0)
pr, pl = su.getPeriodicOrbit2(ICr, Tr, yminr, ICl, Tl, yminl, m, startParams=mpr)
ws1.cslip.ICp_r = ICr
ws1.cslip.ICp_l = ICl
ws1.cslip.ICp_r_doc = "initial conditions (r) for periodic solution (mean-IC periodic)"
ws1.cslip.ICp_l_doc = "initial conditions (l) for periodic solution (mean-IC periodic)"
# change last element to be energy input - this is consistent with the format used in the rest of the code
#Pp_r = array(Pp_r)
#Pp_l = array(Pp_l)
#Pp_r[4] = dE_r
#ws1.cslip.Pp_r = Pp_r[:5].copy()
#Pp_l[4] = dE_l
#Pp_l = Pp_l[:5]
#ws1.cslip.Pp_l = Pp_l[:5].copy()
ws1.cslip.Pp_r = pr.copy()
ws1.cslip.Pp_l = pl.copy()
ws1.cslip.Pp_r_doc = "parameters (r) for periodic solution (mean-IC periodic)"
ws1.cslip.Pp_l_doc = "parameters (l) for periodic solution (mean-IC periodic)"
#verify periodicity
if not conf.quiet:
sr1 = su.finalState(ws1.cslip.ICp_r, ws1.cslip.Pp_r, addDict = {'m' : ws1.cslip.mass, 'g' : -9.81})
sl1 = su.finalState(ws1.cslip.ICp_l, ws1.cslip.Pp_l, addDict = {'m' : ws1.cslip.mass, 'g' : -9.81})
print "\n S A N I T Y C H E C K:\n"
print "difference between identified periodic orbit and simulation results:"
print "left apex:", sr1 - ws1.cslip.ICp_l
print "right apex:", sl1 - ws1.cslip.ICp_r
print '\n===========================================================\n'
print "Deviation from periodic orbit to average apex state, in units of standard deviations:"
print " [height, horizontal speed, lateral speed]"
print "left step: ", (ws1.cslip.ICp_l - mean(ws1.dataset_full.all_IC_l, axis=0)) / std(
ws1.dataset_full.all_IC_l, axis=0)
print "right step:", (ws1.cslip.ICp_r - mean(ws1.dataset_full.all_IC_r, axis=0)) / std(
ws1.dataset_full.all_IC_r, axis=0)
print '\n===========================================================\n'
print "relative deviation from periodic parameters to mean parameters,"
print "in units of standard deviation of parameters"
print "[k, alpha, l0, beta, deltaE]"
print (ws1.cslip.Pp_l - mean(ws1.dataset_full.all_param_l, axis=0)) / std(ws1.dataset_full.all_param_l, axis=0)
print (ws1.cslip.Pp_r - mean(ws1.dataset_full.all_param_r, axis=0)) / std(ws1.dataset_full.all_param_r, axis=0)
In [13]:
# create factor-model (new code)
# cell_ID 4.2a
states_r = hstack([ws1.dataset.all_IC_r, ws1.fscore_r ])
states_l = hstack([ws1.dataset.all_IC_l, ws1.fscore_l ])
# indices for regression (bypasses bootstrap)
if not 'indices_' in locals():
indices_ = None
if not conf.po_average_over_IC:
raise NotImplementedError("FACTPR-SLIP for (average-param) periodic solution not yet implemented!")
# create controlled ankle slip manually (see "fac_slip" above)
# the getControlMaps - function internally computes the reference solution based upon average ICs!
cmaps, nmaps, refs = getControlMaps(states_r, states_l, ws1.dataset, conf, indices=indices_)
ws1.cslip.fac_slip_ctrl_r = cmaps[0].copy()
ws1.cslip.fac_slip_ctrl_l = cmaps[1].copy()
ws1.cslip.fac_slip_refstate_r = hstack([ws1.cslip.ICp_r, zeros(states_r.shape[1] - 3)])
ws1.cslip.fac_slip_refstate_l = hstack([ws1.cslip.ICp_l, zeros(states_r.shape[1] - 3)])
ws1.cslip.fac_slip_ctrl_r_doc = "control map: full state (minus reference) -> slip parameter (right)"
ws1.cslip.fac_slip_ctrl_l_doc = "control map: full state (minus reference) -> slip parameter (left)"
ws1.cslip.fac_slip_refstate_r_doc = "reference state for factor slip (right)"
ws1.cslip.fac_slip_refstate_l_doc = "reference state for factor slip (left)"
ws1.cslip.fac_slip_dprop_r = nmaps[0].copy()
ws1.cslip.fac_slip_dprop_l = nmaps[1].copy()
ws1.cslip.fac_slip_dprop_r_doc = ("discrete 1-step propagator of additional states (right) " +
"(here: slip parameter themselves)")
ws1.cslip.fac_slip_dprop_l_doc = ("discrete 1-step propagator of additional states (left) " +
"(here: slip parameter themselves)")
ws1.cslip.fac_slip = get_auto_sys(ws1.cslip.Pp_r, ws1.cslip.Pp_l, ws1.cslip.fac_slip_refstate_r,
ws1.cslip.fac_slip_refstate_l, ws1.cslip.fac_slip_ctrl_r, ws1.cslip.fac_slip_ctrl_l,
ws1.cslip.fac_slip_dprop_r, ws1.cslip.fac_slip_dprop_l, {'m': ws1.cslip.mass, 'g':-9.81})
ws1.cslip.fac_slip_doc = "stride function of the autonomous controlled SLIP model ('factor slip')"
In [9]:
# create factor-model
# old code
# cell_ID old_4.2a
if False: # do not run old code
# indices for regression (bypasses bootstrap)
if not 'indices_' in locals():
indices_ = None
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])
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)
ws1.cslip.fac_slip_ctrl_r = fda.meanMat(all_ctrl_r)
ws1.cslip.fac_slip_ctrl_l = fda.meanMat(all_ctrl_l)
ws1.cslip.fac_slip_ctrl_r_doc = "mapping from full model state (minus reference) to SLIP parameters (right)"
ws1.cslip.fac_slip_ctrl_l_doc = "mapping from full model state (minus reference) to SLIP parameters (left)"
# 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)
ws1.cslip.fac_slip_dprop_r = fda.meanMat(all_facprop_r)
ws1.cslip.fac_slip_dprop_l = fda.meanMat(all_facprop_l)
ws1.cslip.fac_slip_dprop_r_doc = "discrete 1-step propagator of additional states (right) (here: factors)"
ws1.cslip.fac_slip_dprop_l_doc = "discrete 1-step propagator of additional states (left) (here: factors)"
# set this to zero to force "zero" as reference for the fscores!
allow_nonzero_ref = 0. if conf.cslip_forceZeroRef else 1.
ws1.cslip.fac_slip_refstate_r = hstack(
[ws1.cslip.ICp_r, allow_nonzero_ref * mean(ws1.fscore_r, axis=0)]) # the latter should be zero anyway
ws1.cslip.fac_slip_refstate_l = hstack(
[ws1.cslip.ICp_l, allow_nonzero_ref * mean(ws1.fscore_l, axis=0)]) # the latter should be zero anyway
ws1.cslip.fac_slip_refstate_r_doc = "reference state (periodic): [IC; scores on factos] (right)"
ws1.cslip.fac_slip_refstate_l_doc = "reference state (periodic): [IC; scores on factos] (left)"
# create autonomous "factor-based" SLIP model
ws1.cslip.fac_slip = get_auto_sys(ws1.cslip.Pp_r, ws1.cslip.Pp_l,
ws1.cslip.fac_slip_refstate_r, ws1.cslip.fac_slip_refstate_l, ws1.cslip.fac_slip_ctrl_r,
ws1.cslip.fac_slip_ctrl_l, ws1.cslip.fac_slip_dprop_r, ws1.cslip.fac_slip_dprop_l,
{'m': ws1.cslip.mass, 'g':-9.81})
ws1.cslip.fac_slip_doc = "stride function of the controlled SLIP model ('factors')"
In [10]:
std(dt_fstate_r, axis=0)
Out[10]:
In [11]:
ws1.cslip.ICp_l
Out[11]:
In [12]:
getControlMaps(ws1.dataset_full)
In [15]:
IC = ws1.cslip.fac_slip_refstate_r
t = ws1.cslip.fac_slip
def multirun(fun, IC, n):
res = array(IC)
for rep in range(n):
res = fun(res)
return res
print multirun(t, IC, 20) - IC
In [16]:
cmaps, nmaps, refs = getControlMaps(dt_fstate_r, dt_fstate_l, dataset_r, conf)
ws1.cslip.ankle_slip_ctrl_r = cmaps[0].copy()
ws1.cslip.ankle_slip_ctrl_l = cmaps[1].copy()
ws1.cslip.ankle_slip_refstate_r = hstack([ws1.cslip.ICp_r, zeros(states_r.shape[1] - 3)])
ws1.cslip.ankle_slip_refstate_l = hstack([ws1.cslip.ICp_l, zeros(states_r.shape[1] - 3)])
ws1.cslip.ankle_slip_ctrl_r_doc = "control map: full state (minus reference) -> slip parameter (right)"
ws1.cslip.ankle_slip_ctrl_l_doc = "control map: full state (minus reference) -> slip parameter (left)"
ws1.cslip.ankle_slip_refstate_r_doc = "reference state for ankle slip (right)"
ws1.cslip.ankle_slip_refstate_l_doc = "reference state for ankle slip (left)"
ws1.cslip.ankle_slip_dprop_r = nmaps[0].copy()
ws1.cslip.ankle_slip_dprop_l = nmaps[1].copy()
ws1.cslip.ankle_slip_dprop_r_doc = "discrete 1-step propagator of additional states (right) (here: ankles)"
ws1.cslip.ankle_slip_dprop_l_doc = "discrete 1-step propagator of additional states (left) (here: ankles)"
In [42]:
std(ws1.fscore_r, axis=0)
Out[42]:
In [17]:
# cell_ID 4.2b
# indices for regression (bypasses bootstrap)
if not 'indices_' in locals():
indices_ = None
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)]
# --- add belt speed
vbelt = .5 *( mean(ws1.dataset_full.all_IC_r[:,1]) - mean(states_r[:,1]) +
mean(ws1.dataset_full.all_IC_l[:,1]) - mean(states_l[:,1]))
states_r[:,1] += vbelt
states_l[:,1] += vbelt
if not conf.po_average_over_IC:
raise NotImplementedError("Ankle-SLIP for (average-param) periodic solution not yet implemented!")
# create controlled ankle slip manually (see "fac_slip" above)
# the getControlMaps - function internally computes the reference solution based upon average ICs!
cmaps, nmaps, refs = getControlMaps(states_r, states_l, dataset_r, conf, indices=indices_)
ws1.cslip.ankle_slip_ctrl_r = cmaps[0].copy()
ws1.cslip.ankle_slip_ctrl_l = cmaps[1].copy()
ws1.cslip.ankle_slip_refstate_r = hstack([ws1.cslip.ICp_r, zeros(states_r.shape[1] - 3)])
ws1.cslip.ankle_slip_refstate_l = hstack([ws1.cslip.ICp_l, zeros(states_r.shape[1] - 3)])
ws1.cslip.ankle_slip_ctrl_r_doc = "control map: full state (minus reference) -> slip parameter (right)"
ws1.cslip.ankle_slip_ctrl_l_doc = "control map: full state (minus reference) -> slip parameter (left)"
ws1.cslip.ankle_slip_refstate_r_doc = "reference state for ankle slip (right)"
ws1.cslip.ankle_slip_refstate_l_doc = "reference state for ankle slip (left)"
ws1.cslip.ankle_slip_dprop_r = nmaps[0].copy()
ws1.cslip.ankle_slip_dprop_l = nmaps[1].copy()
ws1.cslip.ankle_slip_dprop_r_doc = "discrete 1-step propagator of additional states (right) (here: ankles)"
ws1.cslip.ankle_slip_dprop_l_doc = "discrete 1-step propagator of additional states (left) (here: ankles)"
ws1.cslip.ankle_slip = get_auto_sys(ws1.cslip.Pp_r, ws1.cslip.Pp_l, ws1.cslip.ankle_slip_refstate_r,
ws1.cslip.ankle_slip_refstate_l, ws1.cslip.ankle_slip_ctrl_r, ws1.cslip.ankle_slip_ctrl_l,
ws1.cslip.ankle_slip_dprop_r, ws1.cslip.ankle_slip_dprop_l, {'m': ws1.cslip.mass, 'g':-9.81})
ws1.cslip.ankle_slip_doc = "stride function of the autonomous controlled SLIP model ('ankle slip')"
In [18]:
# strictly speaking, the first element is invalid ...
# cell_ID 4.2c
states_r = hstack([ws1.dataset.all_IC_r, roll(ws1.dataset.s_param_l, 1, axis=0) ])
states_l = hstack([ws1.dataset.all_IC_l, ws1.dataset.s_param_r ])
if not 'indices_' in locals():
indices_ = None
if not conf.po_average_over_IC:
raise NotImplementedError("Augmented-SLIP for (average-param) periodic solution not yet implemented!")
# create controlled ankle slip manually (see "fac_slip" above)
# the getControlMaps - function internally computes the reference solution based upon average ICs!
cmaps, nmaps, refs = getControlMaps(states_r, states_l, ws1.dataset, conf, indices=indices_)
ws1.cslip.aug_slip_ctrl_r = cmaps[0].copy()
ws1.cslip.aug_slip_ctrl_l = cmaps[1].copy()
ws1.cslip.aug_slip_refstate_r = hstack([ws1.cslip.ICp_r, zeros(states_r.shape[1] - 3)])
ws1.cslip.aug_slip_refstate_l = hstack([ws1.cslip.ICp_l, zeros(states_r.shape[1] - 3)])
ws1.cslip.aug_slip_ctrl_r_doc = "control map: full state (minus reference) -> slip parameter (right)"
ws1.cslip.aug_slip_ctrl_l_doc = "control map: full state (minus reference) -> slip parameter (left)"
ws1.cslip.aug_slip_refstate_r_doc = "reference state for augmented slip (right)"
ws1.cslip.aug_slip_refstate_l_doc = "reference state for augmented slip (left)"
ws1.cslip.aug_slip_dprop_r = nmaps[0].copy()
ws1.cslip.aug_slip_dprop_l = nmaps[1].copy()
ws1.cslip.aug_slip_dprop_r_doc = ("discrete 1-step propagator of additional states (right) " +
"(here: slip parameter themselves)")
ws1.cslip.aug_slip_dprop_l_doc = ("discrete 1-step propagator of additional states (left) " +
"(here: slip parameter themselves)")
ws1.cslip.aug_slip = get_auto_sys(ws1.cslip.Pp_r, ws1.cslip.Pp_l, ws1.cslip.aug_slip_refstate_r,
ws1.cslip.aug_slip_refstate_l, ws1.cslip.aug_slip_ctrl_r, ws1.cslip.aug_slip_ctrl_l,
ws1.cslip.aug_slip_dprop_r, ws1.cslip.aug_slip_dprop_l, {'m': ws1.cslip.mass, 'g':-9.81})
ws1.cslip.aug_slip_doc = "stride function of the autonomous controlled SLIP model ('augmented slip')"
In [ ]:
# cell_ID 4.2d
# indices for regression (bypasses bootstrap)
if not 'indices_' in locals():
indices_ = None
ws1.k.selection = ws1.full_markers
dataset_full = 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 full markers"
states_r = dataset_full.all_kin_r[:, reord_idx(dataset_full.kin_labels)]
states_l = dataset_full.all_kin_l[:, reord_idx(dataset_full.kin_labels)]
# --- add belt speed
vbelt = .5 *( mean(ws1.dataset_full.all_IC_r[:,1]) - mean(states_r[:,1]) +
mean(ws1.dataset_full.all_IC_l[:,1]) - mean(states_l[:,1]))
states_r[:, 1] += vbelt
states_l[:, 1] += vbelt
if not conf.po_average_over_IC:
raise NotImplementedError("Ankle-SLIP for (average-param) periodic solution not yet implemented!")
# create controlled ankle slip manually (see "fac_slip" above)
# the getControlMaps - function internally computes the reference solution based upon average ICs!
cmaps, nmaps, refs = getControlMaps(states_r, states_l, dataset_r, conf, indices=indices_)
ws1.cslip.full_slip_ctrl_r = cmaps[0].copy()
ws1.cslip.full_slip_ctrl_l = cmaps[1].copy()
ws1.cslip.full_slip_refstate_r = hstack([ws1.cslip.ICp_r, zeros(states_r.shape[1] - 3)])
ws1.cslip.full_slip_refstate_l = hstack([ws1.cslip.ICp_l, zeros(states_r.shape[1] - 3)])
ws1.cslip.full_slip_ctrl_r_doc = "control map: full state (minus reference) -> slip parameter (right)"
ws1.cslip.full_slip_ctrl_l_doc = "control map: full state (minus reference) -> slip parameter (left)"
ws1.cslip.full_slip_refstate_r_doc = "reference state for fullstate slip (right)"
ws1.cslip.full_slip_refstate_l_doc = "reference state for fullstate slip (left)"
ws1.cslip.full_slip_dprop_r = nmaps[0].copy()
ws1.cslip.full_slip_dprop_l = nmaps[1].copy()
ws1.cslip.full_slip_dprop_r_doc = "discrete 1-step propagator of additional states (right) (here: fs)"
ws1.cslip.full_slip_dprop_l_doc = "discrete 1-step propagator of additional states (left) (here: fs)"
ws1.cslip.full_slip = get_auto_sys(ws1.cslip.Pp_r, ws1.cslip.Pp_l, ws1.cslip.full_slip_refstate_r,
ws1.cslip.full_slip_refstate_l, ws1.cslip.full_slip_ctrl_r, ws1.cslip.full_slip_ctrl_l,
ws1.cslip.full_slip_dprop_r, ws1.cslip.full_slip_dprop_l, {'m': ws1.cslip.mass, 'g':-9.81})
ws1.cslip.ankle_slip_doc = "stride function of the autonomous controlled SLIP model ('fullstate slip')"
step 4: create controlled SLIP
to content
also, visualize eigenvalues and eigenvalues of the "reduced" discrete model
In [19]:
# skip floquet model (re-)computation?
quick_run = True # skip floquet model eigenvalue computation if already exists in workspace
# cell_ID 4.3
# eigenvalues for factor slip
J1 = []
h = .001
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))
J1 = vstack(J1).T
# eigenvalues for 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))
J2 = vstack(J2).T
# eigenvalues for augmented slip
J3 = []
h = .003
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]))
J3 = vstack(J3).T
# for comparision: compute eigenvalues of (reduced) Floquet model
_, maps_1, _ = fda.fitData(ws1.reddat_r, ws1.reddat_l, nps=1, nrep=200, rcond=1e-7)
_, maps_2, _ = fda.fitData(ws1.reddat_l[:-1, :], ws1.reddat_r[1:, :], nps=1, nrep=200, rcond=1e-7, )
maps_red = [dot(m1, m2) for m1, m2 in zip(maps_1, maps_2)]
if not(quick_run and 'fmaps' in locals()):
# for comparison: compute eigenvalues of full state with multiple sections (try n=7!)
nps = 5
ws1.k.selection = ws1.full_markers
ws1.k_doc = ("KinData object for s:" + str(ws1.k.subject) +
" t:" + str(ws1.k.ttype) + " with 'all markers' selection")
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:]
# compute multiple-section map
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
_, maps1, _ = fda.fitData(d1d[:skip,rep::nps], d1d[start:,erep::nps], nps=1, rcond=1e-7, nrep=80)
pt_maps.append(maps1)
fmaps = [reduce(dot, ptm) for ptm in zip(*pt_maps[::-1])]
print "done!"
In [20]:
# visualize
# cell_ID 4.3.1
vi_red = ut.VisEig(127, False)
for A in maps_red:
vi_red.add(eig(A)[0])
vi_full = ut.VisEig(127, False)
for A in ws1.maps_full:
vi_full.add(eig(A)[0])
vi_full_multi = ut.VisEig(127, False)
for A in fmaps:
vi_full_multi.add(eig(A)[0])
figure(figsize=(9,9))
# factor
lbl = 'factor-SLIP'
for ev in eig(J1)[0]:
plot(ev.real, ev.imag, 'cd', label=lbl)
lbl = None
# ankle
lbl = 'ankle-SLIP'
for ev in eig(J2)[0]:
plot(ev.real, ev.imag, 'kd', label=lbl)
lbl=None
# augmented
lbl = 'augmented SLIP'
for ev in eig(J3)[0]:
plot(ev.real, ev.imag, 'ro', label=lbl)
lbl=None
axis('equal')
v = vi_red.vis1(N=3)
v[0].set_cmap(cm.spring)
vf = vi_full.vis1(N=5, rlabels=[.25, .5])
vf[0].set_cmap(cm.gray)
vfm = vi_full_multi.vis1(N=3, rlabels=[])
vfm[0].set_cmap(cm.autumn)
vfm[0].set_cmap(cm.winter)
#xlim(-.6,.6)
#ylim(-.6,.6)
xlim(-1.1, 1.1)
ylim(-1.1, 1.1)
legend()
title('\n'.join(['eigenvalues of reduced model vs. factors-controlled SLIP',
'(gray: full state Floquet model) Subject: {}'.format(conf.subject)]))
print "USE LESS CONTOUR LEVEL -> REMOVE CLUTTER"
print "SEE QUESTION ON TOP"
Requirements: Note the next cell automatically (re-)intializes the notebook
Cell 2.01 (here)
Cell 3.1 and 3.2
In [1]:
%%capture capt
# *NOTE* this cell prints two figures when it finished; other output is suppressed by the cell magic command
# --- run required cells automatically? :) (notebook must be stored for this!)
import mutils.io as mio
import mutils.misc as mi
ws5subject = 3
if not 'ws5' in locals():
ws5 = mio.saveable()
ws5.store_figs = True
# load basic config
mi.run_nbcells('AnkleSLIP - find minimal model- Version 3.ipynb', ['0'])
conf.subject = ws5subject
conf.startup_n_full_maps = 10
conf.startup_compute_PCA
mi.run_nbcells('AnkleSLIP - find minimal model- Version 3.ipynb',
['0.1','1', '3', '3.1', '3.2','4', '4.2a', '4.2b', '4.2c'])
In [2]:
# create continuous model
# exclude first apex value because "augmented SLIP" is not defined for this
# these files contain the content of the corresponding cells.
# set this to "True" to re-run everything
ws5.skip_mdls = [4,5] # models *not* to recompute! (1-6; can be empty)
ws5.subj = conf.subject
ws5.out_of_sample = True
ws5.out_of_sample_doc = 'perform out of sample prediction (for CoM)'
ws5.nboot = 100
ws5.nboot_doc = 'number of bootstrap repetitions for predictions'
ws5.use_cached_fullmdl = False # use cached file or not
ws5.nps = 50
print "building models ..."
ws1.k.selection = ['com_x', 'com_y', 'com_z']
odat_c = ws1.k.make_1D(nps=ws5.nps, phases_list=ws1.dataset_full.all_phases_r)[:, 2*ws5.nps:]
ws5.odat = fda.dt_movingavg(odat_c, conf.dt_window, conf.dt_medfilter)[1:, :]
ws5.odat_doc = 'com z, vx, vy, vz with ' + str(ws5.nps) + ' nps'
pr = mod(hstack(ws1.dataset_full.all_phases_r), 2*pi)
pr = fda.dt_movingavg(pr[:, newaxis], conf.dt_window, conf.dt_medfilter)
# create floquet models
if not 1 in ws5.skip_mdls:
ws5.idat1 = ws1.dataset_full.n_kin_r[1:, :]
ws5.idat1_doc = 'model 1: full dataset at (right) apex'
if not 2 in ws5.skip_mdls:
ws5.idat2 = ws1.reddat_r[1:, :]
addstr = 'including IC' if not conf.exclude_IC_from_factors else 'excluding IC'
ws5.idat2_doc = 'model 2: CoM + factors ' + addstr
if not 3 in ws5.skip_mdls:
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'r_anl_y - com_y',
'r_anl_z - com_z', 'r_anl_x - com_x',] #'r_sia_y - l_sia_y'] # add hip rotation indicator?
tmp_dsr = build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
ws5.idat3 = fda.dt_movingavg(tmp_dsr.all_kin_r, conf.dt_window, conf.dt_medfilter)[1:, :]
ws5.idat3_doc = 'model 3: CoM + ankles '
if not 4 in ws5.skip_mdls:
ws5.idat4 = fda.dt_movingavg(hstack([ws1.dataset_full.all_IC_r[1:, :],
pr[1:, :], # phase mod 2pi at apex
hstack(dpr)[1:, newaxis], # phase velocity at apex
]), conf.dt_window, conf.dt_medfilter)
ws5.idat4_doc = 'model 4: CoM + phi + vphi'
if not 5 in ws5.skip_mdls:
ws5.idat5 = fda.dt_movingavg(hstack([ws1.dataset_full.all_IC_r[1:, :], ws1.dataset_full.s_param_l[:-1, :],
pr[1:, :]]), conf.dt_window, conf.dt_medfilter)
ws5.idat5_doc = 'model 5: augmented SLIP (exact apex state + phase info)'
if not 6 in ws5.skip_mdls:
ws5.idat6 = fda.dt_movingavg(hstack([ws1.dataset_full.all_IC_r[1:, :], ws1.dataset_full.s_param_l[:-1, :]]),
conf.dt_window, conf.dt_medfilter)
ws5.idat6_doc = 'model 6: augmented SLIP'
if not 7 in ws5.skip_mdls:
ws1.k.selection = ['com_x', 'com_y', 'com_z', 'r_anl_y - com_y',
'r_anl_z - com_z', 'r_anl_x - com_x', 'cvii_y - sacr_y']
tmp_dsr = build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
ws5.idat7 = fda.dt_movingavg(tmp_dsr.all_kin_r, conf.dt_window, conf.dt_medfilter)[1:, :]
ws5.idat7_doc = 'model 7: CoM + ankles '
print "performing prediction tests"
print "*" * ws5.nps
print "model 1:"
if 1 in ws5.skip_mdls:
print "skipped"
else:
res = None
if ws5.use_cached_fullmdl:
try:
res = mio.mload('tmp/vred_s{}fm.list'.format(conf.subject))
print "loaded mdl from file (median detrended!)"
except IOError:
pass
if res == None: # not loaded; either not found or cache disabled
res = []
for rep in range(ws5.nps):
res.append(vstack(st.predTest(ws5.idat1, ws5.odat[:,rep::ws5.nps], out_of_sample=ws5.out_of_sample,
nboot=ws5.nboot)))
sys.stdout.write('.')
print " done"
mio.msave('tmp/vred_s{}fm.list'.format(conf.subject), res)
ws5.pred_mdl1 = res
ws5.pred_mdl1_doc = ''.join(['list of rem. var. of CoM using mdl1 for prediction; for each of ',
str(ws5.nps), ' phases of gait cyc.'])
print "model 2:"
if 2 in ws5.skip_mdls:
print "skipped"
else:
res = []
for rep in range(ws5.nps):
res.append(vstack(st.predTest(ws5.idat2, ws5.odat[:,rep::ws5.nps], out_of_sample=ws5.out_of_sample,
nboot=ws5.nboot)))
sys.stdout.write('.')
print " done"
ws5.pred_mdl2 = res
ws5.pred_mdl2_doc = ''.join(['list of rem. var. of CoM using mdl2 for prediction; for each of ',
str(ws5.nps), ' phases of gait cyc.'])
print "model 3:"
if 3 in ws5.skip_mdls:
print "skipped"
else:
res = []
for rep in range(ws5.nps):
res.append(vstack(st.predTest(ws5.idat3, ws5.odat[:,rep::ws5.nps], out_of_sample=ws5.out_of_sample,
nboot=ws5.nboot)))
sys.stdout.write('.')
print " done"
ws5.pred_mdl3 = res
ws5.pred_mdl3_doc = ''.join(['list of rem. var. of CoM using mdl3 for prediction; for each of ',
str(ws5.nps), ' phases of gait cyc.'])
print "model 4:"
if 4 in ws5.skip_mdls:
print "skipped"
else:
res = []
for rep in range(ws5.nps):
res.append(vstack(st.predTest(ws5.idat4, ws5.odat[:,rep::ws5.nps], out_of_sample=ws5.out_of_sample,
nboot=ws5.nboot)))
sys.stdout.write('.')
print " done"
ws5.pred_mdl4 = res
ws5.pred_mdl4_doc = ''.join(['list of rem. var. of CoM using mdl4 for prediction; for each of ',
str(ws5.nps),' phases of gait cyc.' ])
print "model 5:"
if 5 in ws5.skip_mdls:
print "skipped"
else:
res = []
for rep in range(ws5.nps):
res.append(vstack(st.predTest(ws5.idat5, ws5.odat[:,rep::ws5.nps], out_of_sample=ws5.out_of_sample,
nboot=ws5.nboot)))
sys.stdout.write('.')
print " done"
ws5.pred_mdl5 = res
ws5.pred_mdl5_doc = ''.join(['list of rem. var. of CoM using mdl5 for prediction; for each of ',
str(ws5.nps), ' phases of gait cyc.'])
print "model 6:"
if 6 in ws5.skip_mdls:
print "skipped"
else:
res = []
for rep in range(ws5.nps):
res.append(vstack(st.predTest(ws5.idat6, ws5.odat[:,rep::ws5.nps], out_of_sample=ws5.out_of_sample,
nboot=ws5.nboot)))
sys.stdout.write('.')
print " done"
ws5.pred_mdl6 = res
ws5.pred_mdl6_doc = ''.join(['list of rem. var. of CoM using mdl5 for prediction; for each of ' ,
str(ws5.nps), ' phases of gait cyc.'])
print "model 7:"
if 7 in ws5.skip_mdls:
print "skipped"
else:
res = []
for rep in range(ws5.nps):
res.append(vstack(st.predTest(ws5.idat7, ws5.odat[:,rep::ws5.nps], out_of_sample=ws5.out_of_sample,
nboot=ws5.nboot)))
sys.stdout.write('.')
print " done"
ws5.pred_mdl7 = res
ws5.pred_mdl7_doc = ''.join(['list of rem. var. of CoM using mdl7 for prediction; for each of ',
str(ws5.nps),' phases of gait cyc.' ])
In [3]:
# --- visualize
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)
means = []
stds = []
for mdl in arange(7) + 1:
if mdl not in ws5.skip_mdls:
code = 'vstack([mean(x, axis=0) for x in ws5.pred_mdl{}])'.format(mdl)
means.append(eval(code))
code = 'vstack([std(x, axis=0) for x in ws5.pred_mdl{}])'.format(mdl)
stds.append(eval(code))
fig = figure(figsize=(6.83,2))
phase = linspace(0,2*pi, 50, endpoint=False)
titles = ['CoM z'.format(conf.subject),
'CoM vx'.format(conf.subject),
'CoM vy'.format(conf.subject),]
colors_ = ['#006767','k', '#009900', 'r','c']
fmt_m = {'linestyle' : '-', 'marker' : '', 'linewidth' : 2}
fmt_s = {'linestyle' : '--', 'linewidth' : 1}
mdl_lbl = ['full state', 'factor SLIP', 'ankle SLIP', 'augmented SLIP', 'ankle* SLIP']
for dim in range(3):
subplot(1, 3, dim+1)
for nr, (m, s) in enumerate(zip(means, stds)):
plot(phase, m[:, dim], color=colors_[nr], **fmt_m)
plot(phase, m[:, dim] + s[:, dim], color=colors_[nr], **fmt_s)
#plot(phase, m[:, dim] + s[:, dim], styles_[nr])
title(titles[dim], fontproperties=fp)
gca().set_xticks([0, 3.14, 6.28])
gca().set_xticklabels(['0', r'$\pi$', r'$2 \pi$'], fontproperties=fp)
gca().set_yticks(arange(6) * .2)
ylim(0, 1.1)
plot([0,6.28], [1, 1], 'k--')
plot([pi, pi],[0, 1], 'k--')
grid('on')
if dim == 1:
xlabel('phase [rad]', fontproperties=fp)
if dim == 0:
gca().set_yticklabels(arange(6) * .2, fontproperties=fp)
ylabel('relative remaining variance', fontproperties=fp)
else:
gca().set_yticklabels([])
lax = axes([.025,.9,.95, .1], frameon=False)
lax.text(.75, 0, 'predictions for subject {}:'.format(conf.subject), va='center', fontproperties=fp)
xpos = [2.15, 2.9, 3.85, 4.75, 5.95]
for mdl in range(5):
lax.plot([xpos[mdl], xpos[mdl] + .15], [0, 0], '-', color=colors_[mdl], **fmt_m)
lax.text(xpos[mdl] + .2 , 0, mdl_lbl[mdl], va='center', fontproperties=fp) # va: equivalent to verticalalignment
lax.set_ylim([-1,1])
lax.set_xlim(.5,7)
lax.set_xticks([])
lax.set_yticks([])
subplots_adjust(left=.075, right=.95, wspace=.1,bottom=.2, top=.8)
if ws5.store_figs:
savefig('img/fig_fmdl_pred_subj{}.pdf'.format(conf.subject))
savefig('img/fig_fmdl_pred_subj{}.svg'.format(conf.subject))
pass # suppress output of last function call
Step 11.1: ankle-SLIP
Step 11.2: augmented-SLIP
Step 11.3: factor-SLIP
Step 11.4: full-state-SLIP
Step 11.5: visualize
Step 11.6: run for all subjects (self-contained cell!)
Approach:
In [1]:
# run required cells (notebook should be saved for this, at least the referred cells shoult be up to data)
import mutils.io as mio
import mutils.misc as mi
# load basic config
nb_name = 'AnkleSLIP - find minimal model- Version 3.ipynb'
mi.run_nbcells(nb_name, ['0'])
conf.subject = 1
conf.quiet = True
mi.run_nbcells(nb_name,
['0.1','1', '3', '3.1', '3.2',]) #'4', '4.2a', '4.2b', '4.2c'])
print "\nnotebook initialized"
In [24]:
array([1,2,3],dtype=int)
Out[24]:
In [ ]:
In [58]:
# --- utility functions etc.
# NOTE: these functions (getRed) rely on external variables! conf, d (ws1.dataset_full), ...
#d = ws1.dataset_full # shortcut
# cell_ID 11.0a
def getControlMaps2(state_r, state_l, pr, pl, Tr, yminr, Tl, yminl, m, indices, conf, d=None, rcond=1e-8):
"""
state_r: median -> reference state r
state_l: median -> reference state l
indices: subset of indices to use for regression
m1: state_r -> (pr - pr_ref)
m2: state_l -> (pl - pl_ref)
m3: state_r -> state_l[3:]
m4: state_l -> state_r[3:]
"""
# remove potential too high indices
idx = array([x for x in indices if x < state_r.shape[0] - 1], dtype=int)
#def getPeriodicOrbit2(ICr, Tr, yminr, ICl, Tl, yminl, m, startParams=[14000.,
#[ICp_r, Pp_r, dE_r], [ICp_l, Pp_l, dE_l] =
ICp_r = mean(state_r[idx, :3], axis=0)
ICp_l = mean(state_l[idx, :3], axis=0)
TR = mean(Tr[idx, :])
TL = mean(Tl[idx, :])
ymin_r = mean(yminr[idx, :])
ymin_l = mean(yminl[idx, :])
Pp_r, Pp_l = su.getPeriodicOrbit2(ICp_r, TR, ymin_r, ICp_l, TL, ymin_l,
m, startParams=mean(pr[idx, :], axis=0)[:5])
# now: detrend data
# non-SLIP state
dt_nss_r = fda.dt_movingavg(state_r[:, 3:], conf.dt_window, conf.dt_medfilter)
dt_nss_l = fda.dt_movingavg(state_l[:, 3:], conf.dt_window, conf.dt_medfilter)
# full state
dt_fs_r = fda.dt_movingavg(state_r, conf.dt_window, conf.dt_medfilter)
dt_fs_l = fda.dt_movingavg(state_l, conf.dt_window, conf.dt_medfilter)
# SLIP parameters
dt_pl = fda.dt_movingavg(d.all_param_l, conf.dt_window, conf.dt_medfilter)
dt_pr = fda.dt_movingavg(d.all_param_r, conf.dt_window, conf.dt_medfilter)
# compute non-SLIP state prediction maps (propagators)
_, all_prop_r, _ = fda.fitData(dt_fs_r[idx, :], dt_nss_l[idx, :], nps=1, nrep=25,
sections=[0,], rcond=rcond)
_, all_prop_l, _ = fda.fitData(dt_fs_l[idx, :], dt_nss_r[idx+1, :], nps=1, nrep=25,
sections=[0,], rcond=rcond)
prop_r = fda.meanMat(all_prop_r)
prop_l = fda.meanMat(all_prop_l)
# compute parameter prediction maps
_, all_ctrl_r, _ = fda.fitData(dt_fs_r[idx, :], dt_pr[idx, :], nps=1, nrep=25,
sections=[0,], rcond=rcond)
_, all_ctrl_l, _ = fda.fitData(dt_fs_l[idx, :], dt_pl[idx, :], nps=1, nrep=25,
sections=[0,], rcond=rcond)
ctrl_r = fda.meanMat(all_ctrl_r)
ctrl_l = fda.meanMat(all_ctrl_l)
return (ctrl_r, ctrl_l), (prop_r, prop_l), (ICp_r, ICp_l, Pp_r, Pp_l)
def getPOs(d):
""" returns a list parameters yielding periodic orbits """
n = d.all_IC_r.shape[0]
print "v" * ((n//25) + 1)
TR = vstack(d.TR).squeeze()
TL = vstack(d.TL).squeeze()
yR = vstack(d.yminR).squeeze()
yL = vstack(d.yminL).squeeze()
m = mean(d.masses)
ICr = d.all_IC_r
ICl = d.all_IC_l
pr = []
pl = []
for nr in range(n):
Pp_r, Pp_l = su.getPeriodicOrbit2(ICr[nr, :], TR[nr], yR[nr], ICl[nr,:], TL[nr], yL[nr],
m, startParams=d.all_param_r[nr, :])
pr.append(Pp_r)
pl.append(Pp_l)
if not mod(nr, 25):
sys.stdout.write('.')
print "done!"
return vstack(pr), vstack(pl)
# --- compute controller and perform predictions
def getRed(states_r_o, states_l_o, is_factor=False, nboot=100):
"""
computes the variance reduction by simulation
:args:
states_r_o, states_l_o: full states of the models. [y,vx,vz, <additional states>]
is_factor (bool): if True, factors are recomputed per bootstrap iteration
nboot (int): number of bootstrap iterations
:returns:
vred_step, vred_stride: bootstrapped lists of achieved variance reduction for CoM
"""
vred_step = []
vred_stride = []
dataset = d
# nboot = 100 # quick for now
print "v" * nboot
for bsrep in range(nboot):
#print "rep:", bsrep + 1
sys.stdout.write('.')
idx = randint(0, states_r_o.shape[0], states_r_o.shape[0])
idx.sort()
# --- select regression and prediction subset
oidx = fda.otheridx(idx, states_r_o.shape[0])[:-1] # skip last index in any case
if is_factor:
# compute factor direction out-of-sample
sr_ = fda.dt_movingavg(states_r_o, conf.dt_window, conf.dt_medfilter)
sl_ = fda.dt_movingavg(states_l_o, conf.dt_window, conf.dt_medfilter)
sr_ = sr_ / std(sr_, axis=0)
sl_ = sl_ / std(sl_, axis=0)
facs_r = st.find_factors(sr_[idx, :].T, d.s_param_r[idx, :].T, k=5)
facs_l = st.find_factors(sl_[idx, :].T, d.s_param_l[idx, :].T, k=5)
# old
#fscore_r = dot(facs_r.T, sr_[idx, :].T).T
#fscore_l = dot(facs_l.T, sl_[idx, :].T).T
#states_r = hstack([states_r_o[idx, :3], fscore_r])
#states_l = hstack([states_l_o[idx, :3], fscore_l])
fscore_r = dot(facs_r.T, sr_.T).T
fscore_l = dot(facs_l.T, sl_.T).T
states_r = hstack([states_r_o[:, :3], fscore_r])
states_l = hstack([states_l_o[:, :3], fscore_l])
else:
states_r = states_r_o
states_l = states_l_o
mass = mean(array(d.masses))
# compute control maps
(ctrl_r, ctrl_l), (prop_r, prop_l), (ICp_r, ICp_l, Pp_r, Pp_l) = getControlMaps2(states_r,
states_l, d.all_param_r, d.all_param_l, vstack(d.TR), vstack(d.yminR),
vstack(d.TL), vstack(d.yminL), mass, idx, conf, dataset)
# --- compute baseline params for "other indices":
ICp_r = mean(states_r[oidx, :3], axis=0)
ICp_l = mean(states_l[oidx, :3], axis=0)
TR = mean(vstack(d.TR)[oidx, :])
TL = mean(vstack(d.TL)[oidx, :])
ymin_r = mean(vstack(d.yminR)[oidx, :])
ymin_l = mean(vstack(d.yminL)[oidx, :])
po_r, po_l = su.getPeriodicOrbit2(ICp_r, TR, ymin_r, ICp_l, TL, ymin_l,
mean(array(d.masses)), startParams=mean(d.all_param_r[idx, :], axis=0)[:5])
refstate_r = mean(states_r[oidx, :], axis=0)
refstate_l = mean(states_l[oidx, :], axis=0)
# ---- simulate steps
idat_r = states_r[oidx, :]
idat_l = states_l[oidx, :]
odat_r = idat_l
odat_l = states_r[oidx+1, :]
# detrended
idat_r_dt = fda.dt_movingavg(states_r, conf.dt_window, conf.dt_medfilter)[oidx, :]
idat_l_dt = fda.dt_movingavg(states_l, conf.dt_window, conf.dt_medfilter)[oidx, :]
odat_r_dt = fda.dt_movingavg(states_l, conf.dt_window, conf.dt_medfilter)[oidx, :]
odat_l_dt = fda.dt_movingavg(states_r, conf.dt_window, conf.dt_medfilter)[oidx+1, :]
# detrended "periodic orbit" parameters (they do not correspond to periodic orbit wrt avg. IC)
# compute periodic orbits for every
#ICp_r = mean(states_r[oidx, :3], axis=0)
#ICp_l = mean(states_l[oidx, :3], axis=0)
TRc = vstack(d.TR).squeeze()
TLc = vstack(d.TL).squeeze()
ymin_rc = vstack(d.yminR).squeeze()
ymin_lc = vstack(d.yminL).squeeze()
#po_r = []
#po_l = []
#for ox in oidx:
# po_rc, po_lc = su.getPeriodicOrbit2(states_r[ox,:], TRc[ox], ymin_rc[ox], states_l[ox,:],
# TLc[ox], ymin_lc[ox],
# mean(array(d.masses)), startParams=mean(d.all_param_r[idx, :], axis=0)[:5])
# po_r.append(po_rc)
# po_l.append(po_lc)
po_r = (d.ppr - fda.dt_movingavg(d.ppr, conf.dt_window, conf.dt_medfilter))[oidx, :]
po_l = (d.ppl - fda.dt_movingavg(d.ppl, conf.dt_window, conf.dt_medfilter))[oidx, :]
#po_r = vstack(po_r)
#po_l = vstack(po_l)
#sys.stdout.write('p')
#po_r = (d.all_param_r - fda.dt_movingavg(d.all_param_r, conf.dt_window, conf.dt_medfilter))[oidx, :]
#po_l = (d.all_param_l - fda.dt_movingavg(d.all_param_l, conf.dt_window, conf.dt_medfilter))[oidx, :]
pred_r_1 = [] # a step
pred_r_2 = [] # a stride
for step in range(len(oidx)):
# --- first step
# --- --- update SLIP params
##p_up_r = dot(ctrl_r, idat_r[step,:] - refstate_r)
p_up_r = dot(ctrl_r, idat_r_dt[step,:])
##k, alpha, L0, beta, dE = po_r + p_up_r
k, alpha, L0, beta, dE = po_r[step, :] + p_up_r
pars = [k, L0, mass, alpha, beta, dE]
try:
IC = idat_r[step, :3]
fs1 = sl.qSLIP_step3D(IC, pars)[1][-1, [1,3,5]]
except (ValueError, TypeError):
# no SLIP computation possible -> use mean value as prediction
sys.stdout.write('x')
fs1 = ICp_l
# --- second step
# --- --- propagate additional states
##augstate = dot(prop_r, idat_r[step,:] - refstate_r)
augstate = dot(prop_r, idat_r_dt[step, :])
##refstate_c = refstate_l[:3]
refstate_c = idat_l[step, :3] - idat_l_dt[step, :3]
p_up_l = dot(ctrl_l, hstack([fs1 - refstate_c, augstate]))
##k, alpha, L0, beta, dE = po_l + p_up_l
k, alpha, L0, beta, dE = po_l[step, :] + p_up_l
pars = [k, L0, mass, alpha, beta, dE]
try:
fs2 = sl.qSLIP_step3D(fs1, pars)[1][-1, [1,3,5]]
except (ValueError, TypeError):
sys.stdout.write('X')
fs2 = ICp_r
pred_r_1.append(fs1)
pred_r_2.append(fs2)
#if not mod(step, 50):
# sys.stdout.write('.')
pred_r_1 = vstack(pred_r_1)
pred_r_2 = vstack(pred_r_2)
#print "done"
#sys.stdout.write('.')
##vred_step.append(var(odat_r[:,:3] - pred_r_1, axis=0) / var(odat_r[:,:3], axis=0))
##vred_stride.append(var(odat_l[:,:3] - pred_r_2, axis=0) / var(odat_l[:,:3], axis=0))
vred_step.append(var(odat_r[:,:3] - pred_r_1, axis=0) / var(odat_r_dt[:,:3], axis=0))
vred_stride.append(var(odat_l[:,:3] - pred_r_2, axis=0) / var(odat_l_dt[:,:3], axis=0))
print " done"
return vred_step, vred_stride
# --- compute prediction directly (linear fit)
def getRed_d(states_r_o, states_l_o, is_factor=False, nboot=100):
"""
computes the variance reduction by direct prediction
:args:
states_r, states_l: full states of the models. [y,vx,vz, <additional states>]
is_factor (bool): If true, comptue the factors from the full state (separately per bootstrap iteration)
nboot (int): number of bootstrap iterations
:returns:
vred_step, vred_stride: bootstrapped lists of achieved variance reduction for CoM
"""
vred_step = []
vred_stride = []
#nboot = 100 # quicker for now, usually: 50
print "V" * nboot
for bsrep in range(nboot):
#print "rep:", bsrep + 1
sys.stdout.write('.')
idx = randint(0, states_r_o.shape[0] - 1, states_r_o.shape[0])
idx.sort()
# --- select regression and prediction subset
oidx = fda.otheridx(idx, states_r_o.shape[0])[:-1] # skip last index in any case
if is_factor:
# compute factor direction out-of-sample
sr_ = fda.dt_movingavg(states_r_o, conf.dt_window, conf.dt_medfilter)
sl_ = fda.dt_movingavg(states_l_o, conf.dt_window, conf.dt_medfilter)
sr_ = sr_ / std(sr_, axis=0)
sl_ = sl_ / std(sl_, axis=0)
facs_r = st.find_factors(sr_[idx, :].T, d.s_param_r[idx, :].T, k=5)
facs_l = st.find_factors(sl_[idx, :].T, d.s_param_l[idx, :].T, k=5)
# old:
#fscore_r = dot(facs_r.T, sr_[idx, :].T).T
#fscore_l = dot(facs_l.T, sl_[idx, :].T).T
#states_r = hstack([states_r_o[idx, :3], fscore_r])
#states_l = hstack([states_l_o[idx, :3], fscore_l])
fscore_r = dot(facs_r.T, sr_.T).T
fscore_l = dot(facs_l.T, sl_.T).T
states_r = hstack([states_r_o[:, :3], fscore_r])
states_l = hstack([states_l_o[:, :3], fscore_l])
else:
states_r = states_r_o
states_l = states_l_o
# --- compute baseline params for "other indices":
ICp_r = mean(states_r[oidx, :3], axis=0)
ICp_l = mean(states_l[oidx, :3], axis=0)
refstate_r = mean(states_r[oidx, :], axis=0)
refstate_l = mean(states_l[oidx, :], axis=0)
# --- compute regression
As1 = []
As2 = []
# for fair comparison: also bootstrap and average regression maps
for irep in range(nboot):
subidx = randint(0, len(idx), len(idx))
sidx = idx[subidx]
idat = fda.dt_movingavg(states_r[sidx, :], conf.dt_window, conf.dt_medfilter)
odat1 = fda.dt_movingavg(states_l[sidx, :3], conf.dt_window, conf.dt_medfilter)
odat2 = fda.dt_movingavg(states_r[sidx + 1, :3], conf.dt_window, conf.dt_medfilter)
As1.append(dot( odat1.T, pinv(idat.T, rcond=1e-7) ))
As2.append(dot( odat2.T, pinv(idat.T, rcond=1e-7) ))
A1 = fda.meanMat(As1)
A2 = fda.meanMat(As2)
idat_r = states_r[oidx, :] - refstate_r
odat_r = states_l[oidx, :3] - refstate_l[:3]
pred_r_1 = dot(A1, idat_r.T).T
odat_l = states_r[oidx+1, :3] - refstate_r[:3]
pred_r_2 = dot(A2, idat_r.T).T
#old
vred_step.append(var(odat_r[:,:3] - pred_r_1, axis=0) / var(odat_r[:,:3], axis=0))
vred_stride.append(var(odat_l[:,:3] - pred_r_2, axis=0) / var(odat_l[:,:3], axis=0))
#
print " done"
return vred_step, vred_stride
In [3]:
#bs1 = ppr - fda.dt_movingavg(ppr, conf.dt_window, conf.dt_medfilter)
#bs2 = d.all_param_r - fda.dt_movingavg(d.all_param_r, conf.dt_window, conf.dt_medfilter)##
#figure(figsize=(11,4))
#for dim in range(5):
# subplot(2,3,dim+1)
# plot(bs1[:,dim], 'r-')
# plot(bs2[:,dim], 'g+')
In [4]:
#std(ppr - d.all_param_r, axis=0)
In [55]:
# cell_ID 11.1a
# --- build dataset
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
if False: # method changed
try:
fname = 'tmp/ppr_ppl_s{}.dict'.format(conf.subject)
p = mio.mload(fname)
ppr = p['ppr']
ppl = p['ppl']
print "data loaded"
except IOError:
#p = {'ppr' : ppr, 'ppl' : ppl}
ppr, ppl = getPOs(d)
mio.msave(fname, p)
d.ppr = d.all_param_r
d.ppl = d.all_param_l
In [56]:
# cell_ID 11.1b
# compute using simulation
vr_step_ank, vr_stride_ank = getRed(states_r, states_l)
# skip:
# compute using affine model with (bootstrap-)averaged prediction matrices
# note: this is not technically sound
vr_step_d_ank, vr_stride_d_ank = getRed_d(states_r, states_l, nboot=10)
# compute using affine model with single prediction matrix per bootstrap iteration
sr_dt = fda.dt_movingavg(states_r, conf.dt_window, conf.dt_medfilter)
sl_dt = fda.dt_movingavg(states_l, conf.dt_window, conf.dt_medfilter)
vr_step_dx_ank = st.predTest(sr_dt, sl_dt[:,:3], out_of_sample=True)
vr_stride_dx_ank = st.predTest(sr_dt[:-1,:], sr_dt[1:,:3], out_of_sample=True)
In [57]:
#vr_step_dx_ank = st.predTest(sr_dt, sl_dt[:,:3], out_of_sample=True)
#vr_stride_dx_ank = st.predTest(sr_dt[:-1,:], sr_dt[1:,:3], out_of_sample=True)
# cell_ID 11.1c
# --- quick visualization
figure()
subplot(1,2,1)
b1 = boxplot(vstack(vr_step_ank), positions = arange(3) + 1, widths=.15)
b2 = boxplot(vstack(vr_step_d_ank), positions = arange(3) + 1.2, widths=.15)
b3 = boxplot(vstack(vr_step_dx_ank), positions = arange(3) + 1.4, widths=.15)
mi.recolor(b1, 'b')
mi.recolor(b2, 'r')
mi.recolor(b3, 'k')
title('subj {} ankle SLIP: 1 step'.format(conf.subject))
ylim(0,1.9)
xlim(0.5,3.8)
subplot(1,2,2)
b1 = boxplot(vstack(vr_stride_ank), positions = arange(3) + 1, widths=.15)
b2 = boxplot(vstack(vr_stride_d_ank), positions = arange(3) + 1.2, widths=.15)
b3 = boxplot(vstack(vr_stride_dx_ank), positions = arange(3) + 1.4, widths=.15)
mi.recolor(b1, 'b')
mi.recolor(b2, 'r')
mi.recolor(b3, 'k')
title('subj {} ankle SLIP: 1 stride'.format(conf.subject))
ylim(0,1.9)
xlim(0.5,3.8)
pass
#vred = var(odat_r[:,:3] - pred_r_1, axis=0) / var(odat_r[:,:3], axis=0)
#print vred
In [65]:
# cell_ID 11.2a
states_r = hstack([d.all_IC_r, roll(d.all_param_l, 1, axis=0) ]) # actually, the first state is invalid
states_l = hstack([d.all_IC_l, d.all_param_r])
vr_step_aug, vr_stride_aug = getRed(states_r, states_l)
# compute using affine model with (bootstrap-)averaged prediction matrices
# NOTE: this is not technically sound!
vr_step_d_aug, vr_stride_d_aug = getRed_d(states_r, states_l, nboot=10)
# compute using affine model with single prediction matrix per bootstrap iteration
sr_dt = fda.dt_movingavg(states_r, conf.dt_window, conf.dt_medfilter)
sl_dt = fda.dt_movingavg(states_l, conf.dt_window, conf.dt_medfilter)
vr_step_dx_aug = st.predTest(sr_dt, sl_dt[:,:3], out_of_sample=True)
vr_stride_dx_aug = st.predTest(sr_dt[:-1,:], sr_dt[1:,:3], out_of_sample=True)
In [66]:
# --- quick visualization
# cell_ID 11.2b
figure()
subplot(1,2,1)
b1 = boxplot(vstack(vr_step_aug), positions = arange(3) + 1, widths=.15)
b2 = boxplot(vstack(vr_step_d_aug), positions = arange(3) + 1.2, widths=.15)
b3 = boxplot(vstack(vr_step_dx_aug), positions = arange(3) + 1.4, widths=.15)
mi.recolor(b1, 'b')
mi.recolor(b2, 'r')
mi.recolor(b3, 'k')
title('aug(r) SLIP:\n 1 step')
ylim(0,1.1)
xlim(0.5,3.8)
subplot(1,2,2)
b1 = boxplot(vstack(vr_stride_aug), positions = arange(3) + 1, widths=.15)
b2 = boxplot(vstack(vr_stride_d_aug), positions = arange(3) + 1.2, widths=.15)
b3 = boxplot(vstack(vr_stride_dx_aug), positions = arange(3) + 1.4, widths=.15)
mi.recolor(b1, 'b')
mi.recolor(b2, 'r')
mi.recolor(b3, 'k')
title('aug(r) SLIP:\n 1 stride')
ylim(0,1.1)
xlim(0.5,3.8)
pass
# --- quick visualization
figure()
subplot(1,2,1)
b1 = boxplot(vstack(vr_step_aug))
b2 = boxplot(vstack(vr_step_ank))
mi.recolor(b1,'r')
mi.recolor(b2,'b')
subplot(1,2,2)
b1 = boxplot(vstack(vr_stride_aug))
b2 = boxplot(vstack(vr_stride_ank))
mi.recolor(b1,'r')
mi.recolor(b2,'b')
pass
In [60]:
# this must change to out-of-sample factor computation!
#ws1.facs_r = st.find_factors(ws1.dataset.s_kin_r[:, sel_idx].T, ws1.dataset.s_param_r.T, k=conf.n_factors)
#ws1.facs_l = st.find_factors(ws1.dataset.s_kin_l[:, sel_idx].T, ws1.dataset.s_param_l.T, k=conf.n_factors)
#ws1.fscore_r = dot(ws1.facs_r.T, ws1.dataset.s_kin_r[:, sel_idx].T).T
#ws1.fscore_r_doc = "scores of data when projected on the factors (right)"
#ws1.fscore_l = dot(ws1.facs_l.T, ws1.dataset.s_kin_l[:, sel_idx].T).T
#ws1.fscore_l_doc = "scores of data when projected on the factors (left)"
# cell_ID 11.3a
#states_r = hstack([d.all_IC_r, ws1.fscore_r ])
#states_l = hstack([d.all_IC_l, ws1.fscore_l])
# state is the same as "full state SLIP", but factors are computed in the getRed functions!
use_factors = True
if use_factors: #out of sample predictions for factors
lbls = [label for label in d.kin_labels if '-' in label and 'v_' not in label]
ws1.k.selection = lbls
k_ap_r = ws1.k.get_kin_apex(d.all_phases_r)
k_ap_r = hstack(k_ap_r).T
k_ap_l = ws1.k.get_kin_apex(d.all_phases_l)
k_ap_l = hstack(k_ap_l).T
states_r = hstack([d.all_IC_r, k_ap_r])
states_l = hstack([d.all_IC_l, k_ap_l])
else:
states_r = hstack([d.all_IC_r, ws1.fscore_r ])
states_l = hstack([d.all_IC_l, ws1.fscore_l])
vr_step_fac, vr_stride_fac = getRed(states_r, states_l, is_factor=use_factors)
# compute using affine model with (bootstrap-)averaged prediction matrices
# note: this is not technically sound (averaging over a sub-sample...)
vr_step_d_fac, vr_stride_d_fac = getRed_d(states_r, states_l, is_factor=use_factors, nboot=10)
# compute using affine model with single prediction matrix per bootstrap iteration
sr_dt = fda.dt_movingavg(states_r, conf.dt_window, conf.dt_medfilter)
sl_dt = fda.dt_movingavg(states_l, conf.dt_window, conf.dt_medfilter)
vr_step_dx_fac = st.predTest(sr_dt, sl_dt[:,:3], out_of_sample=True)
vr_stride_dx_fac = st.predTest(sr_dt[:-1,:], sr_dt[1:,:3], out_of_sample=True)
In [61]:
# cell_ID 11.3b
# --- quick visualization
figure(figsize=(10,5))
subplot(1,2,1)
b1 = boxplot(vstack(vr_step_fac), positions = arange(3) + 1, widths=.15)
b2 = boxplot(vstack(vr_step_d_fac), positions = arange(3) + 1.2, widths=.15)
b3 = boxplot(vstack(vr_step_dx_fac), positions = arange(3) + 1.4, widths=.15)
mi.recolor(b1, 'b')
mi.recolor(b2, 'r')
mi.recolor(b3, 'k')
title('subject {}, factor SLIP:\n 1 step (r)'.format(conf.subject))
ylim(0,1.1)
xlim(0.5,3.8)
subplot(1,2,2)
b1 = boxplot(vstack(vr_stride_fac), positions = arange(3) + 1, widths=.15)
b2 = boxplot(vstack(vr_stride_d_fac), positions = arange(3) + 1.2, widths=.15)
b3 = boxplot(vstack(vr_stride_dx_fac), positions = arange(3) + 1.4, widths=.15)
mi.recolor(b1, 'b')
mi.recolor(b2, 'r')
mi.recolor(b3, 'k')
title('subject {}, factor SLIP:\n 1 stride (r)'.format(conf.subject))
ylim(0,1.1)
xlim(0.5,3.8)
lax = axes([.6, .2, .3, .25], frameon=True)
textprops={'fontsize' : 9}
def plotentry(h, color, ltext, x=0, lw=1):
plot([x, x+1],[h, h], color=color, linewidth=lw)
text(x + 1.2, h, ltext, **textprops)
plotentry(0, 'b','simulation')
plotentry(1,'r', 'affine model [averaged matrix]') # fs affine
plotentry(2, 'k', 'affine model [single matrix]') # factor
lax.set_yticks([])
lax.set_xticks([])
lax.set_xlim(-.5,10)
lax.set_ylim([3, -1])
pass
In [62]:
# cell_ID 11.4a
lbls = [label for label in d.kin_labels if '-' in label and 'v_' not in label]
ws1.k.selection = lbls
k_ap_r = ws1.k.get_kin_apex(d.all_phases_r)
k_ap_r = hstack(k_ap_r).T
k_ap_l = ws1.k.get_kin_apex(d.all_phases_l)
k_ap_l = hstack(k_ap_l).T
states_r = hstack([d.all_IC_r, k_ap_r ])
states_l = hstack([d.all_IC_l, k_ap_l])
vr_step_full, vr_stride_full = getRed(states_r, states_l)
# compute using affine model with (bootstrap-)averaged prediction matrices
# NOTE: this is actually not technically sound!
vr_step_d_full, vr_stride_d_full = getRed_d(states_r, states_l, nboot=10)
# compute using affine model with single prediction matrix per bootstrap iteration
sr_dt = fda.dt_movingavg(states_r, conf.dt_window, conf.dt_medfilter)
sl_dt = fda.dt_movingavg(states_l, conf.dt_window, conf.dt_medfilter)
vr_step_dx_full = st.predTest(sr_dt, sl_dt[:,:3], out_of_sample=True)
vr_stride_dx_full = st.predTest(sr_dt[:-1,:], sr_dt[1:,:3], out_of_sample=True)
In [63]:
# cell_ID 11.4b
# --- quick visualization
#vr_step_dx_full = st.predTest(sr_dt, sl_dt[:,:3], out_of_sample=False)
#vr_stride_dx_full = st.predTest(sr_dt[:-1,:], sr_dt[1:,:3], out_of_sample=False)
figure()
subplot(1,2,1)
b1 = boxplot(vstack(vr_step_full), positions = arange(3) + 1, widths=.15)
b2 = boxplot(vstack(vr_step_d_full), positions = arange(3) + 1.2, widths=.15)
b3 = boxplot(vstack(vr_step_dx_full), positions = arange(3) + 1.4, widths=.15)
mi.recolor(b1, 'b')
mi.recolor(b2, 'r')
mi.recolor(b3, 'k')
title('full state(r) SLIP:\n 1 step')
ylim(0,1.1)
xlim(0.5,3.8)
subplot(1,2,2)
b1 = boxplot(vstack(vr_stride_full), positions = arange(3) + 1, widths=.15)
b2 = boxplot(vstack(vr_stride_d_full), positions = arange(3) + 1.2, widths=.15)
b3 = boxplot(vstack(vr_stride_dx_full), positions = arange(3) + 1.4, widths=.15)
mi.recolor(b1, 'b')
mi.recolor(b2, 'r')
mi.recolor(b3, 'k')
title('full state(r) SLIP:\n 1 stride')
ylim(0,1.1)
xlim(0.5,3.8)
pass
In [ ]:
In [14]:
# --- "full" figure
figure(figsize=(12,5))
# --- one step
subplot(1,2,1)
b0a = boxplot(vstack(vr_step_full), positions = arange(3)*2 + .6, widths=.15)
b0b = boxplot(vstack(vr_step_dx_full), positions = arange(3)*2 + .8, widths=.15)
b1 = boxplot(vstack(vr_step_fac), positions = arange(3)*2 + 1, widths=.15)
b2 = boxplot(vstack(vr_step_dx_fac), positions = arange(3)*2 + 1.2, widths=.15)
b4 = boxplot(vstack(vr_step_ank), positions = arange(3)*2 + 1.4, widths=.15)
b5 = boxplot(vstack(vr_step_dx_ank), positions = arange(3)*2 + 1.6, widths=.15)
b3 = boxplot(vstack(vr_step_aug), positions = arange(3)*2 + 1.8, widths=.15)
b3b = boxplot(vstack(vr_step_dx_aug), positions = arange(3)*2 + 2.0, widths=.15)
mi.recolor(b0a,'#006767', lw=2)
mi.recolor(b0b,'#339a9a')
mi.recolor(b1,'k', lw=2)
mi.recolor(b2,'#797979')
mi.recolor(b3,'r', lw=2)
mi.recolor(b3b,'#ff6767')
mi.recolor(b4,'#009900', lw=2)
mi.recolor(b5,'#55ff55')
plot([.2, 7.],[1., 1.], 'k--')
ylim(0,1.45)
xlim(0.2,7)
gca().set_xticks(arange(3)*2 + 1.5)
gca().set_xticklabels(['z', 'vy', 'vx'])
ylabel('relative remaining variance')
title('subject {}: 1 step (right)'.format(conf.subject))
xlabel('CoM state predicted')
# --- one stride
subplot(1,2,2)
b0a = boxplot(vstack(vr_stride_full), positions = arange(3)*2 + .6, widths=.15)
b0b = boxplot(vstack(vr_stride_dx_full), positions = arange(3)*2 + .8, widths=.15)
b1 = boxplot(vstack(vr_stride_fac), positions = arange(3)*2 + 1, widths=.15)
b2 = boxplot(vstack(vr_stride_dx_fac), positions = arange(3)*2 + 1.2, widths=.15)
b4 = boxplot(vstack(vr_stride_ank), positions = arange(3)*2 + 1.4, widths=.15)
b5 = boxplot(vstack(vr_stride_dx_ank), positions = arange(3)*2 + 1.6, widths=.15)
b3 = boxplot(vstack(vr_stride_aug), positions = arange(3)*2 + 1.8, widths=.15)
b3b = boxplot(vstack(vr_stride_dx_aug), positions = arange(3)*2 + 2.0, widths=.15)
mi.recolor(b0a,'#006767', lw=2)
mi.recolor(b0b,'#339a9a')
mi.recolor(b1,'k', lw=2)
mi.recolor(b2,'#797979')
mi.recolor(b3,'r', lw=2)
mi.recolor(b3b,'#ff6767')
mi.recolor(b4,'#009900', lw=2)
mi.recolor(b5,'#55ff55')
plot([.2, 7.],[1., 1.], 'k--')
ylim(0,1.45)
xlim(0.2,7)
gca().set_xticks(arange(3)*2 + 1.5)
gca().set_xticklabels(['z', 'vy', 'vx'])
ylabel('relative remaining variance')
title('subject {}: 1 stride (right)'.format(conf.subject))
xlabel('CoM state predicted')
# legend axis
lax = axes([.6, .2, .3, .25], frameon=True)
textprops={'fontsize' : 9}
def plotentry(h, color, ltext, x=0, lw=1):
plot([x, x+1],[h, h], color=color, linewidth=lw)
text(x + 1.2, h, ltext, **textprops)
plotentry(0, '#006767','full state [sim]', lw=2)
plotentry(0,'#339a9a', 'full state [affine]', x=5) # fs affine
plotentry(1, 'k', 'factors [sim]', lw=2) # factor
plotentry(1, '#797979', 'factors [affine]', x=5) # fact. affine
plotentry(2,'#009900', 'ankle state [sim]', lw=2) # ank
plotentry(2,'#55ff55', 'ankle state [affine]', x=5,) # ank affine
plotentry(3,'r', 'SLIP param. [sim]', lw=2) # aug
plotentry(3,'#ff6767', 'SLIP param. [affine]', x=5) # aug affine
lax.set_yticks([])
lax.set_xticks([])
lax.set_xlim(-.5,10)
lax.set_ylim([4, -1])
Out[14]:
In [2]:
import mutils.io as mio
import mutils.misc as mi
# load basic config
nb_name = 'AnkleSLIP - find minimal model- Version 3.ipynb'
ws11 = mio.saveable()
mi.run_nbcells(nb_name, ['0'])
conf.subject = 1
conf.quiet = True
ws11.vr_step_full = []
ws11.vr_step_dx_full = []
ws11.vr_step_fac = []
ws11.vr_step_dx_fac = []
ws11.vr_step_ank = []
ws11.vr_step_dx_ank = []
ws11.vr_step_aug = []
ws11.vr_step_dx_aug = []
ws11.vr_stride_full = []
ws11.vr_stride_dx_full = []
ws11.vr_stride_fac = []
ws11.vr_stride_dx_fac = []
ws11.vr_stride_ank = []
ws11.vr_stride_dx_ank = []
ws11.vr_stride_aug = []
ws11.vr_stride_dx_aug = []
print '\ndone'
In [ ]:
# compute data for all subjects
for subject_ in [1, 2, 3, 7]:
conf.subject = subject_
conf.quiet = True
# reload data
mi.run_nbcells(nb_name,
['0.1','1', '3', '3.1', '3.2', '11.0a'])
print "\nnotebook initialized for subject {}".format(subject_)
mi.run_nbcells(nb_name,
['11.1a', '11.1b', '11.2a', '11.3a', '11.4a'])
ws11.vr_step_full.append(vstack(vr_step_full))
ws11.vr_step_dx_full.append(vstack(vr_step_dx_full))
ws11.vr_step_fac.append(vstack(vr_step_fac))
ws11.vr_step_dx_fac.append(vstack(vr_step_dx_fac))
ws11.vr_step_ank.append(vstack(vr_step_ank))
ws11.vr_step_dx_ank.append(vstack(vr_step_dx_ank))
ws11.vr_step_aug.append(vstack(vr_step_aug))
ws11.vr_step_dx_aug.append(vstack(vr_step_dx_aug))
ws11.vr_stride_full.append(vstack(vr_stride_full))
ws11.vr_stride_dx_full.append(vstack(vr_stride_dx_full))
ws11.vr_stride_fac.append(vstack(vr_stride_fac))
ws11.vr_stride_dx_fac.append(vstack(vr_stride_dx_fac))
ws11.vr_stride_ank.append(vstack(vr_stride_ank))
ws11.vr_stride_dx_ank.append(vstack(vr_stride_dx_ank))
ws11.vr_stride_aug.append(vstack(vr_stride_aug))
ws11.vr_stride_dx_aug.append(vstack(vr_stride_dx_aug))
In [5]:
import mutils.plotting as mplt
colors_ = mplt.colorset_distinct
fmt1 = {'marker' : '>',
'markersize' : 3.5,
'markeredgecolor' : 'None'}
fmt2 = {'marker' : '<',
'markersize' : 2.5,
'markeredgecolor' : 'None'}
import matplotlib.font_manager as fmt
FM = fmt.FontManager()
#fnt = FM.findfont('Times New Roman') # for Proc R Soc
fnt = FM.findfont('Arial') # for NComm
fp = fmt.FontProperties(fname=fnt)
fp.set_size(9.0)
# make font in math the same as in text
import matplotlib as mpl
mpl.rcParams['mathtext.default'] = 'regular' # this is the key
In [14]:
# create the figure
vis_subject = 7 # from 1,2,3,7
subjects_ = [1, 2, 3, 7]
plot_subject = find(array(subjects_) == vis_subject)[0]
# horizontal positions of elements
posP = arange(ws11.vr_step_full[0].shape[1]) * 2.5
posPb = posP + .2
fig = figure(figsize=(18./2.56, 3))
# --- left panel
ax0 = fig.add_subplot(1, 2, 1)
# boxplots
b1 = boxplot(ws11.vr_step_full[plot_subject], positions = posP + .1, widths=.14, sym='')
b2 = boxplot(ws11.vr_step_fac[plot_subject], positions = posP + .6, widths=.14, sym='')
b3 = boxplot(ws11.vr_step_ank[plot_subject], positions = posP + 1.1, widths=.14, sym='')
b4 = boxplot(ws11.vr_step_aug[plot_subject], positions = posP + 1.6, widths=.14, sym='')
mi.recolor(b1, colors_[0], lw=1.25)
mi.recolor(b2, colors_[1], lw=1.25)
mi.recolor(b3, colors_[2], lw=1.25)
mi.recolor(b4, colors_[3], lw=1.25)
b1b = boxplot(ws11.vr_step_dx_full[plot_subject], positions = posPb + .1, widths=.14, sym='')
b2b = boxplot(ws11.vr_step_dx_fac[plot_subject], positions = posPb + .6, widths=.14, sym='')
b3b = boxplot(ws11.vr_step_dx_ank[plot_subject], positions = posPb + 1.1, widths=.14, sym='')
b4b = boxplot(ws11.vr_step_dx_aug[plot_subject], positions = posPb + 1.6, widths=.14, sym = '')
mi.recolor(b1b, colors_[0], lw=.5)
mi.recolor(b2b, colors_[1], lw=.5)
mi.recolor(b3b, colors_[2], lw=.5)
mi.recolor(b4b, colors_[3], lw=.5)
for snr, subject_ in enumerate([1, 2, 3, 7]):
if subject_ == vis_subject:
continue # skip triangle for subject for which the boxplots are drawn
for dim in range(ws11.vr_step_full[0].shape[1]):
plot(posP[dim] + .1 - .1, median(ws11.vr_step_full[snr][:, dim]),
color=colors_[0], **fmt1)
plot(posPb[dim] + .1 + .1, median(ws11.vr_step_dx_full[snr][:, dim]),
color=colors_[0], **fmt2)
plot(posP[dim] + .6 - .1, median(ws11.vr_step_fac[snr][:, dim]),
color=colors_[1], **fmt1)
plot(posPb[dim] + .6 + .1, median(ws11.vr_step_dx_fac[snr][:, dim]),
color=colors_[1], **fmt2)
plot(posP[dim] + 1.1 - .1, median(ws11.vr_step_ank[snr][:, dim]),
color=colors_[2], **fmt1)
plot(posPb[dim] + 1.1 + .1, median(ws11.vr_step_dx_ank[snr][:, dim]),
color=colors_[2], **fmt2)
plot(posP[dim] + 1.6 - .1, median(ws11.vr_step_aug[snr][:, dim]),
color=colors_[3], **fmt1)
plot(posPb[dim] + 1.6 + .1, median(ws11.vr_step_dx_aug[snr][:, dim]),
color=colors_[3], **fmt2)
# --- right panel
ax1 = fig.add_subplot(1, 2, 2)
# boxplots
b1 = boxplot(ws11.vr_stride_full[plot_subject], positions = posP + .1, widths=.14, sym='')
b2 = boxplot(ws11.vr_stride_fac[plot_subject], positions = posP + .6, widths=.14, sym='')
b3 = boxplot(ws11.vr_stride_ank[plot_subject], positions = posP + 1.1, widths=.14, sym='')
b4 = boxplot(ws11.vr_stride_aug[plot_subject], positions = posP + 1.6, widths=.14, sym='')
mi.recolor(b1, colors_[0], lw=1.25)
mi.recolor(b2, colors_[1], lw=1.25)
mi.recolor(b3, colors_[2], lw=1.25)
mi.recolor(b4, colors_[3], lw=1.25)
b1b = boxplot(ws11.vr_stride_dx_full[plot_subject], positions = posPb + .1, widths=.14, sym='')
b2b = boxplot(ws11.vr_stride_dx_fac[plot_subject], positions = posPb + .6, widths=.14, sym='')
b3b = boxplot(ws11.vr_stride_dx_ank[plot_subject], positions = posPb + 1.1, widths=.14, sym='')
b4b = boxplot(ws11.vr_stride_dx_aug[plot_subject], positions = posPb + 1.6, widths=.14, sym = '')
mi.recolor(b1b, colors_[0], lw=.5)
mi.recolor(b2b, colors_[1], lw=.5)
mi.recolor(b3b, colors_[2], lw=.5)
mi.recolor(b4b, colors_[3], lw=.5)
for snr, subject_ in enumerate([1, 2, 3, 7]):
if subject_ == vis_subject:
continue # skip triangle for subject for which the boxplots are drawn
for dim in range(ws11.vr_stride_full[0].shape[1]):
plot(posP[dim] + .1 - .1, median(ws11.vr_stride_full[snr][:, dim]),
color=colors_[0], **fmt1)
plot(posPb[dim] + .1 + .1, median(ws11.vr_stride_dx_full[snr][:, dim]),
color=colors_[0], **fmt2)
plot(posP[dim] + .6 - .1, median(ws11.vr_stride_fac[snr][:, dim]),
color=colors_[1], **fmt1)
plot(posPb[dim] + .6 + .1, median(ws11.vr_stride_dx_fac[snr][:, dim]),
color=colors_[1], **fmt2)
plot(posP[dim] + 1.1 - .1, median(ws11.vr_stride_ank[snr][:, dim]),
color=colors_[2], **fmt1)
plot(posPb[dim] + 1.1 + .1, median(ws11.vr_stride_dx_ank[snr][:, dim]),
color=colors_[2], **fmt2)
plot(posP[dim] + 1.6 - .1, median(ws11.vr_stride_aug[snr][:, dim]),
color=colors_[3], **fmt1)
plot(posPb[dim] + 1.6 + .1, median(ws11.vr_stride_dx_aug[snr][:, dim]),
color=colors_[3], **fmt2)
# --- formatting
ax0.set_xticks(posP + .85)
ax0.set_yticks(arange(6) * .2)
ax0.set_yticklabels(arange(6) * .2, fontproperties=fp)
ax0.set_xticklabels(['height', 'horiz. velocity', 'lat.velocity'], fontproperties=fp)
ax0.set_xlabel('predicted CoM state', fontproperties=fp)
ax0.get_xaxis().set_label_coords(1.02, -0.12)
ax0.set_ylabel('relative remaining variance', fontproperties=fp)
ax0.plot([-10, 20], [1, 1], 'k--')
ax0.set_xlim(-1, ws11.vr_step_full[0].shape[1] * 3 - 1)
ax0.set_ylim(0, 1.15)
ax1.set_xticks(posP + .85)
ax1.set_yticks(arange(6) * .2)
ax1.get_yaxis().set_ticks_position('right')
ax1.set_yticklabels(arange(6) * .2, fontproperties=fp)
ax1.get_yaxis().set_ticks_position('both')
ax1.set_xticklabels(['height', 'horiz. velocity', 'lat.velocity'], fontproperties=fp)
ax1.plot([-10, 20], [1, 1], 'k--')
ax1.set_xlim(-1, ws11.vr_step_full[0].shape[1] * 3 - 1)
ax1.set_ylim(0, 1.15)
ax0.set_title('prediction accuracy after one step (right)', fontproperties=fp)
ax1.set_title('prediction accuracy after one stride (right)', fontproperties=fp)
subplots_adjust(right=.95, bottom=.15, top=.92, left=.075, wspace=.025)
# legend axis
axl = axes([.53, .2, .41, .25], frameon=False)
#axl.plot(randn(100))
axl.set_xticks([])
axl.set_yticks([])
axl.plot([0,.2], [0, 0], '-', color=colors_[0], lw=1.25)
axl.text(0.25, 0, 'full state SLIP', fontproperties=fp, va='center')
axl.plot([0,.2], [-1, -1], '-', color=colors_[1], lw=1.25)
axl.text(0.25, -1, 'factor SLIP', fontproperties=fp, va='center')
axl.plot([0,.2], [-2, -2], '-', color=colors_[2], lw=1.25)
axl.text(0.25, -2, 'ankle SLIP', fontproperties=fp, va='center')
axl.plot([0,.2], [-3, -3], '-', color=colors_[3], lw=1.25)
axl.text(0.25, -3, 'augmented SLIP', fontproperties=fp, va='center')
axl.plot([3, 3.05, 3.2], [0, 0, 0], '-', color=colors_[0], lw=.5)
axl.text(3.25, 0, 'full state affine', fontproperties=fp, va='center')
axl.plot([3, 3.05, 3.2], [-1, -1, -1], '-', color=colors_[1], lw=.5)
axl.text(3.25, -1, 'factor state affine', fontproperties=fp, va='center')
axl.plot([3, 3.05, 3.2], [-2, -2, -2], '-', color=colors_[2], lw=.5)
axl.text(3.25, -2, 'ankle SLIP affine', fontproperties=fp, va='center')
axl.plot([3, 3.05, 3.2], [-3, -3, -3], '-', color=colors_[3], lw=.5)
axl.text(3.25, -3, 'augmented SLIP affine', fontproperties=fp, va='center')
axl.set_xlim(-.6,6.5)
axl.set_ylim(-4.5, 0.5)
savefig('img/apex_states_predictions_s{}.pdf'.format(vis_subject))
pass
In [ ]:
In [19]:
# new format
# create the figure
vis_subject = 7 # from 1,2,3,7
subjects_ = [1, 2, 3, 7]
plot_subject = find(array(subjects_) == vis_subject)[0]
# horizontal positions of elements
posP = arange(ws11.vr_step_full[0].shape[1]) * 2.5
posPb = posP + .2
#fig = figure(figsize=(18./2.56, 3.1))
fig = figure(figsize=(18./2.56, 2.7))
# --- left panel
ax0 = fig.add_subplot(2, 1, 1)
# boxplots
b1 = boxplot(ws11.vr_step_full[plot_subject], positions = posP + .1, widths=.14, sym='')
b2 = boxplot(ws11.vr_step_fac[plot_subject], positions = posP + .6, widths=.14, sym='')
b3 = boxplot(ws11.vr_step_ank[plot_subject], positions = posP + 1.1, widths=.14, sym='')
b4 = boxplot(ws11.vr_step_aug[plot_subject], positions = posP + 1.6, widths=.14, sym='')
mi.recolor(b1, colors_[0], lw=1.25)
mi.recolor(b2, colors_[1], lw=1.25)
mi.recolor(b3, colors_[2], lw=1.25)
mi.recolor(b4, colors_[3], lw=1.25)
b1b = boxplot(ws11.vr_step_dx_full[plot_subject], positions = posPb + .1, widths=.14, sym='')
b2b = boxplot(ws11.vr_step_dx_fac[plot_subject], positions = posPb + .6, widths=.14, sym='')
b3b = boxplot(ws11.vr_step_dx_ank[plot_subject], positions = posPb + 1.1, widths=.14, sym='')
b4b = boxplot(ws11.vr_step_dx_aug[plot_subject], positions = posPb + 1.6, widths=.14, sym = '')
mi.recolor(b1b, colors_[0], lw=.5)
mi.recolor(b2b, colors_[1], lw=.5)
mi.recolor(b3b, colors_[2], lw=.5)
mi.recolor(b4b, colors_[3], lw=.5)
for snr, subject_ in enumerate([1, 2, 3, 7]):
if subject_ == vis_subject:
continue # skip triangle for subject for which the boxplots are drawn
for dim in range(ws11.vr_step_full[0].shape[1]):
plot(posP[dim] + .1 - .1, median(ws11.vr_step_full[snr][:, dim]),
color=colors_[0], **fmt1)
plot(posPb[dim] + .1 + .1, median(ws11.vr_step_dx_full[snr][:, dim]),
color=colors_[0], **fmt2)
plot(posP[dim] + .6 - .1, median(ws11.vr_step_fac[snr][:, dim]),
color=colors_[1], **fmt1)
plot(posPb[dim] + .6 + .1, median(ws11.vr_step_dx_fac[snr][:, dim]),
color=colors_[1], **fmt2)
plot(posP[dim] + 1.1 - .1, median(ws11.vr_step_ank[snr][:, dim]),
color=colors_[2], **fmt1)
plot(posPb[dim] + 1.1 + .1, median(ws11.vr_step_dx_ank[snr][:, dim]),
color=colors_[2], **fmt2)
plot(posP[dim] + 1.6 - .1, median(ws11.vr_step_aug[snr][:, dim]),
color=colors_[3], **fmt1)
plot(posPb[dim] + 1.6 + .1, median(ws11.vr_step_dx_aug[snr][:, dim]),
color=colors_[3], **fmt2)
# --- right panel
ax1 = fig.add_subplot(2, 1, 2)
# boxplots
b1 = boxplot(ws11.vr_stride_full[plot_subject], positions = posP + .1, widths=.14, sym='')
b2 = boxplot(ws11.vr_stride_fac[plot_subject], positions = posP + .6, widths=.14, sym='')
b3 = boxplot(ws11.vr_stride_ank[plot_subject], positions = posP + 1.1, widths=.14, sym='')
b4 = boxplot(ws11.vr_stride_aug[plot_subject], positions = posP + 1.6, widths=.14, sym='')
mi.recolor(b1, colors_[0], lw=1.25)
mi.recolor(b2, colors_[1], lw=1.25)
mi.recolor(b3, colors_[2], lw=1.25)
mi.recolor(b4, colors_[3], lw=1.25)
b1b = boxplot(ws11.vr_stride_dx_full[plot_subject], positions = posPb + .1, widths=.14, sym='')
b2b = boxplot(ws11.vr_stride_dx_fac[plot_subject], positions = posPb + .6, widths=.14, sym='')
b3b = boxplot(ws11.vr_stride_dx_ank[plot_subject], positions = posPb + 1.1, widths=.14, sym='')
b4b = boxplot(ws11.vr_stride_dx_aug[plot_subject], positions = posPb + 1.6, widths=.14, sym = '')
mi.recolor(b1b, colors_[0], lw=.5)
mi.recolor(b2b, colors_[1], lw=.5)
mi.recolor(b3b, colors_[2], lw=.5)
mi.recolor(b4b, colors_[3], lw=.5)
for snr, subject_ in enumerate([1, 2, 3, 7]):
if subject_ == vis_subject:
continue # skip triangle for subject for which the boxplots are drawn
for dim in range(ws11.vr_stride_full[0].shape[1]):
plot(posP[dim] + .1 - .1, median(ws11.vr_stride_full[snr][:, dim]),
color=colors_[0], **fmt1)
plot(posPb[dim] + .1 + .1, median(ws11.vr_stride_dx_full[snr][:, dim]),
color=colors_[0], **fmt2)
plot(posP[dim] + .6 - .1, median(ws11.vr_stride_fac[snr][:, dim]),
color=colors_[1], **fmt1)
plot(posPb[dim] + .6 + .1, median(ws11.vr_stride_dx_fac[snr][:, dim]),
color=colors_[1], **fmt2)
plot(posP[dim] + 1.1 - .1, median(ws11.vr_stride_ank[snr][:, dim]),
color=colors_[2], **fmt1)
plot(posPb[dim] + 1.1 + .1, median(ws11.vr_stride_dx_ank[snr][:, dim]),
color=colors_[2], **fmt2)
plot(posP[dim] + 1.6 - .1, median(ws11.vr_stride_aug[snr][:, dim]),
color=colors_[3], **fmt1)
plot(posPb[dim] + 1.6 + .1, median(ws11.vr_stride_dx_aug[snr][:, dim]),
color=colors_[3], **fmt2)
# --- formatting
ax0.set_xticks(posP + .85)
ax0.set_yticks(arange(6) * .2)
ax0.set_yticklabels(arange(6) * .2, fontproperties=fp)
ax0.set_xticklabels([]) #'height', 'horiz. velocity', 'lat.velocity'], fontproperties=fp)
ax0.set_ylabel('rrv', fontproperties=fp)
ax0.plot([-10, 20], [1, 1], 'k--')
ax0.set_xlim(-.7, ws11.vr_step_full[0].shape[1] * 3 - 1.9)
ax0.set_ylim(0, 1.15)
ax0.arrow(-0.5,0.95,0,-0.8, head_width=0.05, head_length=0.1, lw=1, color='#000000')
ax0.text(-0.375,0.5, 'better prediction', ha='left', va='center', rotation=90, fontproperties=fp)
ax1.set_xticks(posP + .85)
ax1.set_yticks(arange(6) * .2)
#ax1.get_yaxis().set_ticks_position('right')
ax1.set_yticklabels(arange(6) * .2, fontproperties=fp)
ax1.get_yaxis().set_ticks_position('both')
ax1.set_xticklabels(['height', 'horizontal velocity', 'lateral velocity'], fontproperties=fp, ha='left')
ax1.text(-0.8, -0.11, 'predicted CoM state:', fontproperties=fp, va='center')
#ax1.set_xlabel('predicted CoM state', fontproperties=fp)
#ax1.get_xaxis().set_label_coords(0.485, -0.20)
ax1.set_ylabel('rrv', fontproperties=fp)
ax1.plot([-10, 20], [1, 1], 'k--')
ax1.set_xlim(-.7, ws11.vr_step_full[0].shape[1] * 3 - 1.9)
ax1.set_ylim(0, 1.15)
ax1.arrow(-0.5,0.95,0,-0.8, head_width=0.05, head_length=0.1, lw=1, color='#000000')
ax1.text(-0.375,0.5, 'better prediction', ha='left', va='center', rotation=90, fontproperties=fp)
#ax0.set_title('prediction accuracy after one step (right)', fontproperties=fp)
#ax1.set_title('prediction accuracy after one stride (right)', fontproperties=fp)
subplots_adjust(right=.985, bottom=.075, top=.98, left=.075, wspace=.025, hspace=0.025)
# legend axis
#axl = axes([.53, .2, .41, .25], frameon=False)
#axl = axes([.075, .91, .91, .08], frameon=False)
axl = axes([.27, .12, .68, .125], frameon=False)
#axl.plot(randn(100))
axl.set_xticks([])
axl.set_yticks([])
axl.plot([0, .2], [0, 0], '-', color=colors_[0], lw=1.25)
axl.text(0.25, 0, 'simulated', fontproperties=fp, va='center')
axl.text(0.1, 1, 'Full state SLIP', fontproperties=fp, va='center')
axl.plot([0, 0.2], [-1, -1], '-', color=colors_[0], lw=.5)
axl.text(0.25, -1, 'linear', fontproperties=fp, va='center')
axl.plot([1, 1.2], [-1, -1], '-', color=colors_[1], lw=.5)
axl.text(1.25, -1, 'linear', fontproperties=fp, va='center')
axl.text(1.1, 1, 'Factor-SLIP', fontproperties=fp, va='center')
axl.plot([1, 1.2], [0, 0], '-', color=colors_[1], lw=1.25)
axl.text(1.25, 0, 'simulated', fontproperties=fp, va='center')
axl.plot([2, 2.2], [-1, -1], '-', color=colors_[2], lw=.5)
axl.text(2.25, -1, 'linear', fontproperties=fp, va='center')
axl.text(2.1, 1, 'Ankle-SLIP', fontproperties=fp, va='center')
axl.plot([2, 2.2], [0, 0], '-', color=colors_[2], lw=1.25)
axl.text(2.25, 0, 'simulated', fontproperties=fp, va='center')
axl.plot([3, 3.2], [-1, -1], '-', color=colors_[3], lw=.5)
axl.text(3.25, -1, 'linear', fontproperties=fp, va='center')
axl.text(3.1, 1, 'Augmented-SLIP', fontproperties=fp, va='center')
axl.plot([3, 3.2], [0, 0], '-', color=colors_[3], lw=1.25)
axl.text(3.25, 0, 'simulated', fontproperties=fp, va='center')
axl.set_xlim(-.1, 4.35)
axl.set_ylim(-1.5, 1.65)
ax0.text(-0., 0.2, 'one step', fontproperties=fp, bbox=dict(facecolor='#ffffff', edgecolor='k', pad=4.0))
ax1.text(-0., 0.2, 'one stride', fontproperties=fp, bbox=dict(facecolor='#ffffff', edgecolor='k', pad=4.0))
savefig('img/apex_states_predictions_s{}_v2.pdf'.format(vis_subject))
print 'figure stored as img/apex_states_predictions_s{}_v2.pdf'.format(vis_subject)
pass
In [7]:
vis_subject
Out[7]:
In [33]:
yax = ax1.get_yaxis().set_ticks_position('left')
In [20]:
# give rrv in numbers
print " rrv of ankle compared to full state models"
print " <apex height> <horiz. vel.> <lat. vel>"
print " ank | full | RR | ank | full | RR | ank | full | RR | "
print "RR: ratio of variance reduction: (1 - rrv_1)/(1 - rrv_2)"
print "=" * 50
for snr in range(4):
print ["{:1.4f}, {:1.4f}, {:1.4f}".format(median(ws11.vr_step_dx_ank[snr][:, dim]),
median(ws11.vr_step_dx_full[snr][:, dim]),
(1. - median(ws11.vr_step_dx_ank[snr][:, dim]))/(1 - median(ws11.vr_step_dx_full[snr][:, dim])) )
for dim in range(3)]
In [9]:
ws11.vr_step_dx_ank[snr].shape
Out[9]:
In [10]:
len(ws11.vr_step_dx_ank)
Out[10]:
In [ ]: