Python/Matlab examples for using the tools

[mitkof6]

I find it a little bit difficult to setup the AnalyzeTool without the .xml file (through the API). It would be great if examples are provided (especially for Python and Matlab users) that demonstrate the simulation pipelines (e.g. gait2354, SO, IK, ID, JRA, etc.). Something like this: IK.

For example, it isn't very intuitive how to append reserve actuators and how to add different analysis. The doxygen is a little bit vague and does not explain the meaning of some functions and how to use them.

[chrisdembia]

Such examples would be great. @jimmyDunne worked on a utility script for appending reserve actuators.

#2060 is related.

[mitkof6]

@chrisdembia it seems that a lot of documentation is located in the .cpp files, it would be nice if these can be accessed in the doxygen.

[chrisdembia]

Yes we've discussed that a lot. Over time a lot of the documentation in the cpp files has become less than correct, and we think that incorrect documentation is worse than no documentation. We definitely want to move the comments from the cpp files to the headers, but they documentation should be corrected in the process.

[mitkof6]

Continuing our discussion, here are some comments for using the AnalyzeTool. This is how I managed to initiate the SO (v3.3) after a lot of effort:

string tempExternalLoadsXML = parameters.resultsDir +
    parameters.subjectName + "_external_loads.xml";
// prepare external loads
ExternalLoads externalLoads(parameters.model,
                            parameters.groundReactionXMLTemplate);
// update template
externalLoads.setExternalLoadsModelKinematicsFileName(
    parameters.inverseKinematicsMotion);
externalLoads.setDataFileName(parameters.groundReactionForces);
externalLoads.setLowpassCutoffFrequencyForLoadKinematics(6);
externalLoads.print(tempExternalLoadsXML);

// add reserve actuators
ForceSet forceSet(parameters.model, parameters.reserveActuators);
for (int i = 0; i < forceSet.getSize(); i++) {
    parameters.model.updForceSet().append(forceSet.get(i));
}

// add static optimization
Storage motion(parameters.inverseKinematicsMotion);
auto staticOptimization = new StaticOptimization();
staticOptimization->setStartTime(motion.getFirstTime());
staticOptimization->setEndTime(motion.getLastTime());
staticOptimization->setUseModelForceSet(true);
staticOptimization->setActivationExponent(2);
staticOptimization->setUseMusclePhysiology(true);
staticOptimization->setConvergenceCriterion(0.0001);
staticOptimization->setMaxIterations(100);
parameters.model.addAnalysis(staticOptimization);

AnalyzeTool tool;
tool.setName(parameters.subjectName);
tool.setModel(parameters.model);
tool.setInitialTime(motion.getFirstTime());
tool.setFinalTime(motion.getLastTime());
tool.setCoordinatesFileName(parameters.inverseKinematicsMotion);
// let createExternalLoads to handle the creation of the external loads
tool.setExternalLoadsFileName(tempExternalLoadsXML);
// _loadModelAndInput be set so that .mot file is loaded correctly
tool.setLoadModelAndInput(true);
tool.setLowpassCutoffFrequency(6);
tool.setResultsDir(parameters.resultsDir);
tool.run();

From a user perspective the tool.setExternalLoads does not work because in the run method of the analyze tool the external forces are created with createExternalLoads

https://github.com/opensim-org/opensim-core/blob/master/OpenSim/Tools/AnalyzeTool.cpp#L531

so the user can't pass the ExternalLoads object but is restricted to set the path to the .xml file. Another problem is the specification of the reserve actuators. The additional force set will not be added to the model even if it is specified by setForceSetFiles, because this is performed in the constructor of the tool rather than before the run method.

https://github.com/opensim-org/opensim-core/blob/master/OpenSim/Tools/AnalyzeTool.cpp#L91

Another problem is the motion vs state input. If the coordinates are specified from a .mot file (e.g. IK) then the _loadModelAndInput must be set to true, which is not well documented. This should be handled automatically when the user sets the corresponding input.

https://github.com/opensim-org/opensim-core/blob/master/OpenSim/Tools/AnalyzeTool.cpp#L540

The tools are an important part of the API, especially with the python and java wrappings. I suggest that the API is refined based on use cases so as to permit easy initialization and invocation, without restricting to the xml approach alone.

[chrisdembia]

Thanks @mitkof6; those are all good points. The documentation for _loadModelAndInput is easy to address.

[mitkof6]

Hi,

Two more comments.

1. The order of invocation of setLowpassCutoffFrequency matters in the AnalyzeTool

2. In OpenSim v3.3/3.2 I can call my python code to execute the JointReaction analysis. In the new version of from github the tool fails with an exception

Windows: access violation - no rtti data

Archlinux: Python segmentation fault (core dumped)

def perform_jra(model_file, ik_file, grf_file, grf_xml, reserve_actuators,
                muscle_forces_file, results_dir):
# 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')

# add reserve actuators
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)
model.addAnalysis(joint_reaction)

# analysis
analysis = opensim.AnalyzeTool(model)
analysis.setName(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()

When I comment joint_reaction.setForcesFileName(muscle_forces_file) the tool runs normally without exception.

Furthermore, force_set.setMemoryOwner(False) is necessary because the model takes ownership and when the execution goes out of scope we have segmentation fault in python.

PS: you can use the file from the gait model to try and reproduce this error.

[romainmartinez]

I get the same segmentation error as @mitkof6 when I run the Joint Reaction analysis with python.
The error also goes away when I comment:

joint_reaction.setForcesFileName('StaticOptimization_force.sto')

I would really be interested to know if PR #2191 could solve our problem.

[mitkof6]

@romainmartinez yes it seems so. You can try to verify this if you have time so we can be 100% sure.

[mitkof6]

When I perform batch joint reaction analysis the RAM is constantly increasing so I have to close python and re-run the batch from where it stopped. Is this a python side effect or some memory leak? This is my code:

def perform_jra(model_file, ik_file, grf_file, grf_xml, reserve_actuators,
                         muscle_forces_file, results_dir, prefix=''):
  # 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')

  # 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)
  model.addAnalysis(joint_reaction)

  # 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()

and the batch loop

print('Performing joint reaction analysis batch ...')
for i, force_file in enumerate(tqdm(force_files)):
    if i < 389:  # due to memory leaks in OpenSim this must be restarted
        continue
    perform_jra(model_file, ik_file, grf_file, grf_xml_file,
                reserve_actuators_file, feasible_set_dir + force_file,
                jra_results_dir, prefix=str(i) + '_')

[gattia]

@mitkof6 did you find a solution to this? Im finding the same problem.

I'm even running my batch processing using python subprocess so that I am calling a script that does my JR analysis (and all steps leading up to it) and once it's done that entire script should be closed, running the next batch in a separate python subprocess. I would think that ending the script should release that memory, but it doesnt seem to be the case.

A related quirk I've noticed is that killing a python script that is running an opensim analysis doesnt kill the opensim analysis - it keeps running in the background (shown as a python process), even though the python script is closed on the command line. In System Monitor, there is still a python process that is running and taking up memory. For the killed python script, it seems to take up memory and CPU resources until the opensim analysis is done (I cant know this for sure, but the time it takes to release makes this seem plausible), and then it gets released.

Im wondering if part of the batch processing problem is that there are still some latent C++ processes being run by opensim that if given some more time before starting the next loop might release that memory. But when not given the time, they just keep piling on. I say this because of the fact that killing a python script running an opensim analysis eventually does release the memory, it just doesnt do so until everything from the underlying C++ code that had been initiated is done.

Sorry, that's a lot of info/ideas. Im not a computer scientist, so these ideas could be out to left field. Hopefully this is helpful in some way.

sudo code example of the python script running subprocesses:

#opensim_batch.py is a file that does all my JR analysis and saves the result to file. 
#This file is something like the code at the following, but includes saving results to file. 
#https://github.com/gattia/randome_code/blob/master/opensim_scale_ik_id_so.py

import os
import sys
import gc
import pandas as pd

#Read in the list of participant IDs that is fed to the subprocess for processing using opensim script
file_with_list_data_to_process = 'path_to_file_for_batch_processing
list_files_to_process = pd.read_csv(file_with_list_data_to_process, index_col=False)

#loop over all participant IDs. 
for file_index in range(list_files_to_process.shape[0]):
    
    #run the python script as a subprocess giving participant ID as a keyword argument
    subprocess.call(['python',
                               'opensim_batch.py',
                               list_files_to_process[0,'Participant'])
    #Finally, run the garbage collector to make sure that memory gets freed. 
    gc.collect()

[gattia]

Just to clarify, I was talking about the memory issue mentioned for batch processing... sorry, I didn't say that explicitly :).

[mitkof6]

Hi @gattia,

Unfortunately, I haven't resolved this issue. It might be a memory leak in the joint reaction tool or anywhere else. I have to close the python process each time to free the memory.

[gattia]

@mitkof6 thanks for the response. I've done some more playing and realized that what I mentioned in my post (using subprocess to call a second python script that performs the analysis) actually does work. It's a bit of a hack, and doesnt solve the underlying problem, but if you are running into memory issues, it might be of interest.

In the event that you are interested in this solution, I've tried to briefly describe what I've done, and the advantages to it below.

I am performing batch processing for participants as well as trials within participants. In the code block I had in my previous comment, the opensim_batch.py batch processes trials for an individual participant and this slowly accumulates memory. However, because each participant is called independently using the subprocess module from another python script it releases memory after completing every participant. Using subprocess in this way is nice because depending on how much RAM and CPU compute you have, you can actually run multiple instances of this at the base script at the same time, effectively getting parallel computation of many analyses - I'm finding that opensim essentially maxes out one thread, so I could in theory run in parallel 2x the number of cores on my machine, or 2xnumber of cores-1 if you want to leave something available for other processes. In addition, you could also extend what I've mentioned so each individual trial, as well as participant, is processed in a subprocess, and this should completely eliminate the memory accumulation.

[aymanhab]

Could you guys provide a minimal C++ example to reproduce the issue so we can find out if this is Bindings (python, Matlab) specific or if it can be fixed for all use cases? Thanks

[aymanhab]

@mitkof6 ExternalLoads had undergone refactoring before 4.0 release so your comment above may or maynot apply now. Just a headsup.

[gattia]

Sorry, I have zero experience with C++, otherwise I'd be happy to help out.

In case it helps, my current version of opensim was built using this commit

The for loop I run is very similar to that described by @mitkof6, except I include IK in my loop as well. I essentially run from line 35 this script and onwards in a loop

[gattia]

I could try creating a self-contained minimal example in python using some of the opensim example data if that will help someone else re-create in C++.

[aymanhab]

Thanks @gattia yes, that would be awesome, a minimal example in any language would be a great first step 👍

[gattia]

@aymanhab I've uploaded a repository that includes some of the example gait data and tests accumulation of memory for batch processing. It records the memory accumulated after every loop for each of IK, ID, SO, and JR. I've included in the README examples of the resulsts printed from the script - these results are the average memory accumulated per loop (20 loops). For the example data I have, it seems that SO is the biggest culprit (40MB/loop), but all of them seem to accumulate at least a little memory.

[mitkof6]

@aymanhab I have prepared a project that reproduces this behavior as part of a demo for TGCS. Currently, only Linux is supported because I have problems porting gmp library for Windows (not tested on Mac). I am working on a Windows version. You will have to build the C++ source code and run CalculateFeasibleMuscleForces first. Then to reproduce this behavior run the

CalculateFeasibleJointReactionLoads 0 300

If you open a task manager or system monitor, you will observe that the RAM constantly increases with each iteration.

https://gitlab.com/mitkof6/feasible_muscle_force_analysis

[Se123R]

From a user perspective the tool.setExternalLoads does not work because in the run method of the analyze tool the external forces are created with createExternalLoads

https://github.com/opensim-org/opensim-core/blob/master/OpenSim/Tools/AnalyzeTool.cpp#L531

so the user can't pass the ExternalLoads object but is restricted to set the path to the .xml file. Another problem is the specification of the reserve actuators. The additional force set will not be added to the model even if it is specified by setForceSetFiles, because this is performed in the constructor of the tool rather than before the run method.

@mitkof6 I have a similar problem in running CMC in Matlab did you found a solution? Here is my code.

cmc = CMCTool('subject01_Setup_CMC.xml');
cmc.setDesiredKinematicsFileName('subject0_walk1_ik.mot'); 
ExLoad=ExternalLoads('subject01_walk_grf.xml',true);
ExLoad.setDataFileName('New_subject01_walk_grf.mot');
cmc.setExternalLoads(ExLoad);
cmc.setStartTime(0.9);
cmc.setFinalTime(1.1);
cmc.setResultsDir('CMC');
cmc.run();

@chrisdembia I can not change the "force data file" directory by using setDataFileName. It considers directory which is in 'subject01_walk_grf.xml'. So when I want to apply CMC on dataset with different external force it would not change. Do you have any idea? Thank you.

[mitkof6]

If the set external force does not work, then you should write the external loads setup file and use the cmc.setExternalLoadsFileName() instead. Here is an example how I avoid this problem:

def perform_so(model_file, ik_file, grf_file, grf_xml, reserve_actuators, results_dir):
    """Performs Static Optimization 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
    results_dir: str
        directory to store the results
    """
    # model
    model = opensim.Model(model_file)

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

    # add reserve actuators
    force_set = opensim.SetForces(reserve_actuators, True)
    force_set.setMemoryOwner(False)  # model will be the owner
    for i in range(0, force_set.getSize()):
        model.updForceSet().append(force_set.get(i))

    # construct static optimization
    motion = opensim.Storage(ik_file)
    static_optimization = opensim.StaticOptimization()
    static_optimization.setStartTime(motion.getFirstTime())
    static_optimization.setEndTime(motion.getLastTime())
    static_optimization.setUseModelForceSet(True)
    static_optimization.setUseMusclePhysiology(True)
    static_optimization.setActivationExponent(2)
    static_optimization.setConvergenceCriterion(0.0001)
    static_optimization.setMaxIterations(100)
    model.addAnalysis(static_optimization)

    # analysis
    analysis = opensim.AnalyzeTool(model)
    analysis.setName(name)
    analysis.setModel(model)
    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()
    so_force_file = results_dir + name + '_so_forces.sto'
    so_activations_file = results_dir + name + '_so_activations.sto'
    return (so_force_file, so_activations_file)

Please note that if you use OpenSim v3.3 you should use ForceSet instead of SetForce. Keep in mind that the setForceSetFile might be ignored if called in wrong order.

[chrisdembia]

Thank you so much for sharing that, @mitkof6 !

[Se123R]

@mitkof6 Thank you. It solved based on your idea, except I changed printToXML to print. Because MATLAB did not recognize printToXML. Here is the code.

ExLoad=ExternalLoads('subject01_walk_grf.xml',true);
ExLoad.setDataFileName('New_subject01_walk_grf.mot');
ExLoad.print('New_subject01_walk_grf.xml')

cmc = CMCTool('subject01_Setup_CMC.xml');
cmc.setDesiredKinematicsFileName('subject0_walk1_ik.mot'); 
cmc.setExternalLoadsFileName('New_subject01_walk_grf.xml');
cmc.setStartTime(0.9);
cmc.setFinalTime(1.1);
cmc.setResultsDir(results_dir);
cmc.run();