The purpose of this notebook is to check whether we get similar gain patterns from artificially created steps based on perfectly repeatable cyclic walking with random noise added to each measurement. It also gives some idea of how much we are perturbing the person by comparing the variation in the joint angles, rates, and torques in normal walking and perturbed walking.
import numpy as np
import pandas
from IPython.core.pylabtools import figsize
from gaitanalysis import gait
from gaitanalysis.controlid import SimpleControlSolver
import sys
sys.path.append('../src')
import utils
from gait_landmark_settings import settings
%matplotlib inline
Load the path to the directory with the experimental data.
trial = utils.Trial('019')
trial.prep_data('First Normal Walking')
Loading pre-cleaned data: /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/cleaned-data/019-first-normal-walking.h5 _write_event_data_frame_to_disk took 0.02 s to execute. Loading pre-computed inverse dynamics from /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gait-data/019-first-normal-walking.h5. _write_inverse_dynamics_to_disk took 0.06 s to execute. Loading pre-computed gait cycles from /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gait-data/019-first-normal-walking.h5. _section_into_gait_cycles took 0.06 s to execute.
trial.prep_data('Longitudinal Perturbation')
Loading pre-cleaned data: /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/cleaned-data/019-longitudinal-perturbation.h5 _write_event_data_frame_to_disk took 0.07 s to execute. Loading pre-computed inverse dynamics from /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gait-data/019-longitudinal-perturbation.h5. _write_inverse_dynamics_to_disk took 0.18 s to execute. Loading pre-computed gait cycles from /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gait-data/019-longitudinal-perturbation.h5. _section_into_gait_cycles took 0.26 s to execute.
Number of valid normal and perturbed steps.
normal_cycles, _ = trial._remove_bad_gait_cycles('First Normal Walking')
perturbed_cycles, _ = trial._remove_bad_gait_cycles('Longitudinal Perturbation')
num_normal_cycles = normal_cycles.shape[0]
num_normal_cycles
53
num_perturbed_cycles = perturbed_cycles.shape[0]
num_perturbed_cycles
460
sensors = trial.sensors
controls = trial.controls
Now compute the standard deviation of each signal across steps, to show the variation in the steps.
normal_std = normal_cycles.std(axis='items')[sensors + controls]
perturbed_std = perturbed_cycles.std(axis='items')[sensors + controls]
This table shows the ratio of the standard deviation of the joint angles in perturbed walking and the standard deviation of joint angles in unperturbed walking.
ratio_std_angle = (perturbed_std / normal_std)[[s for s in sensors if s.endswith('Angle')]]
ratio_std_angle
Right.Ankle.PlantarFlexion.Angle | Right.Knee.Flexion.Angle | Right.Hip.Flexion.Angle | Left.Ankle.PlantarFlexion.Angle | Left.Knee.Flexion.Angle | Left.Hip.Flexion.Angle | |
---|---|---|---|---|---|---|
0.00 | 1.566836 | 2.630389 | 1.477480 | 1.749671 | 1.665014 | 1.215681 |
0.05 | 2.059160 | 3.017920 | 1.464052 | 1.679368 | 1.667757 | 1.183173 |
0.10 | 2.436364 | 2.360532 | 1.371916 | 1.712240 | 1.468306 | 1.244763 |
0.15 | 1.902425 | 1.824315 | 1.201782 | 1.574920 | 1.310895 | 1.295006 |
0.20 | 1.468306 | 1.501638 | 1.176849 | 1.724208 | 1.251143 | 1.288113 |
0.25 | 1.345427 | 1.363272 | 1.424501 | 1.460129 | 1.455903 | 1.291402 |
0.30 | 1.569264 | 1.397347 | 1.677289 | 1.423371 | 1.477394 | 1.244054 |
0.35 | 1.623236 | 1.483223 | 1.795209 | 1.526252 | 1.461988 | 1.351938 |
0.40 | 1.390728 | 1.619975 | 1.913013 | 1.877501 | 1.470230 | 1.311642 |
0.45 | 1.396095 | 1.593145 | 1.881835 | 1.442670 | 1.656805 | 1.074231 |
0.50 | 1.493026 | 1.603355 | 1.719323 | 1.628578 | 1.668514 | 0.821170 |
0.55 | 1.494665 | 1.539993 | 1.442179 | 1.529402 | 1.592058 | 1.186555 |
0.60 | 1.448875 | 1.444366 | 1.451140 | 1.714218 | 1.838389 | 1.455963 |
0.65 | 1.302056 | 1.404191 | 1.395299 | 1.514150 | 1.950900 | 1.417497 |
0.70 | 1.186661 | 1.427134 | 1.384713 | 1.311796 | 1.494022 | 1.407294 |
0.75 | 1.062704 | 1.776564 | 1.296341 | 1.129233 | 1.390375 | 1.492891 |
0.80 | 1.164200 | 1.459187 | 1.200584 | 1.085766 | 1.464684 | 1.651544 |
0.85 | 1.456970 | 1.521140 | 1.323666 | 1.076937 | 1.420984 | 1.670096 |
0.90 | 1.617654 | 1.658391 | 1.560421 | 1.099987 | 1.385267 | 1.516423 |
0.95 | 1.825885 | 1.968238 | 1.608207 | 1.452898 | 1.507854 | 1.356928 |
The minimum and maximum values for each angle give and idea of how much the person is being perturbed. We see a max of 2X the normal walking and lows around 1X (i.e. no change).
ratio_std_angle.min()
Right.Ankle.PlantarFlexion.Angle 1.062704 Right.Knee.Flexion.Angle 1.363272 Right.Hip.Flexion.Angle 1.176849 Left.Ankle.PlantarFlexion.Angle 1.076937 Left.Knee.Flexion.Angle 1.251143 Left.Hip.Flexion.Angle 0.821170 dtype: float64
ratio_std_angle.max()
Right.Ankle.PlantarFlexion.Angle 2.436364 Right.Knee.Flexion.Angle 3.017920 Right.Hip.Flexion.Angle 1.913013 Left.Ankle.PlantarFlexion.Angle 1.877501 Left.Knee.Flexion.Angle 1.950900 Left.Hip.Flexion.Angle 1.670096 dtype: float64
We can also look at the rates.
ratio_std_rate = (perturbed_std / normal_std)[[s for s in sensors if s.endswith('Rate')]]
ratio_std_rate.min()
Right.Ankle.PlantarFlexion.Rate 0.936086 Right.Knee.Flexion.Rate 1.137838 Right.Hip.Flexion.Rate 1.150567 Left.Ankle.PlantarFlexion.Rate 0.925559 Left.Knee.Flexion.Rate 0.994680 Left.Hip.Flexion.Rate 0.807775 dtype: float64
ratio_std_rate.max()
Right.Ankle.PlantarFlexion.Rate 1.878168 Right.Knee.Flexion.Rate 2.496494 Right.Hip.Flexion.Rate 2.529571 Left.Ankle.PlantarFlexion.Rate 2.002819 Left.Knee.Flexion.Rate 2.944835 Left.Hip.Flexion.Rate 3.408380 dtype: float64
These also range 1X to 2X.
ratio_std_torque = (perturbed_std / normal_std)[controls]
ratio_std_torque.min()
Right.Ankle.PlantarFlexion.Moment 1.080816 Right.Knee.Flexion.Moment 1.314786 Right.Hip.Flexion.Moment 1.235382 Left.Ankle.PlantarFlexion.Moment 1.243807 Left.Knee.Flexion.Moment 1.278643 Left.Hip.Flexion.Moment 1.355717 dtype: float64
ratio_std_torque.max()
Right.Ankle.PlantarFlexion.Moment 3.089017 Right.Knee.Flexion.Moment 2.576820 Right.Hip.Flexion.Moment 2.330060 Left.Ankle.PlantarFlexion.Moment 3.180187 Left.Knee.Flexion.Moment 2.605091 Left.Hip.Flexion.Moment 2.055234 dtype: float64
The torques have higher max values, some close to 3X. The min is still around 1X.
Samin gave me some similar numbers from here results from the 2-link inverted pendulum balancing problem.
standing_std_no_external_perturbation = pandas.DataFrame({'Ankle.Angle': [0.0043],
'Ankle.Rate': [0.0122],
'Hip.Angle': [0.0102],
'Hip.Rate': [0.0405],
'Ankle.Torque': [50.11],
'Hip.Torque': [50.04]})
standing_std_min_external_perturbation = pandas.DataFrame({'Ankle.Angle': [0.0107],
'Ankle.Rate': [0.0513],
'Hip.Angle': [0.0157],
'Hip.Rate': [0.0501],
'Ankle.Torque': [54.42],
'Hip.Torque': [50.15]})
standing_std_better_external_perturbation = pandas.DataFrame({'Ankle.Angle': [0.0105],
'Ankle.Rate': [0.0277],
'Hip.Angle': [0.0118],
'Hip.Rate': [0.093],
'Ankle.Torque': [51.228],
'Hip.Torque': [50.07]})
So at the bare minimum perturbation for reasonable control id with the direct method, these are the ratios of the standard deviations:
standing_std_min_external_perturbation / standing_std_no_external_perturbation
Ankle.Angle | Ankle.Rate | Ankle.Torque | Hip.Angle | Hip.Rate | Hip.Torque | |
---|---|---|---|---|---|---|
0 | 2.488372 | 4.204918 | 1.086011 | 1.539216 | 1.237037 | 1.002198 |
So, we seem to have enough perturbations in the walking ID for some of the angles. The minimum perturbations in the walking are all a bit lower than these numbers, with the ankle angle and rate being the poorest.
We can visualize how variable the steps are too.
figsize(14, 14)
axes = trial.gait_data_objs['First Normal Walking'].plot_gait_cycles(
*(sensors[:2] + controls[:2]), mean=False)
axes = trial.gait_data_objs['Longitudinal Perturbation'].plot_gait_cycles(
*(sensors[:2] + controls[:2]), mean=False)
Now we will select one of the steps and build a panel of steps that based on that single step and add random noise to the columns for the sensors and controls.
def create_artificial_steps(base_step_number):
step_data_frames = []
for i in range(num_normal_cycles):
# pick the 20th step
df = normal_cycles.iloc[base_step_number].copy()
for col in sensors + controls:
# normal() doesn't say that you can pass in an array for the scale but I think it works by broadcasting.
df[col] = df[col] + np.random.normal(0.0, perturbed_std[col].values, (len(df),))
step_data_frames.append(df)
artificial_steps = pandas.Panel(dict(zip(range(len(step_data_frames)), step_data_frames)))
return artificial_steps
artificial_steps = create_artificial_steps(20)
axes = gait.plot_gait_cycles(artificial_steps, *(sensors[:2] + controls[:2]))
trial.identify_controller('First Normal Walking', 'joint isolated')
Loading pre-cleaned data: /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/cleaned-data/019-first-normal-walking.h5 _write_event_data_frame_to_disk took 0.02 s to execute. Loading pre-computed inverse dynamics from /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gait-data/019-first-normal-walking.h5. _write_inverse_dynamics_to_disk took 0.04 s to execute. Loading pre-computed gait cycles from /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gait-data/019-first-normal-walking.h5. _section_into_gait_cycles took 0.04 s to execute. Loading pre-computed gains from: /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gains/joint-isolated/019-first-normal-walking.npz /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gains/joint-isolated/019-first-normal-walking.h5 identify_controller took 0.18 s to execute.
figsize(10, 8)
fig, axes = trial.plot_joint_isolated_gains('First Normal Walking', 'joint isolated')
Generating gain plot. plot_joint_isolated_gains took 0.81 s to execute.
trial.identify_controller('Longitudinal Perturbation', 'joint isolated')
Loading pre-cleaned data: /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/cleaned-data/019-longitudinal-perturbation.h5 _write_event_data_frame_to_disk took 0.05 s to execute. Loading pre-computed inverse dynamics from /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gait-data/019-longitudinal-perturbation.h5. _write_inverse_dynamics_to_disk took 0.19 s to execute. Loading pre-computed gait cycles from /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gait-data/019-longitudinal-perturbation.h5. _section_into_gait_cycles took 0.19 s to execute. Loading pre-computed gains from: /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gains/joint-isolated/019-longitudinal-perturbation.npz /home/moorepants/Manuscripts/gait-control-direct-id-paper/processed-data/gains/joint-isolated/019-longitudinal-perturbation.h5 identify_controller took 0.54 s to execute.
fig, axes = trial.plot_joint_isolated_gains('Longitudinal Perturbation', 'joint isolated')
Generating gain plot. plot_joint_isolated_gains took 1.48 s to execute.
solver = SimpleControlSolver(create_artificial_steps(20), sensors, controls)
gain_inclusion_matrix = np.zeros((len(controls), len(sensors))).astype(bool)
for i, row in enumerate(gain_inclusion_matrix):
row[2 * i:2 * i + 2] = True
artificial_result = solver.solve(gain_inclusion_matrix=gain_inclusion_matrix)
fig, axes = utils.plot_joint_isolated_gains(trial.sensors, trial.controls,
artificial_result[0], artificial_result[3])
Generating gain plot. 0.33 s
Solve for the gains from the same base step but with different noise applied to ensure we get different results than the first.
solver = SimpleControlSolver(create_artificial_steps(20), sensors, controls)
artificial_result = solver.solve(gain_inclusion_matrix=gain_inclusion_matrix)
fig, axes = utils.plot_joint_isolated_gains(trial.sensors, trial.controls,
artificial_result[0], artificial_result[3])
Generating gain plot. 0.42 s
So this is good news. We are not able to identifying the same gains from this artificial data. And the gains have no apparent, repeatable patterns. It also seems to give gains of zero, which is what you may expect from this atificial data.
This seems to to show that the gain patterns we get from the real data, both normal walking and perturbed walking, are not some artifact of the nominal cyclic motion and the mathematical method. It may be the case that walking without the longitudinal perturbations is still internally perturbed enough for us to extract some reasonable gains for the controller. This also indicates that the recoveries from perturbations are likely not random, because we would not see the repeatable results in the controller identification.
!git rev-parse HEAD
6af0c229bb0a9edf78f8489d5f0efe02d5685a57
!git --git-dir=/home/moorepants/src/GaitAnalysisToolKit/.git --work-tree=/home/moorepants/src/GaitAnalysisToolKit rev-parse HEAD
9e80dfdcfe0a14b44e0ebcbadb6e9e827d215c3c
%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
Installed version_information.py. To use it, type: %load_ext version_information
%load_ext version_information
%version_information numpy, scipy, pandas, matplotlib, tables, oct2py, dtk, gaitanalysis
Software | Version |
---|---|
Python | 2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] |
IPython | 3.0.0 |
OS | Linux 3.13.0 49 generic x86_64 with debian jessie sid |
numpy | 1.9.2 |
scipy | 0.15.1 |
pandas | 0.16.0 |
matplotlib | 1.4.2 |
tables | 3.1.1 |
oct2py | 3.1.0 |
dtk | 0.4.0 |
gaitanalysis | 0.2.0dev |
Thu Apr 23 17:25:14 2015 PDT |
!conda list
# packages in environment at /home/moorepants/anaconda/envs/gait-direct-id-paper: # backports.ssl-match-hostname 3.4.0.2 <pip> cairo 1.12.18 0 dateutil 2.4.1 py27_0 dynamicisttoolkit 0.4.0 py27_0 fastcache 1.0.2 py27_0 fontconfig 2.11.1 2 freetype 2.4.10 0 gaitanalysistoolkit (/home/moorepants/src/GaitAnalysisToolKit) 0.2.0.dev0 <pip> hdf5 1.8.14 0 ipython 3.0.0 py27_0 ipython-notebook 3.0.0 py27_1 jinja2 2.7.3 py27_1 jsonschema 2.4.0 py27_0 libpng 1.5.13 1 libsodium 0.4.5 0 libxml2 2.9.0 0 markupsafe 0.23 py27_0 matplotlib 1.4.2 np19py27_0 mistune 0.5.1 py27_0 nose 1.3.4 py27_1 numexpr 2.3.1 np19py27_0 numpy 1.9.2 py27_0 oct2py 3.1.0 py27_0 openssl 1.0.1k 1 pandas 0.16.0 np19py27_1 pip 6.1.1 py27_0 pixman 0.26.2 0 ptyprocess 0.4 py27_0 py2cairo 1.10.0 py27_2 pygments 2.0.2 py27_0 pyparsing 2.0.1 py27_0 pyqt 4.10.4 py27_0 pyside 1.2.1 py27_1 pytables 3.1.1 np19py27_2 python 2.7.9 2 python-dateutil 2.4.2 py27_0 pytz 2015.2 py27_0 pyyaml 3.11 py27_0 pyzmq 14.5.0 py27_0 qt 4.8.5 0 readline 6.2 2 scipy 0.15.1 np19py27_0 seaborn 0.5.1 np19py27_0 setuptools 15.0 py27_0 shiboken 1.2.1 py27_0 sip 4.15.5 py27_0 six 1.9.0 py27_0 snakeviz 0.3.1 <pip> sqlite 3.8.4.1 1 ssl_match_hostname 3.4.0.2 py27_0 sympy 0.7.6 py27_0 system 5.8 2 tables 3.1.1 <pip> terminado 0.5 py27_0 tk 8.5.18 0 tornado 4.1 py27_0 yaml 0.1.4 0 zeromq 4.0.4 0 zlib 1.2.8 0
!pip freeze
backports.ssl-match-hostname==3.4.0.2 DynamicistToolKit==0.4.0 fastcache==1.0.2 -e git+git@github.com:csu-hmc/GaitAnalysisToolKit.git@9e80dfdcfe0a14b44e0ebcbadb6e9e827d215c3c#egg=GaitAnalysisToolKit-origin_speedup-inverse-dynamics ipython==3.0.0 Jinja2==2.7.3 jsonschema==2.4.0 MarkupSafe==0.23 matplotlib==1.4.2 mistune==0.5.1 nose==1.3.4 numexpr==2.3.1 numpy==1.9.2 oct2py==3.1.0 pandas==0.16.0 ptyprocess==0.4 Pygments==2.0.2 pyparsing==2.0.1 PySide==1.2.1 python-dateutil==2.4.2 pytz==2015.2 PyYAML==3.11 pyzmq==14.5.0 scipy==0.15.1 seaborn==0.5.1 six==1.9.0 snakeviz==0.3.1 sympy==0.7.6 tables==3.1.1 terminado==0.5 tornado==4.1