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@12538c404b0597f152f30adc2106847ef56fa589#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==2013b pyzmq==13.1.0 scipy==0.12.0 six==1.3.0 sympy==0.7.3 tornado==3.1 wsgiref==0.1.2
%pylab inline
Populating the interactive namespace from numpy and matplotlib
rcParams['figure.figsize'] = 16, 8
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.
obinna['FP2.ForY'].plot()
<matplotlib.axes.AxesSubplot at 0x3628d50>
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.
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.
right_steps = data.split_at('right', num_samples=30)
data.plot_steps('FP2.ForY', 'RKneeFlexion.Ang', 'RKneeFlexion.Rate', 'RKneeFlexion.Mom', marker='o')
array([<matplotlib.axes.AxesSubplot object at 0x4e98fd0>, <matplotlib.axes.AxesSubplot object at 0x369a950>, <matplotlib.axes.AxesSubplot object at 0x36ab7d0>, <matplotlib.axes.AxesSubplot object at 0x459e550>], 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 0x4b05fd0>, <matplotlib.axes.AxesSubplot object at 0x3637090>, <matplotlib.axes.AxesSubplot object at 0x4317b90>, <matplotlib.axes.AxesSubplot object at 0x4334290>], 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.
$$m_m(t) = m(t)_{lc} + K(t) s_e(t) = m(t)_{lc} + K(t) (s_{lc}(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)_{lc} + K(t) s_{lc}(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
(30, 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([ 20.49582551])
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.
solver.plot_estimated_vs_measure_controls(estimated_controls, variance)
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)
axes = solver.plot_gains(gains, gain_var)
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])
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=30)
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.12806411])
It is interesting to note that the gains here are of a very different scale than the previous data.
axes = solver.plot_gains(gains, gain_var)
solver.plot_estimated_vs_measure_controls(estimated_controls, variance)