Implementing musculoskeletal model driven by FES

.github/ISSUE_TEMPLATE/enhancement_report.md

@@ -0,0 +1,17 @@
+---
+name: Enhancement report
+about: Create a report to help us improve
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+### Subject of the enhancement
+Give a title to your enhancement here.
+
+### Proposed enhancement
+Describe your enhancement here.
+
+### Additional information
+Any additional information, like libraries versions, etc.

README.md

@@ -1,17 +1,11 @@
-<p align="center">
-    <img
-      src="docs/cocofest_logo.png"
-      alt="logo"
-      width="200"
-    />
-</p>
+<img align="right" width="400" src="docs/cocofest_logo.png">
 
 # COCOFEST
 
 `Cocofest` : Custom Optimal Control Optimization for Functional Electrical Stimulation, is an optimal control program (OCP) package for functional electrical stimulation (FES).
 It is based on the [bioptim](https://github.com/pyomeca/bioptim) framework for the optimal control construction.
 Bioptim uses [biorbd](https://github.com/pyomeca/biorbd) a biomechanics library and benefits from the powerful algorithmic diff provided by [CasADi](https://web.casadi.org/).
-To solve the OCP, the robust solver [Ipopt](https://github.com/coin-or/Ipopt) has been implemented. 
+To solve the OCP, the robust solver [Ipopt](https://github.com/coin-or/Ipopt) has been implemented.
 
 ## Status
 
@@ -37,11 +31,11 @@ To solve the OCP, the robust solver [Ipopt](https://github.com/coin-or/Ipopt) ha
 [Create your own FES OCP](#create-your-own-fes-ocp)
 
 [Examples](#examples)
+- [Musculoskeletal model driven by FES models](#musculoskeletal-model-driven-by-FES-models)
 
 <details>
 <summary><a href="#other-functionalities">Other functionalities</a></summary>
 
-- [With fatigue](#with-fatigue)
 - [Initial value problem](#initital-value-problem)
 - [Summation truncation](#summation-truncation)
 
@@ -56,11 +50,36 @@ Cloning the repository is the first step to be able to use the package.
 
 ## Dependencies
 `Cocofest` relies on several libraries. 
-Based on `bioptim`, the user is invited to directly download the framework from anaconda or from the [sources](https://github.com/pyomeca/bioptim) by cloning the repository
+Follows the steps to install everything you need to use `Cocofest`.
+</br>
+First, you need to create a new conda environment
+```bash
+conda create -n YOUR_ENV_NAME python=3.10
+```
+
+Then, activate the environment
+```bash
+conda activate YOUR_ENV_NAME
+```
+
+This step will allow you to install the dependencies in the environment much quicker
+```bash:
+conda install -cconda-forge conda-libmamba-solver
+```
+
+After, install the dependencies
 ```bash
-conda install -c conda-forge bioptim
+conda install numpy matplotlib pytest casadi biorbd -cconda-forge --solver=libmamba
 ```
-The other [bioptim dependencies](https://github.com/pyomeca/bioptim#dependencies) must be installed as well.
+
+Finally, install the bioptim setup.py file located in your cocofest/external/bioptim folder
+```bash
+cd <path_up_to_cocofest_file>/external/bioptim
+python setup.py install
+```
+
+You are now ready to use `Cocofest`!
+
 
 # Available FES models
 The available FES models are likely to increase so stay tune.
@@ -96,7 +115,7 @@ the minimum and maximum time between two stimulation pulse, the time bimapping
 (If True, will act like a frequency at constant pulse interval).
 
 ```python
-ocp = OcpFes.prepare_ocp(...,
+ocp = OcpFes().prepare_ocp(...,
                          n_stim=10,
                          n_shooting=20,
                          ...,)
@@ -110,39 +129,41 @@ result = ocp.solve()
 
 # Examples
 You can find all the available examples in the [examples](https://github.com/Kev1CO/cocofest/tree/main/examples) file.
+## Musculoskeletal model driven by FES models
+The following example is a musculoskeletal model driven by the Ding2007 FES model.
+The objective function is to reach a 90° forearm position and 0° arm position at the movement end.
+The stimulation last 1s and the stimulation frequency is 10Hz.
+The optimized parameter are each stimulation pulse width.
 
-# Other functionalities
+<p align="center">
+  <img width="500" src=docs/arm_flexion.gif>
+</p>
 
-## With fatigue
-It is possible to compute the models with their fatigue equations or not.
-For example, the "ding2003" model can be use with fatigue DingModelFrequencyWithFatigue or without DingModelFrequency.
-If no fatigue is applied, the fatigue equation will not be added to the model and the muscle force will remain 
-constant during the simulation regardless of the previous stimulation appearance.
 
-```python
-ocp = OcpFes.prepare_ocp(model=DingModelFrequencyWithFatigue())
-```
+
+# Other functionalities
 
 ## Initital value problem
 You can also compute the models form initial value problem.
 For that, use the IvpFes class to build the computed problem.
 
 ```python
-ocp = IvpFes(
-    ding_model=DingModelFrequency(),
-)
+ocp = IvpFes(model=DingModelFrequency(), ...)
 ```
 
 ## Summation truncation
 The summation truncation is an integer parameter that can be added to the model.
 It will truncate the stimulation apparition list used for the calcium summation.
-The integer number defines the stimulation number to keep for this summation.
+The integer number defines the stimulation number to keep prior this summation calculation (e.g only the 5 past stimulation will be included).
 
 ```python
-ocp = OcpFes.prepare_ocp()
+ocp = OcpFes().prepare_ocp(model=DingModelFrequency(sum_stim_truncation=5))
 ```
 
 
 # Citing
 `Cocofest` is not yet published in a journal.
-But if you use `Cocofest` in your research, please kindly cite this package by giving the repository link.
+But if you use `Cocofest` in your research, please kindly cite this zenodo link [10.5281/zenodo.10427934](https://doi.org/10.5281/zenodo.10427934).
+
+# Acknowledgements
+The Cocofest [logo](docs/cocofest_logo.png) has been design by [MaxMV](https://www.instagram.com/max_mv3/)

cocofest/__init__.py

@@ -5,12 +5,12 @@
 from .models.ding2007_with_fatigue import DingModelPulseDurationFrequencyWithFatigue
 from .models.hmed2018 import DingModelIntensityFrequency
 from .models.hmed2018_with_fatigue import DingModelIntensityFrequencyWithFatigue
-from .optimization.fes_multi_start import FunctionalElectricStimulationMultiStart
+from .models.dynamical_model import FesMskModel
 from .optimization.fes_ocp import OcpFes
 from .optimization.fes_identification_ocp import OcpFesId
+from .optimization.fes_ocp_dynamics import OcpFesMsk
 from .integration.ivp_fes import IvpFes
 from .fourier_approx import FourierSeries
-from .read_data import ExtractData
 from .identification.ding2003_force_parameter_identification import DingModelFrequencyForceParameterIdentification
 from .identification.ding2007_force_parameter_identification import (
     DingModelPulseDurationFrequencyForceParameterIdentification,

cocofest/custom_objectives.py

@@ -2,6 +2,7 @@
 This custom objective is to enable the tracking of a curve by a state at all node. Used for sample data control problems
 such as functional electro stimulation
 """
+
 import numpy as np
 from casadi import MX, SX
 
@@ -68,3 +69,39 @@
             return 1 - (force / controller.states[key].cx)
         else:
             raise RuntimeError(f"Minimization type {minimization_type} not implemented")
+
+    @staticmethod
+    def minimize_overall_muscle_fatigue(controller: PenaltyController) -> MX:
+        """
+        Minimize the overall muscle fatigue.
+
+        Parameters
+        ----------
+        controller: PenaltyController
+            The penalty node elements
+
+        Returns
+        -------
+        The sum of each force scaling factor
+        """
+        muscle_name_list = controller.model.bio_model.muscle_names
+        muscle_fatigue = [controller.states["A_" + muscle_name_list[x]].cx for x in range(len(muscle_name_list))]
+        return sum(muscle_fatigue)
+
+    @staticmethod
+    def minimize_overall_muscle_force_production(controller: PenaltyController) -> MX:
+        """
+        Minimize the overall muscle force production.
+
+        Parameters
+        ----------
+        controller: PenaltyController
+            The penalty node elements
+
+        Returns
+        -------
+        The sum of each force
+        """
+        muscle_name_list = controller.model.bio_model.muscle_names
+        muscle_force = [controller.states["F_" + muscle_name_list[x]].cx ** 3 for x in range(len(muscle_name_list))]
+        return sum(muscle_force)

cocofest/integration/ivp_fes.py

@@ -57,12 +57,14 @@ class IvpFes(OptimalControlProgram):
 
     def __init__(
         self,
-        model: DingModelFrequency
-        | DingModelFrequencyWithFatigue
-        | DingModelPulseDurationFrequency
-        | DingModelPulseDurationFrequencyWithFatigue
-        | DingModelIntensityFrequency
-        | DingModelIntensityFrequencyWithFatigue,
+        model: (
+            DingModelFrequency
+            | DingModelFrequencyWithFatigue
+            | DingModelPulseDurationFrequency
+            | DingModelPulseDurationFrequencyWithFatigue
+            | DingModelIntensityFrequency
+            | DingModelIntensityFrequencyWithFatigue
+        ),
         n_stim: int = None,
         n_shooting: int = None,
         final_time: float = None,
@@ -257,12 +259,14 @@ def build_initial_guess_from_ocp(self, ocp):
     @classmethod
     def from_frequency_and_final_time(
         cls,
-        model: DingModelFrequency
-        | DingModelFrequencyWithFatigue
-        | DingModelPulseDurationFrequency
-        | DingModelPulseDurationFrequencyWithFatigue
-        | DingModelIntensityFrequency
-        | DingModelIntensityFrequencyWithFatigue,
+        model: (
+            DingModelFrequency
+            | DingModelFrequencyWithFatigue
+            | DingModelPulseDurationFrequency
+            | DingModelPulseDurationFrequencyWithFatigue
+            | DingModelIntensityFrequency
+            | DingModelIntensityFrequencyWithFatigue
+        ),
         n_shooting: int,
         final_time: float,
         frequency: int | float = None,
@@ -298,12 +302,14 @@ def from_frequency_and_final_time(
     @classmethod
     def from_frequency_and_n_stim(
         cls,
-        model: DingModelFrequency
-        | DingModelFrequencyWithFatigue
-        | DingModelPulseDurationFrequency
-        | DingModelPulseDurationFrequencyWithFatigue
-        | DingModelIntensityFrequency
-        | DingModelIntensityFrequencyWithFatigue,
+        model: (
+            DingModelFrequency
+            | DingModelFrequencyWithFatigue
+            | DingModelPulseDurationFrequency
+            | DingModelPulseDurationFrequencyWithFatigue
+            | DingModelIntensityFrequency
+            | DingModelIntensityFrequencyWithFatigue
+        ),
         n_stim: int,
         n_shooting: int,
         frequency: int | float = None,