From 97d2b1e6dc963a35c024af64ab8d48a5f7ff72fa Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 19 Aug 2026 06:32:08 -0500 Subject: [PATCH 1/2] Add compiled Python bindings Build a shared library with juliac and expose the MatrixCovers API through ctypes and NumPy. Copy inputs at the ABI boundary and report failures through structured status values. The compiled library includes solvers that do not require package extensions. Unsupported penalty and solver combinations return an error. Assisted-by: Claude Sonnet 5 Assisted-by: Claude Fable 5 --- .gitignore | 3 + lib/Project.toml | 15 ++ lib/build-env/Project.toml | 13 + lib/build.jl | 45 ++++ lib/python/_facade.py | 278 +++++++++++++++++++++ lib/src/matrixcovers.jl | 455 ++++++++++++++++++++++++++++++++++ lib/test/python/test_smoke.py | 107 ++++++++ 7 files changed, 916 insertions(+) create mode 100644 lib/Project.toml create mode 100644 lib/build-env/Project.toml create mode 100644 lib/build.jl create mode 100644 lib/python/_facade.py create mode 100644 lib/src/matrixcovers.jl create mode 100644 lib/test/python/test_smoke.py diff --git a/.gitignore b/.gitignore index 19ab606..5f09308 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,6 @@ /docs/Manifest*.toml /docs/build/ /test/generator/Manifest*.toml +/lib/out/ +/lib/Manifest.toml +/lib/build-env/Manifest.toml diff --git a/lib/Project.toml b/lib/Project.toml new file mode 100644 index 0000000..982870f --- /dev/null +++ b/lib/Project.toml @@ -0,0 +1,15 @@ +name = "matrixcovers" +uuid = "9a0ce11c-0abe-4c5a-899b-d53d9853b6d8" +version = "0.0.1" + +[deps] +JLWInterop = "65e54657-ed21-41a3-96db-71ab7fa6d94b" +MatrixCovers = "727e6139-ff52-4636-a344-ed1d23e73ffc" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +# `build.jl` adds the local MatrixCovers source to a temporary copy. + +[compat] +JLWInterop = "0.1" +MatrixCovers = "1" +julia = "1.13" diff --git a/lib/build-env/Project.toml b/lib/build-env/Project.toml new file mode 100644 index 0000000..47a4a6c --- /dev/null +++ b/lib/build-env/Project.toml @@ -0,0 +1,13 @@ +# Build dependencies for the matrixcovers Python wheel. The runtime +# dependencies are defined in `../Project.toml`. +# +# julia +rc --project=lib/build-env -e 'using Pkg; Pkg.instantiate()' + +[deps] +JuliaC = "acedd4c2-ced6-4a15-accc-2607eb759ba2" +JuliaLibWrapping = "d61f35a8-f6af-436f-bc10-cee6b101f7bd" + +[compat] +JuliaC = "0.3" +JuliaLibWrapping = "0.2" +julia = "1.13" diff --git a/lib/build.jl b/lib/build.jl new file mode 100644 index 0000000..33627b4 --- /dev/null +++ b/lib/build.jl @@ -0,0 +1,45 @@ +# Build the Python wheel from the repository root with Julia 1.13: +# +# julia +rc --project=lib/build-env lib/build.jl +# +# Instantiate the build environment first: +# +# julia +rc --project=lib/build-env -e 'using Pkg; Pkg.instantiate()' +# +# The temporary project supplies the absolute source path required by juliac. + +using TOML: TOML + +const HERE = @__DIR__ +const REPO_ROOT = abspath(joinpath(HERE, "..")) + +function prepare_project() + toml = TOML.parsefile(joinpath(HERE, "Project.toml")) + sources = get(toml, "sources", Dict{String, Any}()) + sources["MatrixCovers"] = Dict("path" => REPO_ROOT) + toml["sources"] = sources + tmp = mktempdir(; prefix = "matrixcovers-lib-project-") + open(joinpath(tmp, "Project.toml"), "w") do io + TOML.print(io, toml; sorted = true) + end + return tmp +end + +const REPO_VERSION = TOML.parsefile(joinpath(REPO_ROOT, "Project.toml"))["version"] + +using JuliaLibWrapping, JuliaC + +result = standard_build(HERE; + libname = "matrixcovers", + python_package = "matrixcovers", + project = prepare_project(), + version = REPO_VERSION, + verbose = true, +) + +# Replace the generated facade with the public API. +cp(joinpath(HERE, "python", "_facade.py"), + joinpath(HERE, "out", "matrixcovers", "_facade.py"); + force = true) + +@info "Built matrixcovers" library=result.library bundle=result.bundle_dir diff --git a/lib/python/_facade.py b/lib/python/_facade.py new file mode 100644 index 0000000..2b7c72f --- /dev/null +++ b/lib/python/_facade.py @@ -0,0 +1,278 @@ +"""Python bindings for MatrixCovers.jl. + +Inputs are converted to `float64` arrays and outputs are newly allocated +`numpy.ndarray`s. Invalid options raise `ValueError`; errors reported by the +compiled library raise `JLWError`. + +Penalties are `"abslog1"`, `"abslog2"`, `"abslinear1"`, and `"abslinear2"`. +The `_min` functions support only `"abslog2"`; other penalties require Julia +package extensions that are not included in the compiled library. +""" +from . import _lowlevel +import numpy as np + +from ._lowlevel import JLWError + +_PENALTY_CODES = { + "abslog1": 1, + "abslog2": 2, + "abslinear1": 3, + "abslinear2": 4, +} + + +def _penalty_code(penalty): + try: + return _PENALTY_CODES[penalty] + except KeyError: + raise ValueError( + f"unknown penalty {penalty!r}; expected one of {sorted(_PENALTY_CODES)}" + ) from None + + +_LINSOLVE_CODES = {"auto": 1, "dense": 2, "lsqr": 3} + + +def _linsolve_code(linsolve): + try: + return _LINSOLVE_CODES[linsolve] + except KeyError: + raise ValueError( + f"unknown linsolve {linsolve!r}; expected one of {sorted(_LINSOLVE_CODES)}" + ) from None + + +def _sentinel_int(value): + return -1 if value is None else int(value) + + +def _sentinel_float(value): + return -1.0 if value is None else float(value) + + +def _as_matrix(A): + return np.asfortranarray(A, dtype=np.float64) + + +def _as_vector(v): + return np.ascontiguousarray(v, dtype=np.float64) + + +def symcover(A, *, maxiter=None): + """Heuristic symmetric hard cover: `a` with `a[i]*a[j] >= abs(A[i, j])`.""" + _A = _as_matrix(A) + a = np.zeros(_A.shape[0], dtype=np.float64) + _lowlevel.mc_symcover( + _lowlevel.CMatrix_Float64.from_numpy(_A), + _sentinel_int(maxiter), + _lowlevel.CVector_Float64.from_numpy(a), + ) + return a + + +def cover(A, *, maxiter=None): + """Heuristic hard cover: `(a, b)` with `a[i]*b[j] >= abs(A[i, j])`.""" + _A = _as_matrix(A) + m, n = _A.shape + a = np.zeros(m, dtype=np.float64) + b = np.zeros(n, dtype=np.float64) + _lowlevel.mc_cover( + _lowlevel.CMatrix_Float64.from_numpy(_A), + _sentinel_int(maxiter), + _lowlevel.CVector_Float64.from_numpy(a), + _lowlevel.CVector_Float64.from_numpy(b), + ) + return a, b + + +def symcover_min(A, *, penalty="abslog2", maxiter=None, linsolve="auto"): + """phi-minimal symmetric hard cover of `A`.""" + _A = _as_matrix(A) + a = np.zeros(_A.shape[0], dtype=np.float64) + _lowlevel.mc_symcover_min( + _penalty_code(penalty), + _lowlevel.CMatrix_Float64.from_numpy(_A), + _sentinel_int(maxiter), + _linsolve_code(linsolve), + _lowlevel.CVector_Float64.from_numpy(a), + ) + return a + + +def cover_min(A, *, penalty="abslog2", maxiter=None, linsolve="auto"): + """phi-minimal hard cover of `A`.""" + _A = _as_matrix(A) + m, n = _A.shape + a = np.zeros(m, dtype=np.float64) + b = np.zeros(n, dtype=np.float64) + _lowlevel.mc_cover_min( + _penalty_code(penalty), + _lowlevel.CMatrix_Float64.from_numpy(_A), + _sentinel_int(maxiter), + _linsolve_code(linsolve), + _lowlevel.CVector_Float64.from_numpy(a), + _lowlevel.CVector_Float64.from_numpy(b), + ) + return a, b + + +def soft_symcover(A, *, penalty="abslinear2", maxiter=None, starts=None, sigma=None, seed=0): + """Symmetric soft cover of `A` minimizing the penalty, with no coverage constraint.""" + _A = _as_matrix(A) + a = np.zeros(_A.shape[0], dtype=np.float64) + _lowlevel.mc_soft_symcover( + _penalty_code(penalty), + _lowlevel.CMatrix_Float64.from_numpy(_A), + _sentinel_int(maxiter), + _sentinel_int(starts), + _sentinel_float(sigma), + int(seed), + _lowlevel.CVector_Float64.from_numpy(a), + ) + return a + + +def soft_cover(A, *, penalty="abslinear2", maxiter=None, starts=None, sigma=None, seed=0): + """Asymmetric soft cover of `A` minimizing the penalty, with no coverage constraint.""" + _A = _as_matrix(A) + m, n = _A.shape + a = np.zeros(m, dtype=np.float64) + b = np.zeros(n, dtype=np.float64) + _lowlevel.mc_soft_cover( + _penalty_code(penalty), + _lowlevel.CMatrix_Float64.from_numpy(_A), + _sentinel_int(maxiter), + _sentinel_int(starts), + _sentinel_float(sigma), + int(seed), + _lowlevel.CVector_Float64.from_numpy(a), + _lowlevel.CVector_Float64.from_numpy(b), + ) + return a, b + + +def soft_symcover_min(A, *, penalty="abslog2", maxiter=None): + """phi-minimal symmetric soft cover of `A`, with no coverage constraint.""" + _A = _as_matrix(A) + a = np.zeros(_A.shape[0], dtype=np.float64) + _lowlevel.mc_soft_symcover_min( + _penalty_code(penalty), + _lowlevel.CMatrix_Float64.from_numpy(_A), + _sentinel_int(maxiter), + _lowlevel.CVector_Float64.from_numpy(a), + ) + return a + + +def soft_cover_min(A, *, penalty="abslog2", maxiter=None): + """phi-minimal asymmetric soft cover of `A`, with no coverage constraint.""" + _A = _as_matrix(A) + m, n = _A.shape + a = np.zeros(m, dtype=np.float64) + b = np.zeros(n, dtype=np.float64) + _lowlevel.mc_soft_cover_min( + _penalty_code(penalty), + _lowlevel.CMatrix_Float64.from_numpy(_A), + _sentinel_int(maxiter), + _lowlevel.CVector_Float64.from_numpy(a), + _lowlevel.CVector_Float64.from_numpy(b), + ) + return a, b + + +def iscover(a, A, b=None, *, rtol=0.0, atol=0.0): + """Whether `a`, `b` cover `A`: `a[i]*b[j] >= abs(A[i, j])*(1 - rtol) - atol`. + + `b=None` (the default) tests the symmetric cover `a*a'`, and requires `A` + to be square. + """ + _A = _as_matrix(A) + _a = _as_vector(a) + if b is None: + _result = _lowlevel.mc_iscover_sym( + _lowlevel.CVector_Float64.from_numpy(_a), + _lowlevel.CMatrix_Float64.from_numpy(_A), + float(rtol), float(atol), + ) + else: + _b = _as_vector(b) + _result = _lowlevel.mc_iscover( + _lowlevel.CVector_Float64.from_numpy(_a), + _lowlevel.CVector_Float64.from_numpy(_b), + _lowlevel.CMatrix_Float64.from_numpy(_A), + float(rtol), float(atol), + ) + return bool(_result.value) + + +def cover_objective(a, A, b=None, *, penalty="abslog2"): + """`sum(phi(abs(A[i,j]) / (a[i]*b[j])))`; `b=None` tests the symmetric cover `a*a'`.""" + _A = _as_matrix(A) + _a = _as_vector(a) + if b is None: + _result = _lowlevel.mc_cover_objective_sym( + _penalty_code(penalty), + _lowlevel.CVector_Float64.from_numpy(_a), + _lowlevel.CMatrix_Float64.from_numpy(_A), + ) + else: + _b = _as_vector(b) + _result = _lowlevel.mc_cover_objective( + _penalty_code(penalty), + _lowlevel.CVector_Float64.from_numpy(_a), + _lowlevel.CVector_Float64.from_numpy(_b), + _lowlevel.CMatrix_Float64.from_numpy(_A), + ) + return float(_result.value) + + +def gramcover(a, b, A, *, w=None, W=None): + """Symmetric cover of a (weighted) Gram matrix of `A`, from a cover `(a, b)` of `A`. + + `G = A'*A` when neither `w` nor `W` is given, `G = A'*diag(w)*A` for a + vector `w`, and `G = A'*W*A` for a matrix `W`. Passing both `w` and `W` + raises `ValueError`. + """ + if w is not None and W is not None: + raise ValueError("pass at most one of `w` or `W`, not both") + _a = _as_vector(a) + _b = _as_vector(b) + _A = _as_matrix(A) + s = np.zeros(_A.shape[1], dtype=np.float64) + if w is not None: + _w = _as_vector(w) + _lowlevel.mc_gramcover_weighted( + _lowlevel.CVector_Float64.from_numpy(_a), + _lowlevel.CVector_Float64.from_numpy(_b), + _lowlevel.CMatrix_Float64.from_numpy(_A), + _lowlevel.CVector_Float64.from_numpy(_w), + _lowlevel.CVector_Float64.from_numpy(s), + ) + elif W is not None: + _W = _as_matrix(W) + _lowlevel.mc_gramcover_matrix( + _lowlevel.CVector_Float64.from_numpy(_a), + _lowlevel.CVector_Float64.from_numpy(_b), + _lowlevel.CMatrix_Float64.from_numpy(_A), + _lowlevel.CMatrix_Float64.from_numpy(_W), + _lowlevel.CVector_Float64.from_numpy(s), + ) + else: + _lowlevel.mc_gramcover( + _lowlevel.CVector_Float64.from_numpy(_a), + _lowlevel.CVector_Float64.from_numpy(_b), + _lowlevel.CMatrix_Float64.from_numpy(_A), + _lowlevel.CVector_Float64.from_numpy(s), + ) + return s + + +__all__ = [ + "JLWError", + "symcover", "cover", + "symcover_min", "cover_min", + "soft_symcover", "soft_cover", + "soft_symcover_min", "soft_cover_min", + "iscover", "cover_objective", + "gramcover", +] diff --git a/lib/src/matrixcovers.jl b/lib/src/matrixcovers.jl new file mode 100644 index 0000000..f4d5bfe --- /dev/null +++ b/lib/src/matrixcovers.jl @@ -0,0 +1,455 @@ +# C ABI for MatrixCovers. Inputs and caller-owned output buffers contain +# `Float64`; matrices are column-major. Inputs are copied before calling +# MatrixCovers, and output sizes are checked before copying results back. +# +# # Penalty enum (`penalty::Int32`) +# +# 1 = AbsLog{1}() 2 = AbsLog{2}() 3 = AbsLinear{1}() 4 = AbsLinear{2}() +# +# Penalties are decoded with concrete branches to keep dispatch trim-safe. +# The `_min` entry points support only `AbsLog{2}` because other penalties +# require package extensions that are not linked. Soft covers and +# `cover_objective` support all four penalties. +# +# # linsolve enum (`linsolve::Int32`) +# +# 1 = :auto 2 = :dense 3 = :lsqr +# +# Used by `mc_symcover_min` and `mc_cover_min`. +# +# # Numeric tuning knobs +# +# Negative `maxiter`, `starts`, or `sigma` values omit that keyword, preserving +# MatrixCovers defaults. `starts`, `sigma`, and `seed` apply only to +# `AbsLinear` soft covers. +# +# # Status codes (`JLWStatus.code`) +# +# 0 ok +# 1 invalid penalty enum +# 2 penalty requires an extension not linked into this library +# 3 invalid linsolve enum +# 4 DimensionMismatch +# 6 ArgumentError +# 99 unexpected internal error +module matrixcovers + +using JLWInterop +using MatrixCovers +using Random: MersenneTwister + +# Results that contain both a status and a value. + +struct MCBoolResult + status::JLWStatus + value::Int32 +end + +struct MCScalarResult + status::JLWStatus + value::Float64 +end + +# Boundary helpers + +function _copyin_vec(v::CVector{Float64}) + out = Vector{Float64}(undef, length(v)) + copyto!(out, v) + return out +end + +function _copyin_mat(A::CMatrix{Float64}) + m, n = size(A) + out = Matrix{Float64}(undef, m, n) + copyto!(out, A) + return out +end + +# Explicit comparisons keep the mapping trim-safe. +function _linsolve_symbol(linsolve::Int32) + linsolve == Int32(1) && return :auto + linsolve == Int32(2) && return :dense + linsolve == Int32(3) && return :lsqr + return :invalid +end + +# Inline exception handling so `trim=:safe` can narrow the exception type. +macro status_from_exception(e) + quote + let ex = $(esc(e)) + if ex isa DimensionMismatch + m = ex.msg + jlw_error(4, m isa String ? m : "dimension mismatch") + elseif ex isa ArgumentError + m = ex.msg + jlw_error(6, m isa String ? m : "invalid argument") + else + jlw_error(99, "unexpected internal error") + end + end + end +end + +# soft_symcover / soft_cover kwarg plumbing +# +# Penalty types and explicit keyword combinations keep calls trim-safe. + +function _soft_symcover_log(::Type{P}, M::Matrix{Float64}, maxiter::Int64) where {P<:MatrixCovers.AbstractCoverPenalty} + return maxiter < 0 ? soft_symcover(P(), M) : soft_symcover(P(), M; maxiter = Int(maxiter)) +end + +function _soft_cover_log(::Type{P}, M::Matrix{Float64}, maxiter::Int64) where {P<:MatrixCovers.AbstractCoverPenalty} + return maxiter < 0 ? soft_cover(P(), M) : soft_cover(P(), M; maxiter = Int(maxiter)) +end + +# `AbsLinear` solvers also accept `starts`, `sigma`, and `rng`. +function _soft_symcover_lin(::Type{P}, M::Matrix{Float64}, maxiter::Int64, starts::Int64, + sigma::Float64, seed::UInt64) where {P<:MatrixCovers.AbstractCoverPenalty} + rng = MersenneTwister(seed) + if maxiter < 0 + if starts < 0 + return sigma < 0 ? soft_symcover(P(), M; rng) : + soft_symcover(P(), M; sigma, rng) + else + return sigma < 0 ? soft_symcover(P(), M; starts = Int(starts), rng) : + soft_symcover(P(), M; starts = Int(starts), sigma, rng) + end + else + if starts < 0 + return sigma < 0 ? soft_symcover(P(), M; maxiter = Int(maxiter), rng) : + soft_symcover(P(), M; maxiter = Int(maxiter), sigma, rng) + else + return sigma < 0 ? soft_symcover(P(), M; maxiter = Int(maxiter), starts = Int(starts), rng) : + soft_symcover(P(), M; maxiter = Int(maxiter), starts = Int(starts), sigma, rng) + end + end +end + +function _soft_cover_lin(::Type{P}, M::Matrix{Float64}, maxiter::Int64, starts::Int64, + sigma::Float64, seed::UInt64) where {P<:MatrixCovers.AbstractCoverPenalty} + rng = MersenneTwister(seed) + if maxiter < 0 + if starts < 0 + return sigma < 0 ? soft_cover(P(), M; rng) : + soft_cover(P(), M; sigma, rng) + else + return sigma < 0 ? soft_cover(P(), M; starts = Int(starts), rng) : + soft_cover(P(), M; starts = Int(starts), sigma, rng) + end + else + if starts < 0 + return sigma < 0 ? soft_cover(P(), M; maxiter = Int(maxiter), rng) : + soft_cover(P(), M; maxiter = Int(maxiter), sigma, rng) + else + return sigma < 0 ? soft_cover(P(), M; maxiter = Int(maxiter), starts = Int(starts), rng) : + soft_cover(P(), M; maxiter = Int(maxiter), starts = Int(starts), sigma, rng) + end + end +end + +# Hard covers + +Base.@ccallable function mc_symcover(A::CMatrix{Float64}, maxiter::Int64, + a::CVector{Float64})::JLWStatus + try + M = _copyin_mat(A) + av = maxiter < 0 ? symcover(M) : symcover(M; maxiter = Int(maxiter)) + length(a) == length(av) || return jlw_error(4, "output length must match matrix size") + copyto!(a, av) + return jlw_ok() + catch e + return @status_from_exception(e) + end +end + +Base.@ccallable function mc_cover(A::CMatrix{Float64}, maxiter::Int64, + a::CVector{Float64}, b::CVector{Float64})::JLWStatus + try + M = _copyin_mat(A) + av, bv = maxiter < 0 ? cover(M) : cover(M; maxiter = Int(maxiter)) + (length(a) == length(av) && length(b) == length(bv)) || + return jlw_error(4, "output lengths must match matrix size") + copyto!(a, av) + copyto!(b, bv) + return jlw_ok() + catch e + return @status_from_exception(e) + end +end + +Base.@ccallable function mc_symcover_min(penalty::Int32, A::CMatrix{Float64}, maxiter::Int64, + linsolve::Int32, a::CVector{Float64})::JLWStatus + try + if penalty == Int32(2) + ls = _linsolve_symbol(linsolve) + ls === :invalid && return jlw_error(3, "linsolve must be 1 (:auto), 2 (:dense), or 3 (:lsqr)") + M = _copyin_mat(A) + av = maxiter < 0 ? symcover_min(AbsLog{2}(), M; linsolve = ls) : + symcover_min(AbsLog{2}(), M; maxiter = Int(maxiter), linsolve = ls) + length(a) == length(av) || return jlw_error(4, "output length must match matrix size") + copyto!(a, av) + return jlw_ok() + elseif penalty == Int32(1) + return jlw_error(2, "penalty AbsLog{1} requires the MatrixCoversJuMPExt extension (JuMP and HiGHS)") + elseif penalty == Int32(3) || penalty == Int32(4) + return jlw_error(2, "penalty AbsLinear requires the MatrixCoversIpoptExt extension (JuMP and Ipopt)") + else + return jlw_error(1, "unknown penalty enum, expected 1 (AbsLog1), 2 (AbsLog2), 3 (AbsLinear1), or 4 (AbsLinear2)") + end + catch e + return @status_from_exception(e) + end +end + +Base.@ccallable function mc_cover_min(penalty::Int32, A::CMatrix{Float64}, maxiter::Int64, + linsolve::Int32, a::CVector{Float64}, + b::CVector{Float64})::JLWStatus + try + if penalty == Int32(2) + ls = _linsolve_symbol(linsolve) + ls === :invalid && return jlw_error(3, "linsolve must be 1 (:auto), 2 (:dense), or 3 (:lsqr)") + M = _copyin_mat(A) + av, bv = maxiter < 0 ? cover_min(AbsLog{2}(), M; linsolve = ls) : + cover_min(AbsLog{2}(), M; maxiter = Int(maxiter), linsolve = ls) + (length(a) == length(av) && length(b) == length(bv)) || + return jlw_error(4, "output lengths must match matrix size") + copyto!(a, av) + copyto!(b, bv) + return jlw_ok() + elseif penalty == Int32(1) + return jlw_error(2, "penalty AbsLog{1} requires the MatrixCoversJuMPExt extension (JuMP and HiGHS)") + elseif penalty == Int32(3) || penalty == Int32(4) + return jlw_error(2, "penalty AbsLinear requires the MatrixCoversIpoptExt extension (JuMP and Ipopt)") + else + return jlw_error(1, "unknown penalty enum, expected 1 (AbsLog1), 2 (AbsLog2), 3 (AbsLinear1), or 4 (AbsLinear2)") + end + catch e + return @status_from_exception(e) + end +end + +# Soft covers + +Base.@ccallable function mc_soft_symcover(penalty::Int32, A::CMatrix{Float64}, maxiter::Int64, + starts::Int64, sigma::Float64, seed::UInt64, + a::CVector{Float64})::JLWStatus + try + M = _copyin_mat(A) + av = if penalty == Int32(1) + _soft_symcover_log(AbsLog{1}, M, maxiter) + elseif penalty == Int32(2) + _soft_symcover_log(AbsLog{2}, M, maxiter) + elseif penalty == Int32(3) + _soft_symcover_lin(AbsLinear{1}, M, maxiter, starts, sigma, seed) + elseif penalty == Int32(4) + _soft_symcover_lin(AbsLinear{2}, M, maxiter, starts, sigma, seed) + else + return jlw_error(1, "unknown penalty enum, expected 1 (AbsLog1), 2 (AbsLog2), 3 (AbsLinear1), or 4 (AbsLinear2)") + end + length(a) == length(av) || return jlw_error(4, "output length must match matrix size") + copyto!(a, av) + return jlw_ok() + catch e + return @status_from_exception(e) + end +end + +Base.@ccallable function mc_soft_cover(penalty::Int32, A::CMatrix{Float64}, maxiter::Int64, + starts::Int64, sigma::Float64, seed::UInt64, + a::CVector{Float64}, b::CVector{Float64})::JLWStatus + try + M = _copyin_mat(A) + av, bv = if penalty == Int32(1) + _soft_cover_log(AbsLog{1}, M, maxiter) + elseif penalty == Int32(2) + _soft_cover_log(AbsLog{2}, M, maxiter) + elseif penalty == Int32(3) + _soft_cover_lin(AbsLinear{1}, M, maxiter, starts, sigma, seed) + elseif penalty == Int32(4) + _soft_cover_lin(AbsLinear{2}, M, maxiter, starts, sigma, seed) + else + return jlw_error(1, "unknown penalty enum, expected 1 (AbsLog1), 2 (AbsLog2), 3 (AbsLinear1), or 4 (AbsLinear2)") + end + (length(a) == length(av) && length(b) == length(bv)) || + return jlw_error(4, "output lengths must match matrix size") + copyto!(a, av) + copyto!(b, bv) + return jlw_ok() + catch e + return @status_from_exception(e) + end +end + +Base.@ccallable function mc_soft_symcover_min(penalty::Int32, A::CMatrix{Float64}, maxiter::Int64, + a::CVector{Float64})::JLWStatus + try + if penalty == Int32(2) + M = _copyin_mat(A) + av = maxiter < 0 ? soft_symcover_min(AbsLog{2}(), M) : + soft_symcover_min(AbsLog{2}(), M; maxiter = Int(maxiter)) + length(a) == length(av) || return jlw_error(4, "output length must match matrix size") + copyto!(a, av) + return jlw_ok() + elseif penalty == Int32(1) + return jlw_error(2, "MatrixCovers does not implement AbsLog{1} for soft_symcover_min") + elseif penalty == Int32(3) || penalty == Int32(4) + return jlw_error(2, "penalty AbsLinear requires the MatrixCoversIpoptExt extension (JuMP and Ipopt)") + else + return jlw_error(1, "unknown penalty enum, expected 1 (AbsLog1), 2 (AbsLog2), 3 (AbsLinear1), or 4 (AbsLinear2)") + end + catch e + return @status_from_exception(e) + end +end + +Base.@ccallable function mc_soft_cover_min(penalty::Int32, A::CMatrix{Float64}, maxiter::Int64, + a::CVector{Float64}, b::CVector{Float64})::JLWStatus + try + if penalty == Int32(2) + M = _copyin_mat(A) + av, bv = maxiter < 0 ? soft_cover_min(AbsLog{2}(), M) : + soft_cover_min(AbsLog{2}(), M; maxiter = Int(maxiter)) + (length(a) == length(av) && length(b) == length(bv)) || + return jlw_error(4, "output lengths must match matrix size") + copyto!(a, av) + copyto!(b, bv) + return jlw_ok() + elseif penalty == Int32(1) + return jlw_error(2, "MatrixCovers does not implement AbsLog{1} for soft_cover_min") + elseif penalty == Int32(3) || penalty == Int32(4) + return jlw_error(2, "penalty AbsLinear requires the MatrixCoversIpoptExt extension (JuMP and Ipopt)") + else + return jlw_error(1, "unknown penalty enum, expected 1 (AbsLog1), 2 (AbsLog2), 3 (AbsLinear1), or 4 (AbsLinear2)") + end + catch e + return @status_from_exception(e) + end +end + +# Predicates and objectives + +Base.@ccallable function mc_iscover_sym(a::CVector{Float64}, A::CMatrix{Float64}, + rtol::Float64, atol::Float64)::MCBoolResult + try + av = _copyin_vec(a) + M = _copyin_mat(A) + ok = iscover(av, M; rtol, atol) + return MCBoolResult(jlw_ok(), ok ? Int32(1) : Int32(0)) + catch e + return MCBoolResult(@status_from_exception(e), Int32(0)) + end +end + +Base.@ccallable function mc_iscover(a::CVector{Float64}, b::CVector{Float64}, + A::CMatrix{Float64}, rtol::Float64, + atol::Float64)::MCBoolResult + try + av = _copyin_vec(a) + bv = _copyin_vec(b) + M = _copyin_mat(A) + ok = iscover(av, bv, M; rtol, atol) + return MCBoolResult(jlw_ok(), ok ? Int32(1) : Int32(0)) + catch e + return MCBoolResult(@status_from_exception(e), Int32(0)) + end +end + +Base.@ccallable function mc_cover_objective_sym(penalty::Int32, a::CVector{Float64}, + A::CMatrix{Float64})::MCScalarResult + try + av = _copyin_vec(a) + M = _copyin_mat(A) + if penalty == Int32(1) + return MCScalarResult(jlw_ok(), cover_objective(AbsLog{1}(), av, M)) + elseif penalty == Int32(2) + return MCScalarResult(jlw_ok(), cover_objective(AbsLog{2}(), av, M)) + elseif penalty == Int32(3) + return MCScalarResult(jlw_ok(), cover_objective(AbsLinear{1}(), av, M)) + elseif penalty == Int32(4) + return MCScalarResult(jlw_ok(), cover_objective(AbsLinear{2}(), av, M)) + else + return MCScalarResult(jlw_error(1, "unknown penalty enum, expected 1 (AbsLog1), 2 (AbsLog2), 3 (AbsLinear1), or 4 (AbsLinear2)"), NaN) + end + catch e + return MCScalarResult(@status_from_exception(e), NaN) + end +end + +Base.@ccallable function mc_cover_objective(penalty::Int32, a::CVector{Float64}, + b::CVector{Float64}, + A::CMatrix{Float64})::MCScalarResult + try + av = _copyin_vec(a) + bv = _copyin_vec(b) + M = _copyin_mat(A) + if penalty == Int32(1) + return MCScalarResult(jlw_ok(), cover_objective(AbsLog{1}(), av, bv, M)) + elseif penalty == Int32(2) + return MCScalarResult(jlw_ok(), cover_objective(AbsLog{2}(), av, bv, M)) + elseif penalty == Int32(3) + return MCScalarResult(jlw_ok(), cover_objective(AbsLinear{1}(), av, bv, M)) + elseif penalty == Int32(4) + return MCScalarResult(jlw_ok(), cover_objective(AbsLinear{2}(), av, bv, M)) + else + return MCScalarResult(jlw_error(1, "unknown penalty enum, expected 1 (AbsLog1), 2 (AbsLog2), 3 (AbsLinear1), or 4 (AbsLinear2)"), NaN) + end + catch e + return MCScalarResult(@status_from_exception(e), NaN) + end +end + +# Gram covers + +Base.@ccallable function mc_gramcover(a::CVector{Float64}, b::CVector{Float64}, + A::CMatrix{Float64}, s::CVector{Float64})::JLWStatus + try + av = _copyin_vec(a) + bv = _copyin_vec(b) + M = _copyin_mat(A) + sv = Vector{Float64}(undef, size(M, 2)) + gramcover!(sv, av, bv, M) + length(s) == length(sv) || return jlw_error(4, "output length must match the number of columns of A") + copyto!(s, sv) + return jlw_ok() + catch e + return @status_from_exception(e) + end +end + +Base.@ccallable function mc_gramcover_weighted(a::CVector{Float64}, b::CVector{Float64}, + A::CMatrix{Float64}, w::CVector{Float64}, + s::CVector{Float64})::JLWStatus + try + av = _copyin_vec(a) + bv = _copyin_vec(b) + M = _copyin_mat(A) + wv = _copyin_vec(w) + sv = Vector{Float64}(undef, size(M, 2)) + gramcover!(sv, av, bv, M, wv) + length(s) == length(sv) || return jlw_error(4, "output length must match the number of columns of A") + copyto!(s, sv) + return jlw_ok() + catch e + return @status_from_exception(e) + end +end + +Base.@ccallable function mc_gramcover_matrix(a::CVector{Float64}, b::CVector{Float64}, + A::CMatrix{Float64}, W::CMatrix{Float64}, + s::CVector{Float64})::JLWStatus + try + av = _copyin_vec(a) + bv = _copyin_vec(b) + M = _copyin_mat(A) + Wm = _copyin_mat(W) + sv = Vector{Float64}(undef, size(M, 2)) + gramcover!(sv, av, bv, M, Wm) + length(s) == length(sv) || return jlw_error(4, "output length must match the number of columns of A") + copyto!(s, sv) + return jlw_ok() + catch e + return @status_from_exception(e) + end +end + +end # module diff --git a/lib/test/python/test_smoke.py b/lib/test/python/test_smoke.py new file mode 100644 index 0000000..abff63e --- /dev/null +++ b/lib/test/python/test_smoke.py @@ -0,0 +1,107 @@ +"""Smoke test for the bundled `matrixcovers` wheel.""" +import sys +import numpy as np +import matrixcovers as mc + +failures = [] + + +def check(name, cond): + if cond: + print(f"ok: {name}") + else: + print(f"FAIL: {name}") + failures.append(name) + + +A_sym = np.array([[4.0, 1.0], [1.0, 4.0]]) +A_asym = np.array([[1.0, 2.0, 3.0], [6.0, 5.0, 4.0]]) + +# symcover / iscover +a = mc.symcover(A_sym) +check("symcover matches reference", np.allclose(a, [2.0, 2.0], rtol=1e-8)) +check("symcover result covers A", mc.iscover(a, A_sym, rtol=1e-8)) + +# symcover_min +a_min = mc.symcover_min(A_sym) +check("symcover_min matches reference", np.allclose(a_min, [2.0, 2.0], rtol=1e-6)) + +# cover / cover_min +a_asym, b_asym = mc.cover(A_asym) +check( + "cover matches reference", + np.allclose(a_asym, [1.2544610775677627, 3.475905976749231], rtol=1e-6) + and np.allclose(b_asym, [1.7261686708831454, 1.621762761307448, 2.3914651906272066], rtol=1e-6), +) + +a_cmin, b_cmin = mc.cover_min(A_asym) +check( + "cover_min matches reference", + np.allclose(a_cmin, [1.1986299952850965, 3.2535823504068366], rtol=1e-6) + and np.allclose(b_cmin, [1.8441211421157804, 1.6685716376905815, 2.502857440802694], rtol=1e-6), +) + +# soft_symcover / soft_cover +a_soft = mc.soft_symcover(A_sym, penalty="abslog2") +check("soft_symcover(abslog2) matches reference", np.allclose(a_soft, [1.414213562373095, 1.414213562373095], rtol=1e-6)) + +# The default seed makes the multistart result reproducible. +a_soft_default = mc.soft_symcover(A_sym) +check("soft_symcover default matches reference", + np.allclose(a_soft_default, [1.8439088914584778, 1.8439088914585735], rtol=1e-6)) + +a_soft_l1 = mc.soft_symcover(A_sym, penalty="abslinear1") +check("soft_symcover(abslinear1) matches reference", np.allclose(a_soft_l1, [2.0, 2.0], rtol=1e-6)) + +a_sc, b_sc = mc.soft_cover(A_asym) +check( + "soft_cover default matches reference", + np.allclose(a_sc, [1.0817791286952234, 2.7823842907260183], rtol=1e-6) + and np.allclose(b_sc, [1.7867556808278862, 1.8232808727324803, 2.3172250172537487], rtol=1e-6), +) + +# cover_objective +for penalty in ("abslog1", "abslog2", "abslinear1", "abslinear2"): + val = mc.cover_objective(a, A_sym, penalty=penalty) + check(f"cover_objective({penalty}) is finite", np.isfinite(val)) + +# gramcover +s = mc.gramcover(a_asym, b_asym, A_asym) +G = A_asym.T @ A_asym +check("gramcover covers A'*A", np.all(s[:, None] * s[None, :] >= np.abs(G) - 1e-8)) + +# error paths +try: + mc.symcover_min(A_sym, penalty="not-a-penalty") + check("unknown penalty string raises ValueError", False) +except ValueError: + check("unknown penalty string raises ValueError", True) + +try: + from matrixcovers import _lowlevel + _lowlevel.mc_symcover_min( + 99, + _lowlevel.CMatrix_Float64.from_numpy(np.asfortranarray(A_sym)), + -1, 1, + _lowlevel.CVector_Float64.from_numpy(np.zeros(2)), + ) + check("bad low-level penalty enum raises JLWError", False) +except mc.JLWError as e: + check("bad low-level penalty enum raises JLWError", e.code == 1) + +try: + mc.iscover(np.zeros(5), A_sym) + check("shape-mismatched call raises cleanly", False) +except mc.JLWError as e: + check("shape-mismatched call raises cleanly", e.code == 4) + +try: + mc.symcover_min(A_sym, penalty="abslinear2") + check("extension-only penalty raises JLWError code 2", False) +except mc.JLWError as e: + check("extension-only penalty raises JLWError code 2", e.code == 2) + +if failures: + print(f"\n{len(failures)} check(s) failed: {failures}") + sys.exit(1) +print("\nALL CHECKS PASSED") From 76cb2307ccc5b90b213a97434eb425338c3d34b1 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 19 Aug 2026 06:32:08 -0500 Subject: [PATCH 2/2] Build Python wheels in release CI Use JuliaLibWrapping's reusable workflow to build wheels for releases and manual runs. Document installation and basic Python usage. Assisted-by: Claude Sonnet 5 --- .github/workflows/python-wheel.yml | 16 ++++++++++++++++ README.md | 25 +++++++++++++++++++++++++ lib/build-env/Project.toml | 2 +- 3 files changed, 42 insertions(+), 1 deletion(-) create mode 100644 .github/workflows/python-wheel.yml diff --git a/.github/workflows/python-wheel.yml b/.github/workflows/python-wheel.yml new file mode 100644 index 0000000..9c9adbc --- /dev/null +++ b/.github/workflows/python-wheel.yml @@ -0,0 +1,16 @@ +# Pin the `@main` ref to a JuliaLibWrapping release tag once one containing +# build-wrappers.yml exists. +name: Python wheel +on: + release: + types: [published] + workflow_dispatch: +permissions: + contents: write +jobs: + wheel: + uses: JuliaInterop/JuliaLibWrapping.jl/.github/workflows/build-wrappers.yml@main + with: + build-dir: lib + cpu-target: "generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)" + smoke-test: lib/test/python/test_smoke.py diff --git a/README.md b/README.md index 77f4fb3..b8d13a6 100644 --- a/README.md +++ b/README.md @@ -74,3 +74,28 @@ true See the [documentation](https://HolyLab.github.io/MatrixCovers.jl/dev/) for motivation, examples, and a full API reference. + +## Python + +A subset of MatrixCovers is available as a compiled Python package that does +not require Julia. Install a wheel from a +[GitHub release](https://github.com/HolyLab/MatrixCovers.jl/releases): + +``` +pip install https://github.com/HolyLab/MatrixCovers.jl/releases/download/vX.Y.Z/matrixcovers-X.Y.Z-py3-none-manylinux_2_35_x86_64.whl +``` + +Replace `X.Y.Z` with a released version. Wheels support Linux x86_64 with +glibc >= 2.35 (Ubuntu 22.04+, +Debian 12+, Fedora 36+; not RHEL/Rocky 9). + +```python +import numpy as np +import matrixcovers as mc + +A = np.array([[4.0, 2.0], [2.0, 16.0]]) +a = mc.symcover(A) # a[i] * a[j] >= abs(A[i, j]) +print(mc.iscover(a, A)) # True +``` + +See `lib/python/_facade.py` for the Python API. diff --git a/lib/build-env/Project.toml b/lib/build-env/Project.toml index 47a4a6c..26b0f62 100644 --- a/lib/build-env/Project.toml +++ b/lib/build-env/Project.toml @@ -9,5 +9,5 @@ JuliaLibWrapping = "d61f35a8-f6af-436f-bc10-cee6b101f7bd" [compat] JuliaC = "0.3" -JuliaLibWrapping = "0.2" +JuliaLibWrapping = "0.1.2" julia = "1.13"