Sohn2013_hinlimb_models.matWe'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
CatsEach 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]
afl95scaling 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]: