diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 86a0bc9..509e71e 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -33,15 +33,15 @@ jobs: # - x86 steps: - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@v2 + - uses: julia-actions/setup-julia@v3 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v5 + - uses: codecov/codecov-action@v6 with: files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index df4a0a7..8250dc2 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -25,8 +25,8 @@ jobs: DOC_TEMPLATE_VERSION: "v0.7.0" # Change this to the specific tag version you want steps: - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@v2 - - uses: julia-actions/cache@v2 + - uses: julia-actions/setup-julia@v3 + - uses: julia-actions/cache@v3 - name: Use Documentation Template run: | ./docs/get_docs_utils.sh ${{ env.DOC_TEMPLATE_VERSION }} diff --git a/.github/workflows/nightly.yml b/.github/workflows/nightly.yml index bb70dea..2626b00 100644 --- a/.github/workflows/nightly.yml +++ b/.github/workflows/nightly.yml @@ -25,10 +25,10 @@ jobs: - x64 steps: - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@v2 + - uses: julia-actions/setup-julia@v3 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 diff --git a/.gitignore b/.gitignore index 2117a39..60d314a 100644 --- a/.gitignore +++ b/.gitignore @@ -33,10 +33,13 @@ Manifest.toml docs/Manifest.toml # Project specific ignores below -# generated example artifacts +# generated example artifacts /examples/**/plots/ /examples/**/trajectories/ +# benchmark output artifacts +/benchmark/results/ + # external pkgs and configs pardiso.lic /.CondaPkg/ diff --git a/Project.toml b/Project.toml index 796bc0e..fdc9bd2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "DirectTrajOpt" uuid = "c823fa1f-8872-4af5-b810-2b9b72bbbf56" -version = "0.9.0" +version = "0.9.2" authors = ["Aaron Trowbridge and contributors"] [deps] @@ -40,10 +40,10 @@ LinearAlgebra = "1.10, 1.11, 1.12" MadNLP = "0.9" MathOptInterface = "1.49" NamedTrajectories = "0.8" -OrdinaryDiffEqTsit5 = "1.9" +OrdinaryDiffEqTsit5 = "1.9, 2" Random = "1.10, 1.11, 1.12" Reexport = "1.2" -SciMLBase = "2.145.0 - 2.145, 2" +SciMLBase = "2.148, 3" SparseArrays = "1.10, 1.11, 1.12" Test = "1.10, 1.11, 1.12" TestItemRunner = "1.1" diff --git a/benchmark/Project.toml b/benchmark/Project.toml new file mode 100644 index 0000000..f08b1cf --- /dev/null +++ b/benchmark/Project.toml @@ -0,0 +1,17 @@ +[deps] +DirectTrajOpt = "c823fa1f-8872-4af5-b810-2b9b72bbbf56" +ExponentialAction = "e24c0720-ea99-47e8-929e-571b494574d3" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +HarmoniqsBenchmarks = "f45d0b76-2d23-4568-9599-481e0da131db" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +NamedTrajectories = "538bc3a1-5ab9-4fc3-b776-35ca1e893e08" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[sources] +DirectTrajOpt = {path = ".."} +# TODO: drop rev pin once HarmoniqsBenchmarks.jl#1 (feat/alloc-profile) merges +HarmoniqsBenchmarks = {url = "https://github.com/harmoniqs/HarmoniqsBenchmarks.jl", rev = "feat/alloc-profile"} diff --git a/benchmark/alloc_profile.jl b/benchmark/alloc_profile.jl new file mode 100644 index 0000000..ae5bf8d --- /dev/null +++ b/benchmark/alloc_profile.jl @@ -0,0 +1,133 @@ +# ============================================================================= +# Ipopt + MadNLP allocation profile — bilinear toy problem +# +# Runs `solve!` once per solver under Profile.Allocs via benchmark_memory! +# from HarmoniqsBenchmarks.jl and saves the sampled trace to +# benchmark/results/allocs/ for hot-path triage. The Piccolissimo alloc- +# profile testitem covers the Altissimo side; this script is the sibling +# for the in-tree NLP solvers. +# +# Uses the same `bilinear_dynamics_and_trajectory` fixture the main test +# suite uses, so the profiled problem is deterministic and small (N=10, +# 4-state × 2-control) — we care about allocation *patterns*, not absolute +# counts on a production-size problem. +# +# Run: +# julia --project=benchmark benchmark/alloc_profile.jl +# ============================================================================= + +using Random +using NamedTrajectories +using SparseArrays +using LinearAlgebra +using DirectTrajOpt +using MathOptInterface +const MOI = MathOptInterface +using Ipopt +using MadNLP +using HarmoniqsBenchmarks + +# Resolve the MadNLPSolverExt extension module so MadNLPOptions is accessible +# (matches the pattern used in Piccolissimo.jl/benchmark/benchmarks.jl). +const MadNLPSolverExt = [ + mod for mod in reverse(Base.loaded_modules_order) + if Symbol(mod) == :MadNLPSolverExt +][1] + +# Pull in the bilinear fixture without duplicating it. +include(joinpath(@__DIR__, "..", "test", "test_utils.jl")) + +Random.seed!(42) + +const RESULTS_DIR = joinpath(@__DIR__, "results", "allocs") +mkpath(RESULTS_DIR) + +# ---------------------------------------------------------------------------- +# Problem builder — wraps the shared fixture with a QuadraticRegularizer-style +# objective so both Ipopt and MadNLP see the same NLP. +# ---------------------------------------------------------------------------- +function build_problem(; N = 10) + G, traj = bilinear_dynamics_and_trajectory(; N = N) + + integrators = [ + BilinearIntegrator(G, :x, :u, traj), + DerivativeIntegrator(:u, :du, traj), + DerivativeIntegrator(:du, :ddu, traj), + ] + + J = TerminalObjective(x -> norm(x - traj.goal.x)^2, :x, traj) + J += QuadraticRegularizer(:u, traj, 1.0) + + prob = DirectTrajOptProblem(traj, J, integrators) + return prob, traj +end + +# ---------------------------------------------------------------------------- +# Profile one solver. Warmup runs on a throwaway deepcopy so JIT/compile +# allocations stay out of the recorded trace. +# ---------------------------------------------------------------------------- +function profile_solver(; solver_name, options_ctor, N = 10, sample_rate = 1.0) + prob_warmup, traj = build_problem(; N = N) + prob_profiled, _ = build_problem(; N = N) + + state_dim = traj.dims[:x] + ctrl_dim = sum(traj.dims[cn] for cn in traj.control_names if cn != traj.timestep; init = 0) + + println("\n[$(solver_name)] JIT warmup on throwaway problem copy...") + DirectTrajOpt.solve!(prob_warmup; options = options_ctor()) + + println("[$(solver_name)] Profiling allocations (sample_rate=$(sample_rate))...") + profile = benchmark_memory!( + package = "DirectTrajOpt", + solver = solver_name, + benchmark_name = "bilinear_N$(N)_$(lowercase(solver_name))", + N = traj.N, + state_dim = state_dim, + control_dim = ctrl_dim, + sample_rate = sample_rate, + warmup = false, + runner = "local", + ) do + DirectTrajOpt.solve!(prob_profiled; options = options_ctor()) + end + + mb = profile.total_bytes / (1024 * 1024) + println("[$(solver_name)] captured $(profile.total_count) samples, $(round(mb; digits=2)) MB total") + + path = save_alloc_profile(RESULTS_DIR, profile.benchmark_name, profile) + println("[$(solver_name)] saved to $(path)") + return profile, path +end + +# ---------------------------------------------------------------------------- +# Entry points +# +# sample_rate default is 0.01 because Ipopt/MadNLP generate orders of magnitude +# more fine-grained allocations than the solve's wall-time budget accommodates +# at sample_rate=1.0 (an N=10 bilinear toy can hang for 15+ minutes at 1.0). +# 0.01 still gives statistically useful traces for hot-path triage. +# ---------------------------------------------------------------------------- +function main(; N = 10, sample_rate = 0.01) + ipopt_profile, ipopt_path = profile_solver(; + solver_name = "Ipopt", + options_ctor = () -> IpoptOptions(max_iter = 50, print_level = 0), + N = N, + sample_rate = sample_rate, + ) + + madnlp_profile, madnlp_path = profile_solver(; + solver_name = "MadNLP", + options_ctor = () -> MadNLPSolverExt.MadNLPOptions(max_iter = 50, print_level = Int(MadNLP.ERROR)), + N = N, + sample_rate = sample_rate, + ) + + println("\nDone.") + println(" Ipopt profile: $(ipopt_path) ($(ipopt_profile.total_count) samples)") + println(" MadNLP profile: $(madnlp_path) ($(madnlp_profile.total_count) samples)") + return (ipopt = ipopt_profile, madnlp = madnlp_profile) +end + +if abspath(PROGRAM_FILE) == @__FILE__ + main() +end diff --git a/benchmark/analyze_allocs.jl b/benchmark/analyze_allocs.jl new file mode 100644 index 0000000..9d1ff53 --- /dev/null +++ b/benchmark/analyze_allocs.jl @@ -0,0 +1,136 @@ +using HarmoniqsBenchmarks +using Printf + +const DEFAULT_RESULTS_DIR = joinpath(@__DIR__, "results", "allocs") +results_dir() = isempty(ARGS) ? DEFAULT_RESULTS_DIR : ARGS[1] + +# Noise filters — frames / types from Profile.Allocs itself or the Julia +# toplevel/runtime that do not tell us anything about user-code hotpaths. +const NOISE_FRAME_PATTERNS = [ + "Profile.Allocs", + "gc-alloc-profiler", + "gc-stock.c", + "gc.c:", + "jl_apply", + "jl_toplevel_", + "ijl_toplevel_", + "jl_interpret_toplevel_thunk", + "jl_repl_entrypoint", + "interpreter.c", + "_include(", + "include_string(", + "loading.jl", + "client.jl", + "_start() at sys.so", + "ip:0x", + "_start at ", + " at Base.jl:", + "true_main at jlapi.c", + "__libc_start_main", + "loader_exe.c", + "jl_system_image_data", + "macro expansion at Allocs.jl", + "boot.jl:", + "jl_f__call_latest", +] + +const WRAPPER_FRAME_PATTERNS = [ + "alloc_profile.jl", + "benchmark_memory!", + "HarmoniqsBenchmarks", +] + +const NOISE_TYPE_PATTERNS = [ + "Profile.Allocs", +] + +_is_noise_frame(f) = any(p -> occursin(p, f), NOISE_FRAME_PATTERNS) +_is_noise_type(t) = any(p -> occursin(p, t), NOISE_TYPE_PATTERNS) + +function _first_user_frame(stack) + for f in stack + _is_noise_frame(f) && continue + any(p -> occursin(p, f), WRAPPER_FRAME_PATTERNS) && continue + return f + end + return isempty(stack) ? "" : stack[end] +end + +_is_wrapper_frame(f) = any(p -> occursin(p, f), WRAPPER_FRAME_PATTERNS) + +function top_frames(profile; k = 25, scale_to_total = true, drop_wrappers = true) + by_frame = Dict{String, Tuple{Int, Int}}() + for s in profile.samples + _is_noise_type(s.type_name) && continue + for frame in s.stacktrace + _is_noise_frame(frame) && continue + drop_wrappers && _is_wrapper_frame(frame) && continue + cnt, bytes = get(by_frame, frame, (0, 0)) + by_frame[frame] = (cnt + 1, bytes + s.size_bytes) + end + end + ranked = sort(collect(by_frame); by = x -> -x[2][2])[1:min(k, length(by_frame))] + scale = scale_to_total ? (1 / profile.sample_rate) : 1.0 + println("\nTop $(length(ranked)) user frames by allocated bytes (scaled ×$(Int(scale))):") + println(rpad(" bytes", 14), rpad("samples", 10), "frame") + for (frame, (cnt, bytes)) in ranked + @printf " %-12s %-8d %s\n" _fmt_bytes(bytes * scale) cnt _truncate(frame, 140) + end +end + +function top_leaf_callsites(profile; k = 25, scale_to_total = true) + by_leaf = Dict{String, Tuple{Int, Int}}() + for s in profile.samples + _is_noise_type(s.type_name) && continue + leaf = _first_user_frame(s.stacktrace) + cnt, bytes = get(by_leaf, leaf, (0, 0)) + by_leaf[leaf] = (cnt + 1, bytes + s.size_bytes) + end + ranked = sort(collect(by_leaf); by = x -> -x[2][2])[1:min(k, length(by_leaf))] + scale = scale_to_total ? (1 / profile.sample_rate) : 1.0 + println("\nTop $(length(ranked)) leaf call sites by allocated bytes (scaled ×$(Int(scale))):") + println(rpad(" bytes", 14), rpad("samples", 10), "leaf") + for (leaf, (cnt, bytes)) in ranked + @printf " %-12s %-8d %s\n" _fmt_bytes(bytes * scale) cnt _truncate(leaf, 140) + end +end + +function top_types(profile; k = 15, scale_to_total = true) + by_type = Dict{String, Tuple{Int, Int}}() + for s in profile.samples + _is_noise_type(s.type_name) && continue + cnt, bytes = get(by_type, s.type_name, (0, 0)) + by_type[s.type_name] = (cnt + 1, bytes + s.size_bytes) + end + ranked = sort(collect(by_type); by = x -> -x[2][2])[1:min(k, length(by_type))] + scale = scale_to_total ? (1 / profile.sample_rate) : 1.0 + println("\nTop $(length(ranked)) allocated types (scaled ×$(Int(scale))):") + println(rpad(" bytes", 14), rpad("samples", 10), "type") + for (t, (cnt, bytes)) in ranked + @printf " %-12s %-8d %s\n" _fmt_bytes(bytes * scale) cnt _truncate(t, 120) + end +end + +_fmt_bytes(b) = b >= 1 << 30 ? @sprintf("%.2f GB", b / (1 << 30)) : + b >= 1 << 20 ? @sprintf("%.2f MB", b / (1 << 20)) : + b >= 1 << 10 ? @sprintf("%.2f KB", b / (1 << 10)) : + @sprintf("%d B", Int(round(b))) + +_truncate(s, n) = length(s) <= n ? s : string(first(s, n - 1), "…") + +function main() + dir = results_dir() + files = sort(filter(f -> endswith(f, "_allocs.jld2"), readdir(dir; join = true))) + isempty(files) && (println("no *_allocs.jld2 files under $dir"); return) + for path in files + profile = load_alloc_profile(path) + println("=" ^ 100) + println(basename(path)) + @printf " solver=%s N=%d sample_rate=%g samples=%d total=%s\n" profile.solver profile.N profile.sample_rate profile.total_count _fmt_bytes(profile.total_bytes) + top_types(profile; k = 10) + top_leaf_callsites(profile; k = 20) + top_frames(profile; k = 20) + end +end + +main() diff --git a/ext/MadNLPSolverExt/MadNLPSolverExt.jl b/ext/MadNLPSolverExt/MadNLPSolverExt.jl index e310eaa..4aceef8 100644 --- a/ext/MadNLPSolverExt/MadNLPSolverExt.jl +++ b/ext/MadNLPSolverExt/MadNLPSolverExt.jl @@ -20,4 +20,182 @@ using DirectTrajOpt.Solvers include("solver.jl") include("utils.jl") + +# Coverage targets: ext/MadNLPSolverExt/ + src/solvers/madnlp_solver/ + +@testitem "MadNLPOptions construction" setup=[DTOTestHelpers] begin + opts = DirectTrajOpt.MadNLPOptions() + @test opts.tol == 1e-8 + @test opts.max_iter == 3000 + @test opts.print_level == 3 + @test opts.hessian_approximation == "exact" + + opts2 = DirectTrajOpt.MadNLPOptions(max_iter = 100, tol = 1e-6) + @test opts2.max_iter == 100 + @test opts2.tol == 1e-6 + @test opts isa Solvers.AbstractSolverOptions +end + +@testitem "MadNLP basic solve" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + traj_before = deepcopy(prob.trajectory.data) + solve!(prob; options = DirectTrajOpt.MadNLPOptions(max_iter = 50), verbose = false) + @test prob.trajectory.data != traj_before +end + +@testitem "MadNLP verbose=false" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + output = capture_stdout() do + solve!(prob; options = DirectTrajOpt.MadNLPOptions(max_iter = 10), verbose = false) + end + @test !contains(output, "initializing optimizer") +end + +@testitem "MadNLP verbose=true" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + output = capture_stdout() do + solve!(prob; options = DirectTrajOpt.MadNLPOptions(max_iter = 10), verbose = true) + end + @test contains(output, "initializing optimizer") + @test contains(output, "evaluator created") + @test contains(output, "optimizer initialization complete") +end + +# @testitem "MadNLP unknown kwargs passthrough" setup=[DTOTestHelpers] begin +# # On this branch the @warn for unknown kwargs is commented out in +# # ext/MadNLPSolverExt/solver.jl. Unknown kwargs are silently collected. +# # When warnings are re-enabled, add @test_logs assertion here. +# prob, _ = make_standard_prob() +# solve!(prob; options=DirectTrajOpt.MadNLPOptions(max_iter=5), verbose=false, totally_fake_option=42) +# @test true +# end + +@testitem "MadNLP eval_hessian kwarg routing" setup=[DTOTestHelpers] begin + # eval_hessian=false routes to hessian_approximation="compact_lbfgs". + # The @warn is commented out on this branch — just verify no error. + prob, _ = make_standard_prob() + solve!( + prob; + options = DirectTrajOpt.MadNLPOptions(max_iter = 5), + verbose = false, + eval_hessian = false, + ) + @test true +end + +@testitem "MadNLP eval_hessian kwarg routing" setup=[DTOTestHelpers] begin + # eval_hessian=false routes to hessian_approximation="compact_lbfgs". + # The @warn is commented out on this branch — just verify no error. + prob, _ = make_standard_prob() + result = _solve_with_kwargs( + prob, + DirectTrajOpt.MadNLPOptions(max_iter = 5); + verbose = false, + eval_hessian = false, + ) + @test true +end + +@testitem "MadNLP compact_lbfgs hessian" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + opts = + DirectTrajOpt.MadNLPOptions(max_iter = 10, hessian_approximation = "compact_lbfgs") + solve!(prob; options = opts, verbose = false) + @test true +end + +@testitem "MadNLP with global variables" setup=[DTOTestHelpers] begin + G, traj = bilinear_dynamics_and_trajectory(add_global = true) + integrators = [ + BilinearIntegrator(G, :x, :u, traj), + DerivativeIntegrator(:u, :du, traj), + DerivativeIntegrator(:du, :ddu, traj), + ] + J = TerminalObjective(x -> norm(x - traj.goal.x)^2, :x, traj) + J += QuadraticRegularizer(:u, traj, 1.0) + J += QuadraticRegularizer(:du, traj, 1.0) + J += MinimumTimeObjective(traj) + J += GlobalObjective(g -> norm(g)^2, :g, traj; Q = 1.0) + + g_ug = NonlinearGlobalKnotPointConstraint( + ug -> begin + u = ug[1:traj.dims[:u]] + g = ug[(traj.dims[:u]+1):end] + return [norm(u) * (1.0 + norm(g)) - 2.0] + end, + [:u], + [:g], + traj; + times = 2:(traj.N-1), + equality = false, + ) + prob = + DirectTrajOptProblem(traj, J, integrators; constraints = AbstractConstraint[g_ug]) + solve!(prob; options = DirectTrajOpt.MadNLPOptions(max_iter = 50), verbose = false) + + for k = 2:(traj.N-1) + u = traj[k][:u] + g = traj.global_data[traj.global_components[:g]] + @test norm(u) * (1.0 + norm(g)) <= 2.0 + 1e-5 + end +end + +@testitem "_solve_with_kwargs with MumpsSolver (default)" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + DirectTrajOpt._solve_with_kwargs( + prob, + DirectTrajOpt.MadNLPOptions(max_iter = 50); + verbose = false, + kkt_system = MadNLP.SparseKKTSystem, + linear_solver = MadNLP.MumpsSolver, + ) + @test true +end + +@testitem "_solve_with_kwargs with LapackCPUSolver" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + DirectTrajOpt._solve_with_kwargs( + prob, + DirectTrajOpt.MadNLPOptions(max_iter = 50); + verbose = false, + kkt_system = MadNLP.SparseUnreducedKKTSystem, + linear_solver = MadNLP.LapackCPUSolver, + ) + @test true +end + +@testitem "_solve_with_kwargs with LOQOUpdate adaptive barrier" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + # LOQOUpdate: adaptive barrier from Nocedal et al. 2009 §3. + # Uses average/min complementarity ratio to set the barrier parameter. + # Falls back to monotone if insufficient progress. + DirectTrajOpt._solve_with_kwargs( + prob, + DirectTrajOpt.MadNLPOptions(max_iter = 50); + verbose = false, + kkt_system = MadNLP.SparseKKTSystem, + linear_solver = MadNLP.MumpsSolver, + barrier = MadNLP.LOQOUpdate(1e-8, 10.0), + ) + @test true +end + +@testitem "_solve_with_kwargs with QualityFunctionUpdate adaptive barrier" setup=[ + DTOTestHelpers, +] begin + prob, _ = make_standard_prob() + # QualityFunctionUpdate: adaptive barrier from Nocedal et al. 2009 §4. + # Minimizes an ℓ1 quality function via golden search; falls back to + # monotone if insufficient progress. + DirectTrajOpt._solve_with_kwargs( + prob, + DirectTrajOpt.MadNLPOptions(max_iter = 50); + verbose = false, + kkt_system = MadNLP.SparseKKTSystem, + linear_solver = MadNLP.MumpsSolver, + barrier = MadNLP.QualityFunctionUpdate(1e-8, 10.0), + ) + @test true +end + end diff --git a/ext/MadNLPSolverExt/solver.jl b/ext/MadNLPSolverExt/solver.jl index 514d29a..ca60e53 100644 --- a/ext/MadNLPSolverExt/solver.jl +++ b/ext/MadNLPSolverExt/solver.jl @@ -208,7 +208,12 @@ function DirectTrajOpt.set_options!(optimizer::AbstractOptimizer, options::MadNL if name in ignored_options continue end - # TODO: allow internal defaults, i.e. do not set the internal options dict unless the user actually specified the associated opt + # `nothing` means "use MadNLP's own default" — don't overwrite the optimizer's + # internal dict in that case. Applies to the pass-through fields + # (linear_solver, array_type, kkt_system, cudss_ordering). + if value === nothing + continue + end if name == :print_level optimizer.options[name] = MadNLP.LogLevels(value) elseif name == :hessian_approximation diff --git a/ext/MadNLPSolverExt/utils.jl b/ext/MadNLPSolverExt/utils.jl index e57f119..4cd983a 100644 --- a/ext/MadNLPSolverExt/utils.jl +++ b/ext/MadNLPSolverExt/utils.jl @@ -104,7 +104,9 @@ function DirectTrajOpt._solve_with_kwargs( # TODO: Verify this is working as expected update_trajectory!(prob, optimizer, variables) - return nothing + # TODO: Make this standard? + # return nothing + return optimizer end diff --git a/scripts/diagnose_phase_bifurcation.jl b/scripts/diagnose_phase_bifurcation.jl new file mode 100644 index 0000000..d2c39a2 --- /dev/null +++ b/scripts/diagnose_phase_bifurcation.jl @@ -0,0 +1,201 @@ +# Diagnose global phase bifurcation between Ipopt and MadNLP +# +# Run this AFTER the main comparison script has populated: +# qcp_ipopt, qcp_madnlp, integrator, qtraj +# +# Or run standalone — it sets everything up. + +import Ipopt +import MadNLP + +using DirectTrajOpt +using Piccolo +using Piccolissimo +using LinearAlgebra + +function _get_MadNLPSolverExt() + cur_mods = reverse(Base.loaded_modules_order) + candidate_mods = [n for n = cur_mods if Symbol(n) == :MadNLPSolverExt] + length(candidate_mods) == 1 || error("Expected 1 MadNLPSolverExt, found $(length(candidate_mods))") + return candidate_mods[1] +end +const MadNLPSolverExt = _get_MadNLPSolverExt() + +# ---------------------------------------------------------------------------- # +# Setup: same problem for both solvers +# ---------------------------------------------------------------------------- # + +H_drift = 0.5 * PAULIS.Z +H_drives = [PAULIS.X, PAULIS.Y] +drive_bounds = [1., 1.] +sys = QuantumSystem(H_drift, H_drives, drive_bounds) +T = 10.0 +N = 100 +U_goal = GATES.X + +qtraj = UnitaryTrajectory(sys, U_goal, T) +integrator = HermitianExponentialIntegrator(qtraj, N) + +# Build two identical problems +qcp = SmoothPulseProblem(qtraj, N; Q=1e2, R=1e-2, ddu_bound=1e0, Δt_bounds=(0.5, 1.5), integrator=integrator) +qcp_ipopt = deepcopy(qcp) +qcp_madnlp = deepcopy(qcp) + +# ---------------------------------------------------------------------------- # +# Solve +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("Solving with Ipopt...") +println("="^70) +solve!(qcp_ipopt; options=IpoptOptions(max_iter=100)) + +println("\n" * "="^70) +println("Solving with MadNLP...") +println("="^70) +solve!(qcp_madnlp; options=MadNLPSolverExt.MadNLPOptions(max_iter=100)) + +# ---------------------------------------------------------------------------- # +# Extract final unitaries from the solver trajectory data +# ---------------------------------------------------------------------------- # + +function get_final_unitary(qcp) + traj = get_trajectory(qcp) + # The state Ũ⃗ at the final knot point + Ũ⃗_final = traj[:Ũ⃗][:, end] + return iso_vec_to_operator(Ũ⃗_final) +end + +U_ipopt = get_final_unitary(qcp_ipopt) +U_madnlp = get_final_unitary(qcp_madnlp) + +# ---------------------------------------------------------------------------- # +# Diagnose: what phase did each solver converge to? +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("PHASE ANALYSIS") +println("="^70) + +# The trace inner product tr(U_goal' * U_final) +# For a perfect solution: tr(X' * (α·X)) = α·tr(I) = 2α +# where α = ±i (from Schrödinger evolution: U = exp(-iHt)) +trace_ipopt = tr(U_goal' * U_ipopt) +trace_madnlp = tr(U_goal' * U_madnlp) + +println("\ntr(U_goal' * U_final):") +println(" Ipopt: $trace_ipopt") +println(" MadNLP: $trace_madnlp") + +println("\nPhase angle (arg of trace):") +println(" Ipopt: $(angle(trace_ipopt)) rad ($(rad2deg(angle(trace_ipopt)))°)") +println(" MadNLP: $(angle(trace_madnlp)) rad ($(rad2deg(angle(trace_madnlp)))°)") + +println("\nPhase sign diagnostic (im * tr / |tr|):") +phase_ipopt = trace_ipopt / abs(trace_ipopt) +phase_madnlp = trace_madnlp / abs(trace_madnlp) +println(" Ipopt: $phase_ipopt") +println(" MadNLP: $phase_madnlp") + +println("\nFidelity (abs2-based, phase-insensitive):") +n = size(U_goal, 1) +fid_ipopt = abs2(trace_ipopt) / n^2 +fid_madnlp = abs2(trace_madnlp) / n^2 +println(" Ipopt: $fid_ipopt") +println(" MadNLP: $fid_madnlp") + +println("\nFidelity via Piccolo API:") +println(" Ipopt: $(fidelity(qcp_ipopt))") +println(" MadNLP: $(fidelity(qcp_madnlp))") + +# ---------------------------------------------------------------------------- # +# Check if solvers converged to different phase basins +# ---------------------------------------------------------------------------- # + +phase_diff = angle(trace_ipopt) - angle(trace_madnlp) +# Normalize to [-π, π] +phase_diff = mod(phase_diff + π, 2π) - π + +println("\n" * "="^70) +println("PHASE DIFFERENCE BETWEEN SOLVERS") +println("="^70) +println(" Δφ = $(phase_diff) rad ($(rad2deg(phase_diff))°)") +if abs(phase_diff) > π/2 + println(" ⚠ SOLVERS CONVERGED TO DIFFERENT PHASE BASINS") + println(" This is the bifurcation Gennadi identified.") + println(" The abs2 in the objective makes ±(im·U_goal) equivalent.") +else + println(" ✓ Solvers converged to the same phase basin this time.") + println(" Re-run to potentially observe bifurcation (depends on numerics).") +end + +# ---------------------------------------------------------------------------- # +# Show the actual final unitaries for inspection +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("FINAL UNITARIES") +println("="^70) + +println("\nU_goal = X gate:") +display(U_goal) + +println("\n\nU_final (Ipopt):") +display(U_ipopt) + +println("\n\nU_final (MadNLP):") +display(U_madnlp) + +println("\n\nU_final / (im * U_goal) — should be ≈ ±1 * I if converged:") +println("\n Ipopt: U / (i·X) =") +display(U_ipopt / (im * U_goal)) +println("\n MadNLP: U / (i·X) =") +display(U_madnlp / (im * U_goal)) + +# ---------------------------------------------------------------------------- # +# Constraint violations (for completeness) +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("CONSTRAINT VIOLATIONS") +println("="^70) + +deltas_ipopt = zeros(integrator.dim) +deltas_madnlp = zeros(integrator.dim) +DirectTrajOpt.evaluate!(deltas_ipopt, integrator, get_trajectory(qcp_ipopt)) +DirectTrajOpt.evaluate!(deltas_madnlp, integrator, get_trajectory(qcp_madnlp)) + +println(" max |constraint| Ipopt: $(maximum(abs.(deltas_ipopt)))") +println(" max |constraint| MadNLP: $(maximum(abs.(deltas_madnlp)))") + +# ---------------------------------------------------------------------------- # +# Datavec comparison +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("TRAJECTORY DIVERGENCE") +println("="^70) + +data_ipopt = get_trajectory(qcp_ipopt).data +data_madnlp = get_trajectory(qcp_madnlp).data +max_diff = maximum(abs.(data_ipopt .- data_madnlp)) +println(" max |data_ipopt - data_madnlp| = $max_diff") + +println("\n" * "="^70) +println("DIAGNOSIS COMPLETE") +println("="^70) +println(""" +The objective `abs2(tr(U_goal' * U)) / n^2` at: + qc2/src/quantum_objectives.jl:49 + +is phase-insensitive: abs2(z) = abs2(-z) = abs2(iz) = ... +Both +im*U_goal and -im*U_goal are equally valid optima. + +When Δt is a free variable, the expanded search space creates +multiple basins of attraction. The solver's early numerical +trajectory (from initialization + step selection) determines +which basin it falls into. + +With fixed Δt, the landscape is more constrained and both +solvers consistently find the same basin. +""") diff --git a/src/constraints/linear/equality_constraint.jl b/src/constraints/linear/equality_constraint.jl index b1768c4..9853b73 100644 --- a/src/constraints/linear/equality_constraint.jl +++ b/src/constraints/linear/equality_constraint.jl @@ -143,6 +143,47 @@ function fix_trajectory_variable!( return constraints end +export fix_global_variable! + +""" + fix_global_variable!(constraints, name, value) + +Pin a global (time-invariant) variable to `value` using a `GlobalEqualityConstraint`. +Removes any existing `BoundsConstraint` or `EqualityConstraint` on the same global +variable to avoid MOI conflicts. Companion to [`fix_trajectory_variable!`](@ref) for +the global-variable case. + +This is the integrator-agnostic mechanism for pinning a calibrated parameter +(e.g. learned `θ` in QILC alternating calibration) into the control NLP — any +globals-aware integrator (`HermitianExponentialIntegrator`, `SplineIntegrator`, +`NonHermitianExponentialIntegrator`, ...) will read the pinned value through the +NLP variable rather than requiring an integrator rebuild. + +# Arguments +- `constraints::Vector{<:AbstractConstraint}`: mutable constraint list +- `name::Symbol`: global variable name to pin +- `value::AbstractVector{Float64}`: pin value (length must equal the global variable dim) +""" +function fix_global_variable!( + constraints::Vector{<:AbstractConstraint}, + name::Symbol, + value::AbstractVector{Float64}, +) + # Remove any existing BoundsConstraint and EqualityConstraint on this global variable + # (per-timestep dedup at line ~134 only filters non-global constraints; this is the + # globals counterpart). Pinning supersedes any prior bounds or pin on the same global. + filter!(constraints) do c + if c isa BoundsConstraint && c.var_names == name && c.is_global + return false + elseif c isa EqualityConstraint && c.var_names == name && c.is_global + return false + end + return true + end + push!(constraints, GlobalEqualityConstraint(name, collect(Float64, value))) + return constraints +end + function Base.show(io::IO, c::EqualityConstraint) print(io, "EqualityConstraint: \"$(c.label)\"") end @@ -300,3 +341,54 @@ end @test prob.trajectory[t][:u] ≈ u_ref[:, t] atol=1e-10 end end + +@testitem "fix_global_variable! pins global to supplied value" begin + include("../../../test/test_utils.jl") + + G, traj = bilinear_dynamics_and_trajectory(add_global = true) + + integrators = [ + BilinearIntegrator(G, :x, :u, traj), + DerivativeIntegrator(:u, :du, traj), + DerivativeIntegrator(:du, :ddu, traj), + ] + + J = TerminalObjective(x -> norm(x - traj.goal.x)^2, :x, traj) + J += QuadraticRegularizer(:u, traj, 1.0) + J += MinimumTimeObjective(traj) + + # Build a normal problem first (so bounds/equality on globals may exist) + prob_orig = DirectTrajOptProblem(traj, J, integrators) + + # Pin global :g to a specific value via fix_global_variable! + g_dim = length(traj.global_components[:g]) + g_value = fill(0.42, g_dim) + constraints = deepcopy(prob_orig.constraints) + fix_global_variable!(constraints, :g, g_value) + + # After fixing, no remaining BoundsConstraint or EqualityConstraint on :g + # other than our new GlobalEqualityConstraint + g_eq_cons = filter( + c -> c isa EqualityConstraint && c.var_names == :g && c.is_global, + constraints, + ) + @test length(g_eq_cons) == 1 + @test !any(c -> c isa BoundsConstraint && c.var_names == :g && c.is_global, constraints) + + # Re-applying should still leave exactly one equality constraint (dedup works) + fix_global_variable!(constraints, :g, fill(0.5, g_dim)) + g_eq_cons2 = filter( + c -> c isa EqualityConstraint && c.var_names == :g && c.is_global, + constraints, + ) + @test length(g_eq_cons2) == 1 + @test g_eq_cons2[1].values ≈ fill(0.5, g_dim) + + # Solve and confirm the global is pinned at the supplied value + fix_global_variable!(constraints, :g, g_value) # pin at 0.42 again + prob = DirectTrajOptProblem(traj, J, integrators, constraints) + solve!(prob; max_iter = 100) + + g_components = traj.global_components[:g] + @test prob.trajectory.global_data[g_components] ≈ g_value atol=1e-6 +end diff --git a/src/solvers/_solvers.jl b/src/solvers/_solvers.jl index 2407bc7..e81bb2d 100644 --- a/src/solvers/_solvers.jl +++ b/src/solvers/_solvers.jl @@ -32,4 +32,31 @@ include("constrain.jl") include("evaluator.jl") include("solve.jl") + +# Coverage targets: src/solvers/_solvers.jl (50% → ~100%) + +@testitem "DefaultSolverOptions get/set" setup=[DTOTestHelpers] begin + original = Solvers._get_DefaultSolverOptions() + @test original == IpoptSolverExt.IpoptOptions + + Solvers._set_DefaultSolverOptions(DirectTrajOpt.MadNLPOptions) + @test Solvers._get_DefaultSolverOptions() == DirectTrajOpt.MadNLPOptions + + Solvers._set_DefaultSolverOptions(Solvers.DefaultSolverOptions) + @test Solvers._get_DefaultSolverOptions() == Solvers.DefaultSolverOptions + + Solvers._set_DefaultSolverOptions(original) + @test Solvers._get_DefaultSolverOptions() == IpoptSolverExt.IpoptOptions +end + +@testitem "AbstractOptimizer alias" setup=[DTOTestHelpers] begin + @test Solvers.AbstractOptimizer === MOI.AbstractOptimizer +end + +@testitem "AbstractSolverOptions hierarchy" setup=[DTOTestHelpers] begin + @test IpoptSolverExt.IpoptOptions <: Solvers.AbstractSolverOptions + @test DirectTrajOpt.MadNLPOptions <: Solvers.AbstractSolverOptions + @test Solvers.DefaultSolverOptions <: Solvers.AbstractSolverOptions +end + end diff --git a/src/solvers/constrain.jl b/src/solvers/constrain.jl index 0571f77..05d500d 100644 --- a/src/solvers/constrain.jl +++ b/src/solvers/constrain.jl @@ -351,3 +351,75 @@ function (con::TimeConsistencyConstraint)( ) end end + + +# Coverage targets: src/solvers/constrain.jl (92% → ~95%) + cross-solver agreement + +@testitem "constrain! verbose output" setup=[DTOTestHelpers] begin + G, traj = bilinear_dynamics_and_trajectory() + integrators = [ + BilinearIntegrator(G, :x, :u, traj), + DerivativeIntegrator(:u, :du, traj), + DerivativeIntegrator(:du, :ddu, traj), + ] + J = QuadraticRegularizer(:u, traj, 1.0) + prob = DirectTrajOptProblem(traj, J, integrators) + + evaluator = Solvers.Evaluator(prob; eval_hessian = true, verbose = false) + nl_cons = Solvers.get_nonlinear_constraints(prob) + block_data = MOI.NLPBlockData(nl_cons, evaluator, true) + + optimizer = Ipopt.Optimizer() + MOI.set(optimizer, MOI.NLPBlock(), block_data) + MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE) + + data_dim = traj.dim * traj.N + variables = MOI.add_variables(optimizer, data_dim + traj.global_dim) + MOI.set( + optimizer, + MOI.VariablePrimalStart(), + variables[1:data_dim], + collect(traj.datavec), + ) + + linear_constraints = AbstractLinearConstraint[filter( + c->c isa AbstractLinearConstraint, + prob.constraints, + )...] + + output = capture_stdout() do + Solvers.constrain!( + optimizer, + variables, + linear_constraints, + prob.trajectory; + verbose = true, + ) + end + @test contains(output, "applying constraint") +end + +@testitem "Cross-solver agreement (Ipopt vs MadNLP)" setup=[DTOTestHelpers] begin + using Random + + include( + joinpath(dirname(dirname(pathof(DirectTrajOpt))), "test", "solver_test_utils.jl"), + ) + + seed = UInt64(42) + + prob_ipopt = get_seeded_prob_solved( + seed, + IpoptSolverExt.IpoptOptions(; max_iter = 100, print_level = 0), + ) + prob_madnlp = + get_seeded_prob_solved(seed, DirectTrajOpt.MadNLPOptions(; max_iter = 100)) + + traj_ipopt = prob_ipopt.trajectory + traj_madnlp = prob_madnlp.trajectory + + traj_dist = (traj_madnlp.data[:, :] .- traj_ipopt.data[:, :]) .^ 2 + traj_dist = sqrt(sum(traj_dist)) / length(traj_dist) + + @test traj_dist < 1e-4 +end diff --git a/src/solvers/evaluator.jl b/src/solvers/evaluator.jl index 8c59bab..5989f0c 100644 --- a/src/solvers/evaluator.jl +++ b/src/solvers/evaluator.jl @@ -769,3 +769,84 @@ end ) @test all(isapprox.(triu(∂²ℒ), triu(sparse(∂²ℒ_finitediff)), atol = 1e-2)) end + + +# Coverage targets: src/solvers/evaluator.jl (61% → ~90%) + +@testitem "features_available branches" setup=[DTOTestHelpers] begin + prob, _, eval_hess = make_evaluator(eval_hessian = true) + @test :Grad in MOI.features_available(eval_hess) + @test :Jac in MOI.features_available(eval_hess) + @test :Hess in MOI.features_available(eval_hess) + + _, _, eval_no_hess = make_evaluator(eval_hessian = false) + @test :Grad in MOI.features_available(eval_no_hess) + @test :Jac in MOI.features_available(eval_no_hess) + @test :Hess ∉ MOI.features_available(eval_no_hess) +end + +@testitem "eval_constraint_jacobian_product" setup=[DTOTestHelpers] begin + _, traj, evaluator = make_evaluator() + + Z⃗ = collect(traj.datavec) + n_vars = length(Z⃗) + n_cons = evaluator.n_constraints + + # Build dense Jacobian from sparse structure + jac_structure = MOI.jacobian_structure(evaluator) + jac_values = zeros(length(jac_structure)) + MOI.eval_constraint_jacobian(evaluator, jac_values, Z⃗) + + J_dense = zeros(n_cons, n_vars) + for (idx, (row, col)) in enumerate(jac_structure) + J_dense[row, col] += jac_values[idx] + end + + w = randn(n_vars) + y_moi = zeros(n_cons) + MOI.eval_constraint_jacobian_product(evaluator, y_moi, Z⃗, w) + + @test isapprox(y_moi, J_dense * w, atol = 1e-10) +end + +@testitem "eval_constraint_jacobian_transpose_product" setup=[DTOTestHelpers] begin + _, traj, evaluator = make_evaluator() + + Z⃗ = collect(traj.datavec) + n_vars = length(Z⃗) + n_cons = evaluator.n_constraints + + jac_structure = MOI.jacobian_structure(evaluator) + jac_values = zeros(length(jac_structure)) + MOI.eval_constraint_jacobian(evaluator, jac_values, Z⃗) + + J_dense = zeros(n_cons, n_vars) + for (idx, (row, col)) in enumerate(jac_structure) + J_dense[row, col] += jac_values[idx] + end + + w = randn(n_cons) + y_moi = zeros(n_vars) + MOI.eval_constraint_jacobian_transpose_product(evaluator, y_moi, Z⃗, w) + + @test isapprox(y_moi, J_dense' * w, atol = 1e-10) +end + +@testitem "Evaluator verbose construction" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + output = capture_stdout() do + Solvers.Evaluator(prob; eval_hessian = true, verbose = true) + end + @test contains(output, "building evaluator") + @test contains(output, "evaluator ready") + @test contains(output, "jacobian structure") + @test contains(output, "hessian structure") +end + +@testitem "Evaluator silent construction" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + output = capture_stdout() do + Solvers.Evaluator(prob; eval_hessian = true, verbose = false) + end + @test output == "" +end diff --git a/src/solvers/ipopt_solver/callbacks.jl b/src/solvers/ipopt_solver/callbacks.jl index 6cb62a3..98d7829 100644 --- a/src/solvers/ipopt_solver/callbacks.jl +++ b/src/solvers/ipopt_solver/callbacks.jl @@ -530,4 +530,136 @@ end @test Callbacks.test_early_stop(deepcopy(prob), 50, 100, 25) # [50 < 100] < 25 end + +# Coverage targets: src/solvers/ipopt_solver/callbacks.jl (59% → ~85%) + +@testitem "Callback trajectory update" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + @test Callbacks.test_update_trajectory(deepcopy(prob), true) + @test Callbacks.test_update_trajectory(deepcopy(prob), false) +end + +@testitem "Callback early stop" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + @test Callbacks.test_early_stop(deepcopy(prob), 50, 100, 25) # stop < min < max + @test Callbacks.test_early_stop(deepcopy(prob), 50, 100, 75) # min < stop < max + @test Callbacks.test_early_stop(deepcopy(prob), 50, 100, 125) # min < max < stop +end + +@testitem "Callback trajectory history" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + @test Callbacks.test_update_trajectory_history(deepcopy(prob); n_iters = 10) +end + +@testitem "callback_rollout_fidelity_factory" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + + call_count = Ref(0) + mock_fid_fn = (traj, sys) -> begin + call_count[] += 1 + return 0.5 + 0.01 * call_count[] + end + + fidelities = Dict{Int32,Float64}() + callback = Callbacks.callback_factory( + Callbacks.callback_update_trajectory_factory(prob), + Callbacks.callback_rollout_fidelity_factory( + prob, + nothing, + mock_fid_fn, + fidelities; + freq = 2, + fid_thresh = nothing, + ), + Callbacks.callback_stop_iteration_factory(10), + ) + + optimizer, variables = IpoptSolverExt.get_optimizer_and_variables( + prob, + IpoptOptions(; max_iter = 20, print_level = 0), + callback, + ) + IpoptSolverExt.MOI.optimize!(optimizer) + + @test length(fidelities) > 0 +end + +@testitem "callback_best_rollout_fidelity_factory" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + + call_count = Ref(0) + mock_fid_fn = (traj, sys) -> begin + call_count[] += 1 + return 0.1 * call_count[] + end + + trajectories = Dict{Int32,Any}() + callback = Callbacks.callback_factory( + Callbacks.callback_update_trajectory_factory(prob), + Callbacks.callback_best_rollout_fidelity_factory( + prob, + nothing, + mock_fid_fn, + trajectories; + max_trajectories = 3, + freq = 1, + fid_thresh = nothing, + ), + Callbacks.callback_stop_iteration_factory(8), + ) + + optimizer, variables = IpoptSolverExt.get_optimizer_and_variables( + prob, + IpoptOptions(; max_iter = 20, print_level = 0), + callback, + ) + IpoptSolverExt.MOI.optimize!(optimizer) + + @test 0 < length(trajectories) <= 3 + for (k, (fid, t)) in trajectories + @test fid isa Number + @test t isa NamedTrajectory + end +end + +@testitem "callback_say_hello_factory" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + + output = capture_stdout() do + callback = Callbacks.callback_factory( + Callbacks.callback_say_hello_factory("TEST_MSG_12345"), + Callbacks.callback_stop_iteration_factory(2), + ) + optimizer, variables = IpoptSolverExt.get_optimizer_and_variables( + prob, + IpoptOptions(; max_iter = 5, print_level = 0), + callback, + ) + IpoptSolverExt.MOI.optimize!(optimizer) + end + @test contains(output, "TEST_MSG_12345") +end + +@testitem "callback_update_optimizer_state_history" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + stats = Callbacks.IpoptOptimizerState[] + + callback = Callbacks.callback_factory( + Callbacks.callback_update_optimizer_state_history_factory(stats), + Callbacks.callback_stop_iteration_factory(5), + ) + + optimizer, variables = IpoptSolverExt.get_optimizer_and_variables( + prob, + IpoptOptions(; max_iter = 20, print_level = 0), + callback, + ) + IpoptSolverExt.MOI.optimize!(optimizer) + + @test length(stats) > 0 + @test haskey(stats[1], :iter_count) + @test haskey(stats[1], :obj_value) +end + + end diff --git a/src/solvers/madnlp_solver/options.jl b/src/solvers/madnlp_solver/options.jl index 9a6fb73..6d3a382 100644 --- a/src/solvers/madnlp_solver/options.jl +++ b/src/solvers/madnlp_solver/options.jl @@ -7,6 +7,13 @@ export MadNLPOptions print_level::Int = 3 # (MadNLP.TRACE::MadNLP.LogLevels = 1, ..., MadNLP.ERROR::MadNLP.LogLevels = 6) hessian_approximation::String = "exact" # (exact = MadNLP.ExactHessian, compact_lbfgs = MadNLP.CompactLBFGS) # no other QN methods supported in conjunction with MadNLP.SparseCallback + # Pass-throughs consumed by MadNLP's MOI layer (not by MadNLP itself); + # leave as `nothing` to use MadNLP defaults. Only forwarded when non-nothing. + linear_solver::Any = nothing # e.g. MadNLPGPU.CUDSSSolver, MadNLP.LapackCPUSolver + array_type::Any = nothing # e.g. CUDA.CuArray for GPU + kkt_system::Any = nothing # e.g. MadNLP.SparseUnreducedKKTSystem + cudss_ordering::Any = nothing # e.g. MadNLPGPU.AMD_ORDERING + # # Only supported by DirectTrajOpt._solve, as an optional kwarg override of `hessian_approximation`; # # `hessian_approximation = eval_hessian ? "exact" : "compact_lbfgs"` # eval_hessian::Bool = true diff --git a/src/solvers/solve.jl b/src/solvers/solve.jl index 95900fa..732dec1 100644 --- a/src/solvers/solve.jl +++ b/src/solvers/solve.jl @@ -203,3 +203,115 @@ function DirectTrajOpt.solve!( return nothing end + + +# Coverage targets: src/solvers/solve.jl (35% → ~80%) + +@testitem "get_num_variables without globals" setup=[DTOTestHelpers] begin + G, traj = bilinear_dynamics_and_trajectory() + integrators = [BilinearIntegrator(G, :x, :u, traj)] + J = QuadraticRegularizer(:u, traj, 1.0) + prob = DirectTrajOptProblem(traj, J, integrators) + + @test traj.global_dim == 0 + @test Solvers.get_num_variables(prob) == traj.dim * traj.N +end + +@testitem "get_num_variables with globals" setup=[DTOTestHelpers] begin + G, traj = bilinear_dynamics_and_trajectory(add_global = true) + integrators = [BilinearIntegrator(G, :x, :u, traj)] + J = QuadraticRegularizer(:u, traj, 1.0) + prob = DirectTrajOptProblem(traj, J, integrators) + + @test traj.global_dim > 0 + @test Solvers.get_num_variables(prob) == traj.dim * traj.N + traj.global_dim +end + +@testitem "get_nonlinear_constraints — dynamics only" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + # Remove the nonlinear constraint to test dynamics-only path + G, traj = bilinear_dynamics_and_trajectory() + integrators = [ + BilinearIntegrator(G, :x, :u, traj), + DerivativeIntegrator(:u, :du, traj), + DerivativeIntegrator(:du, :ddu, traj), + ] + J = QuadraticRegularizer(:u, traj, 1.0) + prob_dyn = DirectTrajOptProblem(traj, J, integrators) + nl_cons = Solvers.get_nonlinear_constraints(prob_dyn) + + # All dynamics → equality bounds (0, 0) + @test all(c -> c.lower == 0.0 && c.upper == 0.0, nl_cons) +end + +@testitem "get_nonlinear_constraints — inequality" setup=[DTOTestHelpers] begin + G, traj = bilinear_dynamics_and_trajectory() + integrators = [ + BilinearIntegrator(G, :x, :u, traj), + DerivativeIntegrator(:u, :du, traj), + DerivativeIntegrator(:du, :ddu, traj), + ] + J = QuadraticRegularizer(:u, traj, 1.0) + g_ineq = NonlinearKnotPointConstraint( + u -> [norm(u) - 1.0], + :u, + traj; + times = 2:(traj.N-1), + equality = false, + ) + prob = + DirectTrajOptProblem(traj, J, integrators; constraints = AbstractConstraint[g_ineq]) + nl_cons = Solvers.get_nonlinear_constraints(prob) + + dynamics_dim = sum(i.dim for i in integrators) + ineq_cons = nl_cons[(dynamics_dim+1):end] + @test all(c -> c.lower == -Inf && c.upper == 0.0, ineq_cons) +end + +@testitem "get_nonlinear_constraints — equality" setup=[DTOTestHelpers] begin + G, traj = bilinear_dynamics_and_trajectory() + integrators = [ + BilinearIntegrator(G, :x, :u, traj), + DerivativeIntegrator(:u, :du, traj), + DerivativeIntegrator(:du, :ddu, traj), + ] + J = QuadraticRegularizer(:u, traj, 1.0) + g_eq = NonlinearKnotPointConstraint( + u -> [norm(u) - 0.5], + :u, + traj; + times = 2:(traj.N-1), + equality = true, + ) + prob = + DirectTrajOptProblem(traj, J, integrators; constraints = AbstractConstraint[g_eq]) + nl_cons = Solvers.get_nonlinear_constraints(prob) + + dynamics_dim = sum(i.dim for i in integrators) + eq_cons = nl_cons[(dynamics_dim+1):end] + @test all(c -> c.lower == 0.0 && c.upper == 0.0, eq_cons) +end + +@testitem "_solve error stubs" setup=[DTOTestHelpers] begin + G, traj = bilinear_dynamics_and_trajectory() + integrators = [BilinearIntegrator(G, :x, :u, traj)] + J = QuadraticRegularizer(:u, traj, 1.0) + prob = DirectTrajOptProblem(traj, J, integrators) + + # Ensure CI knows these are meant to error; maybe wrap in a try block? + + # Non-AbstractSolverOptions type + @test DirectTrajOpt._solve(prob, "not_options") === nothing + + # AbstractSolverOptions subtype with no backend + struct _FakeSolverOptions <: Solvers.AbstractSolverOptions end + @test DirectTrajOpt._solve(prob, _FakeSolverOptions()) === nothing +end + +@testitem "solve! default dispatch uses Ipopt" setup=[DTOTestHelpers] begin + prob, _ = make_standard_prob() + @test Solvers._get_DefaultSolverOptions() == IpoptSolverExt.IpoptOptions + traj_before = deepcopy(prob.trajectory.data) + solve!(prob; max_iter = 2, print_level = 0, verbose = false) + @test prob.trajectory.data != traj_before +end diff --git a/test/runtests.jl b/test/runtests.jl index 9f95075..bebf5d7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using DirectTrajOpt using TestItemRunner +include("test_snippets.jl") # Run all testitem tests in package @run_package_tests diff --git a/test/test_snippets.jl b/test/test_snippets.jl new file mode 100644 index 0000000..8c125bd --- /dev/null +++ b/test/test_snippets.jl @@ -0,0 +1,73 @@ +using TestItemRunner +using TestItems + +# ============================================================================ +# Shared test setup via @testsnippet and @testmodule +# ============================================================================ + +# @testsnippet: inlined into every @testitem that lists it in setup=[...]. +# Use for lightweight helper definitions (functions, constants). +# Re-evaluated per testitem — no shared mutable state. + +@testsnippet DTOTestHelpers begin + using DirectTrajOpt + using DirectTrajOpt: Solvers, IpoptSolverExt, Callbacks + import MathOptInterface as MOI + import MadNLP + import Ipopt + using LinearAlgebra + using SparseArrays + using NamedTrajectories + using Test + + # Pull in bilinear_dynamics_and_trajectory and friends from DTO's own test utils + include(joinpath(dirname(dirname(pathof(DirectTrajOpt))), "test", "test_utils.jl")) + + # Standard problem builder used by most tests + function make_standard_prob(; add_global = false) + G, traj = bilinear_dynamics_and_trajectory(; add_global = add_global) + integrators = [ + BilinearIntegrator(G, :x, :u, traj), + DerivativeIntegrator(:u, :du, traj), + DerivativeIntegrator(:du, :ddu, traj), + ] + J = TerminalObjective(x -> norm(x - traj.goal.x)^2, :x, traj) + J += QuadraticRegularizer(:u, traj, 1.0) + J += QuadraticRegularizer(:du, traj, 1.0) + J += MinimumTimeObjective(traj) + + g_u_norm = NonlinearKnotPointConstraint( + u -> [norm(u) - 1.0], + :u, + traj; + times = 2:(traj.N-1), + equality = false, + ) + + prob = DirectTrajOptProblem( + traj, + J, + integrators; + constraints = AbstractConstraint[g_u_norm], + ) + return prob, traj + end + + # Standard problem + evaluator builder + function make_evaluator(; eval_hessian = true) + prob, traj = make_standard_prob() + evaluator = Solvers.Evaluator(prob; eval_hessian = eval_hessian, verbose = false) + return prob, traj, evaluator + end + + # Pipe-based stdout capture for Julia 1.12+. + # Standalone-driver workaround — drop when a better mechanism is available. + function capture_stdout(f) + pipe = Pipe() + redirect_stdout(pipe) do + f() + end + close(pipe.in) + return read(pipe.out, String) + end +end