[RTR] Parameters uniformisation

bioptim/__init__.py

@@ -87,6 +87,8 @@
 # --- Managing the parameters --- #
 ParameterList
     A list of Parameter
+ParameterConstrainer
+    A placeholder for the ParameterList
 
 
 # --- Managing the boundaries of the variables --- #
@@ -217,7 +219,7 @@
     CyclicMovingHorizonEstimator,
     MultiCyclicNonlinearModelPredictiveControl,
 )
-from .optimization.parameters import ParameterList
+from .optimization.parameters import ParameterList, ParameterContainer
 from .optimization.solution.solution import Solution
 from .optimization.solution.solution_data import SolutionMerge, TimeAlignment
 from .optimization.optimization_variable import OptimizationVariableList

bioptim/dynamics/configure_problem.py

@@ -914,7 +914,7 @@ def configure_dynamics_function(ocp, nlp, dyn_func, **extra_params):
             The function to get the derivative of the states
         """
 
-        DynamicsFunctions.apply_parameters(nlp.parameters.mx, nlp)
+        DynamicsFunctions.apply_parameters(nlp)
 
         if not isinstance(dyn_func, (tuple, list)):
             dyn_func = (dyn_func,)
@@ -924,8 +924,8 @@ def configure_dynamics_function(ocp, nlp, dyn_func, **extra_params):
                 nlp.time_mx,
                 nlp.states.scaled.mx_reduced,
                 nlp.controls.scaled.mx_reduced,
-                nlp.parameters.mx,
-                nlp.algebraic_states.scaled.mx,
+                nlp.parameters.scaled.mx_reduced,
+                nlp.algebraic_states.scaled.mx_reduced,
                 nlp,
                 **extra_params,
             )
@@ -941,8 +941,8 @@ def configure_dynamics_function(ocp, nlp, dyn_func, **extra_params):
                         time_span_sym,
                         nlp.states.scaled.mx_reduced,
                         nlp.controls.scaled.mx_reduced,
-                        nlp.parameters.mx,
-                        nlp.algebraic_states.scaled.mx,
+                        nlp.parameters.scaled.mx_reduced,
+                        nlp.algebraic_states.scaled.mx_reduced,
                     ],
                     [dynamics_dxdt],
                     ["t_span", "x", "u", "p", "a"],
@@ -972,8 +972,8 @@ def configure_dynamics_function(ocp, nlp, dyn_func, **extra_params):
                             time_span_sym,
                             nlp.states.scaled.mx_reduced,
                             nlp.controls.scaled.mx_reduced,
-                            nlp.parameters.mx,
-                            nlp.algebraic_states.scaled.mx,
+                            nlp.parameters.scaled.mx_reduced,
+                            nlp.algebraic_states.scaled.mx_reduced,
                             nlp.states_dot.scaled.mx_reduced,
                         ],
                         [dynamics_eval.defects],
@@ -1016,16 +1016,16 @@ def configure_contact_function(ocp, nlp, dyn_func: Callable, **extra_params):
                 time_span_sym,
                 nlp.states.scaled.mx_reduced,
                 nlp.controls.scaled.mx_reduced,
-                nlp.parameters.mx,
-                nlp.algebraic_states.scaled.mx,
+                nlp.parameters.scaled.mx_reduced,
+                nlp.algebraic_states.scaled.mx_reduced,
             ],
             [
                 dyn_func(
                     time_span_sym,
                     nlp.states.scaled.mx_reduced,
                     nlp.controls.scaled.mx_reduced,
-                    nlp.parameters.mx,
-                    nlp.algebraic_states.scaled.mx,
+                    nlp.parameters.scaled.mx_reduced,
+                    nlp.algebraic_states.scaled.mx_reduced,
                     nlp,
                     **extra_params,
                 )
@@ -1079,9 +1079,16 @@ def configure_soft_contact_function(ocp, nlp):
         global_soft_contact_force_func = nlp.model.soft_contact_forces(
             nlp.states["q"].mapping.to_second.map(q), nlp.states["qdot"].mapping.to_second.map(qdot)
         )
+
+        # TODO: do not declare unuseful functions!
         nlp.soft_contact_forces_func = Function(
             "soft_contact_forces_func",
-            [nlp.time_mx, nlp.states.mx_reduced, nlp.controls.mx_reduced, nlp.parameters.mx],
+            [
+                nlp.time_mx,
+                nlp.states.scaled.mx_reduced,
+                nlp.controls.scaled.mx_reduced,
+                nlp.parameters.scaled.mx_reduced,
+            ],
             [global_soft_contact_force_func],
             ["t", "x", "u", "p"],
             ["soft_contact_forces"],

bioptim/dynamics/dynamics_functions.py

@@ -301,6 +301,8 @@ def stochastic_torque_driven(
         q = DynamicsFunctions.get(nlp.states["q"], states)
         qdot = DynamicsFunctions.get(nlp.states["qdot"], states)
         tau = DynamicsFunctions.get(nlp.controls["tau"], controls)
+        motor_noise = DynamicsFunctions.get(nlp.parameters["motor_noise"], parameters)
+        sensory_noise = DynamicsFunctions.get(nlp.parameters["sensory_noise"], parameters)
 
         tau += nlp.model.compute_torques_from_noise_and_feedback(
             nlp=nlp,
@@ -309,8 +311,8 @@ def stochastic_torque_driven(
             controls=controls,
             parameters=parameters,
             algebraic_states=algebraic_states,
-            sensory_noise=nlp.parameters["sensory_noise"].mx,
-            motor_noise=nlp.parameters["motor_noise"].mx,
+            sensory_noise=sensory_noise,
+            motor_noise=motor_noise,
         )
         tau = tau + nlp.model.friction_coefficients @ qdot if with_friction else tau
 
@@ -366,6 +368,8 @@ def stochastic_torque_driven_free_floating_base(
         qdot_roots = DynamicsFunctions.get(nlp.states["qdot_roots"], states)
         qdot_joints = DynamicsFunctions.get(nlp.states["qdot_joints"], states)
         tau_joints = DynamicsFunctions.get(nlp.controls["tau_joints"], controls)
+        motor_noise = DynamicsFunctions.get(nlp.parameters["motor_noise"], parameters)
+        sensory_noise = DynamicsFunctions.get(nlp.parameters["sensory_noise"], parameters)
 
         q_full = vertcat(q_roots, q_joints)
         qdot_full = vertcat(qdot_roots, qdot_joints)
@@ -378,8 +382,8 @@ def stochastic_torque_driven_free_floating_base(
             controls=controls,
             parameters=parameters,
             algebraic_states=algebraic_states,
-            motor_noise=nlp.parameters["motor_noise"].mx,
-            sensory_noise=nlp.parameters["sensory_noise"].mx,
+            motor_noise=motor_noise,
+            sensory_noise=sensory_noise,
         )
         tau_joints = tau_joints + nlp.model.friction_coefficients @ qdot_joints if with_friction else tau_joints
 
@@ -992,7 +996,7 @@ def get(var: OptimizationVariable, cx: MX | SX):
         return var.mapping.to_second.map(cx[var.index, :])
 
     @staticmethod
-    def apply_parameters(parameters: MX.sym, nlp):
+    def apply_parameters(nlp):
         """
         Apply the parameter variables to the model. This should be called before calling the dynamics
 
@@ -1004,10 +1008,13 @@ def apply_parameters(parameters: MX.sym, nlp):
             The definition of the system
         """
 
-        for param in nlp.parameters:
+        for param_key in nlp.parameters:
             # Call the pre dynamics function
-            if param.function[0]:
-                param.function[0](nlp.model, parameters[param.index], **param.params)
+            if nlp.parameters[param_key].function:
+                param = nlp.parameters[param_key]
+                param_scaling = nlp.parameters[param_key].scaling.scaling
+                param_reduced = nlp.parameters.scaled.mx_reduced[param.index]
+                param.function(nlp.model, param_reduced * param_scaling, **param.kwargs)
 
     @staticmethod
     def compute_qdot(nlp, q: MX | SX, qdot: MX | SX):

bioptim/dynamics/integrator.py

@@ -62,9 +62,8 @@ def __init__(self, ode: dict, ode_opt: dict):
         self.cx = ode_opt["cx"]
         self.t_span_sym = ode["t"]
         self.x_sym = ode["x"]
-        self.u_sym = ode["p"]
+        self.u_sym = ode["p"]  # TODO: @ pariterre could we put ode["u"] instead?
         self.param_sym = ode["param"]
-        self.param_scaling = ode_opt["param_scaling"]
         self.a_sym = ode["a"]
         self.fun = ode["ode"]
         self.implicit_fun = ode["implicit_ode"]
@@ -83,7 +82,7 @@ def __init__(self, ode: dict, ode_opt: dict):
             self.dxdt(
                 states=self.x_sym,
                 controls=self.u_sym,
-                params=self.param_sym * self.param_scaling,
+                params=self.param_sym,
                 algebraic_states=self.a_sym,
             ),
             self._input_names,

bioptim/dynamics/ode_solver.py

@@ -157,7 +157,7 @@ def param_ode(self, nlp) -> MX:
         -------
         The symbolic parameters
         """
-        return nlp.parameters.cx
+        return nlp.parameters.scaled.cx_start
 
     def initialize_integrator(self, ocp, nlp, dynamics_index: int, node_index: int, **extra_opt) -> Callable:
         """
@@ -190,7 +190,6 @@ def initialize_integrator(self, ocp, nlp, dynamics_index: int, node_index: int,
             "cx": nlp.cx,
             "control_type": nlp.control_type,
             "defects_type": self.defects_type,
-            "param_scaling": vertcat(*[nlp.parameters[key].scaling.scaling for key in nlp.parameters.keys()]),
             "ode_index": node_index if nlp.dynamics_func[dynamics_index].size2_out("xdot") > 1 else 0,
             "duplicate_starting_point": self.duplicate_starting_point,
             **extra_opt,

bioptim/examples/getting_started/custom_parameters.py

@@ -29,6 +29,7 @@
     ObjectiveList,
     PhaseDynamics,
     VariableScaling,
+    VariableScalingList,
 )
 
 
@@ -175,7 +176,7 @@ def prepare_ocp(
     u_bounds["tau"][1, :] = 0
 
     # Define the parameter to optimize
-    parameters = ParameterList()
+    parameters = ParameterList(use_sx=use_sx)
     parameter_objectives = ParameterObjectiveList()
     parameter_bounds = BoundsList()
     parameter_init = InitialGuessList()

bioptim/examples/stochastic_optimal_control/arm_reaching_muscle_driven.py

@@ -82,8 +82,8 @@ def stochastic_forward_dynamics(
     motor_noise = 0
     sensory_noise = 0
     if with_noise:
-        motor_noise = nlp.parameters["motor_noise"].mx
-        sensory_noise = nlp.parameters["sensory_noise"].mx
+        motor_noise = DynamicsFunctions.get(nlp.parameters["motor_noise"], parameters)
+        sensory_noise = DynamicsFunctions.get(nlp.parameters["sensory_noise"], parameters)
 
     mus_excitations_fb = mus_excitations
     noise_torque = np.zeros(nlp.model.motor_noise_magnitude.shape)

bioptim/examples/stochastic_optimal_control/common.py

@@ -22,8 +22,8 @@ def dynamics_torque_driven_with_feedbacks(time, states, controls, parameters, al
         k = DynamicsFunctions.get(nlp.algebraic_states["k"], algebraic_states)
         k_matrix = StochasticBioModel.reshape_to_matrix(k, nlp.model.matrix_shape_k)
 
-        motor_noise = nlp.parameters["motor_noise"].mx
-        sensory_noise = nlp.parameters["sensory_noise"].mx
+        motor_noise = DynamicsFunctions.get(nlp.parameters["motor_noise"], parameters)
+        sensory_noise = DynamicsFunctions.get(nlp.parameters["sensory_noise"], parameters)
         end_effector = nlp.model.sensory_reference(time, states, controls, parameters, algebraic_states, nlp)
 
         tau_feedback = get_excitation_with_feedback(k_matrix, end_effector, ref, sensory_noise)

bioptim/examples/stochastic_optimal_control/mass_point_model.py

@@ -81,7 +81,8 @@ def dynamics(self, states, controls, parameters, algebraic_states, nlp, with_noi
         u = DynamicsFunctions.get(nlp.controls["u"], controls)
         motor_noise = 0
         if with_noise:
-            motor_noise = nlp.parameters["motor_noise"].mx
+            motor_noise = DynamicsFunctions.get(nlp.parameters["motor_noise"], parameters)
+
         qddot = -self.kapa * (q - u) - self.beta * qdot * sqrt(qdot[0] ** 2 + qdot[1] ** 2 + self.c**2) + motor_noise
 
         return DynamicsEvaluation(dxdt=vertcat(qdot, qddot), defects=None)

bioptim/examples/stochastic_optimal_control/rockit_model.py

@@ -69,7 +69,7 @@ def dynamics(self, states, controls, parameters, algebraic_states, nlp, with_noi
 
         motor_noise = 0
         if with_noise:
-            motor_noise = nlp.parameters["motor_noise"].mx
+            motor_noise = DynamicsFunctions.get(nlp.parameters["motor_noise"], parameters)
 
         qddot = -0.1 * (1 - q**2) * qdot - q + u + motor_noise
 

bioptim/examples/torque_driven_ocp/minimize_maximum_torque_by_extra_parameter.py

@@ -26,6 +26,7 @@
     BiorbdModel,
     PenaltyController,
     ParameterObjectiveList,
+    VariableScaling,
 )
 from matplotlib import pyplot as plt
 
@@ -58,20 +59,12 @@ def prepare_ocp(
     tau_mappings.add("tau", to_second=[None, 0], to_first=[1])
 
     # Define the parameter to optimize
-    parameters = ParameterList()
+    parameters = ParameterList(use_sx=False)
     parameter_init = InitialGuessList()
     parameter_bounds = BoundsList()
 
-    parameters.add(
-        "max_tau",
-        my_parameter_function,
-        size=1,
-    )
-    parameters.add(
-        "min_tau",
-        my_parameter_function,
-        size=1,
-    )
+    parameters.add("max_tau", my_parameter_function, size=1, scaling=VariableScaling("max_tau", [1]))
+    parameters.add("min_tau", my_parameter_function, size=1, scaling=VariableScaling("min_tau", [1]))
 
     parameter_init["max_tau"] = 1
     parameter_init["min_tau"] = -1

bioptim/gui/check_conditioning.py

@@ -85,8 +85,7 @@ def jacobian_hessian_constraints():
                     for controller in controllers:
                         controller.node_index = constraints.node_idx[0]
 
-                    _, t0, x, u, a = constraints.get_variable_inputs(controllers)
-                    p = nlp.parameters.cx
+                    _, t0, x, u, p, a = constraints.get_variable_inputs(controllers)
 
                     list_constraints.append(
                         jacobian(
@@ -187,8 +186,7 @@ def jacobian_hessian_constraints():
 
                         for controller in controllers:
                             controller.node_index = constraints.node_idx[0]
-                        _, t0, x, u, a = constraints.get_variable_inputs(controllers)
-                        p = nlp.parameters.cx
+                        _, t0, x, u, p, a = constraints.get_variable_inputs(controllers)
 
                         hessian_cas = hessian(
                             constraints.function[node_index](t0, phases_dt, x, u, p, a)[axis], vertcat_obj
@@ -390,14 +388,14 @@ def hessian_objective():
 
                 for controller in controllers:
                     controller.node_index = obj.node_idx[0]
-                _, t0, x, u, s = obj.get_variable_inputs(controllers)
-                params = nlp.parameters.cx
+                _, t0, x, u, p, a = obj.get_variable_inputs(controllers)
+
                 target = PenaltyHelpers.target(obj, obj.node_idx.index(node_index))
 
-                p = obj.weighted_function[node_index](t0, phases_dt, x, u, params, s, obj.weight, target)
+                penalty = obj.weighted_function[node_index](t0, phases_dt, x, u, p, a, obj.weight, target)
 
-                for i in range(p.shape[0]):
-                    objective += p[i] ** 2
+                for i in range(penalty.shape[0]):
+                    objective += penalty[i] ** 2
 
             # create function to build the hessian
             vertcat_obj = vertcat([], *nlp.X_scaled, *nlp.U_scaled)  # time, states, controls

bioptim/gui/graph.py

@@ -306,8 +306,8 @@ def print(self):
             print(f"**********")
             print(f"PARAMETERS: ")
             print("")
-            for parameter in self.ocp.nlp[phase_idx].parameters:
-                parameter, initial_guess, min_bound, max_bound, scaling = self._scaling_parameter(parameter.name)
+            for key in self.ocp.nlp[phase_idx].parameters:
+                parameter, initial_guess, min_bound, max_bound, scaling = self._scaling_parameter(key)
                 print(f"Name: {parameter.name}")
                 print(f"Size: {parameter.size}")
                 print(f"Initial_guess: {initial_guess}")
@@ -673,8 +673,8 @@ def _draw_nlp_cluster(self, g, phase_idx: int):
 
             if len(self.ocp.nlp[phase_idx].parameters) > 0:
                 param_idx = 0
-                for param in self.ocp.nlp[phase_idx].parameters:
-                    self._draw_parameter_node(g, phase_idx, param_idx, param.name)
+                for key in self.ocp.nlp[phase_idx].parameters:
+                    self._draw_parameter_node(g, phase_idx, param_idx, key)
                     param_idx += 1
             else:
                 node_str = "<u><b>Parameters</b></u><br/> No parameter set"

bioptim/gui/plot.py

@@ -768,7 +768,7 @@ def _compute_y_from_plot_func(
                     idx,
                     lambda p_idx, n_idx, sn_idx: u[n_idx][:, sn_idx] if n_idx < len(u) else np.ndarray((0, 1)),
                 )
-                p_node = PenaltyHelpers.parameters(penalty, lambda: np.array(p))
+                p_node = PenaltyHelpers.parameters(penalty, 0, lambda p_idx, n_idx, sn_idx: np.array(p))
                 a_node = PenaltyHelpers.states(
                     penalty,
                     idx,

bioptim/interfaces/acados_interface.py

@@ -164,8 +164,8 @@ def __acados_export_model(self, ocp):
         p = ocp.nlp[0].parameters.cx
         a = ocp.nlp[0].algebraic_states.cx_start
         if ocp.parameters:
-            for param in ocp.parameters:
-                if str(param.cx)[:11] == f"time_phase_":
+            for key in ocp.parameters:
+                if str(ocp.parameters[key].cx)[:11] == f"time_phase_":
                     raise RuntimeError("Time constraint not implemented yet with Acados.")
 
         self.nparams = ocp.nlp[0].parameters.shape