This notebook is intended to address this question:
Does the data processing introduce a bias? The joint moment calculation makes use of marker coordinates. The joint angle calculation uses the same marker coordinates. If there was a random measurement error in a marker, it would affect both and this may result in correlations that show up in the sys id results. We should be able to test this effect and conclude whether it is a problem. Here is an idea: (1) make an average (periodic and unperturbed) gait cycle from raw data (marker coordinates and GRF) and repeat it over and over to create the same amount of data that you usually work with. (2) add small random numbers to all samples, (3) compute angles and torques in the same way you normally do, (4) do the sys id and see what comes out. If this result is the same as what you get from the perturbation experiment, that would be bad news...
import sys
sys.path.append('..')
import operator
import numpy as np
import pandas
import matplotlib.pyplot as plt
from gaitanalysis.motek import markers_for_2D_inverse_dynamics
from gaitanalysis.gait import WalkingData
from gaitanalysis.controlid import SimpleControlSolver
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
Pick a trial.
trial_number = '018'
utils.remove_precomputed_data(utils.tmp_data_dir(), trial_number)
Temporary data directory is set to ../data
First, extract the data from the first normal walking period in the trial.
walking_event_data_frame, meta_data, 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 Cleaning the data. Saving cleaned data: ../data/cleaned-data-018-first-normal-walking.h5 4.92 s
/home/moorepants/anaconda/envs/walk/lib/python2.7/site-packages/tables/path.py:100: NaturalNameWarning: object name is not a valid Python identifier: 'First Normal Walking'; it does not match the pattern ``^[a-zA-Z_][a-zA-Z0-9_]*$``; you will not be able to use natural naming to access this object; using ``getattr()`` will still work, though NaturalNameWarning)
walking_data, walking_data_path = \
utils.write_inverse_dynamics_to_disk(walking_event_data_frame,
meta_data,
walking_event_data_path)
Computing the inverse dynamics. Saving inverse dynamics to ../data/walking-data-018-first-normal-walking.h5. 31.00 s
params = settings[trial_number]
params
(10.0, 30.0, 75, 120)
steps, walking_data = utils.section_signals_into_steps(walking_data, walking_data_path,
filter_frequency=params[0],
threshold=params[1],
num_samples_lower_bound=params[2],
num_samples_upper_bound=params[3],
num_samples=None)
Finding the ground reaction force landmarks. 0.02 s Spliting the data into steps. 2.43 s
steps.shape
(60, 93, 221)
figsize(10, 10)
axes = walking_data.step_data.hist()
Now extract one good step from the bunch. The 21st step looks as good as any:
steps.iloc[20]['FP2.ForY'].plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f7983ea0310>
This gives a list of measurements needed for the 2D inverse dynamics caluclations and also a list of just the marker labels.
measurements = markers_for_2D_inverse_dynamics('full')
measurement_list = reduce(operator.add, measurements)
marker_list = reduce(operator.add, measurements[:2])
Copy the step so it no longer points to the original data arrays.
step = steps.iloc[20][measurement_list + ['Original Time']].copy()
This code fits an nth order Fourier series to each signal needed for the 2D dynamics and then builds a smoothed step data frame.
smoothed_step = step.copy()
time = step['Original Time'].values
freq = 1.0 / (time[-1] - time[0]) # cycle / sec
print('Step frequency = {:0.2f} hz'.format(freq))
omega = 2.0 * np.pi * freq # rad /sec
print('Step frequency = {:0.2f} rad/s'.format(omega))
fourier_order = 20
initial_coeff = np.ones(1 + 2 * fourier_order)
for measurement in measurement_list:
signal = step[measurement].values
popt, pcov = utils.fit_fourier(time, signal, initial_coeff, omega)
eval_fourier = utils.fourier_series(omega)
smoothed_step[measurement] = eval_fourier(time, *popt)
Step frequency = 1.06 hz Step frequency = 6.69 rad/s
This plot shows the Fourier series (green line) vs the original data (blue dots).
figsize(14, 40)
axes = step.plot(style='.', subplots=True)
for i, ax in enumerate(axes):
smoothed_step.iloc[:, i].plot(ax=ax)
Now build some fake data by repeating the smoothed steps..
fakedata = {}
for i in range(500):
fakedata[i] = smoothed_step.copy()
fake_data_df = pandas.concat(fakedata, ignore_index=True)
fake_data_df.index = np.linspace(0.0, len(fake_data_df) * 0.01 - 0.01, len(fake_data_df))
Now add random Gaussian noise with a mean of zero and standard deviation of 5 mm to all of the marker measurements.
shape = fake_data_df[marker_list].shape
fake_data_df[marker_list] += np.random.normal(scale=0.005, size=shape)
This is a plot of the first ~11 gait cycles of the fake data (includes the noise on the markers).
figsize(14, 28)
axes = fake_data_df[['FP2.ForY'] + marker_list][:1000].plot(subplots=True)
Now compute the inverse dynamics on the noisy data.
walking_data = WalkingData(fake_data_df)
subject_mass = float(meta_data['subject']['mass'])
args = list(measurements) + [subject_mass, 6.0]
df = walking_data.inverse_dynamics_2d(*args)
This shows an example of the joint angles, rates, and torques.
figsize(14, 24)
subset = ['FP2.ForY',
'Right.Ankle.PlantarFlexion.Moment',
'Right.Ankle.PlantarFlexion.Angle',
'Right.Ankle.PlantarFlexion.Rate']
axes = walking_data.raw_data[subset][:1000].plot(subplots=True)
Now, to identify the controller the steps need to be found.
null = walking_data.grf_landmarks('FP2.ForY', 'FP1.ForY', threshold=30.0)
fake_steps = walking_data.split_at('right', num_samples=20)
The only variation in the steps is due to the marker coordinates noise.
figsize(12, 12)
axes = walking_data.plot_steps(*subset)
Now identify a joint centric gain scheduled controller as usual.
sensors, controls = utils.load_sensors_and_controls()
solver = SimpleControlSolver(fake_steps, 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
result = solver.solve(gain_omission_matrix=gain_omission_matrix)
fig, axes = utils.plot_joint_isolated_gains(sensors, controls, result[0], result[3])
Generating gain plot. 0.38 s
Interesting that we see a curve similar in shape, but not magnitude, in the angle to moment gain curves. Below I calculate the ratio of the right vertical ground reaction force to the gain curve and then plot the gain curves, normalized to by this ratio with respect to the vertical ground reaction force.
for m, a, seg in zip([0, 1, 2], [0, 2, 4], ['Ankle', 'Knee', 'Hip']):
ratios = fake_steps.iloc[0]['FP2.ForY'] / result[0][:, m, a]
ratio = ratios[:11].mean()
print('GRF to {} Gain ratio = {:0.3f}'.format(seg, ratio))
plt.plot(ratio * result[0][:, m, a], label=seg)
plt.plot(fake_steps.iloc[0]['FP2.ForY'], label='ForY')
plt.legend()
GRF to Ankle Gain ratio = 22.384 GRF to Knee Gain ratio = 7.630 GRF to Hip Gain ratio = 3.779
<matplotlib.legend.Legend at 0x7f797e466310>
!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:37:33 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.2.1 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