From e1a176b59afcf6fb8c1c912bf3b2d3a28437405f Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Wed, 28 Feb 2024 18:16:23 -0500 Subject: [PATCH 01/18] for integrate continuity --- bioptim/optimization/solution/solution.py | 26 +++++++++++++++-------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index 08b53f1dc..ea0f73fb0 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -230,14 +230,14 @@ def from_initial_guess(cls, ocp, sol: list): "an InitialGuess[List] of len 4 (states, controls, parameters, algebraic_states), " "or a None" ) - if sum([len(s) != len(all_ns) if p != 4 else False for p, s in enumerate(sol)]) != 0: - raise ValueError("The InitialGuessList len must match the number of phases") - if n_param != 0: - if len(sol) != 3 and len(sol[3]) != 1 and sol[3][0].shape != (n_param, 1): - raise ValueError( - "The 3rd element is the InitialGuess of the parameter and " - "should be a unique vector of size equal to n_param" - ) + # if sum([len(s) != len(all_ns) if p != 4 else False for p, s in enumerate(sol)]) != 0: + # raise ValueError("The InitialGuessList len must match the number of phases") + # if n_param != 0: + # if len(sol) != 3 and len(sol[3]) != 1 and sol[3][0].shape != (n_param, 1): + # raise ValueError( + # "The 3rd element is the InitialGuess of the parameter and " + # "should be a unique vector of size equal to n_param" + # ) dt, sol_states, sol_controls, sol_params, sol_algebraic_states = sol @@ -281,7 +281,9 @@ def from_initial_guess(cls, ocp, sol: list): # For parameters if n_param: for p, ss in enumerate(sol_params): - vector = np.concatenate((vector, np.repeat(ss.init, all_ns[p] + 1)[:, np.newaxis])) + for key in ss.keys(): + # vector = np.concatenate((vector, np.repeat(ss[key].init, all_ns[p] + 1)[:, np.newaxis])) + vector = np.concatenate((vector, np.repeat(ss[key].init, 1)[:, np.newaxis])) # For algebraic_states variables for p, ss in enumerate(sol_algebraic_states): @@ -734,6 +736,12 @@ def integrate( for ns, sol_ns in enumerate(integrated_sol): out[p][key][ns] = sol_ns[nlp.states[key].index, :] + if shooting_type == Shooting.SINGLE and p < len(self.ocp.nlp) - 1: + state_list = [] + for key in nlp.states.keys(): + state_list.append(out[p][key][-1][0]) + x[p+1][0] = np.array(state_list) + if to_merge: out = SolutionData.from_unscaled(self.ocp, out, "x").to_dict(to_merge=to_merge, scaled=False) From 81ee186a2b723e508d573f093ec94f045c9f9ea1 Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Fri, 1 Mar 2024 15:45:18 -0500 Subject: [PATCH 02/18] integration correction --- bioptim/optimization/solution/solution.py | 50 +++++++++++++++++++---- 1 file changed, 41 insertions(+), 9 deletions(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index ea0f73fb0..d3f25af7e 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -232,6 +232,12 @@ def from_initial_guess(cls, ocp, sol: list): ) # if sum([len(s) != len(all_ns) if p != 4 else False for p, s in enumerate(sol)]) != 0: # raise ValueError("The InitialGuessList len must match the number of phases") + + # p!=3 : pour les paramètres + + # TODO : allow empty num + # Donner autant d'initial guess vide que de phase + # if n_param != 0: # if len(sol) != 3 and len(sol[3]) != 1 and sol[3][0].shape != (n_param, 1): # raise ValueError( @@ -688,6 +694,7 @@ def integrate( shooting_type: Shooting = Shooting.SINGLE, integrator: SolutionIntegrator = SolutionIntegrator.OCP, to_merge: SolutionMerge | list[SolutionMerge] = None, + return_time: bool = False, ): has_direct_collocation = sum([nlp.ode_solver.is_direct_collocation for nlp in self.ocp.nlp]) > 0 if has_direct_collocation and integrator == SolutionIntegrator.OCP: @@ -715,11 +722,14 @@ def integrate( u = self._stepwise_controls.to_dict(to_merge=SolutionMerge.KEYS, scaled=False) a = self._decision_algebraic_states.to_dict(to_merge=SolutionMerge.KEYS, scaled=False) + time_vector = [] out: list = [None] * len(self.ocp.nlp) integrated_sol = None for p, nlp in enumerate(self.ocp.nlp): next_x = self._states_for_phase_integration(shooting_type, p, integrated_sol, x, u, params, a) - integrated_sol = solve_ivp_interface( + + # TODO: Get the actual time vector + phase_times, integrated_sol = solve_ivp_interface( shooting_type=shooting_type, nlp=nlp, t=t_spans[p], @@ -730,22 +740,32 @@ def integrate( method=integrator, ) + time_vector.append(phase_times) out[p] = {} for key in nlp.states.keys(): out[p][key] = [None] * nlp.n_states_nodes for ns, sol_ns in enumerate(integrated_sol): out[p][key][ns] = sol_ns[nlp.states[key].index, :] - if shooting_type == Shooting.SINGLE and p < len(self.ocp.nlp) - 1: - state_list = [] - for key in nlp.states.keys(): - state_list.append(out[p][key][-1][0]) - x[p+1][0] = np.array(state_list) - if to_merge: out = SolutionData.from_unscaled(self.ocp, out, "x").to_dict(to_merge=to_merge, scaled=False) - return out if len(out) > 1 else out[0] + # Make a _remove_integrated_duplicates(out) function + if ( + shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) + and SolutionMerge.NODES in to_merge + ): + for p in range(len(out)): + _, index_duplicate = np.unique(time_vector[p]) + # Remove the duplicate for both merged keys and non-merged keys + # Something like out[p][key] = out[p][key][index_duplicate] + + # Also deal with the Discontinuous and non-discontinuous phases if phase is merged + + if return_time: + return out if len(out) > 1 else out[0], time_vector if len(time_vector) > 1 else time_vector[0] + else: + return out if len(out) > 1 else out[0] def _states_for_phase_integration( self, @@ -791,6 +811,15 @@ def _states_for_phase_integration( if phase_idx == 0: return [decision_states[phase_idx][0]] + # if shooting_type == Shooting.SINGLE and p < len(self.ocp.nlp) - 1: + # state_list = [] + # for key in nlp.states.keys(): + # state_list.append(out[p][key][-1][0]) + # x[p+1][0] = np.array(state_list) + + # TODO A mettre ici + + penalty = self.ocp.phase_transitions[phase_idx - 1] t0 = PenaltyHelpers.t0(penalty, 0, lambda p, n: self._stepwise_times[p][n][0]) @@ -799,6 +828,9 @@ def _states_for_phase_integration( # based on the phase transition objective or constraint function. That is why we need to concatenate # twice the last state x = PenaltyHelpers.states(penalty, 0, lambda p, n, sn: integrated_states[-1]) + + # x = decision_phase x - 1 + u = PenaltyHelpers.controls( penalty, 0, @@ -820,7 +852,7 @@ def _states_for_phase_integration( f"please integrate with Shooting.SINGLE_DISCONTINUOUS_PHASE." ) - return [decision_states[phase_idx][0] + dx] + return [(integrated_states[-1] if shooting_type == Shooting.SINGLE else decision_states[phase_idx]) + dx] def _integrate_stepwise(self) -> None: """ From a55b88889233f6c1585863f16806de6283906569 Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Mon, 4 Mar 2024 17:08:40 -0500 Subject: [PATCH 03/18] PR : integrate shooting single and delete duplicates of merged results --- bioptim/interfaces/solve_ivp_interface.py | 4 +- bioptim/optimization/solution/solution.py | 83 +++++++++++++---------- 2 files changed, 50 insertions(+), 37 deletions(-) diff --git a/bioptim/interfaces/solve_ivp_interface.py b/bioptim/interfaces/solve_ivp_interface.py index 5c254fc44..935dedca6 100644 --- a/bioptim/interfaces/solve_ivp_interface.py +++ b/bioptim/interfaces/solve_ivp_interface.py @@ -48,6 +48,7 @@ def solve_ivp_interface( y = [] control_type = nlp.control_type + integration_time_vector = [] for node in range(nlp.ns): if method == SolutionIntegrator.OCP: @@ -55,6 +56,7 @@ def solve_ivp_interface( else: t_span = t[node] t_eval = np.linspace(float(t_span[0]), float(t_span[1]), nlp.n_states_stepwise_steps(node)) + integration_time_vector.append(t_eval[0]) # If multiple shooting, we need to set the first x0, otherwise use the previous answer x0i = np.array(x[node] if node == 0 or shooting_type == Shooting.MULTIPLE else y[-1][:, -1]) @@ -93,7 +95,7 @@ def solve_ivp_interface( y.append(x[-1] if shooting_type == Shooting.MULTIPLE else y[-1][:, -1][:, np.newaxis]) - return y + return integration_time_vector, y def _solve_ivp_scipy_interface( diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index d3f25af7e..f468924f6 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -230,20 +230,18 @@ def from_initial_guess(cls, ocp, sol: list): "an InitialGuess[List] of len 4 (states, controls, parameters, algebraic_states), " "or a None" ) - # if sum([len(s) != len(all_ns) if p != 4 else False for p, s in enumerate(sol)]) != 0: - # raise ValueError("The InitialGuessList len must match the number of phases") - # p!=3 : pour les paramètres + if len(sol[0]) != len(all_ns): + raise ValueError("The time step dt array len must match the number of phases") + if sum([len(s) != len(all_ns) if p != 3 and len(sol[p + 1].keys()) != 0 else False for p, s in enumerate(sol[:1])]) != 0: + raise ValueError("The InitialGuessList len must match the number of phases") - # TODO : allow empty num - # Donner autant d'initial guess vide que de phase - - # if n_param != 0: - # if len(sol) != 3 and len(sol[3]) != 1 and sol[3][0].shape != (n_param, 1): - # raise ValueError( - # "The 3rd element is the InitialGuess of the parameter and " - # "should be a unique vector of size equal to n_param" - # ) + if n_param != 0: + if len(sol) != 3 and len(sol[3]) != 1 and sol[3][0].shape != (n_param, 1): + raise ValueError( + "The 3rd element is the InitialGuess of the parameter and " + "should be a unique vector of size equal to n_param" + ) dt, sol_states, sol_controls, sol_params, sol_algebraic_states = sol @@ -290,6 +288,8 @@ def from_initial_guess(cls, ocp, sol: list): for key in ss.keys(): # vector = np.concatenate((vector, np.repeat(ss[key].init, all_ns[p] + 1)[:, np.newaxis])) vector = np.concatenate((vector, np.repeat(ss[key].init, 1)[:, np.newaxis])) + # TODO : ask pariterre about this modification (the previous version did not enable the parameters + # to be taken into consideration for each phase) # For algebraic_states variables for p, ss in enumerate(sol_algebraic_states): @@ -728,7 +728,6 @@ def integrate( for p, nlp in enumerate(self.ocp.nlp): next_x = self._states_for_phase_integration(shooting_type, p, integrated_sol, x, u, params, a) - # TODO: Get the actual time vector phase_times, integrated_sol = solve_ivp_interface( shooting_type=shooting_type, nlp=nlp, @@ -749,18 +748,7 @@ def integrate( if to_merge: out = SolutionData.from_unscaled(self.ocp, out, "x").to_dict(to_merge=to_merge, scaled=False) - - # Make a _remove_integrated_duplicates(out) function - if ( - shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) - and SolutionMerge.NODES in to_merge - ): - for p in range(len(out)): - _, index_duplicate = np.unique(time_vector[p]) - # Remove the duplicate for both merged keys and non-merged keys - # Something like out[p][key] = out[p][key][index_duplicate] - - # Also deal with the Discontinuous and non-discontinuous phases if phase is merged + time_vector, out = self._remove_integrated_duplicates(time_vector, out, shooting_type, to_merge) if return_time: return out if len(out) > 1 else out[0], time_vector if len(time_vector) > 1 else time_vector[0] @@ -811,15 +799,6 @@ def _states_for_phase_integration( if phase_idx == 0: return [decision_states[phase_idx][0]] - # if shooting_type == Shooting.SINGLE and p < len(self.ocp.nlp) - 1: - # state_list = [] - # for key in nlp.states.keys(): - # state_list.append(out[p][key][-1][0]) - # x[p+1][0] = np.array(state_list) - - # TODO A mettre ici - - penalty = self.ocp.phase_transitions[phase_idx - 1] t0 = PenaltyHelpers.t0(penalty, 0, lambda p, n: self._stepwise_times[p][n][0]) @@ -829,8 +808,6 @@ def _states_for_phase_integration( # twice the last state x = PenaltyHelpers.states(penalty, 0, lambda p, n, sn: integrated_states[-1]) - # x = decision_phase x - 1 - u = PenaltyHelpers.controls( penalty, 0, @@ -854,6 +831,40 @@ def _states_for_phase_integration( return [(integrated_states[-1] if shooting_type == Shooting.SINGLE else decision_states[phase_idx]) + dx] + def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merge): + if ( + shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) + and SolutionMerge.NODES in to_merge + ): + if SolutionMerge.PHASES in to_merge: + redundant_index = [] + merged_out = {} + for i in range(len(self.ocp.nlp)): + redundant_index.append( + (np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int) + + redundant_index[-1][-1]).tolist() if i != 0 else np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, + self.ocp.nlp[0].ns + 1, dtype=int).tolist()) + redundant_index = np.array(redundant_index).flatten().tolist() + for key in out.keys(): + merged_out[key] = np.delete(out[key][0], redundant_index[:-1])[:-1] + time_vector = [i for time_vector_sub in time_vector for i in time_vector_sub] + + else: + merged_out = [] + for i in range(len(out)): + redundant_index = np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int).tolist() + redundant_index = [0] + redundant_index[:-1] + out_per_phase = {} + for key in out[i].keys(): + out_per_phase[key] = np.delete(out[i][key][0], redundant_index) + merged_out.append(out_per_phase) + + else: + merged_out = out + + return time_vector, merged_out + + def _integrate_stepwise(self) -> None: """ This method integrate to stepwise level the states. That is the states that are used in the dynamics and From 305ff91957c38c773e2f018264aab302d6b782d6 Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Mon, 4 Mar 2024 17:23:52 -0500 Subject: [PATCH 04/18] fixing dimension + black --- bioptim/optimization/solution/solution.py | 34 +++++++++++++++-------- 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index f468924f6..ee14155cd 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -233,7 +233,15 @@ def from_initial_guess(cls, ocp, sol: list): if len(sol[0]) != len(all_ns): raise ValueError("The time step dt array len must match the number of phases") - if sum([len(s) != len(all_ns) if p != 3 and len(sol[p + 1].keys()) != 0 else False for p, s in enumerate(sol[:1])]) != 0: + if ( + sum( + [ + len(s) != len(all_ns) if p != 3 and len(sol[p + 1].keys()) != 0 else False + for p, s in enumerate(sol[:1]) + ] + ) + != 0 + ): raise ValueError("The InitialGuessList len must match the number of phases") if n_param != 0: @@ -832,27 +840,30 @@ def _states_for_phase_integration( return [(integrated_states[-1] if shooting_type == Shooting.SINGLE else decision_states[phase_idx]) + dx] def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merge): - if ( - shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) - and SolutionMerge.NODES in to_merge - ): + if shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) and SolutionMerge.NODES in to_merge: if SolutionMerge.PHASES in to_merge: redundant_index = [] merged_out = {} for i in range(len(self.ocp.nlp)): redundant_index.append( - (np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int) + - redundant_index[-1][-1]).tolist() if i != 0 else np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, - self.ocp.nlp[0].ns + 1, dtype=int).tolist()) + ( + np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int) + + redundant_index[-1][-1] + ).tolist() + if i != 0 + else np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int).tolist() + ) redundant_index = np.array(redundant_index).flatten().tolist() for key in out.keys(): - merged_out[key] = np.delete(out[key][0], redundant_index[:-1])[:-1] + merged_out[key] = np.delete(out[key][0], redundant_index[:-1])[:-1][np.newaxis, :] time_vector = [i for time_vector_sub in time_vector for i in time_vector_sub] else: merged_out = [] for i in range(len(out)): - redundant_index = np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int).tolist() + redundant_index = np.linspace( + 1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int + ).tolist() redundant_index = [0] + redundant_index[:-1] out_per_phase = {} for key in out[i].keys(): @@ -864,7 +875,6 @@ def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merg return time_vector, merged_out - def _integrate_stepwise(self) -> None: """ This method integrate to stepwise level the states. That is the states that are used in the dynamics and @@ -886,7 +896,7 @@ def _integrate_stepwise(self) -> None: unscaled: list = [None] * len(self.ocp.nlp) for p, nlp in enumerate(self.ocp.nlp): - integrated_sol = solve_ivp_interface( + phase_times, integrated_sol = solve_ivp_interface( shooting_type=Shooting.MULTIPLE, nlp=nlp, t=t_spans[p], From 34ceabcec2d31d94f44d723d1e738aea6faae6e8 Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Tue, 5 Mar 2024 17:55:30 -0500 Subject: [PATCH 05/18] changing redundent node value detection (not working yet) --- bioptim/interfaces/solve_ivp_interface.py | 3 +- bioptim/optimization/solution/solution.py | 35 +++++++++++++---------- 2 files changed, 22 insertions(+), 16 deletions(-) diff --git a/bioptim/interfaces/solve_ivp_interface.py b/bioptim/interfaces/solve_ivp_interface.py index 935dedca6..3418c568d 100644 --- a/bioptim/interfaces/solve_ivp_interface.py +++ b/bioptim/interfaces/solve_ivp_interface.py @@ -56,7 +56,7 @@ def solve_ivp_interface( else: t_span = t[node] t_eval = np.linspace(float(t_span[0]), float(t_span[1]), nlp.n_states_stepwise_steps(node)) - integration_time_vector.append(t_eval[0]) + integration_time_vector.append(np.linspace(float(t[node][0]), float(t[node][1]), nlp.n_states_stepwise_steps(node))) # If multiple shooting, we need to set the first x0, otherwise use the previous answer x0i = np.array(x[node] if node == 0 or shooting_type == Shooting.MULTIPLE else y[-1][:, -1]) @@ -94,6 +94,7 @@ def solve_ivp_interface( y.append(result) y.append(x[-1] if shooting_type == Shooting.MULTIPLE else y[-1][:, -1][:, np.newaxis]) + integration_time_vector = np.concatenate(integration_time_vector, axis=0) return integration_time_vector, y diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index ee14155cd..e780002cb 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -840,22 +840,24 @@ def _states_for_phase_integration( return [(integrated_states[-1] if shooting_type == Shooting.SINGLE else decision_states[phase_idx]) + dx] def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merge): - if shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) and SolutionMerge.NODES in to_merge: - if SolutionMerge.PHASES in to_merge: - redundant_index = [] + if shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) and (to_merge == SolutionMerge.NODES or SolutionMerge.NODES in to_merge): + if (isinstance(to_merge, list) and SolutionMerge.PHASES in to_merge) or to_merge == SolutionMerge.PHASES : merged_out = {} - for i in range(len(self.ocp.nlp)): - redundant_index.append( - ( - np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int) - + redundant_index[-1][-1] - ).tolist() - if i != 0 - else np.linspace(1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int).tolist() - ) - redundant_index = np.array(redundant_index).flatten().tolist() + + redundant_index_between_phases = np.cumsum([len(time_vector[i]) + 1 for i in range(len(self.ocp.nlp))]).tolist() + time_vector = [i for time_vector_sub in time_vector for i in time_vector_sub] if len(time_vector) != 1 else time_vector [0] + unique_time_vector, unique_index = np.unique(time_vector, return_index=True) + redundant_index = [i for i in range(len(time_vector)) if i not in unique_index] + for key in out.keys(): - merged_out[key] = np.delete(out[key][0], redundant_index[:-1])[:-1][np.newaxis, :] + merged_out[key] = [None] * len(out[key]) + for i in range(len(out[key])): + temp_merged_out = np.delete(out[key][i], redundant_index_between_phases[:-1])[:-1][np.newaxis, :] + merged_out[key][i] = np.delete(temp_merged_out, redundant_index[:-1])[:-1][np.newaxis, :] + + + merged_out[key] = np.concatenate(merged_out[key], axis=0) + # merged_out[key] = np.delete(out[key][0], redundant_index[:-1])[:-1][np.newaxis, :] time_vector = [i for time_vector_sub in time_vector for i in time_vector_sub] else: @@ -867,7 +869,10 @@ def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merg redundant_index = [0] + redundant_index[:-1] out_per_phase = {} for key in out[i].keys(): - out_per_phase[key] = np.delete(out[i][key][0], redundant_index) + out_per_phase[key] = [None] * len(out[i][key]) + for j in range(len(out[i][key])): + out_per_phase[key][j] = np.delete(out[i][key][j], redundant_index[i])[:-1][np.newaxis, :] + out_per_phase[key] = np.concatenate(out_per_phase[key], axis=0) merged_out.append(out_per_phase) else: From 3ec3cabedbcd3ef05b15493cbe8477cfcf1e2f72 Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Wed, 6 Mar 2024 13:35:05 -0500 Subject: [PATCH 06/18] new redundant node detection for merge pahse/node (need one for node only) --- bioptim/optimization/solution/solution.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index e780002cb..380ee6d3e 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -853,12 +853,10 @@ def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merg merged_out[key] = [None] * len(out[key]) for i in range(len(out[key])): temp_merged_out = np.delete(out[key][i], redundant_index_between_phases[:-1])[:-1][np.newaxis, :] - merged_out[key][i] = np.delete(temp_merged_out, redundant_index[:-1])[:-1][np.newaxis, :] - + merged_out[key][i] = np.delete(temp_merged_out, redundant_index)[np.newaxis, :] merged_out[key] = np.concatenate(merged_out[key], axis=0) - # merged_out[key] = np.delete(out[key][0], redundant_index[:-1])[:-1][np.newaxis, :] - time_vector = [i for time_vector_sub in time_vector for i in time_vector_sub] + time_vector = unique_time_vector else: merged_out = [] From b628a705989878527603fb2f5b3906dd1c153907 Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Wed, 6 Mar 2024 16:06:04 -0500 Subject: [PATCH 07/18] enabled only to_merge = merge.node --- bioptim/optimization/solution/solution.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index 380ee6d3e..cee88845a 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -837,7 +837,7 @@ def _states_for_phase_integration( f"please integrate with Shooting.SINGLE_DISCONTINUOUS_PHASE." ) - return [(integrated_states[-1] if shooting_type == Shooting.SINGLE else decision_states[phase_idx]) + dx] + return [(integrated_states[-1] if shooting_type == Shooting.SINGLE else decision_states[phase_idx][0]) + dx] def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merge): if shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) and (to_merge == SolutionMerge.NODES or SolutionMerge.NODES in to_merge): @@ -854,22 +854,23 @@ def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merg for i in range(len(out[key])): temp_merged_out = np.delete(out[key][i], redundant_index_between_phases[:-1])[:-1][np.newaxis, :] merged_out[key][i] = np.delete(temp_merged_out, redundant_index)[np.newaxis, :] - merged_out[key] = np.concatenate(merged_out[key], axis=0) time_vector = unique_time_vector else: merged_out = [] for i in range(len(out)): - redundant_index = np.linspace( - 1, self.ocp.nlp[0].ns * 2 + 1, self.ocp.nlp[0].ns + 1, dtype=int - ).tolist() - redundant_index = [0] + redundant_index[:-1] + phase_time_vector = time_vector[i] + unique_time_vector, unique_index = np.unique(phase_time_vector, return_index=True) + redundant_index = [i for i in range(len(phase_time_vector)) if i not in unique_index] + out_per_phase = {} for key in out[i].keys(): out_per_phase[key] = [None] * len(out[i][key]) for j in range(len(out[i][key])): - out_per_phase[key][j] = np.delete(out[i][key][j], redundant_index[i])[:-1][np.newaxis, :] + temp_merged_out = np.delete(out[i][key][j], redundant_index)[:-1][np.newaxis, :] + temp_merged_out = temp_merged_out[0][1:] if i > 0 else temp_merged_out + out_per_phase[key][j] = temp_merged_out out_per_phase[key] = np.concatenate(out_per_phase[key], axis=0) merged_out.append(out_per_phase) From 5496db587c0634a09abd7e0a908e0c6dec8e8f9a Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Wed, 6 Mar 2024 16:30:15 -0500 Subject: [PATCH 08/18] black --- bioptim/interfaces/solve_ivp_interface.py | 4 +++- bioptim/optimization/solution/solution.py | 21 +++++++++++++++------ 2 files changed, 18 insertions(+), 7 deletions(-) diff --git a/bioptim/interfaces/solve_ivp_interface.py b/bioptim/interfaces/solve_ivp_interface.py index 3418c568d..1c50abc3f 100644 --- a/bioptim/interfaces/solve_ivp_interface.py +++ b/bioptim/interfaces/solve_ivp_interface.py @@ -56,7 +56,9 @@ def solve_ivp_interface( else: t_span = t[node] t_eval = np.linspace(float(t_span[0]), float(t_span[1]), nlp.n_states_stepwise_steps(node)) - integration_time_vector.append(np.linspace(float(t[node][0]), float(t[node][1]), nlp.n_states_stepwise_steps(node))) + integration_time_vector.append( + np.linspace(float(t[node][0]), float(t[node][1]), nlp.n_states_stepwise_steps(node)) + ) # If multiple shooting, we need to set the first x0, otherwise use the previous answer x0i = np.array(x[node] if node == 0 or shooting_type == Shooting.MULTIPLE else y[-1][:, -1]) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index cee88845a..439664875 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -840,19 +840,29 @@ def _states_for_phase_integration( return [(integrated_states[-1] if shooting_type == Shooting.SINGLE else decision_states[phase_idx][0]) + dx] def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merge): - if shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) and (to_merge == SolutionMerge.NODES or SolutionMerge.NODES in to_merge): - if (isinstance(to_merge, list) and SolutionMerge.PHASES in to_merge) or to_merge == SolutionMerge.PHASES : + if shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) and ( + to_merge == SolutionMerge.NODES or SolutionMerge.NODES in to_merge + ): + if (isinstance(to_merge, list) and SolutionMerge.PHASES in to_merge) or to_merge == SolutionMerge.PHASES: merged_out = {} - redundant_index_between_phases = np.cumsum([len(time_vector[i]) + 1 for i in range(len(self.ocp.nlp))]).tolist() - time_vector = [i for time_vector_sub in time_vector for i in time_vector_sub] if len(time_vector) != 1 else time_vector [0] + redundant_index_between_phases = np.cumsum( + [len(time_vector[i]) + 1 for i in range(len(self.ocp.nlp))] + ).tolist() + time_vector = ( + [i for time_vector_sub in time_vector for i in time_vector_sub] + if len(time_vector) != 1 + else time_vector[0] + ) unique_time_vector, unique_index = np.unique(time_vector, return_index=True) redundant_index = [i for i in range(len(time_vector)) if i not in unique_index] for key in out.keys(): merged_out[key] = [None] * len(out[key]) for i in range(len(out[key])): - temp_merged_out = np.delete(out[key][i], redundant_index_between_phases[:-1])[:-1][np.newaxis, :] + temp_merged_out = np.delete(out[key][i], redundant_index_between_phases[:-1])[:-1][ + np.newaxis, : + ] merged_out[key][i] = np.delete(temp_merged_out, redundant_index)[np.newaxis, :] merged_out[key] = np.concatenate(merged_out[key], axis=0) time_vector = unique_time_vector @@ -871,7 +881,6 @@ def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merg temp_merged_out = np.delete(out[i][key][j], redundant_index)[:-1][np.newaxis, :] temp_merged_out = temp_merged_out[0][1:] if i > 0 else temp_merged_out out_per_phase[key][j] = temp_merged_out - out_per_phase[key] = np.concatenate(out_per_phase[key], axis=0) merged_out.append(out_per_phase) else: From 58ed59a9fc7fac95baca9eaea2a5978639d3c1fc Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Thu, 7 Mar 2024 15:42:14 -0500 Subject: [PATCH 09/18] round the time vector because float(DM) can round the two same DM value differently at e10-8~ --- bioptim/optimization/solution/solution.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index 439664875..06ed79613 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -854,7 +854,10 @@ def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merg if len(time_vector) != 1 else time_vector[0] ) + + time_vector = np.array(time_vector).round(decimals=6).tolist() unique_time_vector, unique_index = np.unique(time_vector, return_index=True) + redundant_index = [i for i in range(len(time_vector)) if i not in unique_index] for key in out.keys(): From 2369cc3e4be3c9a3feb868c1563bf689d2864286 Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Mon, 11 Mar 2024 18:46:16 -0400 Subject: [PATCH 10/18] Applying comments from PR review --- bioptim/interfaces/solve_ivp_interface.py | 7 +- bioptim/optimization/solution/solution.py | 146 ++++++++++++---------- 2 files changed, 79 insertions(+), 74 deletions(-) diff --git a/bioptim/interfaces/solve_ivp_interface.py b/bioptim/interfaces/solve_ivp_interface.py index 1c50abc3f..5c254fc44 100644 --- a/bioptim/interfaces/solve_ivp_interface.py +++ b/bioptim/interfaces/solve_ivp_interface.py @@ -48,7 +48,6 @@ def solve_ivp_interface( y = [] control_type = nlp.control_type - integration_time_vector = [] for node in range(nlp.ns): if method == SolutionIntegrator.OCP: @@ -56,9 +55,6 @@ def solve_ivp_interface( else: t_span = t[node] t_eval = np.linspace(float(t_span[0]), float(t_span[1]), nlp.n_states_stepwise_steps(node)) - integration_time_vector.append( - np.linspace(float(t[node][0]), float(t[node][1]), nlp.n_states_stepwise_steps(node)) - ) # If multiple shooting, we need to set the first x0, otherwise use the previous answer x0i = np.array(x[node] if node == 0 or shooting_type == Shooting.MULTIPLE else y[-1][:, -1]) @@ -96,9 +92,8 @@ def solve_ivp_interface( y.append(result) y.append(x[-1] if shooting_type == Shooting.MULTIPLE else y[-1][:, -1][:, np.newaxis]) - integration_time_vector = np.concatenate(integration_time_vector, axis=0) - return integration_time_vector, y + return y def _solve_ivp_scipy_interface( diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index 06ed79613..fef66baf7 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -233,15 +233,11 @@ def from_initial_guess(cls, ocp, sol: list): if len(sol[0]) != len(all_ns): raise ValueError("The time step dt array len must match the number of phases") - if ( - sum( - [ - len(s) != len(all_ns) if p != 3 and len(sol[p + 1].keys()) != 0 else False - for p, s in enumerate(sol[:1]) - ] - ) - != 0 - ): + + is_right_size = [len(s) != len(all_ns) if p != 3 and len(sol[p + 1].keys()) != 0 else False + for p, s in enumerate(sol[:1])] + + if sum(is_right_size) != 0: raise ValueError("The InitialGuessList len must match the number of phases") if n_param != 0: @@ -294,10 +290,7 @@ def from_initial_guess(cls, ocp, sol: list): if n_param: for p, ss in enumerate(sol_params): for key in ss.keys(): - # vector = np.concatenate((vector, np.repeat(ss[key].init, all_ns[p] + 1)[:, np.newaxis])) vector = np.concatenate((vector, np.repeat(ss[key].init, 1)[:, np.newaxis])) - # TODO : ask pariterre about this modification (the previous version did not enable the parameters - # to be taken into consideration for each phase) # For algebraic_states variables for p, ss in enumerate(sol_algebraic_states): @@ -394,6 +387,7 @@ def stepwise_time( to_merge: SolutionMerge | list[SolutionMerge] = None, time_alignment: TimeAlignment = TimeAlignment.STATES, continuous: bool = True, + duplicated_times: bool = True, ) -> list | np.ndarray: """ Returns the time vector at each node that matches stepwise_states or stepwise_controls @@ -410,6 +404,9 @@ def stepwise_time( continuous: bool If the time should be continuous throughout the whole ocp. If False, then the time is reset at the beginning of each phase. + duplicated_times: bool + If the times should be duplicated for each nodes. + If False, then the returned time vector will not have any duplicated times Returns ------- @@ -421,6 +418,7 @@ def stepwise_time( to_merge=to_merge, time_alignment=time_alignment, continuous=continuous, + duplicated_times=duplicated_times, ) def _process_time_vector( @@ -429,6 +427,7 @@ def _process_time_vector( to_merge: SolutionMerge | list[SolutionMerge], time_alignment: TimeAlignment, continuous: bool, + duplicated_times: bool = True, ): if to_merge is None or isinstance(to_merge, SolutionMerge): to_merge = [to_merge] @@ -483,6 +482,14 @@ def _process_time_vector( else: raise ValueError("time_alignment should be either TimeAlignment.STATES or TimeAlignment.CONTROLS") + if not duplicated_times: + for i in range(len(times)): + for j in range(len(times[i])): + keep_condition = times[i][j].shape[0] == 1 and i == len(times) - 1 + times[i][j] = times[i][j][:] if keep_condition else times[i][j][:-1] + if j == len(times[i]) - 1 and i != len(times) - 1: + del times[i][j] + if continuous: for phase_idx, phase_time in enumerate(times): if phase_idx == 0: @@ -702,6 +709,7 @@ def integrate( shooting_type: Shooting = Shooting.SINGLE, integrator: SolutionIntegrator = SolutionIntegrator.OCP, to_merge: SolutionMerge | list[SolutionMerge] = None, + duplicated_times: bool = True, return_time: bool = False, ): has_direct_collocation = sum([nlp.ode_solver.is_direct_collocation for nlp in self.ocp.nlp]) > 0 @@ -730,13 +738,12 @@ def integrate( u = self._stepwise_controls.to_dict(to_merge=SolutionMerge.KEYS, scaled=False) a = self._decision_algebraic_states.to_dict(to_merge=SolutionMerge.KEYS, scaled=False) - time_vector = [] out: list = [None] * len(self.ocp.nlp) integrated_sol = None for p, nlp in enumerate(self.ocp.nlp): next_x = self._states_for_phase_integration(shooting_type, p, integrated_sol, x, u, params, a) - phase_times, integrated_sol = solve_ivp_interface( + integrated_sol = solve_ivp_interface( shooting_type=shooting_type, nlp=nlp, t=t_spans[p], @@ -747,18 +754,21 @@ def integrate( method=integrator, ) - time_vector.append(phase_times) out[p] = {} for key in nlp.states.keys(): out[p][key] = [None] * nlp.n_states_nodes for ns, sol_ns in enumerate(integrated_sol): - out[p][key][ns] = sol_ns[nlp.states[key].index, :] + if duplicated_times: + out[p][key][ns] = sol_ns[nlp.states[key].index, :] + else: + duplicated_times_condition = p == len(self.ocp.nlp) - 1 and ns == nlp.ns + out[p][key][ns] = sol_ns[nlp.states[key].index, :] if duplicated_times_condition else sol_ns[nlp.states[key].index, :-1] if to_merge: out = SolutionData.from_unscaled(self.ocp, out, "x").to_dict(to_merge=to_merge, scaled=False) - time_vector, out = self._remove_integrated_duplicates(time_vector, out, shooting_type, to_merge) if return_time: + time_vector = self.stepwise_time(to_merge=to_merge, duplicated_times=duplicated_times) return out if len(out) > 1 else out[0], time_vector if len(time_vector) > 1 else time_vector[0] else: return out if len(out) > 1 else out[0] @@ -839,57 +849,57 @@ def _states_for_phase_integration( return [(integrated_states[-1] if shooting_type == Shooting.SINGLE else decision_states[phase_idx][0]) + dx] - def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merge): - if shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) and ( - to_merge == SolutionMerge.NODES or SolutionMerge.NODES in to_merge - ): - if (isinstance(to_merge, list) and SolutionMerge.PHASES in to_merge) or to_merge == SolutionMerge.PHASES: - merged_out = {} - - redundant_index_between_phases = np.cumsum( - [len(time_vector[i]) + 1 for i in range(len(self.ocp.nlp))] - ).tolist() - time_vector = ( - [i for time_vector_sub in time_vector for i in time_vector_sub] - if len(time_vector) != 1 - else time_vector[0] - ) - - time_vector = np.array(time_vector).round(decimals=6).tolist() - unique_time_vector, unique_index = np.unique(time_vector, return_index=True) - - redundant_index = [i for i in range(len(time_vector)) if i not in unique_index] - - for key in out.keys(): - merged_out[key] = [None] * len(out[key]) - for i in range(len(out[key])): - temp_merged_out = np.delete(out[key][i], redundant_index_between_phases[:-1])[:-1][ - np.newaxis, : - ] - merged_out[key][i] = np.delete(temp_merged_out, redundant_index)[np.newaxis, :] - merged_out[key] = np.concatenate(merged_out[key], axis=0) - time_vector = unique_time_vector - - else: - merged_out = [] - for i in range(len(out)): - phase_time_vector = time_vector[i] - unique_time_vector, unique_index = np.unique(phase_time_vector, return_index=True) - redundant_index = [i for i in range(len(phase_time_vector)) if i not in unique_index] - - out_per_phase = {} - for key in out[i].keys(): - out_per_phase[key] = [None] * len(out[i][key]) - for j in range(len(out[i][key])): - temp_merged_out = np.delete(out[i][key][j], redundant_index)[:-1][np.newaxis, :] - temp_merged_out = temp_merged_out[0][1:] if i > 0 else temp_merged_out - out_per_phase[key][j] = temp_merged_out - merged_out.append(out_per_phase) - - else: - merged_out = out - - return time_vector, merged_out + # def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merge): + # if shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) and ( + # to_merge == SolutionMerge.NODES or SolutionMerge.NODES in to_merge + # ): + # if (isinstance(to_merge, list) and SolutionMerge.PHASES in to_merge) or to_merge == SolutionMerge.PHASES: + # merged_out = {} + # + # redundant_index_between_phases = np.cumsum( + # [len(time_vector[i]) + 1 for i in range(len(self.ocp.nlp))] + # ).tolist() + # time_vector = ( + # [i for time_vector_sub in time_vector for i in time_vector_sub] + # if len(time_vector) != 1 + # else time_vector[0] + # ) + # + # time_vector = np.array(time_vector).round(decimals=6).tolist() + # unique_time_vector, unique_index = np.unique(time_vector, return_index=True) + # + # redundant_index = [i for i in range(len(time_vector)) if i not in unique_index] + # + # for key in out.keys(): + # merged_out[key] = [None] * len(out[key]) + # for i in range(len(out[key])): + # temp_merged_out = np.delete(out[key][i], redundant_index_between_phases[:-1])[:-1][ + # np.newaxis, : + # ] + # merged_out[key][i] = np.delete(temp_merged_out, redundant_index)[np.newaxis, :] + # merged_out[key] = np.concatenate(merged_out[key], axis=0) + # time_vector = unique_time_vector + # + # else: + # merged_out = [] + # for i in range(len(out)): + # phase_time_vector = time_vector[i] + # unique_time_vector, unique_index = np.unique(phase_time_vector, return_index=True) + # redundant_index = [i for i in range(len(phase_time_vector)) if i not in unique_index] + # + # out_per_phase = {} + # for key in out[i].keys(): + # out_per_phase[key] = [None] * len(out[i][key]) + # for j in range(len(out[i][key])): + # temp_merged_out = np.delete(out[i][key][j], redundant_index)[:-1][np.newaxis, :] + # temp_merged_out = temp_merged_out[0][1:] if i > 0 else temp_merged_out + # out_per_phase[key][j] = temp_merged_out + # merged_out.append(out_per_phase) + # + # else: + # merged_out = out + # + # return time_vector, merged_out def _integrate_stepwise(self) -> None: """ From e0a635d6ee6d896e7be74acec6dbda4b856203bb Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Mon, 11 Mar 2024 19:08:17 -0400 Subject: [PATCH 11/18] adding assert to existing tests for dim comparision duplicates=False --- bioptim/optimization/solution/solution.py | 2 +- tests/shard4/test_simulate.py | 8 ++++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index fef66baf7..c34b8ed61 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -922,7 +922,7 @@ def _integrate_stepwise(self) -> None: unscaled: list = [None] * len(self.ocp.nlp) for p, nlp in enumerate(self.ocp.nlp): - phase_times, integrated_sol = solve_ivp_interface( + integrated_sol = solve_ivp_interface( shooting_type=Shooting.MULTIPLE, nlp=nlp, t=t_spans[p], diff --git a/tests/shard4/test_simulate.py b/tests/shard4/test_simulate.py index 34c10a351..0b810c569 100644 --- a/tests/shard4/test_simulate.py +++ b/tests/shard4/test_simulate.py @@ -606,6 +606,14 @@ def test_integrate_multiphase(shooting, integrator, ode_solver, phase_dynamics, assert sol_integrated[i][key].shape == (shapes[k], n_shooting[i] * n_steps + 1) assert states[i][key].shape == (shapes[k], n_shooting[i] * n_steps + 1) + opts = {"shooting_type": shooting, "integrator": integrator, "to_merge": [SolutionMerge.PHASES, SolutionMerge.NODES]} + sol_integrated, sol_time = sol.integrate(**opts, duplicated_times=False, return_time=True) + sol_time = np.concatenate(sol_time, axis=0) + assert len(sol_time) == len(np.unique(sol_time)) + for i in range(len(sol_integrated)): + for k, key in enumerate(states[i]): + assert len(sol_integrated[key][k]) == len(sol_time) + def test_check_models_comes_from_same_super_class(): from bioptim.examples.getting_started import example_multiphase as ocp_module From 87801318c01f946f2bc5c31683f8a391e48251f4 Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Mon, 11 Mar 2024 19:09:37 -0400 Subject: [PATCH 12/18] black --- bioptim/optimization/solution/solution.py | 11 ++++++++--- tests/shard4/test_simulate.py | 6 +++++- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index c34b8ed61..226a0c74f 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -234,8 +234,9 @@ def from_initial_guess(cls, ocp, sol: list): if len(sol[0]) != len(all_ns): raise ValueError("The time step dt array len must match the number of phases") - is_right_size = [len(s) != len(all_ns) if p != 3 and len(sol[p + 1].keys()) != 0 else False - for p, s in enumerate(sol[:1])] + is_right_size = [ + len(s) != len(all_ns) if p != 3 and len(sol[p + 1].keys()) != 0 else False for p, s in enumerate(sol[:1]) + ] if sum(is_right_size) != 0: raise ValueError("The InitialGuessList len must match the number of phases") @@ -762,7 +763,11 @@ def integrate( out[p][key][ns] = sol_ns[nlp.states[key].index, :] else: duplicated_times_condition = p == len(self.ocp.nlp) - 1 and ns == nlp.ns - out[p][key][ns] = sol_ns[nlp.states[key].index, :] if duplicated_times_condition else sol_ns[nlp.states[key].index, :-1] + out[p][key][ns] = ( + sol_ns[nlp.states[key].index, :] + if duplicated_times_condition + else sol_ns[nlp.states[key].index, :-1] + ) if to_merge: out = SolutionData.from_unscaled(self.ocp, out, "x").to_dict(to_merge=to_merge, scaled=False) diff --git a/tests/shard4/test_simulate.py b/tests/shard4/test_simulate.py index 0b810c569..0e01d4fba 100644 --- a/tests/shard4/test_simulate.py +++ b/tests/shard4/test_simulate.py @@ -606,7 +606,11 @@ def test_integrate_multiphase(shooting, integrator, ode_solver, phase_dynamics, assert sol_integrated[i][key].shape == (shapes[k], n_shooting[i] * n_steps + 1) assert states[i][key].shape == (shapes[k], n_shooting[i] * n_steps + 1) - opts = {"shooting_type": shooting, "integrator": integrator, "to_merge": [SolutionMerge.PHASES, SolutionMerge.NODES]} + opts = { + "shooting_type": shooting, + "integrator": integrator, + "to_merge": [SolutionMerge.PHASES, SolutionMerge.NODES], + } sol_integrated, sol_time = sol.integrate(**opts, duplicated_times=False, return_time=True) sol_time = np.concatenate(sol_time, axis=0) assert len(sol_time) == len(np.unique(sol_time)) From 49221a0db9bd5757e88c7ac071b829cac1e620f6 Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Tue, 12 Mar 2024 09:32:49 -0400 Subject: [PATCH 13/18] removing useless fun --- bioptim/optimization/solution/solution.py | 52 ----------------------- 1 file changed, 52 deletions(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index 226a0c74f..6f1f73ee5 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -854,58 +854,6 @@ def _states_for_phase_integration( return [(integrated_states[-1] if shooting_type == Shooting.SINGLE else decision_states[phase_idx][0]) + dx] - # def _remove_integrated_duplicates(self, time_vector, out, shooting_type, to_merge): - # if shooting_type in (Shooting.SINGLE, Shooting.SINGLE_DISCONTINUOUS_PHASE) and ( - # to_merge == SolutionMerge.NODES or SolutionMerge.NODES in to_merge - # ): - # if (isinstance(to_merge, list) and SolutionMerge.PHASES in to_merge) or to_merge == SolutionMerge.PHASES: - # merged_out = {} - # - # redundant_index_between_phases = np.cumsum( - # [len(time_vector[i]) + 1 for i in range(len(self.ocp.nlp))] - # ).tolist() - # time_vector = ( - # [i for time_vector_sub in time_vector for i in time_vector_sub] - # if len(time_vector) != 1 - # else time_vector[0] - # ) - # - # time_vector = np.array(time_vector).round(decimals=6).tolist() - # unique_time_vector, unique_index = np.unique(time_vector, return_index=True) - # - # redundant_index = [i for i in range(len(time_vector)) if i not in unique_index] - # - # for key in out.keys(): - # merged_out[key] = [None] * len(out[key]) - # for i in range(len(out[key])): - # temp_merged_out = np.delete(out[key][i], redundant_index_between_phases[:-1])[:-1][ - # np.newaxis, : - # ] - # merged_out[key][i] = np.delete(temp_merged_out, redundant_index)[np.newaxis, :] - # merged_out[key] = np.concatenate(merged_out[key], axis=0) - # time_vector = unique_time_vector - # - # else: - # merged_out = [] - # for i in range(len(out)): - # phase_time_vector = time_vector[i] - # unique_time_vector, unique_index = np.unique(phase_time_vector, return_index=True) - # redundant_index = [i for i in range(len(phase_time_vector)) if i not in unique_index] - # - # out_per_phase = {} - # for key in out[i].keys(): - # out_per_phase[key] = [None] * len(out[i][key]) - # for j in range(len(out[i][key])): - # temp_merged_out = np.delete(out[i][key][j], redundant_index)[:-1][np.newaxis, :] - # temp_merged_out = temp_merged_out[0][1:] if i > 0 else temp_merged_out - # out_per_phase[key][j] = temp_merged_out - # merged_out.append(out_per_phase) - # - # else: - # merged_out = out - # - # return time_vector, merged_out - def _integrate_stepwise(self) -> None: """ This method integrate to stepwise level the states. That is the states that are used in the dynamics and From e9cdff0940326bb27bd05ef050d78bfb3168ef77 Mon Sep 17 00:00:00 2001 From: Pariterre Date: Tue, 12 Mar 2024 11:21:49 -0400 Subject: [PATCH 14/18] Fixed new error message of casadi --- tests/shard6/test_global_stochastic_except_collocation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/shard6/test_global_stochastic_except_collocation.py b/tests/shard6/test_global_stochastic_except_collocation.py index b67f98870..4e9d751b5 100644 --- a/tests/shard6/test_global_stochastic_except_collocation.py +++ b/tests/shard6/test_global_stochastic_except_collocation.py @@ -35,7 +35,7 @@ def test_arm_reaching_muscle_driven(use_sx): match = None else: match = re.escape( - "Error in Function::call for 'tp' [MXFunction] at .../casadi/core/function.cpp:339:\n" + "Error in Function::call for 'tp' [MXFunction] at .../casadi/core/function.cpp:370:\n" ".../casadi/core/linsol_internal.cpp:65: eval_sx not defined for LinsolQr" ) with pytest.raises(RuntimeError, match=match): @@ -282,7 +282,7 @@ def test_arm_reaching_torque_driven_explicit(use_sx): match = None else: match = re.escape( - "Error in Function::call for 'tp' [MXFunction] at .../casadi/core/function.cpp:339:\n" + "Error in Function::call for 'tp' [MXFunction] at .../casadi/core/function.cpp:370:\n" ".../casadi/core/linsol_internal.cpp:65: eval_sx not defined for LinsolQr" ) with pytest.raises(RuntimeError, match=match): From ba394618b43828fe17c026162045b90750e1f25a Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Tue, 12 Mar 2024 12:00:23 -0400 Subject: [PATCH 15/18] including pr remarks --- bioptim/optimization/solution/solution.py | 24 ++++++++++++++++++++++- tests/shard4/test_simulate.py | 15 +++++++------- 2 files changed, 31 insertions(+), 8 deletions(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index 6f1f73ee5..4323e034b 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -713,6 +713,28 @@ def integrate( duplicated_times: bool = True, return_time: bool = False, ): + """ + Create a deepcopy of the Solution + + Parameters + ---------- + shooting_type: Shooting + The integration shooting type to use + integrator: SolutionIntegrator + The type of integrator to use + to_merge: SolutionMerge | list[SolutionMerge, ...] + The type of merge to perform. If None, then no merge is performed. + duplicated_times: bool + If the times should be duplicated for each node. + If False, then the returned time vector will not have any duplicated times. + return_time: bool + If the time vector should be returned + + Returns + ------- + Return the integrated states + """ + has_direct_collocation = sum([nlp.ode_solver.is_direct_collocation for nlp in self.ocp.nlp]) > 0 if has_direct_collocation and integrator == SolutionIntegrator.OCP: raise ValueError( @@ -773,7 +795,7 @@ def integrate( out = SolutionData.from_unscaled(self.ocp, out, "x").to_dict(to_merge=to_merge, scaled=False) if return_time: - time_vector = self.stepwise_time(to_merge=to_merge, duplicated_times=duplicated_times) + time_vector = np.concatenate(self.stepwise_time(to_merge=to_merge, duplicated_times=duplicated_times)) return out if len(out) > 1 else out[0], time_vector if len(time_vector) > 1 else time_vector[0] else: return out if len(out) > 1 else out[0] diff --git a/tests/shard4/test_simulate.py b/tests/shard4/test_simulate.py index 0e01d4fba..b31bd85cc 100644 --- a/tests/shard4/test_simulate.py +++ b/tests/shard4/test_simulate.py @@ -606,13 +606,14 @@ def test_integrate_multiphase(shooting, integrator, ode_solver, phase_dynamics, assert sol_integrated[i][key].shape == (shapes[k], n_shooting[i] * n_steps + 1) assert states[i][key].shape == (shapes[k], n_shooting[i] * n_steps + 1) - opts = { - "shooting_type": shooting, - "integrator": integrator, - "to_merge": [SolutionMerge.PHASES, SolutionMerge.NODES], - } - sol_integrated, sol_time = sol.integrate(**opts, duplicated_times=False, return_time=True) - sol_time = np.concatenate(sol_time, axis=0) + sol_integrated, sol_time = sol.integrate( + shooting_type=shooting, + integrator=integrator, + to_merge=[SolutionMerge.PHASES, SolutionMerge.NODES], + duplicated_times=False, + return_time=True, + ) + assert len(sol_time) == len(np.unique(sol_time)) for i in range(len(sol_integrated)): for k, key in enumerate(states[i]): From 448aafbb08e713395dfde871d798431e5768f960 Mon Sep 17 00:00:00 2001 From: Pariterre Date: Tue, 12 Mar 2024 14:26:39 -0400 Subject: [PATCH 16/18] Fixed error messages --- ...st_global_stochastic_except_collocation.py | 20 ++----------------- 1 file changed, 2 insertions(+), 18 deletions(-) diff --git a/tests/shard6/test_global_stochastic_except_collocation.py b/tests/shard6/test_global_stochastic_except_collocation.py index 4e9d751b5..167ddd239 100644 --- a/tests/shard6/test_global_stochastic_except_collocation.py +++ b/tests/shard6/test_global_stochastic_except_collocation.py @@ -30,15 +30,7 @@ def test_arm_reaching_muscle_driven(use_sx): sensory_noise_magnitude = vertcat(wPq_magnitude, wPqdot_magnitude) if use_sx: - if platform.system() == "Windows": - # It is not possible to test the error message on Windows as it uses absolute path - match = None - else: - match = re.escape( - "Error in Function::call for 'tp' [MXFunction] at .../casadi/core/function.cpp:370:\n" - ".../casadi/core/linsol_internal.cpp:65: eval_sx not defined for LinsolQr" - ) - with pytest.raises(RuntimeError, match=match): + with pytest.raises(RuntimeError, match=".*eval_sx not defined for LinsolQr"): ocp = ocp_module.prepare_socp( final_time=final_time, n_shooting=n_shooting, @@ -277,15 +269,7 @@ def test_arm_reaching_torque_driven_explicit(use_sx): bioptim_folder = os.path.dirname(ocp_module.__file__) if use_sx: - if platform.system() == "Windows": - # It is not possible to test the error message on Windows as it uses absolute path - match = None - else: - match = re.escape( - "Error in Function::call for 'tp' [MXFunction] at .../casadi/core/function.cpp:370:\n" - ".../casadi/core/linsol_internal.cpp:65: eval_sx not defined for LinsolQr" - ) - with pytest.raises(RuntimeError, match=match): + with pytest.raises(RuntimeError, match=".*eval_sx not defined for LinsolQr"): ocp = ocp_module.prepare_socp( biorbd_model_path=bioptim_folder + "/models/LeuvenArmModel.bioMod", final_time=final_time, From c79604121467b92d8bd5a73be5b99d5f084399cc Mon Sep 17 00:00:00 2001 From: Pariterre Date: Tue, 12 Mar 2024 14:29:47 -0400 Subject: [PATCH 17/18] Added comment to clarify a peculiar check --- bioptim/optimization/solution/solution.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index 4323e034b..948526316 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -486,6 +486,7 @@ def _process_time_vector( if not duplicated_times: for i in range(len(times)): for j in range(len(times[i])): + # Last node of last phase is always kept keep_condition = times[i][j].shape[0] == 1 and i == len(times) - 1 times[i][j] = times[i][j][:] if keep_condition else times[i][j][:-1] if j == len(times[i]) - 1 and i != len(times) - 1: @@ -784,6 +785,7 @@ def integrate( if duplicated_times: out[p][key][ns] = sol_ns[nlp.states[key].index, :] else: + # Last node of last phase is always kept duplicated_times_condition = p == len(self.ocp.nlp) - 1 and ns == nlp.ns out[p][key][ns] = ( sol_ns[nlp.states[key].index, :] From 4cda59c20832923d484e062e11676f6ea3527d7a Mon Sep 17 00:00:00 2001 From: Kevin CO Date: Tue, 12 Mar 2024 16:32:08 -0400 Subject: [PATCH 18/18] adapting return time to to_merge input --- bioptim/optimization/solution/solution.py | 30 ++++++++++++++++++++++- tests/shard4/test_simulate.py | 13 ++++++++++ 2 files changed, 42 insertions(+), 1 deletion(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index 948526316..188cbe281 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -797,7 +797,7 @@ def integrate( out = SolutionData.from_unscaled(self.ocp, out, "x").to_dict(to_merge=to_merge, scaled=False) if return_time: - time_vector = np.concatenate(self.stepwise_time(to_merge=to_merge, duplicated_times=duplicated_times)) + time_vector = self._return_time_vector(to_merge=to_merge, duplicated_times=duplicated_times) return out if len(out) > 1 else out[0], time_vector if len(time_vector) > 1 else time_vector[0] else: return out if len(out) > 1 else out[0] @@ -918,6 +918,34 @@ def _integrate_stepwise(self) -> None: self._stepwise_states = SolutionData.from_unscaled(self.ocp, unscaled, "x") + def _return_time_vector(self, to_merge: SolutionMerge | list[SolutionMerge], duplicated_times: bool): + """ + Returns the time vector at each node that matches stepwise_states or stepwise_controls + Parameters + ---------- + to_merge: SolutionMerge | list[SolutionMerge, ...] + The merge type to perform. If None, then no merge is performed. + duplicated_times: bool + If the times should be duplicated for each node. + If False, then the returned time vector will not have any duplicated times. + Returns + ------- + The time vector at each node that matches stepwise_states or stepwise_controls + """ + if to_merge is None: + to_merge = [] + if isinstance(to_merge, SolutionMerge): + to_merge = [to_merge] + if SolutionMerge.NODES and SolutionMerge.PHASES in to_merge: + time_vector = np.concatenate(self.stepwise_time(to_merge=to_merge, duplicated_times=duplicated_times)) + elif SolutionMerge.NODES in to_merge: + time_vector = self.stepwise_time(to_merge=to_merge, duplicated_times=duplicated_times) + for i in range(len(self.ocp.nlp)): + time_vector[i] = np.concatenate(time_vector[i]) + else: + time_vector = self.stepwise_time(to_merge=to_merge, duplicated_times=duplicated_times) + return time_vector + def interpolate(self, n_frames: int | list | tuple, scaled: bool = False): """ Interpolate the states diff --git a/tests/shard4/test_simulate.py b/tests/shard4/test_simulate.py index b31bd85cc..9d093a9da 100644 --- a/tests/shard4/test_simulate.py +++ b/tests/shard4/test_simulate.py @@ -619,6 +619,19 @@ def test_integrate_multiphase(shooting, integrator, ode_solver, phase_dynamics, for k, key in enumerate(states[i]): assert len(sol_integrated[key][k]) == len(sol_time) + sol_integrated, sol_time = sol.integrate( + shooting_type=shooting, + integrator=integrator, + to_merge=[SolutionMerge.NODES], + duplicated_times=False, + return_time=True, + ) + for i in range(len(sol_time)): + assert len(sol_time[i]) == len(np.unique(sol_time[i])) + for i in range(len(sol_integrated)): + for k, key in enumerate(states[i]): + assert len(sol_integrated[i][key][k]) == len(sol_time[i]) + def test_check_models_comes_from_same_super_class(): from bioptim.examples.getting_started import example_multiphase as ocp_module