From eae4404849f8682a41c9b8d4d9fd68ae9a79a987 Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Sun, 31 May 2026 17:31:41 -0400 Subject: [PATCH 1/3] reaname args --- julia/src/Marble.jl | 8 ++++---- julia/src/conversion.jl | 14 +++++++------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/julia/src/Marble.jl b/julia/src/Marble.jl index a15a4ff..d7523bd 100644 --- a/julia/src/Marble.jl +++ b/julia/src/Marble.jl @@ -15,16 +15,16 @@ module Marble @initcxx end - function setup!(solver::Marble.Solver,model::Model, ind_cc1, ind_cc2, comp_type; kwargs...) + function setup!(solver::Marble.Solver,model::Model, ind_cc1, ind_cc2, cc_type; kwargs...) Marble.update_settings!(solver; kwargs...) - setup!(solver, MathOptNLPModel(model), ind_cc1, ind_cc2, comp_type; kwargs...) + setup!(solver, MathOptNLPModel(model), ind_cc1, ind_cc2, cc_type; kwargs...) return nothing end - function setup!(solver::Marble.Solver,nlp::AbstractNLPModel, ind_cc1, ind_cc2, comp_type; kwargs...) + function setup!(solver::Marble.Solver,nlp::AbstractNLPModel, ind_cc1, ind_cc2, cc_type; kwargs...) Marble.update_settings!(solver; kwargs...) opts = Marble.options(solver) - data = jump_to_marble(nlp, ind_cc1, ind_cc2, comp_type) + data = jump_to_marble(nlp, ind_cc1, ind_cc2, cc_type) Marble.set_problem!(solver, data.Q, data.q, data.c0, data.J_eq, data.b_eq, data.J_ineq, data.b_ineq, data.L, data.l, data.R, data.r, opts) return nothing end diff --git a/julia/src/conversion.jl b/julia/src/conversion.jl index 276edc5..65fc6e7 100644 --- a/julia/src/conversion.jl +++ b/julia/src/conversion.jl @@ -34,7 +34,7 @@ function jump_to_marble( nlp::AbstractNLPModel, ind_cc1, ind_cc2, - comp_type::Union{ + cc_type::Union{ Tuple{Symbol,Symbol}, AbstractVector{<:Tuple{Symbol,Symbol}}, }, @@ -44,20 +44,20 @@ function jump_to_marble( ncc = length(ind_cc1) - comp_type = comp_type isa Tuple ? fill(comp_type, ncc) : collect(comp_type) + cc_type = cc_type isa Tuple ? fill(cc_type, ncc) : collect(cc_type) nvar = nlp.meta.nvar ncon = nlp.meta.ncon @assert nlp.meta.nnln == 0 "Expected no nonlinear constraints in `nlp`" @assert length(ind_cc2) == ncc "Expected equal number of complementarity endpoints" - @assert length(comp_type) == ncc "Expected one complementarity type per pair" - @assert all(t -> all(in((:var, :con)), t), comp_type) """ - Expected comp_type entries of form (:var, :var), (:var, :con), (:con, :var), or (:con, :con) + @assert length(cc_type) == ncc "Expected one complementarity type per pair" + @assert all(t -> all(in((:var, :con)), t), cc_type) """ + Expected cc_type entries of form (:var, :var), (:var, :con), (:con, :var), or (:con, :con) """ - kind1 = first.(comp_type) - kind2 = last.(comp_type) + kind1 = first.(cc_type) + kind2 = last.(cc_type) valid(k, i) = 1 <= i <= (k === :var ? nvar : ncon) row(k, i) = k === :var ? ncon + i : i From 1b057589ce6b1cf7c6a2b37489d0b81160c774f9 Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Sun, 31 May 2026 18:18:28 -0400 Subject: [PATCH 2/3] update python, add julia and python tests --- .gitignore | 5 +- CMakeLists.txt | 30 ++- bindings/python_bindings.cpp | 10 +- julia/Project.toml | 6 + julia/src/Marble.jl | 46 ++++ julia/test.jl | 222 ------------------ julia/test/runtests.jl | 100 ++++++++ julia/tester.jl | 156 ------------ python/examples/macmpec_lcqp.py | 46 ---- python/examples/simple_test.py | 27 +++ python/marble/__init__.py | 107 +++++++++ .../{typings/marble.pyi => marble/_core.pyi} | 12 +- python/pyproject.toml | 7 +- python/tester.py | 39 --- python/tests/test_marble.py | 99 ++++++++ 15 files changed, 433 insertions(+), 479 deletions(-) delete mode 100644 julia/test.jl create mode 100644 julia/test/runtests.jl delete mode 100644 julia/tester.jl delete mode 100644 python/examples/macmpec_lcqp.py create mode 100644 python/examples/simple_test.py create mode 100644 python/marble/__init__.py rename python/{typings/marble.pyi => marble/_core.pyi} (95%) delete mode 100644 python/tester.py create mode 100644 python/tests/test_marble.py diff --git a/.gitignore b/.gitignore index e8adc9a..5b061e7 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ .vscode/ .claude/ build/* -.venv/ .DS_Store -julia/examples/Manifest.toml \ No newline at end of file +julia/examples/Manifest.toml +python/marble/__pycache__/ +python/.venv/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 1f85223..af8acb9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -91,26 +91,39 @@ if(MARBLE_BUILD_PYTHON) FetchContent_MakeAvailable(pybind11) endif() + # The compiled extension is the private submodule `marble._core`; the pure + # Python `marble` package (python/marble/) wraps it with the high-level API. pybind11_add_module(marble_python MODULE bindings/python_bindings.cpp) target_link_libraries(marble_python PRIVATE marble) set_target_properties(marble_python PROPERTIES - OUTPUT_NAME marble - LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/python + OUTPUT_NAME _core + LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/python/marble ) - # .pyi stub generation (for VSCode / Pylance autocomplete) + # Stage the pure Python package next to the freshly built extension so that + # adding build/python to sys.path is enough to `import marble` from the tree. + add_custom_command(TARGET marble_python POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy_directory + ${CMAKE_SOURCE_DIR}/python/marble ${CMAKE_BINARY_DIR}/python/marble + COMMENT "Staging marble package into build/python/marble" + VERBATIM + ) + + # .pyi stub generation for the compiled core (VSCode / Pylance autocomplete) if(MARBLE_GENERATE_PYI) find_program(PYBIND11_STUBGEN pybind11-stubgen) if(PYBIND11_STUBGEN) add_custom_command(TARGET marble_python POST_BUILD - COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_SOURCE_DIR}/python/typings COMMAND ${CMAKE_COMMAND} -E env PYTHONPATH=${CMAKE_BINARY_DIR}/python - ${PYBIND11_STUBGEN} marble - --output-dir ${CMAKE_SOURCE_DIR}/python/typings + ${PYBIND11_STUBGEN} marble._core + --output-dir ${CMAKE_SOURCE_DIR}/python --numpy-array-remove-parameters --ignore-unresolved-names "^(m|n)$" - COMMENT "Generating marble.pyi stubs" + COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_SOURCE_DIR}/python/marble/_core.pyi + ${CMAKE_BINARY_DIR}/python/marble/_core.pyi + COMMENT "Generating marble/_core.pyi stubs" VERBATIM ) else() @@ -120,7 +133,8 @@ if(MARBLE_BUILD_PYTHON) endif() endif() - install(TARGETS marble_python LIBRARY DESTINATION .) + # Wheel layout: site-packages/marble/{__init__.py, _core.so} + install(TARGETS marble_python LIBRARY DESTINATION marble) endif() # --------------------------------------------------------------------------- diff --git a/bindings/python_bindings.cpp b/bindings/python_bindings.cpp index 8d71dd0..89f446e 100644 --- a/bindings/python_bindings.cpp +++ b/bindings/python_bindings.cpp @@ -38,8 +38,8 @@ static Vec arr_to_vec(py::array_t arr) { return Eigen::Map(arr.data(), arr.size()); } -PYBIND11_MODULE(marble, m) { - m.doc() = "Marble: constrained optimization solver with complementarity constraints"; +PYBIND11_MODULE(_core, m) { + m.doc() = "Marble compiled core: constrained optimization solver with complementarity constraints"; // ------------------------------------------------------------------------- // Problem @@ -78,6 +78,10 @@ PYBIND11_MODULE(marble, m) { .def_readonly("n_eq", &Problem::n_eq) .def_readonly("n_ineq", &Problem::n_ineq) .def_readonly("n_comp", &Problem::n_comp) + .def("obj", [](const Problem& p, const Vec& z) { return p.obj(z); }, py::arg("z")) + .def("residual_eq", [](const Problem& p, const Vec& z) { return p.residual_eq(z); }, py::arg("z")) + .def("residual_ineq", [](const Problem& p, const Vec& z) { return p.residual_ineq(z); }, py::arg("z")) + .def("residual_comp", [](const Problem& p, const Vec& z) { return p.residual_comp(z); }, py::arg("z")) .def_property_readonly("cost_gradient", [](const Problem& p) { return p.cost_gradient; }) .def_readonly("cost_const", &Problem::cost_const) // Constraint matrices (returned as scipy-compatible CSC tuple: rows, cols, colptr, rowval, nzval) @@ -326,6 +330,8 @@ PYBIND11_MODULE(marble, m) { py::arg("niter")) .def("get_problem", &Solver::get_problem, py::return_value_policy::reference_internal) + .def("get_options", [](Solver& s) -> Solver::Options& { return s.options; }, + py::return_value_policy::reference_internal) // Main solve .def("solve", [](Solver& s) { return s.solve(); }) .def("convergence", &Solver::convergence, py::arg("options")) diff --git a/julia/Project.toml b/julia/Project.toml index 13e5ab3..af86599 100644 --- a/julia/Project.toml +++ b/julia/Project.toml @@ -25,3 +25,9 @@ NLPModelsJuMP = "0.13.5" Preferences = "1.5.2" SparseArrays = "1.10.0" julia = "1.10" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/julia/src/Marble.jl b/julia/src/Marble.jl index d7523bd..578f9f1 100644 --- a/julia/src/Marble.jl +++ b/julia/src/Marble.jl @@ -29,6 +29,52 @@ module Marble return nothing end + function setup!(solver::Marble.Solver, Q::AbstractMatrix, q::AbstractVector, c0::Real=0.0; + J_eq=nothing, b_eq=nothing, J_ineq=nothing, b_ineq=nothing, + L=nothing, l=nothing, R=nothing, r=nothing, kwargs...) + Marble.update_settings!(solver; kwargs...) + opts = Marble.options(solver) + + n = length(q) + J_eq = isnothing(J_eq) ? zeros(0, n) : J_eq + b_eq = isnothing(b_eq) ? zeros(0) : b_eq + J_ineq = isnothing(J_ineq) ? zeros(0, n) : J_ineq + b_ineq = isnothing(b_ineq) ? zeros(0) : b_ineq + L = isnothing(L) ? zeros(0, n) : L + l = isnothing(l) ? zeros(0) : l + R = isnothing(R) ? zeros(0, n) : R + r = isnothing(r) ? zeros(0) : r + + _set_problem!(solver, opts, Q, q, Float64(c0), + J_eq, b_eq, J_ineq, b_ineq, L, l, R, r) + return nothing + end + + # Dispatch to the dense or sparse set_problem! binding based on storage. When + # any block is sparse every block is converted to a SparseMatrixCSC so the + # solver sees consistent compressed-sparse data + function _set_problem!(solver::Marble.Solver, opts, Q, q, c0, + J_eq, b_eq, J_ineq, b_ineq, L, l, R, r) + blocks = (Q, J_eq, J_ineq, L, R) + fvec(v) = collect(Float64, v) + if any(b -> b isa AbstractSparseMatrix, blocks) + Qs, Es, Is, Ls, Rs = sparse(Q), sparse(J_eq), sparse(J_ineq), sparse(L), sparse(R) + Marble.set_problem!(solver, size(Q, 2), + Qs.colptr, Qs.rowval, Qs.nzval, fvec(q), c0, + length(b_eq), Es.colptr, Es.rowval, Es.nzval, fvec(b_eq), + length(b_ineq), Is.colptr, Is.rowval, Is.nzval, fvec(b_ineq), + length(l), Ls.colptr, Ls.rowval, Ls.nzval, fvec(l), + Rs.colptr, Rs.rowval, Rs.nzval, fvec(r), opts) + else + fmat(M) = Matrix{Float64}(M) + Marble.set_problem!(solver, + fmat(Q), fvec(q), c0, + fmat(J_eq), fvec(b_eq), fmat(J_ineq), fvec(b_ineq), + fmat(L), fvec(l), fmat(R), fvec(r), opts) + end + return nothing + end + const OPTIONS = [ :convergence_kkt_norm, :convergence_eq_violation, :convergence_ineq_violation, :convergence_comp_violation, :outer_step_kkt_norm, diff --git a/julia/test.jl b/julia/test.jl deleted file mode 100644 index 7185074..0000000 --- a/julia/test.jl +++ /dev/null @@ -1,222 +0,0 @@ -using Pkg; -Pkg.activate(joinpath(@__DIR__)) -using CxxWrap -using Revise -using CILogDomain -using SparseArrays -using LinearAlgebra - -module Marble - using CxxWrap - @wrapmodule(() -> joinpath(@__DIR__, "submodules/Marble/build/libmarble_wrapper")) - - function __init__() - @initcxx - end -end - -name = :hopper -## - -problem = create_problem(name)#, eq = 1, ineq = 1, comp = 1); -solver_options = SolverOptions() -solver_options.verbose = true -solver_options.max_iters = 10_000 -z_ours, solver = solve(problem; options=solver_options); - -# Check eigen to julia and back -println(Marble.test_array([1.0, 2.0, 3.0])) - -arr = [1.0 2.0; 3.0 4.0] -println(Marble.test_matrix(arr, 2, 2)) -println(arr) - -# Set up problem and check sizes -H, g, f0 = Matrix(solver.cost_hess), solver.cost_grad, solver.cost_const -J_eq, J_ineq, J_comp = Matrix(solver.conjac[solver.eq_inds, :]), Matrix(solver.conjac[solver.ineq_inds, :]), Matrix(solver.conjac[solver.comp_inds, :]) -c_eq, c_ineq, c_comp = solver.conrhs[solver.eq_inds], solver.conrhs[solver.ineq_inds], solver.conrhs[solver.comp_inds] - -prob = Marble.Problem(H, g, f0, J_eq, c_eq, J_ineq, c_ineq, - J_comp, c_comp) -@assert Marble.nz(prob) == solver.nz -@assert Marble.n_eq(prob) == solver.n_eq -@assert Marble.n_ineq(prob) == solver.n_ineq -@assert Marble.n_comp(prob) == solver.n_comp - -# Set up solver, test setting and getting problem -marble_solver = Marble.Solver() - -# Check retraction and retraction derivative -s = randn(10) -@assert norm(Marble.retract(marble_solver, s, 1e-2) - solver.r(s, 1e-2)) < 1e-10 -@assert norm(Marble.retract_deriv(marble_solver, s, 1e-2) - solver.d_r(s, 1e-2)) < 1e-10 -@assert norm(Marble.retract_second_deriv(marble_solver, s, 1e-2) - solver.dd_r(s, 1e-2)) < 1e-10 - -# Set problem -Marble.set_problem(marble_solver, prob) -prob = Marble.get_problem(marble_solver) -@assert Marble.nz(prob) == solver.nz -@assert Marble.n_eq(prob) == solver.n_eq -@assert Marble.n_ineq(prob) == solver.n_ineq -@assert Marble.n_comp(prob) == solver.n_comp - -# Get workspace -workspace = Marble.get_workspace(marble_solver) - -# Ruiz scaling check -@assert norm(Marble.scaling(workspace) - diag(solver.ruiz_mat), Inf) < 1e-10 - -# Test indices -Marble.z_inds(marble_solver) == solver.kkt_inds.z .- 1 -Marble.s_ineq_inds(marble_solver) == solver.kkt_inds.v .- 1 -Marble.s_comp_inds(marble_solver) == solver.kkt_inds.σ .- 1 -Marble.m_eq_inds(marble_solver) == solver.kkt_inds.λ .- 1 -Marble.m_ineq_inds(marble_solver) == solver.kkt_inds.μ .- 1 -Marble.m_comp_inds(marble_solver) == solver.kkt_inds.τ .- 1 - -# Check KKT structure -function get_kkt(marble_solver) - nr, nc, colptr, rowval, nzval = Marble.kkt_system(Marble.get_workspace(marble_solver)) - return copy(SparseMatrixCSC(nr, nc, colptr.+1, rowval.+1, nzval)) -end -kkt_sparse = get_kkt(marble_solver) -perm = Marble.amd_perm_vec(workspace) .+ 1 -iperm = Marble.amd_iperm_vec(workspace) .+ 1 - -# Check indices into nzval (can't compare against solver one because we only use upper triangular form) -# Also the matrix is permuted, so we have to check against permuted indices -kkt_sparse.nzval[:] = 1:length(kkt_sparse.nzval) -kkt_full = (kkt_sparse + kkt_sparse' - spdiagm(diag(kkt_sparse)))[iperm, iperm] # Unpermuted version - -norm(diagm(kkt_sparse.nzval[Marble.z_z_inds(marble_solver) .+ 1]) - kkt_full[solver.kkt_inds.z, solver.kkt_inds.z], Inf) -norm(diagm(kkt_sparse.nzval[Marble.s_ineq_s_ineq_inds(marble_solver) .+ 1]) - kkt_full[solver.kkt_inds.v, solver.kkt_inds.v], Inf) -norm(diagm(kkt_sparse.nzval[Marble.s_ineq_m_ineq_inds(marble_solver) .+ 1]) - kkt_full[solver.kkt_inds.v, solver.kkt_inds.μ], Inf) -norm(diagm(kkt_sparse.nzval[Marble.s_comp_s_comp_inds(marble_solver) .+ 1]) - kkt_full[solver.kkt_inds.σ, solver.kkt_inds.σ], Inf) -norm(kkt_sparse.nzval[Marble.s_comp_m_comp_inds(marble_solver) .+ 1] - kkt_full[solver.kkt_inds.σ, solver.kkt_inds.τ].nzval, Inf) - -# Update ineq, comp and penalty terms -iter = solver.iters[2] -Marble.update_KKT_ineq(marble_solver, iter.v, sqrt(iter.κ)) -Marble.update_KKT_comp(marble_solver, iter.σ, iter.τ, sqrt(iter.κ)) -Marble.update_KKT_penalty(marble_solver, 1/iter.ρ) -kkt = get_kkt(marble_solver)[iperm, iperm] # Unpermuted -kkt = (kkt + kkt' - spdiagm(diag(kkt))) - -# Check stationarity rows -scaling = solver.ruiz_mat -hess = scaling*lagrangian_hessian(solver, iter)*scaling -@assert norm(hess[solver.kkt_inds.z, :] - kkt[solver.kkt_inds.z, :], Inf) < 1e-10 -@assert norm(hess[solver.kkt_inds.v, :] - kkt[solver.kkt_inds.v, :], Inf) < 1e-10 -@assert norm(hess[solver.kkt_inds.σ, :] - kkt[solver.kkt_inds.σ, :], Inf) < 1e-10 - -# Check entire matrix -@assert norm(hess - kkt, Inf) < 1e-10 - -# Check across all iterates -for iter in solver.iters - Marble.update_KKT_ineq(marble_solver, iter.v, sqrt(iter.κ)) - Marble.update_KKT_comp(marble_solver, iter.σ, iter.τ, sqrt(iter.κ)) - Marble.update_KKT_penalty(marble_solver, 1/iter.ρ) - hess = (scaling*lagrangian_hessian(solver, iter)*scaling)[perm, perm] - kkt = get_kkt(marble_solver) - @assert norm(triu(hess) - kkt, Inf) < 1e-10 -end - -# Check writing to workspace -function set_from_iter(marble_solver, iter) - workspace = Marble.get_workspace(marble_solver) - Marble.z(workspace) .= iter.z - Marble.s_ineq(workspace) .= iter.v - Marble.s_comp(workspace) .= iter.σ - Marble.m_eq(workspace) .= iter.λ - Marble.m_ineq(workspace) .= iter.μ - Marble.m_comp(workspace) .= iter.τ - Marble.m_eq_est(workspace) .= iter.α - Marble.m_ineq_est(workspace) .= iter.β - Marble.m_comp_est(workspace) .= iter.γ - Marble.update_KKT_residual(marble_solver, sqrt(iter.κ), 1/iter.ρ) - Marble.update_KKT_system(marble_solver, sqrt(iter.κ), 1/iter.ρ) - Marble.update_KKT_primal_regularizer(marble_solver, iter.reg) -end -set_from_iter(marble_solver, solver.iters[2]) -residual = Marble.kkt_residual(workspace) -kkt = get_kkt(marble_solver) -@assert norm(lagrangian_gradient(solver, iter) - residual, Inf) < 1e-10 -@assert norm(triu((scaling*(lagrangian_hessian(solver, iter) + iter.reg*solver.reg_mat)*scaling)[perm, perm]) - kkt, Inf) < 1e-10 - -# Check across all iteraters -for iter in solver.iters - # Write current solution guess - set_from_iter(marble_solver, iter) - - # Check against iterate - residual = Marble.kkt_residual(workspace) - kkt = get_kkt(marble_solver) - @assert norm(lagrangian_gradient(solver, iter) - residual, Inf) < 1e-10 - @assert norm(triu((scaling*(lagrangian_hessian(solver, iter) + iter.reg*solver.reg_mat)*scaling)[perm, perm]) - kkt, Inf) < 1e-10 -end - -# Check regularizer updating -Marble.update_KKT_primal_regularizer(marble_solver, 0.0) -kkt = get_kkt(marble_solver) -Marble.update_KKT_primal_regularizer(marble_solver, 1e-4) -kkt_reg = get_kkt(marble_solver) -@assert norm((kkt_reg - kkt) - 1e-4*(scaling*solver.reg_mat*scaling)[perm, perm], Inf) < 1e-12 -Marble.update_KKT_primal_regularizer(marble_solver, 1.23e-7) -kkt_reg = copy(get_kkt(marble_solver)) -@assert norm((kkt_reg - kkt) - 1.23e-7*(scaling*solver.reg_mat*scaling)[perm, perm], Inf) < 1e-12 -Marble.update_KKT_primal_regularizer(marble_solver, 0.0) -kkt_reg = copy(get_kkt(marble_solver)) -@assert norm(kkt_reg - kkt, Inf) < 1e-14 - -# Test factorization and solve -iter = solver.iters[5] -set_from_iter(marble_solver, iter) -res = copy(Marble.kkt_residual(Marble.get_workspace(marble_solver))) -reg = 1e-7 -Marble.update_KKT_primal_regularizer(marble_solver, reg) -kkt = get_kkt(marble_solver)[iperm, iperm] # Unpermuted -kkt = (kkt + kkt' - spdiagm(diag(kkt))) - -@assert Marble.analytical_factorization(marble_solver) -@assert Marble.numerical_factorization(marble_solver) -Marble.backsolve(marble_solver) -qdldl_soln = Marble.newton_step(Marble.get_workspace(marble_solver)) -@assert(norm(kkt*(scaling \ qdldl_soln) + (scaling*res), Inf) < 1e-9) - -# Test for each iter -using QDLDL -for iter in solver.iters - if iter.type != :Newton - continue - end - set_from_iter(marble_solver, iter) - - # Get C++ terms - res = copy(Marble.kkt_residual(Marble.get_workspace(marble_solver))) - kkt = get_kkt(marble_solver)[iperm, iperm] # Unpermuted - kkt = (kkt + kkt' - spdiagm(diag(kkt))) - - # Solve C++ system - @assert Marble.analytical_factorization(marble_solver) - @assert Marble.numerical_factorization(marble_solver) - @assert Marble.check_inertia(marble_solver) - Marble.backsolve(marble_solver) - qdldl_soln = Marble.newton_step(Marble.get_workspace(marble_solver)) - - # Check against Julia - hess = lagrangian_hessian(solver, iter) + iter.reg*solver.reg_mat - grad = lagrangian_gradient(solver, iter) - F = QDLDL.qdldl(hess) - qdldl_soln_julia = QDLDL.solve(F, -grad) - @assert norm(hess*qdldl_soln + grad, Inf) < 1e-5 - - # Check iterate passes filter - marble_filter_viol, marble_filter_obj = Marble.entry_from_solution(marble_solver, sqrt(iter.κ), 1/iter.ρ) - iter_filter_obj, iter_filter_viol = filter_criteria(solver, iter) - - @assert abs(marble_filter_viol - iter_filter_viol) < 1e-10 - @assert abs(marble_filter_obj - iter_filter_obj) < 1e-10 -end - -println("All tests passed!") diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl new file mode 100644 index 0000000..31e5208 --- /dev/null +++ b/julia/test/runtests.jl @@ -0,0 +1,100 @@ +# Marble test suite +# +# Mirror of the Python suite (python/tests/test_marble.py): the same example +# problem written directly as matrices (from julia/examples/simple_test.jl) and +# the same eight test cases in the same order. + +using Test +using LinearAlgebra +using SparseArrays +using Marble + +# Example problem (julia/examples/simple_test.jl), written directly as matrices: +# min x'x +# s.t. x1 = 1, x2 >= 1, 0 <= (x3 + 1) ⟂ (x4 - 1) >= 0 +# optimum x* = [1, 1, 0, 1] +const Q = 2.0 * Matrix(1.0I, 4, 4) +const q = zeros(4) +const C0 = 0.0 +const J_EQ, B_EQ = [1.0 0.0 0.0 0.0], [-1.0] +const J_INEQ, B_INEQ = [0.0 1.0 0.0 0.0], [-1.0] +const L, EL = [0.0 0.0 1.0 0.0], [1.0] +const R, ER = [0.0 0.0 0.0 1.0], [-1.0] +const ZSTAR = [1.0, 1.0, 0.0, 1.0] + +# Build the solver, set up the example problem (dense or sparse), and solve. +# Option settings are forwarded to setup! as keyword arguments. +function setup_and_solve(; sparse_problem = false, opts...) + conv = sparse_problem ? sparse : identity + solver = Marble.Solver() + Marble.setup!(solver, conv(Q), q, C0; + J_eq = conv(J_EQ), b_eq = B_EQ, J_ineq = conv(J_INEQ), b_ineq = B_INEQ, + L = conv(L), l = EL, R = conv(R), r = ER, opts...) + return solver, Marble.solve!(solver) +end + +@testset "Marble" begin + + @testset "problem dimensions" begin + solver, _ = setup_and_solve() + p = solver.problem + @test Marble.nz(p) == 4 && Marble.n_eq(p) == 1 && Marble.n_ineq(p) == 1 && Marble.n_comp(p) == 1 + end + + @testset "solves to known optimum" begin + _, res = setup_and_solve() + z = collect(Marble.z(res)) + @test Marble.converged(res) && isapprox(z, ZSTAR, atol = 1e-4) + end + + @testset "equality constraint satisfied" begin + # x1 = 1 + _, res = setup_and_solve() + @test abs(collect(Marble.z(res))[1] - 1.0) < 1e-4 + end + + @testset "inequality constraint satisfied" begin + # x2 >= 1 + _, res = setup_and_solve() + @test collect(Marble.z(res))[2] >= 1.0 - 1e-4 + end + + @testset "complementarity satisfied" begin + # 0 <= (x3 + 1) ⟂ (x4 - 1) >= 0 + _, res = setup_and_solve() + z = collect(Marble.z(res)) + a = z[3] + 1.0 + b = z[4] - 1.0 + @test a >= -1e-4 && b >= -1e-4 && abs(a * b) < 1e-4 + end + + @testset "objective value" begin + # x'x = 3 at the optimum + solver, res = setup_and_solve() + @test abs(Marble.obj(solver, collect(Marble.z(res))) - 3.0) < 1e-3 + end + + @testset "dense and sparse solve agree" begin + _, rd = setup_and_solve(sparse_problem = false) + _, rs = setup_and_solve(sparse_problem = true) + @test Marble.converged(rd) && Marble.converged(rs) && + isapprox(collect(Marble.z(rd)), collect(Marble.z(rs)), atol = 1e-6) + end + + @testset "complementarity-only QPCC" begin + # same cost x'x, only the complementarity 0 <= (x3 + 1) ⟂ (x4 - 1) >= 0 + # x1, x2 are unconstrained -> 0; x3 = 0, x4 = 1 => x* = [0, 0, 0, 1] + solver = Marble.Solver() + Marble.setup!(solver, Q, q, C0; L = L, l = EL, R = R, r = ER) + res = Marble.solve!(solver) + z = collect(Marble.z(res)) + @test Marble.converged(res) && isapprox(z, [0.0, 0.0, 0.0, 1.0], atol = 1e-4) + end + + @testset "options are honored" begin + s2, r2 = setup_and_solve(max_iters = 123, convergence_kkt_norm = 1e-8) + @test Marble.max_iters(s2.options) == 123 && + isapprox(Marble.convergence_kkt_norm(s2.options), 1e-8) && + Marble.converged(r2) + end +end diff --git a/julia/tester.jl b/julia/tester.jl deleted file mode 100644 index ee18efe..0000000 --- a/julia/tester.jl +++ /dev/null @@ -1,156 +0,0 @@ -using CxxWrap -using LinearAlgebra -using SparseArrays - -# --------------------------------------------------------------------------- -# Load the module -# --------------------------------------------------------------------------- -module Marble - using CxxWrap - @wrapmodule(() -> joinpath(@__DIR__, "../build/libmarble_julia")) - function __init__() - @initcxx - end -end - -# --------------------------------------------------------------------------- -# Helper: reconstruct a Julia SparseMatrixCSC from the kkt_system tuple -# --------------------------------------------------------------------------- -function get_kkt(solver) - ws = Marble.get_workspace(solver) - nr, nc, colptr, rowval, nzval = Marble.kkt_system(ws) - return SparseMatrixCSC(nr, nc, colptr .+ 1, rowval .+ 1, copy(nzval)) -end - -# --------------------------------------------------------------------------- -# Test 1 – Problem construction and dimension accessors -# --------------------------------------------------------------------------- -println("Test 1: Problem construction") - -nz = 2 -H = Matrix(1.0 * I(nz)) # Identity Hessian -g = [1.0, 2.0] # Linear cost -J_empty = zeros(0, nz) -c_empty = zeros(0) - -prob = Marble.Problem(H, g, 0.0, J_empty, c_empty, J_empty, c_empty, J_empty, c_empty) -@assert Marble.nz(prob) == 2 -@assert Marble.n_eq(prob) == 0 -@assert Marble.n_ineq(prob) == 0 -@assert Marble.n_comp(prob) == 0 -println(" PASSED") - -# --------------------------------------------------------------------------- -# Test 2 – SolverOptions getters and setters -# --------------------------------------------------------------------------- -println("Test 2: SolverOptions") - -opts = Marble.SolverOptions() -Marble.max_iters!(opts, 500) -Marble.convergence_kkt_norm!(opts, 1e-6) -Marble.output_dir!(opts, "/dev/null") -@assert Marble.max_iters(opts) == 500 -@assert Marble.convergence_kkt_norm(opts) ≈ 1e-6 -@assert Marble.output_dir(opts) == "/dev/null" -println(" PASSED") - -# --------------------------------------------------------------------------- -# Test 3 – FilterEntry construction and accessors -# --------------------------------------------------------------------------- -println("Test 3: FilterEntry") - -entry = Marble.FilterEntry(0.5, 1.2) -@assert Marble.feas(entry) ≈ 0.5 -@assert Marble.merit(entry) ≈ 1.2 -Marble.feas!(entry, 0.1) -@assert Marble.feas(entry) ≈ 0.1 -println(" PASSED") - -# --------------------------------------------------------------------------- -# Test 4 – Filter update and acceptability -# --------------------------------------------------------------------------- -println("Test 4: Filter") - -filt = Marble.Filter() -@assert Marble.num_entries(filt) == 0 -Marble.update(filt, Marble.FilterEntry(1.0, 1.0)) -@assert Marble.num_entries(filt) == 1 -# entries() returns a flat [feas, merit, ...] vector -flat = Marble.entries(filt) -@assert flat[1] ≈ 1.0 && flat[2] ≈ 1.0 -# A candidate strictly better on both axes should be acceptable -@assert Marble.acceptable(filt, Marble.FilterEntry(0.5, 0.5)) -Marble.clear(filt) -@assert Marble.num_entries(filt) == 0 -println(" PASSED") - -# --------------------------------------------------------------------------- -# Test 5 – Retraction maps -# --------------------------------------------------------------------------- -println("Test 5: Retraction maps") - -solver = Marble.Solver(opts) -Marble.set_problem(solver, prob, ones(nz)) - -s = [0.1, 0.2] -r = Marble.retract(solver, s, 0.1) -dr = Marble.retract_deriv(solver, s, 0.1) -ddr = Marble.retract_second_deriv(solver, s, 0.1) -@assert length(r) == nz -@assert length(dr) == nz -@assert length(ddr) == nz -println(" PASSED") - -# --------------------------------------------------------------------------- -# Test 6 – Solve and check solution -# Unconstrained QP: min 0.5 * x' * I * x + [1,2]' * x -# Optimal solution: x* = [-1, -2] -# --------------------------------------------------------------------------- -println("Test 6: Solve unconstrained QP") - -converged = Marble.solve(solver, opts) -@assert converged - -ws = Marble.get_workspace(solver) -z = copy(Marble.z(ws)) -println(" Solution z = $z") -println(" Expected = [-1.0, -2.0]") -@assert norm(z - [-1.0, -2.0], Inf) < 1e-3 -println(" PASSED") - -# --------------------------------------------------------------------------- -# Test 7 – Workspace residuals -# --------------------------------------------------------------------------- -println("Test 7: Workspace residuals") - -@assert norm(Marble.kkt_residual(ws), Inf) < 1e-5 -@assert norm(Marble.residual_eq(ws), Inf) < 1e-8 -@assert norm(Marble.residual_ineq(ws), Inf) < 1e-8 -@assert norm(Marble.residual_comp(ws), Inf) < 1e-8 -println(" PASSED") - -# --------------------------------------------------------------------------- -# Test 8 – KKT index vectors are non-empty and correctly sized -# --------------------------------------------------------------------------- -println("Test 8: KKT index vectors") - -@assert length(Marble.z_inds(solver)) == nz -@assert length(Marble.m_eq_inds(solver)) == 0 # no equality constraints -println(" n_primals = $(Marble.n_primals(solver))") -println(" n_duals = $(Marble.n_duals(solver))") -println(" n_vars = $(Marble.n_vars(solver))") -println(" PASSED") - -# --------------------------------------------------------------------------- -# Test 9 – KKT system sparse matrix reconstruction -# --------------------------------------------------------------------------- -println("Test 9: KKT system sparse matrix") - -kkt = get_kkt(solver) -@assert size(kkt, 1) == size(kkt, 2) -@assert size(kkt, 1) == Marble.n_vars(solver) -println(" KKT size: $(size(kkt))") -println(" PASSED") - -# --------------------------------------------------------------------------- -println("\nAll tests passed!") diff --git a/python/examples/macmpec_lcqp.py b/python/examples/macmpec_lcqp.py deleted file mode 100644 index 166cc85..0000000 --- a/python/examples/macmpec_lcqp.py +++ /dev/null @@ -1,46 +0,0 @@ -import h5py -import numpy as np -import marble -import json -import logging -import os - - -def load_problem(path): - with h5py.File(path, "r") as f: - return {k: np.asarray(f[k]).T for k in f.keys()} - - -def solve_macmpec_lcqp(name): - data = load_problem(os.path.join(os.path.dirname(__file__), f"../../examples/data/macmpec/h5/{name}.h5")) - logger.info(f"{data['n_eq']} equality constraints, {data['n_ineq']} inequality constraints, {data['n_comp']} complementarity constraints") - - H, q, f0 = data["H"], data["q"], float(data["f0"].item()) - J_eq, b_eq = data["J_eq"], data["b_eq"] - J_ineq, b_ineq = data["J_ineq"], data["b_ineq"] - J_comp, b_comp = data["J_comp"], data["b_comp"] - - obj = lambda z: 0.5 * z.T @ H @ z + q.T @ z + f0 - - solver = marble.Solver() - problem = marble.Problem(H, q, f0, - J_eq, b_eq, - J_ineq, b_ineq, - J_comp, b_comp) - solver.set_problem(problem) - assert solver.solve(marble.SolverOptions()), "Failed to solve problem" - - logger.info(f"\tObjective value: {obj(solver.get_workspace().z)}") - - -if __name__ == "__main__": - logging.basicConfig(level=logging.INFO) - logger = logging.getLogger(__name__) - - json_path = os.path.join(os.path.dirname(__file__), "../../examples/data/macmpec/macmpec_lcqp.json") - with open(json_path, "r") as f: - json_data = json.load(f) - - for key in sorted(json_data.keys()): - logger.info(f"Solving MacMPEC problem: {key}") - solve_macmpec_lcqp(key) diff --git a/python/examples/simple_test.py b/python/examples/simple_test.py new file mode 100644 index 0000000..5f737c6 --- /dev/null +++ b/python/examples/simple_test.py @@ -0,0 +1,27 @@ +# Construct a simple test problem minimizing x'x +# where x[1] = 1 +# x[2] >= 1 +# 0 <= (x[3] + 1) ⟂ (x[4] - 1) >= 0 --> solution is x[3] = 0, x[4] = 1 +import numpy as np +import marble + +Q = 2.0 * np.eye(4) +q = np.zeros(4) + +J_eq, b_eq = [[1.0, 0.0, 0.0, 0.0]], [-1.0] # x[1] == 1 +J_ineq, b_ineq = [[0.0, 1.0, 0.0, 0.0]], [-1.0] # x[2] >= 1 +L, l = [[0.0, 0.0, 1.0, 0.0]], [1.0] # x[3] + 1 +R, r = [[0.0, 0.0, 0.0, 1.0]], [-1.0] # x[4] - 1 + +solver = marble.Solver() +solver.setup(Q, q, 0.0, + J_eq=J_eq, b_eq=b_eq, J_ineq=J_ineq, b_ineq=b_ineq, + L=L, l=l, R=R, r=r, verbosity=1) +results = solver.solve() + +z = np.asarray(results.z) +print(z) +print(solver.problem.obj(z)) +print(solver.problem.residual_eq(z)) +print(solver.problem.residual_ineq(z)) +print(solver.problem.residual_comp(z)) diff --git a/python/marble/__init__.py b/python/marble/__init__.py new file mode 100644 index 0000000..3ea14d8 --- /dev/null +++ b/python/marble/__init__.py @@ -0,0 +1,107 @@ +from __future__ import annotations + +import warnings + +import numpy as np + +from . import _core +from ._core import Filter, FilterEntry, Problem, SolveResult, SolverOptions, Workspace + + +def _scipy_sparse(): + """Return the ``scipy.sparse`` module, or ``None`` if SciPy is unavailable.""" + try: + import scipy.sparse as sp + except ImportError: + return None + return sp + + +def _is_sparse(M) -> bool: + sp = _scipy_sparse() + return sp is not None and sp.issparse(M) + + +def _vec(v) -> np.ndarray: + """Coerce to a contiguous 1-D float64 array.""" + return np.ascontiguousarray(np.asarray(v, dtype=np.float64)).reshape(-1) + + +def _dense(A, ncols: int) -> np.ndarray: + """Coerce a matrix to a dense 2-D float64 array, promoting a 1-D row.""" + if _is_sparse(A): + return np.ascontiguousarray(A.toarray(), dtype=np.float64) + a = np.asarray(A, dtype=np.float64) + if a.ndim == 1: + a = a.reshape(0, ncols) if a.size == 0 else a.reshape(1, -1) + return np.ascontiguousarray(a) + + +def _csc(A, ncols: int): + """Coerce any matrix to a scipy CSC float64 matrix.""" + sp = _scipy_sparse() + if _is_sparse(A): + return A.tocsc().astype(np.float64, copy=False) + return sp.csc_matrix(_dense(A, ncols)) + + +class Solver(_core.Solver): + """Solver for ``min 1/2 x'Q x + q'x + c0`` subject to + + ``J_eq x + b_eq == 0``, ``J_ineq x + b_ineq >= 0`` and the complementarity + ``0 <= L x + l ⟂ R x + r >= 0``. + """ + + def setup(self, Q, q, c0=0.0, *, J_eq=None, b_eq=None, J_ineq=None, b_ineq=None, + L=None, l=None, R=None, r=None, **options): + """Set up the problem from matrices and vectors. + + Only ``Q`` and ``q`` are required; any omitted constraint block defaults + to an empty block of the right shape. Solver options (``max_iters``, + ``convergence_kkt_norm``, ...) are passed as keyword arguments; unknown + names are ignored with a warning. + """ + opts = SolverOptions() + for name, value in options.items(): + if hasattr(opts, name): + setattr(opts, name, value) + else: + warnings.warn(f"Unknown solver option {name!r}; ignoring") + + n = _vec(q).shape[0] + J_eq = np.zeros((0, n)) if J_eq is None else J_eq + b_eq = np.zeros(0) if b_eq is None else b_eq + J_ineq = np.zeros((0, n)) if J_ineq is None else J_ineq + b_ineq = np.zeros(0) if b_ineq is None else b_ineq + L = np.zeros((0, n)) if L is None else L + l = np.zeros(0) if l is None else l + R = np.zeros((0, n)) if R is None else R + r = np.zeros(0) if r is None else r + + # If any block is sparse build the whole problem from CSC, else dense + blocks = (Q, J_eq, J_ineq, L, R) + if any(_is_sparse(b) for b in blocks): + if _scipy_sparse() is None: + raise ImportError("SciPy is required to set up a sparse problem") + conv = _csc + else: + conv = _dense + self.set_problem( + conv(Q, n), _vec(q), float(c0), + conv(J_eq, n), _vec(b_eq), + conv(J_ineq, n), _vec(b_ineq), + conv(L, n), _vec(l), + conv(R, n), _vec(r), + opts, + ) + return self + + @property + def options(self) -> SolverOptions: + """The solver's current options.""" + return self.get_options() + + @property + def problem(self) -> Problem: + """The problem currently set on the solver.""" + return self.get_problem() diff --git a/python/typings/marble.pyi b/python/marble/_core.pyi similarity index 95% rename from python/typings/marble.pyi rename to python/marble/_core.pyi index b997393..f6578e4 100644 --- a/python/typings/marble.pyi +++ b/python/marble/_core.pyi @@ -1,5 +1,5 @@ """ -Marble: constrained optimization solver with complementarity constraints +Marble compiled core: constrained optimization solver with complementarity constraints """ from __future__ import annotations import numpy @@ -45,6 +45,14 @@ class Problem: @typing.overload def __init__(self, cost_hessian: scipy.sparse.csc_matrix, cost_gradient: numpy.ndarray, cost_const: float, J_eq: scipy.sparse.csc_matrix, c_eq: numpy.ndarray, J_ineq: scipy.sparse.csc_matrix, c_ineq: numpy.ndarray, L: scipy.sparse.csc_matrix, l: numpy.ndarray, R: scipy.sparse.csc_matrix, r: numpy.ndarray) -> None: ... + def obj(self, z: numpy.ndarray) -> float: + ... + def residual_comp(self, z: numpy.ndarray) -> numpy.ndarray: + ... + def residual_eq(self, z: numpy.ndarray) -> numpy.ndarray: + ... + def residual_ineq(self, z: numpy.ndarray) -> numpy.ndarray: + ... @property def J_comp(self) -> tuple: ... @@ -145,6 +153,8 @@ class Solver: ... def get_filter(self) -> Filter: ... + def get_options(self) -> SolverOptions: + ... def get_problem(self) -> Problem: ... def get_workspace(self) -> Workspace: diff --git a/python/pyproject.toml b/python/pyproject.toml index 8aa16c4..ac93ef0 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -5,9 +5,9 @@ build-backend = "scikit_build_core.build" [project] name = "marble" version = "0.1.0" -description = "Marble: constrained optimization solver with complementarity constraints" +description = "Marble: constrained optimization solver for complementarity constraints" requires-python = ">=3.8" -dependencies = ["numpy"] +dependencies = ["numpy", "scipy"] [tool.scikit-build] cmake.source-dir = ".." @@ -17,4 +17,5 @@ cmake.args = [ "-DMARBLE_BUILD_JULIA=OFF", "-DMARBLE_GENERATE_PYI=OFF", ] -wheel.packages = [] +# Ship the pure Python package; the compiled _core is installed into it by CMake +wheel.packages = ["marble"] diff --git a/python/tester.py b/python/tester.py deleted file mode 100644 index c8782b9..0000000 --- a/python/tester.py +++ /dev/null @@ -1,39 +0,0 @@ -import marble -import numpy as np - -s = marble.Solver() - -# Simple 2-variable unconstrained QP: -# min 0.5 * x' * I * x + [1, 2]' * x -# Optimal solution: x* = [-1, -2] - -nz = 2 -H = np.eye(nz) -g = np.array([1.0, 2.0]) -J_empty = np.zeros((0, nz)) -c_empty = np.zeros(0) - -prob = marble.Problem( - H, g, 0.0, - J_empty, c_empty, # equality - J_empty, c_empty, # inequality - J_empty, c_empty, # complementarity -) - -print(f"Problem: nz={prob.nz}, n_eq={prob.n_eq}, n_ineq={prob.n_ineq}, n_comp={prob.n_comp}") - -opts = marble.SolverOptions() -opts.max_iters = 500 -opts.convergence_kkt_norm = 1e-6 -print(f"Options: {opts}") - -solver = marble.Solver(opts) -solver.set_problem(prob, np.ones(nz)) - -converged = solver.solve(opts) -print(f"Converged: {converged}") - -ws = solver.get_workspace() -print(f"Solution z: {ws.z}") -print(f"Expected: [-1. -2.]") -print(f"KKT residual norm: {np.linalg.norm(ws.kkt_residual):.2e}") diff --git a/python/tests/test_marble.py b/python/tests/test_marble.py new file mode 100644 index 0000000..c2a7c3e --- /dev/null +++ b/python/tests/test_marble.py @@ -0,0 +1,99 @@ +"""Marble Python test suite. + +Mirror of the Julia suite (julia/test/runtests.jl): the same example problem +written directly as matrices (from julia/examples/simple_test.jl) and the same +eight test cases in the same order. Each test sets up and runs the solver fresh. +Run with ``python -m unittest`` or ``python tests/test_marble.py``. +""" +import unittest + +import numpy as np +import scipy.sparse as sp + +import marble + + +# Example problem (julia/examples/simple_test.jl), written directly as matrices: +# min x'x +# s.t. x1 = 1, x2 >= 1, 0 <= (x3 + 1) ⟂ (x4 - 1) >= 0 +# optimum x* = [1, 1, 0, 1] +Q = 2.0 * np.eye(4) +q = np.zeros(4) +C0 = 0.0 +J_EQ, B_EQ = [[1.0, 0.0, 0.0, 0.0]], [-1.0] +J_INEQ, B_INEQ = [[0.0, 1.0, 0.0, 0.0]], [-1.0] +L, EL = [[0.0, 0.0, 1.0, 0.0]], [1.0] +R, ER = [[0.0, 0.0, 0.0, 1.0]], [-1.0] +ZSTAR = np.array([1.0, 1.0, 0.0, 1.0]) + + +# Build the solver, set up the example problem (dense or sparse), and solve. +# Option settings are forwarded to setup as keyword arguments. +def setup_and_solve(sparse_problem=False, **options): + conv = (lambda M: sp.csc_matrix(M)) if sparse_problem else (lambda M: M) + m = marble.Solver() + m.setup(conv(Q), q, C0, + J_eq=conv(J_EQ), b_eq=B_EQ, J_ineq=conv(J_INEQ), b_ineq=B_INEQ, + L=conv(L), l=EL, R=conv(R), r=ER, **options) + return m, m.solve() + + +class TestMarble(unittest.TestCase): + + def test_problem_dimensions(self): + m, _ = setup_and_solve() + p = m.problem + self.assertTrue(p.nz == 4 and p.n_eq == 1 and p.n_ineq == 1 and p.n_comp == 1) + + def test_solves_to_known_optimum(self): + _, res = setup_and_solve() + z = np.asarray(res.z) + self.assertTrue(res.converged and np.allclose(z, ZSTAR, atol=1e-4)) + + def test_equality_constraint_satisfied(self): + # x1 = 1 + _, res = setup_and_solve() + self.assertTrue(abs(np.asarray(res.z)[0] - 1.0) < 1e-4) + + def test_inequality_constraint_satisfied(self): + # x2 >= 1 + _, res = setup_and_solve() + self.assertTrue(np.asarray(res.z)[1] >= 1.0 - 1e-4) + + def test_complementarity_satisfied(self): + # 0 <= (x3 + 1) ⟂ (x4 - 1) >= 0 + _, res = setup_and_solve() + z = np.asarray(res.z) + a = z[2] + 1.0 + b = z[3] - 1.0 + self.assertTrue(a >= -1e-4 and b >= -1e-4 and abs(a * b) < 1e-4) + + def test_objective_value(self): + # x'x = 3 at the optimum + m, res = setup_and_solve() + self.assertTrue(abs(m.problem.obj(np.asarray(res.z)) - 3.0) < 1e-3) + + def test_dense_and_sparse_solve_agree(self): + _, rd = setup_and_solve(sparse_problem=False) + _, rs = setup_and_solve(sparse_problem=True) + self.assertTrue(rd.converged and rs.converged + and np.allclose(np.asarray(rd.z), np.asarray(rs.z), atol=1e-6)) + + def test_complementarity_only_qpcc(self): + # same cost x'x, only the complementarity 0 <= (x3 + 1) ⟂ (x4 - 1) >= 0 + # x1, x2 are unconstrained -> 0; x3 = 0, x4 = 1 => x* = [0, 0, 0, 1] + m = marble.Solver() + m.setup(Q, q, C0, L=L, l=EL, R=R, r=ER) + res = m.solve() + z = np.asarray(res.z) + self.assertTrue(res.converged and np.allclose(z, [0.0, 0.0, 0.0, 1.0], atol=1e-4)) + + def test_options_are_honored(self): + m, res = setup_and_solve(max_iters=123, convergence_kkt_norm=1e-8) + self.assertTrue(m.options.max_iters == 123 + and m.options.convergence_kkt_norm == 1e-8 + and res.converged) + + +if __name__ == "__main__": + unittest.main() From f6a85002df8551787423a8bc646a30abcc492899 Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Sun, 31 May 2026 18:23:53 -0400 Subject: [PATCH 3/3] update readme --- README.md | 116 ++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 91 insertions(+), 25 deletions(-) diff --git a/README.md b/README.md index eccaf40..210739d 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,9 @@ # Marble [![CMake Build](https://github.com/MarbleSolver/RCQP/actions/workflows/cmake-build.yml/badge.svg)](https://github.com/MarbleSolver/RCQP/actions/workflows/cmake-build.yml) +[![](https://img.shields.io/badge/docs-dev-blue.svg)](https://roboticexplorationlab.org/Marble/) A C++ solver for quadratic programs with linear complementarity constraints, with Python and Julia bindings. -This repository is currently a work-in-progress and will receive some significant clean up. If you want to get started we recommend checking out the Getting Started with [Julia](https://roboticexplorationlab.org/Marble/getting_started/julia/) or [Python](https://roboticexplorationlab.org/Marble/getting_started/python/) pages. We will try to make sure that the interfaces (Julia and Python) do not change, just the internals. - ## Prerequisites | Dependency | Required for | @@ -30,7 +29,7 @@ The project uses CMake presets defined in `CMakePresets.json`. Execute the comma cmake --preset python cmake --build --preset python ``` -Output: `build/python/marble.cpython--.so` +Output: `build/python/marble/_core.cpython--.so` (the compiled core inside the `marble` package) ### Julia bindings ```bash @@ -60,11 +59,59 @@ Inside the container: - `python` uses the venv at `/opt/venv`; `import marble` works immediately - Julia shared library is at `build/lib/libmarble_julia.so` -## Python: pip install (editable) -For a proper install into your virtual environment: +## Julia: use from your own environment + +The Julia bindings in `julia/` form a package named `Marble`. It is not registered, so +add it to your environment with `Pkg.develop`. -**1. Install system dependencies** (CMake, Eigen3, and nlohmann-json are still required at build time): +**1. Build the Julia shared library** (this needs `CxxWrap` available to the Julia that +CMake uses, see [Julia bindings](#julia-bindings) above): +```bash +julia -e 'using Pkg; Pkg.add("CxxWrap")' +cmake --preset julia +cmake --build --preset julia +``` +This writes `build/lib/libmarble_julia.{dylib,so}`, which `Marble.jl` loads on import. +By default it looks in this repo's `build/lib/`, so an in-repo build needs no extra +configuration. + +**2. Add the package to your environment and install its dependencies:** +```julia +using Pkg +Pkg.develop(path="/path/to/RCQP/julia") +Pkg.instantiate() # pulls in CxxWrap and the other dependencies +``` + +**3. Use it:** build a solver, set up the problem matrices, then solve. Matrices may be +dense or sparse, and solver options are passed to `setup!` as keyword arguments. +```julia +using Marble +using LinearAlgebra + +# min 1/2 x'x s.t. x1 + x2 = 2 -> x* = [1, 1] +solver = Marble.Solver() +Marble.setup!(solver, Matrix(1.0I, 2, 2), [0.0, 0.0]; J_eq = [1.0 1.0], b_eq = [-2.0]) +res = Marble.solve!(solver) +println(Marble.converged(res), " ", collect(Marble.z(res))) # true ~[1.0, 1.0] +``` + +If you build or move the shared library elsewhere, point the package at it with a +Preference, then restart Julia so the new path takes effect: +```julia +using Preferences +set_preferences!("Marble", "libmarble_julia_path" => "/abs/path/to/libmarble_julia") +``` +The path omits the file extension; `Marble.jl` appends the platform's `.dylib` / `.so`. + + +## Python: use from your own virtual environment + +marble is not on PyPI, so install it from a local checkout of this repository. This +works with any virtual environment, a fresh one or an existing project's. + +**1. Install the system build dependencies** (CMake, Eigen3, and nlohmann-json are +needed at build time): ```bash # macOS brew install cmake eigen nlohmann-json @@ -73,59 +120,78 @@ brew install cmake eigen nlohmann-json sudo apt install cmake libeigen3-dev nlohmann-json3-dev ``` -**2. Install Python build dependencies into your venv:** +**2. Create and activate a virtual environment:** ```bash -pip install scikit-build-core pybind11 +python -m venv .venv # or activate an existing environment +source .venv/bin/activate # Windows: .venv\Scripts\activate ``` -**3. Install marble:** +**3. Install marble from the repo's `python/` directory:** ```bash -pip install -e . --no-build-isolation +pip install /path/to/RCQP/python ``` +pip builds the C++ extension in an isolated environment, fetching `scikit-build-core` +and `pybind11` automatically, and pulls in `numpy` and `scipy`. Nothing else needs to +be installed first. -After this, `import marble` works from anywhere in that environment. +After this, `import marble` works wherever that environment is active: +```python +import numpy as np +import marble -## Python: usage without installing +# min 1/2 x'x + q'x -> x* = -q +solver = marble.Solver() +solver.setup(np.eye(2), np.array([1.0, 2.0])) +print(solver.solve().z) # [-1. -2.] +``` + +**Developing marble itself:** install the build tools into your environment and use an +editable install so source edits are picked up without reinstalling: +```bash +pip install scikit-build-core pybind11 +pip install -e /path/to/RCQP/python --no-build-isolation +``` + +## Python: usage without installing + +After building (`cmake --build --preset python`), the `marble` package is staged at +`build/python/marble/`. Add `build/python/` to your path at runtime, then import it: -Add `build/python/` to your path at runtime: ```python import sys -from pathlib import Path -# Insert the path to the `build` directory in PYTHONPATH -path_to_marble = ... -sys.path.insert(0, str(path_to_marble / "build" / "python")) +rcqp_root = "/path/to/RCQP" # your checkout of this repo +sys.path.insert(0, rcqp_root + "/build/python") import marble ``` ## VSCode: Python autocomplete -Autocomplete and type-checking are driven by the `.pyi` stubs in `typings/`, which are regenerated automatically each time you build `marble_python` (requires `pybind11-stubgen`). - -Two config files in this repo wire everything up for VSCode automatically: +The high-level API in `python/marble/__init__.py` is type-annotated, and the +compiled core ships a `python/marble/_core.pyi` stub that is regenerated each time +you build `marble_python` (requires `pybind11-stubgen`). Pylance reads both from the +`marble` package, so a single setting wires everything up: ### `.vscode/settings.json` ```json { "python.analysis.extraPaths": ["${workspaceFolder}/build/python"], - "python.analysis.stubPath": "${workspaceFolder}/typings", "python.analysis.useLibraryCodeForTypes": true } ``` -- `extraPaths`: tells Pylance where the compiled `.so` lives so it can be imported -- `stubPath`: tells Pylance where to find the `.pyi` stub file for type info +- `extraPaths`: tells Pylance where the staged `marble` package (with `_core.so` and `_core.pyi`) lives so it can be imported ### Install the Pylance extension -Install the **Pylance** extension in VSCode (`ms-python.vscode-pylance`). It picks up the above settings automatically. +Install the **Pylance** extension in VSCode (`ms-python.vscode-pylance`). It picks up the above setting automatically. If stubs are stale or missing, rebuild: ```bash cmake --build --preset python --target marble_python ``` -If `pybind11-stubgen` is not installed, CMake will warn but the build still succeeds. Ensure `pybind11-stubgen` is installed with: +If `pybind11-stubgen` is not installed, CMake will warn but the build still succeeds. Install it with: ```bash pip install pybind11-stubgen ```