BREAKING CHANGE: Move ode_solver from OptimalControlProgram to DynamicsList.add

README.md

@@ -707,7 +707,6 @@ OptimalControlProgram(
     objective_functions: [Objective, ObjectiveList],
     constraints: [Constraint, ConstraintList],
     parameters: ParameterList,
-    ode_solver: OdeSolver,
     control_type: [ControlType, list],
     all_generalized_mapping: BiMapping,
     q_mapping: BiMapping,

bioptim/__init__.py

@@ -161,46 +161,25 @@
 """
 
 from .dynamics.configure_problem import ConfigureProblem, DynamicsFcn, DynamicsList, Dynamics
-from .dynamics.configure_problem import ConfigureProblem, DynamicsFcn, DynamicsList, Dynamics
-from .dynamics.dynamics_evaluation import DynamicsEvaluation
 from .dynamics.dynamics_evaluation import DynamicsEvaluation
 from .dynamics.dynamics_functions import DynamicsFunctions
-from .dynamics.dynamics_functions import DynamicsFunctions
 from .dynamics.fatigue.effort_perception import EffortPerception, TauEffortPerception
-from .dynamics.fatigue.effort_perception import EffortPerception, TauEffortPerception
-from .dynamics.fatigue.fatigue_dynamics import FatigueList
 from .dynamics.fatigue.fatigue_dynamics import FatigueList
 from .dynamics.fatigue.michaud_fatigue import MichaudFatigue, MichaudTauFatigue
-from .dynamics.fatigue.michaud_fatigue import MichaudFatigue, MichaudTauFatigue
 from .dynamics.fatigue.xia_fatigue import XiaFatigue, XiaTauFatigue, XiaFatigueStabilized
-from .dynamics.fatigue.xia_fatigue import XiaFatigue, XiaTauFatigue, XiaFatigueStabilized
-from .dynamics.ode_solvers import OdeSolver, OdeSolverBase
 from .dynamics.ode_solvers import OdeSolver, OdeSolverBase
 from .gui.online_callback_server import PlottingServer
-from .gui.online_callback_server import PlottingServer
 from .gui.plot import CustomPlot
-from .gui.plot import CustomPlot
-from .interfaces import Solver
 from .interfaces import Solver
 from .limits.constraints import ConstraintFcn, ConstraintList, Constraint, ParameterConstraintList
-from .limits.constraints import ConstraintFcn, ConstraintList, Constraint, ParameterConstraintList
-from .limits.fatigue_path_conditions import FatigueBounds, FatigueInitialGuess
 from .limits.fatigue_path_conditions import FatigueBounds, FatigueInitialGuess
 from .limits.multinode_constraint import MultinodeConstraintFcn, MultinodeConstraintList, MultinodeConstraint
-from .limits.multinode_constraint import MultinodeConstraintFcn, MultinodeConstraintList, MultinodeConstraint
 from .limits.multinode_objective import MultinodeObjectiveFcn, MultinodeObjectiveList, MultinodeObjective
-from .limits.multinode_objective import MultinodeObjectiveFcn, MultinodeObjectiveList, MultinodeObjective
-from .limits.objective_functions import ObjectiveFcn, ObjectiveList, Objective, ParameterObjectiveList
 from .limits.objective_functions import ObjectiveFcn, ObjectiveList, Objective, ParameterObjectiveList
 from .limits.path_conditions import BoundsList, InitialGuessList, Bounds, InitialGuess
-from .limits.path_conditions import BoundsList, InitialGuessList, Bounds, InitialGuess
-from .limits.penalty_controller import PenaltyController
 from .limits.penalty_controller import PenaltyController
 from .limits.penalty_helpers import PenaltyHelpers
-from .limits.penalty_helpers import PenaltyHelpers
 from .limits.phase_transition import PhaseTransitionFcn, PhaseTransitionList, PhaseTransition
-from .limits.phase_transition import PhaseTransitionFcn, PhaseTransitionList, PhaseTransition
-from .misc.__version__ import __version__
 from .misc.__version__ import __version__
 from .misc.enums import (
     Axis,

bioptim/dynamics/configure_new_variable.py

@@ -220,7 +220,7 @@ def _declare_cx_and_plot(self):
             for node_index in range(
                 self.nlp.n_states_nodes if self.nlp.phase_dynamics == PhaseDynamics.ONE_PER_NODE else 1
             ):
-                n_cx = self.nlp.ode_solver.n_required_cx + 2
+                n_cx = self.nlp.dynamics_type.ode_solver.n_required_cx + 2
                 cx_scaled = self.define_cx_scaled(n_col=n_cx, node_index=node_index)
                 cx = self.define_cx_unscaled(cx_scaled, self.nlp.x_scaling[self.name].scaling)
                 self.nlp.states.append(
@@ -280,7 +280,7 @@ def _declare_cx_and_plot(self):
             for node_index in range(
                 self.nlp.n_states_nodes if self.nlp.phase_dynamics == PhaseDynamics.ONE_PER_NODE else 1
             ):
-                n_cx = self.nlp.ode_solver.n_required_cx + 2
+                n_cx = self.nlp.dynamics_type.ode_solver.n_required_cx + 2
                 cx_scaled = self.define_cx_scaled(n_col=n_cx, node_index=node_index)
                 cx = self.define_cx_unscaled(cx_scaled, self.nlp.xdot_scaling[self.name].scaling)
                 self.nlp.states_dot.append(
@@ -295,7 +295,7 @@ def _declare_cx_and_plot(self):
             for node_index in range(
                 self.nlp.n_states_nodes if self.nlp.phase_dynamics == PhaseDynamics.ONE_PER_NODE else 1
             ):
-                n_cx = self.nlp.ode_solver.n_required_cx + 2
+                n_cx = self.nlp.dynamics_type.ode_solver.n_required_cx + 2
                 cx_scaled = self.define_cx_scaled(n_col=n_cx, node_index=node_index)
                 cx = self.define_cx_unscaled(cx_scaled, self.nlp.a_scaling[self.name].scaling)
                 self.nlp.algebraic_states.append(

bioptim/dynamics/configure_problem.py

@@ -6,7 +6,7 @@
 from .configure_new_variable import NewVariableConfiguration
 from .dynamics_functions import DynamicsFunctions
 from .fatigue.fatigue_dynamics import FatigueList
-from .ode_solvers import OdeSolver
+from .ode_solvers import OdeSolver, OdeSolverBase
 from ..gui.plot import CustomPlot
 from ..limits.constraints import ImplicitConstraintFcn
 from ..misc.enums import (
@@ -1341,7 +1341,11 @@ def configure_integrated_value(
 
         # TODO: compute values at collocation points
         # but for now only cx_start can be used
-        n_cx = nlp.ode_solver.n_cx - 1 if isinstance(nlp.ode_solver, OdeSolver.COLLOCATION) else 3
+        n_cx = (
+            nlp.dynamics_type.ode_solver.n_cx - 1
+            if isinstance(nlp.dynamics_type.ode_solver, OdeSolver.COLLOCATION)
+            else 3
+        )
         if n_cx < 3:
             n_cx = 3
 
@@ -1955,6 +1959,7 @@ def __init__(
         skip_continuity: bool = False,
         state_continuity_weight: float | None = None,
         phase_dynamics: PhaseDynamics = PhaseDynamics.SHARED_DURING_THE_PHASE,
+        ode_solver: OdeSolver | OdeSolverBase = OdeSolver.RK4(),
         numerical_data_timeseries: dict[str, np.ndarray] = None,
         **extra_parameters: Any,
     ):
@@ -1976,6 +1981,8 @@ def __init__(
             otherwise it is an objective
         phase_dynamics: PhaseDynamics
             If the dynamics should be shared between the nodes or not
+        ode_solver: OdeSolver
+            The integrator to use to integrate this dynamics.
         numerical_data_timeseries: dict[str, np.ndarray]
             The numerical timeseries at each node. ex: the experimental external forces data should go here.
         """
@@ -1994,6 +2001,9 @@ def __init__(
             dynamic_function = extra_parameters["dynamic_function"]
             del extra_parameters["dynamic_function"]
 
+        if not isinstance(ode_solver, OdeSolverBase):
+            raise RuntimeError("ode_solver should be built an instance of OdeSolver")
+
         super(Dynamics, self).__init__(type=dynamics_type, **extra_parameters)
         self.dynamic_function = dynamic_function
         self.configure = configure
@@ -2002,6 +2012,7 @@ def __init__(
         self.skip_continuity = skip_continuity
         self.state_continuity_weight = state_continuity_weight
         self.phase_dynamics = phase_dynamics
+        self.ode_solver = ode_solver
         self.numerical_data_timeseries = numerical_data_timeseries
 
 
@@ -2028,7 +2039,6 @@ def add(self, dynamics_type: Callable | Dynamics | DynamicsFcn, **extra_paramete
         extra_parameters: dict
             Any parameters to pass to Dynamics
         """
-
         if isinstance(dynamics_type, Dynamics):
             self.copy(dynamics_type)
 

bioptim/dynamics/dynamics_functions.py

@@ -154,7 +154,7 @@ def torque_driven(
 
         defects = None
         # TODO: contacts and fatigue to be handled with implicit dynamics
-        if nlp.ode_solver.defects_type == DefectType.IMPLICIT:
+        if nlp.dynamics_type.ode_solver.defects_type == DefectType.IMPLICIT:
             if not with_contact and fatigue is None:
                 qddot = DynamicsFunctions.get(nlp.states_dot["qddot"], nlp.states_dot.scaled.cx)
                 tau_id = DynamicsFunctions.inverse_dynamics(nlp, q, qdot, qddot, with_contact, external_forces)

bioptim/dynamics/ode_solver_base.py

@@ -229,7 +229,7 @@ def initialize_integrator(
             "implicit_ode": nlp.implicit_dynamics_func,
         }
 
-        return nlp.ode_solver.integrator(ode, ode_opt)
+        return nlp.dynamics_type.ode_solver.integrator(ode, ode_opt)
 
     def prepare_dynamic_integrator(self, ocp, nlp):
         """
@@ -244,26 +244,32 @@ def prepare_dynamic_integrator(self, ocp, nlp):
         """
 
         # Primary dynamics
-        dynamics = [nlp.ode_solver.initialize_integrator(ocp, nlp, dynamics_index=0, node_index=0)]
+        dynamics = [nlp.dynamics_type.ode_solver.initialize_integrator(ocp, nlp, dynamics_index=0, node_index=0)]
         if nlp.phase_dynamics == PhaseDynamics.SHARED_DURING_THE_PHASE:
             dynamics = dynamics * nlp.ns
         else:
             for node_index in range(1, nlp.ns):
-                dynamics.append(nlp.ode_solver.initialize_integrator(ocp, nlp, dynamics_index=0, node_index=node_index))
+                dynamics.append(
+                    nlp.dynamics_type.ode_solver.initialize_integrator(
+                        ocp, nlp, dynamics_index=0, node_index=node_index
+                    )
+                )
         nlp.dynamics = dynamics
 
         # Extra dynamics
         extra_dynamics = []
         for i in range(len(nlp.extra_dynamics_func)):
             extra_dynamics += [
-                nlp.ode_solver.initialize_integrator(ocp, nlp, dynamics_index=i, node_index=0, is_extra_dynamics=True)
+                nlp.dynamics_type.ode_solver.initialize_integrator(
+                    ocp, nlp, dynamics_index=i, node_index=0, is_extra_dynamics=True
+                )
             ]
             if nlp.phase_dynamics == PhaseDynamics.SHARED_DURING_THE_PHASE:
                 extra_dynamics = extra_dynamics * nlp.ns
             else:
                 for node_index in range(1, nlp.ns):
                     extra_dynamics += [
-                        nlp.ode_solver.initialize_integrator(
+                        nlp.dynamics_type.ode_solver.initialize_integrator(
                             ocp, nlp, dynamics_index=i, node_index=node_index, is_extra_dynamics=True
                         )
                     ]

bioptim/examples/acados/cube.py

@@ -24,7 +24,7 @@ def prepare_ocp(biorbd_model_path, n_shooting, tf, ode_solver=OdeSolver.RK4(), u
     bio_model = BiorbdModel(biorbd_model_path)
 
     # Dynamics
-    dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics)
+    dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics)
 
     # Path constraint
     x_bounds = BoundsList()
@@ -43,7 +43,6 @@ def prepare_ocp(biorbd_model_path, n_shooting, tf, ode_solver=OdeSolver.RK4(), u
         tf,
         x_bounds=x_bounds,
         u_bounds=u_bounds,
-        ode_solver=ode_solver,
         use_sx=use_sx,
     )
 

bioptim/examples/custom_model/main.py

@@ -73,6 +73,7 @@ def prepare_ocp(
     dynamics.add(
         configure_dynamics,
         dynamic_function=dynamics,
+        ode_solver=ode_solver,
         expand_dynamics=expand_dynamics,
         phase_dynamics=phase_dynamics,
     )
@@ -111,7 +112,6 @@ def prepare_ocp(
         x_init=x_init,
         u_init=u_init,
         objective_functions=objective_functions,
-        ode_solver=ode_solver,
         use_sx=False,
         n_threads=n_threads,
     )

bioptim/examples/fatigue/static_arm_with_fatigue.py

@@ -156,6 +156,7 @@ def prepare_ocp(
         expand_dynamics=expand_dynamics,
         fatigue=fatigue_dynamics,
         with_residual_torque=torque_level > 0,
+        ode_solver=ode_solver,
         phase_dynamics=phase_dynamics,
     )
 
@@ -208,7 +209,6 @@ def prepare_ocp(
         u_init=u_init,
         objective_functions=objective_functions,
         constraints=constraint,
-        ode_solver=ode_solver,
         use_sx=False,
         n_threads=n_threads,
     )

bioptim/examples/getting_started/custom_constraint.py

@@ -101,7 +101,9 @@ def prepare_ocp(
     objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau", weight=100)
 
     # Dynamics
-    dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics)
+    dynamics = Dynamics(
+        DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics
+    )
 
     # Constraints
     constraints = ConstraintList()
@@ -129,7 +131,6 @@ def prepare_ocp(
         u_bounds=u_bounds,
         objective_functions=objective_functions,
         constraints=constraints,
-        ode_solver=ode_solver,
     )
 
 

bioptim/examples/getting_started/custom_dynamics.py

@@ -161,13 +161,15 @@ def prepare_ocp(
             custom_configure,
             dynamic_function=custom_dynamics,
             my_additional_factor=1,
+            ode_solver=ode_solver,
             expand_dynamics=expand_dynamics,
             phase_dynamics=phase_dynamics,
         )
     else:
         dynamics.add(
             DynamicsFcn.TORQUE_DRIVEN,
             dynamic_function=custom_dynamics,
+            ode_solver=ode_solver,
             expand_dynamics=expand_dynamics,
             phase_dynamics=phase_dynamics,
         )
@@ -201,7 +203,6 @@ def prepare_ocp(
         u_bounds=u_bounds,
         objective_functions=objective_functions,
         constraints=constraints,
-        ode_solver=ode_solver,
         use_sx=use_sx,
     )
 

bioptim/examples/getting_started/custom_initial_guess.py

@@ -126,7 +126,9 @@ def prepare_ocp(
     objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau", weight=100)
 
     # Dynamics
-    dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics)
+    dynamics = Dynamics(
+        DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics
+    )
 
     # Constraints
     constraints = ConstraintList()
@@ -248,7 +250,6 @@ def prepare_ocp(
         u_init=u_init,
         objective_functions=objective_functions,
         constraints=constraints,
-        ode_solver=ode_solver,
         x_scaling=x_scaling,
         xdot_scaling=xdot_scaling,
         u_scaling=u_scaling,

bioptim/examples/getting_started/custom_objectives.py

@@ -120,7 +120,9 @@ def prepare_ocp(
     )
 
     # Dynamics
-    dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics)
+    dynamics = Dynamics(
+        DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics
+    )
 
     # Path constraint
     x_bounds = BoundsList()
@@ -145,7 +147,6 @@ def prepare_ocp(
         x_bounds=x_bounds,
         u_bounds=u_bounds,
         objective_functions=objective_functions,
-        ode_solver=ode_solver,
     )
 
 

bioptim/examples/getting_started/custom_parameters.py

@@ -221,7 +221,9 @@ def prepare_ocp(
     objective_functions.add(ObjectiveFcn.Lagrange.MINIMIZE_STATE, key="q", weight=1)
 
     # Dynamics
-    dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics)
+    dynamics = Dynamics(
+        DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics
+    )
 
     # Path constraint
     x_bounds = BoundsList()
@@ -249,7 +251,6 @@ def prepare_ocp(
         parameter_objectives=parameter_objectives,
         parameter_bounds=parameter_bounds,
         parameter_init=parameter_init,
-        ode_solver=ode_solver,
         use_sx=use_sx,
     )
 

bioptim/examples/getting_started/custom_phase_transitions.py

@@ -115,10 +115,18 @@ def prepare_ocp(
 
     # Dynamics
     dynamics = DynamicsList()
-    dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics)
-    dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics)
-    dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics)
-    dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics)
+    dynamics.add(
+        DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics
+    )
+    dynamics.add(
+        DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics
+    )
+    dynamics.add(
+        DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics
+    )
+    dynamics.add(
+        DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics
+    )
 
     # Constraints
     constraints = ConstraintList()
@@ -201,7 +209,6 @@ def prepare_ocp(
         x_init=x_init,
         objective_functions=objective_functions,
         constraints=constraints,
-        ode_solver=ode_solver,
         phase_transitions=phase_transitions,
     )