Structure of Sohn2013_hinlimb_models.mat

We've got a .mat file with data that we want to explore. As specified by Hongchul Sohn, the author of the data (Sohn et al 2013):

http://www.sciencedirect.com/science/article/pii/S0021929013000614

Cats

Each cell, ({1} birdy, {2} nikey, {3} russl), contains:

  • name (string): cat's name
  • R (7x31): moment arm matrix
  • J (6x7): geometric Jacobian, mapping endpoing wrench (6x1) and joint torque (7x1)
  • q (7x1): model posture matched to experimental cat configuration at normal preferred stance

DoFs {1,3}

DoF names

  • HF (hip flexion)
  • HAd (hip adduction)
  • HR (hip rotation)
  • KE (knee extension)
  • KA (knee adduction)
  • AE (ankle extension)
  • AAd (ankle adduction)

muscles {1x31}

muscle names (abbreviations in paper)

fmax (1x31)

maximum isometric muscle force, descending from original cat model [Burkholder and Nichols, 2004]

afl95

scaling factor from active force-length relationship curve, at 95% of optimal fiber length

cosa95 (1x31)

cosine of pennation at given posture.

Note

  • The above two variables are computed from specific Hill-type muscle model using Neuromechanic.
  • Tendon slack lengths, which we regard as free variable, are set so that all muscles are at 95% of its optimum fiber lengths.
  • Because of this, pennation is same across all postures (or cats).

Example

For model construction, you only need fmax, afl95, cos95 and R for a given cat. Choose birdy (cat = 1), for example:

RFm = Cats{1,cat}.R*diag(afl95*(fmax.*cosa95))

Linear equation for static equilibrium then becomes:

RFm*e == (Cats{1,cat}.J)'*Fend

where e (31x1) is the muscle excitation vector you are solving for and Fend (6x1) is an endpoint wrench vector you specify as the task.

Extraction using numpy & scipy


In [1]:
from scipy.io import loadmat
from numpy.linalg import pinv
from numpy import matrix, diag
mat_data_struct = loadmat('Sohn2013_hinlimb_models.mat')

We've loaded the data, but we don't know anything about it


In [2]:
type(mat_data_struct)


Out[2]:
dict

So it's a dictionary


In [3]:
mat_data_struct.keys()


Out[3]:
['afl95',
 'muscles',
 'fmax',
 '__header__',
 'cosa95',
 '__globals__',
 'Cats',
 'DoFs',
 '__version__']

Right, we are interested in the 'Cats' part (according to the file readme.m)


In [4]:
mat_cats_struct = mat_data_struct['Cats']

In [5]:
type(mat_cats_struct)


Out[5]:
numpy.ndarray

In [6]:
mat_cats_struct.shape


Out[6]:
(1L, 3L)

So it's a 3-element row-vector, just as described in the readme. Let's check out the first cat.


In [7]:
type(mat_cats_struct[0,0])


Out[7]:
numpy.ndarray

In [8]:
mat_cats_struct[0,0].shape


Out[8]:
(1L, 1L)

Ok, that's a little weird. Apparently it's a 1-element array. Whatever. Let's check it out.


In [9]:
mat_cats_struct[0,0]


Out[9]:
array([[ ([[-0.02783489912161, -0.01882911604966, -0.02868287010288, -0.03071652959358, 0.0, 0.0, 0.0, -0.00208368875446, 0.0003485538953436, 0.001456767089118, -0.02927431092181, 0.0, 0.0, 0.0, -0.009245564245646, 0.0, 0.0, 0.003935755549574, 0.0, 0.000896679288098, -0.007143540946358, 0.003742031416251, 0.02282982547408, -0.02761199801176, 0.0, -0.03582234054497, 0.0, 0.0, 0.0, 0.0, 0.0], [0.003990983832939, 0.006648899690929, -0.0006355796399535, -0.001548459275539, 0.0, 0.0, 0.0, -0.002845565498727, -0.007179741520794, -0.005218254438697, 0.005843580316898, 0.0, 0.0, 0.0, 0.001783503698058, 0.0, 0.0, 0.001548346656565, 0.0, -0.008898295416249, -0.004060503947911, -0.001227457615269, 0.0149204158403, 0.001019405134613, 0.0, -0.002623144767741, 0.0, 0.0, 0.0, 0.0, 0.0], [0.003157945210088, 0.0006708891567457, 0.002124034575455, 0.00428315647766, 0.0, 0.0, 0.0, -0.008643568860864, -0.01505518451558, -0.01550768474396, 0.001688019636424, 0.0, 0.0, 0.0, 0.0003056909509738, 0.0, 0.0, -0.001282149883146, 0.0, 0.001198455486097, 0.007721092926651, 0.0003367388778529, -0.001019161859921, -0.001814524756842, 0.0, -0.002837083717716, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, -0.03380243501439, 0.002147091016874, 0.0, 0.0, 0.0, 0.0, 0.0, -0.02812536410839, -0.007968300290057, -0.008778953315791, 0.0, 0.0, 0.0, -0.008528042676866, 0.0, 0.0, 0.0, 0.0, 0.009775322426243, 0.0, -0.004371980070236, 0.0, -0.03422503356838, 0.0, 0.0, 0.009102669381603, 0.008852135067189, 0.008996176966654], [0.0, 0.0, 0.0, 0.01241406697032, -0.005422976628189, 0.0, 0.0, 0.0, 0.0, 0.0, -0.001405103544474, -0.00410598689128, 0.001116558764212, 0.0, 0.0, 0.0, -0.005144334278681, 0.0, 0.0, 0.0, 0.0, -0.004771203141274, 0.0, -0.005108715145826, 0.0, -0.006672177313841, 0.0, 0.0, -0.003267709530142, -0.001944106483436, -0.007106260389893], [0.0, 0.0, 0.0, 0.0, -0.009731983330464, 0.001566067711558, 0.004328204415756, 0.0, 0.0, 0.0, 0.0, 0.009682440308431, 0.011401749685, 0.001190070599856, 0.0, 0.0007996029287431, 0.01459020049161, 0.0, 0.0005579586433572, 0.0, 0.0, 0.0, 0.0, 0.0, 0.01060796892582, 0.0, -0.008676970477828, 0.0002876872276216, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, -0.003332472366159, 0.00995602056047, 0.00143587047278, 0.0, 0.0, 0.0, 0.0, 0.004402845592835, 0.002849124518352, -0.002282577142891, 0.0, -0.005968695563936, 0.00507288443447, 0.0, -0.002801703213964, 0.0, 0.0, 0.0, 0.0, 0.0, 0.004311411851769, 0.0, -0.001902162143581, 0.009734985378671, 0.0, 0.0, 0.0]], [[0.2234976451696, 0.008432122451868, -0.0146817143015, 0.1471455529382, -0.003810460243922, -0.07060299184861, -0.01825125389556], [-0.02835970686496, -0.009624710493188, -0.007787963353863, -0.09014052087962, -0.001723371151422, 0.009531332381853, -0.009033533902514], [0.0, -0.149420018465, -0.1685661873535, -0.01579992548833, -0.03583576980603, 0.0007766049408956, -0.07430851346857], [0.0, 0.7521698653937, 0.6588003950711, 0.1909253966513, -0.2868878091614, -0.01283173881681, 0.9108096163933], [0.0, 0.6589692660463, -0.7519771103363, 0.1414048924364, -0.9549105796903, -0.01358570047598, -0.3722639648047], [1.0, 0.0, -0.02263768946677, 0.9713661252631, 0.0764275457546, -0.9998253723634, -0.1784527477887]], [u'birdy'], [[-48.7786895377033], [1.2971548715395], [-8.67990806283839], [-84.0721897678482], [18.6757078585793], [-31.6202992915975], [-11.8246234522819]])]], 
      dtype=[('R', 'O'), ('J', 'O'), ('name', 'O'), ('q', 'O')])

Alright, that kind of makes sense. Let's unwrap it.


In [10]:
cat_data = mat_cats_struct[0,0][0][0]
R = cat_data[0]
J = cat_data[1]
name = cat_data[2][0]
q = cat_data[3]

Let's make sure the variables match the specs:


In [11]:
print("R dimensions, {0}, are valid? {1}".format(R.shape, R.shape == (7,31)))
print("J dimensions, {0}, are valid? {1}".format(J.shape, J.shape == (6,7)))
print("This should be a string with the cat's name: {0}".format(name))
print("q dimensions, {0}, are valid? {1}".format(q.shape, q.shape == (7,1)))


R dimensions, (7L, 31L), are valid? True
J dimensions, (6L, 7L), are valid? True
This should be a string with the cat's name: birdy
q dimensions, (7L, 1L), are valid? True

Model construction for the first cat (birdy)

We've already set the variables R and J, for the moment arm matrix and Jacobian, respectively. Let's also extract the other globals mentioned above.

First, we set up the cat-global constants:


In [12]:
cat_const_fmax = mat_data_struct['fmax']
cat_const_cosa95 = mat_data_struct['cosa95']
cat_const_afl95 = mat_data_struct['afl95'][0,0]

The following is specific to birdy


In [13]:
F = (cat_const_afl95 * cat_const_fmax * cat_const_cosa95)[0,:]

The equilibrium is now given by

R*F*x == J.T * w

where x is the activation vector, and w is the corresponding endpoint wrench. Thus, we get our matrix of generators:


In [14]:
W = pinv(J.T) * matrix(R) * matrix(diag(F))

W is a 6x31 matrix, where the $i$-th column represents the contribution of the $i$-th muscle.

Putting it all together and saving everything


In [15]:
def compute_cat_generators(cat_index):
    cat_data = mat_cats_struct[0, cat_index][0][0]
    R = cat_data[0]
    J = cat_data[1]
    name = cat_data[2][0]
    q = cat_data[3]
    F = (cat_const_afl95 * cat_const_fmax * cat_const_cosa95)[0,:]
    W = pinv(J.T) * matrix(R) * matrix(diag(F))
    return name, W

In [16]:
from numpy import savetxt

num_cats = 3

for cat_index in range(num_cats):
    name, W = compute_cat_generators(cat_index)
    savetxt("cat.{0}.wrenches.txt".format(name), W.T)
    savetxt("cat.{0}.forces.txt".format(name), W[0:3,:].T)
    savetxt("cat.{0}.torques.txt".format(name), W[3:,:].T)

Here's the list of force generators for the cat hindlimb.


In [17]:
W[0:3,:].T


Out[17]:
matrix([[ -2.44144597e+01,  -1.57391323e+01,  -8.51755859e-01],
        [ -1.76362162e+00,  -1.09434708e+00,  -2.00808351e-01],
        [ -1.12608949e+01,  -7.44499529e+00,   2.91354054e-01],
        [ -2.59086729e+01,   3.93167551e+01,  -6.77622942e+00],
        [ -1.25073288e+00,   1.07304395e+00,   5.56023623e-01],
        [  5.24475736e-01,  -2.21835121e-01,   9.10821010e-01],
        [  2.19437869e+00,  -2.96502640e+00,  -3.89770846e-01],
        [ -1.31759252e-01,  -2.02506505e-01,   4.96595569e-01],
        [  1.64008905e+00,  -9.46539758e-01,   8.15345592e+00],
        [  1.42642854e-01,   1.16865462e-02,   3.36018769e-01],
        [ -4.40270999e+00,   6.89681815e+00,  -5.15604778e-01],
        [  8.76852537e+00,   3.23926608e+00,   5.35209885e-01],
        [  7.89934504e+00,   1.72344139e+00,  -2.03038106e+00],
        [ -4.53775661e-02,  -2.69847066e-01,  -6.74883219e-01],
        [ -8.26081623e-01,  -5.24387291e-01,  -5.19867076e-02],
        [ -1.70624218e-01,  -3.83630242e-02,  -5.24256556e-01],
        [  8.26939533e+00,   9.34346097e-01,   9.53473514e-02],
        [  4.17226442e+00,   2.96652197e+00,  -7.51142707e-01],
        [ -8.98606838e-02,  -1.64605778e-02,  -2.68383755e-01],
        [  1.46375228e-01,  -1.92752547e-01,   9.96999568e-01],
        [ -2.59818579e+00,  -1.72096997e+00,  -1.53072481e-01],
        [ -3.53353536e-01,  -1.21411311e+01,   1.12912012e+00],
        [  3.95114528e+00,   3.10692829e+00,  -1.88074914e+00],
        [ -1.68179590e+01,  -7.42670812e+00,   2.61098965e+00],
        [  1.26808887e+00,  -1.70416356e+00,  -2.06128582e-01],
        [ -1.64609946e+01,   2.15930819e+01,   4.75846817e+00],
        [ -1.31423849e+00,   1.86092751e+00,   4.09018472e-01],
        [  7.25453460e-01,   2.40309334e-02,   1.94219424e+00],
        [ -1.14870538e+00,  -4.47679801e+00,   2.27106826e-01],
        [ -3.87209406e+00,  -1.51137800e+01,   2.23495839e-01],
        [ -1.63985926e+00,  -6.36124470e+00,   1.01701157e+00]])