New dynamic optimization method, NMPC, external force, cycling

.github/.gitmodules

@@ -0,0 +1,3 @@
+[submodule "external/bioptim"]
+    path = external/bioptim
+	url = https://github.com/Kev1CO/bioptim

.github/workflows/run_tests_win.yml

@@ -30,16 +30,15 @@ jobs:
       - name: Setup environment
         uses: conda-incubator/setup-miniconda@v2
         with:
-          miniforge-variant: Mambaforge
           miniforge-version: latest
           use-mamba: true
           activate-environment: cocofest
           environment-file: environment.yml
 
       - name: Print mamba info
         run: |
-          mamba info
-          mamba list
+          conda info
+          conda list
 
       - name: Install bioptim on Windows
         run: |
@@ -49,19 +48,21 @@ jobs:
           cd ..
 
       - name: Install extra dependencies
-        run: mamba install pytest-cov black pytest pytest-cov codecov packaging -cconda-forge
+        run: conda install pytest-cov black pytest pytest-cov codecov packaging -cconda-forge
 
       - name: Run tests with code coverage
         run: pytest -v --color=yes --cov-report term-missing --cov=cocofest --cov-report=xml:coverage.xml tests/shard${{ matrix.shard }}
 
       - name: Archive coverage report
         id: archive
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: coverage${{ matrix.shard }}
           path: |
             coverage.xml
             .coverage
+          if-no-files-found: error
+          include-hidden-files: true
 
   merge-coverage:
     needs: build
@@ -81,15 +82,15 @@ jobs:
       - name: Setup environment
         uses: conda-incubator/setup-miniconda@v2
         with:
-          miniforge-variant: Mambaforge
           miniforge-version: latest
           use-mamba: true
           activate-environment: cocofest
           environment-file: environment.yml
 
       - name: Print mamba info
         run: |
-          mamba info
+          conda info
+          conda list
 
       - name: Install bioptim on Windows
         run: |
@@ -103,7 +104,7 @@ jobs:
 
       - name: Download all workflow run artifacts
         id: download
-        uses: actions/download-artifact@v3
+        uses: actions/download-artifact@v4
 
       - name: Rename coverage files
         run: |
@@ -126,11 +127,9 @@ jobs:
           coverage xml
           coverage report -m
 
-      - name: Upload coverage to Codecov
-        uses: codecov/codecov-action@v2
+      - uses: codecov/codecov-action@v5
         with:
           file: ./coverage.xml
           flags: unittests
           fail_ci_if_error: true
           verbose: true
-          token: ${{ secrets.CODECOV_TOKEN }}

.gitignore

@@ -6,6 +6,16 @@ __pycache__/
 .pytest_cache
 
 
-cocofest/examples/data/debug.py
+examples/data/debug.py
 
 external/bioptim
+
+example/*.pkl
+example/fes_multibody_cycling/*.pkl
+example/reaching_task/*.pkl
+example/minimize_fatigue/*.pkl
+
+example/minimize_fatigue/result_file/*.pkl
+examples/nmpc/results/*.pkl
+
+*.pkl

.gitmodules

@@ -1,3 +1,3 @@
 [submodule "external/bioptim"]
-    path = external/bioptim
+	path = external/bioptim
 	url = https://github.com/Kev1CO/bioptim

README.md

@@ -54,7 +54,7 @@ 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
+conda create -n YOUR_ENV_NAME python=3.11
 ```
 
 Then, activate the environment

cocofest/__init__.py

@@ -10,21 +10,12 @@
 from .models.dynamical_model import FesMskModel
 from .models.model_maker import ModelMaker
 from .optimization.fes_ocp import OcpFes
-from .optimization.fes_identification_ocp import OcpFesId
-from .optimization.fes_ocp_dynamics import OcpFesMsk
-from .optimization.fes_ocp_nmpc_cyclic import NmpcFes
-from .optimization.fes_ocp_dynamics_nmpc_cyclic import NmpcFesMsk
+from .optimization.fes_id_ocp import OcpFesId
+from .optimization.fes_ocp_multibody import OcpFesMsk
+from .optimization.fes_nmpc import FesNmpc
+from .optimization.fes_nmpc_multibody import FesNmpcMsk
 from .integration.ivp_fes import IvpFes
 from .fourier_approx import FourierSeries
-from .identification.ding2003_force_parameter_identification import (
-    DingModelFrequencyForceParameterIdentification,
-)
-from .identification.ding2007_force_parameter_identification import (
-    DingModelPulseWidthFrequencyForceParameterIdentification,
-)
-from .identification.hmed2018_force_parameter_identification import (
-    DingModelPulseIntensityFrequencyForceParameterIdentification,
-)
 from .dynamics.inverse_kinematics_and_dynamics import (
     get_circle_coord,
     inverse_kinematics_cycling,
@@ -33,3 +24,4 @@
 from .result.plot import PlotCyclingResult
 from .result.pickle import SolutionToPickle
 from .result.animate import PickleAnimate
+from .result.graphics import FES_plot

cocofest/custom_constraints.py

@@ -3,96 +3,28 @@
 """
 
 from casadi import MX, SX, vertcat
-
 from bioptim import PenaltyController
 
-from .models.hmed2018 import DingModelPulseIntensityFrequency
+from .models.dynamical_model import FesMskModel
 
 
 class CustomConstraint:
     @staticmethod
-    def cn_sum(controller: PenaltyController, stim_time: list, model_idx: int = None) -> MX | SX:
-        model = controller.model.muscles_dynamics_model[model_idx] if isinstance(model_idx, int) else controller.model
-        cn_sum_key = model.cn_sum_name
-        km_key = model.km_name
-        intensity_in_model = True if isinstance(model, DingModelPulseIntensityFrequency) else False
-        pulse_intensity_key = model.pulse_intensity_name if intensity_in_model else None
-        pulse_intensity = controller.parameters[pulse_intensity_key].cx if intensity_in_model else None
-        lambda_i = model.get_lambda_i(nb_stim=len(stim_time), pulse_intensity=pulse_intensity)
-        km = controller.states[km_key].cx if model._with_fatigue else model.km_rest
-        r0 = model.get_r0(km=km)
-
-        return controller.controls[cn_sum_key].cx - model.cn_sum_fun(
-            r0=r0, t=controller.time.cx, t_stim_prev=stim_time, lambda_i=lambda_i
-        )
-
-    @staticmethod
-    def cn_sum_identification(controller: PenaltyController, stim_time: list, stim_index: list) -> MX | SX:
-        intensity_in_model = True if isinstance(controller.model, DingModelPulseIntensityFrequency) else False
-        ar, bs, Is, cr = None, None, None, None
-        if intensity_in_model:
-            ar = controller.parameters["ar"].cx if "ar" in controller.parameters.keys() else controller.model.ar
-            bs = controller.parameters["bs"].cx if "bs" in controller.parameters.keys() else controller.model.bs
-            Is = controller.parameters["Is"].cx if "Is" in controller.parameters.keys() else controller.model.Is
-            cr = controller.parameters["cr"].cx if "cr" in controller.parameters.keys() else controller.model.cr
-        lambda_i = (
-            [
-                controller.model.lambda_i_calculation_identification(
-                    controller.parameters["pulse_intensity"].cx[i], ar, bs, Is, cr
-                )
-                for i in stim_index
-            ]
-            if intensity_in_model
-            else [1 for _ in range(len(stim_time))]
-        )
-        km = (
-            controller.parameters["km_rest"].cx
-            if "km_rest" in controller.parameters.keys()
-            else controller.model.km_rest
-        )
-        r0 = km + controller.model.r0_km_relationship
-        return controller.controls["Cn_sum"].cx - controller.model.cn_sum_fun(
-            r0=r0, t=controller.time.cx, t_stim_prev=stim_time, lambda_i=lambda_i
-        )
-
-    @staticmethod
-    def a_calculation(controller: PenaltyController, last_stim_index: int) -> MX | SX:
-        a = controller.states["A"].cx if controller.model.with_fatigue else controller.model.a_scale
-        last_stim_index = 0 if controller.parameters["pulse_width"].cx.shape == (1, 1) else last_stim_index
-        a_calculation = controller.model.a_calculation(
-            a_scale=a,
-            pulse_width=controller.parameters["pulse_width"].cx[last_stim_index],
-        )
-        return controller.controls["A_calculation"].cx - a_calculation
-
-    @staticmethod
-    def a_calculation_msk(controller: PenaltyController, last_stim_index: int, model_idx: int) -> MX | SX:
-        model = controller.model.muscles_dynamics_model[model_idx]
-        muscle_name = model.muscle_name
-        a = controller.states["A_" + muscle_name].cx if model.with_fatigue else model.a_scale
-        last_stim_index = (
-            0 if controller.parameters["pulse_width_" + muscle_name].cx.shape == (1, 1) else last_stim_index
-        )
-        a_calculation = model.a_calculation(
-            a_scale=a,
-            pulse_width=controller.parameters["pulse_width_" + muscle_name].cx[last_stim_index],
-        )
-        return controller.controls["A_calculation_" + muscle_name].cx - a_calculation
-
-    @staticmethod
-    def a_calculation_identification(controller: PenaltyController, last_stim_index: int) -> MX | SX:
-        a = (
-            controller.parameters["a_scale"].cx
-            if "a_scale" in controller.parameters.keys()
-            else controller.model.a_scale
-        )
-        pd0 = controller.parameters["pd0"].cx if "pd0" in controller.parameters.keys() else controller.model.pd0
-        pdt = controller.parameters["pdt"].cx if "pdt" in controller.parameters.keys() else controller.model.pdt
-        last_stim_index = 0 if controller.parameters["pulse_width"].cx.shape == (1, 1) else last_stim_index
-        a_calculation = controller.model.a_calculation_identification(
-            a_scale=a,
-            pulse_width=controller.parameters["pulse_width"].cx[last_stim_index],
-            pd0=pd0,
-            pdt=pdt,
-        )
-        return controller.controls["A_calculation"].cx - a_calculation
+    def pulse_intensity_sliding_window_constraint(
+        controller: PenaltyController, last_stim_idx: int, muscle_name: str = ""
+    ) -> MX | SX:
+        key = "pulse_intensity" + "_" + str(muscle_name) if muscle_name else "pulse_intensity"
+        parameters = [controller.parameters[key].cx[i] for i in range(last_stim_idx + 1)]
+        if isinstance(controller.model, FesMskModel):
+            model = controller.model.muscles_dynamics_model[0]
+        else:
+            model = controller.model
+
+        while len(parameters) < controller.controls[key].cx.shape[0]:
+            min_intensity = model.min_pulse_intensity() if isinstance(model.min_pulse_intensity(), int | float) else 0
+            parameters.insert(0, min_intensity)
+        if len(parameters) > controller.controls[key].cx.shape[0]:
+            size_diff = len(parameters) - controller.controls[key].cx.shape[0]
+            parameters = parameters[size_diff:]
+
+        return controller.controls[key].cx - vertcat(*parameters)

cocofest/custom_objectives.py

@@ -1,78 +1,12 @@
 """
-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
+This custom objective class regroups all the custom objectives that are used in the optimization problem.
 """
 
-import numpy as np
-from casadi import MX, SX, sum1, horzcat
-
+from casadi import MX, vertcat
 from bioptim import PenaltyController
-from .fourier_approx import FourierSeries
 
 
 class CustomObjective:
-    @staticmethod
-    def track_state_from_time(controller: PenaltyController, fourier_coeff: np.ndarray, key: str) -> MX | SX:
-        """
-        Minimize the states variables.
-        By default, this function is quadratic, meaning that it minimizes towards the target.
-        Targets (default=np.zeros()) and indices (default=all_idx) can be specified.
-
-        Parameters
-        ----------
-        controller: PenaltyController
-            The penalty node elements
-        fourier_coeff: np.ndarray
-            The values to aim for
-        key: str
-            The name of the state to minimize
-
-        Returns
-        -------
-        The difference between the two keys
-        """
-        # get the approximated force value from the fourier series at the node time
-        value_from_fourier = FourierSeries().fit_func_by_fourier_series_with_real_coeffs(
-            controller.ocp.node_time(phase_idx=controller.phase_idx, node_idx=controller.t[0]),
-            fourier_coeff,
-            mode="casadi",
-        )
-        return value_from_fourier - controller.states[key].cx
-
-    @staticmethod
-    def track_state_from_time_interpolate(
-        controller: PenaltyController,
-        force: np.ndarray,
-        key: str,
-        minimization_type: str = "least square",
-    ) -> MX:
-        """
-        Minimize the states variables.
-        This function least square.
-        Targets (default=np.zeros()) and indices (default=all_idx) can be specified.
-
-        Parameters
-        ----------
-        controller: PenaltyController
-            The penalty node elements
-        force: np.ndarray
-            The force vector
-        key: str
-            The name of the state to minimize
-        minimization_type: str
-            The type of minimization to perform. Either "least square" or "best fit"
-
-        Returns
-        -------
-        The difference between the two keys
-        """
-        if minimization_type == "least square":
-            return force - controller.states[key].cx
-        elif minimization_type == "best fit":
-            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:
         """
@@ -88,13 +22,10 @@
         The sum of each force scaling factor
         """
         muscle_name_list = controller.model.bio_model.muscle_names
-        muscle_fatigue_rest = horzcat(
-            *[controller.model.bio_stim_model[x].a_rest for x in range(1, len(muscle_name_list) + 1)]
-        )
-        muscle_fatigue = horzcat(
+        muscle_fatigue = vertcat(
             *[controller.states["A_" + muscle_name_list[x]].cx for x in range(len(muscle_name_list))]
         )
-        return sum1(muscle_fatigue_rest) / sum1(muscle_fatigue)
+        return muscle_fatigue
 
     @staticmethod
     def minimize_overall_muscle_force_production(controller: PenaltyController) -> MX:
@@ -111,7 +42,7 @@
         The sum of each force
         """
         muscle_name_list = controller.model.bio_model.muscle_names
-        muscle_force = horzcat(
+        muscle_force = vertcat(
             *[controller.states["F_" + muscle_name_list[x]].cx for x in range(len(muscle_name_list))]
         )
-        return sum1(muscle_force)
+        return muscle_force