--- a
+++ b/Python_tool/MuscleParOptTool/optimMuscleParams.py
@@ -0,0 +1,252 @@
+#!/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
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+