Inplace (#297)

* args and setstateequation for trapeze

* path for trapeze

* todo: compute whole dynamics

* undefined reference

* try extended dynaùics first

* try vector of structs

* tests ok

* with dynamics matrix for trapeze

* with lagrange vector for trapeze

* comment

* prepare constraints main loop for multithread

* try c_block views

* tests ok for c_block (except mintf with lagrange cost equal to 1)

* cleanup

* compact setPathConstraints

* loop test

* multi-thread always fails on main constraints loop

* extended dynamics with vectorize; use with new separate (vector) getters for x, u

* midpoint prof

* simple args for midpoint

* shift index for get variables

* grid in DOCP

* trapeze, work array for dynamics

* tests ok

* sync from main

* tests ok

* prepare for inplace

* better prof

* btime

* total allocs do not match .mem values

* more tests

* tests ok

* standard tests ok with inplace

* inplace for midpoint

* check allocs for inplace case

* mem files for inplace/outplace comparison

* grmbl

* prepare for vectorized ocp functions and extended dynamics

* prepare dynamics extended

* todo: inplace version, then optimize

* fbll

* todo: check tests

* tests ok

* still ok

* removed some allocs

* work arrays in disc struct are worse, even with inplace getters

* code warntype

* prof

* test dyn

* inlined path constraints

* getters for v, t0, tf and time grid now defined as part of DOCP

* getter for x defined in discretization struct

* WORSE with getters in disc and objective in docp...

* OK with multiple definition of objective in docp :-)

* cleanup

* clean then sync and merge

* cleaned

* cleanup

* cleanup

* sync 2

* sync ok

* compat for docs

* jet

* ignore

* check getters

* simplification: dynamics_ext is now always inplace even if OCP is not

* views do seem slightly better than slices for getters

* check t_i type

* redefine type_stable ti/tf in docp

* miminalist objective

* getters without allocs; todo: improve _x by passing full argument

* allocs tests for local mayer

* reintroduced scal/vec getters that seem to give less allocs than scalarization functions

* grr

* Val in external getter does not fix type

* parametric types in DOCP for scal/vect variables seem promising

* in progress: constraints with param getters

* test param

* tests param ok except the 3 midpoint

* check block

* fixing regression on t_i types; still some errors in tests

* all tests (except midpoint) ok

* allocs and runtime dispatch for all ocp functions, which show as type Any

* mayer test

* prof

.gitignore

@@ -1,3 +1,6 @@
+# Files generated by profiling
+*.pb.gz
+
 # Files generated by invoking Julia with --code-coverage
 *.jl.cov
 *.jl.*.cov

Project.toml

@@ -1,7 +1,7 @@
 name = "CTDirect"
 uuid = "790bbbee-bee9-49ee-8912-a9de031322d5"
 authors = ["Olivier Cots <[email protected]>"]
-version = "0.12.1"
+version = "0.13.0"
 
 [deps]
 ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
@@ -24,7 +24,7 @@ CTSolveExtMadNLP = ["MadNLP"]
 
 [compat]
 ADNLPModels = "0.8"
-CTBase = "0.13"
+CTBase = "0.14"
 DocStringExtensions = "0.9"
 HSL = "0.4"
 JLD2 = "0.4" # 0.5 has ugly warnings for functions

benchmark/benchmark.jl

@@ -1,10 +1,15 @@
 # Benchmark
-include("../test/deps.jl")
+using CTDirect
+using CTBase
+
+using LinearAlgebra
+using NLPModelsIpopt
+
 using Printf
 
 using MKL # Replace OpenBLAS with Intel MKL +++ should be an option
 
-using MadNLPMumps
+#using MadNLPMumps
 
 #######################################################
 # load examples library
@@ -21,7 +26,7 @@ function bench(;
     precompile = true,
     display = false,
     verbose = true,
-    discretization = :trapeze,
+    disc_method = :trapeze
 )
 
     #######################################################
@@ -67,7 +72,7 @@ function bench(;
                 linear_solver = linear_solver,
                 max_iter = 0,
                 display = display,
-                discretization = discretization,
+                disc_method = disc_method
             )
             t_precomp += t
         end
@@ -86,7 +91,7 @@ function bench(;
             linear_solver = linear_solver,
             grid_size = grid_size,
             tol = tol,
-            discretization = discretization,
+            disc_method = disc_method
         )
         if !isnothing(problem[:obj]) && !isapprox(sol.objective, problem[:obj], rtol = 5e-2)
             error(

benchmark/prof.jl

@@ -1,106 +1,184 @@
-# +++ TODO: make function with bools as args ?
-# Profiling
+# tests to check allocations in particular
 using CTDirect
 using CTBase
 
 using LinearAlgebra
 using NLPModelsIpopt
 using BenchmarkTools
-#using Traceur
-#using Profile
-#using PProf
+using Profile
+using PProf
 using JET
 
-precompile = true
-test_objective = true
-test_constraints = true
-test_transcription = true
-test_solve = true
+include("../test/problems/goddard.jl")
+#include("../test/problems/double_integrator.jl")
 
-test_code_warntype = false
-test_jet = false
+# local version of mayer cost
+function local_mayer(obj, x0, xf, v)
+    obj[1] = xf[3]
+    return
+end
 
-# define OCP
-include("../test/problems/goddard.jl")
-prob = goddard_all()
-ocp = prob[:ocp]
-grid_size = 100
-println("Load problem ", prob[:name])
-
-if precompile
-    println("Precompilation")
-    if test_transcription
-        docp, nlp = direct_transcription(ocp, grid_size = grid_size)
+function init(;in_place, grid_size, disc_method)
+    if in_place
+        prob = goddard_all_inplace()
+        #prob = goddard_a()
+        #prob = double_integrator_a()
+    else
+        prob = goddard_all()
+        #prob = goddard()       
+        #prob = double_integrator_mintf()
     end
-    if test_solve
-        direct_solve(ocp, grid_size = grid_size, display = false, max_iter = 2)
+    ocp = prob[:ocp]
+    docp = CTDirect.DOCP(ocp, grid_size=grid_size, time_grid=CTDirect.__time_grid(), disc_method=string(disc_method))
+    xu = CTDirect.DOCP_initial_guess(docp)
+    return prob, docp, xu
+end
+
+
+function test_unit(;test_get=false, test_dyn=false, test_unit_cons=false, test_mayer=false, test_obj=true, test_block=false, test_cons=true, test_trans=true, test_solve=true, warntype=false, jet=false, profile=false, grid_size=100, disc_method=:trapeze, in_place=true)
+
+    # define problem and variables
+    prob, docp, xu = init(in_place=in_place, grid_size=grid_size, disc_method=disc_method)
+    disc = docp.discretization
+    #= OK, same as calling the functions with docp
+    NLP_objective = (xu) -> CTDirect.DOCP_objective(xu, docp)
+    NLP_constraints! = (c, xu) -> CTDirect.DOCP_constraints!(c, xu, docp) =#
+    c = fill(666.666, docp.dim_NLP_constraints)
+    work = similar(xu, docp.dim_NLP_x)
+
+    # getters
+    if test_get
+        println("Getters")
+        print("t"); @btime CTDirect.get_final_time($xu, $docp)
+        print("x"); @btime CTDirect.get_OCP_state_at_time_step($xu, $docp, 1)
+        print("u"); @btime CTDirect.get_OCP_control_at_time_step($xu, $docp, 1)
+        print("v"); @btime CTDirect.get_OCP_variable($xu, $docp)
+        if warntype
+            @code_warntype CTDirect.get_final_time(xu, docp)
+            @code_warntype CTDirect.get_time_grid(xu, docp)
+            @code_warntype CTDirect.get_OCP_state_at_time_step(xu, docp, 1)
+            @code_warntype CTDirect.get_OCP_control_at_time_step(xu, docp, 1)
+            @code_warntype CTDirect.get_OCP_variable(xu, docp)
+        end
     end
-    if test_objective
-        CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp)
+
+    f = similar(xu, docp.dim_NLP_x)
+
+    if test_dyn
+        if in_place
+            print("dynamics_ext"); @btime $docp.dynamics_ext($f, $t, $x, $u, $v)
+            warntype && @code_warntype docp.dynamics_ext(f, t, x, u, v)
+            Profile.clear_malloc_data()
+            docp.dynamics_ext(f, t, x, u, v)
+        else
+            print("dynamics_ext"); @btime $docp.dynamics_ext($t, $x, $u, $v)
+            warntype && @code_warntype docp.dynamics_ext(t, x, u, v)
+            Profile.clear_malloc_data()
+            docp.dynamics_ext(t, x, u, v)
+        end
     end
-    if test_constraints
-        CTDirect.DOCP_constraints!(
-            zeros(docp.dim_NLP_constraints),
-            CTDirect.DOCP_initial_guess(docp),
-            docp,
-        )
+
+    if test_unit_cons
+        if in_place
+            print("u cons"); @btime $docp.control_constraints[2]($c, $t, $u, $v)
+            print("x cons"); @btime $docp.state_constraints[2]($c, $t, $x, $v)
+            print("xu cons"); @btime $docp.mixed_constraints[2]($c, $t, $x, $u, $v)
+            if warntype
+                @code_warntype docp.control_constraints[2](c, t, u, v)
+                @code_warntype docp.state_constraints[2](c, t, x, v)
+                @code_warntype docp.mixed_constraints[2](c, t, x, u, v)
+            end
+        else
+            print("u cons"); @btime $docp.control_constraints[2]($t, $u, $v)
+            print("x cons"); @btime $docp.state_constraints[2]($t, $x, $v)
+            print("xu cons"); @btime $docp.mixed_constraints[2]($t, $x, $u, $v)
+            if warntype
+                @code_warntype docp.control_constraints[2](t, u, v)
+                @code_warntype docp.state_constraints[2](t, x, v)
+                @code_warntype docp.mixed_constraints[2](t, x, u, v)
+            end
+        end
     end
-end
 
-# evaluation
-if test_objective
-    println("Timed objective")
-    @btime CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp)
-end
-if test_constraints
-    println("Timed constraints")
-    @btime CTDirect.DOCP_constraints!(
-        zeros(docp.dim_NLP_constraints),
-        CTDirect.DOCP_initial_guess(docp),
-        docp,
-    )
-end
+    # objective
+    if test_mayer
+        n = docp.dim_OCP_x
+        nx = docp.dim_NLP_x
+        m = docp.dim_NLP_u
+        N = docp.dim_NLP_steps
+        x0 = CTDirect.get_OCP_state_at_time_step(xu, docp, 1)
+        xf = CTDirect.get_OCP_state_at_time_step(xu, docp, N+1)
+        v = CTDirect.get_OCP_variable(xu, docp)
+        obj = similar(xu,1)
+
+        # local mayer
+        println("")
+        print("Local Mayer: views for x0/xf and scalar v"); @btime local_mayer($obj, (@view $xu[1:$n]), (@view $xu[($nx + $m) * $N + 1: ($nx + $m) * $N + $n]), $xu[end]) # OK
+        print("Local Mayer: param scal/vec getters"); @btime local_mayer($obj, $x0, $xf, $v) # OK
+        print("OCP Mayer: param scal/vec getters"); @btime $docp.ocp.mayer($obj, $x0, $xf, $v) # 3 allocs (112)
 
-# transcription
-if test_transcription
-    println("Timed transcription")
-    @btime docp, nlp = direct_transcription(ocp, grid_size = grid_size)
-end
+        warntype && @code_warntype docp.ocp.mayer(obj, x0, xf, v)
+        jet && display(@report_opt docp.ocp.mayer(obj, x0, xf, v))
+        if profile
+            Profile.Allocs.@profile sample_rate=1.0 docp.ocp.mayer(obj, x0, xf, v)
+            results = Profile.Allocs.fetch()
+            PProf.Allocs.pprof()
+        end
+    end
 
-# full solve
-if test_solve
-    println("Timed full solve")
-    @btime sol = direct_solve(ocp, grid_size = grid_size, display = false)
-end
+    if test_obj
+        print("Objective"); @btime CTDirect.DOCP_objective($xu, $docp)
+        warntype && @code_warntype CTDirect.DOCP_objective(xu, docp) # quasi OK (inplace/outplace for ocp.mayer return ?)
+        jet && display(@report_opt CTDirect.DOCP_objective(xu, docp))
+        if profile
+            Profile.Allocs.@profile sample_rate=1.0 CTDirect.DOCP_objective(xu, docp)
+            results = Profile.Allocs.fetch()
+            PProf.Allocs.pprof()
+        end
+    end
 
-if test_code_warntype
-    if test_objective
-        # NB. Pb with the mayer part: obj is type unstable (Any) because ocp.mayer is Union(Mayer,nothing), even for mayer problems (also, we should not even enter this code part for lagrange problems since has_mayer us defined as const in DOCP oO ...).
-        @code_warntype CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp)
+    if test_block
+        times = CTDirect.get_time_grid(xu, docp) # type OK
+        i = 1
+        v = CTDirect.get_OCP_variable(xu, docp)
+        work = CTDirect.setWorkArray(docp, xu, times, v)
+        print("Constraints block")
+        @btime CTDirect.setConstraintBlock!($docp, $c, $xu, $v, $times, $i, $work)
+        warntype && @code_warntype CTDirect.setConstraintBlock!(docp, c, xu, v, times, i, work)
+        jet && display(@report_opt CTDirect.setConstraintBlock!(docp, c, xu, v, times, i, work))
+        if profile
+            Profile.Allocs.@profile sample_rate=1.0 CTDirect.setConstraintBlock!(docp, c, xu, v, times, i, work)
+            results = Profile.Allocs.fetch()
+            PProf.Allocs.pprof()
+        end        
     end
-    if test_constraints
-        # OK !
-        @code_warntype CTDirect.DOCP_constraints!(
-            zeros(docp.dim_NLP_constraints),
-            CTDirect.DOCP_initial_guess(docp),
-            docp,
-        )
+
+    # DOCP_constraints
+    if test_cons
+        print("Constraints"); @btime CTDirect.DOCP_constraints!($c, $xu, $docp)
+        any(c.==666.666) && error("undefined values in constraints ",c)
+        warntype && @code_warntype CTDirect.DOCP_constraints!(c, xu, docp)
+        jet && display(@report_opt CTDirect.DOCP_constraints!(c, xu, docp))
+        if profile
+            Profile.Allocs.@profile sample_rate=1.0 CTDirect.DOCP_constraints!(c, xu, docp)
+            results = Profile.Allocs.fetch()
+            PProf.Allocs.pprof()
+        end
     end
-end
 
-if test_jet
-    if test_objective
-        # 4 possible errors
-        # due to the ocp.mayer type problem cf above
-        @report_opt CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp)
+    # transcription
+    if test_trans
+        print("Transcription"); @btime direct_transcription($prob.ocp, grid_size=$grid_size)
     end
-    if test_constraints
-        # 50 possible errors: some getindex (Integer vs Int...)
-        # all variables x,u,v
-        @report_opt CTDirect.DOCP_constraints!(
-            zeros(docp.dim_NLP_constraints),
-            CTDirect.DOCP_initial_guess(docp),
-            docp,
-        )
+
+    # solve
+    if test_solve
+        sol = direct_solve(prob.ocp, display=false, grid_size=grid_size)
+        if !isapprox(sol.objective, prob.obj, rtol=1e-2)
+            error("objective mismatch: ", sol.objective, " vs ", prob.obj)
+        end
+        print("Solve"); @btime direct_solve($prob.ocp, display=false, grid_size=$grid_size)
     end
+
 end
+

docs/Project.toml

@@ -10,8 +10,8 @@ NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
-CTBase = "0.13"
-CTDirect = "0.12"
+CTBase = "0.14"
+CTDirect = "0.13"
 Documenter = "1"
 DocumenterMermaid = "0.1"
 HSL = "0.4"

src/default.jl

@@ -6,7 +6,7 @@ $(TYPEDSIGNATURES)
 Used to set the default discretization method.
 The default value is `trapeze`.
 """
-__discretization() = "trapeze"
+__disc_method() = "trapeze"
 
 """
 $(TYPEDSIGNATURES)