pyomeca / bioptim

An optimization framework that links CasADi, Ipopt, ACADOS and biorbd for Optimal Control Problem

Pyomeca is a python library allowing you to carry out a complete biomechanical analysis; in a simple, logical and concise way.

Pyomeca documentation

See Pyomeca's documentation site.

Example

Pyomeca implements specialized functionalities commonly used in biomechanics. As an example, let's process the electromyographic data contained in this c3d file.

You can follow along without installing anything by using our binder server:

from pyomeca import Analogs

data_path = "../tests/data/markers_analogs.c3d"
muscles = [
    "Delt_ant",
    "Delt_med",
    "Delt_post",
    "Supra",
    "Infra",
    "Subscap",
]
emg = Analogs.from_c3d(data_path, suffix_delimiter=".", usecols=muscles)
emg.plot(x="time", col="channel", col_wrap=3)

emg_processed = (
    emg.meca.band_pass(order=2, cutoff=[10, 425])
    .meca.center()
    .meca.abs()
    .meca.low_pass(order=4, cutoff=5, freq=emg.rate)
    .meca.normalize()
)

emg_processed.plot(x="time", col="channel", col_wrap=3)

import matplotlib.pyplot as plt

fig, axes = plt.subplots(ncols=2, figsize=(10, 4))

emg_processed.mean("channel").plot(ax=axes[0])
axes[0].set_title("Mean EMG activation")

emg_processed.plot.hist(ax=axes[1], bins=50)
axes[1].set_title("EMG activation distribution")

See the documentation for more details and examples.

Features

• Signal processing routine commonly used in biomechanics such as filters, normalization, onset detection, outliers detection, derivatives, etc.
• Common matrix manipulation routines implemented such as getting Euler angles to/from a rototranslation matrix, creating a system of axes, setting a rotation or translation, transpose or inverse, etc.
• Easy reading and writing interface to common files in biomechanics (c3d, csv, xlsx,mat, trc, sto, mot)
• All of xarray's awesome features

The following illustration shows all of pyomeca's public API. An interactive version is available in the documentation.

Installation

Pyomeca itself is a pure Python package, but its dependencies are not. The easiest way to get everything installed is to use conda.

To install pyomeca with its recommended dependencies using the conda command line tool:

conda install -c conda-forge pyomeca

Now that you have installed pyomeca, you should be able to import it:

import pyomeca

Integration with other modules

Pyomeca is designed to work well with other libraries that we have developed:

• pyosim: interface between OpenSim and pyomeca to perform batch musculoskeletal analyses
• ezc3d: Easy to use C3D reader/writer in C++, Python and Matlab
• biorbd: C++ interface and add-ons to the Rigid Body Dynamics Library, with Python and Matlab binders.

Bug reports & questions

Pyomeca is Apache-licensed and the source code is available on GitHub. If any questions or issues come up as you use pyomeca, please get in touch via GitHub issues. We welcome any input, feedback, bug reports, and contributions.

Citing Pyomeca

If you use pyomeca in your academic work, please consider citing our paper as:

@article{Martinez2020,
  doi = {10.21105/joss.02431},
  url = {https://doi.org/10.21105/joss.02431},
  year = {2020},
  publisher = {The Open Journal},
  volume = {5},
  number = {53},
  pages = {2431},
  author = {Romain Martinez and Benjamin Michaud and Mickael Begon},
  title = {`pyomeca`: An Open-Source Framework for Biomechanical Analysis},
  journal = {Journal of Open Source Software}
}

Please consider citing the xarray project, which pyomeca is based on:

@article{jors-xarray,
title = {Xarray: N-D labeled arrays and datasets in Python},
author = {Joe Hamman and Stephan Hoyer},
year = {2017},
journal = {Journal of Open Research Software}
}

Contributors

• Romain Martinez
• Benjamin Michaud
• Mickael Begon
• Jenn Dowling-Medley
• Ariane Dang

Pyomeca is an open-source project created and supported by the S2M lab.

bioptim's Issues

custom objective/problem/penalty/etc. can't be save with a .bo

Add post-processing for lagrange multipliers

Add different interface for post analyses

show movement for multiple phase

Cannot use result.animate() for multiple phases.

Track data in parameters

Since we will move to SX, it becomes more interesting to store tracking data in fixed parameters instead of directly into the objective function.
One example of this is:
https://github.com/MMehrez/MPC-and-MHE-implementation-in-MATLAB-using-Casadi/blob/master/workshop_github/Codes_casadi_v3_4_5/MHE_code/MHE_Robot_PS_mul_shooting_v2.m

Make sure plot does not slow down ipopt

Save to a non existing folder

When saving if the folder doesn't exist it fails

ConditionPath as a function of time

Instead of giving all the frame, one could send arbitrary values as a function of time splined

Allow for easy accessible Contraints bounds

Instead of assuming min_bound = max_bounds = 0, this should be passable to the add_penalty function

RK4 steps from outside

ode_opt["number_of_finite_elements"] = 5 should be accessible from outside

nonlinear predictive model control

Add a nonlinear predictive model control

Write a proper readme

The readme is not filled properly but should

Single shooting

Implement a single shooting using the normal controls but sequential states

Plot integration when time is optimized

The plot doesn't work for the integration when the time is variant

Custom plot for user

Add an interface for the user to add new plots

Moving Horizon Estimation

Add an MHE possibility

Custom multiphase constraint

Give the opportunity to have multiphase constraints

Multi-trial

Make multi-trial available

New variable time constraint with time inequality

Add possibility of tracking one state at one node

https://github.com/pyomeca/BiorbdOptim/blob/e904785389195467739539f8bf78910092d7e234/biorbd_optim/penalty.py#L15
If I am not mistaken, currently, the penalty MINIMIZE_STATE only allows one to track one state during the full trajectory.
It would be more general to add the possibility of tracking one state at one particular node.
E.g. when relaxing a path constraint by passing it as a cost function.

Plotting Objective.X.MINIMIZE_STATE

This doesn't work anymore

Plot integration

The integration plot is now only init and final (without intermediate points). It should have more points

Proper integration of the Lagrange objective

The integration of the Lagrange objective function is rectangle. It would make sense to allow a RK integration (or trapezius).

phase time

phase_time = [t_phase[0], t_phase[0] + t_phase[1]]
It could be easier to put directly the duration of the phase instead of cumulative time

Plot objective values

Objective functions (with their respective weight) to compute the objective value should be plotted somehow

Plot tracking data

Data to track should be plotted on their respective figures

Modify min and max at once for Bounds

By adding an [] accessor, we should be able to modify Bounds like so:

x_bounds = BoundsOption([bound_min]*n_q, [bound_max]*n_q)
x_bounds[:n_q, 0] = new_equality_constraint

instead of

x_bounds = BoundsOption([bound_min]*n_q, [bound_max]*n_q)
x_bounds.min[:n_q, 0] = new_equality_constraint
x_bounds.max[:n_q, 0] = new_equality_constraint

Objective minimize_markers_displacement with ACADOS

Error message says there are free variables at the fist node (at least that's how I interpret it).
I guess the problem is that the Lagrange term must be on the entire interval (N points to sum), but finite difference give N-1 points.
Minimize_marker_velocity works fine, so it would be possible to 'fix' the problem by introducing the coordinates_system_idx option in this function Minimize_marker_velocity (to do so I still miss the RT velocity though :( )

Update setup.py

Error when running a program from command line
ModuleNotFoundError: No module named 'biorbd_optim.misc'

Lagrange objective function when optimizing the time

Since dt is an MX, does the ocp.J += casadi.dot(val, val) * weight * nlp["dt"] * nlp["dt"] function that adds the value to the objective scalar change of weighting in function the time of the phase? This may be a bug or a feature... To discuss

Unit test for all constraints, objective function, penalty, get_data

All the options of get data are not unit tested. We should fix that!

Undesired behavior of _check_and_fill_tracking_data_size

https://github.com/pyomeca/BiorbdOptim/blob/e904785389195467739539f8bf78910092d7e234/biorbd_optim/penalty.py#L493

This test crashes when data_to_track is a np.ndarray of the right shape but filled with zeros.
The first part of the logical test is False, so it evaluates the second part which is not defined with a np.ndarray:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Although I understand that only a 0 was expected when tracking a null vector, I believe this behavior to be undesired. E.g. in my case, I am creating a lot of different penalties in a for loop, and some times, the data_to_track can be a null vector.

3D visualisation during optimization

New graph?
BiordOptim?

Allow to instantiate OCP without X_bounds nor U_bounds

Piecewise linear control

Add the capability to control using piecewise controls

Plot parameters

Plot the parameters evolution over iterations

Docstring

Docstrings must deeply be reviewed

node initialization

node initialization

Add flexibility option for constraints [ready do be close]

It would be something like this :

constraints.add(
Constraint.ALIGN_MARKERS,
tol=10e-4
instant=Instant.START,
first_marker_idx=Bow.frog_marker,
second_marker_idx=violon_string.bridge_marker,
)

for g :
max_bound=0+tol, min_bound=0-tol

Plotting bug when an entire phase doesn't have an element

Custom dynamics

Add the capability to create a custom dynamics without changing the dynamics.py file

Add a Lagrange cost on the derivative of the control

I could use this feature, in general it is a nice way to regularize the problems.

Cyclic constraint

Fix the cyclic constraint

Plot bounds on graph

Bounds should be the real ylim of the plot if asked

Deal with quaternions

Change the dict for classes

Dict are fun but confusing. Let's be more verbose on what everything is needed by using classes instead of dict!

Select referential in which minimize_markers_displacement is expressed

Type problem when no data_to_track in minimize_states

File "/home/user/Documents/Programmation/BiorbdOptim/biorbd_optim/penalty.py", line 26, in minimize_states
val = v[states_idx] - data_to_track[t[i], states_idx]
File "/home/user/anaconda3/envs/Trampo_biorbdOptim/lib/python3.7/site-packages/casadi/casadi.py", line 6497, in array
+ "Use an equivalent CasADi function instead of that numpy function.")
Exception: Implicit conversion of symbolic CasADi type to numeric matrix not supported.
This may occur when you pass a CasADi object to a numpy function.
Use an equivalent CasADi function instead of that numpy function.

show_meshes=False

b = BiorbdViz(loaded_model=ocp.nlp[0]["model"], show_meshes=False) b.load_movement(x.T) b.exec()
Casadi cannot open biorbdviz even if show_meshes = False

Allow for modifying constraint and objective function after having created the OCP

Cyclic objective

Fix the cyclic objective