some cleanup (#106)

* basic load/save solution in julia JLD2 format

* added functions for load/save; move to CTBase once finished

* add new struct for OCP solution in discrete form for export

* load / save in JSON format (interpolated solution)

* adjusted split in constructor for OCP solution

* more work needed to make the raw constructor more generic :(

src/solution.jl

@@ -6,63 +6,70 @@ Build OCP functional solution from DOCP vector solution (given as a GenericExecu
 function OCPSolutionFromDOCP(docp, docp_solution_ipopt)
 
     # could pass some status info too (get_status ?)
-    return OCPSolutionFromDOCP_raw(docp, docp_solution_ipopt.solution, objective=docp_solution_ipopt.objective, constraints_violation=docp_solution_ipopt.primal_feas, iterations=docp_solution_ipopt.iter,multipliers_constraints=docp_solution_ipopt.multipliers, multipliers_LB=docp_solution_ipopt.multipliers_L, multipliers_UB=docp_solution_ipopt.multipliers_U, message=docp_solution_ipopt.solver_specific[:internal_msg])
-end
-
-
-"""
-$(TYPEDSIGNATURES)
+    solution = docp_solution_ipopt.solution
 
-Build OCP functional solution from DOCP vector solution (given as raw variables and multipliers plus some optional infos)
-"""
-function OCPSolutionFromDOCP_raw(docp, solution; objective=nothing, constraints_violation=nothing, iterations=0, multipliers_constraints=nothing, multipliers_LB=nothing, multipliers_UB=nothing, message=nothing)
-
-    # NB. still missing: stopping and success info...
-
-    # set objective if needed
-    if objective==nothing
-        objective = DOCP_objective(solution, docp)
-        println("Recomputed raw objective ", objective)
-    end
+    # time grid
+    N = docp.dim_NLP_steps
+    t0 = get_initial_time(solution, docp)
+    tf = max(get_final_time(solution, docp), t0 + 1e-9)
+    T = collect(LinRange(t0, tf, N+1))
+
     # adjust objective sign for maximization problems
-    if !is_min(docp.ocp)
-        objective = - objective
+    if is_min(docp.ocp)
+        objective = docp_solution_ipopt.objective
+    else        
+        objective = - docp_solution_ipopt.objective
     end
 
     # recompute value of constraints at solution
     # NB. the constraint formulation is LB <= C <= UB
     constraints = zeros(docp.dim_NLP_constraints)
     DOCP_constraints!(constraints, solution, docp)
-    # set constraint violation if needed
-    # +++ is not saved in OCP solution currently...
-    if constraints_violation==nothing
-        constraints_check = zeros(docp.dim_NLP_constraints)
-        DOCP_constraints_check!(constraints_check, constraints, docp)
-        println("Recomputed constraints violation ", norm(constraints_check, Inf))
-        variables_check = zeros(docp.dim_NLP_variables)
-        DOCP_variables_check!(variables_check, solution, docp)
-        println("Recomputed variable bounds violation ", norm(variables_check, Inf))
-        constraints_violation = norm(append!(variables_check, constraints_check), Inf)
 
-    end
-
     # parse NLP variables, constraints and multipliers 
-    X, U, v, P, sol_control_constraints, sol_state_constraints, sol_mixed_constraints, sol_variable_constraints, mult_control_constraints, mult_state_constraints, mult_mixed_constraints, mult_variable_constraints, mult_state_box_lower, mult_state_box_upper, mult_control_box_lower, mult_control_box_upper, mult_variable_box_lower, mult_variable_box_upper = parse_DOCP_solution(docp, solution, multipliers_constraints, multipliers_LB, multipliers_UB, constraints)
+    X, U, v, P, sol_control_constraints, sol_state_constraints, sol_mixed_constraints, sol_variable_constraints, mult_control_constraints, mult_state_constraints, mult_mixed_constraints, mult_variable_constraints, mult_state_box_lower, mult_state_box_upper, mult_control_box_lower, mult_control_box_upper, mult_variable_box_lower, mult_variable_box_upper = parse_DOCP_solution(docp, solution, docp_solution_ipopt.multipliers, docp_solution_ipopt.multipliers_L, docp_solution_ipopt.multipliers_U, constraints)
 
-    # time grid
-    N = docp.dim_NLP_steps
-    t0 = get_initial_time(solution, docp)
-    tf = max(get_final_time(solution, docp), t0 + 1e-9)
-    T = collect(LinRange(t0, tf, N+1))
+    # build and return OCP solution
+    return OCPSolutionFromDOCP_raw(docp, T, X, U, v, P,
+    objective=objective, iterations=docp_solution_ipopt.iter,constraints_violation=docp_solution_ipopt.primal_feas, 
+    message=String(docp_solution_ipopt.solver_specific[:internal_msg]),
+    sol_control_constraints=sol_control_constraints, sol_state_constraints=sol_state_constraints, sol_mixed_constraints=sol_mixed_constraints, sol_variable_constraints=sol_variable_constraints, mult_control_constraints=mult_control_constraints, mult_state_constraints=mult_state_constraints, mult_mixed_constraints=mult_mixed_constraints, mult_variable_constraints=mult_variable_constraints, mult_state_box_lower=mult_state_box_lower, mult_state_box_upper=mult_state_box_upper, mult_control_box_lower=mult_control_box_lower, mult_control_box_upper=mult_control_box_upper, mult_variable_box_lower=mult_variable_box_lower, mult_variable_box_upper=mult_variable_box_upper)
+
+end
+
+
+"""
+$(TYPEDSIGNATURES)
+    
+Build OCP functional solution from DOCP vector solution (given as raw variables and multipliers plus some optional infos)
+"""
+# +++ use tuples for more compact arguments
+
+# +++ try to reuse this for the discrete json solution !
+
+# +++ need to remove docp from this one !
+# +++ this means dimensions, 
+# +++ ocp for copy! (-_-) ie dimensions also ?
+# +++ boolean indicators for various constraints types
+# (could check for nothing values instead ?)
+# add a tuple for dimensions
+# a tupple for indicators
+# and do the copy manually ?
+
+function OCPSolutionFromDOCP_raw(docp, T, X, U, v, P;
+    objective=0, iterations=0, constraints_violation=0,
+    message="No msg", stopping=nothing, success=nothing,
+    sol_control_constraints=nothing, sol_state_constraints=nothing, sol_mixed_constraints=nothing, sol_variable_constraints=nothing, mult_control_constraints=nothing, mult_state_constraints=nothing, mult_mixed_constraints=nothing, mult_variable_constraints=nothing, mult_state_box_lower=nothing, 
+    mult_state_box_upper=nothing, mult_control_box_lower=nothing, mult_control_box_upper=nothing, mult_variable_box_lower=nothing, mult_variable_box_upper=nothing)
 
     # variables: remove additional state for lagrange cost
     x = ctinterpolate(T, matrix2vec(X[:,1:docp.ocp.state_dimension], 1))
     p = ctinterpolate(T[1:end-1], matrix2vec(P[:,1:docp.ocp.state_dimension], 1))
     u = ctinterpolate(T, matrix2vec(U, 1))
 
     # generate ocp solution
-    sol = OptimalControlSolution() # +++ constructor with ocp as argument ?
-    copy!(sol, docp.ocp)
+    sol = OptimalControlSolution()
+    copy!(sol, docp.ocp) # +++ use constructor with ocp as argument instead of this ?
     sol.times = T
     # use scalar output for x,u,v,p if dim=1
     sol.state = (sol.state_dimension==1) ? deepcopy(t -> x(t)[1]) : deepcopy(t -> x(t)) 
@@ -71,9 +78,10 @@ function OCPSolutionFromDOCP_raw(docp, solution; objective=nothing, constraints_
     sol.variable = (sol.variable_dimension==1) ? v[1] : v
     sol.objective = objective
     sol.iterations = iterations
-    sol.stopping = nothing #+++
-    sol.message = isnothing(message) ? "No msg" : String(message)
-    sol.success  = nothing #+++
+    sol.stopping = stopping
+    sol.message = message
+    sol.success  = success
+    sol.infos[:constraints_violation] = constraints_violation
 
     # nonlinear constraints and multipliers
     if docp.has_state_constraints
@@ -123,7 +131,6 @@ function OCPSolutionFromDOCP_raw(docp, solution; objective=nothing, constraints_
     end
 
     return sol
-
 end
 
 
@@ -136,7 +143,6 @@ function parse_DOCP_solution(docp, solution, multipliers_constraints, multiplier
 
     # states and controls variables, with box multipliers
     N = docp.dim_NLP_steps
-
     X = zeros(N+1,docp.dim_NLP_state)
     U = zeros(N+1,docp.ocp.control_dimension)
     v = get_variable(solution, docp)
@@ -244,3 +250,36 @@ function parse_DOCP_solution(docp, solution, multipliers_constraints, multiplier
 
     return X, U, v, P, sol_control_constraints, sol_state_constraints, sol_mixed_constraints, sol_variable_constraints, mult_control_constraints, mult_state_constraints, mult_mixed_constraints, mult_variable_constraints, mult_state_box_lower, mult_state_box_upper, mult_control_box_lower, mult_control_box_upper, mult_variable_box_lower, mult_variable_box_upper
 end
+
+
+    #return OCPSolutionFromDOCP_raw(docp, docp_solution_ipopt.solution, objective=docp_solution_ipopt.objective, constraints_violation=docp_solution_ipopt.primal_feas, iterations=docp_solution_ipopt.iter,multipliers_constraints=docp_solution_ipopt.multipliers, multipliers_LB=docp_solution_ipopt.multipliers_L, multipliers_UB=docp_solution_ipopt.multipliers_U, message=docp_solution_ipopt.solver_specific[:internal_msg])
+
+#= OLD
+    # NB. still missing: stopping and success info...
+    # set objective if needed
+    if objective==nothing
+        objective = DOCP_objective(solution, docp)
+        println("Recomputed raw objective ", objective)
+    end
+    # adjust objective sign for maximization problems
+    if !is_min(docp.ocp)
+        objective = - objective
+    end
+
+    # recompute value of constraints at solution
+    # NB. the constraint formulation is LB <= C <= UB
+    constraints = zeros(docp.dim_NLP_constraints)
+    DOCP_constraints!(constraints, solution, docp)
+    # set constraint violation if needed
+    # +++ is not saved in OCP solution currently...
+    if constraints_violation==nothing
+        constraints_check = zeros(docp.dim_NLP_constraints)
+        DOCP_constraints_check!(constraints_check, constraints, docp)
+        println("Recomputed constraints violation ", norm(constraints_check, Inf))
+        variables_check = zeros(docp.dim_NLP_variables)
+        DOCP_variables_check!(variables_check, solution, docp)
+        println("Recomputed variable bounds violation ", norm(variables_check, Inf))
+        constraints_violation = norm(append!(variables_check, constraints_check), Inf)
+
+    end
+=#

src/utils.jl

@@ -212,8 +212,6 @@ end
 # - constructor to recreate OptimalControlSolution from this one
 mutable struct OCP_Solution_discrete
 
-    grid_size
-    objective
     times
     #initial_time_name::Union{String, Nothing}=nothing
     #final_time_name::Union{String, Nothing}=nothing
@@ -231,7 +229,7 @@ mutable struct OCP_Solution_discrete
     #variable_name::Union{String, Nothing}=nothing
     variable
     costate
-    #objective::Union{Nothing, ctNumber}=nothing
+    objective
     #iterations::Union{Nothing, Integer}=nothing
     #stopping::Union{Nothing, Symbol}=nothing # the stopping criterion
     #message::Union{Nothing, String}=nothing # the message corresponding to the stopping criterion
@@ -242,7 +240,7 @@ mutable struct OCP_Solution_discrete
     function OCP_Solution_discrete(solution::OptimalControlSolution)
         solution_d = new()
 
-        # raw copy
+        # raw copy +++ reuse the copy! in CTBase ?
         solution_d.objective = solution.objective
         solution_d.times = solution.times
         solution_d.state_dimension = solution.state_dimension
@@ -251,16 +249,9 @@ mutable struct OCP_Solution_discrete
         solution_d.variable = solution.variable
 
         # interpolate functions into vectors
-        # +++ ther *must* be a quicker way to do this -_-
-        solution_d.grid_size = length(solution_d.times) - 1
-        solution_d.state = zeros(solution_d.grid_size+1, solution_d.state_dimension)
-        solution_d.control = zeros(solution_d.grid_size+1, solution_d.control_dimension)
-        solution_d.costate = zeros(solution_d.grid_size+1, solution_d.state_dimension)
-        for i in 1:solution_d.grid_size
-            solution_d.state[i,:] .= solution.state(solution_d.times[i])
-            solution_d.control[i,:] .= solution.control(solution_d.times[i])
-            solution_d.costate[i,:] .= solution.costate(solution_d.times[i])
-        end
+        solution_d.state = solution.state.(solution_d.times)
+        solution_d.control = solution.control.(solution_d.times)
+        solution_d.costate = solution.costate.(solution_d.times)
         return solution_d
     end
 end
@@ -294,8 +285,19 @@ function load_OCP_solution(filename_prefix="solution"; format="JLD2")
     elseif format == "JSON"
         json_string = read(filename_prefix * ".json", String)
         return JSON3.read(json_string)
+        #+++ parse JSON object containing interpolated solution
+        #+++ then call raw constructor for OCP solution
     else
         println("ERROR: save_OCP_solution: format should be JLD2 or JSON, received ", format)
         return nothing
     end
 end
+
+"""
+$(TYPEDSIGNATURES)
+ 
+Parse interpolated OCP solution saved in JSON format
+"""
+function parse_JSON_solution(json_solution)
+end
+

test/suite/misc.jl

@@ -39,8 +39,8 @@ docp = directTranscription(ocp, grid_size=100)
 nlp = getNLP(docp)
 
 # test OCPSolutionFromDOCP_raw
-dsol = solve(docp, print_level=0, tol=1e-12)
-sol_raw = OCPSolutionFromDOCP_raw(docp, dsol.solution)
+#dsol = solve(docp, print_level=0, tol=1e-12)
+#sol_raw = OCPSolutionFromDOCP_raw(docp, dsol.solution)
 
 # test save / load solution in JLD2 format
 @testset verbose = true showtiming = true ":save_load :JLD2" begin

test/test_misc.jl

@@ -38,8 +38,8 @@ docp2 = directTranscription(ocp, grid_size=100, init=sol)
 dsol2 = solve(docp2, print_level=5, tol=1e-12)
 
 # test OCPSolutionFromDOCP_raw
-println("\nRebuild OCP solution from raw vector")
-sol3 = OCPSolutionFromDOCP_raw(docp2, dsol2.solution)
+#println("\nRebuild OCP solution from raw vector")
+#sol3 = OCPSolutionFromDOCP_raw(docp2, dsol2.solution)
 
 # save / load solution in JLD2 format
 save_OCP_solution(sol, filename_prefix="solution_test")