Show the versions of all the software that is installed into the virutalenv.

In [1]:
!pip freeze
Cython==0.19.1
-e [email protected]:moorepants/[email protected]#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
In [2]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [3]:
rcParams['figure.figsize'] = 16, 8
In [4]:
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.

In [5]:
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.

In [6]:
obinna
Out[6]:
<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.

In [7]:
[col for col in obinna.columns if col.endswith('.Ang') or col.endswith('ForY') or col.endswith('.Mom')]
Out[7]:
['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.

In [8]:
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.

In [9]:
obinna['FP2.ForY'].plot()
Out[9]:
<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.

In [10]:
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.

In [11]:
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]
Out[11]:
<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.

In [12]:
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.

In [13]:
right_steps = data.split_at('right', num_samples=30)
data.plot_steps('FP2.ForY', 'RKneeFlexion.Ang', 'RKneeFlexion.Rate', 'RKneeFlexion.Mom', marker='o')
Out[13]:
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:

In [14]:
mean_right_steps = right_steps.mean(axis='items')
mean_right_steps[['FP2.ForY', 'RKneeFlexion.Ang', 'RKneeFlexion.Rate', 'RKneeFlexion.Mom']].plot(subplots=True)
Out[14]:
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:

In [15]:
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)$:

In [16]:
gains, controls_star, variance, gain_var, control_var, estimated_controls = solver.solve()
In [17]:
gains.shape
Out[17]:
(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.

In [18]:
axes = solver.plot_gains(gains, gain_var)

The variance in the fit is:

In [19]:
variance
Out[19]:
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.

In [20]:
solver.plot_estimated_vs_measure_controls(estimated_controls, variance)
In [21]:
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)
In [22]:
gain_omission_matrix
Out[22]:
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)
In [23]:
axes = solver.plot_gains(gains, gain_var)