analysistool, jointreaction analysis, muscle foces file, matlab

[Martin Mössner]

Up to now I succeeded to wite Matlab code to call the AnalysisTool performing a JointReaction analysis.
In the configuration file JointReaction.xml I have to set the forces_file tag to supplie the muscle forces,
which were priorly calculated ba a call of the AnalysisTool performing a StaticOptimization.
For a single call this works fine

Now I want to make a loop:

for i=1:3
prefix = sprintf('mu%1d',i);
SOfname = [prefix,'_StaticOptimization_force.sto']
AT = AnalyzeTool('JointReaction.xml')
AT.setName=prefix;
AT.setForcesFileName(SOfname);
AT.run()
end

but Matlab claims that setForcesFileName is not known to AT.
Probably JointReaction is a subobject/derived object of Analyze tool.
But how to pass the name of the muscle forces to this subobjece?

[Dimitar Stanev]

The utilization of these tools other than the .xml setup file is a bit tricky. I have done something similar with Python for OpenSim 3.3. However, you will have to adapt it for OpenSim 4.0 since the API has changed a bit. In any case, I hope the following snippet of code can provide you some directions:

Code: Select all

def perform_jra(model_file, ik_file, grf_file, grf_xml, reserve_actuators,
                muscle_forces_file, results_dir, prefix='',
                joint_names=['All'],
                apply_on_bodies=['parent'],
                express_in_frame=['ground']):
    """Performs joint reactions analysis using OpenSim.

    Parameters
    ----------
    model_file: str
        OpenSim model (.osim)
    ik_file: str
        kinematics calculated from Inverse Kinematics
    grf_file: str
        the ground reaction forces
    grf_xml: str
        xml description containing how to apply the GRF forces
    reserve_actuators: str
        path to the reserve actuator .xml file
    muscle_forces_file: str
        path to the file containing the muscle forces from SO
    results_dir: str
        directory to store the results
    prefix: str
        prefix of the resultant joint reaction loads
    joint_names: list
        joint names of interest
    apply_on_bodies: list
        apply on child or parent
    express_in_frame: list
        frame to express results
    """
    assert(len(joint_names) == len(apply_on_bodies) == len(express_in_frame))

    # model
    model = opensim.Model(model_file)

    # prepare external forces xml file
    name = os.path.basename(grf_file)[:-8]
    external_loads = opensim.ExternalLoads(model, grf_xml)
    external_loads.setExternalLoadsModelKinematicsFileName(ik_file)
    external_loads.setDataFileName(grf_file)
    external_loads.setLowpassCutoffFrequencyForLoadKinematics(6)
    external_loads.printToXML(results_dir + name + '.xml')

    # TODO this may not be needed
    # add reserve actuators (must not be appended when performing JRA)
    # force_set = opensim.ForceSet(model, reserve_actuators)
    # force_set.setMemoryOwner(False)  # model will be the owner
    # for i in range(0, force_set.getSize()):
    #     model.updForceSet().append(force_set.get(i))
    #     # model.addForce(force_set.get(i))

    # construct joint reaction analysis
    motion = opensim.Storage(ik_file)
    joint_reaction = opensim.JointReaction(model)
    joint_reaction.setName('JointReaction')
    joint_reaction.setStartTime(motion.getFirstTime())
    joint_reaction.setEndTime(motion.getLastTime())
    joint_reaction.setForcesFileName(muscle_forces_file)
    joint_names_arr = opensim.ArrayStr()
    apply_on_bodies_arr = opensim.ArrayStr()
    express_in_frame_arr = opensim.ArrayStr()
    for j, b, f in zip(joint_names, apply_on_bodies, express_in_frame):
        joint_names_arr.append(j)
        apply_on_bodies_arr.append(b)
        express_in_frame_arr.append(f)

    joint_reaction.setJointNames(joint_names_arr)
    joint_reaction.setOnBody(apply_on_bodies_arr)
    joint_reaction.setInFrame(express_in_frame_arr)
    # model.addAnalysis(joint_reaction)
    model.updAnalysisSet().adoptAndAppend(joint_reaction)
    model.initSystem()

    # analysis
    analysis = opensim.AnalyzeTool(model)
    analysis.setName(prefix + name)
    analysis.setModel(model)
    analysis.setModelFilename(model_file)
    analysis.setInitialTime(motion.getFirstTime())
    analysis.setFinalTime(motion.getLastTime())
    analysis.setLowpassCutoffFrequency(6)
    analysis.setCoordinatesFileName(ik_file)
    analysis.setExternalLoadsFileName(results_dir + name + '.xml')
    analysis.setLoadModelAndInput(True)
    analysis.setResultsDir(results_dir)
    analysis.run()
    jra_file = results_dir + name + '_JointReaction_ReactionLoads.sto'
    return jra_file

Please also be aware of this issue: https://github.com/opensim-org/opensim-core/issues/2076

[Martin Mössner]

The python script above gives some hints, but I was not able to succeed

Following code fragment produces a proper JointReaction configuration:

JR = JointReaction()
JR.setName('JointReaction')
JR.setStartTime(0.0)
JR.setEndTime(1.0)
JR.setForcesFileName([prefix,'_',muscle_forces_from_static_optimization])
JR.print(fullfile(cwd,[prefix,'_jr.xml']))

which can be checked by looking at the printed xml file.
But if I try to change the configuration template, code is running but does not work

AT = AnalyzeTool('cfg_JointReaction_template.xml')
AT.setName(prefix)
% both lines do not change the configuration!
AT.getModel().updAnalysisSet().adoptAndAppend(JR)
% AT.getModel().addAnalysis(JR)
AT.print(fullfile(cwd,[prefix,'_cfg_JointReaction.xml']));
AT.run()

[Martin Mössner]

Why isn't there a function that allows to set the desired analysis in the AnalyzeTool object?
Just the same as in the GUI ... this is very disturbing, curious, time consuming, hard to understand ...

[Aaron Fox]

Hi Martin,

Not sure if this will work for you given I use Matlab instead of Python, but I think recently I used a slightly different command to 'update' the analysis set of a tool (it was CMC not the Analyze tool but I figure it still might hold true). I used the 'getAnalysisSet()' command and tacked on a 'cloneAndAppend(analysis)' command to this. It seemed to add the appropriate analysis to where I wanted it to.

Aaron

[James Cowburn]

Hi Martin,

Following on from what Aaron said, I thought I would provide some example code to show what Aaron was saying. Again it's using Matlab for OpenSim 3.3, but hopefully it'll help point you in the right direction.

Code: Select all

%-----Create JRA-----
JRA = JointReaction();
%Assign JRA Setup parameters
JRA.setName('JointReaction');
JRA.setStartTime(startTime);
JRA.setEndTime(endTime);
JRA.setForcesFileName(forcesFile);
%Set joints for JRA
joints = ArrayStr();
joints.append('walker_knee_r');
JRA.setJointNames(joints);
%Assign bodies JRF are applied to
bodies = ArrayStr();
bodies.append('child');
JRA.setOnBody(bodies);
%Assign bodies to express JRF in
expIn = ArrayStr();
expIn.append('child');
JRA.setInFrame(expIn);

%-----Create Analyze Tool-----
bkTool = AnalyzeTool();
%Append JRA
AS = bkTool.getAnalysisSet();
AS.adoptAndAppend(JRA);
%Assign AnalysisTool setup parameters

James

[Martin Mössner]

Thanks @Aaron, thanks @James
your comments lead to a working version ...
I already figured out that getAnalysisSet() gives a list of analyses

Here is my code. There are some questions, improvements I would like to add.

name0 = 'JointReaction'
cfg1 = ['cfg_',name0,'.xml']
a = AnalyzeTool(cfg1);
a.setName(prefix);

jr = JointReaction()
jr.setName('JointReaction')
jr.setStartTime(0.0)
jr.setEndTime(10.0)
muscleforces
jr.setForcesFileName(muscleforces)
jr.print(fullfile(cwd,[prefix,'_jr.xml']))
joints = ArrayStr();
joints.append('all');
jr.setJointNames(joints);
bodies = ArrayStr();
bodies.append('child');
jr.setOnBody(bodies);
express = ArrayStr();
express.append('ground');
jr.setInFrame(express);

as = a.getAnalysisSet()
% necessary because jointreaction.xml contains a sample jointreaction analysis
as.remove(0)
as.adoptAndAppend(jr)

cfg2 = [prefix,'_',cfg1]
a.print(fullfile(cwd,cfg2));

a = AnalyzeTool(cfg2);
a.run()
fname = [prefix,'_',name0,'_ReactionLoads.sto']

My JointReaction.xml template contains a working template for JointReaction (1st iteration).
In particular I would like to just change the muscleforces file in each iteration.

Therefore I would prefer to copy the jr object from the as list and move it back (replace)

I had to write the new configuration to a xml-file and reconstruct the AnalyzeTool by reopening.
Is this a freature? Why?

[Aaron Fox]

Hi Martin,

So from what I understand you want to keep the original joint reactions analysis ('jr') and just modify the muscle forces file each time, but this 'jr' object gets taken up by the analysis set?

I would potentially take a different approach than having a 'template' xml file that is loaded in and then edited, but rather create a blank AnalyzeTool and set the parameters and create your 'baseline' joint reactions analysis (perhaps even without the muscle forces file to begin with). Then you could possibly generate a loop that on each iteration creates a copy of these two objects, adds the appropriate muscle forces file and runs the analyze tool with these paramters; and then on each fresh loop iteration go back to starting with your 'baseline' objects. Make sense?

Aaron

[Martin Mössner]

Aaron I understand your approach ... but

I start with a working configuration file that works with the GUI in a single case.
In the next step I write Matlab code that does the same.
And in the third step I write a loop to do parameter studies.

In my case the simplest way would be just to change the muscle-forces file in each loop ...
and for this it would be logical just to open the AnalyzeTool in Matlab with the configuration
file of the GUI and then change the muscle forces file in each iteration, just as it is possible, for instance,
with the motion data.

Of course I could bild up the AnalyzeTool with the JointReaction Analysis from scratch, but then I have to
set all options by Matlab code. This is much more complicated (for a beginner), since there do not exist
working examples in Matlab and for OpenSim 4.0. One allways has to upgrade from 3.3 or rewrite from
python or c++ samples

But anyway, I now have a working solution , thanks to all

[Aaron Fox]

Fair enough Martin - I figured the XML load time for each iteration might increase computational time, but if you've got a working solution roll with it!

Aaron