table_generation.py

from gait.generate_table import generate_table as gait_table
from jumper.generate_table import generate_table as jumper_table
from mhe.generate_table import generate_table as mhe_table
from pendulum.generate_table import generate_table as pendulum_table
from pointing.generate_table import generate_table as pointing_table
from somersault.generate_table import generate_table as somersault_table

import numpy as np
from bioptim import Shooting, OdeSolver, SolutionIntegrator


divergence_threshold = 10


class TableOCP:
    def __init__(self):
        self.cols = None

    def add(self, name):
        if not self.cols:
            self.cols = [TableOCP.OCP(name)]
        else:
            self.cols.append(TableOCP.OCP(name))

    def __getitem__(self, item_name):
        return self.cols[[col.name for col in self.cols].index(item_name)]

    def print(self):
        for col in self.cols:
            col.print()

    class OCP:
        def __init__(self, name):
            self.name = name
            self.solver = []

        def print(self):
            print(f"task = {self.name}")
            for solver in self.solver:
                solver.print()

        class Solver:
            def __init__(self, name):
                self.nx = -1
                self.nu = -1
                self.ns = -1
                self.name = name
                self.ode_solver = None
                self.n_iteration = -1
                self.cost = 0
                self.convergence_time = -1
                self.single_shoot_error_t = -1
                self.single_shoot_error_r = -1
                self.single_shoot_divergence_time = None

            def print(self):
                print(f"\t\tsolver = {self.name}")
                steps = f"internal steps = {self.ode_solver.steps}" if isinstance(self.ode_solver, OdeSolver.RK4) \
                    else f"polynomial degree = {self.ode_solver.polynomial_degree}"
                print(f"\t\t\tns = {type(self.ode_solver).__name__}, {steps}")
                print(f"\t\t\tns = {self.ns}")
                print(f"\t\t\tnx = {self.nx}")
                print(f"\t\t\tnu = {self.nu}")
                print(f"\t\t\titerations = {self.n_iteration}")
                print(f"\t\t\tcost = {self.cost}")
                print(f"\t\t\tconvergence_time (s) = {self.convergence_time}")
                print(f"\t\t\tsingle_shoot_error translation (mm) = {self.single_shoot_error_t}")
                print(f"\t\t\tsingle_shoot_error rotation (°) = {self.single_shoot_error_r}")
                print(f"\t\t\tsingle_shoot_time before divergence of {divergence_threshold}(°) = {self.single_shoot_divergence_time}")

            def compute_error_single_shooting(self, sol):
                sol_merged = sol.merge_phases()

                trans_idx = []
                rot_idx = []
                for i in sol.ocp.nlp[0].states["q"].mapping.to_second.map_idx:
                    if i is not None:
                        if sol.ocp.nlp[0].model.nameDof()[i].to_string()[-4:-1] == "Rot":
                            rot_idx += [i]
                        else:
                            trans_idx += [i]
                rot_idx = np.array(list(set(rot_idx)))
                trans_idx = np.array(list(set(trans_idx)))

                sol_int = sol.integrate(shooting_type=Shooting.SINGLE_CONTINUOUS, merge_phases=True, integrator=SolutionIntegrator.SCIPY_RK45)
                jumps = int((sol_merged.states["q"].shape[1] - 1) / (sol_int.states["q"].shape[1] - 1))
                if len(rot_idx) > 0:
                    error_r = sol_int.states["q"][rot_idx, :] - sol_merged.states["q"][rot_idx, ::jumps]
                    self.single_shoot_error_r = (np.sqrt(np.mean(error_r ** 2)) * 180 / np.pi)
                    for i in range(error_r.shape[1]):
                        if np.any(np.abs(error_r[:, i]) > divergence_threshold * np.pi / 180):
                            self.single_shoot_divergence_time = sol_int.phase_time[-1] / sol_int.ns[0] * i
                            break
                else:
                    self.single_shoot_error_r = "N.A."
                if len(trans_idx) > 0:
                    self.single_shoot_error_t = (
                        np.sqrt(
                            np.mean(
                                (sol_int.states["q"][trans_idx, :] - sol_merged.states["q"][trans_idx, ::jumps]) ** 2
                            )
                        )
                        / 1000
                    )
                else:
                    self.single_shoot_error_t = "N.A."


table = TableOCP()

table.add("gait")
table.add("jumper")
table.add("mhe")
table.add("pendulum")
table.add("pointing")
table.add("somersault")

# gait_table(table["gait"])  # requires 16Gb of RAM and takes time to converge (~1h)
jumper_table(table["jumper"])
mhe_table(table["mhe"])
pendulum_table(table["pendulum"])
pointing_table(table["pointing"])
somersault_table(table["somersault"])

table.print()