Skip to content

Commit

Permalink
cleaning up
Browse files Browse the repository at this point in the history
  • Loading branch information
@antoinefalisse
antoinefalisse committed Sep 29, 2022
1 parent b73343d commit 00e428c
Show file tree
Hide file tree
Showing 13 changed files with 1,248 additions and 1,160 deletions.
23 changes: 23 additions & 0 deletions OpenSimPipeline/JointReaction/Setup_JointReaction.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
<?xml version="1.0" encoding="UTF-8" ?>
<OpenSimDocument Version="40000">
<JointReaction name="JointReactionAnalysis">
<!--Flag (true or false) specifying whether on. True by default.-->
<on>true</on>
<!--Start time.-->
<start_time>1</start_time>
<!--End time.-->
<end_time>Inf</end_time>
<!--Specifies how often to store results during a simulation. More specifically, the interval (a positive integer) specifies how many successful integration steps should be taken before results are recorded again.-->
<step_interval>1</step_interval>
<!--Flag (true or false) indicating whether the results are in degrees or not.-->
<in_degrees>true</in_degrees>
<!--The name of a file containing forces storage. If a file name is provided, the forces for all actuators will be applied according to values specified in the forces_file instead of being computed from the states. This option should be used to calculate joint reactions from static optimization results.-->
<forces_file />
<!--Names of the joints on which to perform the analysis. The key word 'All' indicates that the analysis should be performed for all joints.-->
<joint_names> all </joint_names>
<!--Choice of body ('parent' or 'child') for which the reaction loads are calculated. Child body is default. The array must either have one entry or the same number of entries as joints specified above. If the array has one entry only, that selection is applied to all chosen joints.-->
<apply_on_bodies> child</apply_on_bodies>
<!--Names of frames in which the calculated reactions are expressed, or the keyword 'child' or 'parent' to indicate the joint's 'child' or 'parent' Frame. ground is default. If a Frame named 'child' or 'parent' exists and the keyword 'child' or 'parent' is used, the analysis will use that Frame. The array must either have one entry or the same number of entries as joints specified above. If the array has one entry only, that selection is applied to all chosen joints.-->
<express_in_frame> child</express_in_frame>
</JointReaction>
</OpenSimDocument>
677 changes: 677 additions & 0 deletions OpenSimPipeline/JointReaction/computeJointLoading.py
View file

Large diffs are not rendered by default.

137 changes: 28 additions & 109 deletions UtilsDynamicSimulations/OpenSimAD/boundsOpenSimAD.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,21 @@
'''
---------------------------------------------------------------------------
OpenCap processing: boundsOpenSimAD.py
---------------------------------------------------------------------------
Copyright 2022 Stanford University and the Authors
Author(s): Antoine Falisse, Scott Uhlrich
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.
'''

import scipy.interpolate as interpolate
from scipy import signal
import pandas as pd
Expand All @@ -9,10 +27,7 @@ def __init__(self, Qs, joints, muscles):

self.Qs = Qs
self.joints = joints
# self.targetSpeed = targetSpeed
self.muscles = muscles
# self.armJoints = armJoints
# self.mtpJoints = mtpJoints

def splineQs(self):

Expand All @@ -35,11 +50,14 @@ def splineQs(self):
order = 4
w = fc / (fs / 2) # Normalize the frequency
b, a = signal.butter(order/2, w, 'low')
output = signal.filtfilt(b, a, self.Qdotdots_spline.loc[:, self.Qdotdots_spline.columns != 'time'], axis=0,
padtype='odd', padlen=3*(max(len(b),len(a))-1))
output = signal.filtfilt(
b, a,
self.Qdotdots_spline.loc[:, self.Qdotdots_spline.columns!='time'],
axis=0, padtype='odd', padlen=3*(max(len(b),len(a))-1))
output = pd.DataFrame(data=output, columns=self.joints)
self.Qdotdots_spline_filter = pd.concat([pd.DataFrame(data=self.Qdotdots_spline['time'], columns=['time']),
output], axis=1)
self.Qdotdots_spline_filter = pd.concat(
[pd.DataFrame(data=self.Qdotdots_spline['time'], columns=['time']),
output], axis=1)

def getBoundsPosition(self):
self.splineQs()
Expand Down Expand Up @@ -218,7 +236,7 @@ def getBoundsAccelerationFiltered(self):
return (upperBoundsAcceleration, lowerBoundsAcceleration,
scalingAcceleration)

def getBoundsActivation(self, lb_activation=0.05):
def getBoundsActivation(self, lb_activation=0.01):
lb = [lb_activation]
lb_vec = lb * len(self.muscles)
ub = [1]
Expand Down Expand Up @@ -274,9 +292,8 @@ def getBoundsForce(self):

return upperBoundsForce, lowerBoundsForce, scalingForce

def getBoundsActivationDerivative(self):
activationTimeConstant = 0.015
deactivationTimeConstant = 0.06
def getBoundsActivationDerivative(self, activationTimeConstant=0.015,
deactivationTimeConstant=0.06):
lb = [-1 / deactivationTimeConstant]
lb_vec = lb * len(self.muscles)
ub = [1 / activationTimeConstant]
Expand Down Expand Up @@ -442,14 +459,6 @@ def getBoundsLumbarActivation(self, lumbarJoints):

def getBoundsReserveActuators(self, joint, value):

# max_value = {}
# max_value['hip_flexion_l'] = 20
# max_value['hip_flexion_r'] = 20
# max_value['hip_adduction_l'] = 20
# max_value['hip_adduction_r'] = 20
# max_value['hip_rotation_l'] = 20
# max_value['hip_rotation_r'] = 20

lb = [-1]
lb_vec = lb
ub = [1]
Expand All @@ -470,93 +479,3 @@ def getBoundsOffset(self, scaling):
lowerBoundsOffset = pd.DataFrame([-0.5 / scaling], columns=['offset_y'])

return upperBoundsOffset, lowerBoundsOffset

def getBoundsContactParameters(self, NContactSpheres,
parameter_to_optimize):


# lbRadius = 0.022
# ubRadius = 0.042

lbRadius = 0.01
ubRadius = 0.04

if NContactSpheres == 6:

lbLocation_s1 = np.array([-0.01, -0.015]) - 0.005
lbLocation_s2 = np.array([0.14, -0.035]) - 0.005
lbLocation_s3 = np.array([0.12, 0]) - 0.005
lbLocation_s4 = np.array([0.04, -0.03]) - 0.005
lbLocation_s5 = np.array([0.0, -0.03]) - 0.005
lbLocation_s6 = np.array([0.0, -0.02]) - 0.005
# lbLocation_s5 = np.array([0.0, -0.03])
# lbLocation_s6 = np.array([0.05, -0.02])

ubLocation_s1 = np.array([0.05, 0.01]) + 0.005
ubLocation_s2 = np.array([0.185, 0.01]) + 0.005
ubLocation_s3 = np.array([0.165, 0.05]) + 0.005
ubLocation_s4 = np.array([0.12, 0.03]) + 0.005
ubLocation_s5 = np.array([0.07, 0.03]) + 0.005
ubLocation_s6 = np.array([0.052, 0.05]) + 0.005
# ubLocation_s5 = np.array([0.04, 0.03])
# ubLocation_s6 = np.array([0.07, 0.02])

# lbLocation_s1 = np.array([0.00215773306688716053, -0.00434269152195360195]) - 0.03
# lbLocation_s2 = np.array([0.16841223157345971972, -0.03258850869005603529]) - 0.03
# lbLocation_s3 = np.array([0.15095065283989317351, 0.05860493716970469752]) - 0.03
# lbLocation_s4 = np.array([0.07517351958454182581, 0.02992219727974926649]) - 0.03
# lbLocation_s5 = np.array([0.06809743951165971032, -0.02129214951175864221]) - 0.03
# lbLocation_s6 = np.array([0.05107307963374478621, 0.07020500618327656095]) - 0.03

# ubLocation_s1 = np.array([0.00215773306688716053, -0.00434269152195360195]) + 0.03
# ubLocation_s2 = np.array([0.16841223157345971972, -0.03258850869005603529]) + 0.03
# ubLocation_s3 = np.array([0.15095065283989317351, 0.05860493716970469752]) + 0.03
# ubLocation_s4 = np.array([0.07517351958454182581, 0.02992219727974926649]) + 0.03
# ubLocation_s5 = np.array([0.06809743951165971032, -0.02129214951175864221]) + 0.03
# ubLocation_s6 = np.array([0.05107307963374478621, 0.07020500618327656095]) + 0.03

if parameter_to_optimize == 'option1':

lbContactParameters_unsc = np.concatenate((lbLocation_s1, lbLocation_s2, lbLocation_s3, lbLocation_s4, lbLocation_s5, lbLocation_s6, [lbRadius], [lbRadius], [lbRadius], [lbRadius], [lbRadius], [lbRadius]))
ubContactParameters_unsc = np.concatenate((ubLocation_s1, ubLocation_s2, ubLocation_s3, ubLocation_s4, ubLocation_s5, ubLocation_s6, [ubRadius], [ubRadius], [ubRadius], [ubRadius], [ubRadius], [ubRadius]))

scalingContactParameters_v = 1 / (ubContactParameters_unsc - lbContactParameters_unsc)
scalingContactParameters_r = 0.5 - ubContactParameters_unsc / (ubContactParameters_unsc - lbContactParameters_unsc)

lowerBoundsContactParameters = -0.5 * np.ones((1, len(lbContactParameters_unsc)))
upperBoundsContactParameters = 0.5 * np.ones((1, len(ubContactParameters_unsc)))

return upperBoundsContactParameters, lowerBoundsContactParameters, scalingContactParameters_v, scalingContactParameters_r

def getBoundsGR(self, GR, headers):
upperBoundsGR = pd.DataFrame()
lowerBoundsGR = pd.DataFrame()
scalingGR = pd.DataFrame()
for count, header in enumerate(headers):
if (header[0] == 'R' or header[0] == 'L'):
ub = max(max(GR['R' + header[1:]]),
max(GR['L' + header[1:]]))
lb = min(min(GR['R' + header[1:]]),
min(GR['L' + header[1:]]))
else:
raise ValueError("Problem bounds GR")
r = abs(ub - lb)
ub = ub + r
lb = lb - r
upperBoundsGR.insert(count, header, [ub])
lowerBoundsGR.insert(count, header, [lb])

# Scaling
s = pd.concat([abs(upperBoundsGR[header]),
abs(lowerBoundsGR[header])]).max(level=0)
scalingGR.insert(count, header, s)
lowerBoundsGR[header] /= scalingGR[header]
upperBoundsGR[header] /= scalingGR[header]

return upperBoundsGR, lowerBoundsGR, scalingGR

# def getBoundsFinalTime(self):
# upperBoundsFinalTime = pd.DataFrame([1], columns=['time'])
# lowerBoundsFinalTime = pd.DataFrame([0.1], columns=['time'])

# return upperBoundsFinalTime, lowerBoundsFinalTime
Loading

0 comments on commit 00e428c

Please sign in to comment.