Show the versions of all the software that is installed into the virutalenv.
!pip freeze
Cython==0.19.1 -e git+git@github.com:moorepants/DynamicistToolKit.git@fa8c1ff510918ba199b993007c07967c5a5525cb#egg=DynamicistToolKit-dev Jinja2==2.7.1 MarkupSafe==0.18 Pygments==1.6 Sphinx==1.2b1 argparse==1.2.1 coverage==3.6 docutils==0.11 ipython==1.0.0 matplotlib==1.3.0 nose==1.3.0 numpy==1.7.1 numpydoc==0.4 pandas==0.12.0 pyparsing==2.0.1 python-dateutil==2.1 pytz==2013d pyzmq==13.1.0 scipy==0.12.0 six==1.4.1 sympy==0.7.3 tornado==3.1 wsgiref==0.1.2
%pylab inline
Populating the interactive namespace from numpy and matplotlib
import pandas
from numpy import deg2rad
from dtk import walk, process
Obinna provided me with some normal walking data. We exported the data from the D-Flow Motion Capture Module and included the HMB outputs. D-Flow sets missing values from HBM data to 0.000000
in the CSV file so I replace them with NA
when loading. I also set the resulting data frame's index to the TimeStamp
column.
obinna = pandas.read_csv('../data/obinna-walking.txt', delimiter='\t',
index_col="TimeStamp", na_values='0.000000')
The data frame has quite a few time series which are all floats except for the FrameNumber
column which is an integer.
obinna
<class 'pandas.core.frame.DataFrame'> Index: 5559 entries, 534.133513 to 589.712359 Columns: 617 entries, FrameNumber to HBM.COM.Z dtypes: float64(616), int64(1)
For this analysis we are only interested in the vertical ground reaction forces, the joint angles, the joint angular rates, and the joint torques. The D-Flow HBM computations give the angles and the torques. We will compute the rates ourselves. Here is a list of all of the vertical GRF, joint angle, and joint moment data.
[col for col in obinna.columns if col.endswith('.Ang') or col.endswith('ForY') or col.endswith('.Mom')]
['FP1.ForY', 'FP2.ForY', 'PelvisX.Ang', 'PelvisX.Mom', 'PelvisY.Ang', 'PelvisY.Mom', 'PelvisZ.Ang', 'PelvisZ.Mom', 'PelvisYaw.Ang', 'PelvisYaw.Mom', 'PelvisForwardPitch.Ang', 'PelvisForwardPitch.Mom', 'PelvisRightRoll.Ang', 'PelvisRightRoll.Mom', 'TrunkFlexion.Ang', 'TrunkFlexion.Mom', 'TrunkRightBend.Ang', 'TrunkRightBend.Mom', 'TrunkLeftTwist.Ang', 'TrunkLeftTwist.Mom', 'HeadFlexion.Ang', 'HeadFlexion.Mom', 'HeadRightBend.Ang', 'HeadRightBend.Mom', 'HeadLeftTwist.Ang', 'HeadLeftTwist.Mom', 'RShoulderUp.Ang', 'RShoulderUp.Mom', 'LShoulderUp.Ang', 'LShoulderUp.Mom', 'RShoulderForward.Ang', 'RShoulderForward.Mom', 'LShoulderForward.Ang', 'LShoulderForward.Mom', 'RShoulderInward.Ang', 'RShoulderInward.Mom', 'LShoulderInward.Ang', 'LShoulderInward.Mom', 'RShoulderFlexion.Ang', 'RShoulderFlexion.Mom', 'LShoulderFlexion.Ang', 'LShoulderFlexion.Mom', 'RShoulderAbduction.Ang', 'RShoulderAbduction.Mom', 'LShoulderAbduction.Ang', 'LShoulderAbduction.Mom', 'RShoulderInternalRotation.Ang', 'RShoulderInternalRotation.Mom', 'LShoulderInternalRotation.Ang', 'LShoulderInternalRotation.Mom', 'RElbowFlexion.Ang', 'RElbowFlexion.Mom', 'LElbowFlexion.Ang', 'LElbowFlexion.Mom', 'RForeArmPronation.Ang', 'RForeArmPronation.Mom', 'LForeArmPronation.Ang', 'LForeArmPronation.Mom', 'RWristFlexion.Ang', 'RWristFlexion.Mom', 'LWristFlexion.Ang', 'LWristFlexion.Mom', 'RHandAbduction.Ang', 'RHandAbduction.Mom', 'LHandAbduction.Ang', 'LHandAbduction.Mom', 'RHipFlexion.Ang', 'RHipFlexion.Mom', 'LHipFlexion.Ang', 'LHipFlexion.Mom', 'RHipAbduction.Ang', 'RHipAbduction.Mom', 'LHipAbduction.Ang', 'LHipAbduction.Mom', 'RHipInternalRotation.Ang', 'RHipInternalRotation.Mom', 'LHipInternalRotation.Ang', 'LHipInternalRotation.Mom', 'RKneeFlexion.Ang', 'RKneeFlexion.Mom', 'LKneeFlexion.Ang', 'LKneeFlexion.Mom', 'RAnklePlantarFlexion.Ang', 'RAnklePlantarFlexion.Mom', 'LAnklePlantarFlexion.Ang', 'LAnklePlantarFlexion.Mom', 'RFootPronation.Ang', 'RFootPronation.Mom', 'LFootPronation.Ang', 'LFootPronation.Mom', 'RToeFlexion.Ang', 'RToeFlexion.Mom', 'LToeFlexion.Ang', 'LToeFlexion.Mom']
The angles are in degrees, so lets convert them to radians.
for col in obinna.columns:
if col.endswith('.Ang'):
obinna[col] = deg2rad(obinna[col])
The vertical ground reaction forces are stored as variables with the extension .ForY
and the right and left plates are designated by FP2
and FP1
respectively. The following plot shows the right foot vertical ground reaction force in Newtons.
rcParams['figure.figsize'] = 14, 8
obinna['FP2.ForY'].plot()
<matplotlib.axes.AxesSubplot at 0x5105290>
We will first use the ground reaction forces to find the heel strike and toe off times during gait. I'll select a portion of the data which has suitable GRF profiles for extracting the landmarks.
start = 500
stop = 3500
data = walk.WalkingData(obinna.iloc[start:stop].copy())
Now I compute all of the joint angular rates numerically with central differencing.
to_diff = [col for col in obinna.columns if col.endswith('.Ang')]
new_names = [name.split('.')[0] + '.Rate' for name in to_diff]
data.time_derivative(to_diff, new_names)
data.raw_data[new_names]
<class 'pandas.core.frame.DataFrame'> Index: 3000 entries, 539.133881 to 569.12273 Data columns (total 46 columns): PelvisX.Rate 2882 non-null values PelvisY.Rate 2882 non-null values PelvisZ.Rate 2882 non-null values PelvisYaw.Rate 2882 non-null values PelvisForwardPitch.Rate 2882 non-null values PelvisRightRoll.Rate 2882 non-null values TrunkFlexion.Rate 2882 non-null values TrunkRightBend.Rate 2882 non-null values TrunkLeftTwist.Rate 2882 non-null values HeadFlexion.Rate 2882 non-null values HeadRightBend.Rate 2882 non-null values HeadLeftTwist.Rate 2882 non-null values RShoulderUp.Rate 2882 non-null values LShoulderUp.Rate 2882 non-null values RShoulderForward.Rate 2882 non-null values LShoulderForward.Rate 2882 non-null values RShoulderInward.Rate 2882 non-null values LShoulderInward.Rate 2882 non-null values RShoulderFlexion.Rate 2882 non-null values LShoulderFlexion.Rate 2882 non-null values RShoulderAbduction.Rate 2882 non-null values LShoulderAbduction.Rate 2882 non-null values RShoulderInternalRotation.Rate 2882 non-null values LShoulderInternalRotation.Rate 2882 non-null values RElbowFlexion.Rate 2882 non-null values LElbowFlexion.Rate 2882 non-null values RForeArmPronation.Rate 2882 non-null values LForeArmPronation.Rate 2882 non-null values RWristFlexion.Rate 2882 non-null values LWristFlexion.Rate 2882 non-null values RHandAbduction.Rate 2882 non-null values LHandAbduction.Rate 2882 non-null values RHipFlexion.Rate 2882 non-null values LHipFlexion.Rate 2882 non-null values RHipAbduction.Rate 2882 non-null values LHipAbduction.Rate 2882 non-null values RHipInternalRotation.Rate 2882 non-null values LHipInternalRotation.Rate 2882 non-null values RKneeFlexion.Rate 2882 non-null values LKneeFlexion.Rate 2882 non-null values RAnklePlantarFlexion.Rate 2882 non-null values LAnklePlantarFlexion.Rate 2882 non-null values RFootPronation.Rate 2882 non-null values LFootPronation.Rate 2882 non-null values RToeFlexion.Rate 2882 non-null values LToeFlexion.Rate 2882 non-null values dtypes: float64(46)
The strikes and toe offs can then be computed for this data.
rcParams['figure.figsize'] = 14, 10
right_strikes, left_strikes, right_offs, left_offs = \
data.grf_landmarks('FP2.ForY', 'FP1.ForY', do_plot=True, threshold=28.0)
Now we can section the gait cycle based on the right foot and interpolate at 10 distinct equally spaced points in the gait cycle.
rcParams['figure.figsize'] = 14, 10
right_steps = data.split_at('right', num_samples=15)
data.plot_steps('FP2.ForY', 'RKneeFlexion.Ang', 'RKneeFlexion.Rate', 'RKneeFlexion.Mom', marker='o')
array([<matplotlib.axes.AxesSubplot object at 0x50f0910>, <matplotlib.axes.AxesSubplot object at 0x4568690>, <matplotlib.axes.AxesSubplot object at 0x4294c50>, <matplotlib.axes.AxesSubplot object at 0x42b0e10>], dtype=object)
We can see the mean gait cycle easily:
mean_right_steps = right_steps.mean(axis='items')
mean_right_steps[['FP2.ForY', 'RKneeFlexion.Ang', 'RKneeFlexion.Rate', 'RKneeFlexion.Mom']].plot(subplots=True)
array([<matplotlib.axes.AxesSubplot object at 0x423fa90>, <matplotlib.axes.AxesSubplot object at 0x50f8590>, <matplotlib.axes.AxesSubplot object at 0x4932110>, <matplotlib.axes.AxesSubplot object at 0x494b7d0>], dtype=object)
At this point we have the data in the form needed to identify a simple linear controller where we assume the sensors are fed back and subtract from a set point in the walking gait cycle, to give the error in the gait. This error vector is then multiplied by a gain matrix and added to the limit cycle control torques to power the plant in walking.
The following gives the equation that generates the controls given the error in the sensors at a single time step during a single foot step. $t$ here represents the time in the gait cycle, not time across all steps.
$$m_m(t) = m(t)_0 + K(t) s_e(t) = m(t)_0 + K(t) [s_0(t) - s_m(t)]$$Now rearrange the equations such that we have one linear in both the gains, $K(t)$, and in $m^*(t)$:
$$m_m(t) = m(t)_0 + K(t) s_0(t) - K(t) s_m(t) = m^*(t) - K(t) s_m(t)$$Here I specify the hypothesized sensors and controls, and the data to initialized the solver:
sensors = ['RKneeFlexion.Ang',
'RKneeFlexion.Rate',
'RHipFlexion.Ang',
'RHipFlexion.Rate',
'LKneeFlexion.Ang',
'LKneeFlexion.Rate',
'LHipFlexion.Ang',
'LHipFlexion.Rate']
controls = ['RKneeFlexion.Mom',
'RHipFlexion.Mom',
'LKneeFlexion.Mom',
'LHipFlexion.Mom']
solver = walk.SimpleControlSolver(right_steps, sensors, controls)
Now we can solve the linear system for the gains at each time step and the $m^*(t)$:
gains, controls_star, variance, gain_var, control_var, estimated_controls = solver.solve()
gains.shape
(15, 4, 8)
The following plot shows the identified gains as a function of the percentage of the gait cycle. The rows correspond to the controller outputs and the columns to the measured sensors. The shaded region depicts the standard deviation in the estimates.
axes = solver.plot_gains(gains, gain_var)
The variance in the fit is:
variance
array([ 22.0162182])
And the following plot shows the measured controls at each point in the gain cycle (blue line) and the controller generated controls given the measured sensors (green dots). The error bars on the green dots show the variance in the fit.
rcParams['figure.figsize'] = 14, 14
solver.plot_estimated_vs_measure_controls(estimated_controls, variance)
array([<matplotlib.axes.AxesSubplot object at 0x3c72850>, <matplotlib.axes.AxesSubplot object at 0x3e09d50>, <matplotlib.axes.AxesSubplot object at 0x70deb50>, <matplotlib.axes.AxesSubplot object at 0x7cdb4d0>, <matplotlib.axes.AxesSubplot object at 0x7b64790>, <matplotlib.axes.AxesSubplot object at 0x7cc5b10>, <matplotlib.axes.AxesSubplot object at 0x85b9bd0>, <matplotlib.axes.AxesSubplot object at 0x85dc990>], dtype=object)
To get an idea of how each of the controller inputs affects the outputs I plot the contribution from each sensor error and gain to the different controls. To do this I make the assumption that the commanded angles and rates of the joints are the mean over a set of steps.
$$s_0(t) = \bar{s}(t)$$Then I can solve for $m_0(t)$:
$$m_0(t) = m(t) - K(t) [\bar{s}(t) - s(t)]$$And the contribution of $K(t) [\bar{s}(t) - s(t)]$ to each of the controls $m(t)$ can be computed:
$$m_1(t) = m_{0_1}(t) + K_{11} s_{e_1} + \ldots + + K_{1p} s_{e_p}$$I plot the contributions from each sensor error and the nominal moments as stacked bar charts to get an idea of the contributions at each step. The first four steps are shown.
The line plots below the bar charts show the same contributions but as a mean with respect to the step. Notice that the mean of the controller contributions is zero. This seems odd to me and believe it is an artifact of assuming the nominal sensors to be the mean of the sensors over the steps.
I was expecting to see patterns in the controller behavior across the steps but if the majority of the control moments are due to $m_0(t)$ then the controller is just responsible for applying any moments due to deviations from $m_0(t)$.
solver.plot_control_contributions(estimated_controls)
gain_omission_matrix = np.ones((len(controls), len(sensors))).astype(bool)
gain_omission_matrix[2:, :4] = False
gain_omission_matrix[:2, 4:] = False
gains, controls, variance, gain_var, control_var, estimated_controls = solver.solve(gain_omission_matrix=gain_omission_matrix)
gain_omission_matrix
array([[ True, True, True, True, False, False, False, False], [ True, True, True, True, False, False, False, False], [False, False, False, False, True, True, True, True], [False, False, False, False, True, True, True, True]], dtype=bool)
rcParams['figure.figsize'] = 14, 10
axes = solver.plot_gains(gains, gain_var)
rcParams['figure.figsize'] = 14, 14
solver.plot_estimated_vs_measure_controls(estimated_controls, variance)
array([<matplotlib.axes.AxesSubplot object at 0x4a1cbd0>, <matplotlib.axes.AxesSubplot object at 0xe7ec490>, <matplotlib.axes.AxesSubplot object at 0x692fad0>, <matplotlib.axes.AxesSubplot object at 0x67bc650>, <matplotlib.axes.AxesSubplot object at 0x8156710>, <matplotlib.axes.AxesSubplot object at 0x6910890>, <matplotlib.axes.AxesSubplot object at 0x8149490>, <matplotlib.axes.AxesSubplot object at 0x6485150>], dtype=object)
rcParams['figure.figsize'] = 14, 14
solver.plot_control_contributions(estimated_controls)
Ton also gave me some data he collected for subject doing normal walking. Here is explanation:
Here is some data that was collected for a paper which is about to be published. Proofs and supplementary documentation are >attached.
This is from a 77 kg male, 30 seconds of walking preceded by 1 second of standing. The treadmill speed increases gradually >in the beginning, at t=6.0 it is at its final constant speed to best to ignore anything that happens before that time.
The dof file contains the generalized coordinates, dofvel the generalized velocities. The jointmoment file the corresponding >generalized forces. First line has the names of the generalized coordinates. The grf file has the ground reaction forces, >which you would use to determine where you are in the gait cycle (from heelstrike to heelstrike).
CAREN_Trial04.txt is the raw data (mocap and force plate). You could actually load that file in D-Flow to produce the >variables that are on the other files, if you wanted to.
I'll do the same analysis I did for the previous data but limit everything to the right leg because the left leg force plate measurements are botched.
coordinates = pandas.read_csv('../data/CAREN_Trial04_dof.txt', delimiter='\t', index_col="TimeStamp")
speeds = pandas.read_csv('../data/CAREN_Trial04_dofvel.txt', delimiter='\t', index_col="TimeStamp")
generalized_forces = pandas.read_csv('../data/CAREN_Trial04_jointmoment.txt', delimiter='\t', index_col="TimeStamp")
ground_reaction_forces = pandas.read_csv('../data/CAREN_Trial04_grf.txt', delimiter='\t', index_col="TimeStamp")
raw_data = pandas.read_csv('../data/CAREN_Trial04.txt', delimiter='\t', index_col="TimeStamp")
Load in the separate data files and rename the columns to reflate the type of data so they can be merged into one big data frame.
coordinates.columns = [name + '.Ang' for name in coordinates.columns]
speeds.columns = [name + '.Rate' for name in speeds.columns]
generalized_forces.columns = [name + '.Mom' for name in generalized_forces.columns]
all_data = coordinates.join(speeds).join(generalized_forces).join(ground_reaction_forces)
Take the data after 15 seconds to avoid the error in the right force plate data.
data = walk.WalkingData(all_data[all_data.index > 15.0])
rcParams['figure.figsize'] = 14, 10
right_strikes, left_strikes, right_offs, left_offs = \
data.grf_landmarks('RFy','LFy', do_plot=True, threshold=5.0)
right_steps = data.split_at('right', num_samples=15)
axes = data.plot_steps('RFy', 'RKneeFlexion.Ang', 'RKneeFlexion.Rate', 'RKneeFlexion.Mom', marker='o')
sensors = ['RKneeFlexion.Ang',
'RKneeFlexion.Rate',
'RHipFlexion.Ang',
'RHipFlexion.Rate']
controls = ['RKneeFlexion.Mom',
'RHipFlexion.Mom']
solver = walk.SimpleControlSolver(right_steps, sensors, controls)
gains, controls, variance, gain_var, control_var, estimated_controls = solver.solve()
variance
array([ 7.55235945])
It is interesting to note that the gains here are of a very different scale than the previous data.
rcParams['figure.figsize'] = 14, 10
axes = solver.plot_gains(gains, gain_var)
solver.plot_estimated_vs_measure_controls(estimated_controls, variance)
array([<matplotlib.axes.AxesSubplot object at 0xf1b9f90>, <matplotlib.axes.AxesSubplot object at 0x6ddc550>, <matplotlib.axes.AxesSubplot object at 0x11562390>, <matplotlib.axes.AxesSubplot object at 0x7948c90>], dtype=object)
rcParams['figure.figsize'] = 14, 14
solver.plot_control_contributions(estimated_controls)