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 nameR
(7x31): moment arm matrixJ
(6x7): geometric Jacobian, mapping endpoing wrench (6x1) and joint torque (7x1)q
(7x1): model posture matched to experimental cat configuration at normal preferred stanceDoFs
{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.
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.
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]:
So it's a dictionary
In [3]:
mat_data_struct.keys()
Out[3]:
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]:
In [6]:
mat_cats_struct.shape
Out[6]:
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]:
In [8]:
mat_cats_struct[0,0].shape
Out[8]:
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]:
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)))
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.
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]: