Ton simulated the 7-link planar walking model with an open loop control solution + a constant gain feedback controller for 600 seconds. He then provided data files with the time histories of the angles, angular rates, and ground reaction forces. I load the data into the system identification pipeline we are using for the real walking data and if the method works it shoud spit out the constant gain values that Ton chose for his simulation. He chose -300 for all the angle gains and -50 for all the angular rate gains (note his gain signs are likley defined opposite of mine so we should be looking for 300 and 50).
%matplotlib inline
from IPython.core.pylabtools import figsize
figsize(12, 8)
from IPython.display import Image
import matplotlib.pyplot as plt
import sys
sys.path.append('..')
import numpy as np
import pandas
from src import utils
from gaitanalysis.gait import GaitData
from gaitanalysis.controlid import SimpleControlSolver
Load the simulation data.
df = pandas.read_csv('../data/simulated-data-with-controller-noise-and-perturbations-600-seconds.txt',
delimiter='\t', index_col='TimeStamp')
Now change the column headers from Ton's names to the names I've been using.
df.rename(columns=utils.simulated_data_header_map(), inplace=True)
Now split the data into gait cycles and sample at 20 samples per gait cycle.
walk_data = GaitData(df)
null = walk_data.grf_landmarks('FP2.ForY', 'FP1.ForY', threshold=3.0)
cycles = walk_data.split_at('right', num_samples=20)
axes = walk_data.plot_gait_cycles('FP2.ForY', 'Right.Ankle.PlantarFlexion.Angle')
axes = walk_data.gait_cycle_stats.hist()
The preset gait cycle duration is 1.1396 and we can see that is recovered as the mean gait cycle duration.
walk_data.gait_cycle_stats.mean()
Number of Samples 114.961977 Stride Duration 1.139620 Stride Frequency 0.877701 dtype: float64
Now solve for a joint isolated controller, like we've done for the real walking data.
sensors, controls = utils.load_sensors_and_controls()
solver = SimpleControlSolver(cycles.iloc[:400], sensors, controls, cycles.iloc[400:])
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
result = solver.solve(gain_inclusion_matrix=gain_inclusion_matrix)
fig, axes = utils.plot_joint_isolated_gains(sensors, controls, result[0], result[3])
Generating gain plot. 0.33 s
From the plots, the gains seem to be in the range of the values we expected. You can look at it and probably make yourself belive that it is getting the correct answer. But the variablity is very high and the mean's aren't that close. But they seem to be in the ballpark. These are positive because I think I define my feedback term opposite sign as Ton does.
np.set_printoptions(precision=0)
print(result[0].mean(axis=0))
[[ 203. 31. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 159. 27. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 188. 37. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 182. 32. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 163. 29. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 199. 36.]]
print(result[0].std(axis=0))
[[ 321. 16. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 259. 17. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 200. 16. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 330. 17. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 280. 16. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 187. 17.]]
total_time = float(df.index[-1] - df.index[0])
sample_rate = len(df) / total_time
num_samples_in_cycle = int(sample_rate * 1.1396)
num_samples_in_cycle = 115
num_cycles = 525
interp_percent_gait = np.linspace(0.0, 0.95, num=20)
cycles_panel = {}
df['Original Time'] = df.index.values.astype(float)
cycle_columns = df.columns
from scipy.interpolate import interp1d
for i in range(num_cycles - 10):
cycle_values = df.iloc[i * num_samples_in_cycle : num_samples_in_cycle * (i + 1)].values
percent_gait = np.linspace(0.0, 1.0 * (1.0 - 1.0 / cycle_values.shape[0]), cycle_values.shape[0])
f = interp1d(percent_gait, cycle_values, axis=0)
cycles_panel[i] = pandas.DataFrame(f(interp_percent_gait),
index=interp_percent_gait,
columns=cycle_columns)
cycles_panel = pandas.Panel(cycles_panel)
np.set_printoptions(precision=5)
percent_gait
array([ 0. , 0.0087 , 0.01739, 0.02609, 0.03478, 0.04348, 0.05217, 0.06087, 0.06957, 0.07826, 0.08696, 0.09565, 0.10435, 0.11304, 0.12174, 0.13043, 0.13913, 0.14783, 0.15652, 0.16522, 0.17391, 0.18261, 0.1913 , 0.2 , 0.2087 , 0.21739, 0.22609, 0.23478, 0.24348, 0.25217, 0.26087, 0.26957, 0.27826, 0.28696, 0.29565, 0.30435, 0.31304, 0.32174, 0.33043, 0.33913, 0.34783, 0.35652, 0.36522, 0.37391, 0.38261, 0.3913 , 0.4 , 0.4087 , 0.41739, 0.42609, 0.43478, 0.44348, 0.45217, 0.46087, 0.46957, 0.47826, 0.48696, 0.49565, 0.50435, 0.51304, 0.52174, 0.53043, 0.53913, 0.54783, 0.55652, 0.56522, 0.57391, 0.58261, 0.5913 , 0.6 , 0.6087 , 0.61739, 0.62609, 0.63478, 0.64348, 0.65217, 0.66087, 0.66957, 0.67826, 0.68696, 0.69565, 0.70435, 0.71304, 0.72174, 0.73043, 0.73913, 0.74783, 0.75652, 0.76522, 0.77391, 0.78261, 0.7913 , 0.8 , 0.8087 , 0.81739, 0.82609, 0.83478, 0.84348, 0.85217, 0.86087, 0.86957, 0.87826, 0.88696, 0.89565, 0.90435, 0.91304, 0.92174, 0.93043, 0.93913, 0.94783, 0.95652, 0.96522, 0.97391, 0.98261, 0.9913 ])
from gaitanalysis.gait import plot_gait_cycles
axes = plot_gait_cycles(cycles_panel, 'FP2.ForY', 'Right.Ankle.PlantarFlexion.Angle')
solver = SimpleControlSolver(cycles_panel.iloc[:400], sensors, controls, cycles_panel.iloc[400:])
result = solver.solve(gain_inclusion_matrix=gain_inclusion_matrix)
fig, axes = utils.plot_joint_isolated_gains(sensors, controls, result[0], result[3])
Generating gain plot. 0.43 s
df = pandas.read_csv('../data/simulated-data-with-controller-noise-600-seconds.txt',
delimiter='\t', index_col='TimeStamp')
df.rename(columns=utils.simulated_data_header_map(), inplace=True)
walk_data = GaitData(df)
null = walk_data.grf_landmarks('FP2.ForY', 'FP1.ForY', threshold=3.0)
cycles = walk_data.split_at('right', num_samples=20)
axes = walk_data.plot_gait_cycles('FP2.ForY', 'Right.Ankle.PlantarFlexion.Angle')
axes = walk_data.gait_cycle_stats.hist()
sensors, controls = utils.load_sensors_and_controls()
solver = SimpleControlSolver(cycles.iloc[:400], sensors, controls, cycles.iloc[400:])
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
result = solver.solve(gain_inclusion_matrix=gain_inclusion_matrix)
fig, axes = utils.plot_joint_isolated_gains(sensors, controls, result[0], result[3])
Generating gain plot. 0.43 s
np.set_printoptions(precision=0)
print(result[0].mean(axis=0))
print(result[0].std(axis=0))
[[ 240. 34. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 160. 31. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 209. 39. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 209. 34. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 171. 33. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 216. 39.]] [[ 303. 16. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 266. 14. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 187. 14. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 296. 16. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 295. 15. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 180. 16.]]
df = pandas.read_csv('../data/simulated-data-with-perturbations-600-seconds.txt',
delimiter='\t', index_col='TimeStamp')
df.rename(columns=utils.simulated_data_header_map(), inplace=True)
walk_data = GaitData(df)
null = walk_data.grf_landmarks('FP2.ForY', 'FP1.ForY', threshold=3.0)
cycles = walk_data.split_at('right', num_samples=20)
axes = walk_data.plot_gait_cycles('FP2.ForY', 'Right.Ankle.PlantarFlexion.Angle')
axes = walk_data.gait_cycle_stats.hist()
sensors, controls = utils.load_sensors_and_controls()
solver = SimpleControlSolver(cycles.iloc[:400], sensors, controls, cycles.iloc[400:])
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
result = solver.solve(gain_inclusion_matrix=gain_inclusion_matrix)
fig, axes = utils.plot_joint_isolated_gains(sensors, controls, result[0], result[3])
Generating gain plot. 0.32 s
np.set_printoptions(precision=0)
print(result[0].mean(axis=0))
print(result[0].std(axis=0))
[[ 100. 4. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 133. 4. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 75. 24. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 126. 5. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 131. 4. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 105. 14.]] [[ 394. 32. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 187. 18. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 286. 22. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 374. 30. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 157. 13. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 159. 21.]]
Ton took this same data and did a simple identification on the hip joint. The only apparent difference is that he used the known gait average time to section the data into gait cycles, where as I am using the vertical ground reaction force to find the variable gait cycle length for each gait cycle. I also use 3/4 the data for the fit and 1/4 the data for the validation. He recovered the gains much better than I did:
Image('tons-hip-id.png')
sensors = ['Right.Hip.Flexion.Angle', 'Right.Hip.Flexion.Rate']
controls = ['Right.Hip.Flexion.Moment']
solver = SimpleControlSolver(cycles.iloc[:400], sensors, controls, cycles.iloc[400:])
result = solver.solve()
fig, axes = plt.subplots(3, 1)
axes[0].plot(result[1])
axes[0].set_ylabel('m*')
axes[1].plot(result[0][:, 0, 0])
axes[1].set_ylabel('kp')
axes[2].plot(result[0][:, 0, 1])
axes[2].set_ylabel('kd')
<matplotlib.text.Text at 0x7fd03af45d50>
axes = solver.plot_estimated_vs_measure_controls(result[-1], result[2])
plt.xlim((0, 500))
(0, 500)
Both of the data sets that include controller noise give similar solutions and with some stretch you can imagine then numbers are somewhat close to the true solution, but the variability throughout the gait cycle is high. The results from the data without controller noise are odd and don't seem to be close at all.
!git rev-parse HEAD
930978e889c077562c324c986c3f2069b3c8fd40
!git --git-dir=/home/moorepants/src/GaitAnalysisToolKit/.git --work-tree=/home/moorepants/src/GaitAnalysisToolKit rev-parse HEAD
694c97bab8d2bc0b6bf254c88e38329634d5873b
%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.8 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] |
IPython | 2.3.1 |
OS | Linux 3.13.0 40 generic x86_64 with debian jessie sid |
numpy | 1.9.1 |
scipy | 0.14.0 |
pandas | 0.15.1 |
matplotlib | 1.4.2 |
tables | 3.1.1 |
oct2py | 2.4.0 |
dtk | 0.3.5 |
gaitanalysis | 0.1.1 |
Mon Dec 08 15:40:23 2014 EST |
!conda list
# packages in environment at /home/moorepants/anaconda/envs/walk: # backports.ssl-match-hostname 3.4.0.2 <pip> cairo 1.12.18 0 dateutil 2.1 py27_2 dynamicisttoolkit 0.3.5 <pip> fastcache 1.0.1 py27_0 freetype 2.4.10 0 gaitanalysistoolkit 0.1.1 <pip> hdf5 1.8.13 1 ipython 2.3.1 py27_0 ipython-notebook 2.3.1 py27_0 jinja2 2.7.3 py27_1 libpng 1.5.13 1 libsodium 0.4.5 0 markupsafe 0.23 py27_0 matplotlib 1.4.2 np19py27_0 numexpr 2.3.1 np19py27_0 numpy 1.9.1 py27_0 oct2py 2.4.0 <pip> openssl 1.0.1h 1 pandas 0.15.1 np19py27_0 patsy 0.3.0 np19py27_0 pip 1.5.6 py27_0 pixman 0.26.2 0 py2cairo 1.10.0 py27_2 pyparsing 2.0.1 py27_0 pyqt 4.10.4 py27_0 pytables 3.1.1 np19py27_1 python 2.7.8 1 python-dateutil 1.5 <pip> pytz 2014.9 py27_0 pyyaml 3.11 py27_0 pyzmq 14.4.1 py27_0 qt 4.8.5 0 readline 6.2 2 scipy 0.14.0 np19py27_0 setuptools 7.0 py27_0 sip 4.15.5 py27_0 six 1.8.0 py27_0 sqlite 3.8.4.1 0 ssl_match_hostname 3.4.0.2 py27_0 statsmodels 0.6.1 np19py27_0 sympy 0.7.6 py27_0 system 5.8 1 tables 3.1.1 <pip> tk 8.5.15 0 tornado 4.0.2 py27_0 uncertainties 2.4.6.1 <pip> wsgiref 0.1.2 <pip> yaml 0.1.4 0 zeromq 4.0.4 0 zlib 1.2.7 0
!pip freeze
DynamicistToolKit==0.3.5 GaitAnalysisToolKit==0.1.1 Jinja2==2.7.3 MarkupSafe==0.23 PyYAML==3.11 backports.ssl-match-hostname==3.4.0.2 fastcache==1.0.1 ipython==2.3.1 matplotlib==1.4.2 numexpr==2.3.1 numpy==1.9.1 oct2py==2.4.0 pandas==0.15.1 patsy==0.3.0 pyparsing==2.0.1 python-dateutil==1.5 pytz==2014.9 pyzmq==14.4.1 scipy==0.14.0 six==1.8.0 statsmodels==0.6.1 sympy==0.7.6 tables==3.1.1 tornado==4.0.2 uncertainties==2.4.6.1 wsgiref==0.1.2