#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 25 15:50:03 2021
@author: emi
"""
# ------------------------------------------------------------------------ #
# Copyright (c) 2019 Modenese L., Ceseracciu, E., Reggiani M., Lloyd, D.G. #
# #
# Licensed under the Apache License, Version 2.0 (the "License"); #
# you may not use this file except in compliance with the License. #
# You may obtain a copy of the License at #
# http://www.apache.org/licenses/LICENSE-2.0. #
# #
# Unless required by applicable law or agreed to in writing, software #
# distributed under the License is distributed on an "AS IS" BASIS, #
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or #
# implied. See the License for the specific language governing #
# permissions and limitations under the License. #
# #
# Author: Luca Modenese, January 2015 #
# email: l.modenese@imperial.ac.uk #
# ------------------------------------------------------------------------ #
#
# This function optimizes the muscle parameters as described in Modenese L,
# Ceseracciu E, Reggiani M, Lloyd DG (2015). Estimation of
# musculotendon parameters for scaled and subject specific musculoskeletal
# models using an optimization technique. Journal of Biomechanics (submitted)
# and prints the results to command window.
# Also it stores information about the optimization in the structure SimInfo
# Written by Emiliano Ravera emiliano.ravera@uner.edu.ar as part of the
# Python version of work by Luca Modenese in the parameterisation of muscle
# tendon properties.
#------ import packages ---------
import numpy as np
from pathlib import Path
from scipy import linalg, optimize
from sklearn.metrics import mean_squared_error
from time import time
import logging
# adding OpenSim to python script
import sys
sys.path.append('/opt/opensim-core/lib/python3.6/site-packages/')
import opensim
# adding tool function to python script
from Private_functions import sampleMuscleQuantities
#--------------------------------
def optimMuscleParams(osimModel_ref_filepath, osimModel_targ_filepath, N_eval, log_folder):
# results file identifier
res_file_id_exp = '_N' + str(N_eval)
# import models
osimModel_ref = opensim.Model(osimModel_ref_filepath)
osimModel_targ = opensim.Model(osimModel_targ_filepath)
# models details
name = Path(osimModel_targ_filepath).stem
ext = Path(osimModel_targ_filepath).suffix
# assigning new name to the model
osimModel_opt_name = name + '_opt' + res_file_id_exp + ext
osimModel_targ.setName(osimModel_opt_name)
# initializing log file
log_folder = Path(log_folder)
log_folder.mkdir(parents=True, exist_ok=True)
logging.basicConfig(filename = str(log_folder) + '/' + name + '_opt' + res_file_id_exp +'.log', filemode = 'w', format = '%(levelname)s:%(message)s', level = logging.INFO)
# get muscles
muscles = osimModel_ref.getMuscles()
muscles_scaled = osimModel_targ.getMuscles()
# initialize with recognizable values
LmOptLts_opt = -1000*np.ones((muscles.getSize(),2))
SimInfo = {}
for n_mus in range(0, muscles.getSize()):
tic = time()
# current muscle name (here so that it is possible to choose a single muscle when developing).
curr_mus_name = muscles.get(n_mus).getName()
print('processing mus ' + str(n_mus+1) + ': ' + curr_mus_name)
# import muscles
curr_mus = muscles.get(curr_mus_name)
curr_mus_scaled = muscles_scaled.get(curr_mus_name)
# extracting the muscle parameters from reference model
LmOptLts = [curr_mus.getOptimalFiberLength(), curr_mus.getTendonSlackLength()]
PenAngleOpt = curr_mus.getPennationAngleAtOptimalFiberLength()
Mus_ref = sampleMuscleQuantities(osimModel_ref,curr_mus,'all',N_eval)
# calculating minimum fiber length before having pennation 90 deg
# acos(0.1) = 1.47 red = 84 degrees, chosen as in OpenSim
limitPenAngle = np.arccos(0.1)
# this is the minimum length the fiber can be for geometrical reasons.
LfibNorm_min = np.sin(PenAngleOpt) / np.sin(limitPenAngle)
# LfibNorm as calculated above can be shorter than the minimum length
# at which the fiber can generate force (taken to be 0.5 Zajac 1989)
if LfibNorm_min < 0.5:
LfibNorm_min = 0.5
# muscle-tendon paramenters value
MTL_ref = [musc_param_iter[0] for musc_param_iter in Mus_ref]
LfibNorm_ref = [musc_param_iter[1] for musc_param_iter in Mus_ref]
LtenNorm_ref = [musc_param_iter[2]/LmOptLts[1] for musc_param_iter in Mus_ref]
penAngle_ref = [musc_param_iter[4] for musc_param_iter in Mus_ref]
# LfibNomrOnTen_ref = LfibNorm_ref.*cos(penAngle_ref)
LfibNomrOnTen_ref = [(musc_param_iter[1]*np.cos(musc_param_iter[4])) for musc_param_iter in Mus_ref]
# checking the muscle configuration that do not respect the condition.
okList = [pos for pos, value in enumerate(LfibNorm_ref) if value > LfibNorm_min]
# keeping only acceptable values
MTL_ref = np.array([MTL_ref[index] for index in okList])
LfibNorm_ref = np.array([LfibNorm_ref[index] for index in okList])
LtenNorm_ref = np.array([LtenNorm_ref[index] for index in okList])
penAngle_ref = np.array([penAngle_ref[index] for index in okList])
LfibNomrOnTen_ref = np.array([LfibNomrOnTen_ref[index] for index in okList])
# in the target only MTL is needed for all muscles
MTL_targ = sampleMuscleQuantities(osimModel_targ,curr_mus_scaled,'MTL',N_eval)
evalTotPoints = len(MTL_targ)
MTL_targ = np.array([MTL_targ[index] for index in okList])
evalOkPoints = len(MTL_targ)
# The problem to be solved is:
# [LmNorm*cos(penAngle) LtNorm]*[Lmopt Lts]' = MTL;
# written as Ax = b or their equivalent (A^T A) x = (A^T b)
A = np.array([LfibNomrOnTen_ref , LtenNorm_ref]).T
b = MTL_targ
# ===== LINSOL =======
# solving the problem to calculate the muscle param
x = linalg.solve(np.dot(A.T , A) , np.dot(A.T , b))
LmOptLts_opt[n_mus] = x
# checking the results
if np.min(x) <= 0:
# informing the user
line0 = ' '
line1 = 'Negative value estimated for muscle parameter of muscle ' + curr_mus_name + '\n'
line2 = ' Lm Opt Lts' + '\n'
line3 = 'Template model : ' + str(LmOptLts) + '\n'
line4 ='Optimized param : ' + str(LmOptLts_opt[n_mus]) + '\n'
# ===== IMPLEMENTING CORRECTIONS IF ESTIMATION IS NOT CORRECT =======
x = optimize.nnls(np.dot(A.T , A) , np.dot(A.T , b))
x = x[0]
LmOptLts_opt[n_mus] = x
line5 = 'Opt params (optimize.nnls): ' + str(LmOptLts_opt[n_mus])
logging.info(line0 + line1 + line2 + line3 + line4 + line5 + '\n')
# In our tests, if something goes wrong is generally tendon slack
# length becoming negative or zero because tendon length doesn't change
# throughout the range of motion, so lowering the rank of A.
if np.min(x) <= 0:
# analyzes of Lten behaviour
Lten_ref = [musc_param_iter[2] for musc_param_iter in Mus_ref]
Lten_ref = np.array([Lten_ref[index] for index in okList])
if (np.max(Lten_ref) - np.min(Lten_ref)) < 0.0001:
logging.warning(' Tendon length not changing throughout range of motion')
# calculating proportion of tendon and fiber
Lten_fraction = Lten_ref/MTL_ref
Lten_targ = Lten_fraction*MTL_targ
# first round: optimizing Lopt maintaing the proportion of
# tendon as in the reference model
A1 = np.array([LfibNomrOnTen_ref , LtenNorm_ref*0]).T
b1 = MTL_targ - Lten_targ
x1 = optimize.nnls(np.dot(A1.T , A1) , np.dot(A1.T , b1))
x[0] = x1[0][0]
# second round: using the optimized Lopt to recalculate Lts
A2 = np.array([LfibNomrOnTen_ref*0 , LtenNorm_ref]).T
b2 = MTL_targ - np.dot(A1,x1[0])
x2 = optimize.nnls(np.dot(A2.T , A2) , np.dot(A2.T , b2))
x[1] = x2[0][1]
LmOptLts_opt[n_mus] = x
# Here tests about/against optimizers were implemented
# calculating the error (mean squared errors)
fval = mean_squared_error(b, np.dot(A,x), squared=False)
# update muscles from scaled model
curr_mus_scaled.setOptimalFiberLength(LmOptLts_opt[n_mus][0])
curr_mus_scaled.setTendonSlackLength(LmOptLts_opt[n_mus][1])
# PRINT LOGS
toc = time() - tic
line0 = ' '
line1 = 'Calculated optimized muscle parameters for ' + curr_mus.getName() + ' in ' + str(toc) + ' seconds.' + '\n'
line2 = ' Lm Opt Lts' + '\n'
line3 = 'Template model : ' + str(LmOptLts) + '\n'
line4 = 'Optimized param : ' + str(LmOptLts_opt[n_mus]) + '\n'
line5 = 'Nr of eval points : ' + str(evalOkPoints) + '/' + str(evalTotPoints) + ' used' + '\n'
line6 = 'fval : ' + str(fval) + '\n'
line7 = 'var from template [%]: ' + str(100*(np.abs(LmOptLts - LmOptLts_opt[n_mus])) / LmOptLts) + '%' + '\n'
logging.info(line0 + line1 + line2 + line3 + line4 + line5 + line6 + line7 + '\n')
# SIMULATION INFO AND RESULTS
SimInfo[n_mus] = {}
SimInfo[n_mus]['colheader'] = curr_mus.getName()
SimInfo[n_mus]['LmOptLts_ref'] = LmOptLts
SimInfo[n_mus]['LmOptLts_opt'] = LmOptLts_opt[n_mus]
SimInfo[n_mus]['varPercLmOptLts'] = 100*(np.abs(LmOptLts - LmOptLts_opt[n_mus])) / LmOptLts
SimInfo[n_mus]['sampledEvalPoints'] = evalOkPoints
SimInfo[n_mus]['sampledEvalPoints'] = evalTotPoints
SimInfo[n_mus]['fval'] = fval
# assigning optimized model as output
osimModel_opt = osimModel_targ;
return osimModel_opt, SimInfo
Datasets
Models
