diff --git a/.github/workflows/run_tests_win.yml b/.github/workflows/run_tests_win.yml index a7a4a31a..c868d29f 100644 --- a/.github/workflows/run_tests_win.yml +++ b/.github/workflows/run_tests_win.yml @@ -40,6 +40,9 @@ jobs: conda info conda list + - name: Install extra dependencies + run: conda install pip pytest-cov black pytest pytest-cov codecov packaging -cconda-forge + - name: Install bioptim on Windows run: | pwd @@ -47,9 +50,6 @@ jobs: ./bioptim_install_windows.sh 4 ${{ env.PREFIX_WINDOWS }} cd .. - - name: Install extra dependencies - 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 }} diff --git a/cocofest/custom_objectives.py b/cocofest/custom_objectives.py index af4acd68..4360def6 100644 --- a/cocofest/custom_objectives.py +++ b/cocofest/custom_objectives.py @@ -23,7 +23,7 @@ def minimize_overall_muscle_fatigue(controller: PenaltyController) -> MX: ------- The sum of each force scaling factor """ - muscle_name_list = controller.model.bio_model.muscle_names + muscle_name_list = controller.model.muscle_names muscle_model = controller.model.muscles_dynamics_model muscle_fatigue = vertcat( *[ @@ -47,7 +47,7 @@ def minimize_overall_muscle_force_production(controller: PenaltyController) -> M ------- The sum of each force """ - muscle_name_list = controller.model.bio_model.muscle_names + muscle_name_list = controller.model.muscle_names muscle_model = controller.model.muscles_dynamics_model muscle_force = vertcat( *[ @@ -72,7 +72,7 @@ def minimize_overall_stimulation_charge(controller: PenaltyController) -> MX: The sum of each stimulation control """ if isinstance(controller.model, FesMskModel): - muscle_name_list = controller.model.bio_model.muscle_names + muscle_name_list = controller.model.muscle_names if isinstance(controller.model.muscles_dynamics_model[0], DingModelPulseWidthFrequency): stim_charge = vertcat( *[ diff --git a/cocofest/integration/ivp_fes.py b/cocofest/integration/ivp_fes.py index 53a5fc7b..dd3dc838 100644 --- a/cocofest/integration/ivp_fes.py +++ b/cocofest/integration/ivp_fes.py @@ -1,12 +1,10 @@ import numpy as np from bioptim import ( ControlType, - DynamicsList, InitialGuessList, OdeSolver, OptimalControlProgram, ParameterList, - PhaseDynamics, BoundsList, InterpolationType, Solution, @@ -20,6 +18,7 @@ from cocofest.models.ding2007.ding2007_with_fatigue import DingModelPulseWidthFrequencyWithFatigue from cocofest.models.hmed2018.hmed2018 import DingModelPulseIntensityFrequency from cocofest.models.hmed2018.hmed2018_with_fatigue import DingModelPulseIntensityFrequencyWithFatigue +from cocofest.optimization.fes_ocp import OcpFes class IvpFes: @@ -90,9 +89,10 @@ def __init__( numerical_data_time_series, stim_idx_at_node_list = self.model.get_numerical_data_time_series( self.n_shooting, self.final_time ) - self.ode_solver = self.ivp_parameters["ode_solver"] - self._declare_dynamics(numerical_data_time_series) + self.dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=self.ode_solver + ) ( self.x_init, @@ -285,7 +285,7 @@ def _prepare_fake_ocp(self): return OptimalControlProgram( bio_model=[self.model], - dynamics=self.dynamics, + dynamics=self.dynamics_options, n_shooting=self.n_shooting, phase_time=self.final_time, x_init=self.x_init, @@ -318,20 +318,6 @@ def integrate( duplicated_times=duplicated_times, ) - def _declare_dynamics(self, numerical_data_time_series=None): - - self.dynamics = DynamicsList() - self.dynamics.add( - self.model.declare_ding_variables, - dynamic_function=self.model.dynamics, - expand_dynamics=True, - expand_continuity=False, - phase=0, - phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, - numerical_data_timeseries=numerical_data_time_series, - ode_solver=self.ode_solver, - ) - def build_initial_guess_from_ocp(self, ocp, stim_idx_at_node_list=None): """ Build a state, control, parameters and stochastic initial guesses for each phases from a given ocp @@ -341,9 +327,7 @@ def build_initial_guess_from_ocp(self, ocp, stim_idx_at_node_list=None): p = InitialGuessList() s = InitialGuessList() - muscle_name = "_" + ocp.model.muscle_name if ocp.model.muscle_name else "" - for j in range(len(self.model.name_dof)): - key = ocp.model.name_dof[j] + muscle_name + for j, key in enumerate(self.model.name_dofs): x.add(key=key, initial_guess=ocp.model.standard_rest_values()[j], phase=0) if ocp.controls_keys: diff --git a/cocofest/models/ding2003/ding2003.py b/cocofest/models/ding2003/ding2003.py index 2b07193d..a2b753ed 100644 --- a/cocofest/models/ding2003/ding2003.py +++ b/cocofest/models/ding2003/ding2003.py @@ -1,4 +1,4 @@ -from typing import Callable +from typing import Callable, List from math import gcd from fractions import Fraction @@ -6,17 +6,19 @@ from casadi import MX, exp, vertcat from bioptim import ( - ConfigureProblem, DynamicsEvaluation, NonLinearProgram, - OptimalControlProgram, + StateDynamics, + DynamicsFunctions, + OdeSolver, + States, ) from cocofest.models.state_configure import StateConfigure from cocofest.models.fes_model import FesModel -class DingModelFrequency(FesModel): +class DingModelFrequency(FesModel, StateDynamics): """ This is a custom model of the Bioptim package. As CustomModel, some methods are mandatory and must be implemented. to make it work with bioptim. @@ -38,8 +40,9 @@ def __init__( stim_time: list[float] = None, previous_stim: dict = None, sum_stim_truncation: int = 20, + **kwargs, ): - super().__init__() + super().__init__(name=model_name, **kwargs) self._model_name = model_name self._muscle_name = muscle_name self.sum_stim_truncation = sum_stim_truncation @@ -68,6 +71,30 @@ def __init__( self.km_rest = KM_REST_DEFAULT self.fmax = 315.5 # Maximum force (N) at 100 Hz + # ---- Muscle relationship ---- # + self.fes_model = None + self.force_length_relationship = 1 + self.force_velocity_relationship = 1 + self.passive_force_relationship = 0 + + # --- Configure variables --- # + @property + def state_configuration_functions(self) -> List[States | Callable]: + return [StateConfigure().configure_all_muscle_states] + + @property + def control_configuration_functions(self) -> List[States | Callable]: + return [] + + @property + def algebraic_configuration_functions(self) -> List[States | Callable]: + return [] + + @property + def extra_configuration_functions(self) -> List[States | Callable]: + return [] + + # --- Set model parameters --- # def set_a_rest(self, model, a_rest: MX | float): # models is required for bioptim compatibility self.a_rest = a_rest @@ -106,8 +133,8 @@ def serialize(self) -> tuple[Callable, dict]: # ---- Needed for the example ---- # @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return ["Cn" + muscle_name, "F" + muscle_name] @property @@ -155,38 +182,34 @@ def get_lambda_i(nb_stim: int, pulse_intensity: MX | float) -> list[MX | float]: # ---- Model's dynamics ---- # def system_dynamics( self, - cn: MX, - f: MX, - t: MX = None, - t_stim_prev: list[MX] = None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics is the function that describes the models. Parameters ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: list[MX] - The time list of the previous stimulations (s) - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) + time: MX + The system's current node time + states: MX + The state of the system CN, F + controls: MX + The controls of the system, none + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) f_dot = self.f_dot_fun( cn, @@ -194,9 +217,6 @@ def system_dynamics( self.a_rest, self.tau1_rest, self.km_rest, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, ) # Equation n°2 return vertcat(cn_dot, f_dot) @@ -285,9 +305,6 @@ def f_dot_fun( a: MX | float, tau1: MX | float, km: MX | float, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, ) -> MX | float: """ Parameters @@ -302,23 +319,16 @@ def f_dot_fun( The previous step value of time_state_force_no_cross_bridge (s) km: MX | float The previous step value of cross_bridges (unitless) - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) - Returns ------- The value of the derivative force (N) """ return (a * (cn / (km + cn)) - (f / (tau1 + self.tau2 * (cn / (km + cn))))) * ( - force_length_relationship * force_velocity_relationship + passive_force_relationship + self.force_length_relationship * self.force_velocity_relationship + self.passive_force_relationship ) # Equation n°2 - @staticmethod def dynamics( + self, time: MX, states: MX, controls: MX, @@ -326,10 +336,6 @@ def dynamics( algebraic_states: MX, numerical_timeseries: MX, nlp: NonLinearProgram, - fes_model=None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, ) -> DynamicsEvaluation: """ Functional electrical stimulation dynamic @@ -350,57 +356,31 @@ def dynamics( The numerical timeseries of the system nlp: NonLinearProgram A reference to the phase - fes_model: DingModelFrequency - The current phase fes model - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) Returns ------- The derivative of the states in the tuple[MX] format """ - model = fes_model if fes_model else nlp.model + model = self.fes_model if self.fes_model else nlp.model dxdt_fun = model.system_dynamics + dxdt = dxdt_fun( + time=time, + states=states, + controls=controls, + numerical_timeseries=numerical_timeseries, + ) + + defects = None + if isinstance(nlp.dynamics_type.ode_solver, OdeSolver.COLLOCATION): + states_dot_list = [] + for key in model.name_dofs: + states_dot_list.append(DynamicsFunctions.get(nlp.states_dot[key], nlp.states_dot.scaled.cx)) + defects = vertcat(*states_dot_list) - dxdt return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - t=time, - t_stim_prev=numerical_timeseries, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, - ), + dxdt=dxdt, + defects=defects, ) - def declare_ding_variables( - self, - ocp: OptimalControlProgram, - nlp: NonLinearProgram, - numerical_data_timeseries: dict[str, np.ndarray] = None, - contact_type: list = (), - ): - """ - Tell the program which variables are states and controls. - The user is expected to use the ConfigureProblem.configure_xxx functions. - Parameters - ---------- - ocp: OptimalControlProgram - A reference to the ocp - nlp: NonLinearProgram - A reference to the phase - numerical_data_timeseries: dict[str, np.ndarray] - A list of values to pass to the dynamics at each node. Experimental external forces should be included here. - contact_type: list - A list of contact types. This is used to define the contact forces in the dynamics. Not used in this model. - """ - StateConfigure().configure_all_fes_model_states(ocp, nlp, fes_model=self) - ConfigureProblem.configure_dynamics_function(ocp, nlp, dyn_func=self.dynamics) - def _get_additional_previous_stim_time(self): while len(self.previous_stim["time"]) < self.sum_stim_truncation: self.previous_stim["time"].insert(0, -10000000) diff --git a/cocofest/models/ding2003/ding2003_with_fatigue.py b/cocofest/models/ding2003/ding2003_with_fatigue.py index 13dd4038..9cfc111f 100644 --- a/cocofest/models/ding2003/ding2003_with_fatigue.py +++ b/cocofest/models/ding2003/ding2003_with_fatigue.py @@ -3,15 +3,7 @@ from casadi import MX, vertcat import numpy as np -from bioptim import ( - ConfigureProblem, - DynamicsEvaluation, - NonLinearProgram, - OptimalControlProgram, -) - from cocofest.models.ding2003.ding2003 import DingModelFrequency -from cocofest.models.state_configure import StateConfigure class DingModelFrequencyWithFatigue(DingModelFrequency): @@ -101,8 +93,8 @@ def serialize(self) -> tuple[Callable, dict]: # ---- Needed for the example ---- # @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "Cn" + muscle_name, "F" + muscle_name, @@ -139,47 +131,37 @@ def identifiable_parameters(self): # ---- Model's dynamics ---- # def system_dynamics( self, - cn: MX, - f: MX, - a: MX = None, - tau1: MX = None, - km: MX = None, - t: MX = None, - t_stim_prev: MX | float = None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics is the function that describes the models. Parameters ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - a: MX - The value of the scaling factor (unitless) - tau1: MX - The value of the time_state_force_no_cross_bridge (s) - km: MX - The value of the cross_bridges (unitless) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: MX | float - The previous time at which the stimulation was applied (s) - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) + time: MX + The system's current node time + states: MX + The state of the system CN, F, A, Tau1, Km + controls: MX + The controls of the system, none + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + a = states[2] + tau1 = states[3] + km = states[4] + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) # Equation n°1 f_dot = self.f_dot_fun( cn, @@ -187,9 +169,6 @@ def system_dynamics( a, tau1, km, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, ) # Equation n°2 a_dot = self.a_dot_fun(a, f) # Equation n°5 tau1_dot = self.tau1_dot_fun(tau1, f) # Equation n°9 @@ -240,91 +219,3 @@ def km_dot_fun(self, km: MX, f: MX) -> MX | float: The value of the derivative cross_bridges (unitless) """ return -(km - self.km_rest) / self.tau_fat + self.alpha_km * f # Equation n°11 - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F, A, Tau1, Km - controls: MX - The controls of the system, none - parameters: MX - The parameters acting on the system, final time of each phase - algebraic_states: MX - The stochastic variables of the system, none - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: DingModelFrequencyWithFatigue - The current phase fes model - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - a=states[2], - tau1=states[3], - km=states[4], - t=time, - t_stim_prev=numerical_timeseries, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, - ), - defects=None, - ) - - def declare_ding_variables( - self, - ocp: OptimalControlProgram, - nlp: NonLinearProgram, - numerical_data_timeseries: dict[str, np.ndarray] = None, - contact_type: list = (), - ): - """ - Tell the program which variables are states and controls. - The user is expected to use the ConfigureProblem.configure_xxx functions. - Parameters - ---------- - ocp: OptimalControlProgram - A reference to the ocp - nlp: NonLinearProgram - A reference to the phase - numerical_data_timeseries: dict[str, np.ndarray] - A list of values to pass to the dynamics at each node. Experimental external forces should be included here. - contact_type: list - A list of contact types. This is used to define the contact forces in the dynamics. Not used in this model. - """ - StateConfigure().configure_all_fes_model_states(ocp, nlp, fes_model=self) - ConfigureProblem.configure_dynamics_function(ocp, nlp, dyn_func=self.dynamics) diff --git a/cocofest/models/ding2007/ding2007.py b/cocofest/models/ding2007/ding2007.py index 3f2223a1..1e41a478 100644 --- a/cocofest/models/ding2007/ding2007.py +++ b/cocofest/models/ding2007/ding2007.py @@ -1,14 +1,9 @@ -from typing import Callable +from typing import Callable, List -import numpy as np from casadi import MX, vertcat, exp -from bioptim import ( - ConfigureProblem, - DynamicsEvaluation, - NonLinearProgram, - OptimalControlProgram, -) +from bioptim import States + from cocofest.models.ding2003.ding2003 import DingModelFrequency from cocofest.models.state_configure import StateConfigure @@ -66,6 +61,10 @@ def __init__( self.tauc = TAUC_DEFAULT self.fmax = 248 # Maximum force (N) at 100 Hz and 600 us + @property + def control_configuration_functions(self) -> List[States | Callable]: + return [StateConfigure().configure_last_pulse_width] + @property def identifiable_parameters(self): return { @@ -108,41 +107,35 @@ def serialize(self) -> tuple[Callable, dict]: def system_dynamics( self, - cn: MX, - f: MX, - t: MX = None, - t_stim_prev: list[float] | list[MX] = None, - pulse_width: MX = None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics is the function that describes the models. Parameters ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: list[float] | list[MX] - The time list of the previous stimulations (s) - pulse_width: MX - The pulsation duration of the current stimulation (s) - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force coefficient of the muscle (unitless) + time: MX + The system's current node time + states: MX + The state of the system CN, F + controls: MX + The controls of the system, pulse_width + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + pulse_width = controls[0] + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) a_scale = self.a_calculation(a_scale=self.a_scale, pulse_width=pulse_width) f_dot = self.f_dot_fun( @@ -151,9 +144,6 @@ def system_dynamics( a_scale, self.tau1_rest, self.km_rest, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, ) # Equation n°2 from Ding's 2003 article return vertcat(cn_dot, f_dot) @@ -210,90 +200,3 @@ def set_impulse_width(self, value: list[MX]): The pulsation duration list (s) """ self.pulse_width = value - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F, A, Tau1, Km - controls: MX - The controls of the system, none - parameters: MX - The parameters acting on the system, final time of each phase - algebraic_states: MX - The stochastic variables of the system, none - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: DingModelPulseWidthFrequency - The current phase fes model - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force coefficient of the muscle (unitless) - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - t=time, - t_stim_prev=numerical_timeseries, - pulse_width=controls[0], - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, - ), - defects=None, - ) - - def declare_ding_variables( - self, - ocp: OptimalControlProgram, - nlp: NonLinearProgram, - numerical_data_timeseries: dict[str, np.ndarray] = None, - contact_type: list = (), - ): - """ - Tell the program which variables are states and controls. - The user is expected to use the ConfigureProblem.configure_xxx functions. - Parameters - ---------- - ocp: OptimalControlProgram - A reference to the ocp - nlp: NonLinearProgram - A reference to the phase - numerical_data_timeseries: dict[str, np.ndarray] - A list of values to pass to the dynamics at each node. Experimental external forces should be included here. - contact_type: list - A list of contact types. This is used to define the contact forces in the dynamics. Not used in this model. - """ - StateConfigure().configure_all_fes_model_states(ocp, nlp, fes_model=self) - StateConfigure().configure_last_pulse_width(ocp, nlp) - ConfigureProblem.configure_dynamics_function(ocp, nlp, dyn_func=self.dynamics) diff --git a/cocofest/models/ding2007/ding2007_with_fatigue.py b/cocofest/models/ding2007/ding2007_with_fatigue.py index c930f8aa..595e3950 100644 --- a/cocofest/models/ding2007/ding2007_with_fatigue.py +++ b/cocofest/models/ding2007/ding2007_with_fatigue.py @@ -3,14 +3,7 @@ from casadi import MX, vertcat import numpy as np -from bioptim import ( - ConfigureProblem, - DynamicsEvaluation, - NonLinearProgram, - OptimalControlProgram, -) from cocofest.models.ding2007.ding2007 import DingModelPulseWidthFrequency -from cocofest.models.state_configure import StateConfigure class DingModelPulseWidthFrequencyWithFatigue(DingModelPulseWidthFrequency): @@ -59,8 +52,8 @@ def __init__( # ---- Absolutely needed methods ---- # @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "Cn" + muscle_name, "F" + muscle_name, @@ -124,50 +117,38 @@ def serialize(self) -> tuple[Callable, dict]: def system_dynamics( self, - cn: MX, - f: MX, - a: MX = None, - tau1: MX = None, - km: MX = None, - t: MX = None, - t_stim_prev: MX = None, - pulse_width: MX = None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics is the function that describes the models. Parameters ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - a: MX - The value of the scaling factor (unitless) - tau1: MX - The value of the time_state_force_no_cross_bridge (ms) - km: MX - The value of the cross_bridges (unitless) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: MX - The time of the previous stimulation (s) - pulse_width: MX - The time of the impulse (s) - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) + time: MX + The system's current node time + states: MX + The state of the system CN, F, A, Tau1, Km + controls: MX + The controls of the system, pulse_width + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + a = states[2] + tau1 = states[3] + km = states[4] + pulse_width = controls[0] + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) a_scale = self.a_calculation(a_scale=a, pulse_width=pulse_width) @@ -177,9 +158,6 @@ def system_dynamics( a_scale, tau1, km, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, ) # Equation n°2 from Ding's 2003 article a_dot = self.a_dot_fun(a, f) @@ -231,93 +209,3 @@ def km_dot_fun(self, km: MX, f: MX) -> MX | float: The value of the derivative cross_bridges (unitless) """ return -(km - self.km_rest) / self.tau_fat + self.alpha_km * f # Equation n°11 - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F, A, Tau1, Km - controls: MX - The controls of the system, none - parameters: MX - The parameters acting on the system, final time of each phase - algebraic_states: MX - The stochastic variables of the system, none - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: DingModelPulseWidthFrequencyWithFatigue - The current phase fes model - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - a=states[2], - tau1=states[3], - km=states[4], - t=time, - t_stim_prev=numerical_timeseries, - pulse_width=controls[0], - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, - ), - defects=None, - ) - - def declare_ding_variables( - self, - ocp: OptimalControlProgram, - nlp: NonLinearProgram, - numerical_data_timeseries: dict[str, np.ndarray] = None, - contact_type: list = (), - ): - """ - Tell the program which variables are states and controls. - The user is expected to use the ConfigureProblem.configure_xxx functions. - Parameters - ---------- - ocp: OptimalControlProgram - A reference to the ocp - nlp: NonLinearProgram - A reference to the phase - numerical_data_timeseries: dict[str, np.ndarray] - A list of values to pass to the dynamics at each node. Experimental external forces should be included here. - contact_type: list - A list of contact types. This is used to define the contact forces in the dynamics. Not used in this model. - """ - StateConfigure().configure_all_fes_model_states(ocp, nlp, fes_model=self) - StateConfigure().configure_last_pulse_width(ocp, nlp) - ConfigureProblem.configure_dynamics_function(ocp, nlp, dyn_func=self.dynamics) diff --git a/cocofest/models/dynamical_model.py b/cocofest/models/dynamical_model.py index f1ed9f5e..aaef74f9 100644 --- a/cocofest/models/dynamical_model.py +++ b/cocofest/models/dynamical_model.py @@ -1,18 +1,20 @@ -from typing import Callable +from typing import Callable, List import numpy as np from casadi import vertcat, MX, SX, Function from bioptim import ( BiorbdModel, + ConfigureVariables, + ContactType, + DynamicsEvaluation, + DynamicsFunctions, ExternalForceSetTimeSeries, - OptimalControlProgram, NonLinearProgram, - ConfigureProblem, - DynamicsFunctions, - DynamicsEvaluation, - ParameterList, - ContactType, OdeSolver, + OptimalControlProgram, + ParameterList, + StateDynamics, + States, ) from ..models.fes_model import FesModel @@ -27,7 +29,7 @@ ) -class FesMskModel(BiorbdModel): +class FesMskModel(BiorbdModel, StateDynamics): def __init__( self, name: str = None, @@ -41,6 +43,7 @@ def __init__( activate_residual_torque: bool = False, parameters: ParameterList = None, external_force_set: ExternalForceSetTimeSeries = None, + with_contact=False, ): """ The custom model that will be used in the optimal control program for the FES-MSK models @@ -64,8 +67,14 @@ def __init__( parameters: ParameterList The parameters that will be used in the model """ - super().__init__(biorbd_path, parameters=parameters, external_force_set=external_force_set) - self.bio_model = BiorbdModel(biorbd_path, parameters=parameters, external_force_set=external_force_set) + self.with_contact = with_contact + super().__init__( + bio_model=biorbd_path, + parameters=parameters, + external_force_set=external_force_set, + contact_types=[ContactType.RIGID_EXPLICIT] if self.with_contact else [], + ) + self._name = name self.biorbd_path = biorbd_path @@ -84,14 +93,53 @@ def __init__( ) self.muscles_dynamics_model[i].all_stim = self.muscles_dynamics_model[i].previous_stim["time"] + stim_time - self.bio_stim_model = [self.bio_model] + self.muscles_dynamics_model - + self.bio_stim_model = [self] + self.muscles_dynamics_model self.activate_force_length_relationship = activate_force_length_relationship self.activate_force_velocity_relationship = activate_force_velocity_relationship self.activate_passive_force_relationship = activate_passive_force_relationship self.activate_residual_torque = activate_residual_torque self.parameters_list = parameters self.external_forces_set = external_force_set + self.muscles_model = muscles_model + + @property + def state_configuration_functions(self) -> List[States | Callable]: + return [StateConfigure().configure_all_muscle_msk_states, States.Q, States.QDOT] + + @property + def control_configuration_functions(self) -> list: + return [FesMskModel.declare_model_control] + + @property + def algebraic_configuration_functions(self) -> list: + return [] + + @property + def extra_configuration_functions(self) -> list: + return [] + + @property + def contact_types(self) -> List[ContactType]: + return [ContactType.RIGID_EXPLICIT] if self.with_contact else [] + + def get_rigid_contact_forces( + self, + time, + states, + controls, + parameters, + algebraic_states, + numerical_timeseries, + nlp, + ): + q = DynamicsFunctions.get(nlp.states["q"], states) + qdot = DynamicsFunctions.get(nlp.states["qdot"], states) + tau = DynamicsFunctions.get(nlp.controls["tau"], controls) if "tau" in nlp.controls.keys() else MX() + external_forces = nlp.get_external_forces( + "external_forces", states, controls, algebraic_states, numerical_timeseries + ) + + return nlp.model.rigid_contact_forces()(q, qdot, tau, external_forces, nlp.parameters.cx) # ---- Absolutely needed methods ---- # def serialize(self) -> tuple[Callable, dict]: @@ -99,7 +147,7 @@ def serialize(self) -> tuple[Callable, dict]: FesMskModel, { "name": self._name, - "biorbd_path": self.bio_model.path, + "biorbd_path": self.path, "muscles_model": self.muscles_dynamics_model, "stim_time": self.muscles_dynamics_model[0].stim_time, "previous_stim": self.muscles_dynamics_model[0].previous_stim, @@ -113,10 +161,6 @@ def serialize(self) -> tuple[Callable, dict]: ) # ---- Needed for the example ---- # - @property - def name_dof(self) -> tuple[str]: - return self.bio_model.name_dof - def muscle_name_dof(self, index: int = 0) -> list[str]: return self.muscles_dynamics_model[index].name_dof(with_muscle_name=True) @@ -125,14 +169,14 @@ def nb_state(self) -> int: nb_state = 0 for muscle_model in self.muscles_dynamics_model: nb_state += muscle_model.nb_state - nb_state += self.bio_model.nb_q + nb_state += self.nb_q return nb_state @property def name(self) -> None | str: return self._name - def muscle_dynamic( + def dynamics( self, time: MX | SX, states: MX | SX, @@ -141,8 +185,6 @@ def muscle_dynamic( algebraic_states: MX | SX, numerical_data_timeseries: MX | SX, nlp: NonLinearProgram, - muscle_models: list[FesModel], - state_name_list=None, ) -> DynamicsEvaluation: """ The custom dynamics function that provides the derivative of the states: dxdt = f(t, x, u, p, s) @@ -163,15 +205,10 @@ def muscle_dynamic( A list of values to pass to the dynamics at each node. Experimental external forces should be included here. nlp: NonLinearProgram A reference to the phase - muscle_models: list[FesModel] - The list of the muscle models - state_name_list: list[str] - The states names list Returns ------- The derivative of the states in the tuple[MX | SX] format """ - q = DynamicsFunctions.get(nlp.states["q"], states) qdot = DynamicsFunctions.get(nlp.states["qdot"], states) tau = DynamicsFunctions.get(nlp.controls["tau"], controls) if "tau" in nlp.controls.keys() else 0 @@ -184,8 +221,6 @@ def muscle_dynamic( algebraic_states, numerical_data_timeseries, nlp, - muscle_models, - state_name_list, q, qdot, ) @@ -193,11 +228,11 @@ def muscle_dynamic( # You can directly call biorbd function (as for ddq) or call bioptim accessor (as for dq) dq = DynamicsFunctions.compute_qdot(nlp, q, qdot) total_torque = muscles_tau + tau if self.activate_residual_torque else muscles_tau - external_forces = nlp.get_external_forces(states, controls, algebraic_states, numerical_data_timeseries) - with_contact = ( - True if nlp.model.bio_model.contact_names != () else False - ) # TODO: Add a better way of with_contact=True - ddq = nlp.model.forward_dynamics(with_contact=with_contact)( + external_forces = nlp.get_external_forces( + "external_forces", states, controls, algebraic_states, numerical_data_timeseries + ) + + ddq = nlp.model.forward_dynamics(with_contact=self.with_contact)( q, qdot, total_torque, external_forces, parameters ) # q, qdot, tau, external_forces, parameters @@ -205,9 +240,10 @@ def muscle_dynamic( defects = None if isinstance(nlp.dynamics_type.ode_solver, OdeSolver.COLLOCATION): - slope_q = DynamicsFunctions.get(nlp.states_dot["q"], nlp.states_dot.scaled.cx) - slope_qdot = DynamicsFunctions.get(nlp.states_dot["qdot"], nlp.states_dot.scaled.cx) - defects = vertcat(slope_q, slope_qdot) * nlp.dt - vertcat(dq, ddq) * nlp.dt + states_dot_list = [] + for key in nlp.states_dot.keys(): + states_dot_list.append(DynamicsFunctions.get(nlp.states_dot[key], nlp.states_dot.scaled.cx)) + defects = vertcat(*states_dot_list) * nlp.dt - dxdt * nlp.dt return DynamicsEvaluation(dxdt=dxdt, defects=defects) @@ -220,31 +256,29 @@ def muscles_joint_torque( algebraic_states: MX | SX, numerical_data_timeseries: MX | SX, nlp: NonLinearProgram, - muscle_models: list[FesModel], - state_name_list=None, q: MX | SX = None, qdot: MX | SX = None, ): - dxdt_muscle_list = vertcat() muscle_forces = vertcat() muscle_idx_list = [] + state_name_list = nlp.states.keys() - Q = nlp.model.bio_model.q - Qdot = nlp.model.bio_model.qdot + Q = nlp.model.q + Qdot = nlp.model.qdot - updatedModel = nlp.model.bio_model.model.UpdateKinematicsCustom(Q, Qdot) - nlp.model.bio_model.model.updateMuscles(updatedModel, Q, Qdot) - updated_muscle_length_jacobian = nlp.model.bio_model.model.musclesLengthJacobian(updatedModel, Q, False).to_mx() + updatedModel = nlp.model.model.UpdateKinematicsCustom(Q, Qdot) + nlp.model.model.updateMuscles(updatedModel, Q, Qdot) + updated_muscle_length_jacobian = nlp.model.model.musclesLengthJacobian(updatedModel, Q, False).to_mx() updated_muscle_length_jacobian = Function("musclesLengthJacobian", [Q, Qdot], [updated_muscle_length_jacobian])( q, qdot ) bio_muscle_names_at_index = [] - for i in range(len(nlp.model.bio_model.model.muscles())): - bio_muscle_names_at_index.append(nlp.model.bio_model.model.muscle(i).name().to_string()) + for i in range(len(nlp.model.model.muscles())): + bio_muscle_names_at_index.append(nlp.model.model.muscle(i).name().to_string()) - for muscle_model in muscle_models: + for muscle_model in nlp.model.muscles_dynamics_model: muscle_states_idxs = [ i for i in range(len(state_name_list)) if muscle_model.muscle_name in state_name_list[i] ] @@ -268,10 +302,10 @@ def muscles_joint_torque( muscle_force_length_coeff = ( muscle_force_length_coefficient( model=updatedModel, - muscle=nlp.model.bio_model.model.muscle(muscle_idx), + muscle=nlp.model.model.muscle(muscle_idx), q=Q, ) - if nlp.model.activate_force_velocity_relationship + if nlp.model.activate_force_length_relationship else 1 ) muscle_force_length_coeff = Function("muscle_force_length_coeff", [Q, Qdot], [muscle_force_length_coeff])( @@ -281,7 +315,7 @@ def muscles_joint_torque( muscle_force_velocity_coeff = ( muscle_force_velocity_coefficient( model=updatedModel, - muscle=nlp.model.bio_model.model.muscle(muscle_idx), + muscle=nlp.model.model.muscle(muscle_idx), q=Q, qdot=Qdot, ) @@ -295,7 +329,7 @@ def muscles_joint_torque( muscle_passive_force_coeff = ( muscle_passive_force_coefficient( model=updatedModel, - muscle=nlp.model.bio_model.model.muscle(muscle_idx), + muscle=nlp.model.model.muscle(muscle_idx), q=Q, ) if nlp.model.activate_passive_force_relationship @@ -313,6 +347,12 @@ def muscles_joint_torque( if external_force_in_numerical_data_timeseries else numerical_data_timeseries ) + + muscle_model.fes_model = muscle_model + muscle_model.force_length_relationship = muscle_force_length_coeff + muscle_model.force_velocity_relationship = muscle_force_velocity_coeff + muscle_model.passive_force_relationship = muscle_passive_force_coeff + muscle_dxdt = muscle_model.dynamics( time, muscle_states, @@ -321,10 +361,6 @@ def muscles_joint_torque( algebraic_states, fes_numerical_data_timeseries, nlp, - fes_model=muscle_model, - force_length_relationship=muscle_force_length_coeff, - force_velocity_relationship=muscle_force_velocity_coeff, - passive_force_relationship=muscle_passive_force_coeff, ).dxdt dxdt_muscle_list = vertcat(dxdt_muscle_list, muscle_dxdt) @@ -400,8 +436,6 @@ def forces_from_fes_driven( algebraic_states, numerical_timeseries, nlp, - nlp.model.muscles_dynamics_model, - state_name_list, q, qdot, ) @@ -414,12 +448,15 @@ def forces_from_fes_driven( return nlp.model.rigid_contact_forces()(q, qdot, tau, external_forces, nlp.parameters.cx) - def declare_model_variables( - self, + @staticmethod + def declare_model_control( ocp: OptimalControlProgram, nlp: NonLinearProgram, numerical_data_timeseries: dict[str, np.ndarray] = None, contact_type: list = (), + as_states: bool = False, + as_controls: bool = True, + as_algebraic_states: bool = False, ): """ Tell the program which variables are states and controls. @@ -434,35 +471,22 @@ def declare_model_variables( A list of values to pass to the dynamics at each node. Experimental external forces should be included here. contact_type: list The type of contact to be used + as_states: bool + If the variables are states + as_controls: bool + If the variables are controls + as_algebraic_states: bool + If the variables are algebraic states """ - - state_name_list = StateConfigure().configure_all_muscle_states(self.muscles_dynamics_model, ocp, nlp) - ConfigureProblem.configure_q(ocp, nlp, as_states=True, as_controls=False) - state_name_list.append("q") - ConfigureProblem.configure_qdot(ocp, nlp, as_states=True, as_controls=False) - state_name_list.append("qdot") - for muscle_model in self.muscles_dynamics_model: + for muscle_model in nlp.model.muscles_dynamics_model: if isinstance(muscle_model, DingModelPulseWidthFrequency): StateConfigure().configure_last_pulse_width(ocp, nlp, muscle_name=str(muscle_model.muscle_name)) if isinstance(muscle_model, DingModelPulseIntensityFrequency): StateConfigure().configure_pulse_intensity( ocp, nlp, muscle_name=str(muscle_model.muscle_name), truncation=muscle_model.sum_stim_truncation ) - if self.activate_residual_torque: - ConfigureProblem.configure_tau(ocp, nlp, as_states=False, as_controls=True) - - ConfigureProblem.configure_dynamics_function( - ocp, - nlp, - dyn_func=self.muscle_dynamic, - muscle_models=self.muscles_dynamics_model, - state_name_list=state_name_list, - ) - - if ContactType.RIGID_EXPLICIT in contact_type: - ConfigureProblem.configure_rigid_contact_function( - ocp, nlp, self.forces_from_fes_driven, state_name_list=state_name_list - ) + if nlp.model.activate_residual_torque: + ConfigureVariables.configure_tau(ocp, nlp, as_states=False, as_controls=True, as_algebraic_states=False) @staticmethod def _model_sanity( diff --git a/cocofest/models/fes_model.py b/cocofest/models/fes_model.py index bc14b5e5..b0218078 100644 --- a/cocofest/models/fes_model.py +++ b/cocofest/models/fes_model.py @@ -1,14 +1,19 @@ from abc import ABC, abstractmethod -import numpy as np from casadi import MX -from bioptim import NonLinearProgram, OptimalControlProgram +from bioptim import NonLinearProgram class FesModel(ABC): - def __init__(self): + def __init__(self, name, **kwargs): + super().__init__(**kwargs) self.stim_time = None self.previous_stim = None + self._name = name + + @property + def name(self): + return self._name @abstractmethod def set_a_rest(self, model, a_rest: MX | float): @@ -64,15 +69,6 @@ def serialize(self): """ - @abstractmethod - def name_dof(self): - """ - - Returns - ------- - - """ - @abstractmethod def nb_state(self): """ @@ -112,12 +108,10 @@ def with_fatigue(self): @abstractmethod def system_dynamics( self, - cn: MX, - f: MX, - cn_sum: MX, - force_length_relationship: MX | float, - force_velocity_relationship: MX | float, - passive_force_relationship: MX | float, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ): """ @@ -169,9 +163,6 @@ def f_dot_fun( a: MX | float, tau1: MX | float, km: MX | float, - force_length_relationship: MX | float, - force_velocity_relationship: MX | float, - passive_force_relationship: MX | float, ): """ @@ -180,35 +171,16 @@ def f_dot_fun( """ - @staticmethod @abstractmethod def dynamics( + self, time: MX, states: MX, controls: MX, parameters: MX, algebraic_states: MX, - numerical_data_timeseries: MX, - nlp: NonLinearProgram, - fes_model, - force_length_relationship: MX | float, - force_velocity_relationship: MX | float, - passive_force_relationship: MX | float, - ): - """ - - Returns - ------- - - """ - - @abstractmethod - def declare_ding_variables( - self, - ocp: OptimalControlProgram, + numerical_timeseries: MX, nlp: NonLinearProgram, - numerical_data_timeseries: dict[str, np.ndarray] = None, - contact_type: tuple = (), ): """ diff --git a/cocofest/models/hmed2018/hmed2018.py b/cocofest/models/hmed2018/hmed2018.py index 82f970ae..ab0a6866 100644 --- a/cocofest/models/hmed2018/hmed2018.py +++ b/cocofest/models/hmed2018/hmed2018.py @@ -1,14 +1,10 @@ -from typing import Callable +from typing import Callable, List from casadi import MX, vertcat, tanh import numpy as np -from bioptim import ( - ConfigureProblem, - DynamicsEvaluation, - NonLinearProgram, - OptimalControlProgram, -) +from bioptim import States + from cocofest.models.ding2003.ding2003 import DingModelFrequency from cocofest.models.state_configure import StateConfigure @@ -62,6 +58,10 @@ def __init__( self.cr = CR_DEFAULT self.impulse_intensity = None + @property + def control_configuration_functions(self) -> List[States | Callable]: + return [StateConfigure().configure_pulse_intensity] + @property def identifiable_parameters(self): return { @@ -77,7 +77,7 @@ def identifiable_parameters(self): @property def pulse_intensity_name(self): - muscle_name = "_" + self.muscle_name if self.muscle_name else "" + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return "pulse_intensity" + muscle_name def set_ar(self, model, ar: MX | float): @@ -117,41 +117,35 @@ def serialize(self) -> tuple[Callable, dict]: def system_dynamics( self, - cn: MX, - f: MX, - t: MX = None, - t_stim_prev: list[MX] | list[float] = None, - pulse_intensity: list[MX] | list[float] = None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics is the function that describes the models. Parameters ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: list[MX] | list[float] - The previous stimulation time (s) - pulse_intensity: list[MX] | list[float] - The pulsation intensity of the current stimulation (mA) - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) + time: MX + The system's current node time + states: MX + The state of the system CN, F + controls: MX + The controls of the system, pulse_intensity + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + pulse_intensity = controls + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev, pulse_intensity) f_dot = self.f_dot_fun( cn, @@ -159,9 +153,6 @@ def system_dynamics( self.a_rest, self.tau1_rest, self.km_rest, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, ) # Equation n°2 return vertcat(cn_dot, f_dot) @@ -215,93 +206,6 @@ def set_impulse_intensity(self, value: MX): for i in range(value.shape[0]): self.impulse_intensity.append(value[i]) - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model: NonLinearProgram = None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F, A, Tau1, Km - controls: MX - The controls of the system, none - parameters: MX - The parameters acting on the system, final time of each phase - algebraic_states: MX - The stochastic variables of the system, none - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: DingModelPulseIntensityFrequency - The current phase fes model - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - t=time, - t_stim_prev=numerical_timeseries, - pulse_intensity=controls, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, - ), - defects=None, - ) - - def declare_ding_variables( - self, - ocp: OptimalControlProgram, - nlp: NonLinearProgram, - numerical_data_timeseries: dict[str, np.ndarray] = None, - contact_type: list = (), - ): - """ - Tell the program which variables are states and controls. - The user is expected to use the ConfigureProblem.configure_xxx functions. - Parameters - ---------- - ocp: OptimalControlProgram - A reference to the ocp - nlp: NonLinearProgram - A reference to the phase - numerical_data_timeseries: dict[str, np.ndarray] - A list of values to pass to the dynamics at each node. Experimental external forces should be included here. - contact_type: list - A list of contact types. This is used to define the contact forces in the dynamics. Not used in this model. - """ - StateConfigure().configure_all_fes_model_states(ocp, nlp, fes_model=self) - StateConfigure().configure_pulse_intensity(ocp, nlp, truncation=self.sum_stim_truncation) - ConfigureProblem.configure_dynamics_function(ocp, nlp, dyn_func=self.dynamics) - def min_pulse_intensity(self): """ Returns diff --git a/cocofest/models/hmed2018/hmed2018_with_fatigue.py b/cocofest/models/hmed2018/hmed2018_with_fatigue.py index f26b8379..964a4e20 100644 --- a/cocofest/models/hmed2018/hmed2018_with_fatigue.py +++ b/cocofest/models/hmed2018/hmed2018_with_fatigue.py @@ -3,14 +3,7 @@ from casadi import MX, vertcat import numpy as np -from bioptim import ( - ConfigureProblem, - DynamicsEvaluation, - NonLinearProgram, - OptimalControlProgram, -) from cocofest.models.hmed2018.hmed2018 import DingModelPulseIntensityFrequency -from cocofest.models.state_configure import StateConfigure class DingModelPulseIntensityFrequencyWithFatigue(DingModelPulseIntensityFrequency): @@ -60,8 +53,8 @@ def __init__( # ---- Absolutely needed methods ---- # @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "Cn" + muscle_name, "F" + muscle_name, @@ -123,50 +116,38 @@ def serialize(self) -> tuple[Callable, dict]: def system_dynamics( self, - cn: MX, - f: MX, - a: MX = None, - tau1: MX = None, - km: MX = None, - t: MX = None, - t_stim_prev: list[MX] = None, - pulse_intensity: list[MX] | list[float] = None, - force_length_relationship: float | MX = 1, - force_velocity_relationship: float | MX = 1, - passive_force_relationship: MX | float = 0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics is the function that describes the models. Parameters - ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - a: MX - The value of the scaling factor (unitless) - tau1: MX - The value of the time_state_force_no_cross_bridge (ms) - km: MX - The value of the cross_bridges (unitless) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: list[MX] - The previous time at which the stimulation was applied (s) - pulse_intensity: list[MX] | list[float] - The intensity of the stimulations (mA) - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) + ---------- + time: MX + The system's current node time + states: MX + The state of the system CN, F, A, Tau1, Km + controls: MX + The controls of the system, pulse_intensity + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + a = states[2] + tau1 = states[3] + km = states[4] + pulse_intensity = controls # Check + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev, pulse_intensity) f_dot = self.f_dot_fun( cn, @@ -174,9 +155,6 @@ def system_dynamics( a, tau1, km, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, ) # Equation n°2 a_dot = self.a_dot_fun(a, f) # Equation n°5 tau1_dot = self.tau1_dot_fun(tau1, f) # Equation n°9 @@ -227,93 +205,3 @@ def km_dot_fun(self, km: MX, f: MX) -> MX | float: The value of the derivative cross_bridges (unitless) """ return -(km - self.km_rest) / self.tau_fat + self.alpha_km * f # Equation n°11 - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - force_length_relationship: float | MX = 1, - force_velocity_relationship: float | MX = 1, - passive_force_relationship: MX | float = 0, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F, A, Tau1, Km - controls: MX - The controls of the system, none - parameters: MX - The parameters acting on the system, final time of each phase - algebraic_states: MX - The stochastic variables of the system, none - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: DingModelPulseIntensityFrequencyWithFatigue - The current phase fes model - force_length_relationship: MX | float - The force length relationship value (unitless) - force_velocity_relationship: MX | float - The force velocity relationship value (unitless) - passive_force_relationship: MX | float - The passive force relationship value (unitless) - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - a=states[2], - tau1=states[3], - km=states[4], - t=time, - t_stim_prev=numerical_timeseries, - pulse_intensity=controls, - force_length_relationship=force_length_relationship, - force_velocity_relationship=force_velocity_relationship, - passive_force_relationship=passive_force_relationship, - ), - defects=None, - ) - - def declare_ding_variables( - self, - ocp: OptimalControlProgram, - nlp: NonLinearProgram, - numerical_data_timeseries: dict[str, np.ndarray] = None, - contact_type: list = (), - ): - """ - Tell the program which variables are states and controls. - The user is expected to use the ConfigureProblem.configure_xxx functions. - Parameters - ---------- - ocp: OptimalControlProgram - A reference to the ocp - nlp: NonLinearProgram - A reference to the phase - numerical_data_timeseries: dict[str, np.ndarray] - A list of values to pass to the dynamics at each node. Experimental external forces should be included here. - contact_type: list - A list of contact types. This is used to define the contact forces in the dynamics. Not used in this model. - """ - StateConfigure().configure_all_fes_model_states(ocp, nlp, fes_model=self) - StateConfigure().configure_pulse_intensity(ocp, nlp, truncation=self.sum_stim_truncation) - ConfigureProblem.configure_dynamics_function(ocp, nlp, dyn_func=self.dynamics) diff --git a/cocofest/models/marion2009/marion2009.py b/cocofest/models/marion2009/marion2009.py index 1c5d99c4..6152fe31 100644 --- a/cocofest/models/marion2009/marion2009.py +++ b/cocofest/models/marion2009/marion2009.py @@ -2,11 +2,6 @@ from casadi import MX, vertcat - -from bioptim import ( - DynamicsEvaluation, - NonLinearProgram, -) from cocofest.models.ding2007.ding2007 import DingModelFrequency @@ -103,32 +98,35 @@ def angle_scaling_factor(self, theta: MX) -> MX: def system_dynamics( self, - cn: MX, - f: MX, - t: MX = None, - t_stim_prev: list[float] | list[MX] = None, - theta: MX = None, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics incorporating angle dependency. Parameters ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: list[float] | list[MX] - The time list of the previous stimulations (s) - theta: MX - The current knee angle in degrees + time: MX + The system's current node time + states: MX + The state of the system CN, F + controls: MX + The controls of the system, theta + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + theta = controls[0] if controls.shape[0] > 0 else 90 + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) # Apply angle scaling @@ -144,63 +142,3 @@ def system_dynamics( ) return vertcat(cn_dot, f_dot) - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic including angle dependency - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F - controls: MX - The controls of the system: pulse_width, theta - parameters: MX - The parameters acting on the system, final time of each phase - algebraic_states: MX - The stochastic variables of the system, none - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: Marion2009ModelFrequency - The current phase fes model - force_length_relationship: MX | float - The force length relationship value (unitless), not considered for this model - force_velocity_relationship: MX | float - The force velocity relationship value (unitless), not considered for this model - passive_force_relationship: MX | float - The passive force coefficient of the muscle (unitless), not considered for this model - - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - t=time, - t_stim_prev=numerical_timeseries, - theta=controls[0] if controls.shape[0] > 0 else 90, - ), - defects=None, - ) diff --git a/cocofest/models/marion2009/marion2009_modified.py b/cocofest/models/marion2009/marion2009_modified.py index 741d19f0..acef62d5 100644 --- a/cocofest/models/marion2009/marion2009_modified.py +++ b/cocofest/models/marion2009/marion2009_modified.py @@ -2,10 +2,6 @@ from casadi import MX, vertcat -from bioptim import ( - DynamicsEvaluation, - NonLinearProgram, -) from cocofest.models.ding2007.ding2007 import DingModelPulseWidthFrequency @@ -108,35 +104,36 @@ def angle_scaling_factor(self, theta: MX) -> MX: def system_dynamics( self, - cn: MX, - f: MX, - t: MX = None, - t_stim_prev: list[float] | list[MX] = None, - pulse_width: MX = None, - theta: MX = None, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics incorporating angle dependency. Parameters ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: list[float] | list[MX] - The time list of the previous stimulations (s) - pulse_width: MX - The pulsation duration of the current stimulation (s) - theta: MX - The current knee angle in degrees + time: MX + The system's current node time + states: MX + The state of the system CN, F + controls: MX + The controls of the system, pulse_width, theta + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + pulse_width = controls[0] + theta = controls[1] if controls.shape[0] > 1 else 90 + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) # Calculate base a_scale with pulse width dependency @@ -155,64 +152,3 @@ def system_dynamics( ) return vertcat(cn_dot, f_dot) - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - force_length_relationship: MX | float = 1, - force_velocity_relationship: MX | float = 1, - passive_force_relationship: MX | float = 0, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic including angle dependency - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F - controls: MX - The controls of the system: pulse_width, theta - parameters: MX - The parameters acting on the system, final time of each phase - algebraic_states: MX - The stochastic variables of the system, none - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: Marion2009ModelPulseWidthFrequency - The current phase fes model - force_length_relationship: MX | float - The force length relationship value (unitless), not considered for this model - force_velocity_relationship: MX | float - The force velocity relationship value (unitless), not considered for this model - passive_force_relationship: MX | float - The passive force coefficient of the muscle (unitless), not considered for this model - - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - t=time, - t_stim_prev=numerical_timeseries, - pulse_width=controls[0], - theta=controls[1] if controls.shape[0] > 1 else 90, - ), - defects=None, - ) diff --git a/cocofest/models/marion2009/marion2009_modified_with_fatigue.py b/cocofest/models/marion2009/marion2009_modified_with_fatigue.py index aa7e9e3e..427b4b24 100644 --- a/cocofest/models/marion2009/marion2009_modified_with_fatigue.py +++ b/cocofest/models/marion2009/marion2009_modified_with_fatigue.py @@ -3,10 +3,6 @@ import numpy as np from casadi import MX, vertcat -from bioptim import ( - DynamicsEvaluation, - NonLinearProgram, -) from .marion2009_modified import Marion2009ModelPulseWidthFrequency @@ -58,8 +54,8 @@ def __init__( self.alpha_km = ALPHA_KM_DEFAULT @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self, with_muscle_name: bool = False) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "Cn" + muscle_name, "F" + muscle_name, @@ -107,44 +103,39 @@ def serialize(self) -> tuple[Callable, dict]: def system_dynamics( self, - cn: MX, - f: MX, - a: MX = None, - tau1: MX = None, - km: MX = None, - t: MX = None, - t_stim_prev: list[float] | list[MX] = None, - pulse_width: MX = None, - theta: MX = None, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics incorporating both angle dependency and fatigue states. Parameters ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - a: MX - The value of the scaling factor (unitless) - tau1: MX - The value of the time_state_force_no_cross_bridge (ms) - km: MX - The value of the cross_bridges (unitless) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: list[float] | list[MX] - The time list of the previous stimulations (s) - pulse_width: MX - The pulsation duration of the current stimulation (s) - theta: MX - The current knee angle in degrees + time: MX + The system's current node time + states: MX + The state of the system CN, F + controls: MX + The controls of the system, pulse_width, theta + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + a = states[2] + tau1 = states[3] + km = states[4] + pulse_width = controls[0] + theta = controls[1] if controls.shape[0] > 1 else 90 + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) # Calculate base a_scale with pulse width dependency @@ -213,58 +204,3 @@ def km_dot_fun(self, km: MX, f: MX) -> MX | float: The value of the derivative cross_bridges (unitless) """ return -(km - self.km_rest) / self.tau_fat + self.alpha_km * f - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic including angle dependency and fatigue - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F, A, Tau1, Km - controls: MX - The controls of the system: pulse_width, theta - parameters: MX - The parameters acting on the system, final time of each phase - algebraic_states: MX - The stochastic variables of the system, none - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: Marion2009ModelPulseWidthFrequencyWithFatigue - The current phase fes model - - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - a=states[2], - tau1=states[3], - km=states[4], - t=time, - t_stim_prev=numerical_timeseries, - pulse_width=controls[0], - theta=controls[1] if controls.shape[0] > 1 else 90, - ), - defects=None, - ) diff --git a/cocofest/models/marion2009/marion2009_with_fatigue.py b/cocofest/models/marion2009/marion2009_with_fatigue.py index 948c385d..c4a919df 100644 --- a/cocofest/models/marion2009/marion2009_with_fatigue.py +++ b/cocofest/models/marion2009/marion2009_with_fatigue.py @@ -3,10 +3,6 @@ import numpy as np from casadi import MX, vertcat -from bioptim import ( - DynamicsEvaluation, - NonLinearProgram, -) from .marion2009 import Marion2009ModelFrequency @@ -55,8 +51,8 @@ def __init__( self.alpha_km = ALPHA_KM_DEFAULT @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self, with_muscle_name: bool = False) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "Cn" + muscle_name, "F" + muscle_name, @@ -104,41 +100,38 @@ def serialize(self) -> tuple[Callable, dict]: def system_dynamics( self, - cn: MX, - f: MX, - a: MX = None, - tau1: MX = None, - km: MX = None, - t: MX = None, - t_stim_prev: list[float] | list[MX] = None, - theta: MX = None, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics incorporating both angle dependency and fatigue states. Parameters ---------- - cn: MX - The value of the ca_troponin_complex (unitless) - f: MX - The value of the force (N) - a: MX - The value of the scaling factor (unitless) - tau1: MX - The value of the time_state_force_no_cross_bridge (ms) - km: MX - The value of the cross_bridges (unitless) - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: list[float] | list[MX] - The time list of the previous stimulations (s) - theta: MX - The current knee angle in degrees + time: MX + The system's current node time + states: MX + The state of the system CN, F, A, Tau1, Km + controls: MX + The controls of the system, theta + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + a = states[2] + tau1 = states[3] + km = states[4] + theta = controls[0] if controls.shape[0] > 0 else 90 + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) # Apply angle scaling @@ -204,57 +197,3 @@ def km_dot_fun(self, km: MX, f: MX) -> MX | float: The value of the derivative cross_bridges (unitless) """ return -(km - self.km_rest) / self.tau_fat + self.alpha_km * f - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic including angle dependency and fatigue - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F, A, Tau1, Km - controls: MX - The controls of the system: pulse_width, theta - parameters: MX - The parameters acting on the system, final time of each phase - algebraic_states: MX - The stochastic variables of the system, none - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: Marion2009ModelFrequencyWithFatigue - The current phase fes model - - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - a=states[2], - tau1=states[3], - km=states[4], - t=time, - t_stim_prev=numerical_timeseries, - theta=controls[0] if controls.shape[0] > 0 else 90, - ), - defects=None, - ) diff --git a/cocofest/models/marion2013/marion2013.py b/cocofest/models/marion2013/marion2013.py index 0ac91cf5..45295f39 100644 --- a/cocofest/models/marion2013/marion2013.py +++ b/cocofest/models/marion2013/marion2013.py @@ -2,10 +2,6 @@ from casadi import MX, vertcat, exp, cos, pi import numpy as np -from bioptim import ( - DynamicsEvaluation, - NonLinearProgram, -) from cocofest.models.marion2009.marion2009 import Marion2009ModelFrequency @@ -64,8 +60,8 @@ def __init__( self.FM = FM_DEFAULT @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self, with_muscle_name: bool = False) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "Cn" + muscle_name, "F" + muscle_name, @@ -126,38 +122,37 @@ def calculate_acceleration(self, theta: MX, lambda_angle: MX, f: MX, Fload: MX) def system_dynamics( self, - cn: MX, - f: MX, - theta: MX, - dtheta_dt: MX, - t: MX = None, - t_stim_prev: MX = None, - Fload: MX = 0.0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics that describes the model. Parameters ---------- - cn: MX - The normalized Ca2+-troponin complex concentration - f: MX - The instantaneous force - theta: MX - The flexion angle - dtheta_dt: MX - The angular velocity - t: MX - The current time at which the dynamics is evaluated (ms) - t_stim_prev: MX - The time of the previous stimulation (ms) - Fload: MX - External load force (N) + time: MX + The system's current node time + states: MX + The state of the system CN, F, theta, dtheta + controls: MX + The controls of the system, Fload + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + theta = states[2] + dtheta_dt = states[3] + Fload = controls[0] if controls.shape[0] > 0 else 0.0 + t_stim_prev = numerical_timeseries + cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) # Similar to Ding's model calculation # Calculate Marion-specific force scaling terms @@ -182,56 +177,3 @@ def system_dynamics( Fload, ) return vertcat(cn_dot, f_dot, dtheta_dt, d2theta_dt2) - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F, theta, dtheta_dt - controls: MX - The controls of the system - parameters: MX - The parameters acting on the system - algebraic_states: MX - The stochastic variables of the system - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: Marion2013ModelFrequency - The current phase fes model - - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - theta=states[2], - dtheta_dt=states[3], - t=time, - t_stim_prev=numerical_timeseries, - Fload=controls[0] if controls.shape[0] > 0 else 0.0, - ), - defects=None, - ) diff --git a/cocofest/models/marion2013/marion2013_modified.py b/cocofest/models/marion2013/marion2013_modified.py index 3004f758..ede8afb4 100644 --- a/cocofest/models/marion2013/marion2013_modified.py +++ b/cocofest/models/marion2013/marion2013_modified.py @@ -68,8 +68,8 @@ def __init__( self.FM = FM_DEFAULT @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self, with_muscle_name: bool = False) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "Cn" + muscle_name, "F" + muscle_name, @@ -130,41 +130,38 @@ def calculate_acceleration(self, theta: MX, lambda_angle: MX, f: MX, Fload: MX) def system_dynamics( self, - cn: MX, - f: MX, - theta: MX, - dtheta_dt: MX, - t: MX = None, - t_stim_prev: MX = None, - pulse_width: MX = None, - Fload: MX = 0.0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics that describes the model. Parameters ---------- - cn: MX - The normalized Ca2+-troponin complex concentration - f: MX - The instantaneous force - theta: MX - The flexion angle - dtheta_dt: MX - The angular velocity - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: MX - The time of the previous stimulation (s) - pulse_width: MX - The pulse width of the stimulation (s) - Fload: MX - External load force (N) + time: MX + The system's current node time + states: MX + The state of the system CN, F, theta, dtheta + controls: MX + The controls of the system, pulse_width, Fload + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + theta = states[2] + dtheta_dt = states[3] + pulse_width = controls[0] + Fload = controls[1] if controls.shape[0] > 1 else 0.0 + t_stim_prev = numerical_timeseries + # Get CN dynamics from Ding model cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) @@ -194,57 +191,3 @@ def system_dynamics( ) return vertcat(cn_dot, f_dot, dtheta_dt, d2theta_dt2) - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system CN, F, theta, dtheta_dt - controls: MX - The controls of the system - parameters: MX - The parameters acting on the system - algebraic_states: MX - The stochastic variables of the system - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: Marion2013ModelPulseWidthFrequency - The current phase fes model - - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - theta=states[2], - dtheta_dt=states[3], - t=time, - t_stim_prev=numerical_timeseries, - pulse_width=controls[0], - Fload=controls[1] if controls.shape[0] > 1 else 0.0, - ), - defects=None, - ) diff --git a/cocofest/models/marion2013/marion2013_modified_with_fatigue.py b/cocofest/models/marion2013/marion2013_modified_with_fatigue.py index 282b7842..4d639d70 100644 --- a/cocofest/models/marion2013/marion2013_modified_with_fatigue.py +++ b/cocofest/models/marion2013/marion2013_modified_with_fatigue.py @@ -2,10 +2,6 @@ from casadi import MX, vertcat import numpy as np -from bioptim import ( - DynamicsEvaluation, - NonLinearProgram, -) from .marion2013_modified import Marion2013ModelPulseWidthFrequency @@ -62,8 +58,8 @@ def __init__( self.beta_a = BETA_A_DEFAULT @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self, with_muscle_name: bool = False) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "Cn" + muscle_name, "F" + muscle_name, @@ -166,50 +162,41 @@ def km_dot_fun(self, km: MX, f: MX, velocity: MX) -> MX | float: def system_dynamics( self, - cn: MX, - f: MX, - theta: MX, - dtheta_dt: MX, - a: MX, - tau1: MX, - km: MX, - t: MX = None, - t_stim_prev: MX = None, - pulse_width: MX = None, - Fload: MX = 0.0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics including fatigue effects. Parameters ---------- - cn: MX - The normalized Ca2+-troponin complex concentration - f: MX - The instantaneous force near ankle - theta: MX - The knee flexion angle - dtheta_dt: MX - The angular velocity - a: MX - The current a value at 90° affected by fatigue - tau1: MX - The current tau1 value affected by fatigue - km: MX - The current Km value affected by fatigue - t: MX - The current time at which the dynamics is evaluated (s) - t_stim_prev: MX - The time of the previous stimulation (s) - pulse_width: MX - The pulse width of the stimulation (s) - Fload: MX - External load force (N) + time: MX + The system's current node time + states: MX + The state of the system CN, F, theta, dtheta, A, Tau1, Km + controls: MX + The controls of the system, pulse_width, Fload + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + theta = states[2] + dtheta_dt = states[3] + a = states[4] + tau1 = states[5] + km = states[6] + pulse_width = controls[0] + Fload = controls[1] if controls.shape[0] > 1 else 0.0 + t_stim_prev = numerical_timeseries + # Get CN dynamics from Ding model cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) @@ -244,60 +231,3 @@ def system_dynamics( tau1_dot = self.tau1_dot_fun(tau1, f, dtheta_dt) return vertcat(cn_dot, f_dot, dtheta_dt, d2theta_dt2, a_dot, tau1_dot, km_dot) - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic with fatigue - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system including fatigue states - controls: MX - The controls of the system - parameters: MX - The parameters acting on the system - algebraic_states: MX - The stochastic variables of the system - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: Marion2013ModelFrequencyWithFatigue - The current phase fes model - - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - theta=states[2], - dtheta_dt=states[3], - a=states[4], - tau1=states[5], - km=states[6], - t=time, - t_stim_prev=numerical_timeseries, - pulse_width=controls[0], - Fload=controls[1] if controls.shape[0] > 1 else 0.0, - ), - defects=None, - ) diff --git a/cocofest/models/marion2013/marion2013_with_fatigue.py b/cocofest/models/marion2013/marion2013_with_fatigue.py index 2cf11dce..beeeb0af 100644 --- a/cocofest/models/marion2013/marion2013_with_fatigue.py +++ b/cocofest/models/marion2013/marion2013_with_fatigue.py @@ -2,10 +2,6 @@ from casadi import MX, vertcat import numpy as np -from bioptim import ( - DynamicsEvaluation, - NonLinearProgram, -) from .marion2013 import Marion2013ModelFrequency @@ -59,8 +55,8 @@ def __init__( self.beta_a = BETA_A_DEFAULT @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self, with_muscle_name: bool = False) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "Cn" + muscle_name, "F" + muscle_name, @@ -163,47 +159,40 @@ def km_dot_fun(self, km: MX, f: MX, velocity: MX) -> MX | float: def system_dynamics( self, - cn: MX, - f: MX, - theta: MX, - dtheta_dt: MX, - a: MX, - tau1: MX, - km: MX, - t: MX = None, - t_stim_prev: MX = None, - Fload: MX = 0.0, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics including fatigue effects. Parameters ---------- - cn: MX - The normalized Ca2+-troponin complex concentration - f: MX - The instantaneous force near ankle - theta: MX - The knee flexion angle - dtheta_dt: MX - The angular velocity - a: MX - The current a value at 90° affected by fatigue - km: MX - The current Km value affected by fatigue - tau1: MX - The current tau1 value affected by fatigue - t: MX - The current time at which the dynamics is evaluated (ms) - t_stim_prev: MX - The time of the previous stimulation (ms) - Fload: MX - External load force (N) + time: MX + The system's current node time + states: MX + The state of the system CN, F, theta, dtheta, A, Tau1, Km + controls: MX + The controls of the system, Fload + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The value of the derivative of each state dx/dt at the current time t """ + t = time + cn = states[0] + f = states[1] + theta = states[2] + dtheta_dt = states[3] + a = states[4] + tau1 = states[5] + km = states[6] + Fload = controls[0] if controls.shape[0] > 0 else 0.0 + t_stim_prev = numerical_timeseries + # Get CN dynamics from Ding model cn_dot = self.calculate_cn_dot(cn, t, t_stim_prev) @@ -235,59 +224,3 @@ def system_dynamics( tau1_dot = self.tau1_dot_fun(tau1, f, dtheta_dt) return vertcat(cn_dot, f_dot, dtheta_dt, d2theta_dt2, a_dot, tau1_dot, km_dot) - - @staticmethod - def dynamics( - time: MX, - states: MX, - controls: MX, - parameters: MX, - algebraic_states: MX, - numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - ) -> DynamicsEvaluation: - """ - Functional electrical stimulation dynamic with fatigue - - Parameters - ---------- - time: MX - The system's current node time - states: MX - The state of the system including fatigue states - controls: MX - The controls of the system - parameters: MX - The parameters acting on the system - algebraic_states: MX - The stochastic variables of the system - numerical_timeseries: MX - The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: Marion2013ModelFrequencyWithFatigue - The current phase fes model - - Returns - ------- - The derivative of the states in the tuple[MX] format - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - cn=states[0], - f=states[1], - theta=states[2], - dtheta_dt=states[3], - a=states[4], - tau1=states[5], - km=states[6], - t=time, - t_stim_prev=numerical_timeseries, - Fload=controls[0] if controls.shape[0] > 0 else 0.0, - ), - defects=None, - ) diff --git a/cocofest/models/state_configure.py b/cocofest/models/state_configure.py index cd594343..065863df 100644 --- a/cocofest/models/state_configure.py +++ b/cocofest/models/state_configure.py @@ -1,5 +1,5 @@ from bioptim import ( - ConfigureProblem, + ConfigureVariables, NonLinearProgram, OptimalControlProgram, ) @@ -49,7 +49,7 @@ def configure_ca_troponin_complex( muscle_name = "_" + muscle_name if muscle_name else "" name = "Cn" + muscle_name name_cn = [name] - ConfigureProblem.configure_new_variable( + ConfigureVariables.configure_new_variable( name, name_cn, ocp, @@ -89,7 +89,7 @@ def configure_force( muscle_name = "_" + muscle_name if muscle_name else "" name = "F" + muscle_name name_f = [name] - ConfigureProblem.configure_new_variable( + ConfigureVariables.configure_new_variable( name, name_f, ocp, @@ -129,7 +129,7 @@ def configure_scaling_factor( muscle_name = "_" + muscle_name if muscle_name else "" name = "A" + muscle_name name_a = [name] - return ConfigureProblem.configure_new_variable( + return ConfigureVariables.configure_new_variable( name, name_a, ocp, @@ -169,7 +169,7 @@ def configure_time_state_force_no_cross_bridge( muscle_name = "_" + muscle_name if muscle_name else "" name = "Tau1" + muscle_name name_tau1 = [name] - return ConfigureProblem.configure_new_variable( + return ConfigureVariables.configure_new_variable( name, name_tau1, ocp, @@ -209,7 +209,7 @@ def configure_cross_bridges( muscle_name = "_" + muscle_name if muscle_name else "" name = "Km" + muscle_name name_km = [name] - return ConfigureProblem.configure_new_variable( + return ConfigureVariables.configure_new_variable( name, name_km, ocp, @@ -236,7 +236,7 @@ def configure_cn_sum(ocp, nlp, muscle_name: str = None): muscle_name = "_" + muscle_name if muscle_name else "" name = "Cn_sum" + muscle_name name_cn_sum = [name] - return ConfigureProblem.configure_new_variable(name, name_cn_sum, ocp, nlp, as_states=False, as_controls=True) + return ConfigureVariables.configure_new_variable(name, name_cn_sum, ocp, nlp, as_states=False, as_controls=True) @staticmethod def configure_a_calculation(ocp, nlp, muscle_name: str = None): @@ -255,7 +255,7 @@ def configure_a_calculation(ocp, nlp, muscle_name: str = None): muscle_name = "_" + muscle_name if muscle_name else "" name = "A_calculation" + muscle_name name_cn_sum = [name] - return ConfigureProblem.configure_new_variable(name, name_cn_sum, ocp, nlp, as_states=False, as_controls=True) + return ConfigureVariables.configure_new_variable(name, name_cn_sum, ocp, nlp, as_states=False, as_controls=True) @staticmethod def configure_muscle_activation( @@ -287,7 +287,7 @@ def configure_muscle_activation( muscle_name = "_" + muscle_name if muscle_name else "" name = "a" + muscle_name name_a = [name] - ConfigureProblem.configure_new_variable( + ConfigureVariables.configure_new_variable( name, name_a, ocp, @@ -327,7 +327,7 @@ def configure_fatigue_state( muscle_name = "_" + muscle_name if muscle_name else "" name = "mu" + muscle_name name_mu = [name] - ConfigureProblem.configure_new_variable( + ConfigureVariables.configure_new_variable( name, name_mu, ocp, @@ -367,7 +367,7 @@ def configure_angle( muscle_name = "_" + muscle_name if muscle_name else "" name = "theta" + muscle_name name_theta = [name] - ConfigureProblem.configure_new_variable( + ConfigureVariables.configure_new_variable( name, name_theta, ocp, @@ -407,7 +407,7 @@ def configure_angular_velocity( muscle_name = "_" + muscle_name if muscle_name else "" name = "dtheta_dt" + muscle_name name_dtheta_dt = [name] - ConfigureProblem.configure_new_variable( + ConfigureVariables.configure_new_variable( name, name_dtheta_dt, ocp, @@ -432,12 +432,12 @@ def configure_last_pulse_width(ocp, nlp, muscle_name: str = None): muscle_name = "_" + muscle_name if muscle_name else "" name = "last_pulse_width" + muscle_name last_pulse_width = [name] - return ConfigureProblem.configure_new_variable( + return ConfigureVariables.configure_new_variable( name, last_pulse_width, ocp, nlp, as_states=False, as_controls=True ) @staticmethod - def configure_pulse_intensity(ocp, nlp, muscle_name: str = None, truncation: int = 20): + def configure_pulse_intensity(ocp, nlp, muscle_name: str = None, truncation: int = None): """ Configure the pulse intensity control for the Ding model @@ -448,15 +448,17 @@ def configure_pulse_intensity(ocp, nlp, muscle_name: str = None, truncation: int nlp: NonLinearProgram A reference to the phase """ + muscle_name = nlp.model.muscle_name if muscle_name is None else muscle_name + truncation = nlp.model.sum_stim_truncation if truncation is None else truncation muscle_name = "_" + muscle_name if muscle_name else "" name = "pulse_intensity" + muscle_name pulse_intensity = [str(i) for i in range(truncation)] - return ConfigureProblem.configure_new_variable( + return ConfigureVariables.configure_new_variable( name, pulse_intensity, ocp, nlp, as_states=False, as_controls=True ) @staticmethod - def configure_intensity(ocp, nlp, muscle_name: str = None): + def configure_intensity(ocp, nlp): """ Configure the intensity control for the Veltink1992 model @@ -467,31 +469,43 @@ def configure_intensity(ocp, nlp, muscle_name: str = None): nlp: NonLinearProgram A reference to the phase """ + muscle_name = nlp.model.muscle_name muscle_name = "_" + muscle_name if muscle_name else "" name = "I" + muscle_name pulse_intensity = [name] - return ConfigureProblem.configure_new_variable( + return ConfigureVariables.configure_new_variable( name, pulse_intensity, ocp, nlp, as_states=False, as_controls=True ) - def configure_all_muscle_states(self, muscles_dynamics_model, ocp, nlp): - state_name_list = [] - for muscle_dynamics_model in muscles_dynamics_model: - for state_key in muscle_dynamics_model.name_dof: - if state_key in self.state_dictionary.keys(): - self.state_dictionary[state_key]( + @staticmethod + def configure_all_muscle_states(ocp, nlp): + for state_key in nlp.model.name_dofs: + if state_key in StateConfigure().state_dictionary.keys(): + StateConfigure().state_dictionary[state_key]( + ocp=ocp, + nlp=nlp, + as_states=True, + as_controls=False, + muscle_name=nlp.model.muscle_name, + ) + + @staticmethod + def configure_all_muscle_msk_states(ocp, nlp): + for muscle_dynamics_model in nlp.model.muscles_dynamics_model: + for state_key in muscle_dynamics_model.name_dofs: + separator = "_" + key = state_key.split(separator, 1)[0] + if key in StateConfigure().state_dictionary.keys(): + StateConfigure().state_dictionary[key]( ocp=ocp, nlp=nlp, as_states=True, as_controls=False, muscle_name=muscle_dynamics_model.muscle_name, ) - state_name_list.append(state_key + "_" + muscle_dynamics_model.muscle_name) - - return state_name_list def configure_all_fes_model_states(self, ocp, nlp, fes_model): - for state_key in fes_model.name_dof: + for state_key in fes_model.name_dofs: if state_key in self.state_dictionary.keys(): self.state_dictionary[state_key]( ocp=ocp, diff --git a/cocofest/models/veltink1992/veltink1992.py b/cocofest/models/veltink1992/veltink1992.py index b5f8d412..80243479 100644 --- a/cocofest/models/veltink1992/veltink1992.py +++ b/cocofest/models/veltink1992/veltink1992.py @@ -1,18 +1,14 @@ -from typing import Callable +from typing import Callable, List from casadi import MX, vertcat import numpy as np -from bioptim import ( - ConfigureProblem, - DynamicsEvaluation, - NonLinearProgram, - OptimalControlProgram, -) +from bioptim import StateDynamics, States, OdeSolver, DynamicsEvaluation, NonLinearProgram, DynamicsFunctions from cocofest.models.state_configure import StateConfigure +from cocofest.models.fes_model import FesModel -class VeltinkModelPulseIntensity: +class VeltinkModelPulseIntensity(StateDynamics): """ This is a custom model implementing the muscle activation dynamics from: @@ -44,9 +40,32 @@ def __init__( self.I_threshold = I_threshold if I_threshold is not None else I_THRESHOLD_DEFAULT self.I_saturation = I_saturation if I_saturation is not None else I_SATURATION_DEFAULT + self.contact_types = () + + @property + def name(self): + return self._name + + # --- Configure variables --- # + @property + def state_configuration_functions(self) -> List[States | Callable]: + return [StateConfigure().configure_all_muscle_states] + + @property + def control_configuration_functions(self) -> List[States | Callable]: + return [StateConfigure().configure_intensity] + + @property + def algebraic_configuration_functions(self) -> List[States | Callable]: + return [] + + @property + def extra_configuration_functions(self) -> List[States | Callable]: + return [] + @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self, with_muscle_name: bool = False) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return ["a" + muscle_name] # Only muscle activation state @property @@ -128,30 +147,38 @@ def get_muscle_activation(self, a: MX, u: MX) -> MX: def system_dynamics( self, - a: MX, - I: MX, + time: MX, + states: MX, + controls: MX, + numerical_timeseries: MX, ) -> MX: """ The system dynamics implementing equation (4) for muscle activation. Parameters ---------- - a: MX - Muscle activation state (unitless) - I: MX - Stimulation current amplitude (mA) + time: MX + The system's current node time + states: MX + The state of the system a + controls: MX + The controls of the system, I + numerical_timeseries: MX + The numerical timeseries of the system Returns ------- The derivative of muscle activation state """ + a = states[0] + I = controls[0] u = self.normalize_current(I) a_dot = self.get_muscle_activation(a=a, u=u) return vertcat(a_dot) - @staticmethod def dynamics( + self, time: MX, states: MX, controls: MX, @@ -159,64 +186,46 @@ def dynamics( algebraic_states: MX, numerical_timeseries: MX, nlp: NonLinearProgram, - fes_model=None, ) -> DynamicsEvaluation: """ + Functional electrical stimulation dynamic + Parameters ---------- time: MX The system's current node time states: MX - The state of the system (muscle activation) + The state of the system CN, F, A, Tau1, Km controls: MX - The controls of the system (stimulation current) + The controls of the system, none parameters: MX - The parameters acting on the system + The parameters acting on the system, final time of each phase algebraic_states: MX - The stochastic variables of the system + The stochastic variables of the system, none numerical_timeseries: MX The numerical timeseries of the system nlp: NonLinearProgram A reference to the phase - fes_model: VeltinkModelPulseIntensity - The current phase fes model - Returns ------- - The derivative of the states + The derivative of the states in the tuple[MX] format """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - a=states[0], - I=controls[0], - ), - defects=None, + dxdt_fun = nlp.model.system_dynamics + dxdt = dxdt_fun( + time=time, + states=states, + controls=controls, + numerical_timeseries=numerical_timeseries, ) - def declare_model_variables( - self, - ocp: OptimalControlProgram, - nlp: NonLinearProgram, - numerical_data_timeseries: dict[str, np.ndarray] = None, - contact_type: tuple = (), - ): - """ - Tell the program which variables are states and controls. + defects = None + if isinstance(nlp.dynamics_type.ode_solver, OdeSolver.COLLOCATION): + states_dot_list = [] + for key in nlp.model.name_dofs: + states_dot_list.append(DynamicsFunctions.get(nlp.states_dot[key], nlp.states_dot.scaled.cx)) + defects = vertcat(*states_dot_list) - dxdt - Parameters - ---------- - ocp: OptimalControlProgram - A reference to the ocp - nlp: NonLinearProgram - A reference to the phase - numerical_data_timeseries: dict[str, np.ndarray] - A list of values to pass to the dynamics at each node - contact_type: tuple - The type of contact to use for the model - """ - StateConfigure().configure_all_fes_model_states(ocp, nlp, fes_model=self) - StateConfigure().configure_intensity(ocp, nlp) - ConfigureProblem.configure_dynamics_function(ocp, nlp, dyn_func=self.dynamics) + return DynamicsEvaluation( + dxdt=dxdt, + defects=defects, + ) diff --git a/cocofest/models/veltink1992/veltink1992_and_riener1998.py b/cocofest/models/veltink1992/veltink1992_and_riener1998.py index 1ad277df..8b5aaf4d 100644 --- a/cocofest/models/veltink1992/veltink1992_and_riener1998.py +++ b/cocofest/models/veltink1992/veltink1992_and_riener1998.py @@ -2,15 +2,7 @@ from casadi import MX, vertcat import numpy as np -from bioptim import ( - ConfigureProblem, - DynamicsEvaluation, - NonLinearProgram, - OptimalControlProgram, -) - from cocofest.models.veltink1992.veltink1992 import VeltinkModelPulseIntensity -from cocofest.models.state_configure import StateConfigure class VeltinkRienerModelPulseIntensityWithFatigue(VeltinkModelPulseIntensity): @@ -57,8 +49,8 @@ def __init__( self.T_rec = T_rec if T_rec is not None else T_REC_DEFAULT @property - def name_dof(self, with_muscle_name: bool = False) -> list[str]: - muscle_name = "_" + self.muscle_name if self.muscle_name and with_muscle_name else "" + def name_dofs(self, with_muscle_name: bool = False) -> list[str]: + muscle_name = "_" + self.muscle_name if self.muscle_name is not None else "" return [ "a" + muscle_name, # Muscle activation "mu" + muscle_name, # Fatigue state @@ -118,100 +110,34 @@ def get_mu_dot(self, a: MX, mu: MX) -> MX: def system_dynamics( self, - a: MX, - mu: MX, - I: MX, - ) -> MX: - """ - The system dynamics including fatigue effects. - - Parameters - ---------- - a: MX - Muscle activation state (unitless) - mu: MX - Fatigue state (unitless) - I: MX - Stimulation current amplitude (mA) - - Returns - ------- - The value of the derivative of each state - """ - u = self.normalize_current(I) - a_dot = self.get_muscle_activation(a=a, u=u) - mu_dot = self.get_mu_dot(a=a, mu=mu) - - return vertcat(a_dot, mu_dot) - - @staticmethod - def dynamics( time: MX, states: MX, controls: MX, - parameters: MX, - algebraic_states: MX, numerical_timeseries: MX, - nlp: NonLinearProgram, - fes_model=None, - ) -> DynamicsEvaluation: + ) -> MX: """ + The system dynamics including fatigue effects. + Parameters ---------- time: MX The system's current node time states: MX - The state of the system (muscle activation, fatigue) + The state of the system a, mu controls: MX - The controls of the system (stimulation current) - parameters: MX - The parameters acting on the system - algebraic_states: MX - The stochastic variables of the system + The controls of the system, I numerical_timeseries: MX The numerical timeseries of the system - nlp: NonLinearProgram - A reference to the phase - fes_model: VeltinkRienerModelPulseIntensityWithFatigue - The current phase fes model Returns ------- - The derivative of the states - """ - model = fes_model if fes_model else nlp.model - dxdt_fun = model.system_dynamics - - return DynamicsEvaluation( - dxdt=dxdt_fun( - a=states[0], - mu=states[1], - I=controls[0], - ), - defects=None, - ) - - def declare_model_variables( - self, - ocp: OptimalControlProgram, - nlp: NonLinearProgram, - numerical_data_timeseries: dict[str, np.ndarray] = None, - contact_type: tuple = (), - ): + The value of the derivative of each state """ - Tell the program which variables are states and controls. + a = states[0] + mu = states[1] + I = controls[0] + u = self.normalize_current(I) + a_dot = self.get_muscle_activation(a=a, u=u) + mu_dot = self.get_mu_dot(a=a, mu=mu) - Parameters - ---------- - ocp: OptimalControlProgram - A reference to the ocp - nlp: NonLinearProgram - A reference to the phase - numerical_data_timeseries: dict[str, np.ndarray] - A list of values to pass to the dynamics at each node - contact_type: tuple - The type of contact to use for the model - """ - StateConfigure().configure_all_fes_model_states(ocp, nlp, fes_model=self) - StateConfigure().configure_intensity(ocp, nlp) - ConfigureProblem.configure_dynamics_function(ocp, nlp, dyn_func=self.dynamics) + return vertcat(a_dot, mu_dot) diff --git a/cocofest/optimization/fes_id_ocp.py b/cocofest/optimization/fes_id_ocp.py index 8216c8a2..29ceed1c 100644 --- a/cocofest/optimization/fes_id_ocp.py +++ b/cocofest/optimization/fes_id_ocp.py @@ -37,7 +37,7 @@ def set_x_bounds( # Sets the bound for all the phases x_bounds = BoundsList() - variable_bound_list = model.name_dof + variable_bound_list = model.name_dofs starting_bounds, min_bounds, max_bounds = ( model.standard_rest_values(), model.standard_rest_values(), diff --git a/cocofest/optimization/fes_ocp.py b/cocofest/optimization/fes_ocp.py index 0e8b2754..7e2b6e1f 100644 --- a/cocofest/optimization/fes_ocp.py +++ b/cocofest/optimization/fes_ocp.py @@ -3,7 +3,8 @@ from bioptim import ( BoundsList, ConstraintList, - DynamicsList, + DynamicsOptionsList, + DynamicsOptions, InitialGuessList, InterpolationType, Node, @@ -12,7 +13,6 @@ ParameterList, PhaseDynamics, VariableScaling, - OdeSolver, ) from ..fourier_approx import FourierSeries @@ -72,17 +72,18 @@ def set_constraints(model, n_shooting, stim_idx_at_node_list): return constraints @staticmethod - def declare_dynamics(model, numerical_data_timeseries=None, ode_solver=OdeSolver.RK4(n_integration_steps=10)): - dynamics = DynamicsList() - dynamics.add( - model.declare_ding_variables, - dynamic_function=model.dynamics, - expand_dynamics=True, - phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, - numerical_data_timeseries=numerical_data_timeseries, - ode_solver=ode_solver, + def declare_dynamics_options(numerical_time_series, ode_solver): + dynamics_options = DynamicsOptionsList() + dynamics_options.add( + DynamicsOptions( + expand_dynamics=True, + expand_continuity=False, + phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, + ode_solver=ode_solver, + numerical_data_timeseries=numerical_time_series, + ) ) - return dynamics + return dynamics_options @staticmethod def set_x_bounds(model): @@ -99,7 +100,7 @@ def set_x_bounds(model): # Sets the bound for all the phases x_bounds = BoundsList() - variable_bound_list = model.name_dof + variable_bound_list = model.name_dofs starting_bounds, min_bounds, max_bounds = ( model.standard_rest_values(), model.standard_rest_values(), @@ -131,7 +132,7 @@ def set_x_bounds(model): @staticmethod def set_x_init(model): - variable_bound_list = model.name_dof + variable_bound_list = model.name_dofs x_init = InitialGuessList() for j in range(len(variable_bound_list)): x_init.add(variable_bound_list[j], model.standard_rest_values()[j]) diff --git a/cocofest/optimization/fes_ocp_multibody.py b/cocofest/optimization/fes_ocp_multibody.py index b76beefa..e0d79ad9 100644 --- a/cocofest/optimization/fes_ocp_multibody.py +++ b/cocofest/optimization/fes_ocp_multibody.py @@ -3,12 +3,10 @@ from bioptim import ( BoundsList, ConstraintList, - DynamicsList, ExternalForceSetTimeSeries, InitialGuessList, InterpolationType, ParameterList, - PhaseDynamics, VariableScaling, ) @@ -38,22 +36,6 @@ def get_numerical_time_series_for_external_forces(n_shooting, external_force_dic return numerical_time_series, external_force_set - @staticmethod - def declare_dynamics(bio_models, numerical_time_series, ode_solver, contact_type): - dynamics = DynamicsList() - dynamics.add( - bio_models.declare_model_variables, - dynamic_function=bio_models.muscle_dynamic, - expand_dynamics=True, - expand_continuity=False, - phase=0, - phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, - numerical_data_timeseries=numerical_time_series, - contact_type=contact_type, - ode_solver=ode_solver, - ) - return dynamics - @staticmethod def build_parameters( model: FesMskModel, @@ -127,7 +109,7 @@ def set_x_bounds_fes(bio_models): x_init = InitialGuessList() for model in bio_models.muscles_dynamics_model: muscle_name = model.muscle_name - variable_bound_list = [model.name_dof[i] + "_" + muscle_name for i in range(len(model.name_dof))] + variable_bound_list = model.name_dofs starting_bounds, min_bounds, max_bounds = ( model.standard_rest_values(), @@ -286,6 +268,7 @@ def update_model(model, parameters, external_force_set): previous_stim=model.muscles_dynamics_model[0].previous_stim, activate_force_length_relationship=model.activate_force_length_relationship, activate_force_velocity_relationship=model.activate_force_velocity_relationship, + activate_passive_force_relationship=model.activate_passive_force_relationship, activate_residual_torque=model.activate_residual_torque, parameters=parameters, external_force_set=external_force_set, diff --git a/examples/fes_multibody/cycling/cycling_bayesian_mhe.py b/examples/fes_multibody/cycling/cycling_bayesian_mhe.py index 06b6069f..bfb14520 100644 --- a/examples/fes_multibody/cycling/cycling_bayesian_mhe.py +++ b/examples/fes_multibody/cycling/cycling_bayesian_mhe.py @@ -117,7 +117,10 @@ def prepare_nmpc_bo( numerical_time_series.update(time_series2) # --- Dynamics & states --- # - dynamics = base.set_dynamics(model=model, numerical_time_series=numerical_time_series, ode_solver=ode_solver) + dynamics_options = base.set_dynamics_options( + numerical_time_series=numerical_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) + ) + x_init = base.set_q_qdot_init( n_shooting=window_n_shooting, pedal_config=cycling_info["pedal_config"], @@ -148,7 +151,7 @@ def prepare_nmpc_bo( return base.MyCyclicNMPC( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, cycle_len=cycle_len, cycle_duration=cycle_duration, n_cycles_simultaneous=n_cycles_simultaneous, diff --git a/examples/fes_multibody/cycling/cycling_pulse_width_mhe.py b/examples/fes_multibody/cycling/cycling_pulse_width_mhe.py index d5fe7c34..d5fbb1b2 100644 --- a/examples/fes_multibody/cycling/cycling_pulse_width_mhe.py +++ b/examples/fes_multibody/cycling/cycling_pulse_width_mhe.py @@ -18,9 +18,7 @@ BoundsList, ConstraintList, ConstraintFcn, - ContactType, CostType, - DynamicsList, ExternalForceSetTimeSeries, InitialGuessList, InterpolationType, @@ -29,7 +27,6 @@ ObjectiveFcn, ObjectiveList, OdeSolver, - PhaseDynamics, SolutionMerge, Solution, Solver, @@ -333,7 +330,7 @@ def prepare_nmpc( ) numerical_time_series.update(numerical_data_time_series) # --- Dynamics --- # - dynamics = set_dynamics(model=model, numerical_time_series=numerical_time_series, ode_solver=ode_solver) + dynamics_options = set_dynamics_options(numerical_time_series=numerical_time_series, ode_solver=ode_solver) # --- Set states --- # # --- Set q (position and speed) initial guesses --- # @@ -377,7 +374,7 @@ def prepare_nmpc( return MyCyclicNMPC( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, cycle_len=cycle_len, cycle_duration=cycle_duration, n_cycles_simultaneous=n_cycles_simultaneous, @@ -406,19 +403,11 @@ def set_external_forces(n_shooting, external_force_dict, force_name): return numerical_time_series, external_force_set -def set_dynamics(model, numerical_time_series, ode_solver): - dynamics = DynamicsList() - dynamics.add( - dynamics_type=model.declare_model_variables, - dynamic_function=model.muscle_dynamic, - expand_dynamics=True, - phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, - numerical_data_timeseries=numerical_time_series, - contact_type=[ContactType.RIGID_EXPLICIT], - phase=0, - ode_solver=ode_solver, +def set_dynamics_options(numerical_time_series, ode_solver): + dynamics_options = OcpFesMsk.declare_dynamics_options( + numerical_time_series=numerical_time_series, ode_solver=ode_solver ) - return dynamics + return dynamics_options def set_q_qdot_init( @@ -665,6 +654,7 @@ def updating_model(model: FesMskModel, external_force_set, parameters=None) -> F activate_residual_torque=model.activate_residual_torque, parameters=parameters, external_force_set=external_force_set, + with_contact=True, ) return model @@ -728,6 +718,7 @@ def set_fes_model(model_path, stim_time): activate_passive_force_relationship=True, activate_residual_torque=False, external_force_set=None, # External forces will be added later (resistive_torque) + with_contact=True, ) return fes_model @@ -1023,9 +1014,9 @@ def main( main( stimulation_frequency=30, - n_total_cycle=1, + n_total_cycle=5, n_cycles_simultaneous=[2], # [2, 3, 4, 5] - resistive_torque=-0.3, # (N.m) + resistive_torque=-0.2, # (N.m) cost_fun_weight=[(1, 0, 0)], # (min_force, min_fatigue, min_control) init_guess=False, save=False, diff --git a/examples/fes_multibody/cycling/cycling_with_different_driven_methods.py b/examples/fes_multibody/cycling/cycling_with_different_driven_methods.py index 016099f9..b6c00403 100644 --- a/examples/fes_multibody/cycling/cycling_with_different_driven_methods.py +++ b/examples/fes_multibody/cycling/cycling_with_different_driven_methods.py @@ -14,8 +14,6 @@ ConstraintFcn, ControlType, CostType, - DynamicsFcn, - DynamicsList, ExternalForceSetTimeSeries, InitialGuessList, InterpolationType, @@ -28,7 +26,8 @@ PhaseDynamics, Solver, VariableScalingList, - ContactType, + DynamicsOptionsList, + DynamicsOptions, ) from cocofest import ( @@ -92,9 +91,11 @@ def update_model( previous_stim=model.muscles_dynamics_model[0].previous_stim, activate_force_length_relationship=model.activate_force_length_relationship, activate_force_velocity_relationship=model.activate_force_velocity_relationship, + activate_passive_force_relationship=model.activate_force_velocity_relationship, activate_residual_torque=model.activate_residual_torque, parameters=parameters, external_force_set=external_force_set, + with_contact=model.with_contact, ) else: model = BiorbdModel(model.path, external_force_set=external_force_set) @@ -102,55 +103,55 @@ def update_model( return model -def set_dynamics( - model: BiorbdModel | FesMskModel, - numerical_time_series: dict, - dynamics_type_str: str = "torque_driven", - ode_solver: OdeSolver = OdeSolver.RK4(n_integration_steps=10), -) -> DynamicsList: - """ - Set the dynamics of the optimal control program based on the chosen dynamics type. - - Parameters - ---------- - model: BiorbdModel | FesMskModel - The biomechanical model. - numerical_time_series: dict - External numerical data (e.g., external forces). - dynamics_type_str: str - Type of dynamics ("torque_driven", "muscle_driven", or "fes_driven"). - - Returns - ------- - A DynamicsList configured for the problem. - """ - dynamics_type = ( - DynamicsFcn.TORQUE_DRIVEN - if dynamics_type_str == "torque_driven" - else ( - DynamicsFcn.MUSCLE_DRIVEN - if dynamics_type_str == "muscle_driven" - else model.declare_model_variables if dynamics_type_str == "fes_driven" else None - ) - ) - if dynamics_type is None: - raise ValueError("Dynamics type not recognized") - - dynamics = DynamicsList() - dynamics.add( - dynamics_type=dynamics_type, - dynamic_function=( - None if dynamics_type in (DynamicsFcn.TORQUE_DRIVEN, DynamicsFcn.MUSCLE_DRIVEN) else model.muscle_dynamic - ), - expand_dynamics=True, - expand_continuity=False, - phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, - numerical_data_timeseries=numerical_time_series, - phase=0, - ode_solver=ode_solver, - contact_type=[ContactType.RIGID_EXPLICIT], - ) - return dynamics +# def set_dynamics( +# model: BiorbdModel | FesMskModel, +# numerical_time_series: dict, +# dynamics_type_str: str = "torque_driven", +# ode_solver: OdeSolver = OdeSolver.RK4(n_integration_steps=10), +# ) -> DynamicsList: +# """ +# Set the dynamics of the optimal control program based on the chosen dynamics type. +# +# Parameters +# ---------- +# model: BiorbdModel | FesMskModel +# The biomechanical model. +# numerical_time_series: dict +# External numerical data (e.g., external forces). +# dynamics_type_str: str +# Type of dynamics ("torque_driven", "muscle_driven", or "fes_driven"). +# +# Returns +# ------- +# A DynamicsList configured for the problem. +# """ +# dynamics_type = ( +# DynamicsFcn.TORQUE_DRIVEN +# if dynamics_type_str == "torque_driven" +# else ( +# DynamicsFcn.MUSCLE_DRIVEN +# if dynamics_type_str == "muscle_driven" +# else model.declare_model_variables if dynamics_type_str == "fes_driven" else None +# ) +# ) +# if dynamics_type is None: +# raise ValueError("Dynamics type not recognized") +# +# dynamics = DynamicsList() +# dynamics.add( +# dynamics_type=dynamics_type, +# dynamic_function=( +# None if dynamics_type in (DynamicsFcn.TORQUE_DRIVEN, DynamicsFcn.MUSCLE_DRIVEN) else model.muscle_dynamic +# ), +# expand_dynamics=True, +# expand_continuity=False, +# phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, +# numerical_data_timeseries=numerical_time_series, +# phase=0, +# ode_solver=ode_solver, +# contact_type=[ContactType.RIGID_EXPLICIT], +# ) +# return dynamics def set_objective_functions(model: BiorbdModel | FesMskModel, dynamics_type: str, init_x) -> ObjectiveList: @@ -452,11 +453,21 @@ def prepare_ocp( numerical_time_series.update(numerical_data_time_series) # Set dynamics based on the chosen dynamics type - dynamics = set_dynamics( - model, - numerical_time_series, - dynamics_type_str=dynamics_type, - ode_solver=ode_solver, + # dynamics = set_dynamics( + # model, + # numerical_time_series, + # dynamics_type_str=dynamics_type, + # ode_solver=ode_solver, + # ) + dynamics_options = DynamicsOptionsList() + dynamics_options.add( + DynamicsOptions( + expand_dynamics=True, + expand_continuity=False, + phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, + ode_solver=ode_solver, + numerical_data_timeseries=numerical_time_series, + ) ) # Set initial guess for state variables @@ -486,9 +497,9 @@ def prepare_ocp( return OptimalControlProgram( [model], - dynamics, n_shooting, final_time, + dynamics_options, x_bounds=x_bounds, u_bounds=u_bounds, x_init=x_init, @@ -563,6 +574,7 @@ def main( activate_passive_force_relationship=True, activate_residual_torque=False, external_force_set=None, # External forces will be added later + with_contact=True, ) # Adjust n_shooting based on the stimulation time n_shooting = model.muscles_dynamics_model[0].get_n_shooting(final_time) @@ -578,7 +590,7 @@ def main( dynamics_type=dynamics_type, use_sx=False, ode_solver=OdeSolver.COLLOCATION(polynomial_degree=3, method="radau"), - # ode_solver=OdeSolver.RK4(n_integration_steps=5) + # ode_solver=OdeSolver.RK4(n_integration_steps=5), torque=-0.3, initial_guess_model_path=initial_guess_model_path, ) diff --git a/examples/fes_multibody/elbow_flexion/elbow_flexion_task.py b/examples/fes_multibody/elbow_flexion/elbow_flexion_task.py index 18674678..0ba091e7 100644 --- a/examples/fes_multibody/elbow_flexion/elbow_flexion_task.py +++ b/examples/fes_multibody/elbow_flexion/elbow_flexion_task.py @@ -42,12 +42,8 @@ def prepare_ocp( numerical_data_time_series, stim_idx_at_node_list = muscle_model.get_numerical_data_time_series( n_shooting, final_time ) - - dynamics = OcpFesMsk.declare_dynamics( - model, - numerical_time_series=numerical_data_time_series, - ode_solver=OdeSolver.RK4(n_integration_steps=10), - contact_type=[], + dynamics_options = OcpFesMsk.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) x_bounds, x_init = OcpFesMsk.set_x_bounds(model, msk_info) @@ -99,7 +95,7 @@ def prepare_ocp( return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, @@ -139,6 +135,7 @@ def main(plot=True, model_path="../../msk_models/Arm26/arm26_biceps_triceps.bioM activate_passive_force_relationship=True, activate_residual_torque=False, stim_time=stim_time, + with_contact=False, ) final_time = 1 diff --git a/examples/fes_multibody/elbow_flexion/frequency_optimization_multibody.py b/examples/fes_multibody/elbow_flexion/frequency_optimization_multibody.py index 2535a648..d5492244 100644 --- a/examples/fes_multibody/elbow_flexion/frequency_optimization_multibody.py +++ b/examples/fes_multibody/elbow_flexion/frequency_optimization_multibody.py @@ -21,12 +21,8 @@ def prepare_ocp(model, final_time: float, resistive_torque, msk_info): n_shooting, final_time ) numerical_time_series.update(numerical_data_time_series) - - dynamics = OcpFesMsk.declare_dynamics( - model, - numerical_time_series=numerical_time_series, - contact_type=[], - ode_solver=OdeSolver.RK4(n_integration_steps=10), + dynamics_options = OcpFesMsk.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) # --- Set initial guesses and bounds for states and controls --- # @@ -48,11 +44,12 @@ def prepare_ocp(model, final_time: float, resistive_torque, msk_info): activate_force_velocity_relationship=model.activate_force_velocity_relationship, activate_residual_torque=model.activate_residual_torque, external_force_set=external_force_set, + with_contact=False, ) return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, diff --git a/examples/fes_multibody/reaching/reaching_task.py b/examples/fes_multibody/reaching/reaching_task.py index 3a131ba0..1ffe5a16 100644 --- a/examples/fes_multibody/reaching/reaching_task.py +++ b/examples/fes_multibody/reaching/reaching_task.py @@ -52,6 +52,7 @@ def initialize_model(final_time, model_path): activate_force_velocity_relationship=True, activate_passive_force_relationship=True, activate_residual_torque=False, + with_contact=False, ) return model @@ -66,12 +67,9 @@ def prepare_ocp( numerical_data_time_series, stim_idx_at_node_list = muscle_model.get_numerical_data_time_series( n_shooting, final_time ) - - dynamics = OcpFesMsk.declare_dynamics( - model, + dynamics_options = OcpFesMsk.declare_dynamics_options( numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.COLLOCATION(polynomial_degree=3, method="radau"), - contact_type=[], ) # --- Initialize default FES bounds and initial guess --- # @@ -124,7 +122,7 @@ def prepare_ocp( return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, @@ -152,11 +150,3 @@ def main(plot=True, model_path="../../msk_models/Arm26/arm26_with_reaching_targe if __name__ == "__main__": main() - -# [1] Dahmane, R., Djordjevič, S., Šimunič, B., & Valenčič, V. (2005). -# Spatial fiber type distribution in normal human muscle: histochemical and tensiomyographical evaluation. -# Journal of biomechanics, 38(12), 2451-2459. - -# [2] Klein, C. S., Allman, B. L., Marsh, G. D., & Rice, C. L. (2002). -# Muscle size, strength, and bone geometry in the upper limbs of young and old men. -# The Journals of Gerontology Series A: Biological Sciences and Medical Sciences, 57(7), M455-M459. diff --git a/examples/getting_started/identification/muscle_model_id.py b/examples/getting_started/identification/muscle_model_id.py index d3910f7e..889ccace 100644 --- a/examples/getting_started/identification/muscle_model_id.py +++ b/examples/getting_started/identification/muscle_model_id.py @@ -7,6 +7,7 @@ DingModelPulseIntensityFrequency, IvpFes, OcpFesId, + OcpFes, ) from cocofest.identification.identification_method import DataExtraction @@ -71,10 +72,8 @@ def prepare_ocp( ) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFesId.declare_dynamics( - model=model, - numerical_data_timeseries=numerical_data_time_series, - ode_solver=OdeSolver.RK4(n_integration_steps=10), + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) x_bounds, x_init = OcpFesId.set_x_bounds( @@ -107,7 +106,7 @@ def prepare_ocp( return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, x_init=x_init, diff --git a/examples/getting_started/multibody/pulse_intensity_optimization_multibody.py b/examples/getting_started/multibody/pulse_intensity_optimization_multibody.py index a429804f..a4dff8f7 100644 --- a/examples/getting_started/multibody/pulse_intensity_optimization_multibody.py +++ b/examples/getting_started/multibody/pulse_intensity_optimization_multibody.py @@ -34,11 +34,9 @@ def prepare_ocp(model: FesMskModel, final_time: float, external_force: dict, msk else: numerical_time_series = numerical_data_time_series - dynamics = OcpFesMsk.declare_dynamics( - model, + dynamics_options = OcpFesMsk.declare_dynamics_options( numerical_time_series=numerical_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10), - contact_type=[], ) x_bounds, x_init = OcpFesMsk.set_x_bounds(model, msk_info) @@ -77,7 +75,7 @@ def prepare_ocp(model: FesMskModel, final_time: float, external_force: dict, msk return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, diff --git a/examples/getting_started/multibody/pulse_width_optimization_multibody.py b/examples/getting_started/multibody/pulse_width_optimization_multibody.py index 10dc8110..cd67a5d8 100644 --- a/examples/getting_started/multibody/pulse_width_optimization_multibody.py +++ b/examples/getting_started/multibody/pulse_width_optimization_multibody.py @@ -35,11 +35,9 @@ def prepare_ocp(model: FesMskModel, final_time: float, external_force: dict, msk else: numerical_time_series = numerical_data_time_series - dynamics = OcpFesMsk.declare_dynamics( - model, + dynamics_options = OcpFesMsk.declare_dynamics_options( numerical_time_series=numerical_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10), - contact_type=[], ) x_bounds, x_init = OcpFesMsk.set_x_bounds(model, msk_info) @@ -57,7 +55,7 @@ def prepare_ocp(model: FesMskModel, final_time: float, external_force: dict, msk return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, @@ -83,6 +81,7 @@ def main(plot=True, biorbd_path="../../msk_models/Arm26/arm26_biceps_1dof.bioMod activate_passive_force_relationship=True, activate_residual_torque=False, external_force_set=None, # External forces will be added later + with_contact=False, ) resistive_torque = { diff --git a/examples/getting_started/optimization/force_tracking.py b/examples/getting_started/optimization/force_tracking.py index 77de4c06..2c81c9ea 100644 --- a/examples/getting_started/optimization/force_tracking.py +++ b/examples/getting_started/optimization/force_tracking.py @@ -35,13 +35,11 @@ def prepare_ocp(model: FesModel, final_time: float, pw_max: float, force_trackin time_series, stim_idx_at_node_list = model.get_numerical_data_time_series( n_shooting, final_time ) # Retrieve time and indexes at which occurs the stimulation for the FES dynamic - dynamics = OcpFes.declare_dynamics( - model, - time_series, + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10), - # ode_solver=OdeSolver.COLLOCATION(polynomial_degree=3, method="radau"), # Possibility to use a different solver + # ode_solver=OdeSolver.COLLOCATION(polynomial_degree=3, method="radau"), # Possibility to use a different solver) ) - # --- Set initial guesses and bounds for states and controls --- # x_bounds = OcpFes.set_x_bounds(model) x_init = OcpFes.set_x_init(model) @@ -61,7 +59,7 @@ def prepare_ocp(model: FesModel, final_time: float, pw_max: float, force_trackin return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, diff --git a/examples/getting_started/optimization/frequency_optimization.py b/examples/getting_started/optimization/frequency_optimization.py index dc5f4da5..233c71cd 100644 --- a/examples/getting_started/optimization/frequency_optimization.py +++ b/examples/getting_started/optimization/frequency_optimization.py @@ -4,7 +4,14 @@ No optimization will be done on the stimulation, the frequency is fixed to 1Hz. """ -from bioptim import ObjectiveList, ObjectiveFcn, OptimalControlProgram, ControlType, OdeSolver, Node +from bioptim import ( + ObjectiveList, + ObjectiveFcn, + OptimalControlProgram, + ControlType, + OdeSolver, + Node, +) from cocofest import OcpFes, ModelMaker @@ -12,8 +19,9 @@ def prepare_ocp(model, final_time): # --- Set dynamics --- # n_shooting = model.get_n_shooting(final_time=final_time) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFes.declare_dynamics( - model, numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) + + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) # --- Set initial guesses and bounds for states --- # @@ -28,7 +36,7 @@ def prepare_ocp(model, final_time): return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, diff --git a/examples/getting_started/optimization/pulse_intensity_optimization.py b/examples/getting_started/optimization/pulse_intensity_optimization.py index 40307edc..344e1c68 100644 --- a/examples/getting_started/optimization/pulse_intensity_optimization.py +++ b/examples/getting_started/optimization/pulse_intensity_optimization.py @@ -17,8 +17,8 @@ def prepare_ocp(model, final_time, pi_max): # --- Set dynamics --- # n_shooting = model.get_n_shooting(final_time=final_time) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFes.declare_dynamics( - model, numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) # --- Set initial guesses and bounds for states and controls --- # @@ -51,7 +51,7 @@ def prepare_ocp(model, final_time, pi_max): return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, diff --git a/examples/getting_started/optimization/pulse_width_optimization.py b/examples/getting_started/optimization/pulse_width_optimization.py index 0cc00f22..1b0ceed1 100644 --- a/examples/getting_started/optimization/pulse_width_optimization.py +++ b/examples/getting_started/optimization/pulse_width_optimization.py @@ -13,8 +13,8 @@ def prepare_ocp(model, final_time, pw_max): # --- Set dynamics --- # n_shooting = model.get_n_shooting(final_time=final_time) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFes.declare_dynamics( - model, numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) # --- Set initial guesses and bounds for states and controls --- # @@ -32,7 +32,7 @@ def prepare_ocp(model, final_time, pw_max): return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, diff --git a/examples/getting_started/optimization/pulse_width_optimization_nmpc.py b/examples/getting_started/optimization/pulse_width_optimization_nmpc.py index ffc3911c..a92f343d 100644 --- a/examples/getting_started/optimization/pulse_width_optimization_nmpc.py +++ b/examples/getting_started/optimization/pulse_width_optimization_nmpc.py @@ -50,7 +50,10 @@ def prepare_nmpc( total_cycle_len, total_cycle_duration ) - dynamics = OcpFes.declare_dynamics(model, numerical_data_time_series, ode_solver) + # dynamics = OcpFes.declare_dynamics(model, numerical_data_time_series, ode_solver) + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=ode_solver + ) x_bounds = OcpFes.set_x_bounds(model) x_init_fes = OcpFes.set_x_init(model) @@ -94,7 +97,7 @@ def prepare_nmpc( return FesNmpc( bio_model=model, - dynamics=dynamics, + dynamics=dynamics_options, cycle_len=cycle_len, cycle_duration=cycle_duration, n_cycles_simultaneous=n_cycles_simultaneous, @@ -138,7 +141,8 @@ def main(plot=True): minimize_force=True, minimize_fatigue=False, use_sx=False, - ode_solver=OdeSolver.COLLOCATION(polynomial_degree=5, method="radau"), + ode_solver=OdeSolver.RK4(n_integration_steps=10), + # ode_solver=OdeSolver.COLLOCATION(polynomial_degree=5, method="radau"), ) def update_functions(_nmpc: MultiCyclicNonlinearModelPredictiveControl, cycle_idx: int, _sol: Solution): diff --git a/examples/identification/force_model/ding2003_model_id.py b/examples/identification/force_model/ding2003_model_id.py index 285ce005..bdee3084 100644 --- a/examples/identification/force_model/ding2003_model_id.py +++ b/examples/identification/force_model/ding2003_model_id.py @@ -26,6 +26,7 @@ OcpFesId, ModelMaker, FES_plot, + OcpFes, ) from cocofest.identification.identification_method import DataExtraction @@ -66,10 +67,8 @@ def prepare_ocp( ) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFesId.declare_dynamics( - model=model, - numerical_data_timeseries=numerical_data_time_series, - ode_solver=OdeSolver.RK4(n_integration_steps=10), + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) x_bounds, x_init = OcpFesId.set_x_bounds( @@ -97,7 +96,7 @@ def prepare_ocp( return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, x_init=x_init, diff --git a/examples/identification/force_model/ding2007_model_id.py b/examples/identification/force_model/ding2007_model_id.py index d5144c91..f2c6acda 100644 --- a/examples/identification/force_model/ding2007_model_id.py +++ b/examples/identification/force_model/ding2007_model_id.py @@ -4,10 +4,9 @@ Finally, we use the data to identify the model parameters. """ -import matplotlib.pyplot as plt import numpy as np -from bioptim import SolutionMerge, OdeSolver, OptimalControlProgram, ObjectiveFcn, Node, ControlType, ObjectiveList +from bioptim import OdeSolver, OptimalControlProgram, ObjectiveFcn, Node, ControlType, ObjectiveList from cocofest import ( DingModelPulseWidthFrequency, @@ -15,6 +14,7 @@ ModelMaker, OcpFesId, FES_plot, + OcpFes, ) from cocofest.identification.identification_method import DataExtraction @@ -56,10 +56,8 @@ def prepare_ocp( ) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFesId.declare_dynamics( - model=model, - numerical_data_timeseries=numerical_data_time_series, - ode_solver=OdeSolver.RK4(n_integration_steps=10), + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) x_bounds, x_init = OcpFesId.set_x_bounds( @@ -92,7 +90,7 @@ def prepare_ocp( return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, x_init=x_init, diff --git a/examples/identification/force_model/hmed2018_model_id.py b/examples/identification/force_model/hmed2018_model_id.py index 3e247556..0fd90ca8 100644 --- a/examples/identification/force_model/hmed2018_model_id.py +++ b/examples/identification/force_model/hmed2018_model_id.py @@ -1,6 +1,5 @@ import numpy as np -import matplotlib.pyplot as plt -from bioptim import OdeSolver, ObjectiveList, ObjectiveFcn, Node, OptimalControlProgram, ControlType, SolutionMerge +from bioptim import OdeSolver, ObjectiveList, ObjectiveFcn, Node, OptimalControlProgram, ControlType from cocofest import ( ModelMaker, @@ -8,6 +7,7 @@ IvpFes, OcpFesId, FES_plot, + OcpFes, ) from cocofest.identification.identification_method import DataExtraction @@ -50,10 +50,8 @@ def prepare_ocp( ) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFesId.declare_dynamics( - model=model, - numerical_data_timeseries=numerical_data_time_series, - ode_solver=OdeSolver.RK4(n_integration_steps=10), + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) x_bounds, x_init = OcpFesId.set_x_bounds( @@ -86,7 +84,7 @@ def prepare_ocp( return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, x_init=x_init, @@ -138,7 +136,7 @@ def main(plot=True): default_model = DingModelPulseIntensityFrequency() FES_plot(data=sol).plot( - title="Identification of HMed 2018 parameters", + title="Identification of Hmed 2018 parameters", tracked_data=sim_data, default_model=default_model, show_bounds=False, diff --git a/examples/other_fes_models/marion2009_example.py b/examples/other_fes_models/marion2009_example.py index 53d5ebeb..a0fd10d5 100644 --- a/examples/other_fes_models/marion2009_example.py +++ b/examples/other_fes_models/marion2009_example.py @@ -9,6 +9,7 @@ Solver, ObjectiveList, ObjectiveFcn, + OdeSolver, OptimalControlProgram, ControlType, Node, @@ -23,7 +24,9 @@ def prepare_ocp(model, final_time, pw_max=0.0006): # --- Set dynamics --- # n_shooting = model.get_n_shooting(final_time=final_time) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFes.declare_dynamics(model, numerical_data_time_series) + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) + ) # --- Set initial guesses and bounds for states and controls --- # x_bounds = OcpFes.set_x_bounds(model) @@ -44,7 +47,7 @@ def prepare_ocp(model, final_time, pw_max=0.0006): return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, diff --git a/examples/other_fes_models/marion2013_example.py b/examples/other_fes_models/marion2013_example.py index f5afc01c..59917adf 100644 --- a/examples/other_fes_models/marion2013_example.py +++ b/examples/other_fes_models/marion2013_example.py @@ -11,6 +11,7 @@ Solver, ObjectiveList, ObjectiveFcn, + OdeSolver, OptimalControlProgram, ControlType, Node, @@ -26,7 +27,9 @@ def prepare_ocp(model, final_time, pw_max=0.0006): # --- Set dynamics --- # n_shooting = model.get_n_shooting(final_time=final_time) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFes.declare_dynamics(model, numerical_data_time_series) + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) + ) # --- Set initial guesses and bounds for states and controls --- # x_bounds = OcpFes.set_x_bounds(model) @@ -59,7 +62,7 @@ def prepare_ocp(model, final_time, pw_max=0.0006): return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, diff --git a/examples/other_fes_models/veltink1992_example.py b/examples/other_fes_models/veltink1992_example.py index aed6a0cb..1f7de026 100644 --- a/examples/other_fes_models/veltink1992_example.py +++ b/examples/other_fes_models/veltink1992_example.py @@ -1,12 +1,8 @@ -from cocofest import ( - VeltinkModelPulseIntensity, - VeltinkRienerModelPulseIntensityWithFatigue, -) +from cocofest import VeltinkModelPulseIntensity, VeltinkRienerModelPulseIntensityWithFatigue, OcpFes from bioptim import ( OptimalControlProgram, - PhaseDynamics, - DynamicsList, + OdeSolver, Solver, ObjectiveList, ObjectiveFcn, @@ -17,13 +13,8 @@ def prepare_ocp(model, final_time, n_shooting, fmax): # --- Set dynamics --- # - dynamics = DynamicsList() - dynamics.add( - model.declare_model_variables, - dynamic_function=model.dynamics, - expand_dynamics=True, - phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, - numerical_data_timeseries=None, + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=None, ode_solver=OdeSolver.RK4(n_integration_steps=10) ) # --- Set initial guesses and bounds for states and controls --- # @@ -42,7 +33,7 @@ def prepare_ocp(model, final_time, n_shooting, fmax): return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, objective_functions=objective_functions, diff --git a/examples/sensitivity/force_length_velocity/muscle_force_length_and_force_velocity_relationships_comparison.py b/examples/sensitivity/force_length_velocity/muscle_force_length_and_force_velocity_relationships_comparison.py deleted file mode 100644 index dbcf8940..00000000 --- a/examples/sensitivity/force_length_velocity/muscle_force_length_and_force_velocity_relationships_comparison.py +++ /dev/null @@ -1,65 +0,0 @@ -""" -This example is used to compare the effect of the muscle force-length and force-velocity relationships -on the joint angle. -""" - -import numpy as np - -import matplotlib.pyplot as plt - -from bioptim import SolutionMerge - -from cocofest import DingModelPulseWidthFrequencyWithFatigue, OcpFesMsk, FesMskModel - -minimum_pulse_width = DingModelPulseWidthFrequencyWithFatigue().pd0 - -sol_list = [] -sol_time = [] -activate_force_length_relationship = [False, True] - -for i in range(2): - - model = FesMskModel( - name=None, - biorbd_path="../../msk_models/Arm26/arm26_biceps_1dof.bioMod", - muscles_model=[DingModelPulseWidthFrequencyWithFatigue(muscle_name="BIClong")], - activate_force_length_relationship=activate_force_length_relationship[i], - activate_force_velocity_relationship=activate_force_length_relationship[i], - activate_residual_torque=False, - ) - - ocp = OcpFesMsk.prepare_ocp( - model=model, - stim_time=np.linspace(0, 1, 11)[:-1], - final_time=1, - pulse_width={"fixed": 0.00025}, - msk_info={ - "bound_type": "start", - "bound_data": [0], - "with_residual_torque": False, - }, - use_sx=False, - ) - sol = ocp.solve() - sol_list.append(sol.decision_states(to_merge=[SolutionMerge.NODES, SolutionMerge.PHASES])) - time = np.concatenate( - sol.stepwise_time(to_merge=[SolutionMerge.NODES, SolutionMerge.PHASES], duplicated_times=False), - axis=0, - ) - index = 0 - for j in range(len(sol.ocp.nlp) - 1): - index = index + 1 + sol.ocp.nlp[j].ns - time = np.insert(time, index, time[index - 1]) - - sol_time.append(time) - -plt.plot(sol_time[0], np.degrees(sol_list[0]["q"][0]), label="without relationships") -plt.plot(sol_time[1], np.degrees(sol_list[1]["q"][0]), label="with relationships") - -plt.xlabel("Time (s)") -plt.ylabel("Angle (°)") -plt.legend() -plt.show() - -joint_overestimation = np.degrees(sol_list[0]["q"][0][-1]) - np.degrees(sol_list[1]["q"][0][-1]) -print(f"Joint overestimation: {joint_overestimation} degrees") diff --git a/examples/sensitivity/force_length_velocity/muscle_relationships_comparison.py b/examples/sensitivity/force_length_velocity/muscle_relationships_comparison.py new file mode 100644 index 00000000..0663a7e9 --- /dev/null +++ b/examples/sensitivity/force_length_velocity/muscle_relationships_comparison.py @@ -0,0 +1,131 @@ +""" +This example is used to compare the effect of the muscle force-length, force-velocity, and passive-force relationships. +""" + +import numpy as np + +import matplotlib.pyplot as plt + +from bioptim import ( + SolutionMerge, + OdeSolver, + ObjectiveList, + ObjectiveFcn, + ParameterList, + OptimalControlProgram, + ControlType, + ConstraintFcn, + ConstraintList, + Node, +) + +from cocofest import DingModelPulseWidthFrequencyWithFatigue, OcpFesMsk, FesMskModel, CustomObjective, OcpFes + + +def prepare_ocp(model: FesMskModel, final_time: float, msk_info: dict, fixed_pw): + + muscle_model = model.muscles_dynamics_model[0] + n_shooting = muscle_model.get_n_shooting(final_time) + numerical_data_time_series, stim_idx_at_node_list = muscle_model.get_numerical_data_time_series( + n_shooting, final_time + ) + dynamics_options = OcpFes.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10) + ) + + x_bounds, x_init = OcpFesMsk.set_x_bounds(model, msk_info) + u_bounds, u_init = OcpFesMsk.set_u_bounds(model, msk_info["with_residual_torque"], max_bound=0.0006) + + objective_functions = ObjectiveList() + objective_functions.add( + CustomObjective.minimize_overall_muscle_force_production, + custom_type=ObjectiveFcn.Lagrange, + weight=1, + quadratic=True, + ) + + model = OcpFesMsk.update_model(model, parameters=ParameterList(use_sx=True), external_force_set=None) + + constraints = ConstraintList() + constraints.add( + ConstraintFcn.BOUND_CONTROL, + key="last_pulse_width_BIClong", + node=Node.ALL_SHOOTING, + min_bound=np.array([fixed_pw]), + max_bound=np.array([fixed_pw]), + ) + + return OptimalControlProgram( + bio_model=[model], + dynamics=dynamics_options, + n_shooting=n_shooting, + phase_time=final_time, + objective_functions=objective_functions, + x_init=x_init, + x_bounds=x_bounds, + u_init=u_init, + u_bounds=u_bounds, + constraints=constraints, + control_type=ControlType.CONSTANT, + use_sx=True, + n_threads=20, + ) + + +# --- Main --- # +minimum_pulse_width = DingModelPulseWidthFrequencyWithFatigue().pd0 + +sol_list = [] +sol_time = [] + +relationship_activation_dict = { + "no_relationship": [False, False, False], + "length": [True, False, False], + "velocity": [False, True, False], + "passive_force": [False, False, True], + "length_velocity": [True, True, False], + "length_passive_force": [True, False, True], + "velocity_passive_force": [False, True, True], + "all_relationship": [True, True, True], +} + +keys = list(relationship_activation_dict.keys()) +for key in keys: + stim_time = np.linspace(0, 1, 11)[:-1] + model = FesMskModel( + name=None, + biorbd_path="../../msk_models/Arm26/arm26_biceps_1dof.bioMod", + muscles_model=[DingModelPulseWidthFrequencyWithFatigue(muscle_name="BIClong")], + stim_time=list(stim_time), + activate_force_length_relationship=relationship_activation_dict[key][0], + activate_force_velocity_relationship=relationship_activation_dict[key][1], + activate_passive_force_relationship=relationship_activation_dict[key][2], + activate_residual_torque=False, + ) + + ocp = prepare_ocp( + model=model, + final_time=1, + msk_info={ + "bound_type": "start", + "bound_data": [0], + "with_residual_torque": False, + }, + fixed_pw=0.00025, + ) + + sol = ocp.solve() + sol_list.append(sol.stepwise_states(to_merge=[SolutionMerge.NODES, SolutionMerge.PHASES])) + sol_time.append(sol.stepwise_time(to_merge=[SolutionMerge.NODES, SolutionMerge.PHASES])) + +for i in range(len(sol_time)): + key = keys[i] + plt.plot(sol_time[i], np.degrees(sol_list[i]["q"][0]), label=key) + + joint_error = np.degrees(sol_list[-1]["q"][0][-1]) - np.degrees(sol_list[i]["q"][0][-1]) + print(f"Joint error for {key}: {joint_error} degrees") + +plt.xlabel("Time (s)") +plt.ylabel("Angle (°)") +plt.legend() +plt.show() diff --git a/external/bioptim b/external/bioptim index d0d8d4b7..e23300dd 160000 --- a/external/bioptim +++ b/external/bioptim @@ -1 +1 @@ -Subproject commit d0d8d4b7fb8e0fc67584ec53c1ae8dfeac02bba7 +Subproject commit e23300dd14438922f63fa6c3b3d44747083af79b diff --git a/external/bioptim_install_windows.sh b/external/bioptim_install_windows.sh index c6ae70cc..07dfeb05 100644 --- a/external/bioptim_install_windows.sh +++ b/external/bioptim_install_windows.sh @@ -9,4 +9,4 @@ git submodule update --recursive --init cd bioptim # Installing bioptim -python setup.py install +pip install . diff --git a/tests/shard1/test_models_dynamics_without_bioptim.py b/tests/shard1/test_models_dynamics_without_bioptim.py index 4de8d2f8..e85287e4 100644 --- a/tests/shard1/test_models_dynamics_without_bioptim.py +++ b/tests/shard1/test_models_dynamics_without_bioptim.py @@ -9,7 +9,7 @@ def test_ding2003_dynamics(): model = ModelMaker.create_model("ding2003_with_fatigue") assert model.nb_state == 5 - assert model.name_dof == ["Cn", "F", "A", "Tau1", "Km"] + assert model.name_dofs == ["Cn", "F", "A", "Tau1", "Km"] np.testing.assert_almost_equal( model.standard_rest_values(), np.array([[0], [0], [model.a_rest], [model.tau1_rest], [model.km_rest]]), @@ -47,13 +47,10 @@ def test_ding2003_dynamics(): np.testing.assert_almost_equal( np.array( model.system_dynamics( - cn=5, - f=100, - a=3009, - tau1=0.050957, - km=0.103, - t=0.11, - t_stim_prev=np.array([0, 0.1]), + time=0.11, + states=[5, 100, 3009, 0.050957, 0.103], + controls=None, + numerical_timeseries=np.array([0, 0.1]), ) ).squeeze(), np.array(DM([-219.4399, 2037.0703, -40, 0.0021, 0.0019])).squeeze(), @@ -76,7 +73,7 @@ def test_ding2003_dynamics(): def test_ding2007_dynamics(): model = ModelMaker.create_model("ding2007_with_fatigue") assert model.nb_state == 5 - assert model.name_dof == [ + assert model.name_dofs == [ "Cn", "F", "A", @@ -126,14 +123,10 @@ def test_ding2007_dynamics(): np.testing.assert_almost_equal( np.array( model.system_dynamics( - cn=5, - f=100, - a=4920, - tau1=0.050957, - km=0.103, - t=0.11, - t_stim_prev=np.array([0, 0.1]), - pulse_width=0.0002, + time=0.11, + states=[5, 100, 4920, 0.050957, 0.103], + controls=[0.0002], + numerical_timeseries=np.array([0, 0.1]), ) ).squeeze(), np.array(DM([-4.179e02, -4.905e02, -40.0, 2.1759e-03, 2.1677e-03])).squeeze(), @@ -163,7 +156,7 @@ def test_ding2007_dynamics(): def test_hmed2018_dynamics(): model = ModelMaker.create_model("hmed2018_with_fatigue") assert model.nb_state == 5 - assert model.name_dof == [ + assert model.name_dofs == [ "Cn", "F", "A", @@ -215,14 +208,10 @@ def test_hmed2018_dynamics(): np.testing.assert_almost_equal( np.array( model.system_dynamics( - cn=5, - f=100, - a=3009, - tau1=0.050957, - km=0.103, - t=0.11, - t_stim_prev=np.array([0, 0.1]), - pulse_intensity=[30, 50], + time=0.11, + states=[5, 100, 3009, 0.050957, 0.103], + controls=[30, 50], + numerical_timeseries=np.array([0, 0.1]), ) ).squeeze(), np.array(DM([-241, 2037.07, -40, 0.0021, 1.9e-03])).squeeze(), @@ -256,7 +245,7 @@ def test_hmed2018_dynamics(): def test_veltink1992_dynamics(): model = ModelMaker.create_model("veltink_and_riener1998") assert model.nb_state == 2 - assert model.name_dof == [ + assert model.name_dofs == [ "a", "mu", ] @@ -289,9 +278,10 @@ def test_veltink1992_dynamics(): np.testing.assert_almost_equal( np.array( model.system_dynamics( - a=0.5, - mu=0.1, - I=50, + time=None, + states=[0.5, 0.1], + controls=[50], + numerical_timeseries=None, ) ).squeeze(), np.array(DM([0.96153846, 0.01066667])).squeeze(), @@ -306,7 +296,7 @@ def test_veltink1992_dynamics(): def test_marion2009_dynamics(): model = ModelMaker.create_model("marion2009_with_fatigue") assert model.nb_state == 5 - assert model.name_dof == [ + assert model.name_dofs == [ "Cn", "F", "A", @@ -358,14 +348,10 @@ def test_marion2009_dynamics(): np.testing.assert_almost_equal( np.array( model.system_dynamics( - cn=5, - f=100, - a=1473, - tau1=0.04298, - km=0.128, - t=0.11, - t_stim_prev=np.array([0, 0.1]), - theta=20, + time=0.11, + states=[5, 100, 1473, 0.04298, 0.128], + controls=np.array([20]), + numerical_timeseries=np.array([0, 0.1]), ) ).squeeze(), np.array(DM([-2.19408644e02, 1.04853099e03, -2.32015336e01, 1.56300000e-02, 6.26900000e-03])).squeeze(), @@ -376,7 +362,7 @@ def test_marion2009_dynamics(): def test_marion2009_modified_dynamics(): model = ModelMaker.create_model("marion2009_modified_with_fatigue") assert model.nb_state == 5 - assert model.name_dof == [ + assert model.name_dofs == [ "Cn", "F", "A", @@ -432,15 +418,10 @@ def test_marion2009_modified_dynamics(): np.testing.assert_almost_equal( np.array( model.system_dynamics( - cn=5, - f=100, - a=1473, - tau1=0.04298, - km=0.128, - t=0.11, - t_stim_prev=np.array([0, 0.1]), - pulse_width=0.0004, - theta=20, + time=0.11, + states=[5, 100, 1473, 0.04298, 0.128], + controls=np.array([0.0004, 20]), + numerical_timeseries=np.array([0, 0.1]), ) ).squeeze(), np.array(DM([-2.19408644e02, 6.13622448e02, -2.00600000e01, 1.56300000e-02, 6.26900000e-03])).squeeze(), @@ -451,7 +432,7 @@ def test_marion2009_modified_dynamics(): def test_marion2013_dynamics(): model = ModelMaker.create_model("marion2013_with_fatigue") assert model.nb_state == 7 - assert model.name_dof == [ + assert model.name_dofs == [ "Cn", "F", "theta", @@ -517,16 +498,10 @@ def test_marion2013_dynamics(): np.testing.assert_almost_equal( np.array( model.system_dynamics( - cn=5, - f=100, - theta=90, - dtheta_dt=0, - a=2100, - tau1=0.0361, - km=0.352, - t=0.11, - t_stim_prev=np.array([0, 0.1]), - Fload=0, + time=0.11, + states=[5, 100, 90, 0, 2100, 0.0361, 0.352], + controls=np.array([0]), + numerical_timeseries=np.array([0, 0.1]), ) ).squeeze(), np.array( @@ -549,7 +524,7 @@ def test_marion2013_dynamics(): def test_marion2013_modified_dynamics(): model = ModelMaker.create_model("marion2013_modified_with_fatigue") assert model.nb_state == 7 - assert model.name_dof == [ + assert model.name_dofs == [ "Cn", "F", "theta", @@ -619,17 +594,10 @@ def test_marion2013_modified_dynamics(): np.testing.assert_almost_equal( np.array( model.system_dynamics( - cn=5, - f=100, - theta=90, - dtheta_dt=0, - a=2100, - tau1=0.0361, - km=0.352, - t=0.11, - t_stim_prev=np.array([0, 0.1]), - pulse_width=0.0004, - Fload=0, + time=0.11, + states=[5, 100, 90, 0, 2100, 0.0361, 0.352], + controls=np.array([0.0004, 0]), + numerical_timeseries=np.array([0, 0.1]), ) ).squeeze(), np.array( diff --git a/tests/shard1/test_ocp_id.py b/tests/shard1/test_ocp_id.py index 65e8edbc..4185f36a 100644 --- a/tests/shard1/test_ocp_id.py +++ b/tests/shard1/test_ocp_id.py @@ -69,9 +69,8 @@ def prepare_ocp_ding2003( ) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFesId.declare_dynamics( - model=model, - numerical_data_timeseries=numerical_data_time_series, + dynamics_options = OcpFesId.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10), ) @@ -100,7 +99,7 @@ def prepare_ocp_ding2003( return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, x_init=x_init, @@ -238,9 +237,8 @@ def prepare_ocp_ding2007( ) numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFesId.declare_dynamics( - model=model, - numerical_data_timeseries=numerical_data_time_series, + dynamics_options = OcpFesId.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10), ) @@ -274,7 +272,7 @@ def prepare_ocp_ding2007( return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, x_init=x_init, @@ -422,9 +420,8 @@ def prepare_ocp_hmed2018( numerical_data_time_series, stim_idx_at_node_list = model.get_numerical_data_time_series(n_shooting, final_time) - dynamics = OcpFesId.declare_dynamics( - model=model, - numerical_data_timeseries=numerical_data_time_series, + dynamics_options = OcpFesId.declare_dynamics_options( + numerical_time_series=numerical_data_time_series, ode_solver=OdeSolver.RK4(n_integration_steps=10), ) @@ -458,7 +455,7 @@ def prepare_ocp_hmed2018( return OptimalControlProgram( bio_model=[model], - dynamics=dynamics, + dynamics=dynamics_options, n_shooting=n_shooting, phase_time=final_time, x_init=x_init,