feat: pyorerun interfaced

bioptim/examples/getting_started/pendulum.py

@@ -152,10 +152,14 @@ def main():
     # --- Solve the ocp. Please note that online graphics only works with the Linux operating system --- #
     sol = ocp.solve(Solver.IPOPT(show_online_optim=platform.system() == "Linux"))
 
-    # --- Show the results (graph or animation) --- #
+    # --- Show the results graph --- #
     sol.print_cost()
     # sol.graphs(show_bounds=True, save_name="results.png")
-    sol.animate(n_frames=100)
+
+    # --- Animate the solution --- #
+    viewer = "bioviz"
+    # viewer = "pyorerun"
+    sol.animate(n_frames=0, viewer=viewer, show_now=True)
 
     # # --- Saving the solver's output after the optimization --- #
     # Here is an example of how we recommend to save the solution. Please note that sol.ocp is not picklable and that sol will be loaded using the current bioptim version, not the version at the time of the generation of the results.

bioptim/examples/holonomic_constraints/two_pendulums.py

@@ -6,7 +6,6 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
-
 from casadi import MX, Function
 
 from bioptim import (
@@ -213,11 +212,21 @@ def main():
     # --- Show results --- #
     q, qdot, qddot, lambdas = compute_all_states(sol, bio_model)
 
-    import bioviz
+    viewer = "pyorerun"
+    if viewer == "bioviz":
+        import bioviz
+
+        viz = bioviz.Viz(model_path)
+        viz.load_movement(q)
+        viz.exec()
+
+    if viewer == "pyorerun":
+        import pyorerun
+
+        viz = pyorerun.PhaseRerun(t_span=np.concatenate(sol.decision_time()).squeeze())
+        viz.add_animated_model(pyorerun.BiorbdModel(model_path), q=q)
 
-    viz = bioviz.Viz(model_path)
-    viz.load_movement(q)
-    viz.exec()
+        viz.rerun("double_pendulum")
 
     time = sol.decision_time(to_merge=SolutionMerge.NODES)
     plt.title("Lagrange multipliers of the holonomic constraint")

bioptim/examples/moving_horizon_estimation/mhe.py

@@ -15,11 +15,12 @@
 
 from copy import copy
 
-import matplotlib.pyplot as plt
-from scipy.integrate import solve_ivp
+import biorbd_casadi as biorbd
 import casadi as cas
+import matplotlib.pyplot as plt
 import numpy as np
-import biorbd_casadi as biorbd
+from scipy.integrate import solve_ivp
+
 from bioptim import (
     BioModel,
     BiorbdModel,
@@ -232,12 +233,12 @@ def target(i: int):
     sol = mhe.solve(update_functions, **get_solver_options(solver))
     sol_states = sol.decision_states(to_merge=SolutionMerge.NODES)
 
-    print("ACADOS with Bioptim")
+    print(f"{solver} with Bioptim")
     print(f"Window size of MHE : {window_duration} s.")
     print(f"New measurement every : {1 / n_shoot_per_second} s.")
     print(f"Average time per iteration of MHE : {sol.solver_time_to_optimize / (n_frames_total - 1)} s.")
     print(f"Average real time per iteration of MHE : {sol.real_time_to_optimize / (n_frames_total - 1)} s.")
-    print(f"Norm of the error on q = {np.linalg.norm(states[:bio_model.nb_q, :n_frames_total] - sol_states['q'])}")
+    print(f"Norm of the error on q = {np.linalg.norm(states[:bio_model.nb_q, :n_frames_total + 1] - sol_states['q'])}")
 
     markers_estimated = states_to_markers(bio_model, sol_states["q"])
 
@@ -256,7 +257,7 @@ def target(i: int):
     plt.figure()
     plt.plot(sol_states["q"].T, "--", label="States estimate (q)")
     plt.plot(sol_states["qdot"].T, "--", label="States estimate (qdot)")
-    plt.plot(sol_states[:, :n_frames_total].T, label="State truth")
+    plt.plot(states[:, :n_frames_total].T, label="State truth")
     plt.legend()
     plt.show()
 

bioptim/examples/moving_horizon_estimation/models/cart_pendulum.bioMod

@@ -41,7 +41,7 @@ segment    Arm
 	    0.0000	0.0335	-0.0032
 	    0.0000	-0.0032	0.0090
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 endsegment
 
     // Markers

bioptim/examples/muscle_driven_with_contact/muscle_activations_contacts_tracker.py

@@ -1,14 +1,14 @@
 """
-# TODO: Remove all the examples/muscle_driven_with_contact and make sure everything is properly tested
 All the examples in muscle_driven_with_contact are merely to show some dynamics and prepare some OCP for the tests.
 It is not really relevant and will be removed when unitary tests for the dynamics will be implemented
 """
 
 import importlib.util
-from pathlib import Path
 import platform
+from pathlib import Path
 
 import numpy as np
+
 from bioptim import (
     BiorbdModel,
     OptimalControlProgram,
@@ -21,6 +21,7 @@
     OdeSolver,
     Solver,
     SolutionMerge,
+    Node,
 )
 
 # Load track_segment_on_rt
@@ -41,8 +42,12 @@ def prepare_ocp(
 
     # Add objective functions
     objective_functions = ObjectiveList()
-    objective_functions.add(ObjectiveFcn.Lagrange.TRACK_CONTROL, key="muscles", target=muscle_activations_ref)
-    objective_functions.add(ObjectiveFcn.Lagrange.TRACK_CONTACT_FORCES, target=contact_forces_ref)
+    objective_functions.add(
+        ObjectiveFcn.Lagrange.TRACK_CONTROL, key="muscles", target=muscle_activations_ref, node=Node.ALL_SHOOTING
+    )
+    objective_functions.add(
+        ObjectiveFcn.Lagrange.TRACK_CONTACT_FORCES, target=contact_forces_ref, node=Node.ALL_SHOOTING
+    )
     objective_functions.add(ObjectiveFcn.Lagrange.MINIMIZE_STATE, key="qdot", weight=0.001)
     objective_functions.add(ObjectiveFcn.Lagrange.MINIMIZE_STATE, key="q", weight=0.001)
     objective_functions.add(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="muscles", weight=0.001)
@@ -53,40 +58,38 @@ def prepare_ocp(
     dynamics.add(DynamicsFcn.MUSCLE_DRIVEN, with_residual_torque=True, with_contact=True)
 
     # Path constraint
-    n_q = bio_model.nb_q
-    n_qdot = n_q
-    pose_at_first_node = [0, 0, -0.75, 0.75]
+    q_at_first_node = [0, 0, -0.75, 0.75]
 
     # Initialize x_bounds
     x_bounds = BoundsList()
-    x_bounds.add(bounds=bio_model.bounds_from_ranges(["q", "qdot"]))
-    x_bounds[0][:, 0] = pose_at_first_node + [0] * n_qdot
+    x_bounds["q"] = bio_model.bounds_from_ranges("q")
+    x_bounds["qdot"] = bio_model.bounds_from_ranges("qdot")
+    x_bounds["q"][:, 0] = q_at_first_node
 
     # Initial guess
     x_init = InitialGuessList()
-    x_init.add(pose_at_first_node + [0] * n_qdot)
+    x_init["q"] = q_at_first_node
 
     # Define control path constraint
     u_bounds = BoundsList()
-    u_bounds.add(
-        [tau_min] * bio_model.nb_tau + [activation_min] * bio_model.nb_muscles,
-        [tau_max] * bio_model.nb_tau + [activation_max] * bio_model.nb_muscles,
-    )
+    u_bounds["tau"] = ([tau_min] * bio_model.nb_tau, [tau_max] * bio_model.nb_tau)
+    u_bounds["muscles"] = ([activation_min] * bio_model.nb_muscles, [activation_max] * bio_model.nb_muscles)
 
     u_init = InitialGuessList()
-    u_init.add([tau_init] * bio_model.nb_tau + [activation_init] * bio_model.nb_muscles)
+    u_init["tau"] = [tau_init] * bio_model.nb_tau
+    u_init["muscles"] = [activation_init] * bio_model.nb_muscles
 
     # ------------- #
 
     return OptimalControlProgram(
-        bio_model,
-        dynamics,
-        n_shooting,
-        phase_time,
-        x_init,
-        u_init,
-        x_bounds,
-        u_bounds,
+        bio_model=bio_model,
+        dynamics=dynamics,
+        n_shooting=n_shooting,
+        phase_time=phase_time,
+        x_init=x_init,
+        u_init=u_init,
+        x_bounds=x_bounds,
+        u_bounds=u_bounds,
         objective_functions=objective_functions,
         ode_solver=ode_solver,
     )
@@ -112,17 +115,21 @@ def main():
     controls = sol.decision_controls(to_merge=SolutionMerge.NODES)
     q, qdot, tau, mus = states["q"], states["qdot"], controls["tau"], controls["muscles"]
 
-    x = np.concatenate((q, qdot))
-    u = np.concatenate((tau, mus))
-    contact_forces_ref = np.array(ocp_to_track.nlp[0].contact_forces_func(x[:, :-1], u[:, :-1], []))
+    x = np.concatenate((q, qdot), axis=0)
+    u = np.concatenate((tau, mus), axis=0)
+    contact_forces_ref = (
+        np.array([ocp_to_track.nlp[0].contact_forces_func([], x[:, i], u[:, i], [], [], []) for i in range(ns)])
+        .squeeze()
+        .T
+    )
     muscle_activations_ref = mus
 
     # Track these data
     ocp = prepare_ocp(
         biorbd_model_path=model_path,
         phase_time=final_time,
         n_shooting=ns,
-        muscle_activations_ref=muscle_activations_ref[:, :-1],
+        muscle_activations_ref=muscle_activations_ref,
         contact_forces_ref=contact_forces_ref,
     )
 

bioptim/examples/torque_driven_ocp/example_minimize_segment_velocity.py

@@ -5,6 +5,7 @@
 """
 
 import numpy as np
+
 from bioptim import (
     BiorbdModel,
     OptimalControlProgram,

bioptim/examples/torque_driven_ocp/example_multi_biorbd_model.py

@@ -85,16 +85,19 @@ def main():
 
     sol.graphs(show_bounds=True)
 
-    # --- Show results --- #
-    show_solution_animation = False
-    if show_solution_animation:
-        states = sol.decision_states(to_merge=SolutionMerge.NODES)
-        q = states["q"]
-        import bioviz
-
-        b = bioviz.Viz("models/triple_pendulum_both_inertia.bioMod")
-        b.load_movement(q)
-        b.exec()
+    # --- Animate results with bioviz --- #
+    # show_solution_animation = False
+    # if show_solution_animation:
+    #     states = sol.decision_states(to_merge=SolutionMerge.NODES)
+    #     q = states["q"]
+    #     import bioviz
+    #
+    #     b = bioviz.Viz("models/triple_pendulum_both_inertia.bioMod")
+    #     b.load_movement(q)
+    #     b.exec()
+
+    # --- Animate results with pyorerun --- #
+    sol.animate(viewer="pyorerun")
 
 
 if __name__ == "__main__":

bioptim/examples/torque_driven_ocp/example_quaternions.py

@@ -7,9 +7,10 @@
 
 import platform
 
-import numpy as np
 import biorbd_casadi as biorbd
+import numpy as np
 from casadi import MX, Function
+
 from bioptim import (
     BiorbdModel,
     OptimalControlProgram,

bioptim/examples/torque_driven_ocp/models/double_pendulum.bioMod

@@ -10,7 +10,7 @@ segment Seg1
 	    0.0000	0.0335	-0.0032
 	    0.0000	-0.0032	0.0090
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 endsegment
 
     // Marker 1
@@ -38,7 +38,7 @@ segment Seg2
 	    0.0000	0.0335	-0.0032
 	    0.0000	-0.0032	0.0090
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 endsegment
 
     // Marker 3

bioptim/examples/torque_driven_ocp/models/double_pendulum_both_inertia.bioMod

@@ -11,7 +11,7 @@ segment Seg1_1
 	    0.0000	0.0335	-0.0032
 	    0.0000	-0.0032	0.0090
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 	meshcolor 1 0 0
 endsegment
 
@@ -40,7 +40,7 @@ segment Seg2_1
 	    0.0000	0.0335	-0.0032
 	    0.0000	-0.0032	0.0090
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 	meshcolor 1 0 0
 endsegment
 
@@ -66,7 +66,7 @@ segment Seg1_2
 	    0.0000	0.03	-0.0032
 	    0.0000	-0.0032	0.02
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 	meshcolor 0 0 1
 endsegment
 
@@ -95,7 +95,7 @@ segment Seg2_2
 	    0.0000	0.03	-0.0032
 	    0.0000	-0.0032	0.02
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 	meshcolor 0 0 1
 endsegment
 

bioptim/examples/torque_driven_ocp/models/double_pendulum_modified_inertia.bioMod

@@ -10,7 +10,7 @@ segment Seg1
 	    0.0000	0.03	-0.0032
 	    0.0000	-0.0032	0.02
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 endsegment
 
     // Marker 1
@@ -38,7 +38,7 @@ segment Seg2
 	    0.0000	0.03	-0.0032
 	    0.0000	-0.0032	0.02
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 endsegment
 
     // Marker 3

bioptim/examples/torque_driven_ocp/models/double_pendulum_with_translations.bioMod

@@ -12,7 +12,7 @@ segment Seg1
 	    0.0000	0.0335	-0.0032
 	    0.0000	-0.0032	0.0090
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 endsegment
 
     // Marker 1
@@ -40,7 +40,7 @@ segment Seg2
 	    0.0000	0.0335	-0.0032
 	    0.0000	-0.0032	0.0090
     com  -0.0005 0.0688 -0.9542
-    meshfile mesh/pendulum.vtp
+    meshfile mesh/pendulum.STL
 endsegment
 
     // Marker 3