Skip to content

Commit 7c95008

Browse files
antoinefalisseantoinefalisse
committed
no assumption of symmetricity for polynomials
1 parent ff8f17b commit 7c95008

File tree

8 files changed

+5118
-5066
lines changed

8 files changed

+5118
-5066
lines changed

‎OpenSimModel/Hamner_modified/Model/Hamner_modified_scaled_polynomial_l.npy

16.8 KB
Binary file not shown.

‎OpenSimModel/Hamner_modified/Model/Hamner_modified_scaled_polynomial_r.npy

-3.01 KB
Binary file not shown.

‎OpenSimModel/Hamner_modified/Model/Hamner_modified_scaled_polynomial_trunk.npy

2.51 KB
Binary file not shown.

‎OpenSimModel/templates/MuscleAnalysis/dummy_motion.mot

Lines changed: 5000 additions & 5000 deletions
Large diffs are not rendered by default.

‎main.py

Lines changed: 107 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -145,7 +145,7 @@
145145
from muscleData import getBodyMass
146146
bodyMass = getBodyMass(pathModelFolder, modelName, loadBodyMass)
147147

148-
# %% Muscles.
148+
# %% Muscles (very messy).
149149
# This section is very specific to the OpenSim model being used.
150150
# The list 'muscles' includes all right leg muscles, as well both side
151151
# trunk muscles.
@@ -163,6 +163,9 @@
163163
rightSideMuscles = muscles[:-3]
164164
leftSideMuscles = [muscle[:-1] + 'l' for muscle in rightSideMuscles]
165165
bothSidesMuscles = leftSideMuscles + rightSideMuscles
166+
leftTrunkMuscles = muscles[-3:]
167+
rightTrunkMuscles = muscles[-6:-3]
168+
trunkMuscles = leftTrunkMuscles + rightTrunkMuscles
166169
nMuscles = len(bothSidesMuscles)
167170
nSideMuscles = len(rightSideMuscles)
168171

@@ -391,56 +394,73 @@
391394
from casadiFunctions import polynomialApproximation
392395
leftPolynomialJoints = [
393396
'hip_flexion_l', 'hip_adduction_l', 'hip_rotation_l', 'knee_angle_l',
394-
'ankle_angle_l', 'subtalar_angle_l', 'mtp_angle_l',
395-
'lumbar_extension', 'lumbar_bending', 'lumbar_rotation']
397+
'ankle_angle_l', 'subtalar_angle_l', 'mtp_angle_l']
396398
rightPolynomialJoints = [
397399
'hip_flexion_r', 'hip_adduction_r', 'hip_rotation_r', 'knee_angle_r',
398-
'ankle_angle_r', 'subtalar_angle_r', 'mtp_angle_r',
400+
'ankle_angle_r', 'subtalar_angle_r', 'mtp_angle_r']
401+
trunkPolynomialJoints = [
399402
'lumbar_extension', 'lumbar_bending', 'lumbar_rotation']
400403
if not withMTP:
401404
leftPolynomialJoints.remove(mtpJoints[0])
402405
rightPolynomialJoints.remove(mtpJoints[1])
403-
nPolynomialJoints = len(leftPolynomialJoints)
406+
# nPolynomialJoints = len(leftPolynomialJoints)
404407

405408
loadPolynomialData = False
406409
# Support loading/saving polynomial data such that they do not needed to
407410
# be recomputed every time.
408-
if os.path.exists(os.path.join(
409-
pathModelFolder,
410-
'{}_polynomial_{}.npy'.format(modelName, 'r'))):
411+
if (os.path.exists(
412+
os.path.join(pathModelFolder,
413+
'{}_polynomial_{}.npy'.format(modelName, 'r'))) and
414+
os.path.exists(
415+
os.path.join(pathModelFolder,
416+
'{}_polynomial_{}.npy'.format(modelName, 'l')))):
411417
loadPolynomialData = True
412418
from muscleData import getPolynomialData
413-
polynomialData = getPolynomialData(
419+
polynomialData = {}
420+
421+
polynomialData['r'] = getPolynomialData(
422+
loadPolynomialData, pathModelFolder, modelName,
423+
pathMotionFile4Polynomials, rightPolynomialJoints,
424+
rightSideMuscles[:-3], side='r')
425+
polynomialData['l'] = getPolynomialData(
426+
loadPolynomialData, pathModelFolder, modelName,
427+
pathMotionFile4Polynomials, leftPolynomialJoints,
428+
leftSideMuscles[:-3], side='l')
429+
polynomialData['trunk'] = getPolynomialData(
414430
loadPolynomialData, pathModelFolder, modelName,
415-
pathMotionFile4Polynomials, rightPolynomialJoints, muscles, side='r')
431+
pathMotionFile4Polynomials, trunkPolynomialJoints, trunkMuscles,
432+
side='trunk')
416433
if loadPolynomialData:
417-
polynomialData = polynomialData.item()
434+
polynomialData['r'] = polynomialData['r'].item()
435+
polynomialData['l'] = polynomialData['l'].item()
436+
polynomialData['trunk'] = polynomialData['trunk'].item()
418437

419438
# The function f_polynomial takes as inputs joint positions and velocities
420439
# from one side (trunk included), and returns muscle-tendon lengths,
421440
# velocities, and moments for the muscle of that side (trunk included).
422-
f_polynomial = polynomialApproximation(muscles, polynomialData,
423-
nPolynomialJoints)
441+
f_polynomial = {}
442+
f_polynomial['r'] = polynomialApproximation(
443+
rightSideMuscles[:-3], polynomialData['r'], len(rightPolynomialJoints))
444+
f_polynomial['l'] = polynomialApproximation(
445+
leftSideMuscles[:-3], polynomialData['l'], len(rightPolynomialJoints))
446+
f_polynomial['trunk'] = polynomialApproximation(
447+
trunkMuscles, polynomialData['trunk'], len(trunkPolynomialJoints))
424448
leftPolJointIdx = getJointIndices(joints, leftPolynomialJoints)
425449
rightPolJointIdx = getJointIndices(joints, rightPolynomialJoints)
450+
trunkPolJointIdx = getJointIndices(joints, trunkPolynomialJoints)
426451

427-
# The left and right polynomialMuscleIndices below are used to identify
428-
# the left and right muscles in the output of f_polynomial. Since
429-
# f_polynomial return data from side muscles (trunk included), we have the
430-
# side leg muscles + all trunk muscles as output. Here we make sure we
431-
# only include the side trunk muscles when identifying all side muscles.
432-
# This is pretty sketchy I know.
433-
rightPolMuscleIdx = [muscles.index(i) for i in rightSideMuscles]
434-
rightTrunkMuscles = ['ercspn_r', 'intobl_r', 'extobl_r']
435-
leftTrunkMuscles = ['ercspn_l', 'intobl_l', 'extobl_l']
436-
leftPolMuscleIdx = (
437-
[muscles.index(i) for i in rightSideMuscles
438-
if i not in rightTrunkMuscles] +
439-
[muscles.index(i) for i in leftTrunkMuscles])
452+
# Helper indices.
453+
sidePolMuscleIdx = [muscles.index(i) for i in rightSideMuscles[:-3]]
454+
leftTrunkPolMuscleIdx = [trunkMuscles.index(i) for i in leftTrunkMuscles]
455+
rightTrunkPolMuscleIdx = [trunkMuscles.index(i) for i in rightTrunkMuscles]
440456
from utilities import getMomentArmIndices
441457
momentArmIndices = getMomentArmIndices(
442458
rightSideMuscles, leftPolynomialJoints, rightPolynomialJoints,
443-
polynomialData)
459+
polynomialData['r'])
460+
for joint in trunkPolynomialJoints:
461+
momentArmIndices[joint] = getJointIndices(bothSidesMuscles,
462+
trunkMuscles)
463+
444464
trunkMomentArmPolynomialIndices = (
445465
[muscles.index(i) for i in leftTrunkMuscles] +
446466
[muscles.index(i) for i in rightTrunkMuscles])
@@ -453,10 +473,12 @@
453473
pathModelFolder, 'data4PolynomialFitting_{}.npy'.format(modelName))
454474
data4PolynomialFitting = np.load(path_data4PolynomialFitting,
455475
allow_pickle=True).item()
456-
momentArms = testPolynomials(
457-
data4PolynomialFitting, rightPolynomialJoints, muscles,
458-
f_polynomial, polynomialData, momentArmIndices,
459-
trunkMomentArmPolynomialIndices)
476+
_ = testPolynomials(data4PolynomialFitting, rightPolynomialJoints,
477+
rightSideMuscles[:-3], f_polynomial['r'],
478+
polynomialData['r'], momentArmIndices)
479+
_ = testPolynomials(data4PolynomialFitting, leftPolynomialJoints,
480+
leftSideMuscles[:-3], f_polynomial['l'],
481+
polynomialData['l'], momentArmIndices)
460482

461483
# %% External function.
462484
# The external function is written in C++ and compiled as a library, which
@@ -901,28 +923,36 @@
901923
# Left side.
902924
Qsinj_l = Qskj_nsc[leftPolJointIdx, j+1]
903925
Qdsinj_l = Qdskj_nsc[leftPolJointIdx, j+1]
904-
[lMTj_l, vMTj_l, dMj_l] = f_polynomial(Qsinj_l, Qdsinj_l)
926+
[lMTj_l, vMTj_l, dMj_l] = f_polynomial['l'](Qsinj_l, Qdsinj_l)
905927
# Right side.
906928
Qsinj_r = Qskj_nsc[rightPolJointIdx, j+1]
907929
Qdsinj_r = Qdskj_nsc[rightPolJointIdx, j+1]
908-
[lMTj_r, vMTj_r, dMj_r] = f_polynomial(Qsinj_r, Qdsinj_r)
909-
# Muscle-tendon lengths and velocities.
910-
lMTj_lr = ca.vertcat(lMTj_l[leftPolMuscleIdx],
911-
lMTj_r[rightPolMuscleIdx])
912-
vMTj_lr = ca.vertcat(vMTj_l[leftPolMuscleIdx],
913-
vMTj_r[rightPolMuscleIdx])
930+
[lMTj_r, vMTj_r, dMj_r] = f_polynomial['r'](Qsinj_r, Qdsinj_r)
931+
# Trunk.
932+
Qsinj_trunk = Qskj_nsc[trunkPolJointIdx, j+1]
933+
Qdsinj_trunk = Qdskj_nsc[trunkPolJointIdx, j+1]
934+
[lMTj_trunk, vMTj_trunk, dMj_trunk] = f_polynomial['trunk'](
935+
Qsinj_trunk, Qdsinj_trunk)
936+
# Muscle-tendon lengths and velocities.
937+
lMTj_lr = ca.vertcat(lMTj_l[sidePolMuscleIdx],
938+
lMTj_trunk[leftTrunkPolMuscleIdx],
939+
lMTj_r[sidePolMuscleIdx],
940+
lMTj_trunk[rightTrunkPolMuscleIdx])
941+
vMTj_lr = ca.vertcat(vMTj_l[sidePolMuscleIdx],
942+
vMTj_trunk[leftTrunkPolMuscleIdx],
943+
vMTj_r[sidePolMuscleIdx],
944+
vMTj_trunk[rightTrunkPolMuscleIdx],)
914945
# Moment arms.
915946
dMj = {}
916947
# Left side.
917-
for joint in leftPolynomialJoints:
948+
for c_j, joint in enumerate(leftPolynomialJoints):
918949
if ((joint != 'mtp_angle_l') and
919950
(joint != 'lumbar_extension') and
920951
(joint != 'lumbar_bending') and
921952
(joint != 'lumbar_rotation')):
922-
dMj[joint] = dMj_l[momentArmIndices[joint],
923-
leftPolynomialJoints.index(joint)]
953+
dMj[joint] = dMj_l[momentArmIndices[joint], c_j]
924954
# Right side.
925-
for joint in rightPolynomialJoints:
955+
for c_j, joint in enumerate(rightPolynomialJoints):
926956
if ((joint != 'mtp_angle_r') and
927957
(joint != 'lumbar_extension') and
928958
(joint != 'lumbar_bending') and
@@ -932,12 +962,10 @@
932962
# subtract by the number of side muscles.
933963
c_ma = [
934964
i - nSideMuscles for i in momentArmIndices[joint]]
935-
dMj[joint] = dMj_r[c_ma,
936-
rightPolynomialJoints.index(joint)]
965+
dMj[joint] = dMj_r[c_ma, c_j]
937966
# Trunk.
938-
for joint in trunkJoints:
939-
dMj[joint] = dMj_l[trunkMomentArmPolynomialIndices,
940-
leftPolynomialJoints.index(joint)]
967+
for c_j, joint in enumerate(trunkJoints):
968+
dMj[joint] = dMj_trunk[:, c_j]
941969

942970
###################################################################
943971
# Hill-equilibrium.
@@ -1462,16 +1490,25 @@
14621490
# Left leg.
14631491
Qsk_GC_l = Qs_GC_rad[leftPolJointIdx, k]
14641492
Qdsk_GC_l = Qds_GC_rad[leftPolJointIdx, k]
1465-
[lMTk_GC_l, vMTk_GC_l, _] = f_polynomial(Qsk_GC_l, Qdsk_GC_l)
1493+
[lMTk_GC_l, vMTk_GC_l, _] = f_polynomial['l'](Qsk_GC_l, Qdsk_GC_l)
14661494
# Right leg.
14671495
Qsk_GC_r = Qs_GC_rad[rightPolJointIdx, k]
14681496
Qdsk_GC_r = Qds_GC_rad[rightPolJointIdx, k]
1469-
[lMTk_GC_r, vMTk_GC_r, _] = f_polynomial(Qsk_GC_r, Qdsk_GC_r)
1497+
[lMTk_GC_r, vMTk_GC_r, _] = f_polynomial['r'](Qsk_GC_r, Qdsk_GC_r)
1498+
# Trunk.
1499+
Qsk_GC_trunk = Qs_GC_rad[trunkPolJointIdx, k]
1500+
Qdsk_GC_trunk = Qds_GC_rad[trunkPolJointIdx, k]
1501+
[lMTk_GC_trunk, vMTk_GC_trunk, _] = f_polynomial['trunk'](
1502+
Qsk_GC_trunk, Qdsk_GC_trunk)
14701503
# Both leg.
1471-
lMTk_GC_lr = ca.vertcat(lMTk_GC_l[leftPolMuscleIdx],
1472-
lMTk_GC_r[rightPolMuscleIdx])
1473-
vMTk_GC_lr = ca.vertcat(vMTk_GC_l[leftPolMuscleIdx],
1474-
vMTk_GC_r[rightPolMuscleIdx])
1504+
lMTk_GC_lr = ca.vertcat(lMTk_GC_l[sidePolMuscleIdx],
1505+
lMTk_GC_trunk[leftTrunkPolMuscleIdx],
1506+
lMTk_GC_r[sidePolMuscleIdx],
1507+
lMTk_GC_trunk[rightTrunkPolMuscleIdx])
1508+
vMTk_GC_lr = ca.vertcat(vMTk_GC_l[sidePolMuscleIdx],
1509+
vMTk_GC_trunk[leftTrunkPolMuscleIdx],
1510+
vMTk_GC_r[sidePolMuscleIdx],
1511+
vMTk_GC_trunk[rightTrunkPolMuscleIdx])
14751512

14761513
###################################################################
14771514
# Derive Hill-equilibrium.
@@ -1621,18 +1658,29 @@
16211658
# Left leg.
16221659
Qsinj_opt_l = Qskj_opt_nsc[leftPolJointIdx, j+1]
16231660
Qdsinj_opt_l = Qdskj_opt_nsc[leftPolJointIdx, j+1]
1624-
[lMTj_opt_l, vMTj_opt_l, _] = f_polynomial(Qsinj_opt_l,
1661+
[lMTj_opt_l, vMTj_opt_l, _] = f_polynomial['l'](Qsinj_opt_l,
16251662
Qdsinj_opt_l)
16261663
# Right leg.
16271664
Qsinj_opt_r = Qskj_opt_nsc[rightPolJointIdx, j+1]
16281665
Qdsinj_opt_r = Qdskj_opt_nsc[rightPolJointIdx, j+1]
1629-
[lMTj_opt_r, vMTj_opt_r, _] = f_polynomial(Qsinj_opt_r,
1666+
[lMTj_opt_r, vMTj_opt_r, _] = f_polynomial['r'](Qsinj_opt_r,
16301667
Qdsinj_opt_r)
1668+
# Trunk.
1669+
Qsinj_opt_trunk = Qskj_opt_nsc[trunkPolJointIdx, j+1]
1670+
Qdsinj_opt_trunk = Qdskj_opt_nsc[trunkPolJointIdx, j+1]
1671+
[lMTj_opt_trunk, vMTj_opt_trunk, _] = f_polynomial['trunk'](
1672+
Qsinj_opt_trunk, Qdsinj_opt_trunk)
16311673
# Both legs .
1632-
lMTj_opt_lr = ca.vertcat(lMTj_opt_l[leftPolMuscleIdx],
1633-
lMTj_opt_r[rightPolMuscleIdx])
1634-
vMTj_opt_lr = ca.vertcat(vMTj_opt_l[leftPolMuscleIdx],
1635-
vMTj_opt_r[rightPolMuscleIdx])
1674+
lMTj_opt_lr = ca.vertcat(
1675+
lMTj_opt_l[sidePolMuscleIdx],
1676+
lMTj_opt_trunk[leftTrunkPolMuscleIdx],
1677+
lMTj_opt_r[sidePolMuscleIdx],
1678+
lMTj_opt_trunk[rightTrunkPolMuscleIdx])
1679+
vMTj_opt_lr = ca.vertcat(
1680+
vMTj_opt_l[sidePolMuscleIdx],
1681+
vMTj_opt_trunk[leftTrunkPolMuscleIdx],
1682+
vMTj_opt_r[sidePolMuscleIdx],
1683+
vMTj_opt_trunk[rightTrunkPolMuscleIdx],)
16361684

16371685
###############################################################
16381686
# Derive Hill-equilibrium.

‎plotResults.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@
99
import matplotlib.pyplot as plt
1010

1111
# %% User inputs
12-
cases = ['0', '3']
12+
cases = ['0']
1313

1414
# %% Paths
1515
pathMain = os.getcwd()

‎polynomials.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -476,7 +476,7 @@ def testPolynomials(data4PolynomialFitting, joints, muscles,
476476
fig = plt.figure()
477477
fig.suptitle('Muscle-tendon lengths')
478478
for i in range(len(muscles)):
479-
muscle_obj = muscles[i][:-1] + 'r'
479+
muscle_obj = muscles[i]
480480
if polynomialData[muscle_obj]['dimension'] == 1:
481481
temp = polynomialData[muscle_obj]['spanning']==1
482482
y = (i for i,v in enumerate(temp) if v == True)
@@ -515,7 +515,7 @@ def testPolynomials(data4PolynomialFitting, joints, muscles,
515515
for j in range(NMomentarms):
516516
if joint[-1] == 'r' or joint[-1] == 'l':
517517
muscle_obj_r = (
518-
muscles[momentArmIndices[joint[:-1] + 'l'][j]][:-1] + 'r')
518+
muscles[momentArmIndices[joint[:-1] + 'l'][j]])
519519
muscle_obj = muscles[momentArmIndices[joint[:-1] + 'l'][j]]
520520
else:
521521
muscle_obj_r = muscles[trunkMomentArmPolynomialIndices[j]]

‎utilities.py

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -190,23 +190,27 @@ def getJointIndices(joints, selectedJoints):
190190
return jointIndices
191191

192192
# %% Get moment arm indices.
193-
def getMomentArmIndices(rightMuscles, leftPolynomialJoints,
193+
def getMomentArmIndices(muscles, leftPolynomialJoints,
194194
rightPolynomialJoints, polynomialData):
195195

196196
momentArmIndices = {}
197-
for count, muscle in enumerate(rightMuscles):
197+
for count, muscle in enumerate(muscles):
198+
if not muscle in polynomialData:
199+
continue
198200
spanning = polynomialData[muscle]['spanning']
199201
for i in range(len(spanning)):
200202
if (spanning[i] == 1):
201203
momentArmIndices.setdefault(
202204
leftPolynomialJoints[i], []).append(count)
203-
for count, muscle in enumerate(rightMuscles):
205+
for count, muscle in enumerate(muscles):
206+
if not muscle in polynomialData:
207+
continue
204208
spanning = polynomialData[muscle]['spanning']
205209
for i in range(len(spanning)):
206210
if (spanning[i] == 1):
207211
momentArmIndices.setdefault(
208212
rightPolynomialJoints[i], []).append(
209-
count + len(rightMuscles))
213+
count + len(muscles))
210214

211215
return momentArmIndices
212216

0 commit comments

Comments
 (0)