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 sys
sys.path.append('..')
from src import utils
from src.grf_landmark_settings import settings
%matplotlib inline
from IPython.core.pylabtools import figsize
Load the path to the directory with the experimental data.
trials_dir = utils.trial_data_dir()
Trials data directory is set to /home/moorepants/Data/human-gait/gait-control-identification
Trial #68 is a good trial.
trial_number = '068'
First, extract the data from the first normal walking period in the trial.
normal_walking_event_data_frame, meta_data, normal_walking_event_data_path = \
utils.write_event_data_frame_to_disk(trial_number, event='First Normal Walking')
Trials data directory is set to /home/moorepants/Data/human-gait/gait-control-identification Temporary data directory is set to ../data Loading pre-cleaned data: ../data/cleaned-data-068-first-normal-walking.h5 0.08 s
normal_walking_data, normal_walking_data_path = \
utils.write_inverse_dynamics_to_disk(normal_walking_event_data_frame,
meta_data,
normal_walking_event_data_path)
Loading pre-computed inverse dynamics from ../data/walking-data-068-first-normal-walking.h5. 0.04 s
params = settings[trial_number]
normal_steps, normal_walking_data = \
utils.section_signals_into_steps(normal_walking_data,
normal_walking_data_path,
filter_frequency=params[0],
threshold=params[1],
num_samples_lower_bound=params[2],
num_samples_upper_bound=params[3])
Loading pre-computed steps from ../data/walking-data-068-first-normal-walking.h5. 0.045748
And extract the perturbed walking section from the same trial.
perturbed_event_data_frame, meta_data, perturbed_event_data_path = \
utils.write_event_data_frame_to_disk(trial_number)
Trials data directory is set to /home/moorepants/Data/human-gait/gait-control-identification Temporary data directory is set to ../data Loading pre-cleaned data: ../data/cleaned-data-068-longitudinal-perturbation.h5 0.07 s
perturbed_walking_data, perturbed_walking_data_path = \
utils.write_inverse_dynamics_to_disk(perturbed_event_data_frame, meta_data, perturbed_event_data_path)
Loading pre-computed inverse dynamics from ../data/walking-data-068-longitudinal-perturbation.h5. 0.17 s
perturbed_steps, perturbed_walking_data = \
utils.section_signals_into_steps(perturbed_walking_data,
perturbed_walking_data_path,
filter_frequency=params[0],
threshold=params[1],
num_samples_lower_bound=params[2],
num_samples_upper_bound=params[3])
Loading pre-computed steps from ../data/walking-data-068-longitudinal-perturbation.h5. 0.111402
Number of valid normal and perturbed steps.
num_normal_steps = normal_steps.shape[0]
num_normal_steps
50
perturbed_steps.shape[0]
429
sensors, controls = utils.load_sensors_and_controls()
Now compute the standard deviation of each signal across steps, to show the variation in the steps.
normal_std = normal_steps.std(axis='items')[sensors + controls]
perturbed_std = perturbed_steps.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.618319 | 2.247994 | 1.288912 | 1.926326 | 1.075319 | 1.511918 |
0.05 | 2.192838 | 1.501089 | 1.202325 | 2.552773 | 1.223832 | 1.519456 |
0.10 | 1.781349 | 1.679455 | 1.322374 | 2.916951 | 1.402097 | 1.254302 |
0.15 | 2.095995 | 1.827481 | 1.289481 | 2.780918 | 1.396114 | 1.156941 |
0.20 | 2.158633 | 1.543787 | 1.309101 | 1.085201 | 1.331573 | 1.150347 |
0.25 | 1.955547 | 1.453090 | 1.457115 | 1.279312 | 1.367076 | 1.153786 |
0.30 | 1.825708 | 1.629841 | 1.624708 | 1.287720 | 1.640474 | 1.278629 |
0.35 | 1.618211 | 1.767320 | 1.778991 | 1.205845 | 1.550288 | 1.440512 |
0.40 | 1.309958 | 1.735724 | 1.819850 | 1.065881 | 1.546013 | 1.293779 |
0.45 | 1.229039 | 1.714844 | 1.869774 | 1.130605 | 1.482120 | 1.078247 |
0.50 | 1.481159 | 1.733850 | 1.852123 | 1.231980 | 1.320667 | 0.908879 |
0.55 | 1.842465 | 1.774532 | 1.509914 | 1.305433 | 1.569905 | 0.999879 |
0.60 | 2.089173 | 1.756863 | 1.447760 | 1.975813 | 1.805059 | 1.270518 |
0.65 | 1.858328 | 1.738585 | 1.504956 | 2.374454 | 2.050551 | 1.389473 |
0.70 | 1.578017 | 1.528670 | 1.340027 | 2.533431 | 1.717937 | 1.293544 |
0.75 | 1.658846 | 1.204797 | 1.236146 | 1.987912 | 1.422549 | 1.262855 |
0.80 | 1.514951 | 1.632808 | 1.282370 | 1.673203 | 1.308990 | 1.496922 |
0.85 | 1.209626 | 1.726183 | 1.370029 | 1.489646 | 1.321812 | 1.786014 |
0.90 | 1.253828 | 1.604310 | 1.482364 | 1.694609 | 1.314824 | 1.994658 |
0.95 | 1.525031 | 1.705935 | 1.420937 | 2.176375 | 1.205337 | 1.687014 |
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.209626 Right.Knee.Flexion.Angle 1.204797 Right.Hip.Flexion.Angle 1.202325 Left.Ankle.PlantarFlexion.Angle 1.065881 Left.Knee.Flexion.Angle 1.075319 Left.Hip.Flexion.Angle 0.908879 dtype: float64
ratio_std_angle.max()
Right.Ankle.PlantarFlexion.Angle 2.192838 Right.Knee.Flexion.Angle 2.247994 Right.Hip.Flexion.Angle 1.869774 Left.Ankle.PlantarFlexion.Angle 2.916951 Left.Knee.Flexion.Angle 2.050551 Left.Hip.Flexion.Angle 1.994658 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 1.065337 Right.Knee.Flexion.Rate 1.227428 Right.Hip.Flexion.Rate 1.195008 Left.Ankle.PlantarFlexion.Rate 0.876602 Left.Knee.Flexion.Rate 0.992043 Left.Hip.Flexion.Rate 1.124041 dtype: float64
ratio_std_rate.max()
Right.Ankle.PlantarFlexion.Rate 2.293435 Right.Knee.Flexion.Rate 2.092252 Right.Hip.Flexion.Rate 2.087967 Left.Ankle.PlantarFlexion.Rate 3.148746 Left.Knee.Flexion.Rate 1.932574 Left.Hip.Flexion.Rate 1.987495 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.586617 Right.Knee.Flexion.Moment 1.124440 Right.Hip.Flexion.Moment 0.698268 Left.Ankle.PlantarFlexion.Moment 1.103677 Left.Knee.Flexion.Moment 1.006484 Left.Hip.Flexion.Moment 0.917101 dtype: float64
ratio_std_torque.max()
Right.Ankle.PlantarFlexion.Moment 2.587869 Right.Knee.Flexion.Moment 2.437394 Right.Hip.Flexion.Moment 2.819450 Left.Ankle.PlantarFlexion.Moment 2.814005 Left.Knee.Flexion.Moment 2.479371 Left.Hip.Flexion.Moment 2.008800 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.
import pandas
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 = normal_walking_data.plot_steps(*(sensors[:2] + controls[:2]), mean=False)
axes = perturbed_walking_data.plot_steps(*(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.
import numpy as np
def create_artificial_steps(base_step_number):
step_data_frames = []
for i in range(num_normal_steps):
# pick the 20th step
df = normal_steps.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)
from gaitanalysis.gait import WalkingData
artificial_walking_data = WalkingData(normal_walking_event_data_frame)
artificial_walking_data.steps = artificial_steps
axes = artificial_walking_data.plot_steps(*(sensors[:2] + controls[:2]))
sensor_labels, control_labels, normal_result, solver = \
utils.find_joint_isolated_controller(normal_steps, normal_walking_event_data_path)
Identifying the controller. Loading pre-computed gains from: ../data/joint-isolated-gain-data-068-first-normal-walking.npz ../data/joint-isolated-gain-data-068-first-normal-walking.h5 0.03 s
figsize(10, 8)
fig, axes = utils.plot_joint_isolated_gains(sensor_labels, control_labels, normal_result[0], normal_result[3])
Generating gain plot. 0.64 s
sensor_labels, control_labels, perturbed_result, solver = \
utils.find_joint_isolated_controller(perturbed_steps.iloc[:normal_steps.shape[0]], perturbed_event_data_path)
Identifying the controller. Loading pre-computed gains from: ../data/joint-isolated-gain-data-068-longitudinal-perturbation.npz ../data/joint-isolated-gain-data-068-longitudinal-perturbation.h5 0.01 s
fig, axes = utils.plot_joint_isolated_gains(sensor_labels, control_labels,
perturbed_result[0], perturbed_result[3])
Generating gain plot. 0.37 s
from gaitanalysis.controlid import SimpleControlSolver
solver = SimpleControlSolver(create_artificial_steps(20), sensors, controls)
gain_omission_matrix = np.zeros((len(controls), len(sensors))).astype(bool)
for i, row in enumerate(gain_omission_matrix):
row[2 * i:2 * i + 2] = True
artificial_result = solver.solve(gain_omission_matrix=gain_omission_matrix)
fig, axes = utils.plot_joint_isolated_gains(sensor_labels, control_labels,
artificial_result[0], artificial_result[3])
Generating gain plot. 0.38 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_omission_matrix=gain_omission_matrix)
fig, axes = utils.plot_joint_isolated_gains(sensor_labels, control_labels,
artificial_result[0], artificial_result[3])
Generating gain plot. 0.39 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
7ba68f0160c23a61204291ca107ad570ac6f6e5a
!git --git-dir=/home/moorepants/src/GaitAnalysisToolKit/.git --work-tree=/home/moorepants/src/GaitAnalysisToolKit rev-parse HEAD
a3732352747bc03ca839df9ff02ddcbd889e636d
%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 gaitanalysis, numpy, scipy, pandas, matplotlib, tables, oct2py
Software | Version |
---|---|
Python | 2.7.8 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] |
IPython | 2.3.0 |
OS | Linux 3.13.0 37 generic x86_64 with debian jessie sid |
gaitanalysis | 0.1.0dev |
numpy | 1.8.2 |
scipy | 0.14.0 |
pandas | 0.12.0 |
matplotlib | 1.4.0 |
tables | 3.1.1 |
oct2py | 2.4.0 |
Fri Oct 17 12:14:24 2014 EDT |
!pip freeze
DynamicistToolKit==0.3.5 -e git+git@github.com:csu-hmc/GaitAnalysisToolKit.git@a3732352747bc03ca839df9ff02ddcbd889e636d#egg=GaitAnalysisToolKit-origin/HEAD Jinja2==2.7.2 MarkupSafe==0.18 PyYAML==3.11 backports.ssl-match-hostname==3.4.0.2 ipython==2.3.0 matplotlib==1.4.0 numexpr==2.3.1 numpy==1.8.2 oct2py==2.4.0 pandas==0.12.0 patsy==0.3.0 pyparsing==2.0.1 python-dateutil==1.5 pytz==2014.7 pyzmq==14.3.0 scipy==0.14.0 six==1.8.0 statsmodels==0.5.0 tables==3.1.1 tornado==3.2.1 uncertainties==2.4.6.1 wsgiref==0.1.2