From 72037727a2e2411a3ce3d275068eede6fb028b78 Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Wed, 10 Nov 2021 16:14:36 -0500
Subject: [PATCH 01/12] bug fix when using show_options

---
 bioptim/interfaces/SolverOptions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/bioptim/interfaces/SolverOptions.py b/bioptim/interfaces/SolverOptions.py
index 5fc960230..d85a68998 100644
--- a/bioptim/interfaces/SolverOptions.py
+++ b/bioptim/interfaces/SolverOptions.py
@@ -369,7 +369,7 @@ def as_dict(self, solver):
             solver_options = self.__dict__
             options = {}
             for key in solver_options:
-                if key != "_c_compile" and key != "type" and key != "show_online_optim":
+                if key != "_c_compile" and key != "type" and key != "show_online_optim" and key != "show_options":
                     ipopt_key = "ipopt." + key[1:]
                     options[ipopt_key] = solver_options[key]
             return {**options, **solver.options_common}

From f283307bb60d33fe6efa39c641acd8b2323fe07f Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Wed, 10 Nov 2021 16:59:51 -0500
Subject: [PATCH 02/12] cleaning as_dict for Ipopt Options

---
 bioptim/interfaces/SolverOptions.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/bioptim/interfaces/SolverOptions.py b/bioptim/interfaces/SolverOptions.py
index d85a68998..b40d8f99d 100644
--- a/bioptim/interfaces/SolverOptions.py
+++ b/bioptim/interfaces/SolverOptions.py
@@ -368,8 +368,9 @@ def set_initialization_options(self, val: float):
         def as_dict(self, solver):
             solver_options = self.__dict__
             options = {}
+            non_python_options = ["_c_compile", "type", "show_online_optim", "show_options"]
             for key in solver_options:
-                if key != "_c_compile" and key != "type" and key != "show_online_optim" and key != "show_options":
+                if key not in non_python_options:
                     ipopt_key = "ipopt." + key[1:]
                     options[ipopt_key] = solver_options[key]
             return {**options, **solver.options_common}

From 2bf1d66147d914da6170fc5234cadb7957e3691f Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Thu, 11 Nov 2021 15:25:16 -0500
Subject: [PATCH 03/12] Add a new enum to choose the integrator of the solution

---
 bioptim/gui/plot.py                           | 14 ++---
 bioptim/misc/enums.py                         | 14 +++++
 .../optimization/optimal_control_program.py   | 10 ++--
 bioptim/optimization/solution.py              | 26 ++++++----
 tests/test_simulate.py                        | 51 ++++++++++---------
 tests/utils.py                                |  3 +-
 6 files changed, 72 insertions(+), 46 deletions(-)

diff --git a/bioptim/gui/plot.py b/bioptim/gui/plot.py
index 9f5eb2eca..bcfcf89ef 100644
--- a/bioptim/gui/plot.py
+++ b/bioptim/gui/plot.py
@@ -10,7 +10,7 @@
 from casadi import Callback, nlpsol_out, nlpsol_n_out, Sparsity, DM
 
 from ..limits.path_conditions import Bounds
-from ..misc.enums import PlotType, ControlType, InterpolationType, Shooting
+from ..misc.enums import PlotType, ControlType, InterpolationType, Shooting, SolutionIntegrator
 from ..misc.mapping import Mapping
 from ..optimization.solution import Solution
 
@@ -197,7 +197,7 @@ def __init__(
         automatically_organize: bool = True,
         show_bounds: bool = False,
         shooting_type: Shooting = Shooting.MULTIPLE,
-        use_scipy_integrator: bool = False,
+        integrator: SolutionIntegrator = False,
     ):
         """
         Prepares the figures during the simulation
@@ -212,8 +212,9 @@ def __init__(
             If the axes should fit the bounds (True) or the data (False)
         shooting_type: Shooting
             The type of integration method
-        use_scipy_integrator: bool
-            Use the scipy solve_ivp integrator for RungeKutta 45 instead of currently defined integrator
+        integrator: SolutionIntegrator
+            use the ode defined by OCP or use a separate integrator provided by scipy
+
         """
         for i in range(1, ocp.n_phases):
             if len(ocp.nlp[0].states["q"]) != len(ocp.nlp[i].states["q"]):
@@ -235,7 +236,8 @@ def __init__(
 
         self.t = []
         self.t_integrated = []
-        self.use_scipy_integrator = use_scipy_integrator
+        self.integrator = integrator
+        self.use_scipy_integrator = self.integrator.value is not None
         if isinstance(self.ocp.original_phase_time, (int, float)):
             self.tf = [self.ocp.original_phase_time]
         else:
@@ -585,7 +587,7 @@ def update_data(self, v: dict):
             continuous=False,
             shooting_type=self.shooting_type,
             keep_intermediate_points=True,
-            use_scipy_integrator=self.use_scipy_integrator,
+            integrator=self.integrator,
         ).states
         data_controls = sol.controls
         data_params = sol.parameters
diff --git a/bioptim/misc/enums.py b/bioptim/misc/enums.py
index 54c548587..eaa546055 100644
--- a/bioptim/misc/enums.py
+++ b/bioptim/misc/enums.py
@@ -103,3 +103,17 @@ class VariableType(Enum):
 
     STATES = "states"
     CONTROLS = "controls"
+
+
+class SolutionIntegrator(Enum):
+    """
+    Selection of integrator to use integrate function
+    """
+    DEFAULT = None
+    SCIPY_RK23 = "RK23"
+    SCIPY_RK45 = 'RK45'
+    SCIPY_DOP853 = "DOP853"
+    SCIPY_BDF = 'BDF'
+    SCIPY_LSODA = "LSODA"
+
+
diff --git a/bioptim/optimization/optimal_control_program.py b/bioptim/optimization/optimal_control_program.py
index 5bc256552..ffe149857 100644
--- a/bioptim/optimization/optimal_control_program.py
+++ b/bioptim/optimization/optimal_control_program.py
@@ -28,7 +28,7 @@
 from ..limits.penalty import PenaltyOption
 from ..limits.objective_functions import ObjectiveFunction
 from ..misc.__version__ import __version__
-from ..misc.enums import ControlType, SolverType, Shooting, PlotType, CostType
+from ..misc.enums import ControlType, SolverType, Shooting, PlotType, CostType, SolutionIntegrator
 from ..misc.mapping import BiMappingList, Mapping
 from ..misc.utils import check_version
 from ..optimization.parameters import ParameterList, Parameter
@@ -749,7 +749,7 @@ def prepare_plots(
         automatically_organize: bool = True,
         show_bounds: bool = False,
         shooting_type: Shooting = Shooting.MULTIPLE,
-        use_scipy_integrator: bool = False,
+        integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
     ) -> PlotOcp:
         """
         Create all the plots associated with the OCP
@@ -762,8 +762,8 @@ def prepare_plots(
             If the ylim should fit the bounds
         shooting_type: Shooting
             What type of integration
-        use_scipy_integrator: bool
-            Use the scipy solve_ivp integrator for RungeKutta 45 instead of currently defined integrator
+        integrator: SolutionIntegrator
+            use the ode defined by OCP or use a separate integrator provided by scipy
 
         Returns
         -------
@@ -775,7 +775,7 @@ def prepare_plots(
             automatically_organize=automatically_organize,
             show_bounds=show_bounds,
             shooting_type=shooting_type,
-            use_scipy_integrator=use_scipy_integrator,
+            integrator=integrator,
         )
 
     def solve(
diff --git a/bioptim/optimization/solution.py b/bioptim/optimization/solution.py
index 4cadc1526..1d8256474 100644
--- a/bioptim/optimization/solution.py
+++ b/bioptim/optimization/solution.py
@@ -9,7 +9,7 @@
 from matplotlib import pyplot as plt
 
 from ..limits.path_conditions import InitialGuess, InitialGuessList
-from ..misc.enums import ControlType, CostType, Shooting, InterpolationType, SolverType
+from ..misc.enums import ControlType, CostType, Shooting, InterpolationType, SolverType, SolutionIntegrator
 from ..misc.utils import check_version
 from ..optimization.non_linear_program import NonLinearProgram
 from ..optimization.optimization_variable import OptimizationVariableList, OptimizationVariable
@@ -492,7 +492,7 @@ def integrate(
         keep_intermediate_points: bool = False,
         merge_phases: bool = False,
         continuous: bool = True,
-        use_scipy_integrator: bool = False,
+        integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
     ) -> Any:
         """
         Integrate the states
@@ -509,8 +509,8 @@ def integrate(
         continuous: bool
             If the arrival value of a node should be discarded [True] or kept [False]. The value of an integrated
             arrival node and the beginning of the next one are expected to be almost equal when the problem converged
-        use_scipy_integrator: bool
-            Ignore the dynamics defined by OCP and use a separate integrator provided by scipy
+        integrator: SolutionIntegrator
+            use the ode defined by OCP or use a separate integrator provided by scipy
 
         Returns
         -------
@@ -535,7 +535,7 @@ def integrate(
                 "Shooting.SINGLE_CONTINUOUS and continuous=False cannot be used simultaneously it is a contradiction"
             )
 
-        out = self.__perform_integration(shooting_type, keep_intermediate_points, continuous, use_scipy_integrator)
+        out = self.__perform_integration(shooting_type, keep_intermediate_points, continuous, integrator)
 
         if merge_phases:
             if continuous:
@@ -548,9 +548,12 @@ def integrate(
         return out
 
     def __perform_integration(
-        self, shooting_type: Shooting, keep_intermediate_points: bool, continuous: bool, use_scipy_integrator: bool
-    ):
+        self, shooting_type: Shooting, keep_intermediate_points: bool, continuous: bool, integrator: SolutionIntegrator):
         n_direct_collocation = sum([nlp.ode_solver.is_direct_collocation for nlp in self.ocp.nlp])
+
+        scipy_integrator = integrator.value
+        use_scipy_integrator = scipy_integrator is not None
+
         if n_direct_collocation > 0 and not use_scipy_integrator:
             if continuous:
                 raise RuntimeError(
@@ -613,8 +616,9 @@ def __perform_integration(
                     n_points = n_steps + 1 if continuous else n_steps
                     t_eval = np.linspace(t_init, t_end, n_points) if keep_intermediate_points else [t_init, t_end]
                     integrated = solve_ivp(
-                        lambda t, x: np.array(nlp.dynamics_func(x, u, params))[:, 0], [t_init, t_end], x0, t_eval=t_eval
-                    ).y
+                        lambda t, x: np.array(nlp.dynamics_func(x, u, params))[:, 0], [t_init, t_end], x0,
+                        t_eval=t_eval,
+                        method=scipy_integrator, atol=1e-10, rtol=1e-10).y
 
                     next_state_col = (
                         (s + 1) * (nlp.ode_solver.steps + 1) if nlp.ode_solver.is_direct_collocation else s + 1
@@ -853,7 +857,7 @@ def graphs(
         show_bounds: bool = False,
         show_now: bool = True,
         shooting_type: Shooting = Shooting.MULTIPLE,
-        use_scipy_integrator: bool = False,
+        integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
     ):
         """
         Show the graphs of the simulation
@@ -875,7 +879,7 @@ def graphs(
         if self.is_merged or self.is_interpolated or self.is_integrated:
             raise NotImplementedError("It is not possible to graph a modified Solution yet")
 
-        plot_ocp = self.ocp.prepare_plots(automatically_organize, show_bounds, shooting_type, use_scipy_integrator)
+        plot_ocp = self.ocp.prepare_plots(automatically_organize, show_bounds, shooting_type, integrator)
         plot_ocp.update_data(self.vector)
         if show_now:
             plt.show()
diff --git a/tests/test_simulate.py b/tests/test_simulate.py
index 382475542..68a73ad22 100644
--- a/tests/test_simulate.py
+++ b/tests/test_simulate.py
@@ -4,7 +4,7 @@
 import pytest
 
 import numpy as np
-from bioptim import Shooting, OdeSolver
+from bioptim import Shooting, OdeSolver, SolutionIntegrator
 
 
 def test_merge_phases_one_phase():
@@ -167,8 +167,8 @@ def test_interpolate_multiphases_merge_phase():
 
 
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
-@pytest.mark.parametrize("use_scipy", [True, False])
-def test_integrate(use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.SCIPY_RK45, SolutionIntegrator.DEFAULT])
+def test_integrate(integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import pendulum as ocp_module
 
@@ -185,7 +185,7 @@ def test_integrate(use_scipy, ode_solver):
 
     sol = ocp.solve()
 
-    opts = {"shooting_type": Shooting.MULTIPLE, "keep_intermediate_points": False, "use_scipy_integrator": use_scipy}
+    opts = {"shooting_type": Shooting.MULTIPLE, "keep_intermediate_points": False, "integrator": integrator}
     with pytest.raises(
         ValueError,
         match="Shooting.MULTIPLE and keep_intermediate_points=False "
@@ -194,7 +194,7 @@ def test_integrate(use_scipy, ode_solver):
         _ = sol.integrate(**opts)
 
     opts["keep_intermediate_points"] = True
-    if ode_solver == OdeSolver.COLLOCATION and not use_scipy:
+    if ode_solver == OdeSolver.COLLOCATION and integrator.value is None:
         with pytest.raises(RuntimeError, match="Integration with direct collocation must be not continuous"):
             sol.integrate(**opts)
         return
@@ -202,7 +202,7 @@ def test_integrate(use_scipy, ode_solver):
     sol_integrated = sol.integrate(**opts)
     shapes = (4, 2, 2)
 
-    decimal = 5 if use_scipy else 8
+    decimal = 5 if integrator.value is not None else 8
     for i, key in enumerate(sol.states):
         np.testing.assert_almost_equal(
             sol_integrated.states[key][:, [0, -1]], sol.states[key][:, [0, -1]], decimal=decimal
@@ -241,7 +241,7 @@ def test_integrate_single_shoot(keep_intermediate_points, ode_solver):
 
     sol = ocp.solve()
 
-    opts = {"keep_intermediate_points": keep_intermediate_points, "use_scipy_integrator": False}
+    opts = {"keep_intermediate_points": keep_intermediate_points, "integrator": SolutionIntegrator.DEFAULT}
     if ode_solver == OdeSolver.COLLOCATION:
         with pytest.raises(RuntimeError, match="Integration with direct collocation must be not continuous"):
             sol.integrate(**opts)
@@ -296,7 +296,7 @@ def test_integrate_single_shoot_use_scipy(keep_intermediate_points, ode_solver):
 
     sol = ocp.solve()
 
-    opts = {"keep_intermediate_points": keep_intermediate_points, "use_scipy_integrator": True}
+    opts = {"keep_intermediate_points": keep_intermediate_points, "integrator": SolutionIntegrator.SCIPY_RK45}
 
     sol_integrated = sol.integrate(**opts)
     shapes = (4, 2, 2)
@@ -436,12 +436,13 @@ def test_integrate_single_shoot_use_scipy(keep_intermediate_points, ode_solver):
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
 @pytest.mark.parametrize("merge", [False, True])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_non_continuous(shooting, merge, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_non_continuous(shooting, merge, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import pendulum as ocp_module
 
     bioptim_folder = os.path.dirname(ocp_module.__file__)
+    use_scipy = integrator is not None
 
     n_shooting = 10
 
@@ -458,7 +459,7 @@ def test_integrate_non_continuous(shooting, merge, use_scipy, ode_solver):
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": False,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
     if shooting == Shooting.SINGLE_CONTINUOUS:
         with pytest.raises(
@@ -515,8 +516,8 @@ def test_integrate_non_continuous(shooting, merge, use_scipy, ode_solver):
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
 @pytest.mark.parametrize("keep_intermediate_points", [True, False])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_multiphase(shooting, keep_intermediate_points, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import example_multiphase as ocp_module
 
@@ -527,11 +528,12 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, use_scipy, ode
     sol = ocp.solve()
     n_shooting = [20, 30, 20]
 
+    use_scipy = integrator is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": keep_intermediate_points,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
 
     if shooting == Shooting.SINGLE_CONTINUOUS:
@@ -603,8 +605,8 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, use_scipy, ode
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
 @pytest.mark.parametrize("keep_intermediate_points", [True, False])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_multiphase_merged(shooting, keep_intermediate_points, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_multiphase_merged(shooting, keep_intermediate_points, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import example_multiphase as ocp_module
 
@@ -617,11 +619,12 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, use_sci
 
     sol = ocp.solve()
 
+    use_scipy= integrator.value is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": keep_intermediate_points,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
 
     if shooting == Shooting.SINGLE_CONTINUOUS:
@@ -698,8 +701,8 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, use_sci
 
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_multiphase_non_continuous(shooting, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_multiphase_non_continuous(shooting, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import example_multiphase as ocp_module
 
@@ -710,11 +713,12 @@ def test_integrate_multiphase_non_continuous(shooting, use_scipy, ode_solver):
     sol = ocp.solve()
     n_shooting = [20, 30, 20]
 
+    use_scipy = integrator.value is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": True,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
 
     if shooting == Shooting.SINGLE_CONTINUOUS:
@@ -768,8 +772,8 @@ def test_integrate_multiphase_non_continuous(shooting, use_scipy, ode_solver):
 
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_multiphase_merged_non_continuous(shooting, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_multiphase_merged_non_continuous(shooting, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import example_multiphase as ocp_module
 
@@ -779,11 +783,12 @@ def test_integrate_multiphase_merged_non_continuous(shooting, use_scipy, ode_sol
 
     sol = ocp.solve()
 
+    use_scipy = integrator.value is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": False,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
     if shooting == Shooting.SINGLE_CONTINUOUS:
         with pytest.raises(
diff --git a/tests/utils.py b/tests/utils.py
index dfa1a47a0..2526a58e7 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -18,6 +18,7 @@
     InitialGuess,
     Shooting,
     Solver,
+    SolutionIntegrator,
 )
 
 
@@ -117,7 +118,7 @@ def simulate(sol, decimal_value=7):
             merge_phases=True,
             shooting_type=Shooting.SINGLE_CONTINUOUS,
             keep_intermediate_points=True,
-            use_scipy_integrator=False,
+            integrator=SolutionIntegrator.DEFAULT,
         )
 
         # Evaluate the final error of the single shooting integration versus the finale node

From 2d2b82472b2aee821f2be82391603384607d742f Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Thu, 11 Nov 2021 15:25:58 -0500
Subject: [PATCH 04/12] modified init

---
 bioptim/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/bioptim/__init__.py b/bioptim/__init__.py
index f534d65bc..493295338 100644
--- a/bioptim/__init__.py
+++ b/bioptim/__init__.py
@@ -165,7 +165,7 @@
 from .limits.path_conditions import BoundsList, Bounds, InitialGuessList, InitialGuess, QAndQDotBounds
 from .limits.fatigue_path_conditions import FatigueBounds, FatigueInitialGuess
 from .limits.penalty_node import PenaltyNode, PenaltyNodeList
-from .misc.enums import Axis, Node, InterpolationType, PlotType, ControlType, CostType, Shooting, VariableType
+from .misc.enums import Axis, Node, InterpolationType, PlotType, ControlType, CostType, Shooting, VariableType, SolutionIntegrator
 from .misc.mapping import BiMappingList, BiMapping, Mapping
 from .optimization.non_linear_program import NonLinearProgram
 from .optimization.optimal_control_program import OptimalControlProgram

From 7ac711392da9155e40228d589b002a3d7a193ce5 Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Thu, 11 Nov 2021 15:59:19 -0500
Subject: [PATCH 05/12] readme update

---
 README.md | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/README.md b/README.md
index 8a1385792..872d85776 100644
--- a/README.md
+++ b/README.md
@@ -88,6 +88,7 @@ As a tour guide that uses this binder, you can watch the `bioptim` workshop that
   - [InterpolationType](#enum-interpolationtype)
   - [Shooting](#enum-shooting)
   - [CostType](#enum-costtype)
+  - [SolutionIntegrator](#enum-solutionintegrator)
         
 [Examples](#examples)
 - [Run examples](#run-examples)
@@ -1396,6 +1397,16 @@ The type of cost
 - CONSTRAINTS: The constraints
 - ALL: All the previously described cost type
 
+### Enum: SolutionIntegrator
+The type of integrator used to integrate the solution of the optimal control problem
+- DEFAULT: The default integrator initially chosen with [OdeSolver](#class-odesolver)
+- SCIPY_RK23: The scipy integrator RK23
+- SCIPY_RK45: The scipy integrator RK45
+- SCIPY_DOP853: The scipy integrator DOP853
+- SCIPY_BDF: The scipy integrator BDF
+- SCIPY_LSODA: The scipy integrator LSODA
+
+
 # Examples
 In this section, you will find the description of all the examples implemented with bioptim. They are ordered in 
 separate files. Each subsection corresponds to the different files, dealing with different examples and topics.

From eb15185e0b9f338b9441a3ec67078100cb1c7b5f Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Thu, 11 Nov 2021 16:48:37 -0500
Subject: [PATCH 06/12] readme and i took into account review

---
 README.md                                     |  1 -
 bioptim/__init__.py                           | 12 ++++-
 bioptim/gui/plot.py                           | 24 ++++++---
 bioptim/misc/enums.py                         |  7 ++-
 .../optimization/optimal_control_program.py   |  2 +-
 bioptim/optimization/solution.py              | 37 +++++++------
 tests/test_simulate.py                        | 53 +++++++++----------
 7 files changed, 80 insertions(+), 56 deletions(-)

diff --git a/README.md b/README.md
index 872d85776..9bd008b2b 100644
--- a/README.md
+++ b/README.md
@@ -1384,7 +1384,6 @@ The accepted values are:
 - SPLINE: Requires five columns. It performs a cubic spline to interpolate between the nodes.
 - CUSTOM: User defined interpolation function
 
-
 ### Enum: Shooting
 The type of integration to perform
 - MULTIPLE: resets the state at each node
diff --git a/bioptim/__init__.py b/bioptim/__init__.py
index 493295338..d95de9b61 100644
--- a/bioptim/__init__.py
+++ b/bioptim/__init__.py
@@ -165,7 +165,17 @@
 from .limits.path_conditions import BoundsList, Bounds, InitialGuessList, InitialGuess, QAndQDotBounds
 from .limits.fatigue_path_conditions import FatigueBounds, FatigueInitialGuess
 from .limits.penalty_node import PenaltyNode, PenaltyNodeList
-from .misc.enums import Axis, Node, InterpolationType, PlotType, ControlType, CostType, Shooting, VariableType, SolutionIntegrator
+from .misc.enums import (
+    Axis,
+    Node,
+    InterpolationType,
+    PlotType,
+    ControlType,
+    CostType,
+    Shooting,
+    VariableType,
+    SolutionIntegrator,
+)
 from .misc.mapping import BiMappingList, BiMapping, Mapping
 from .optimization.non_linear_program import NonLinearProgram
 from .optimization.optimal_control_program import OptimalControlProgram
diff --git a/bioptim/gui/plot.py b/bioptim/gui/plot.py
index bcfcf89ef..967e4f489 100644
--- a/bioptim/gui/plot.py
+++ b/bioptim/gui/plot.py
@@ -197,7 +197,7 @@ def __init__(
         automatically_organize: bool = True,
         show_bounds: bool = False,
         shooting_type: Shooting = Shooting.MULTIPLE,
-        integrator: SolutionIntegrator = False,
+        integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
     ):
         """
         Prepares the figures during the simulation
@@ -213,7 +213,7 @@ def __init__(
         shooting_type: Shooting
             The type of integration method
         integrator: SolutionIntegrator
-            use the ode defined by OCP or use a separate integrator provided by scipy
+             Use the ode defined by OCP or use a separate integrator provided by scipy
 
         """
         for i in range(1, ocp.n_phases):
@@ -237,7 +237,7 @@ def __init__(
         self.t = []
         self.t_integrated = []
         self.integrator = integrator
-        self.use_scipy_integrator = self.integrator.value is not None
+
         if isinstance(self.ocp.original_phase_time, (int, float)):
             self.tf = [self.ocp.original_phase_time]
         else:
@@ -296,12 +296,14 @@ def __update_time_vector(self):
         self.t_integrated = []
         last_t = 0
         for phase_idx, nlp in enumerate(self.ocp.nlp):
-            n_int_steps = nlp.ode_solver.steps_scipy if self.use_scipy_integrator else nlp.ode_solver.steps
+            n_int_steps = (
+                nlp.ode_solver.steps_scipy if self.integrator != SolutionIntegrator.DEFAULT else nlp.ode_solver.steps
+            )
             dt_ns = self.tf[phase_idx] / nlp.ns
             time_phase_integrated = []
             last_t_int = copy(last_t)
             for _ in range(nlp.ns):
-                if nlp.ode_solver.is_direct_collocation and not self.use_scipy_integrator:
+                if nlp.ode_solver.is_direct_collocation and self.integrator == SolutionIntegrator.DEFAULT:
                     time_phase_integrated.append(np.array(nlp.dynamics[0].step_time) * dt_ns + last_t_int)
                 else:
                     time_phase_integrated.append(np.linspace(last_t_int, last_t_int + dt_ns, n_int_steps + 1))
@@ -428,7 +430,11 @@ def legend_without_duplicate_labels(ax):
                         )
                     elif plot_type == PlotType.INTEGRATED:
                         plots_integrated = []
-                        n_int_steps = nlp.ode_solver.steps_scipy if self.use_scipy_integrator else nlp.ode_solver.steps
+                        n_int_steps = (
+                            nlp.ode_solver.steps_scipy
+                            if self.integrator != SolutionIntegrator.DEFAULT
+                            else nlp.ode_solver.steps
+                        )
                         zero = np.zeros(n_int_steps + 1)
                         color = self.plot_func[variable][i].color if self.plot_func[variable][i].color else "tab:brown"
                         for cmp in range(nlp.ns):
@@ -600,7 +606,11 @@ def update_data(self, v: dict):
             self.__update_xdata()
 
         for i, nlp in enumerate(self.ocp.nlp):
-            step_size = nlp.ode_solver.steps_scipy + 1 if self.use_scipy_integrator else nlp.ode_solver.steps + 1
+            step_size = (
+                nlp.ode_solver.steps_scipy + 1
+                if self.integrator != SolutionIntegrator.DEFAULT
+                else nlp.ode_solver.steps + 1
+            )
             n_elements = nlp.ns * step_size + 1
 
             state = np.ndarray((0, n_elements))
diff --git a/bioptim/misc/enums.py b/bioptim/misc/enums.py
index eaa546055..22cd0e2fc 100644
--- a/bioptim/misc/enums.py
+++ b/bioptim/misc/enums.py
@@ -109,11 +109,10 @@ class SolutionIntegrator(Enum):
     """
     Selection of integrator to use integrate function
     """
+
     DEFAULT = None
     SCIPY_RK23 = "RK23"
-    SCIPY_RK45 = 'RK45'
+    SCIPY_RK45 = "RK45"
     SCIPY_DOP853 = "DOP853"
-    SCIPY_BDF = 'BDF'
+    SCIPY_BDF = "BDF"
     SCIPY_LSODA = "LSODA"
-
-
diff --git a/bioptim/optimization/optimal_control_program.py b/bioptim/optimization/optimal_control_program.py
index ffe149857..71cc7c496 100644
--- a/bioptim/optimization/optimal_control_program.py
+++ b/bioptim/optimization/optimal_control_program.py
@@ -763,7 +763,7 @@ def prepare_plots(
         shooting_type: Shooting
             What type of integration
         integrator: SolutionIntegrator
-            use the ode defined by OCP or use a separate integrator provided by scipy
+            Use the ode defined by OCP or use a separate integrator provided by scipy
 
         Returns
         -------
diff --git a/bioptim/optimization/solution.py b/bioptim/optimization/solution.py
index 1d8256474..b20a4d67b 100644
--- a/bioptim/optimization/solution.py
+++ b/bioptim/optimization/solution.py
@@ -510,7 +510,7 @@ def integrate(
             If the arrival value of a node should be discarded [True] or kept [False]. The value of an integrated
             arrival node and the beginning of the next one are expected to be almost equal when the problem converged
         integrator: SolutionIntegrator
-            use the ode defined by OCP or use a separate integrator provided by scipy
+            Use the ode defined by OCP or use a separate integrator provided by scipy
 
         Returns
         -------
@@ -548,27 +548,25 @@ def integrate(
         return out
 
     def __perform_integration(
-        self, shooting_type: Shooting, keep_intermediate_points: bool, continuous: bool, integrator: SolutionIntegrator):
+        self, shooting_type: Shooting, keep_intermediate_points: bool, continuous: bool, integrator: SolutionIntegrator
+    ):
         n_direct_collocation = sum([nlp.ode_solver.is_direct_collocation for nlp in self.ocp.nlp])
 
-        scipy_integrator = integrator.value
-        use_scipy_integrator = scipy_integrator is not None
-
-        if n_direct_collocation > 0 and not use_scipy_integrator:
+        if n_direct_collocation > 0 and integrator != SolutionIntegrator.DEFAULT:
             if continuous:
                 raise RuntimeError(
-                    "Integration with direct collocation must be not continuous if not use_scipy_integrator"
+                    "Integration with direct collocation must be not continuous if a scipy integrator is used"
                 )
 
             if shooting_type != Shooting.MULTIPLE:
                 raise RuntimeError(
                     "Integration with direct collocation must using shooting_type=Shooting.MULTIPLE "
-                    "if not use_scipy_integrator"
+                    "if not a scipy integrator is not used"
                 )
 
         # Copy the data
         out = self.copy(skip_data=True)
-        out.recomputed_time_steps = use_scipy_integrator
+        out.recomputed_time_steps = integrator != SolutionIntegrator.DEFAULT
         out._states = []
         for _ in range(len(self._states)):
             out._states.append({})
@@ -578,7 +576,7 @@ def __perform_integration(
         for p, nlp in enumerate(self.ocp.nlp):
             param_scaling = nlp.parameters.scaling
             n_states = self._states[p]["all"].shape[0]
-            n_steps = nlp.ode_solver.steps_scipy if use_scipy_integrator else nlp.ode_solver.steps
+            n_steps = nlp.ode_solver.steps_scipy if integrator == SolutionIntegrator.DEFAULT else nlp.ode_solver.steps
             if not continuous:
                 n_steps += 1
             if keep_intermediate_points:
@@ -599,7 +597,11 @@ def __perform_integration(
                         )
                     x0 += np.array(val)[:, 0]
             else:
-                col = slice(0, n_steps) if nlp.ode_solver.is_direct_collocation and not use_scipy_integrator else 0
+                col = (
+                    slice(0, n_steps)
+                    if nlp.ode_solver.is_direct_collocation and integrator != SolutionIntegrator.DEFAULT
+                    else 0
+                )
                 x0 = self._states[p]["all"][:, col]
 
             for s in range(self.ns[p]):
@@ -610,15 +612,20 @@ def __perform_integration(
                 else:
                     raise NotImplementedError(f"ControlType {nlp.control_type} " f"not yet implemented in integrating")
 
-                if use_scipy_integrator:
+                if integrator != SolutionIntegrator.DEFAULT:
                     t_init = sum(out.phase_time[:p]) / nlp.ns
                     t_end = sum(out.phase_time[: (p + 2)]) / nlp.ns
                     n_points = n_steps + 1 if continuous else n_steps
                     t_eval = np.linspace(t_init, t_end, n_points) if keep_intermediate_points else [t_init, t_end]
                     integrated = solve_ivp(
-                        lambda t, x: np.array(nlp.dynamics_func(x, u, params))[:, 0], [t_init, t_end], x0,
+                        lambda t, x: np.array(nlp.dynamics_func(x, u, params))[:, 0],
+                        [t_init, t_end],
+                        x0,
                         t_eval=t_eval,
-                        method=scipy_integrator, atol=1e-10, rtol=1e-10).y
+                        method=integrator.value,
+                        atol=1e-10,
+                        rtol=1e-10,
+                    ).y
 
                     next_state_col = (
                         (s + 1) * (nlp.ode_solver.steps + 1) if nlp.ode_solver.is_direct_collocation else s + 1
@@ -872,7 +879,7 @@ def graphs(
             If the show method should be called. This is blocking
         shooting_type: Shooting
             The type of interpolation
-        use_scipy_integrator: bool
+        integrator: SolutionIntegrator
             Use the scipy solve_ivp integrator for RungeKutta 45 instead of currently defined integrator
         """
 
diff --git a/tests/test_simulate.py b/tests/test_simulate.py
index 68a73ad22..0849710bb 100644
--- a/tests/test_simulate.py
+++ b/tests/test_simulate.py
@@ -202,7 +202,7 @@ def test_integrate(integrator, ode_solver):
     sol_integrated = sol.integrate(**opts)
     shapes = (4, 2, 2)
 
-    decimal = 5 if integrator.value is not None else 8
+    decimal = 5 if integrator != SolutionIntegrator.DEFAULT else 8
     for i, key in enumerate(sol.states):
         np.testing.assert_almost_equal(
             sol_integrated.states[key][:, [0, -1]], sol.states[key][:, [0, -1]], decimal=decimal
@@ -442,7 +442,6 @@ def test_integrate_non_continuous(shooting, merge, integrator, ode_solver):
     from bioptim.examples.getting_started import pendulum as ocp_module
 
     bioptim_folder = os.path.dirname(ocp_module.__file__)
-    use_scipy = integrator is not None
 
     n_shooting = 10
 
@@ -478,7 +477,11 @@ def test_integrate_non_continuous(shooting, merge, integrator, ode_solver):
 
     opts["keep_intermediate_points"] = True
     opts["merge_phases"] = merge
-    if ode_solver == OdeSolver.COLLOCATION and shooting != Shooting.MULTIPLE and not use_scipy:
+    if (
+        ode_solver == OdeSolver.COLLOCATION
+        and shooting != Shooting.MULTIPLE
+        and integrator != SolutionIntegrator.DEFAULT
+    ):
         with pytest.raises(
             RuntimeError,
             match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE",
@@ -489,14 +492,14 @@ def test_integrate_non_continuous(shooting, merge, integrator, ode_solver):
     sol_integrated = sol.integrate(**opts)
     shapes = (4, 2, 2)
 
-    decimal = 1 if use_scipy or ode_solver == OdeSolver.COLLOCATION else 8
+    decimal = 1 if integrator != SolutionIntegrator.DEFAULT or ode_solver == OdeSolver.COLLOCATION else 8
     for i, key in enumerate(sol.states):
         np.testing.assert_almost_equal(
             sol_integrated.states[key][:, [0, -1]], sol.states[key][:, [0, -1]], decimal=decimal
         )
 
         if ode_solver == OdeSolver.COLLOCATION:
-            if use_scipy:
+            if integrator != SolutionIntegrator.DEFAULT:
                 assert sol_integrated.states[key].shape == (shapes[i], n_shooting * (5 + 1) + 1)
             else:
                 assert sol_integrated.states[key].shape == (shapes[i], n_shooting * (4 + 1) + 1)
@@ -528,7 +531,6 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, od
     sol = ocp.solve()
     n_shooting = [20, 30, 20]
 
-    use_scipy = integrator is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
@@ -545,7 +547,7 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, od
         return
 
     if ode_solver == OdeSolver.COLLOCATION:
-        if not use_scipy:
+        if integrator != SolutionIntegrator.DEFAULT:
             if shooting != Shooting.MULTIPLE:
                 with pytest.raises(
                     RuntimeError, match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE"
@@ -564,7 +566,7 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, od
             return
 
     opts["continuous"] = True
-    if ode_solver == OdeSolver.COLLOCATION and not use_scipy:
+    if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
         with pytest.raises(
             RuntimeError,
             match="Integration with direct collocation must be not continuous",
@@ -575,10 +577,10 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, od
     sol_integrated = sol.integrate(**opts)
     shapes = (6, 3, 3)
 
-    decimal = 1 if use_scipy else 8
+    decimal = 1 if integrator != SolutionIntegrator.DEFAULT else 8
     for i in range(len(sol_integrated.states)):
         for k, key in enumerate(sol.states[i]):
-            if not use_scipy or shooting == Shooting.MULTIPLE:
+            if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
                 np.testing.assert_almost_equal(
                     sol_integrated.states[i][key][:, [0, -1]], sol.states[i][key][:, [0, -1]], decimal=decimal
                 )
@@ -586,7 +588,7 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, od
             if keep_intermediate_points:
                 assert sol_integrated.states[i][key].shape == (shapes[k], n_shooting[i] * 5 + 1)
             else:
-                if not use_scipy or shooting == Shooting.MULTIPLE:
+                if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
                     np.testing.assert_almost_equal(sol_integrated.states[i][key], sol.states[i][key])
                 assert sol_integrated.states[i][key].shape == (shapes[k], n_shooting[i] + 1)
             if ode_solver == OdeSolver.COLLOCATION:
@@ -619,7 +621,6 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, integra
 
     sol = ocp.solve()
 
-    use_scipy= integrator.value is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
@@ -635,7 +636,7 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, integra
             _ = sol.integrate(**opts)
         return
 
-    if ode_solver == OdeSolver.COLLOCATION and not use_scipy:
+    if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
         if shooting != Shooting.MULTIPLE:
             with pytest.raises(
                 RuntimeError, match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE"
@@ -655,7 +656,7 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, integra
 
     opts["merge_phases"] = True
     opts["continuous"] = True
-    if ode_solver == OdeSolver.COLLOCATION and not use_scipy:
+    if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
         with pytest.raises(
             RuntimeError,
             match="Integration with direct collocation must be not continuous",
@@ -667,17 +668,17 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, integra
     sol_integrated = sol.integrate(**opts)
     shapes = (6, 3, 3)
 
-    decimal = 0 if use_scipy else 8
+    decimal = 0 if integrator != SolutionIntegrator.DEFAULT else 8
     for k, key in enumerate(sol.states[0]):
         expected = np.array([sol.states[0][key][:, 0], sol.states[-1][key][:, -1]]).T
-        if not use_scipy or shooting == Shooting.MULTIPLE:
+        if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
             np.testing.assert_almost_equal(sol_integrated.states[key][:, [0, -1]], expected, decimal=decimal)
 
         if keep_intermediate_points:
             assert sol_integrated.states[key].shape == (shapes[k], sum(n_shooting) * 5 + 1)
         else:
             # The interpolation prevents from comparing all points
-            if not use_scipy or shooting == Shooting.MULTIPLE:
+            if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
                 expected = np.concatenate(
                     (sol.states[0][key][:, 0:1], sol.states[-1][key][:, -1][:, np.newaxis]), axis=1
                 )
@@ -713,7 +714,6 @@ def test_integrate_multiphase_non_continuous(shooting, integrator, ode_solver):
     sol = ocp.solve()
     n_shooting = [20, 30, 20]
 
-    use_scipy = integrator.value is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
@@ -730,7 +730,7 @@ def test_integrate_multiphase_non_continuous(shooting, integrator, ode_solver):
             _ = sol.integrate(**opts)
         return
 
-    if ode_solver == OdeSolver.COLLOCATION and not use_scipy:
+    if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
         if shooting != Shooting.MULTIPLE:
             with pytest.raises(
                 RuntimeError, match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE"
@@ -741,10 +741,10 @@ def test_integrate_multiphase_non_continuous(shooting, integrator, ode_solver):
     sol_integrated = sol.integrate(**opts)
     shapes = (6, 3, 3)
 
-    decimal = 1 if use_scipy or ode_solver == OdeSolver.COLLOCATION else 8
+    decimal = 1 if integrator != SolutionIntegrator.DEFAULT or ode_solver == OdeSolver.COLLOCATION else 8
     for i in range(len(sol_integrated.states)):
         for k, key in enumerate(sol.states[i]):
-            if not use_scipy or shooting == Shooting.MULTIPLE:
+            if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
                 np.testing.assert_almost_equal(
                     sol_integrated.states[i][key][:, [0, -1]], sol.states[i][key][:, [0, -1]], decimal=decimal
                 )
@@ -752,7 +752,7 @@ def test_integrate_multiphase_non_continuous(shooting, integrator, ode_solver):
                     sol_integrated.states[i][key][:, [0, -2]], sol.states[i][key][:, [0, -1]], decimal=decimal
                 )
 
-            if ode_solver == OdeSolver.COLLOCATION and not use_scipy:
+            if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
                 assert sol_integrated.states[i][key].shape == (shapes[k], n_shooting[i] * (4 + 1) + 1)
             else:
                 assert sol_integrated.states[i][key].shape == (shapes[k], n_shooting[i] * (5 + 1) + 1)
@@ -783,7 +783,6 @@ def test_integrate_multiphase_merged_non_continuous(shooting, integrator, ode_so
 
     sol = ocp.solve()
 
-    use_scipy = integrator.value is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
@@ -806,7 +805,7 @@ def test_integrate_multiphase_merged_non_continuous(shooting, integrator, ode_so
             _ = sol.integrate(**opts)
 
     elif ode_solver == OdeSolver.COLLOCATION:
-        if not use_scipy:
+        if integrator != SolutionIntegrator.DEFAULT:
             with pytest.raises(
                 RuntimeError,
                 match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE",
@@ -820,11 +819,11 @@ def test_integrate_multiphase_merged_non_continuous(shooting, integrator, ode_so
     sol_integrated = sol.integrate(**opts)
     shapes = (6, 3, 3)
 
-    decimal = 0 if use_scipy or ode_solver == OdeSolver.COLLOCATION else 8
-    steps = 4 if not use_scipy and ode_solver == OdeSolver.COLLOCATION else 5
+    decimal = 0 if integrator != SolutionIntegrator.DEFAULT or ode_solver == OdeSolver.COLLOCATION else 8
+    steps = 4 if integrator != SolutionIntegrator.DEFAULT and ode_solver == OdeSolver.COLLOCATION else 5
     for k, key in enumerate(sol.states[0]):
         expected = np.array([sol.states[0][key][:, 0], sol.states[-1][key][:, -1]]).T
-        if not use_scipy or shooting == Shooting.MULTIPLE:
+        if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
             np.testing.assert_almost_equal(sol_integrated.states[key][:, [0, -1]], expected, decimal=decimal)
             np.testing.assert_almost_equal(sol_integrated.states[key][:, [0, -2]], expected, decimal=decimal)
 

From 7a22ca7e8d0990d11df3d362fe66326c5815a426 Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Thu, 11 Nov 2021 18:04:22 -0500
Subject: [PATCH 07/12] fix push

---
 bioptim/optimization/solution.py | 6 +++---
 tests/test_simulate.py           | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/bioptim/optimization/solution.py b/bioptim/optimization/solution.py
index b20a4d67b..32696dda0 100644
--- a/bioptim/optimization/solution.py
+++ b/bioptim/optimization/solution.py
@@ -552,7 +552,7 @@ def __perform_integration(
     ):
         n_direct_collocation = sum([nlp.ode_solver.is_direct_collocation for nlp in self.ocp.nlp])
 
-        if n_direct_collocation > 0 and integrator != SolutionIntegrator.DEFAULT:
+        if n_direct_collocation > 0 and integrator == SolutionIntegrator.DEFAULT:
             if continuous:
                 raise RuntimeError(
                     "Integration with direct collocation must be not continuous if a scipy integrator is used"
@@ -576,7 +576,7 @@ def __perform_integration(
         for p, nlp in enumerate(self.ocp.nlp):
             param_scaling = nlp.parameters.scaling
             n_states = self._states[p]["all"].shape[0]
-            n_steps = nlp.ode_solver.steps_scipy if integrator == SolutionIntegrator.DEFAULT else nlp.ode_solver.steps
+            n_steps = nlp.ode_solver.steps_scipy if integrator != SolutionIntegrator.DEFAULT else nlp.ode_solver.steps
             if not continuous:
                 n_steps += 1
             if keep_intermediate_points:
@@ -599,7 +599,7 @@ def __perform_integration(
             else:
                 col = (
                     slice(0, n_steps)
-                    if nlp.ode_solver.is_direct_collocation and integrator != SolutionIntegrator.DEFAULT
+                    if nlp.ode_solver.is_direct_collocation and integrator == SolutionIntegrator.DEFAULT
                     else 0
                 )
                 x0 = self._states[p]["all"][:, col]
diff --git a/tests/test_simulate.py b/tests/test_simulate.py
index 0849710bb..0705827ce 100644
--- a/tests/test_simulate.py
+++ b/tests/test_simulate.py
@@ -194,7 +194,7 @@ def test_integrate(integrator, ode_solver):
         _ = sol.integrate(**opts)
 
     opts["keep_intermediate_points"] = True
-    if ode_solver == OdeSolver.COLLOCATION and integrator.value is None:
+    if ode_solver == OdeSolver.COLLOCATION and integrator == SolutionIntegrator.DEFAULT:
         with pytest.raises(RuntimeError, match="Integration with direct collocation must be not continuous"):
             sol.integrate(**opts)
         return

From ff7a8bfd23cf040a3a62bc8915553396853113c4 Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Fri, 12 Nov 2021 09:54:22 -0500
Subject: [PATCH 08/12] unmixing the use of scipy or not

---
 bioptim/optimization/solution.py |  2 +-
 tests/test_simulate.py           | 31 ++++++++++++++++---------------
 2 files changed, 17 insertions(+), 16 deletions(-)

diff --git a/bioptim/optimization/solution.py b/bioptim/optimization/solution.py
index 32696dda0..e40b4d364 100644
--- a/bioptim/optimization/solution.py
+++ b/bioptim/optimization/solution.py
@@ -561,7 +561,7 @@ def __perform_integration(
             if shooting_type != Shooting.MULTIPLE:
                 raise RuntimeError(
                     "Integration with direct collocation must using shooting_type=Shooting.MULTIPLE "
-                    "if not a scipy integrator is not used"
+                    "if a scipy integrator is not used"
                 )
 
         # Copy the data
diff --git a/tests/test_simulate.py b/tests/test_simulate.py
index 0705827ce..d90d67712 100644
--- a/tests/test_simulate.py
+++ b/tests/test_simulate.py
@@ -480,7 +480,7 @@ def test_integrate_non_continuous(shooting, merge, integrator, ode_solver):
     if (
         ode_solver == OdeSolver.COLLOCATION
         and shooting != Shooting.MULTIPLE
-        and integrator != SolutionIntegrator.DEFAULT
+        and integrator == SolutionIntegrator.DEFAULT
     ):
         with pytest.raises(
             RuntimeError,
@@ -547,7 +547,7 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, od
         return
 
     if ode_solver == OdeSolver.COLLOCATION:
-        if integrator != SolutionIntegrator.DEFAULT:
+        if integrator == SolutionIntegrator.DEFAULT:
             if shooting != Shooting.MULTIPLE:
                 with pytest.raises(
                     RuntimeError, match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE"
@@ -566,7 +566,7 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, od
             return
 
     opts["continuous"] = True
-    if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
+    if ode_solver == OdeSolver.COLLOCATION and integrator == SolutionIntegrator.DEFAULT:
         with pytest.raises(
             RuntimeError,
             match="Integration with direct collocation must be not continuous",
@@ -580,7 +580,7 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, od
     decimal = 1 if integrator != SolutionIntegrator.DEFAULT else 8
     for i in range(len(sol_integrated.states)):
         for k, key in enumerate(sol.states[i]):
-            if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
+            if integrator == SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
                 np.testing.assert_almost_equal(
                     sol_integrated.states[i][key][:, [0, -1]], sol.states[i][key][:, [0, -1]], decimal=decimal
                 )
@@ -588,7 +588,7 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, od
             if keep_intermediate_points:
                 assert sol_integrated.states[i][key].shape == (shapes[k], n_shooting[i] * 5 + 1)
             else:
-                if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
+                if integrator == SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
                     np.testing.assert_almost_equal(sol_integrated.states[i][key], sol.states[i][key])
                 assert sol_integrated.states[i][key].shape == (shapes[k], n_shooting[i] + 1)
             if ode_solver == OdeSolver.COLLOCATION:
@@ -636,7 +636,7 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, integra
             _ = sol.integrate(**opts)
         return
 
-    if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
+    if ode_solver == OdeSolver.COLLOCATION and integrator == SolutionIntegrator.DEFAULT:
         if shooting != Shooting.MULTIPLE:
             with pytest.raises(
                 RuntimeError, match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE"
@@ -656,7 +656,7 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, integra
 
     opts["merge_phases"] = True
     opts["continuous"] = True
-    if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
+    if ode_solver == OdeSolver.COLLOCATION and integrator == SolutionIntegrator.DEFAULT:
         with pytest.raises(
             RuntimeError,
             match="Integration with direct collocation must be not continuous",
@@ -671,14 +671,14 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, integra
     decimal = 0 if integrator != SolutionIntegrator.DEFAULT else 8
     for k, key in enumerate(sol.states[0]):
         expected = np.array([sol.states[0][key][:, 0], sol.states[-1][key][:, -1]]).T
-        if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
+        if integrator == SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
             np.testing.assert_almost_equal(sol_integrated.states[key][:, [0, -1]], expected, decimal=decimal)
 
         if keep_intermediate_points:
             assert sol_integrated.states[key].shape == (shapes[k], sum(n_shooting) * 5 + 1)
         else:
             # The interpolation prevents from comparing all points
-            if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
+            if integrator == SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
                 expected = np.concatenate(
                     (sol.states[0][key][:, 0:1], sol.states[-1][key][:, -1][:, np.newaxis]), axis=1
                 )
@@ -730,7 +730,7 @@ def test_integrate_multiphase_non_continuous(shooting, integrator, ode_solver):
             _ = sol.integrate(**opts)
         return
 
-    if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
+    if ode_solver == OdeSolver.COLLOCATION and integrator == SolutionIntegrator.DEFAULT:
         if shooting != Shooting.MULTIPLE:
             with pytest.raises(
                 RuntimeError, match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE"
@@ -744,7 +744,7 @@ def test_integrate_multiphase_non_continuous(shooting, integrator, ode_solver):
     decimal = 1 if integrator != SolutionIntegrator.DEFAULT or ode_solver == OdeSolver.COLLOCATION else 8
     for i in range(len(sol_integrated.states)):
         for k, key in enumerate(sol.states[i]):
-            if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
+            if integrator == SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
                 np.testing.assert_almost_equal(
                     sol_integrated.states[i][key][:, [0, -1]], sol.states[i][key][:, [0, -1]], decimal=decimal
                 )
@@ -752,7 +752,7 @@ def test_integrate_multiphase_non_continuous(shooting, integrator, ode_solver):
                     sol_integrated.states[i][key][:, [0, -2]], sol.states[i][key][:, [0, -1]], decimal=decimal
                 )
 
-            if ode_solver == OdeSolver.COLLOCATION and integrator != SolutionIntegrator.DEFAULT:
+            if ode_solver == OdeSolver.COLLOCATION and integrator == SolutionIntegrator.DEFAULT:
                 assert sol_integrated.states[i][key].shape == (shapes[k], n_shooting[i] * (4 + 1) + 1)
             else:
                 assert sol_integrated.states[i][key].shape == (shapes[k], n_shooting[i] * (5 + 1) + 1)
@@ -808,7 +808,8 @@ def test_integrate_multiphase_merged_non_continuous(shooting, integrator, ode_so
         if integrator != SolutionIntegrator.DEFAULT:
             with pytest.raises(
                 RuntimeError,
-                match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE",
+                match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE"
+                "if a scipy integrator is not used",
             ):
                 _ = sol.integrate(**opts)
             return
@@ -820,10 +821,10 @@ def test_integrate_multiphase_merged_non_continuous(shooting, integrator, ode_so
     shapes = (6, 3, 3)
 
     decimal = 0 if integrator != SolutionIntegrator.DEFAULT or ode_solver == OdeSolver.COLLOCATION else 8
-    steps = 4 if integrator != SolutionIntegrator.DEFAULT and ode_solver == OdeSolver.COLLOCATION else 5
+    steps = 4 if integrator == SolutionIntegrator.DEFAULT and ode_solver == OdeSolver.COLLOCATION else 5
     for k, key in enumerate(sol.states[0]):
         expected = np.array([sol.states[0][key][:, 0], sol.states[-1][key][:, -1]]).T
-        if integrator != SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
+        if integrator == SolutionIntegrator.DEFAULT or shooting == Shooting.MULTIPLE:
             np.testing.assert_almost_equal(sol_integrated.states[key][:, [0, -1]], expected, decimal=decimal)
             np.testing.assert_almost_equal(sol_integrated.states[key][:, [0, -2]], expected, decimal=decimal)
 

From 7edd14117fafca3094cdaa06781bc9cc08791d85 Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Fri, 12 Nov 2021 09:56:33 -0500
Subject: [PATCH 09/12] fixing a last test.

---
 tests/test_simulate.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/test_simulate.py b/tests/test_simulate.py
index d90d67712..6619c20b6 100644
--- a/tests/test_simulate.py
+++ b/tests/test_simulate.py
@@ -805,10 +805,10 @@ def test_integrate_multiphase_merged_non_continuous(shooting, integrator, ode_so
             _ = sol.integrate(**opts)
 
     elif ode_solver == OdeSolver.COLLOCATION:
-        if integrator != SolutionIntegrator.DEFAULT:
+        if integrator == SolutionIntegrator.DEFAULT:
             with pytest.raises(
                 RuntimeError,
-                match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE"
+                match="Integration with direct collocation must using shooting_type=Shooting.MULTIPLE "
                 "if a scipy integrator is not used",
             ):
                 _ = sol.integrate(**opts)

From dd15a5a51fc5f0e1bd6d60dacfcc960d8f91b731 Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Fri, 12 Nov 2021 10:00:32 -0500
Subject: [PATCH 10/12] Default solution integrator as an empty string

---
 bioptim/misc/enums.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/bioptim/misc/enums.py b/bioptim/misc/enums.py
index 22cd0e2fc..cb03658ea 100644
--- a/bioptim/misc/enums.py
+++ b/bioptim/misc/enums.py
@@ -110,7 +110,7 @@ class SolutionIntegrator(Enum):
     Selection of integrator to use integrate function
     """
 
-    DEFAULT = None
+    DEFAULT = ""
     SCIPY_RK23 = "RK23"
     SCIPY_RK45 = "RK45"
     SCIPY_DOP853 = "DOP853"

From b6fe295a4158a08eb25fab8da17e63888736a270 Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Fri, 12 Nov 2021 10:46:48 -0500
Subject: [PATCH 11/12] ready to merge

---
 bioptim/optimization/solution.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/bioptim/optimization/solution.py b/bioptim/optimization/solution.py
index e40b4d364..2ba27127b 100644
--- a/bioptim/optimization/solution.py
+++ b/bioptim/optimization/solution.py
@@ -623,8 +623,6 @@ def __perform_integration(
                         x0,
                         t_eval=t_eval,
                         method=integrator.value,
-                        atol=1e-10,
-                        rtol=1e-10,
                     ).y
 
                     next_state_col = (

From e16441c257cb55f0ecef20ecbb9919657933916a Mon Sep 17 00:00:00 2001
From: Ipuch <puchaud.pierre@gmail.com>
Date: Fri, 12 Nov 2021 11:18:22 -0500
Subject: [PATCH 12/12] non empty str for SolutionIntegrator.DEFAULT

---
 bioptim/misc/enums.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/bioptim/misc/enums.py b/bioptim/misc/enums.py
index cb03658ea..41143d875 100644
--- a/bioptim/misc/enums.py
+++ b/bioptim/misc/enums.py
@@ -110,7 +110,7 @@ class SolutionIntegrator(Enum):
     Selection of integrator to use integrate function
     """
 
-    DEFAULT = ""
+    DEFAULT = "DEFAULT"
     SCIPY_RK23 = "RK23"
     SCIPY_RK45 = "RK45"
     SCIPY_DOP853 = "DOP853"