From 40dafc0adbb638067bda7c005ba058cb1da3392b Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Sat, 4 Apr 2026 01:29:39 -0400 Subject: [PATCH 1/9] first version python bindings --- .gitignore | 3 +- .python-version | 1 + CMakeLists.txt | 193 +++++++++++++++++++++++++-------- CMakePresets.json | 42 +++++++ README.md | 153 ++++++++++++++++++++++++-- julia_bindings.cpp | 259 ++++++++++++++++++++++++++++++++++++++++++++ python_bindings.cpp | 242 +++++++++++++++++++++++++++++++++++++++++ tester.jl | 156 ++++++++++++++++++++++++++ tester.py | 47 ++++++++ typings/rcqp.pyi | 224 ++++++++++++++++++++++++++++++++++++++ wrapper.cpp | 216 ------------------------------------ 11 files changed, 1264 insertions(+), 272 deletions(-) create mode 100644 .python-version create mode 100644 CMakePresets.json create mode 100644 julia_bindings.cpp create mode 100644 python_bindings.cpp create mode 100644 tester.jl create mode 100644 tester.py create mode 100644 typings/rcqp.pyi delete mode 100644 wrapper.cpp diff --git a/.gitignore b/.gitignore index ab24a97..5cb03a8 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ .vscode -build/* \ No newline at end of file +build/* +.venv/ \ No newline at end of file diff --git a/.python-version b/.python-version new file mode 100644 index 0000000..bd28b9c --- /dev/null +++ b/.python-version @@ -0,0 +1 @@ +3.9 diff --git a/CMakeLists.txt b/CMakeLists.txt index 190e3bd..5352fa1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,67 +1,166 @@ -cmake_minimum_required(VERSION 3.10.2) -project(rcqp) +cmake_minimum_required(VERSION 3.28) # 3.28 needed for EXCLUDE_FROM_ALL in FetchContent_Declare + +project(rcqp LANGUAGES C CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) + set(CMAKE_POSITION_INDEPENDENT_CODE ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) +# --------------------------------------------------------------------------- +# Build options +# --------------------------------------------------------------------------- +option(RCQP_BUILD_PYTHON "Build Python bindings" ON) +option(RCQP_BUILD_JULIA "Build Julia bindings" ON) +option(RCQP_GENERATE_PYI "Generate .pyi stubs via pybind11-stubgen" ON) + +# --------------------------------------------------------------------------- +# Output directories (keeps build/ tidy) +# --------------------------------------------------------------------------- +set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) +set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) + +# --------------------------------------------------------------------------- # Dependencies -find_package(Eigen3 3.3 REQUIRED NO_MODULE) # provides Eigen3::Eigen -find_package(QDLDL) +# --------------------------------------------------------------------------- +include(FetchContent) + +find_package(Eigen3 3.3 REQUIRED NO_MODULE) find_package(nlohmann_json CONFIG REQUIRED) -if(NOT QDLDL_FOUND) - include(FetchContent) +# QDLDL – prefer a system install, otherwise fetch as static-only. +# EXCLUDE_FROM_ALL suppresses QDLDL's cmake --install rules so they don't +# interfere with pip install / scikit-build-core (requires CMake 3.28). +find_package(qdldl CONFIG QUIET) +if(NOT qdldl_FOUND) + set(QDLDL_BUILD_STATIC_LIB ON CACHE BOOL "" FORCE) + set(QDLDL_BUILD_SHARED_LIB OFF CACHE BOOL "" FORCE) + set(QDLDL_BUILD_DEMO_EXE OFF CACHE BOOL "" FORCE) + set(QDLDL_UNITTESTS OFF CACHE BOOL "" FORCE) FetchContent_Declare( - qdldl - GIT_REPOSITORY https://github.com/oxfordcontrol/QDLDL.git - GIT_TAG master + qdldl + GIT_REPOSITORY https://github.com/osqp/qdldl.git + GIT_TAG master + EXCLUDE_FROM_ALL ) FetchContent_MakeAvailable(qdldl) - # If the fetched project provides a target, use it; otherwise adjust below. endif() -# Create alias if needed for consistent naming -if(NOT TARGET qdldl::qdldl) - if(TARGET qdldl) - add_library(qdldl::qdldl ALIAS qdldl) - endif() +# Normalise QDLDL to a single interface target regardless of how it was found +add_library(rcqp_qdldl INTERFACE) +if(TARGET qdldl::qdldlstatic) + target_link_libraries(rcqp_qdldl INTERFACE qdldl::qdldlstatic) +elseif(TARGET qdldl::qdldl) + target_link_libraries(rcqp_qdldl INTERFACE qdldl::qdldl) +elseif(TARGET qdldlstatic) + target_link_libraries(rcqp_qdldl INTERFACE qdldlstatic) +elseif(TARGET qdldl) + target_link_libraries(rcqp_qdldl INTERFACE qdldl) +else() + message(FATAL_ERROR "Could not find a usable QDLDL target") endif() -find_package(JlCxx) -get_target_property(JlCxx_location JlCxx::cxxwrap_julia LOCATION) -get_filename_component(JlCxx_location ${JlCxx_location} DIRECTORY) -set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib;${JlCxx_location}") +# --------------------------------------------------------------------------- +# Core library (static → gets bundled into each wrapper, no rpath issues) +# --------------------------------------------------------------------------- +add_library(rcqp STATIC + src/solver.cpp + src/problem.cpp + src/workspace.cpp + src/filter.cpp +) +target_include_directories(rcqp PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include) +target_link_libraries(rcqp PUBLIC Eigen3::Eigen nlohmann_json::nlohmann_json rcqp_qdldl) -message(STATUS "Found JlCxx at ${JlCxx_location}") +# --------------------------------------------------------------------------- +# Python bindings (cmake -DRCQP_BUILD_PYTHON=ON — default ON) +# --------------------------------------------------------------------------- +if(RCQP_BUILD_PYTHON) + find_package(Python COMPONENTS Interpreter Development.Module REQUIRED) -# Setup library -add_library(rcqp SHARED - ${CMAKE_CURRENT_SOURCE_DIR}/src/solver.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/problem.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/workspace.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/filter.cpp -) + find_package(pybind11 CONFIG QUIET) + if(NOT pybind11_FOUND) + FetchContent_Declare( + pybind11 + GIT_REPOSITORY https://github.com/pybind/pybind11.git + GIT_TAG v2.13.6 + ) + FetchContent_MakeAvailable(pybind11) + endif() -target_include_directories(rcqp - PUBLIC - ${CMAKE_CURRENT_SOURCE_DIR}/include -) + pybind11_add_module(rcqp_python MODULE python_bindings.cpp) + target_link_libraries(rcqp_python PRIVATE rcqp) + set_target_properties(rcqp_python PROPERTIES + OUTPUT_NAME rcqp + LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/python + ) + + # .pyi stub generation (for VSCode / Pylance autocomplete) + if(RCQP_GENERATE_PYI) + find_program(PYBIND11_STUBGEN pybind11-stubgen) + if(PYBIND11_STUBGEN) + add_custom_command(TARGET rcqp_python POST_BUILD + COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_SOURCE_DIR}/typings + COMMAND ${CMAKE_COMMAND} -E env + PYTHONPATH=${CMAKE_BINARY_DIR}/python + ${PYBIND11_STUBGEN} rcqp + --output-dir ${CMAKE_SOURCE_DIR}/typings + --numpy-array-remove-parameters + --ignore-unresolved-names "^(m|n)$" + COMMENT "Generating rcqp.pyi stubs" + VERBATIM + ) + else() + message(WARNING + "pybind11-stubgen not found — .pyi stubs won't be regenerated. " + "Run: pip install pybind11-stubgen") + endif() + endif() + + install(TARGETS rcqp_python LIBRARY DESTINATION .) +endif() -target_compile_features(rcqp PUBLIC cxx_std_17) - -# Wrapper library for Julia -add_library(rcqp_wrapper SHARED wrapper.cpp) -target_link_libraries(rcqp_wrapper - PRIVATE - rcqp - JlCxx::cxxwrap_julia) - -# Link dependencies -target_link_libraries(rcqp - PUBLIC - Eigen3::Eigen - qdldl - nlohmann_json::nlohmann_json -) \ No newline at end of file +# --------------------------------------------------------------------------- +# Julia bindings (cmake -DRCQP_BUILD_JULIA=ON) +# --------------------------------------------------------------------------- +if(RCQP_BUILD_JULIA) + # Auto-detect JlCxx prefix by asking Julia where CxxWrap installed it. + # This works on any machine that has Julia + CxxWrap.jl without any + # hardcoded paths. Users can still override with -DJlCxx_DIR=. + if(NOT DEFINED JlCxx_DIR) + find_program(JULIA_EXECUTABLE julia) + if(JULIA_EXECUTABLE) + execute_process( + COMMAND ${JULIA_EXECUTABLE} --startup-file=no -e + "using CxxWrap; print(CxxWrap.prefix_path())" + OUTPUT_VARIABLE _jlcxx_prefix + ERROR_QUIET + OUTPUT_STRIP_TRAILING_WHITESPACE + ) + if(_jlcxx_prefix) + set(JlCxx_DIR "${_jlcxx_prefix}/lib/cmake/JlCxx" + CACHE PATH "Path to JlCxxConfig.cmake (auto-detected from Julia)") + message(STATUS "Auto-detected JlCxx_DIR: ${JlCxx_DIR}") + else() + message(FATAL_ERROR + "Could not detect JlCxx prefix from Julia.\n" + "Make sure CxxWrap.jl is installed: julia -e 'using Pkg; Pkg.add(\"CxxWrap\")'\n" + "Or set JlCxx_DIR manually: cmake -DJlCxx_DIR=") + endif() + else() + message(FATAL_ERROR + "julia executable not found.\n" + "Install Julia from https://julialang.org/downloads or set JlCxx_DIR manually.") + endif() + endif() + + find_package(JlCxx REQUIRED) + add_library(rcqp_julia SHARED julia_bindings.cpp) + target_link_libraries(rcqp_julia PRIVATE rcqp JlCxx::cxxwrap_julia) + set_target_properties(rcqp_julia PROPERTIES + LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib + ) +endif() diff --git a/CMakePresets.json b/CMakePresets.json new file mode 100644 index 0000000..5fa18f0 --- /dev/null +++ b/CMakePresets.json @@ -0,0 +1,42 @@ +{ + "version": 3, + "configurePresets": [ + { + "name": "python", + "displayName": "Python bindings only", + "binaryDir": "${sourceDir}/build", + "cacheVariables": { + "CMAKE_BUILD_TYPE": "Release", + "RCQP_BUILD_PYTHON": "ON", + "RCQP_BUILD_JULIA": "OFF", + "RCQP_GENERATE_PYI": "ON" + } + }, + { + "name": "julia", + "displayName": "Julia bindings only (requires Julia + CxxWrap.jl)", + "binaryDir": "${sourceDir}/build", + "cacheVariables": { + "CMAKE_BUILD_TYPE": "Release", + "RCQP_BUILD_PYTHON": "OFF", + "RCQP_BUILD_JULIA": "ON" + } + }, + { + "name": "all", + "displayName": "Python + Julia bindings (requires Julia + CxxWrap.jl)", + "binaryDir": "${sourceDir}/build", + "cacheVariables": { + "CMAKE_BUILD_TYPE": "Release", + "RCQP_BUILD_PYTHON": "ON", + "RCQP_BUILD_JULIA": "ON", + "RCQP_GENERATE_PYI": "ON" + } + } + ], + "buildPresets": [ + { "name": "python", "configurePreset": "python" }, + { "name": "julia", "configurePreset": "julia" }, + { "name": "all", "configurePreset": "all" } + ] +} diff --git a/README.md b/README.md index e694b1c..5791026 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,154 @@ # RCQP -A C++ solver for quadratic programs with linear complementarity constraints. +A C++ solver for quadratic programs with linear complementarity constraints, with Python and Julia bindings. -## Installation & Compilation +## Prerequisites -### MacOS and Ubuntu -First install the CxxWrap Julia package: +| Dependency | Required for | Install | +|---|---|---| +| CMake ≥ 3.28 | all | [cmake.org](https://cmake.org/download/) or `brew install cmake` | +| Eigen3 | all | `brew install eigen` / `apt install libeigen3-dev` | +| nlohmann-json | all | `brew install nlohmann-json` / `apt install nlohmann-json3-dev` | +| Python ≥ 3.8 + dev headers | Python bindings | [python.org](https://www.python.org/) | +| pybind11-stubgen | `.pyi` autocomplete stubs | `pip install pybind11-stubgen` | +| Julia ≥ 1.9 + CxxWrap.jl | Julia bindings | see below | + +**Julia setup** (only needed for Julia bindings): ```bash julia -e 'using Pkg; Pkg.add("CxxWrap")' ``` -We require at least v0.17.5, which can be checked after installation using `Pkg.status("CxxWrap")` -Inside the repo root, run the following commands: +## Building + +The project uses CMake presets defined in `CMakePresets.json`. + +### Python bindings (default) +```bash +cmake --preset python +cmake --build --preset python +``` +Output: `build/python/rcqp.cpython--.so` + +### Julia bindings +```bash +cmake --preset julia +cmake --build --preset julia +``` +CMake auto-detects JlCxx by calling `julia -e 'using CxxWrap; print(CxxWrap.prefix_path())'` — no manual paths required. + +Output: `build/lib/librcqp_julia.dylib` (macOS) / `build/lib/librcqp_julia.so` (Linux) + +### Both +```bash +cmake --preset all +cmake --build --preset all +``` + + + +## Python — pip install (editable) + +For a proper install into your virtual environment (no `sys.path` hacks): +```bash +pip install scikit-build-core pybind11-stubgen +pip install -e . --no-build-isolation +``` + +After this, `import rcqp` works from anywhere in that environment. + +## Python — usage without installing + +Add `build/python/` to your path at runtime: +```python +import sys +from pathlib import Path +sys.path.insert(0, str(Path(__file__).resolve().parent / "build" / "python")) + +import rcqp +import numpy as np + +opts = rcqp.SolverOptions() +opts.max_iters = 500 + +solver = rcqp.Solver(opts) +solver.set_problem(prob, scaling) +converged = solver.solve(opts) + +ws = solver.get_workspace() +print(ws.z) # primal solution +``` + +## VSCode — Python autocomplete + +Autocomplete and type-checking are driven by the `.pyi` stubs in `typings/`, which are regenerated automatically each time you build `rcqp_python` (requires `pybind11-stubgen`). + +Two config files in this repo wire everything up for VSCode automatically: + +### `.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 + +### Install the Pylance extension +Install the **Pylance** extension in VSCode (`ms-python.vscode-pylance`). It picks up the above settings automatically. + +If stubs are stale or missing, rebuild: +```bash +cmake --build --preset python --target rcqp_python +``` + +If `pybind11-stubgen` is not installed, CMake will warn but the build still succeeds — you just won't get updated stubs. +```bash +pip install pybind11-stubgen +``` + +## Julia — usage + +```julia +using CxxWrap + +module RCQP + using CxxWrap + @wrapmodule(() -> joinpath(@__DIR__, "build/lib/librcqp_julia")) + function __init__() + @initcxx + end +end + +prob = RCQP.Problem(H, g, 0.0, J_eq, c_eq, J_ineq, c_ineq, J_comp, c_comp) + +opts = RCQP.SolverOptions() +RCQP.max_iters!(opts, 500) + +solver = RCQP.Solver(opts) +RCQP.set_problem(solver, prob, ones(nz)) +converged = RCQP.solve(solver, opts) + +ws = RCQP.get_workspace(solver) +z = copy(RCQP.z(ws)) +``` + +See `tester.jl` for a complete working example. diff --git a/julia_bindings.cpp b/julia_bindings.cpp new file mode 100644 index 0000000..f34355e --- /dev/null +++ b/julia_bindings.cpp @@ -0,0 +1,259 @@ +#include +#include "solver.h" + +// Filter::Entry is a plain struct; tell CxxWrap not to treat it as a mirrored type +// so we can register it with add_type and attach methods to it. +namespace jlcxx { + template<> struct IsMirroredType : std::false_type {}; +} + +// --------------------------------------------------------------------------- +// Eigen ↔ Julia array helpers +// --------------------------------------------------------------------------- + +template +Eigen::Matrix +to_eigen(jlcxx::ArrayRef& arr, int rows, int cols) { + if (rows == 0 || cols == 0) + return Eigen::Matrix(rows, cols); + return Eigen::Matrix::Map(arr.data(), rows, cols); +} + +template +Eigen::Matrix +to_eigen(jlcxx::ArrayRef& arr) { + if (arr.size() == 0) + return Eigen::Matrix(0); + return Eigen::Matrix::Map(arr.data(), arr.size()); +} + +template +jlcxx::ArrayRef to_julia(Eigen::Matrix& vec) { + return jlcxx::make_julia_array(vec.data(), vec.size()); +} + +template +jlcxx::ArrayRef to_julia(Eigen::Map>& vec) { + return jlcxx::make_julia_array(vec.data(), vec.size()); +} + +// --------------------------------------------------------------------------- +// Module definition +// --------------------------------------------------------------------------- + +JLCXX_MODULE define_julia_module(jlcxx::Module& mod) { + + // ----------------------------------------------------------------------- + // Problem + // ----------------------------------------------------------------------- + mod.add_type("Problem") + .constructor([](jlcxx::ArrayRef cost_hessian, + jlcxx::ArrayRef cost_gradient, + double cost_const, + jlcxx::ArrayRef J_eq, jlcxx::ArrayRef c_eq, + jlcxx::ArrayRef J_ineq, jlcxx::ArrayRef c_ineq, + jlcxx::ArrayRef J_comp, jlcxx::ArrayRef c_comp) { + int nz = cost_gradient.size(); + int n_eq = c_eq.size(); + int n_ineq = c_ineq.size(); + int n_comp = c_comp.size(); + return new Problem( + to_eigen(cost_hessian, nz, nz), to_eigen(cost_gradient), cost_const, + to_eigen(J_eq, n_eq, nz), to_eigen(c_eq), + to_eigen(J_ineq, n_ineq, nz), to_eigen(c_ineq), + to_eigen(J_comp, n_comp, nz), to_eigen(c_comp)); + }) + .method("nz", [](const Problem& p) { return p.nz; }) + .method("n_eq", [](const Problem& p) { return p.n_eq; }) + .method("n_ineq", [](const Problem& p) { return p.n_ineq; }) + .method("n_comp", [](const Problem& p) { return p.n_comp; }); + + // ----------------------------------------------------------------------- + // SolverOptions + // ----------------------------------------------------------------------- + #define OPTION_RW(name, T) \ + .method(#name, [](const Solver::Options& o) { return o.name; }) \ + .method(#name "!", [](Solver::Options& o, T v) { o.name = v; }) + + mod.add_type("SolverOptions") + .constructor() + OPTION_RW(convergence_kkt_norm, double) + OPTION_RW(convergence_eq_violation, double) + OPTION_RW(convergence_ineq_violation, double) + OPTION_RW(convergence_comp_violation, double) + OPTION_RW(outer_step_kkt_norm, double) + OPTION_RW(penalty_initial, double) + OPTION_RW(penalty_max, double) + OPTION_RW(penalty_scaling, double) + OPTION_RW(relaxation_initial, double) + OPTION_RW(relaxation_min, double) + OPTION_RW(relaxation_scaling, double) + OPTION_RW(max_iters, int) + OPTION_RW(max_iters_linesearch, int) + OPTION_RW(gamma_objective, double) + OPTION_RW(gamma_constraint, double) + .method("output_dir", [](const Solver::Options& o) { return o.output_dir.string(); }) + .method("output_dir!", [](Solver::Options& o, const std::string& v) { o.output_dir = v; }); + + #undef OPTION_RW + + // ----------------------------------------------------------------------- + // FilterEntry + // ----------------------------------------------------------------------- + mod.add_type("FilterEntry") + .constructor() + .constructor([](double feas, double merit) { + Filter::Entry* e = new Filter::Entry(); e->feas = feas; e->merit = merit; return e; + }) + .method("feas", [](const Filter::Entry& e) { return e.feas; }) + .method("feas!", [](Filter::Entry& e, double v) { e.feas = v; }) + .method("merit", [](const Filter::Entry& e) { return e.merit; }) + .method("merit!", [](Filter::Entry& e, double v) { e.merit = v; }); + + // ----------------------------------------------------------------------- + // Filter + // ----------------------------------------------------------------------- + mod.add_type("Filter") + .constructor() + .constructor([](double gamma_objective, double gamma_constraint) { + return new Filter(gamma_objective, gamma_constraint); + }) + .method("clear", &Filter::clear) + .method("size", [](const Filter& f) { return f.entries.size(); }) + // Returns a flat vector [feas0, merit0, feas1, merit1, ...] + // Use reshape(entries(f), 2, :) in Julia to get a 2×n matrix + .method("entries", [](const Filter& f) { + std::vector out; + out.reserve(f.entries.size() * 2); + for (const auto& e : f.entries) { + out.push_back(e.feas); + out.push_back(e.merit); + } + return out; + }) + // sufficient_progress returns std::pair which CxxWrap can't wrap directly, + // so expose each component as its own method + .method("sufficient_feas_progress", [](Filter& f, const Filter::Entry& c, const Filter::Entry& e) { + return f.sufficient_progress(c, e).first; + }) + .method("sufficient_merit_progress", [](Filter& f, const Filter::Entry& c, const Filter::Entry& e) { + return f.sufficient_progress(c, e).second; + }) + .method("candidate_acceptable", &Filter::candidate_acceptable) + .method("candidate_dominated", &Filter::candidate_dominated) + .method("acceptable", &Filter::acceptable) + .method("update", &Filter::update); + + // ----------------------------------------------------------------------- + // Workspace (read-only accessors) + // ----------------------------------------------------------------------- + mod.add_type("Workspace") + // Full solution and decomposed views + .method("solution", [](Workspace& w) { return to_julia(w.solution); }) + .method("z", [](Workspace& w) { return to_julia(w.z); }) + .method("s_ineq", [](Workspace& w) { return to_julia(w.s_ineq); }) + .method("s_comp", [](Workspace& w) { return to_julia(w.s_comp); }) + .method("m_eq", [](Workspace& w) { return to_julia(w.m_eq); }) + .method("m_ineq", [](Workspace& w) { return to_julia(w.m_ineq); }) + .method("m_comp", [](Workspace& w) { return to_julia(w.m_comp); }) + // Multiplier estimates + .method("m_eq_est", [](Workspace& w) { return to_julia(w.m_eq_est); }) + .method("m_ineq_est", [](Workspace& w) { return to_julia(w.m_ineq_est); }) + .method("m_comp_est", [](Workspace& w) { return to_julia(w.m_comp_est); }) + // Residuals + .method("kkt_residual", [](Workspace& w) { return to_julia(w.kkt_residual); }) + .method("residual_eq", [](Workspace& w) { return to_julia(w.residual_eq); }) + .method("residual_ineq", [](Workspace& w) { return to_julia(w.residual_ineq); }) + .method("residual_comp", [](Workspace& w) { return to_julia(w.residual_comp); }) + // Scalar parameters + .method("relax_param", [](const Workspace& w) { return w.relax_param; }) + .method("penalty_param", [](const Workspace& w) { return w.penalty_param; }) + // Newton step + .method("newton_step", [](Workspace& w) { return to_julia(w.newton_step); }) + // KKT system: returns (rows, cols, colptr, rowval, nzval) tuple + .method("kkt_system", [](Workspace& w) { + SMat& kkt = w.kkt_system; + kkt.makeCompressed(); + auto colptr = jlcxx::make_julia_array(kkt.outerIndexPtr(), kkt.outerSize() + 1); + auto rowval = jlcxx::make_julia_array(kkt.innerIndexPtr(), kkt.nonZeros()); + auto nzval = jlcxx::make_julia_array(kkt.valuePtr(), kkt.nonZeros()); + return std::make_tuple((int)kkt.rows(), (int)kkt.cols(), colptr, rowval, nzval); + }) + // QDLDL factorization data + .method("D", [](Workspace& w) { return to_julia(w.D); }) + .method("amd_perm_vec", [](Workspace& w) { return to_julia(w.amd_perm_vec); }) + .method("amd_iperm_vec",[](Workspace& w) { return to_julia(w.amd_iperm_vec); }); + + // ----------------------------------------------------------------------- + // Solver + // ----------------------------------------------------------------------- + mod.add_type("Solver") + .constructor() + .constructor([](const Solver::Options& opts) { return new Solver(opts); }) + // Problem setup + .method("set_problem", [](Solver& s, Problem& prob, jlcxx::ArrayRef scaling) { + s.set_problem(prob, to_eigen(scaling)); + }) + .method("get_problem", &Solver::get_problem) + // Main solve + .method("solve", &Solver::solve) + .method("convergence", &Solver::convergence) + // Workspace / filter access + .method("get_workspace", &Solver::get_workspace) + .method("get_filter", &Solver::get_filter) + // Retraction maps + .method("retract", [](Solver& s, jlcxx::ArrayRef v, double sqrt_relax_param) { + Vec r = s.retract(to_eigen(v), sqrt_relax_param); + return std::vector(r.data(), r.data() + r.size()); + }) + .method("retract_deriv", [](Solver& s, jlcxx::ArrayRef v, double sqrt_relax_param) { + Vec r = s.retract_deriv(to_eigen(v), sqrt_relax_param); + return std::vector(r.data(), r.data() + r.size()); + }) + .method("retract_second_deriv", [](Solver& s, jlcxx::ArrayRef v, double sqrt_relax_param) { + Vec r = s.retract_second_deriv(to_eigen(v), sqrt_relax_param); + return std::vector(r.data(), r.data() + r.size()); + }) + // KKT updates + .method("update_KKT_residual", &Solver::update_KKT_residual) + .method("update_KKT_system", &Solver::update_KKT_system) + .method("update_KKT_ineq", [](Solver& s, jlcxx::ArrayRef s_ineq, double sqrt_relax_param) { + s.update_KKT_ineq(to_eigen(s_ineq), sqrt_relax_param); + }) + .method("update_KKT_comp", [](Solver& s, jlcxx::ArrayRef s_comp, jlcxx::ArrayRef m_comp, double sqrt_relax_param) { + s.update_KKT_comp(to_eigen(s_comp), to_eigen(m_comp), sqrt_relax_param); + }) + .method("update_KKT_penalty", &Solver::update_KKT_penalty) + .method("update_KKT_primal_regularizer", &Solver::update_KKT_primal_regularizer) + // Factorization and backsolve + .method("initialize_kkt_sparsity", &Solver::initialize_kkt_sparsity) + .method("compute_amd_ordering", &Solver::compute_amd_ordering) + .method("analytical_factorization", &Solver::analytical_factorization) + .method("numerical_factorization", &Solver::numerical_factorization) + .method("check_inertia", &Solver::check_inertia) + .method("backsolve", &Solver::backsolve) + // Linesearch + .method("filter_linesearch", &Solver::filter_linesearch) + .method("entry_from_solution", [](const Solver& s, double sqrt_relax_param, double inv_penalty_param) { + Filter::Entry e = s.entry_from_solution(sqrt_relax_param, inv_penalty_param); + return std::make_tuple(e.feas, e.merit); + }) + // Dimensions + .method("n_primals", [](const Solver& s) { return s.n_primals; }) + .method("n_duals", [](const Solver& s) { return s.n_duals; }) + .method("n_vars", [](const Solver& s) { return s.n_vars; }) + // KKT index vectors + .method("z_inds", [](Solver& s) { return to_julia(s.z_inds); }) + .method("s_ineq_inds", [](Solver& s) { return to_julia(s.s_ineq_inds); }) + .method("s_comp_inds", [](Solver& s) { return to_julia(s.s_comp_inds); }) + .method("m_eq_inds", [](Solver& s) { return to_julia(s.m_eq_inds); }) + .method("m_ineq_inds", [](Solver& s) { return to_julia(s.m_ineq_inds); }) + .method("m_comp_inds", [](Solver& s) { return to_julia(s.m_comp_inds); }) + .method("comp_L_inds", [](Solver& s) { return to_julia(s.comp_L_inds); }) + .method("comp_R_inds", [](Solver& s) { return to_julia(s.comp_R_inds); }) + .method("z_z_inds", [](Solver& s) { return to_julia(s.z_z_inds); }) + .method("s_ineq_s_ineq_inds", [](Solver& s) { return to_julia(s.s_ineq_s_ineq_inds); }) + .method("s_ineq_m_ineq_inds", [](Solver& s) { return to_julia(s.s_ineq_m_ineq_inds); }) + .method("s_comp_s_comp_inds", [](Solver& s) { return to_julia(s.s_comp_s_comp_inds); }) + .method("s_comp_m_comp_inds", [](Solver& s) { return to_julia(s.s_comp_m_comp_inds); }); +} diff --git a/python_bindings.cpp b/python_bindings.cpp new file mode 100644 index 0000000..c195d1c --- /dev/null +++ b/python_bindings.cpp @@ -0,0 +1,242 @@ +#include +#include +#include +#include "solver.h" + +namespace py = pybind11; + +// Helper: copy Eigen::Map → Vec so pybind11/eigen can convert to numpy +static Vec copy_map(const Eigen::Map& m) { return Vec(m); } + +PYBIND11_MODULE(rcqp, m) { + m.doc() = "RCQP: constrained optimization solver with complementarity constraints"; + + // ------------------------------------------------------------------------- + // Problem + // ------------------------------------------------------------------------- + py::class_(m, "Problem") + // Dense numpy-array constructor (mirrors Julia wrapper) + .def(py::init([](const Mat& cost_hessian, const Vec& cost_gradient, double cost_const, + const Mat& J_eq, const Vec& c_eq, + const Mat& J_ineq, const Vec& c_ineq, + const Mat& J_comp, const Vec& c_comp) { + return new Problem(cost_hessian, cost_gradient, cost_const, + J_eq, c_eq, J_ineq, c_ineq, J_comp, c_comp); + }), + py::arg("cost_hessian"), py::arg("cost_gradient"), py::arg("cost_const"), + py::arg("J_eq"), py::arg("c_eq"), + py::arg("J_ineq"), py::arg("c_ineq"), + py::arg("J_comp"), py::arg("c_comp")) + // Sparse constructor – accepts Eigen::SparseMatrix which pybind11 converts + // from scipy.sparse.csc_matrix automatically + .def(py::init([](const SMat& cost_hessian, const Vec& cost_gradient, double cost_const, + const SMat& J_eq, const Vec& c_eq, + const SMat& J_ineq, const Vec& c_ineq, + const SMat& J_comp, const Vec& c_comp) { + return new Problem(cost_hessian, cost_gradient, cost_const, + J_eq, c_eq, J_ineq, c_ineq, J_comp, c_comp); + }), + py::arg("cost_hessian"), py::arg("cost_gradient"), py::arg("cost_const"), + py::arg("J_eq"), py::arg("c_eq"), + py::arg("J_ineq"), py::arg("c_ineq"), + py::arg("J_comp"), py::arg("c_comp")) + .def_readonly("nz", &Problem::nz) + .def_readonly("n_eq", &Problem::n_eq) + .def_readonly("n_ineq", &Problem::n_ineq) + .def_readonly("n_comp", &Problem::n_comp) + .def_property_readonly("cost_gradient", [](const Problem& p) { return p.cost_gradient; }) + .def_readonly("cost_const", &Problem::cost_const); + + // ------------------------------------------------------------------------- + // Solver::Options + // ------------------------------------------------------------------------- + py::class_(m, "SolverOptions") + .def(py::init<>()) + .def_readwrite("convergence_kkt_norm", &Solver::Options::convergence_kkt_norm) + .def_readwrite("convergence_eq_violation", &Solver::Options::convergence_eq_violation) + .def_readwrite("convergence_ineq_violation", &Solver::Options::convergence_ineq_violation) + .def_readwrite("convergence_comp_violation", &Solver::Options::convergence_comp_violation) + .def_readwrite("outer_step_kkt_norm", &Solver::Options::outer_step_kkt_norm) + .def_readwrite("penalty_initial", &Solver::Options::penalty_initial) + .def_readwrite("penalty_max", &Solver::Options::penalty_max) + .def_readwrite("penalty_scaling", &Solver::Options::penalty_scaling) + .def_readwrite("relaxation_initial", &Solver::Options::relaxation_initial) + .def_readwrite("relaxation_min", &Solver::Options::relaxation_min) + .def_readwrite("relaxation_scaling", &Solver::Options::relaxation_scaling) + .def_readwrite("max_iters", &Solver::Options::max_iters) + .def_readwrite("max_iters_linesearch", &Solver::Options::max_iters_linesearch) + .def_readwrite("gamma_objective", &Solver::Options::gamma_objective) + .def_readwrite("gamma_constraint", &Solver::Options::gamma_constraint) + .def_property("output_dir", + [](const Solver::Options& o) { return o.output_dir.string(); }, + [](Solver::Options& o, const std::string& v) { o.output_dir = v; }) + .def("__repr__", [](const Solver::Options& o) { + return ""; + }); + + // ------------------------------------------------------------------------- + // Filter::Entry + // ------------------------------------------------------------------------- + py::class_(m, "FilterEntry") + .def(py::init<>()) + .def_readwrite("feas", &Filter::Entry::feas) + .def_readwrite("merit", &Filter::Entry::merit) + .def("__repr__", [](const Filter::Entry& e) { + return ""; + }); + + // ------------------------------------------------------------------------- + // Filter + // ------------------------------------------------------------------------- + py::class_(m, "Filter") + .def(py::init<>()) + .def(py::init(), + py::arg("gamma_objective"), py::arg("gamma_constraint")) + .def("clear", &Filter::clear) + .def_property_readonly("size", + [](const Filter& f) { return f.entries.size(); }) + .def_property_readonly("entries", + [](const Filter& f) { + std::vector> out; + out.reserve(f.entries.size()); + for (const auto& e : f.entries) + out.emplace_back(e.feas, e.merit); + return out; + }) + .def("sufficient_progress", &Filter::sufficient_progress, + py::arg("candidate"), py::arg("entry")) + .def("candidate_acceptable", &Filter::candidate_acceptable, + py::arg("candidate"), py::arg("entry")) + .def("candidate_dominated", &Filter::candidate_dominated, + py::arg("candidate"), py::arg("entry")) + .def("acceptable", &Filter::acceptable, py::arg("candidate")) + .def("update", &Filter::update, py::arg("new_entry")); + + // ------------------------------------------------------------------------- + // Workspace (read-only view into solver internals) + // ------------------------------------------------------------------------- + py::class_(m, "Workspace") + // Solution components – copy Maps to plain Vec so numpy owns the data + .def_property_readonly("solution", + [](const Workspace& w) { return w.solution; }) + .def_property_readonly("z", + [](const Workspace& w) { return copy_map(w.z); }) + .def_property_readonly("s_ineq", + [](const Workspace& w) { return copy_map(w.s_ineq); }) + .def_property_readonly("s_comp", + [](const Workspace& w) { return copy_map(w.s_comp); }) + .def_property_readonly("m_eq", + [](const Workspace& w) { return copy_map(w.m_eq); }) + .def_property_readonly("m_ineq", + [](const Workspace& w) { return copy_map(w.m_ineq); }) + .def_property_readonly("m_comp", + [](const Workspace& w) { return copy_map(w.m_comp); }) + // Multiplier estimates + .def_property_readonly("m_eq_est", + [](const Workspace& w) { return w.m_eq_est; }) + .def_property_readonly("m_ineq_est", + [](const Workspace& w) { return w.m_ineq_est; }) + .def_property_readonly("m_comp_est", + [](const Workspace& w) { return w.m_comp_est; }) + // Residuals + .def_property_readonly("kkt_residual", + [](const Workspace& w) { return w.kkt_residual; }) + .def_property_readonly("residual_eq", + [](const Workspace& w) { return w.residual_eq; }) + .def_property_readonly("residual_ineq", + [](const Workspace& w) { return w.residual_ineq; }) + .def_property_readonly("residual_comp", + [](const Workspace& w) { return w.residual_comp; }) + // Scalar parameters + .def_property_readonly("relax_param", + [](const Workspace& w) { return w.relax_param; }) + .def_property_readonly("penalty_param", + [](const Workspace& w) { return w.penalty_param; }) + // Newton step + .def_property_readonly("newton_step", + [](const Workspace& w) { return w.newton_step; }); + + // ------------------------------------------------------------------------- + // Solver + // ------------------------------------------------------------------------- + py::class_(m, "Solver") + .def(py::init<>()) + .def(py::init(), py::arg("options")) + // Problem setup + .def("set_problem", [](Solver& s, Problem& prob, const Vec& scaling) { + s.set_problem(prob, scaling); + }, + py::arg("problem"), py::arg("scaling")) + .def("get_problem", &Solver::get_problem, + py::return_value_policy::reference_internal) + // Main solve + .def("solve", &Solver::solve, py::arg("options")) + .def("convergence", &Solver::convergence, py::arg("options")) + // Workspace / filter access + .def("get_workspace", &Solver::get_workspace, + py::return_value_policy::reference_internal) + .def("get_filter", &Solver::get_filter, + py::return_value_policy::reference_internal) + // Retraction maps + .def("retract", &Solver::retract, + py::arg("s"), py::arg("sqrt_relax_param")) + .def("retract_deriv", &Solver::retract_deriv, + py::arg("s"), py::arg("sqrt_relax_param")) + .def("retract_second_deriv", &Solver::retract_second_deriv, + py::arg("s"), py::arg("sqrt_relax_param")) + // KKT updates + .def("update_KKT_residual", &Solver::update_KKT_residual, + py::arg("sqrt_relax_param"), py::arg("inv_penalty_param")) + .def("update_KKT_system", &Solver::update_KKT_system, + py::arg("sqrt_relax_param"), py::arg("inv_penalty_param")) + .def("update_KKT_ineq", [](Solver& s, const Vec& s_ineq, double sqrt_relax_param) { + s.update_KKT_ineq(s_ineq, sqrt_relax_param); + }, + py::arg("s_ineq"), py::arg("sqrt_relax_param")) + .def("update_KKT_comp", [](Solver& s, const Vec& s_comp, const Vec& m_comp, double sqrt_relax_param) { + s.update_KKT_comp(s_comp, m_comp, sqrt_relax_param); + }, + py::arg("s_comp"), py::arg("m_comp"), py::arg("sqrt_relax_param")) + .def("update_KKT_penalty", &Solver::update_KKT_penalty, + py::arg("inv_penalty_param")) + .def("update_KKT_primal_regularizer", &Solver::update_KKT_primal_regularizer, + py::arg("reg")) + // Factorization and solve + .def("initialize_kkt_sparsity", &Solver::initialize_kkt_sparsity) + .def("compute_amd_ordering", &Solver::compute_amd_ordering) + .def("analytical_factorization",&Solver::analytical_factorization) + .def("numerical_factorization", &Solver::numerical_factorization) + .def("check_inertia", &Solver::check_inertia) + .def("backsolve", &Solver::backsolve) + // Linesearch + .def("filter_linesearch", &Solver::filter_linesearch, + py::arg("sqrt_relax_param"), py::arg("inv_penalty_param"), py::arg("max_iters")) + .def("entry_from_solution", [](const Solver& s, double sqrt_relax_param, double inv_penalty_param) { + Filter::Entry e = s.entry_from_solution(sqrt_relax_param, inv_penalty_param); + return py::make_tuple(e.feas, e.merit); + }, + py::arg("sqrt_relax_param"), py::arg("inv_penalty_param")) + // Dimension info + .def_readonly("n_primals", &Solver::n_primals) + .def_readonly("n_duals", &Solver::n_duals) + .def_readonly("n_vars", &Solver::n_vars) + // KKT index vectors (returned as numpy int arrays) + .def_property_readonly("z_inds", + [](const Solver& s) { return s.z_inds; }) + .def_property_readonly("s_ineq_inds", + [](const Solver& s) { return s.s_ineq_inds; }) + .def_property_readonly("s_comp_inds", + [](const Solver& s) { return s.s_comp_inds; }) + .def_property_readonly("m_eq_inds", + [](const Solver& s) { return s.m_eq_inds; }) + .def_property_readonly("m_ineq_inds", + [](const Solver& s) { return s.m_ineq_inds; }) + .def_property_readonly("m_comp_inds", + [](const Solver& s) { return s.m_comp_inds; }) + .def_property_readonly("comp_L_inds", + [](const Solver& s) { return s.comp_L_inds; }) + .def_property_readonly("comp_R_inds", + [](const Solver& s) { return s.comp_R_inds; }); +} diff --git a/tester.jl b/tester.jl new file mode 100644 index 0000000..a046486 --- /dev/null +++ b/tester.jl @@ -0,0 +1,156 @@ +using CxxWrap +using LinearAlgebra +using SparseArrays + +# --------------------------------------------------------------------------- +# Load the module +# --------------------------------------------------------------------------- +module RCQP + using CxxWrap + @wrapmodule(() -> joinpath(@__DIR__, "build/librcqp_julia")) + function __init__() + @initcxx + end +end + +# --------------------------------------------------------------------------- +# Helper: reconstruct a Julia SparseMatrixCSC from the kkt_system tuple +# --------------------------------------------------------------------------- +function get_kkt(solver) + ws = RCQP.get_workspace(solver) + nr, nc, colptr, rowval, nzval = RCQP.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 = RCQP.Problem(H, g, 0.0, J_empty, c_empty, J_empty, c_empty, J_empty, c_empty) +@assert RCQP.nz(prob) == 2 +@assert RCQP.n_eq(prob) == 0 +@assert RCQP.n_ineq(prob) == 0 +@assert RCQP.n_comp(prob) == 0 +println(" PASSED") + +# --------------------------------------------------------------------------- +# Test 2 – SolverOptions getters and setters +# --------------------------------------------------------------------------- +println("Test 2: SolverOptions") + +opts = RCQP.SolverOptions() +RCQP.max_iters!(opts, 500) +RCQP.convergence_kkt_norm!(opts, 1e-6) +RCQP.output_dir!(opts, "/dev/null") +@assert RCQP.max_iters(opts) == 500 +@assert RCQP.convergence_kkt_norm(opts) ≈ 1e-6 +@assert RCQP.output_dir(opts) == "/dev/null" +println(" PASSED") + +# --------------------------------------------------------------------------- +# Test 3 – FilterEntry construction and accessors +# --------------------------------------------------------------------------- +println("Test 3: FilterEntry") + +entry = RCQP.FilterEntry(0.5, 1.2) +@assert RCQP.feas(entry) ≈ 0.5 +@assert RCQP.merit(entry) ≈ 1.2 +RCQP.feas!(entry, 0.1) +@assert RCQP.feas(entry) ≈ 0.1 +println(" PASSED") + +# --------------------------------------------------------------------------- +# Test 4 – Filter update and acceptability +# --------------------------------------------------------------------------- +println("Test 4: Filter") + +filt = RCQP.Filter() +@assert RCQP.size(filt) == 0 +RCQP.update(filt, RCQP.FilterEntry(1.0, 1.0)) +@assert RCQP.size(filt) == 1 +# entries() returns a flat [feas, merit, ...] vector +flat = RCQP.entries(filt) +@assert flat[1] ≈ 1.0 && flat[2] ≈ 1.0 +# A candidate strictly better on both axes should be acceptable +@assert RCQP.acceptable(filt, RCQP.FilterEntry(0.5, 0.5)) +RCQP.clear(filt) +@assert RCQP.size(filt) == 0 +println(" PASSED") + +# --------------------------------------------------------------------------- +# Test 5 – Retraction maps +# --------------------------------------------------------------------------- +println("Test 5: Retraction maps") + +solver = RCQP.Solver(opts) +RCQP.set_problem(solver, prob, ones(nz)) + +s = [0.1, 0.2] +r = RCQP.retract(solver, s, 0.1) +dr = RCQP.retract_deriv(solver, s, 0.1) +ddr = RCQP.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 = RCQP.solve(solver, opts) +@assert converged + +ws = RCQP.get_workspace(solver) +z = copy(RCQP.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(RCQP.kkt_residual(ws), Inf) < 1e-5 +@assert norm(RCQP.residual_eq(ws), Inf) < 1e-8 +@assert norm(RCQP.residual_ineq(ws), Inf) < 1e-8 +@assert norm(RCQP.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(RCQP.z_inds(solver)) == nz +@assert length(RCQP.m_eq_inds(solver)) == 0 # no equality constraints +println(" n_primals = $(RCQP.n_primals(solver))") +println(" n_duals = $(RCQP.n_duals(solver))") +println(" n_vars = $(RCQP.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) == RCQP.n_vars(solver) +println(" KKT size: $(size(kkt))") +println(" PASSED") + +# --------------------------------------------------------------------------- +println("\nAll tests passed!") diff --git a/tester.py b/tester.py new file mode 100644 index 0000000..67bceac --- /dev/null +++ b/tester.py @@ -0,0 +1,47 @@ +from pathlib import Path +import sys + +sys.path.insert(0, str(Path(__file__).resolve().parent / "build" / "python")) + +import rcqp + +print("module:", rcqp) +print("members:", dir(rcqp)) + +import numpy as np + +s = rcqp.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 = rcqp.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 = rcqp.SolverOptions() +opts.max_iters = 500 +opts.convergence_kkt_norm = 1e-6 +print(f"Options: {opts}") + +solver = rcqp.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/typings/rcqp.pyi b/typings/rcqp.pyi new file mode 100644 index 0000000..1cea473 --- /dev/null +++ b/typings/rcqp.pyi @@ -0,0 +1,224 @@ +""" +RCQP: constrained optimization solver with complementarity constraints +""" +from __future__ import annotations +import numpy +import scipy.sparse +import typing +__all__: list[str] = ['Filter', 'FilterEntry', 'Problem', 'Solver', 'SolverOptions', 'Workspace'] +class Filter: + @typing.overload + def __init__(self) -> None: + ... + @typing.overload + def __init__(self, gamma_objective: float, gamma_constraint: float) -> None: + ... + def acceptable(self, candidate: FilterEntry) -> bool: + ... + def candidate_acceptable(self, candidate: FilterEntry, entry: FilterEntry) -> bool: + ... + def candidate_dominated(self, candidate: FilterEntry, entry: FilterEntry) -> bool: + ... + def clear(self) -> None: + ... + def sufficient_progress(self, candidate: FilterEntry, entry: FilterEntry) -> tuple[bool, bool]: + ... + def update(self, new_entry: FilterEntry) -> None: + ... + @property + def entries(self) -> list[tuple[float, float]]: + ... + @property + def size(self) -> int: + ... +class FilterEntry: + feas: float + merit: float + def __init__(self) -> None: + ... + def __repr__(self) -> str: + ... +class Problem: + @typing.overload + def __init__(self, cost_hessian: numpy.ndarray, cost_gradient: numpy.ndarray, cost_const: float, J_eq: numpy.ndarray, c_eq: numpy.ndarray, J_ineq: numpy.ndarray, c_ineq: numpy.ndarray, J_comp: numpy.ndarray, c_comp: numpy.ndarray) -> None: + ... + @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, J_comp: scipy.sparse.csc_matrix, c_comp: numpy.ndarray) -> None: + ... + @property + def cost_const(self) -> float: + ... + @property + def cost_gradient(self) -> numpy.ndarray: + ... + @property + def n_comp(self) -> int: + ... + @property + def n_eq(self) -> int: + ... + @property + def n_ineq(self) -> int: + ... + @property + def nz(self) -> int: + ... +class Solver: + @typing.overload + def __init__(self) -> None: + ... + @typing.overload + def __init__(self, options: SolverOptions) -> None: + ... + def analytical_factorization(self) -> bool: + ... + def backsolve(self) -> None: + ... + def check_inertia(self) -> bool: + ... + def compute_amd_ordering(self) -> None: + ... + def convergence(self, options: SolverOptions) -> bool: + ... + def entry_from_solution(self, sqrt_relax_param: float, inv_penalty_param: float) -> tuple: + ... + def filter_linesearch(self, sqrt_relax_param: float, inv_penalty_param: float, max_iters: int) -> bool: + ... + def get_filter(self) -> Filter: + ... + def get_problem(self) -> Problem: + ... + def get_workspace(self) -> Workspace: + ... + def initialize_kkt_sparsity(self) -> None: + ... + def numerical_factorization(self) -> bool: + ... + def retract(self, s: numpy.ndarray, sqrt_relax_param: float) -> numpy.ndarray: + ... + def retract_deriv(self, s: numpy.ndarray, sqrt_relax_param: float) -> numpy.ndarray: + ... + def retract_second_deriv(self, s: numpy.ndarray, sqrt_relax_param: float) -> numpy.ndarray: + ... + def set_problem(self, problem: Problem, scaling: numpy.ndarray) -> None: + ... + def solve(self, options: SolverOptions) -> bool: + ... + def update_KKT_comp(self, s_comp: numpy.ndarray, m_comp: numpy.ndarray, sqrt_relax_param: float) -> None: + ... + def update_KKT_ineq(self, s_ineq: numpy.ndarray, sqrt_relax_param: float) -> None: + ... + def update_KKT_penalty(self, inv_penalty_param: float) -> None: + ... + def update_KKT_primal_regularizer(self, reg: float) -> None: + ... + def update_KKT_residual(self, sqrt_relax_param: float, inv_penalty_param: float) -> None: + ... + def update_KKT_system(self, sqrt_relax_param: float, inv_penalty_param: float) -> None: + ... + @property + def comp_L_inds(self) -> numpy.ndarray: + ... + @property + def comp_R_inds(self) -> numpy.ndarray: + ... + @property + def m_comp_inds(self) -> numpy.ndarray: + ... + @property + def m_eq_inds(self) -> numpy.ndarray: + ... + @property + def m_ineq_inds(self) -> numpy.ndarray: + ... + @property + def n_duals(self) -> int: + ... + @property + def n_primals(self) -> int: + ... + @property + def n_vars(self) -> int: + ... + @property + def s_comp_inds(self) -> numpy.ndarray: + ... + @property + def s_ineq_inds(self) -> numpy.ndarray: + ... + @property + def z_inds(self) -> numpy.ndarray: + ... +class SolverOptions: + convergence_comp_violation: float + convergence_eq_violation: float + convergence_ineq_violation: float + convergence_kkt_norm: float + gamma_constraint: float + gamma_objective: float + max_iters: int + max_iters_linesearch: int + outer_step_kkt_norm: float + output_dir: str + penalty_initial: float + penalty_max: float + penalty_scaling: float + relaxation_initial: float + relaxation_min: float + relaxation_scaling: float + def __init__(self) -> None: + ... + def __repr__(self) -> str: + ... +class Workspace: + @property + def kkt_residual(self) -> numpy.ndarray: + ... + @property + def m_comp(self) -> numpy.ndarray: + ... + @property + def m_comp_est(self) -> numpy.ndarray: + ... + @property + def m_eq(self) -> numpy.ndarray: + ... + @property + def m_eq_est(self) -> numpy.ndarray: + ... + @property + def m_ineq(self) -> numpy.ndarray: + ... + @property + def m_ineq_est(self) -> numpy.ndarray: + ... + @property + def newton_step(self) -> numpy.ndarray: + ... + @property + def penalty_param(self) -> float: + ... + @property + def relax_param(self) -> float: + ... + @property + def residual_comp(self) -> numpy.ndarray: + ... + @property + def residual_eq(self) -> numpy.ndarray: + ... + @property + def residual_ineq(self) -> numpy.ndarray: + ... + @property + def s_comp(self) -> numpy.ndarray: + ... + @property + def s_ineq(self) -> numpy.ndarray: + ... + @property + def solution(self) -> numpy.ndarray: + ... + @property + def z(self) -> numpy.ndarray: + ... diff --git a/wrapper.cpp b/wrapper.cpp deleted file mode 100644 index 9477c49..0000000 --- a/wrapper.cpp +++ /dev/null @@ -1,216 +0,0 @@ -#include -#include "solver.h" - -using jl_VecXd = jlcxx::ArrayRef; -using jl_MatXd = jlcxx::ArrayRef; -using jl_VecXi = jlcxx::ArrayRef; -using jl_MatXi = jlcxx::ArrayRef; - -// Helper functions for converting between Eigen types and Julia types -template -Eigen::Matrix to_eigen(jlcxx::ArrayRef& arr, int rows, int cols) { - if (rows == 0 || cols == 0) - return Eigen::Matrix(rows, cols); - return Eigen::Matrix::Map(arr.data(), rows, cols); -} - -template -Eigen::Matrix to_eigen(jlcxx::ArrayRef& arr) { - if (arr.size() == 0) - return Eigen::Matrix(0); - return Eigen::Matrix::Map(arr.data(), arr.size()); -} - -template -jlcxx::ArrayRef to_julia(Eigen::Matrix& vec) { - return jlcxx::make_julia_array(vec.data(), vec.size()); -} - -template -jlcxx::ArrayRef to_julia(Eigen::Map>& vec) { - return jlcxx::make_julia_array(vec.data(), vec.size()); -} - -template -jlcxx::ArrayRef to_julia(Eigen::Matrix& mat) { - return jlcxx::make_julia_array(mat.data(), mat.rows(), mat.cols()); -} - -Eigen::MatrixXd eigen_mat = Eigen::MatrixXd::Random(3, 3); -Eigen::VectorXd eigen_vec = Eigen::VectorXd::Random(5); - -JLCXX_MODULE define_julia_module(jlcxx::Module& mod) { - // Problem class bindings - mod.add_type("Problem") - .constructor([](jl_MatXd cost_hessian, jl_VecXd cost_gradient, double cost_const, - jl_MatXd J_eq, jl_VecXd c_eq, - jl_MatXd J_ineq, jl_VecXd c_ineq, - jl_MatXd J_comp, jl_VecXd c_comp) { - int nz = cost_gradient.size(); - int n_eq = c_eq.size(); - int n_ineq = c_ineq.size(); - int n_comp = c_comp.size(); - - return new Problem( - to_eigen(cost_hessian, nz, nz), to_eigen(cost_gradient), cost_const, - to_eigen(J_eq, n_eq, nz), to_eigen(c_eq), - to_eigen(J_ineq, n_ineq, nz), to_eigen(c_ineq), - to_eigen(J_comp, n_comp, nz), to_eigen(c_comp) - ); - }) - .method("nz", [](Problem& p) -> const int { return p.nz; }) - .method("n_eq", [](Problem& p) -> const int { return p.n_eq; }) - .method("n_ineq", [](Problem& p) -> const int { return p.n_ineq; }) - .method("n_comp", [](Problem& p) -> const int { return p.n_comp; }); - - #define OPTION_SETTER(name, type) \ - .method(#name "!", [](Solver::Options& o, type v) { o.name = v; }) \ - .method(#name, [](Solver::Options& o) { return o.name; }) - - mod.add_type("SolverOptions") - .constructor() - OPTION_SETTER(convergence_kkt_norm, double) - OPTION_SETTER(convergence_eq_violation, double) - OPTION_SETTER(convergence_ineq_violation, double) - OPTION_SETTER(convergence_comp_violation, double) - OPTION_SETTER(outer_step_kkt_norm, double) - OPTION_SETTER(penalty_initial, double) - OPTION_SETTER(penalty_max, double) - OPTION_SETTER(penalty_scaling, double) - OPTION_SETTER(relaxation_initial, double) - OPTION_SETTER(relaxation_min, double) - OPTION_SETTER(relaxation_scaling, double) - OPTION_SETTER(max_iters, int) - OPTION_SETTER(max_iters_linesearch, int) - OPTION_SETTER(gamma_objective, double) - OPTION_SETTER(gamma_constraint, double) - .method("output_dir!", [](Solver::Options& o, const std::string& v) { o.output_dir = v; }) - .method("output_dir", [](Solver::Options& o) { return o.output_dir.string(); }); - - #undef OPTION_SETTER - - // Filter class bindings - mod.add_type("Filter") - .constructor() - .method("filter_clear", &Filter::clear) - .method("filter_size", []( Filter& f ) { return f.entries.size(); }) - .method("filter_entries", [](Filter& f) { - std::vector entries; - entries.reserve(f.entries.size() * 2); - for (const auto& entry : f.entries) { - entries.push_back(entry.feas); - entries.push_back(entry.merit); - } - return entries; - }); - - // Workspace class bindings - mod.add_type("Workspace") - .constructor([](const Problem& prob) { return new Workspace(prob); }) - .method("solution", [](Workspace& w) { return to_julia(w.solution); }) - .method("z", [](Workspace& w) { return to_julia(w.z); }) - .method("s_ineq", [](Workspace& w) { return to_julia(w.s_ineq); }) - .method("s_comp", [](Workspace& w) { return to_julia(w.s_comp); }) - .method("m_eq", [](Workspace& w) { return to_julia(w.m_eq); }) - .method("m_ineq", [](Workspace& w) { return to_julia(w.m_ineq); }) - .method("m_comp", [](Workspace& w) { return to_julia(w.m_comp); }) - .method("m_eq_est", [](Workspace& w) { return to_julia(w.m_eq_est); }) - .method("m_ineq_est", [](Workspace& w) { return to_julia(w.m_ineq_est); }) - .method("m_comp_est", [](Workspace& w) { return to_julia(w.m_comp_est); }) - .method("kkt_residual", [](Workspace& w) { return to_julia(w.kkt_residual); }) - .method("residual_eq", [](Workspace& w) { return to_julia(w.residual_eq); }) - .method("residual_ineq", [](Workspace& w) { return to_julia(w.residual_ineq); }) - .method("residual_comp", [](Workspace& w) { return to_julia(w.residual_comp); }) - .method("relax_param", [](Workspace& w) { return w.relax_param; }) - .method("penalty_param", [](Workspace& w) { return w.penalty_param; }) - .method("kkt_system", [](Workspace& w) { - SMat& kkt = w.kkt_system; - kkt.makeCompressed(); - jlcxx::ArrayRef colptr = jlcxx::make_julia_array(kkt.outerIndexPtr(), kkt.outerSize() + 1); - jlcxx::ArrayRef rowval = jlcxx::make_julia_array(kkt.innerIndexPtr(), kkt.nonZeros()); - jl_VecXd nzval = jlcxx::make_julia_array(kkt.valuePtr(), kkt.nonZeros()); - - return std::make_tuple(kkt.rows(), kkt.cols(), colptr, rowval, nzval); - }) - .method("newton_step", [](Workspace& w) { return to_julia(w.newton_step); }) - .method("D", [](Workspace& w) { return to_julia(w.D); }) - .method("amd_perm_vec", [](Workspace& w) { return to_julia(w.amd_perm_vec); }) - .method("amd_iperm_vec", [](Workspace& w) { return to_julia(w.amd_iperm_vec); }); - // .method("amd_perm", [](Workspace& w) { return to_julia(w.amd_perm.toDense()); }) -; - // Solver class bindings - mod.add_type("Solver") - .constructor() - .method("retract", [](Solver& solver, jl_VecXd s, double sqrt_relax_param) { - Vec p = solver.retract(to_eigen(s), sqrt_relax_param); - return std::vector(p.data(), p.data() + p.size()); - }) - .method("retract_deriv", [](Solver& solver, jl_VecXd s, double sqrt_relax_param) { - Vec p_deriv = solver.retract_deriv(to_eigen(s), sqrt_relax_param); - return std::vector(p_deriv.data(), p_deriv.data() + p_deriv.size()); - }) - .method("retract_second_deriv", [](Solver& solver, jl_VecXd s, double sqrt_relax_param) { - Vec p_second_deriv = solver.retract_second_deriv(to_eigen(s), sqrt_relax_param); - return std::vector(p_second_deriv.data(), p_second_deriv.data() + p_second_deriv.size()); - }) - .method("update_KKT_residual", &Solver::update_KKT_residual) - .method("update_KKT_system", &Solver::update_KKT_system) - .method("update_KKT_ineq", [](Solver& solver, jl_VecXd s_ineq, double sqrt_relax_param) { - solver.update_KKT_ineq(to_eigen(s_ineq), sqrt_relax_param); - }) - .method("update_KKT_comp", [](Solver& solver, jl_VecXd s_comp, jl_VecXd m_comp, double sqrt_relax_param) { - solver.update_KKT_comp(to_eigen(s_comp), to_eigen(m_comp), sqrt_relax_param); - }) - .method("update_KKT_penalty", [](Solver& solver, double inv_penalty_param) { - solver.update_KKT_penalty(inv_penalty_param); - }) - .method("update_KKT_primal_regularizer", &Solver::update_KKT_primal_regularizer) - .method("analytical_factorization", &Solver::analytical_factorization) - .method("numerical_factorization", &Solver::numerical_factorization) - .method("backsolve", &Solver::backsolve) - .method("check_inertia", &Solver::check_inertia) - .method("compute_amd_ordering", &Solver::compute_amd_ordering) - .method("solve", &Solver::solve) - .method("set_problem", [](Solver& solver, Problem& prob, jl_VecXd scaling) { - solver.set_problem(prob, to_eigen(scaling)); - }) - .method("get_problem", &Solver::get_problem) - .method("z_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.z_inds); }) - .method("s_ineq_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.s_ineq_inds); }) - .method("s_comp_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.s_comp_inds); }) - .method("m_eq_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.m_eq_inds); }) - .method("m_ineq_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.m_ineq_inds); }) - .method("m_comp_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.m_comp_inds); }) - .method("comp_L_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.comp_L_inds); }) - .method("comp_R_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.comp_R_inds); }) - .method("z_z_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.z_z_inds); }) - .method("s_ineq_s_ineq_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.s_ineq_s_ineq_inds); }) - .method("s_ineq_m_ineq_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.s_ineq_m_ineq_inds); }) - .method("s_comp_s_comp_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.s_comp_s_comp_inds); }) - .method("s_comp_m_comp_inds", [](Solver& s) -> jl_VecXi { return to_julia(s.s_comp_m_comp_inds); }) - .method("filter_linesearch", &Solver::filter_linesearch) - .method("entry_from_solution", [](Solver& solver, double sqrt_relax_param, double inv_penalty_param) { - Filter::Entry entry = solver.entry_from_solution(sqrt_relax_param, inv_penalty_param); - return std::make_tuple(entry.feas, entry.merit); - }) - .method("convergence", &Solver::convergence) - .method("get_workspace", &Solver::get_workspace) - .method("get_filter", &Solver::get_filter); - - mod.method("test_array", [](jlcxx::ArrayRef arr) { - std::cout << "here" << std::endl; - eigen_vec = to_eigen(arr); - std::cout << "Received array of size: " << eigen_vec.size() << std::endl; - std::cout << "Contents: " << eigen_vec.transpose() << std::endl; - eigen_vec[2] += 5; - return to_julia(eigen_vec); - }); - mod.method("test_matrix", [](jlcxx::ArrayRef arr, int rows, int cols) { - std::cout << "here" << std::endl; - eigen_mat = to_eigen(arr, rows, cols); - std::cout << "Received matrix of size: " << eigen_mat.rows() << "x" << eigen_mat.cols() << std::endl; - std::cout << "Contents:\n" << eigen_mat << std::endl; - eigen_mat(0, 0) += 5; - return to_julia(eigen_mat); - }); -} \ No newline at end of file From 6d0add5f75a6a044d8889624ccb5160008489133 Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Sat, 4 Apr 2026 01:39:23 -0400 Subject: [PATCH 2/9] update stuff to be able to install rcqp with pip in vevnv --- README.md | 92 ++++++++++++++------------------------------------ pyproject.toml | 19 +++++++++++ tester.py | 5 --- 3 files changed, 45 insertions(+), 71 deletions(-) create mode 100644 pyproject.toml diff --git a/README.md b/README.md index 5791026..2e029c5 100644 --- a/README.md +++ b/README.md @@ -33,7 +33,7 @@ Output: `build/python/rcqp.cpython--.so` cmake --preset julia cmake --build --preset julia ``` -CMake auto-detects JlCxx by calling `julia -e 'using CxxWrap; print(CxxWrap.prefix_path())'` — no manual paths required. +CMake auto-detects JlCxx by calling `julia -e 'using CxxWrap; print(CxxWrap.prefix_path())'` Output: `build/lib/librcqp_julia.dylib` (macOS) / `build/lib/librcqp_julia.so` (Linux) @@ -43,58 +43,46 @@ cmake --preset all cmake --build --preset all ``` - - -## Python — pip install (editable) - -For a proper install into your virtual environment (no `sys.path` hacks): + +**3. Install rcqp:** ```bash -pip install scikit-build-core pybind11-stubgen pip install -e . --no-build-isolation ``` After this, `import rcqp` works from anywhere in that environment. -## Python — usage without installing +## Python: usage without installing Add `build/python/` to your path at runtime: ```python import sys from pathlib import Path -sys.path.insert(0, str(Path(__file__).resolve().parent / "build" / "python")) -import rcqp -import numpy as np - -opts = rcqp.SolverOptions() -opts.max_iters = 500 - -solver = rcqp.Solver(opts) -solver.set_problem(prob, scaling) -converged = solver.solve(opts) +# Insert the path to the `build` directory in PYTHONPATH +path_to_rcqp = ... +sys.path.insert(0, str(path_to_rcqp / "build" / "python")) -ws = solver.get_workspace() -print(ws.z) # primal solution +import rcqp ``` -## VSCode — Python autocomplete +## VSCode: Python autocomplete Autocomplete and type-checking are driven by the `.pyi` stubs in `typings/`, which are regenerated automatically each time you build `rcqp_python` (requires `pybind11-stubgen`). @@ -109,8 +97,8 @@ Two config files in this repo wire everything up for VSCode automatically: } ``` -- `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 compiled `.so` lives so it can be imported +- `stubPath`: tells Pylance where to find the `.pyi` stub file for type info ### Install the Pylance extension Install the **Pylance** extension in VSCode (`ms-python.vscode-pylance`). It picks up the above settings automatically. @@ -120,35 +108,7 @@ If stubs are stale or missing, rebuild: cmake --build --preset python --target rcqp_python ``` -If `pybind11-stubgen` is not installed, CMake will warn but the build still succeeds — you just won't get updated stubs. +If `pybind11-stubgen` is not installed, CMake will warn but the build still succeeds. Ensure `pybind11-stubgen` is installed with: ```bash pip install pybind11-stubgen -``` - -## Julia — usage - -```julia -using CxxWrap - -module RCQP - using CxxWrap - @wrapmodule(() -> joinpath(@__DIR__, "build/lib/librcqp_julia")) - function __init__() - @initcxx - end -end - -prob = RCQP.Problem(H, g, 0.0, J_eq, c_eq, J_ineq, c_ineq, J_comp, c_comp) - -opts = RCQP.SolverOptions() -RCQP.max_iters!(opts, 500) - -solver = RCQP.Solver(opts) -RCQP.set_problem(solver, prob, ones(nz)) -converged = RCQP.solve(solver, opts) - -ws = RCQP.get_workspace(solver) -z = copy(RCQP.z(ws)) -``` - -See `tester.jl` for a complete working example. +``` \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..bc31d28 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,19 @@ +[build-system] +requires = ["scikit-build-core>=0.9", "pybind11>=2.13"] +build-backend = "scikit_build_core.build" + +[project] +name = "rcqp" +version = "0.1.0" +description = "RCQP: constrained optimization solver with complementarity constraints" +requires-python = ">=3.8" +dependencies = ["numpy"] + +[tool.scikit-build] +cmake.build-type = "Release" +cmake.args = [ + "-DRCQP_BUILD_PYTHON=ON", + "-DRCQP_BUILD_JULIA=OFF", + "-DRCQP_GENERATE_PYI=OFF", +] +wheel.packages = [] diff --git a/tester.py b/tester.py index 67bceac..6c003eb 100644 --- a/tester.py +++ b/tester.py @@ -1,8 +1,3 @@ -from pathlib import Path -import sys - -sys.path.insert(0, str(Path(__file__).resolve().parent / "build" / "python")) - import rcqp print("module:", rcqp) From b8b55307db43afc9e49dbf18e2dc071da61a8790 Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Sat, 4 Apr 2026 01:50:14 -0400 Subject: [PATCH 3/9] add docker container --- .dockerignore | 5 ++++ Dockerfile | 70 +++++++++++++++++++++++++++++++++++++++++++++++++++ README.md | 13 ++++++++++ 3 files changed, 88 insertions(+) create mode 100644 .dockerignore create mode 100644 Dockerfile diff --git a/.dockerignore b/.dockerignore new file mode 100644 index 0000000..cbddda0 --- /dev/null +++ b/.dockerignore @@ -0,0 +1,5 @@ +build/ +.venv/ +__pycache__/ +*.pyc +.git/ diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..48de347 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,70 @@ +FROM julia:1.10 + +ENV DEBIAN_FRONTEND=noninteractive + +# Compiler toolchain + C++ headers that the RCQP build needs at compile time: +# git — required by CMake FetchContent to clone qdldl at configure time +# libeigen3-dev — linear algebra (matrix/vector types) +# nlohmann-json3-dev — JSON (used for solver config/output) +# python3-venv — needed to create an isolated Python environment below +RUN apt-get update && apt-get install -y \ + build-essential \ + git \ + python3 \ + python3-venv \ + python3-dev \ + libeigen3-dev \ + nlohmann-json3-dev \ + && rm -rf /var/lib/apt/lists/* + +# CxxWrap.jl is the Julia package that lets Julia call into a compiled C++ library. +# It provides the JlCxx CMake target that julia_bindings.cpp links against. +RUN julia -e 'using Pkg; Pkg.add("CxxWrap")' + +# Create an isolated Python virtual environment at /opt/venv. +# This sidesteps the "externally-managed-environment" restriction on the +# system Python, which blocks pip from installing packages system-wide. +ENV VIRTUAL_ENV=/opt/venv +RUN python3 -m venv $VIRTUAL_ENV + +# Prepend the venv to PATH so every subsequent `python` / `pip` / `cmake` +# command in this file (and in the running container) uses the venv. +ENV PATH="$VIRTUAL_ENV/bin:$PATH" + +# Install cmake and pybind11 into the venv. +# cmake — apt only provides ~3.25; this project requires 3.28+, so we get +# a newer version from PyPI (pip's cmake package bundles the binary) +# pybind11 — the C++/Python bridge; CMake needs its cmake config files to +# find and link pybind11 when building the Python extension module +# numpy — runtime dependency for working with the solver in Python +RUN pip install cmake pybind11 numpy + +WORKDIR /rcqp +COPY . . + +# Build both the Python and Julia bindings in one pass using the "all" preset +# defined in CMakePresets.json. +# +# cmake --preset all +# "configure" step: reads CMakeLists.txt, finds all dependencies, and +# generates the actual build files (Makefiles / Ninja files) in build/. +# Equivalent to: +# cmake -B build \ +# -DCMAKE_BUILD_TYPE=Release \ +# -DRCQP_BUILD_PYTHON=ON \ +# -DRCQP_BUILD_JULIA=ON \ +# -DRCQP_GENERATE_PYI=ON +# +# cmake --build --preset all +# "build" step: runs the compiler using the files generated above. +# Produces: +# build/python/rcqp..so — Python extension module +# build/lib/librcqp_julia.so — Julia shared library +# +RUN cmake --preset all && cmake --build --preset all + +# Add the directory containing the compiled Python extension to PYTHONPATH +# so that `import rcqp` works inside the venv without a separate `pip install`. +ENV PYTHONPATH="/rcqp/build/python" + +CMD ["/bin/bash"] diff --git a/README.md b/README.md index 2e029c5..5d0b46e 100644 --- a/README.md +++ b/README.md @@ -43,6 +43,19 @@ cmake --preset all cmake --build --preset all ``` +## Docker + +The included `Dockerfile` builds both Python and Julia bindings inside a self-contained image. + +```bash +docker build -t rcqp . +docker run -it rcqp +``` + +Inside the container: +- `python` uses the venv at `/opt/venv`; `import rcqp` works immediately +- Julia shared library is at `build/lib/librcqp_julia.so` + ## Python: pip install (editable) For a proper install into your virtual environment: From c6f65c9fdbb0990327976358c467454d940f99d8 Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Sat, 4 Apr 2026 01:51:26 -0400 Subject: [PATCH 4/9] add bindings/ dir --- CMakeLists.txt | 4 +- .../julia_bindings.cpp | 0 .../python_bindings.cpp | 0 typings/rcqp.pyi | 198 +++++++++++++----- 4 files changed, 144 insertions(+), 58 deletions(-) rename julia_bindings.cpp => bindings/julia_bindings.cpp (100%) rename python_bindings.cpp => bindings/python_bindings.cpp (100%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5352fa1..88160b1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -91,7 +91,7 @@ if(RCQP_BUILD_PYTHON) FetchContent_MakeAvailable(pybind11) endif() - pybind11_add_module(rcqp_python MODULE python_bindings.cpp) + pybind11_add_module(rcqp_python MODULE bindings/python_bindings.cpp) target_link_libraries(rcqp_python PRIVATE rcqp) set_target_properties(rcqp_python PROPERTIES OUTPUT_NAME rcqp @@ -158,7 +158,7 @@ if(RCQP_BUILD_JULIA) endif() find_package(JlCxx REQUIRED) - add_library(rcqp_julia SHARED julia_bindings.cpp) + add_library(rcqp_julia SHARED bindings/julia_bindings.cpp) target_link_libraries(rcqp_julia PRIVATE rcqp JlCxx::cxxwrap_julia) set_target_properties(rcqp_julia PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib diff --git a/julia_bindings.cpp b/bindings/julia_bindings.cpp similarity index 100% rename from julia_bindings.cpp rename to bindings/julia_bindings.cpp diff --git a/python_bindings.cpp b/bindings/python_bindings.cpp similarity index 100% rename from python_bindings.cpp rename to bindings/python_bindings.cpp diff --git a/typings/rcqp.pyi b/typings/rcqp.pyi index 1cea473..532a0ec 100644 --- a/typings/rcqp.pyi +++ b/typings/rcqp.pyi @@ -3,6 +3,7 @@ RCQP: constrained optimization solver with complementarity constraints """ from __future__ import annotations import numpy +import numpy.typing import scipy.sparse import typing __all__: list[str] = ['Filter', 'FilterEntry', 'Problem', 'Solver', 'SolverOptions', 'Workspace'] @@ -11,7 +12,7 @@ class Filter: def __init__(self) -> None: ... @typing.overload - def __init__(self, gamma_objective: float, gamma_constraint: float) -> None: + def __init__(self, gamma_objective: typing.SupportsFloat | typing.SupportsIndex, gamma_constraint: typing.SupportsFloat | typing.SupportsIndex) -> None: ... def acceptable(self, candidate: FilterEntry) -> bool: ... @@ -32,24 +33,34 @@ class Filter: def size(self) -> int: ... class FilterEntry: - feas: float - merit: float def __init__(self) -> None: ... def __repr__(self) -> str: ... + @property + def feas(self) -> float: + ... + @feas.setter + def feas(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def merit(self) -> float: + ... + @merit.setter + def merit(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... class Problem: @typing.overload - def __init__(self, cost_hessian: numpy.ndarray, cost_gradient: numpy.ndarray, cost_const: float, J_eq: numpy.ndarray, c_eq: numpy.ndarray, J_ineq: numpy.ndarray, c_ineq: numpy.ndarray, J_comp: numpy.ndarray, c_comp: numpy.ndarray) -> None: + def __init__(self, cost_hessian: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, n]"], cost_gradient: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], cost_const: typing.SupportsFloat | typing.SupportsIndex, J_eq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, n]"], c_eq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], J_ineq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, n]"], c_ineq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], J_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, n]"], c_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"]) -> None: ... @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, J_comp: scipy.sparse.csc_matrix, c_comp: numpy.ndarray) -> None: + def __init__(self, cost_hessian: scipy.sparse.csc_matrix, cost_gradient: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], cost_const: typing.SupportsFloat | typing.SupportsIndex, J_eq: scipy.sparse.csc_matrix, c_eq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], J_ineq: scipy.sparse.csc_matrix, c_ineq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], J_comp: scipy.sparse.csc_matrix, c_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"]) -> None: ... @property def cost_const(self) -> float: ... @property - def cost_gradient(self) -> numpy.ndarray: + def cost_gradient(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property def n_comp(self) -> int: @@ -80,9 +91,9 @@ class Solver: ... def convergence(self, options: SolverOptions) -> bool: ... - def entry_from_solution(self, sqrt_relax_param: float, inv_penalty_param: float) -> tuple: + def entry_from_solution(self, sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex) -> tuple[float, float]: ... - def filter_linesearch(self, sqrt_relax_param: float, inv_penalty_param: float, max_iters: int) -> bool: + def filter_linesearch(self, sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex, max_iters: typing.SupportsInt | typing.SupportsIndex) -> bool: ... def get_filter(self) -> Filter: ... @@ -94,42 +105,42 @@ class Solver: ... def numerical_factorization(self) -> bool: ... - def retract(self, s: numpy.ndarray, sqrt_relax_param: float) -> numpy.ndarray: + def retract(self, s: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... - def retract_deriv(self, s: numpy.ndarray, sqrt_relax_param: float) -> numpy.ndarray: + def retract_deriv(self, s: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... - def retract_second_deriv(self, s: numpy.ndarray, sqrt_relax_param: float) -> numpy.ndarray: + def retract_second_deriv(self, s: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... - def set_problem(self, problem: Problem, scaling: numpy.ndarray) -> None: + def set_problem(self, problem: Problem, scaling: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"]) -> None: ... def solve(self, options: SolverOptions) -> bool: ... - def update_KKT_comp(self, s_comp: numpy.ndarray, m_comp: numpy.ndarray, sqrt_relax_param: float) -> None: + def update_KKT_comp(self, s_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], m_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> None: ... - def update_KKT_ineq(self, s_ineq: numpy.ndarray, sqrt_relax_param: float) -> None: + def update_KKT_ineq(self, s_ineq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> None: ... - def update_KKT_penalty(self, inv_penalty_param: float) -> None: + def update_KKT_penalty(self, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex) -> None: ... - def update_KKT_primal_regularizer(self, reg: float) -> None: + def update_KKT_primal_regularizer(self, reg: typing.SupportsFloat | typing.SupportsIndex) -> None: ... - def update_KKT_residual(self, sqrt_relax_param: float, inv_penalty_param: float) -> None: + def update_KKT_residual(self, sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex) -> None: ... - def update_KKT_system(self, sqrt_relax_param: float, inv_penalty_param: float) -> None: + def update_KKT_system(self, sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex) -> None: ... @property - def comp_L_inds(self) -> numpy.ndarray: + def comp_L_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: ... @property - def comp_R_inds(self) -> numpy.ndarray: + def comp_R_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: ... @property - def m_comp_inds(self) -> numpy.ndarray: + def m_comp_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: ... @property - def m_eq_inds(self) -> numpy.ndarray: + def m_eq_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: ... @property - def m_ineq_inds(self) -> numpy.ndarray: + def m_ineq_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: ... @property def n_duals(self) -> int: @@ -141,59 +152,134 @@ class Solver: def n_vars(self) -> int: ... @property - def s_comp_inds(self) -> numpy.ndarray: + def s_comp_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: ... @property - def s_ineq_inds(self) -> numpy.ndarray: + def s_ineq_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: ... @property - def z_inds(self) -> numpy.ndarray: + def z_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: ... class SolverOptions: - convergence_comp_violation: float - convergence_eq_violation: float - convergence_ineq_violation: float - convergence_kkt_norm: float - gamma_constraint: float - gamma_objective: float - max_iters: int - max_iters_linesearch: int - outer_step_kkt_norm: float output_dir: str - penalty_initial: float - penalty_max: float - penalty_scaling: float - relaxation_initial: float - relaxation_min: float - relaxation_scaling: float def __init__(self) -> None: ... def __repr__(self) -> str: ... + @property + def convergence_comp_violation(self) -> float: + ... + @convergence_comp_violation.setter + def convergence_comp_violation(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def convergence_eq_violation(self) -> float: + ... + @convergence_eq_violation.setter + def convergence_eq_violation(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def convergence_ineq_violation(self) -> float: + ... + @convergence_ineq_violation.setter + def convergence_ineq_violation(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def convergence_kkt_norm(self) -> float: + ... + @convergence_kkt_norm.setter + def convergence_kkt_norm(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def gamma_constraint(self) -> float: + ... + @gamma_constraint.setter + def gamma_constraint(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def gamma_objective(self) -> float: + ... + @gamma_objective.setter + def gamma_objective(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def max_iters(self) -> int: + ... + @max_iters.setter + def max_iters(self, arg0: typing.SupportsInt | typing.SupportsIndex) -> None: + ... + @property + def max_iters_linesearch(self) -> int: + ... + @max_iters_linesearch.setter + def max_iters_linesearch(self, arg0: typing.SupportsInt | typing.SupportsIndex) -> None: + ... + @property + def outer_step_kkt_norm(self) -> float: + ... + @outer_step_kkt_norm.setter + def outer_step_kkt_norm(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def penalty_initial(self) -> float: + ... + @penalty_initial.setter + def penalty_initial(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def penalty_max(self) -> float: + ... + @penalty_max.setter + def penalty_max(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def penalty_scaling(self) -> float: + ... + @penalty_scaling.setter + def penalty_scaling(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def relaxation_initial(self) -> float: + ... + @relaxation_initial.setter + def relaxation_initial(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def relaxation_min(self) -> float: + ... + @relaxation_min.setter + def relaxation_min(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... + @property + def relaxation_scaling(self) -> float: + ... + @relaxation_scaling.setter + def relaxation_scaling(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + ... class Workspace: @property - def kkt_residual(self) -> numpy.ndarray: + def kkt_residual(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def m_comp(self) -> numpy.ndarray: + def m_comp(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def m_comp_est(self) -> numpy.ndarray: + def m_comp_est(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def m_eq(self) -> numpy.ndarray: + def m_eq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def m_eq_est(self) -> numpy.ndarray: + def m_eq_est(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def m_ineq(self) -> numpy.ndarray: + def m_ineq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def m_ineq_est(self) -> numpy.ndarray: + def m_ineq_est(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def newton_step(self) -> numpy.ndarray: + def newton_step(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property def penalty_param(self) -> float: @@ -202,23 +288,23 @@ class Workspace: def relax_param(self) -> float: ... @property - def residual_comp(self) -> numpy.ndarray: + def residual_comp(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def residual_eq(self) -> numpy.ndarray: + def residual_eq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def residual_ineq(self) -> numpy.ndarray: + def residual_ineq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def s_comp(self) -> numpy.ndarray: + def s_comp(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def s_ineq(self) -> numpy.ndarray: + def s_ineq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def solution(self) -> numpy.ndarray: + def solution(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... @property - def z(self) -> numpy.ndarray: + def z(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: ... From 454050d248b593acf0a6cc2f76b155331a790d64 Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Sat, 4 Apr 2026 01:53:44 -0400 Subject: [PATCH 5/9] update readme --- README.md | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 5d0b46e..0f53c36 100644 --- a/README.md +++ b/README.md @@ -3,14 +3,14 @@ A C++ solver for quadratic programs with linear complementarity constraints, wit ## Prerequisites -| Dependency | Required for | Install | -|---|---|---| -| CMake ≥ 3.28 | all | [cmake.org](https://cmake.org/download/) or `brew install cmake` | -| Eigen3 | all | `brew install eigen` / `apt install libeigen3-dev` | -| nlohmann-json | all | `brew install nlohmann-json` / `apt install nlohmann-json3-dev` | -| Python ≥ 3.8 + dev headers | Python bindings | [python.org](https://www.python.org/) | -| pybind11-stubgen | `.pyi` autocomplete stubs | `pip install pybind11-stubgen` | -| Julia ≥ 1.9 + CxxWrap.jl | Julia bindings | see below | +| Dependency | Required for | +|---|---| +| CMake ≥ 3.28 | all | +| Eigen3 | all | +| nlohmann-json | all | +| Python ≥ 3.8 + dev headers | Python bindings | +| pybind11-stubgen | `.pyi` autocomplete stubs | +| Julia ≥ 1.9 + CxxWrap.jl | Julia bindings | **Julia setup** (only needed for Julia bindings): ```bash From 8a55791d78a8c2d063d366afb14e62cac4d2d51b Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Mon, 6 Apr 2026 23:47:44 -0400 Subject: [PATCH 6/9] update python bindings --- bindings/python_bindings.cpp | 48 +++++++- tester.py | 5 +- typings/rcqp.pyi | 205 +++++++++++++---------------------- 3 files changed, 120 insertions(+), 138 deletions(-) diff --git a/bindings/python_bindings.cpp b/bindings/python_bindings.cpp index c195d1c..c6acbfe 100644 --- a/bindings/python_bindings.cpp +++ b/bindings/python_bindings.cpp @@ -67,6 +67,7 @@ PYBIND11_MODULE(rcqp, m) { .def_readwrite("max_iters_linesearch", &Solver::Options::max_iters_linesearch) .def_readwrite("gamma_objective", &Solver::Options::gamma_objective) .def_readwrite("gamma_constraint", &Solver::Options::gamma_constraint) + .def_readwrite("ruiz_iterations", &Solver::Options::ruiz_iterations) .def_property("output_dir", [](const Solver::Options& o) { return o.output_dir.string(); }, [](Solver::Options& o, const std::string& v) { o.output_dir = v; }) @@ -156,7 +157,29 @@ PYBIND11_MODULE(rcqp, m) { [](const Workspace& w) { return w.penalty_param; }) // Newton step .def_property_readonly("newton_step", - [](const Workspace& w) { return w.newton_step; }); + [](const Workspace& w) { return w.newton_step; }) + // KKT system: returns (rows, cols, colptr, rowval, nzval) as a tuple + .def_property_readonly("kkt_system", + [](Workspace& w) { + SMat& kkt = w.kkt_system; + kkt.makeCompressed(); + Eigen::VectorXi colptr(kkt.outerSize() + 1); + Eigen::VectorXi rowval(kkt.nonZeros()); + Vec nzval(kkt.nonZeros()); + std::copy(kkt.outerIndexPtr(), kkt.outerIndexPtr() + kkt.outerSize() + 1, colptr.data()); + std::copy(kkt.innerIndexPtr(), kkt.innerIndexPtr() + kkt.nonZeros(), rowval.data()); + std::copy(kkt.valuePtr(), kkt.valuePtr() + kkt.nonZeros(), nzval.data()); + return py::make_tuple((int)kkt.rows(), (int)kkt.cols(), colptr, rowval, nzval); + }) + // QDLDL factorization data + .def_property_readonly("D", + [](const Workspace& w) { return w.D; }) + .def_property_readonly("amd_perm_vec", + [](const Workspace& w) { return w.amd_perm_vec; }) + .def_property_readonly("amd_iperm_vec", + [](const Workspace& w) { return w.amd_iperm_vec; }) + .def_property_readonly("scaling", + [](const Workspace& w) { return w.scaling; }); // ------------------------------------------------------------------------- // Solver @@ -165,10 +188,15 @@ PYBIND11_MODULE(rcqp, m) { .def(py::init<>()) .def(py::init(), py::arg("options")) // Problem setup - .def("set_problem", [](Solver& s, Problem& prob, const Vec& scaling) { - s.set_problem(prob, scaling); + .def("set_problem", [](Solver& s, Problem& prob) { + s.set_problem(prob); + }, + py::arg("problem")) + .def("ruiz_equilibration", [](Solver& s, int niter) { + s.ruiz_equilibration(niter); + return s.get_workspace().scaling; }, - py::arg("problem"), py::arg("scaling")) + py::arg("niter")) .def("get_problem", &Solver::get_problem, py::return_value_policy::reference_internal) // Main solve @@ -238,5 +266,15 @@ PYBIND11_MODULE(rcqp, m) { .def_property_readonly("comp_L_inds", [](const Solver& s) { return s.comp_L_inds; }) .def_property_readonly("comp_R_inds", - [](const Solver& s) { return s.comp_R_inds; }); + [](const Solver& s) { return s.comp_R_inds; }) + .def_property_readonly("z_z_inds", + [](const Solver& s) { return s.z_z_inds; }) + .def_property_readonly("s_ineq_s_ineq_inds", + [](const Solver& s) { return s.s_ineq_s_ineq_inds; }) + .def_property_readonly("s_ineq_m_ineq_inds", + [](const Solver& s) { return s.s_ineq_m_ineq_inds; }) + .def_property_readonly("s_comp_s_comp_inds", + [](const Solver& s) { return s.s_comp_s_comp_inds; }) + .def_property_readonly("s_comp_m_comp_inds", + [](const Solver& s) { return s.s_comp_m_comp_inds; }); } diff --git a/tester.py b/tester.py index 6c003eb..8613e65 100644 --- a/tester.py +++ b/tester.py @@ -1,8 +1,4 @@ import rcqp - -print("module:", rcqp) -print("members:", dir(rcqp)) - import numpy as np s = rcqp.Solver() @@ -23,6 +19,7 @@ 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 = rcqp.SolverOptions() diff --git a/typings/rcqp.pyi b/typings/rcqp.pyi index 532a0ec..efc3029 100644 --- a/typings/rcqp.pyi +++ b/typings/rcqp.pyi @@ -3,7 +3,6 @@ RCQP: constrained optimization solver with complementarity constraints """ from __future__ import annotations import numpy -import numpy.typing import scipy.sparse import typing __all__: list[str] = ['Filter', 'FilterEntry', 'Problem', 'Solver', 'SolverOptions', 'Workspace'] @@ -12,7 +11,7 @@ class Filter: def __init__(self) -> None: ... @typing.overload - def __init__(self, gamma_objective: typing.SupportsFloat | typing.SupportsIndex, gamma_constraint: typing.SupportsFloat | typing.SupportsIndex) -> None: + def __init__(self, gamma_objective: float, gamma_constraint: float) -> None: ... def acceptable(self, candidate: FilterEntry) -> bool: ... @@ -33,34 +32,24 @@ class Filter: def size(self) -> int: ... class FilterEntry: + feas: float + merit: float def __init__(self) -> None: ... def __repr__(self) -> str: ... - @property - def feas(self) -> float: - ... - @feas.setter - def feas(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: - ... - @property - def merit(self) -> float: - ... - @merit.setter - def merit(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: - ... class Problem: @typing.overload - def __init__(self, cost_hessian: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, n]"], cost_gradient: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], cost_const: typing.SupportsFloat | typing.SupportsIndex, J_eq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, n]"], c_eq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], J_ineq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, n]"], c_ineq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], J_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, n]"], c_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"]) -> None: + def __init__(self, cost_hessian: numpy.ndarray, cost_gradient: numpy.ndarray, cost_const: float, J_eq: numpy.ndarray, c_eq: numpy.ndarray, J_ineq: numpy.ndarray, c_ineq: numpy.ndarray, J_comp: numpy.ndarray, c_comp: numpy.ndarray) -> None: ... @typing.overload - def __init__(self, cost_hessian: scipy.sparse.csc_matrix, cost_gradient: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], cost_const: typing.SupportsFloat | typing.SupportsIndex, J_eq: scipy.sparse.csc_matrix, c_eq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], J_ineq: scipy.sparse.csc_matrix, c_ineq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], J_comp: scipy.sparse.csc_matrix, c_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"]) -> None: + 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, J_comp: scipy.sparse.csc_matrix, c_comp: numpy.ndarray) -> None: ... @property def cost_const(self) -> float: ... @property - def cost_gradient(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def cost_gradient(self) -> numpy.ndarray: ... @property def n_comp(self) -> int: @@ -91,9 +80,9 @@ class Solver: ... def convergence(self, options: SolverOptions) -> bool: ... - def entry_from_solution(self, sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex) -> tuple[float, float]: + def entry_from_solution(self, sqrt_relax_param: float, inv_penalty_param: float) -> tuple: ... - def filter_linesearch(self, sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex, max_iters: typing.SupportsInt | typing.SupportsIndex) -> bool: + def filter_linesearch(self, sqrt_relax_param: float, inv_penalty_param: float, max_iters: int) -> bool: ... def get_filter(self) -> Filter: ... @@ -105,42 +94,44 @@ class Solver: ... def numerical_factorization(self) -> bool: ... - def retract(self, s: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def retract(self, s: numpy.ndarray, sqrt_relax_param: float) -> numpy.ndarray: ... - def retract_deriv(self, s: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def retract_deriv(self, s: numpy.ndarray, sqrt_relax_param: float) -> numpy.ndarray: ... - def retract_second_deriv(self, s: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def retract_second_deriv(self, s: numpy.ndarray, sqrt_relax_param: float) -> numpy.ndarray: ... - def set_problem(self, problem: Problem, scaling: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"]) -> None: + def ruiz_equilibration(self, niter: int) -> numpy.ndarray: + ... + def set_problem(self, problem: Problem) -> None: ... def solve(self, options: SolverOptions) -> bool: ... - def update_KKT_comp(self, s_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], m_comp: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> None: + def update_KKT_comp(self, s_comp: numpy.ndarray, m_comp: numpy.ndarray, sqrt_relax_param: float) -> None: ... - def update_KKT_ineq(self, s_ineq: typing.Annotated[numpy.typing.ArrayLike, numpy.float64, "[m, 1]"], sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex) -> None: + def update_KKT_ineq(self, s_ineq: numpy.ndarray, sqrt_relax_param: float) -> None: ... - def update_KKT_penalty(self, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex) -> None: + def update_KKT_penalty(self, inv_penalty_param: float) -> None: ... - def update_KKT_primal_regularizer(self, reg: typing.SupportsFloat | typing.SupportsIndex) -> None: + def update_KKT_primal_regularizer(self, reg: float) -> None: ... - def update_KKT_residual(self, sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex) -> None: + def update_KKT_residual(self, sqrt_relax_param: float, inv_penalty_param: float) -> None: ... - def update_KKT_system(self, sqrt_relax_param: typing.SupportsFloat | typing.SupportsIndex, inv_penalty_param: typing.SupportsFloat | typing.SupportsIndex) -> None: + def update_KKT_system(self, sqrt_relax_param: float, inv_penalty_param: float) -> None: ... @property - def comp_L_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: + def comp_L_inds(self) -> numpy.ndarray: ... @property - def comp_R_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: + def comp_R_inds(self) -> numpy.ndarray: ... @property - def m_comp_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: + def m_comp_inds(self) -> numpy.ndarray: ... @property - def m_eq_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: + def m_eq_inds(self) -> numpy.ndarray: ... @property - def m_ineq_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: + def m_ineq_inds(self) -> numpy.ndarray: ... @property def n_duals(self) -> int: @@ -152,134 +143,87 @@ class Solver: def n_vars(self) -> int: ... @property - def s_comp_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: - ... - @property - def s_ineq_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: - ... - @property - def z_inds(self) -> typing.Annotated[numpy.typing.NDArray[numpy.int32], "[m, 1]"]: - ... -class SolverOptions: - output_dir: str - def __init__(self) -> None: - ... - def __repr__(self) -> str: - ... - @property - def convergence_comp_violation(self) -> float: - ... - @convergence_comp_violation.setter - def convergence_comp_violation(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: - ... - @property - def convergence_eq_violation(self) -> float: - ... - @convergence_eq_violation.setter - def convergence_eq_violation(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: - ... - @property - def convergence_ineq_violation(self) -> float: - ... - @convergence_ineq_violation.setter - def convergence_ineq_violation(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def s_comp_inds(self) -> numpy.ndarray: ... @property - def convergence_kkt_norm(self) -> float: - ... - @convergence_kkt_norm.setter - def convergence_kkt_norm(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def s_comp_m_comp_inds(self) -> numpy.ndarray: ... @property - def gamma_constraint(self) -> float: - ... - @gamma_constraint.setter - def gamma_constraint(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def s_comp_s_comp_inds(self) -> numpy.ndarray: ... @property - def gamma_objective(self) -> float: - ... - @gamma_objective.setter - def gamma_objective(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def s_ineq_inds(self) -> numpy.ndarray: ... @property - def max_iters(self) -> int: - ... - @max_iters.setter - def max_iters(self, arg0: typing.SupportsInt | typing.SupportsIndex) -> None: + def s_ineq_m_ineq_inds(self) -> numpy.ndarray: ... @property - def max_iters_linesearch(self) -> int: - ... - @max_iters_linesearch.setter - def max_iters_linesearch(self, arg0: typing.SupportsInt | typing.SupportsIndex) -> None: + def s_ineq_s_ineq_inds(self) -> numpy.ndarray: ... @property - def outer_step_kkt_norm(self) -> float: - ... - @outer_step_kkt_norm.setter - def outer_step_kkt_norm(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def z_inds(self) -> numpy.ndarray: ... @property - def penalty_initial(self) -> float: + def z_z_inds(self) -> numpy.ndarray: ... - @penalty_initial.setter - def penalty_initial(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: - ... - @property - def penalty_max(self) -> float: +class SolverOptions: + convergence_comp_violation: float + convergence_eq_violation: float + convergence_ineq_violation: float + convergence_kkt_norm: float + gamma_constraint: float + gamma_objective: float + max_iters: int + max_iters_linesearch: int + outer_step_kkt_norm: float + output_dir: str + penalty_initial: float + penalty_max: float + penalty_scaling: float + relaxation_initial: float + relaxation_min: float + relaxation_scaling: float + ruiz_iterations: float + def __init__(self) -> None: ... - @penalty_max.setter - def penalty_max(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def __repr__(self) -> str: ... +class Workspace: @property - def penalty_scaling(self) -> float: - ... - @penalty_scaling.setter - def penalty_scaling(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def D(self) -> numpy.ndarray: ... @property - def relaxation_initial(self) -> float: - ... - @relaxation_initial.setter - def relaxation_initial(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def amd_iperm_vec(self) -> numpy.ndarray: ... @property - def relaxation_min(self) -> float: - ... - @relaxation_min.setter - def relaxation_min(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def amd_perm_vec(self) -> numpy.ndarray: ... @property - def relaxation_scaling(self) -> float: - ... - @relaxation_scaling.setter - def relaxation_scaling(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None: + def kkt_residual(self) -> numpy.ndarray: ... -class Workspace: @property - def kkt_residual(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def kkt_system(self) -> tuple: ... @property - def m_comp(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def m_comp(self) -> numpy.ndarray: ... @property - def m_comp_est(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def m_comp_est(self) -> numpy.ndarray: ... @property - def m_eq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def m_eq(self) -> numpy.ndarray: ... @property - def m_eq_est(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def m_eq_est(self) -> numpy.ndarray: ... @property - def m_ineq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def m_ineq(self) -> numpy.ndarray: ... @property - def m_ineq_est(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def m_ineq_est(self) -> numpy.ndarray: ... @property - def newton_step(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def newton_step(self) -> numpy.ndarray: ... @property def penalty_param(self) -> float: @@ -288,23 +232,26 @@ class Workspace: def relax_param(self) -> float: ... @property - def residual_comp(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def residual_comp(self) -> numpy.ndarray: + ... + @property + def residual_eq(self) -> numpy.ndarray: ... @property - def residual_eq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def residual_ineq(self) -> numpy.ndarray: ... @property - def residual_ineq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def s_comp(self) -> numpy.ndarray: ... @property - def s_comp(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def s_ineq(self) -> numpy.ndarray: ... @property - def s_ineq(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def scaling(self) -> numpy.ndarray: ... @property - def solution(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def solution(self) -> numpy.ndarray: ... @property - def z(self) -> typing.Annotated[numpy.typing.NDArray[numpy.float64], "[m, 1]"]: + def z(self) -> numpy.ndarray: ... From f6cbef4691c44b5e4d1c76fee485bff706c632b4 Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Mon, 6 Apr 2026 23:56:15 -0400 Subject: [PATCH 7/9] test build req --- .github/workflows/cmake-build.yml | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 .github/workflows/cmake-build.yml diff --git a/.github/workflows/cmake-build.yml b/.github/workflows/cmake-build.yml new file mode 100644 index 0000000..19d7ab1 --- /dev/null +++ b/.github/workflows/cmake-build.yml @@ -0,0 +1,28 @@ +name: CMake Build + +on: + pull_request: + branches: [master] + +jobs: + build: + runs-on: ubuntu-24.04 + + steps: + - uses: actions/checkout@v4 + + - name: Install system dependencies + run: | + sudo apt-get update -q + sudo apt-get install -y libeigen3-dev nlohmann-json3-dev + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.12" + + - name: Configure (Python bindings, no stub generation) + run: cmake --preset python -DRCQP_GENERATE_PYI=OFF + + - name: Build + run: cmake --build --preset python -- -j$(nproc) From 53e6116ab20c0efedae1af86d7942274f3f88821 Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Mon, 6 Apr 2026 23:57:50 -0400 Subject: [PATCH 8/9] update workflow --- .github/workflows/cmake-build.yml | 61 +++++++++++++++++++++++++++++-- 1 file changed, 58 insertions(+), 3 deletions(-) diff --git a/.github/workflows/cmake-build.yml b/.github/workflows/cmake-build.yml index 19d7ab1..d63420b 100644 --- a/.github/workflows/cmake-build.yml +++ b/.github/workflows/cmake-build.yml @@ -5,9 +5,9 @@ on: branches: [master] jobs: - build: + build-python: + name: Python bindings runs-on: ubuntu-24.04 - steps: - uses: actions/checkout@v4 @@ -21,8 +21,63 @@ jobs: with: python-version: "3.12" - - name: Configure (Python bindings, no stub generation) + - name: Configure run: cmake --preset python -DRCQP_GENERATE_PYI=OFF - name: Build run: cmake --build --preset python -- -j$(nproc) + + build-julia: + name: Julia bindings + runs-on: ubuntu-24.04 + steps: + - uses: actions/checkout@v4 + + - name: Install system dependencies + run: | + sudo apt-get update -q + sudo apt-get install -y libeigen3-dev nlohmann-json3-dev + + - name: Set up Julia + uses: julia-actions/setup-julia@v2 + with: + version: "1" + + - name: Install CxxWrap.jl + run: julia --startup-file=no -e 'using Pkg; Pkg.add("CxxWrap")' + + - name: Configure + run: cmake --preset julia + + - name: Build + run: cmake --build --preset julia -- -j$(nproc) + + build-all: + name: Python + Julia bindings + runs-on: ubuntu-24.04 + steps: + - uses: actions/checkout@v4 + + - name: Install system dependencies + run: | + sudo apt-get update -q + sudo apt-get install -y libeigen3-dev nlohmann-json3-dev + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.12" + + - name: Set up Julia + uses: julia-actions/setup-julia@v2 + with: + version: "1" + + - name: Install CxxWrap.jl + run: julia --startup-file=no -e 'using Pkg; Pkg.add("CxxWrap")' + + - name: Configure + run: cmake --preset all -DRCQP_GENERATE_PYI=OFF + + - name: Build + run: cmake --build --preset all -- -j$(nproc) From 8ba84384f12ca22c1eea7bdc9051efda5ac0942c Mon Sep 17 00:00:00 2001 From: Micah Reich Date: Tue, 7 Apr 2026 00:02:56 -0400 Subject: [PATCH 9/9] update readme --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 0f53c36..c245fe1 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,6 @@ # RCQP +[![CMake Build](https://github.com/MarbleSolver/RCQP/actions/workflows/cmake-build.yml/badge.svg?branch=master)](https://github.com/MarbleSolver/RCQP/actions/workflows/cmake-build.yml) + A C++ solver for quadratic programs with linear complementarity constraints, with Python and Julia bindings. ## Prerequisites