From dd74b979da8e85197b1841a3c87c4430a39b25b0 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 19 Aug 2026 06:32:08 -0500 Subject: [PATCH 1/3] 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 4a70e7f95996aefc5dd0eb25b8509a3fd065d699 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 19 Aug 2026 06:32:08 -0500 Subject: [PATCH 2/3] 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 af6dad9..c9683e9 100644 --- a/README.md +++ b/README.md @@ -70,3 +70,28 @@ true See the [documentation](https://HolyLab.github.io/MatrixCovers.jl/dev/) for the algorithm guide and 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" From a3f96c5156b8cc9b2dc9f27343c6bf7112bbb2bf Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 1 Sep 2026 09:40:34 -0500 Subject: [PATCH 3/3] Update to JLW 0.2 --- lib/Project.toml | 2 +- lib/build-env/Project.toml | 2 +- lib/build.jl | 3 + lib/python/_facade.py | 474 +++++++++++++------------ lib/src/matrixcovers.jl | 626 ++++++++++++---------------------- lib/src/trimmability.jl | 165 +++++++++ lib/test/python/test_smoke.py | 32 +- 7 files changed, 652 insertions(+), 652 deletions(-) create mode 100644 lib/src/trimmability.jl diff --git a/lib/Project.toml b/lib/Project.toml index 982870f..86b3df0 100644 --- a/lib/Project.toml +++ b/lib/Project.toml @@ -10,6 +10,6 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" # `build.jl` adds the local MatrixCovers source to a temporary copy. [compat] -JLWInterop = "0.1" +JLWInterop = "0.2" MatrixCovers = "1" julia = "1.13" diff --git a/lib/build-env/Project.toml b/lib/build-env/Project.toml index 26b0f62..47a4a6c 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.1.2" +JuliaLibWrapping = "0.2" julia = "1.13" diff --git a/lib/build.jl b/lib/build.jl index 33627b4..992ddc2 100644 --- a/lib/build.jl +++ b/lib/build.jl @@ -17,6 +17,9 @@ function prepare_project() toml = TOML.parsefile(joinpath(HERE, "Project.toml")) sources = get(toml, "sources", Dict{String, Any}()) sources["MatrixCovers"] = Dict("path" => REPO_ROOT) + # Developer override for building against an unreleased JLWInterop. + jlwinterop_path = get(ENV, "MATRIXCOVERS_JLWINTEROP_PATH", "") + isempty(jlwinterop_path) || (sources["JLWInterop"] = Dict("path" => jlwinterop_path)) toml["sources"] = sources tmp = mktempdir(; prefix = "matrixcovers-lib-project-") open(joinpath(tmp, "Project.toml"), "w") do io diff --git a/lib/python/_facade.py b/lib/python/_facade.py index 2b7c72f..635adee 100644 --- a/lib/python/_facade.py +++ b/lib/python/_facade.py @@ -1,229 +1,276 @@ """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 +Matrix and vector arguments are converted to Fortran-order/contiguous +`float64` arrays; outputs are newly allocated `numpy.ndarray`s. `Penalty` and +`Linsolve` are enum classes: a keyword typed against one accepts a member +(`Penalty.abslog2`), its name as a string (`"abslog2"`), or the underlying +int. An unrecognized name or value raises `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, -} +from . import _lowlevel # noqa: F401 +import numpy as np # noqa: F401 + +from ._lowlevel import ( + JLWStatus, + JLWResult_Bool, + COpt_Int64, + CString_owned, + CVector_owned_Float64, + JLWResult_CVector_owned_Float64, + CVector_borrowed_Float64, + COpt_Float64, + JLWResult_Float64, + CMatrix_borrowed_Float64, + Linsolve, + Penalty, + JLWError, + _enum_coerce, +) -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 +def _as_matrix(A): + return np.asfortranarray(A, dtype=np.float64) -_LINSOLVE_CODES = {"auto": 1, "dense": 2, "lsqr": 3} +def _as_vector(v): + return np.ascontiguousarray(v, dtype=np.float64) -def _linsolve_code(linsolve): +def soft_cover_ab(A, *, penalty=Penalty.abslinear2, maxiter=None, starts=None, sigma=None, seed=0): + """Asymmetric soft cover of `A` minimizing the penalty, with no coverage +constraint, as `vcat(a, b)`; `a` has length `size(A, 1)`. `abslog1` and +`abslog2` accept only `maxiter`; `abslinear1` and `abslinear2` additionally +accept `starts`, `sigma`, and `seed`.""" + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _penalty = _enum_coerce(Penalty, penalty) + _maxiter = COpt_Int64.from_optional(maxiter) + _starts = COpt_Int64.from_optional(starts) + _sigma = COpt_Float64.from_optional(sigma) + _r = _lowlevel.matrixcovers_soft_cover_ab(_A, _penalty, _maxiter, _starts, _sigma, seed) 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) + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def iscover_sym(a, A, *, rtol=0.0, atol=0.0): + """Whether `a` (playing the role of both `a` and `b`) covers `A`: `a[i]*a[j] >= abs(A[i, j])`.""" + _a = CVector_borrowed_Float64.from_numpy(_as_vector(a)) + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _r = _lowlevel.matrixcovers_iscover_sym(_a, _A, rtol, atol) + return _r.value + +def soft_symcover_min(A, *, penalty=Penalty.abslog2, maxiter=None): + """ϕ-minimal symmetric soft cover of `A`, with no coverage constraint. Only +`abslog2` is solved natively; `MatrixCovers` does not implement `abslog1` for +this entrypoint, and the `abslinear` penalties need the JuMP/Ipopt extension.""" + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _penalty = _enum_coerce(Penalty, penalty) + _maxiter = COpt_Int64.from_optional(maxiter) + _r = _lowlevel.matrixcovers_soft_symcover_min(_A, _penalty, _maxiter) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def gramcover_weighted(a, b, A, w): + """Symmetric cover of the weighted Gram matrix `A'*Diagonal(w)*A`, from a cover `(a, b)` of `A`.""" + _a = CVector_borrowed_Float64.from_numpy(_as_vector(a)) + _b = CVector_borrowed_Float64.from_numpy(_as_vector(b)) + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _w = CVector_borrowed_Float64.from_numpy(_as_vector(w)) + _r = _lowlevel.matrixcovers_gramcover_weighted(_a, _b, _A, _w) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def soft_cover_min_ab(A, *, penalty=Penalty.abslog2, maxiter=None): + """ϕ-minimal asymmetric soft cover of `A`, with no coverage constraint, as +`vcat(a, b)`; `a` has length `size(A, 1)`. Only `abslog2` is solved natively; +`MatrixCovers` does not implement `abslog1` for this entrypoint, and the +`abslinear` penalties need the JuMP/Ipopt extension.""" + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _penalty = _enum_coerce(Penalty, penalty) + _maxiter = COpt_Int64.from_optional(maxiter) + _r = _lowlevel.matrixcovers_soft_cover_min_ab(_A, _penalty, _maxiter) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def symcover_min(A, *, penalty=Penalty.abslog2, maxiter=None, linsolve=Linsolve.auto): + """ϕ-minimal symmetric hard cover of `A`. Only `abslog2` is solved natively; +`abslog1` needs the JuMP/HiGHS extension and the `abslinear` penalties need the +JuMP/Ipopt extension, neither linked into this library.""" + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _penalty = _enum_coerce(Penalty, penalty) + _maxiter = COpt_Int64.from_optional(maxiter) + _linsolve = _enum_coerce(Linsolve, linsolve) + _r = _lowlevel.matrixcovers_symcover_min(_A, _penalty, _maxiter, _linsolve) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def cover_objective_ab(a, b, A, *, penalty=Penalty.abslog2): + """`sum(penalty(abs(A[i,j]) / (a[i]*b[j])))`.""" + _a = CVector_borrowed_Float64.from_numpy(_as_vector(a)) + _b = CVector_borrowed_Float64.from_numpy(_as_vector(b)) + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _penalty = _enum_coerce(Penalty, penalty) + _r = _lowlevel.matrixcovers_cover_objective_ab(_a, _b, _A, _penalty) + return _r.value + +def cover_min_ab(A, *, penalty=Penalty.abslog2, maxiter=None, linsolve=Linsolve.auto): + """ϕ-minimal hard cover of `A`, as `vcat(a, b)`; `a` has length `size(A, 1)`. Only +`abslog2` is solved natively; `abslog1` needs the JuMP/HiGHS extension and the +`abslinear` penalties need the JuMP/Ipopt extension, neither linked into this +library.""" + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _penalty = _enum_coerce(Penalty, penalty) + _maxiter = COpt_Int64.from_optional(maxiter) + _linsolve = _enum_coerce(Linsolve, linsolve) + _r = _lowlevel.matrixcovers_cover_min_ab(_A, _penalty, _maxiter, _linsolve) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def soft_symcover(A, *, penalty=Penalty.abslinear2, maxiter=None, starts=None, sigma=None, seed=0): + """Symmetric soft cover of `A` minimizing the penalty, with no coverage +constraint. `abslog1` and `abslog2` accept only `maxiter`; `abslinear1` and +`abslinear2` additionally accept `starts`, `sigma`, and `seed` (which seeds the +multistart perturbation stream).""" + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _penalty = _enum_coerce(Penalty, penalty) + _maxiter = COpt_Int64.from_optional(maxiter) + _starts = COpt_Int64.from_optional(starts) + _sigma = COpt_Float64.from_optional(sigma) + _r = _lowlevel.matrixcovers_soft_symcover(_A, _penalty, _maxiter, _starts, _sigma, seed) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def gramcover_matrix(a, b, A, W): + """Symmetric cover of the weighted Gram matrix `A'*W*A`, from a cover `(a, b)` of `A`.""" + _a = CVector_borrowed_Float64.from_numpy(_as_vector(a)) + _b = CVector_borrowed_Float64.from_numpy(_as_vector(b)) + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _W = CMatrix_borrowed_Float64.from_numpy(_as_matrix(W)) + _r = _lowlevel.matrixcovers_gramcover_matrix(_a, _b, _A, _W) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def iscover_ab(a, b, A, *, rtol=0.0, atol=0.0): + """Whether `a`, `b` cover `A`: `a[i]*b[j] >= abs(A[i, j])`.""" + _a = CVector_borrowed_Float64.from_numpy(_as_vector(a)) + _b = CVector_borrowed_Float64.from_numpy(_as_vector(b)) + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _r = _lowlevel.matrixcovers_iscover_ab(_a, _b, _A, rtol, atol) + return _r.value + +def cover_objective_sym(a, A, *, penalty=Penalty.abslog2): + """`sum(penalty(abs(A[i,j]) / (a[i]*a[j])))`.""" + _a = CVector_borrowed_Float64.from_numpy(_as_vector(a)) + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _penalty = _enum_coerce(Penalty, penalty) + _r = _lowlevel.matrixcovers_cover_objective_sym(_a, _A, _penalty) + return _r.value +def symcover(A, *, maxiter=None): + """Heuristic symmetric hard cover: `a` with `a[i]*a[j] >= abs(A[i, j])`.""" + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _maxiter = COpt_Int64.from_optional(maxiter) + _r = _lowlevel.matrixcovers_symcover(_A, _maxiter) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def cover_ab(A, *, maxiter=None): + """Heuristic hard cover of `A`, as `vcat(a, b)`; `a` has length `size(A, 1)`.""" + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _maxiter = COpt_Int64.from_optional(maxiter) + _r = _lowlevel.matrixcovers_cover_ab(_A, _maxiter) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out + +def _gramcover_plain(a, b, A): + """Symmetric cover of the Gram matrix `A'*A`, from a cover `(a, b)` of `A`.""" + _a = CVector_borrowed_Float64.from_numpy(_as_vector(a)) + _b = CVector_borrowed_Float64.from_numpy(_as_vector(b)) + _A = CMatrix_borrowed_Float64.from_numpy(_as_matrix(A)) + _r = _lowlevel.matrixcovers_gramcover(_a, _b, _A) + try: + _out = np.array(_r.value.as_numpy(), copy=True) + finally: + _r.value.free() + return _out -def _sentinel_float(value): - return -1.0 if value is None else float(value) +# Hand-written public API. Each of these splits or dispatches across the +# generated entrypoints above, whose own names gain a leading underscore +# where reused (`gramcover` -> `_gramcover_plain`) to free the public name. -def _as_matrix(A): - return np.asfortranarray(A, dtype=np.float64) +def cover(A, *, maxiter=None): + """Heuristic hard cover of `A`: `(a, b)` with `a[i]*b[j] >= abs(A[i, j])`.""" + out = cover_ab(A, maxiter=maxiter) + m = np.shape(A)[0] + return out[:m], out[m:] -def _as_vector(v): - return np.ascontiguousarray(v, dtype=np.float64) +def cover_min(A, *, penalty=Penalty.abslog2, maxiter=None, linsolve=Linsolve.auto): + """ϕ-minimal hard cover of `A`: `(a, b)`. Only `abslog2` is solved natively.""" + out = cover_min_ab(A, penalty=penalty, maxiter=maxiter, linsolve=linsolve) + m = np.shape(A)[0] + return out[:m], out[m:] -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 soft_cover(A, *, penalty=Penalty.abslinear2, maxiter=None, starts=None, sigma=None, seed=0): + """Asymmetric soft cover of `A` minimizing the penalty, with no coverage constraint.""" + out = soft_cover_ab(A, penalty=penalty, maxiter=maxiter, starts=starts, sigma=sigma, seed=seed) + m = np.shape(A)[0] + return out[:m], out[m:] -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 soft_cover_min(A, *, penalty=Penalty.abslog2, maxiter=None): + """ϕ-minimal asymmetric soft cover of `A`, with no coverage constraint.""" + out = soft_cover_min_ab(A, penalty=penalty, maxiter=maxiter) + m = np.shape(A)[0] + return out[:m], out[m:] 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`. + """Whether `a`, `b` cover `A`: `a[i]*b[j] >= abs(A[i, j])`. `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) + return iscover_sym(a, A, rtol=rtol, atol=atol) + return iscover_ab(a, b, A, rtol=rtol, atol=atol) + + +def cover_objective(a, A, b=None, *, penalty=Penalty.abslog2): + """`sum(penalty(abs(A[i,j]) / (a[i]*b[j])))`; `b=None` tests the symmetric cover `a*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) + return cover_objective_sym(a, A, penalty=penalty) + return cover_objective_ab(a, b, A, penalty=penalty) def gramcover(a, b, A, *, w=None, W=None): @@ -235,40 +282,15 @@ def gramcover(a, b, A, *, w=None, W=None): """ 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 + return gramcover_weighted(a, b, A, w) + if W is not None: + return gramcover_matrix(a, b, A, W) + return _gramcover_plain(a, b, A) __all__ = [ - "JLWError", + "JLWError", "Penalty", "Linsolve", "symcover", "cover", "symcover_min", "cover_min", "soft_symcover", "soft_cover", diff --git a/lib/src/matrixcovers.jl b/lib/src/matrixcovers.jl index f4d5bfe..e24b793 100644 --- a/lib/src/matrixcovers.jl +++ b/lib/src/matrixcovers.jl @@ -1,455 +1,247 @@ -# 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 +# The binding layer: one `@api` entrypoint per exposed `MatrixCovers` call, +# each with the docstring the generated Python and C interfaces carry. Calls +# `juliac --trim=safe` cannot resolve in their natural form are reshaped in +# `trimmability.jl`; entrypoints that need such a call delegate to it. + +""" + matrixcovers + +The binding layer for [`MatrixCovers`](@ref). `@api` declares which of the +package's cover solvers, objective, and predicate a foreign caller may reach; +[`Penalty`](@ref) and [`Linsolve`](@ref) select the penalty and inner linear +solve for the entrypoints that need them. + +Every `Matrix{Float64}`/`Vector{Float64}` argument arrives as a zero-copy view +of the caller's memory. None of the wrapped `MatrixCovers` functions mutate a +caller-supplied array, so no entrypoint copies its input; a hard-cover +entrypoint that needs scratch storage for its result allocates that storage +itself. +""" module matrixcovers using JLWInterop -using MatrixCovers +using MatrixCovers: MatrixCovers, AbsLog, AbsLinear using Random: MersenneTwister -# Results that contain both a status and a value. +@export_release_entrypoints -struct MCBoolResult - status::JLWStatus - value::Int32 -end +""" + Penalty -struct MCScalarResult - status::JLWStatus - value::Float64 -end +Which cover penalty an entrypoint scores or minimizes against: `abslog1` +(`AbsLog{1}`), `abslog2` (`AbsLog{2}`), `abslinear1` (`AbsLinear{1}`), or +`abslinear2` (`AbsLinear{2}`). +""" +@enum Penalty::Int32 abslog1 = 1 abslog2 = 2 abslinear1 = 3 abslinear2 = 4 -# Boundary helpers +""" + Linsolve -function _copyin_vec(v::CVector{Float64}) - out = Vector{Float64}(undef, length(v)) - copyto!(out, v) - return out -end +The inner linear solve [`symcover_min`](@ref)/[`cover_min`](@ref) uses for +`AbsLog{2}`: `auto` (dense, for a dense matrix), `dense`, or `lsqr` +(matrix-free, intended for large sparse supports). +""" +@enum Linsolve::Int32 auto = 1 dense = 2 lsqr = 3 -function _copyin_mat(A::CMatrix{Float64}) - m, n = size(A) - out = Matrix{Float64}(undef, m, n) - copyto!(out, A) - return out -end +const _EXT_JUMP = "penalty AbsLog{1} requires the MatrixCoversJuMPExt extension (JuMP and HiGHS)" +const _EXT_IPOPT = "penalty AbsLinear requires the MatrixCoversIpoptExt extension (JuMP and Ipopt)" -# 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 +_linsolve_symbol(ls::Linsolve) = ls === auto ? :auto : ls === dense ? :dense : :lsqr -# 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 +include("trimmability.jl") -# soft_symcover / soft_cover kwarg plumbing -# -# Penalty types and explicit keyword combinations keep calls trim-safe. +# --- Hard covers ------------------------------------------------------------ -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 +"Heuristic symmetric hard cover: `a` with `a[i]*a[j] >= abs(A[i, j])`." +symcover(A::Matrix{Float64}; maxiter::Union{Int64, Nothing} = nothing) = _symcover(A, maxiter) -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 +@api symcover(A::Matrix{Float64}; maxiter::Union{Int64, Nothing} = nothing)::Vector{Float64} -# 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 +"Heuristic hard cover of `A`, as `vcat(a, b)`; `a` has length `size(A, 1)`." +cover_ab(A::Matrix{Float64}; maxiter::Union{Int64, Nothing} = nothing) = _cover_ab(A, maxiter) -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 +@api cover_ab(A::Matrix{Float64}; maxiter::Union{Int64, Nothing} = nothing)::Vector{Float64} -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 +""" +ϕ-minimal symmetric hard cover of `A`. Only `abslog2` is solved natively; +`abslog1` needs the JuMP/HiGHS extension and the `abslinear` penalties need the +JuMP/Ipopt extension, neither linked into this library. +""" +function symcover_min(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing, linsolve::Linsolve = auto) + penalty === abslog2 || throw(ArgumentError(penalty === abslog1 ? _EXT_JUMP : _EXT_IPOPT)) + ls = _linsolve_symbol(linsolve) + return isnothing(maxiter) ? MatrixCovers.symcover_min(AbsLog{2}(), A; linsolve = ls) : + MatrixCovers.symcover_min(AbsLog{2}(), A; maxiter = Int(maxiter), linsolve = ls) 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 +@api symcover_min(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing, linsolve::Linsolve = auto)::Vector{Float64} -# 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 +""" +ϕ-minimal hard cover of `A`, as `vcat(a, b)`; `a` has length `size(A, 1)`. Only +`abslog2` is solved natively; `abslog1` needs the JuMP/HiGHS extension and the +`abslinear` penalties need the JuMP/Ipopt extension, neither linked into this +library. +""" +function cover_min_ab(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing, linsolve::Linsolve = auto) + penalty === abslog2 || throw(ArgumentError(penalty === abslog1 ? _EXT_JUMP : _EXT_IPOPT)) + ls = _linsolve_symbol(linsolve) + a, b = isnothing(maxiter) ? MatrixCovers.cover_min(AbsLog{2}(), A; linsolve = ls) : + MatrixCovers.cover_min(AbsLog{2}(), A; maxiter = Int(maxiter), linsolve = ls) + return vcat(a, b) 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 +@api cover_min_ab(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing, linsolve::Linsolve = auto)::Vector{Float64} -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 +# --- Soft covers ------------------------------------------------------------- +# +# The soft-cover entrypoints resolve each optional keyword to a concrete +# sentinel (`-1`/`NaN` for "omitted") and do nothing else: the branch on +# `penalty` and the sentinels lives in the `@noinline` dispatchers in +# `trimmability.jl`. + +""" +Symmetric soft cover of `A` minimizing the penalty, with no coverage +constraint. `abslog1` and `abslog2` accept only `maxiter`; `abslinear1` and +`abslinear2` additionally accept `starts`, `sigma`, and `seed` (which seeds the +multistart perturbation stream). +""" +function soft_symcover(A::Matrix{Float64}; penalty::Penalty = abslinear2, + maxiter::Union{Int64, Nothing} = nothing, + starts::Union{Int64, Nothing} = nothing, + sigma::Union{Float64, Nothing} = nothing, seed::Int64 = 0) + mi = isnothing(maxiter) ? -1 : Int(maxiter) + st = isnothing(starts) ? -1 : Int(starts) + sg = isnothing(sigma) ? NaN : Float64(sigma) + return _soft_symcover_dispatch(A, penalty, mi, st, sg, seed) +end + +@api soft_symcover(A::Matrix{Float64}; penalty::Penalty = abslinear2, + maxiter::Union{Int64, Nothing} = nothing, + starts::Union{Int64, Nothing} = nothing, + sigma::Union{Float64, Nothing} = nothing, seed::Int64 = 0)::Vector{Float64} + +""" +Asymmetric soft cover of `A` minimizing the penalty, with no coverage +constraint, as `vcat(a, b)`; `a` has length `size(A, 1)`. `abslog1` and +`abslog2` accept only `maxiter`; `abslinear1` and `abslinear2` additionally +accept `starts`, `sigma`, and `seed`. +""" +function soft_cover_ab(A::Matrix{Float64}; penalty::Penalty = abslinear2, + maxiter::Union{Int64, Nothing} = nothing, + starts::Union{Int64, Nothing} = nothing, + sigma::Union{Float64, Nothing} = nothing, seed::Int64 = 0) + mi = isnothing(maxiter) ? -1 : Int(maxiter) + st = isnothing(starts) ? -1 : Int(starts) + sg = isnothing(sigma) ? NaN : Float64(sigma) + return _soft_cover_dispatch(A, penalty, mi, st, sg, seed) +end + +@api soft_cover_ab(A::Matrix{Float64}; penalty::Penalty = abslinear2, + maxiter::Union{Int64, Nothing} = nothing, + starts::Union{Int64, Nothing} = nothing, + sigma::Union{Float64, Nothing} = nothing, seed::Int64 = 0)::Vector{Float64} + +""" +ϕ-minimal symmetric soft cover of `A`, with no coverage constraint. Only +`abslog2` is solved natively; `MatrixCovers` does not implement `abslog1` for +this entrypoint, and the `abslinear` penalties need the JuMP/Ipopt extension. +""" +function soft_symcover_min(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing) + mi = isnothing(maxiter) ? -1 : Int(maxiter) + return _soft_symcover_min_dispatch(A, penalty, mi) +end + +@api soft_symcover_min(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing)::Vector{Float64} + +""" +ϕ-minimal asymmetric soft cover of `A`, with no coverage constraint, as +`vcat(a, b)`; `a` has length `size(A, 1)`. Only `abslog2` is solved natively; +`MatrixCovers` does not implement `abslog1` for this entrypoint, and the +`abslinear` penalties need the JuMP/Ipopt extension. +""" +function soft_cover_min_ab(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing) + mi = isnothing(maxiter) ? -1 : Int(maxiter) + return _soft_cover_min_dispatch(A, penalty, mi) +end + +@api soft_cover_min_ab(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing)::Vector{Float64} + +# --- Predicates and objectives ----------------------------------------------- + +"Whether `a` (playing the role of both `a` and `b`) covers `A`: `a[i]*a[j] >= abs(A[i, j])`." +iscover_sym(a::Vector{Float64}, A::Matrix{Float64}; rtol::Float64 = 0.0, atol::Float64 = 0.0) = + MatrixCovers.iscover(a, A; rtol, atol) + +@api iscover_sym(a::Vector{Float64}, A::Matrix{Float64}; rtol::Float64 = 0.0, atol::Float64 = 0.0)::Bool + +"Whether `a`, `b` cover `A`: `a[i]*b[j] >= abs(A[i, j])`." +iscover_ab(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}; rtol::Float64 = 0.0, atol::Float64 = 0.0) = + MatrixCovers.iscover(a, b, A; rtol, atol) + +@api iscover_ab(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}; + rtol::Float64 = 0.0, atol::Float64 = 0.0)::Bool + +_cover_objective(::Type{P}, a::Vector{Float64}, A::Matrix{Float64}) where {P <: MatrixCovers.AbstractCoverPenalty} = + MatrixCovers.cover_objective(P(), a, A) +_cover_objective(::Type{P}, a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}) where {P <: MatrixCovers.AbstractCoverPenalty} = + MatrixCovers.cover_objective(P(), a, b, A) + +"`sum(penalty(abs(A[i,j]) / (a[i]*a[j])))`." +function cover_objective_sym(a::Vector{Float64}, A::Matrix{Float64}; penalty::Penalty = abslog2) + penalty === abslog1 && return _cover_objective(AbsLog{1}, a, A) + penalty === abslog2 && return _cover_objective(AbsLog{2}, a, A) + penalty === abslinear1 && return _cover_objective(AbsLinear{1}, a, A) + return _cover_objective(AbsLinear{2}, a, A) +end + +@api cover_objective_sym(a::Vector{Float64}, A::Matrix{Float64}; penalty::Penalty = abslog2)::Float64 + +"`sum(penalty(abs(A[i,j]) / (a[i]*b[j])))`." +function cover_objective_ab(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}; penalty::Penalty = abslog2) + penalty === abslog1 && return _cover_objective(AbsLog{1}, a, b, A) + penalty === abslog2 && return _cover_objective(AbsLog{2}, a, b, A) + penalty === abslinear1 && return _cover_objective(AbsLinear{1}, a, b, A) + return _cover_objective(AbsLinear{2}, a, b, A) 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 +@api cover_objective_ab(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}; + penalty::Penalty = abslog2)::Float64 -# 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 +# --- Gram covers -------------------------------------------------------------- -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 +"Symmetric cover of the Gram matrix `A'*A`, from a cover `(a, b)` of `A`." +function gramcover(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}) + s = Vector{Float64}(undef, size(A, 2)) + MatrixCovers.gramcover!(s, a, b, A) + return s 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 +@api gramcover(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64})::Vector{Float64} -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 +"Symmetric cover of the weighted Gram matrix `A'*Diagonal(w)*A`, from a cover `(a, b)` of `A`." +function gramcover_weighted(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}, w::Vector{Float64}) + s = Vector{Float64}(undef, size(A, 2)) + MatrixCovers.gramcover!(s, a, b, A, w) + return s 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 +@api gramcover_weighted(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}, w::Vector{Float64})::Vector{Float64} -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 +"Symmetric cover of the weighted Gram matrix `A'*W*A`, from a cover `(a, b)` of `A`." +function gramcover_matrix(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}, W::Matrix{Float64}) + s = Vector{Float64}(undef, size(A, 2)) + MatrixCovers.gramcover!(s, a, b, A, W) + return s 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 +@api gramcover_matrix(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}, W::Matrix{Float64})::Vector{Float64} end # module diff --git a/lib/src/trimmability.jl b/lib/src/trimmability.jl new file mode 100644 index 0000000..c86dd45 --- /dev/null +++ b/lib/src/trimmability.jl @@ -0,0 +1,165 @@ +# Workarounds for `juliac --trim=safe`, which compiles only calls it can +# resolve statically. A few `MatrixCovers` calls are not resolvable in their +# natural form from an `@api` entrypoint; each function here reshapes one of +# them — allocating an output so a mutating method can be called, resolving an +# optional keyword to a concrete sentinel, or isolating a branch behind +# `@noinline` — and carries the constraint that forces the shape. The +# entrypoints themselves are in `matrixcovers.jl`. + +# --- Hard covers ------------------------------------------------------------ + +# `MatrixCovers.symcover` allocates its own output and forwards `maxiter` +# through two layers of `kwargs...` before reaching a concretely-typed callee +# (`tighten_cover!`); `--trim=safe` cannot statically resolve that chain. The +# mutating `symcover!` removes one layer, which is enough to make the call +# resolvable, so the output is allocated here and that method called instead. +function _symcover(A::Matrix{Float64}, maxiter::Union{Int64, Nothing}) + a = Vector{Float64}(undef, size(A, 1)) + return isnothing(maxiter) ? MatrixCovers.symcover!(a, A) : MatrixCovers.symcover!(a, A; maxiter = Int(maxiter)) +end + +# `cover!`'s own `tighten_cover!(a, b, A; kwargs...)` call is not in tail +# position (it balances and inflates the result afterward), and `--trim=safe` +# cannot resolve a `kwargs...`-forwarded call there, even through the mutating +# form. `cover!`'s five-step body is reproduced here instead, with an explicit +# (non-forwarded) `maxiter` on the one step that takes it. +function _cover_ab(A::Matrix{Float64}, maxiter::Union{Int64, Nothing}) + a = Vector{Float64}(undef, size(A, 1)) + b = Vector{Float64}(undef, size(A, 2)) + MatrixCovers.unconstrained_min!(AbsLog{2}(), a, b, A) + MatrixCovers.boost_feasible!(a, b, A) + isnothing(maxiter) ? MatrixCovers.tighten_cover!(a, b, A) : + MatrixCovers.tighten_cover!(a, b, A; maxiter = Int(maxiter)) + MatrixCovers._balance_cover!(a, b, A) + MatrixCovers.inflate_feasible!(a, b, A) + return vcat(a, b) +end + +# --- Soft covers ------------------------------------------------------------- +# +# `MatrixCovers.soft_symcover_min`/`soft_cover_min` (`AbsLog{2}`) reduce to +# `_symcover_min_abslog2`/`_cover_min_abslog2` (the native hard-cover kernel) +# with `κs = ()`, `boost = false` — no penalty continuation, no boost to +# feasibility, giving the unconstrained minimum with no coverage constraint. +# That kernel is not exported, but calling it directly removes two layers of +# `kwargs...` forwarding `--trim=safe` cannot resolve. Three further +# conditions are needed for the call to resolve when reached (even +# transitively, even through further `@noinline` layers) from a function +# whose own signature carries an optional (`Union{T,Nothing}`) keyword — as +# every `@api` entrypoint does: `maxiter` must always be supplied (a +# local copy of the kernel's own default, 40, standing in when the caller +# omits it, since keyword *presence* varying across call sites is itself +# unresolvable); `κs` must be a concretely and consistently typed empty +# container (`Float64[]`, not `()` — the untyped empty tuple forces a second, +# mismatched specialization of the kernel alongside the one its +# `NTuple{4,Float64}` default already requires); and each of `soft_symcover_min`'s +# reachable `@api` entrypoints must resolve the optional keyword to a +# concrete sentinel (`-1` for "omitted") and delegate to its *own* `@noinline` +# dispatcher and kernel, never sharing either with another entrypoint — +# reachability from more than one `@api` root, or a shared kernel across +# entrypoints whose own keyword counts differ, reintroduces the failure. +# +# `soft_symcover`/`soft_cover` have three optional keywords (`maxiter`, +# `starts`, `sigma`), and `--trim=safe` cannot resolve a call into the +# `_soft_*_abslog2`/`AbsLinear` machinery from a function with *that +# many* optional keywords of its own, independent of which branch is taken and +# regardless of intervening `@noinline` layers — unlike `soft_symcover_min` +# (one optional keyword), which resolves the identical call fine. Every +# keyword is therefore resolved to a concrete sentinel value (`-1`/`NaN` for +# "omitted"; a real `maxiter`/`starts` is always non-negative, a real `sigma` +# always finite) in the entrypoint, and the branch on `penalty` and the +# sentinels is pulled into its own `@noinline` dispatcher — the entrypoint's +# own body touches nothing but that resolution, so no optional-keyword type +# reaches any call. + +@noinline function _soft_symcover_dispatch(A::Matrix{Float64}, penalty::Penalty, + mi::Int64, st::Int64, sg::Float64, seed::Int64) + if penalty === abslog1 + return mi < 0 ? MatrixCovers.soft_symcover(AbsLog{1}(), A) : + MatrixCovers.soft_symcover(AbsLog{1}(), A; maxiter = mi) + end + if penalty === abslog2 + a, _ = _soft_symcover_min_dispatch_kernel(A, mi < 0 ? 40 : mi) + return a + end + if penalty === abslinear1 + m = mi < 0 ? 20 : mi + s = st < 0 ? 5 : st + return isnan(sg) ? _soft_symcover_call(AbsLinear{1}(), A, m, s, seed) : + _soft_symcover_call(AbsLinear{1}(), A, m, s, sg, seed) + end + m = mi < 0 ? 32 : mi + s = st < 0 ? 5 : st + return isnan(sg) ? _soft_symcover_call(AbsLinear{2}(), A, m, s, seed) : + _soft_symcover_call(AbsLinear{2}(), A, m, s, sg, seed) +end + +@noinline _soft_symcover_min_dispatch_kernel(A::Matrix{Float64}, mi::Int64) = + MatrixCovers._symcover_min_abslog2(A; maxiter = mi, κs = Float64[], boost = false, fname = :soft_symcover) + +@noinline _soft_symcover_call(::AbsLinear{1}, A::Matrix{Float64}, mi::Int64, st::Int64, seed::Int64) = + MatrixCovers.soft_symcover(AbsLinear{1}(), A; maxiter = mi, starts = st, rng = MersenneTwister(seed)) +@noinline _soft_symcover_call(::AbsLinear{1}, A::Matrix{Float64}, mi::Int64, st::Int64, sg::Float64, seed::Int64) = + MatrixCovers.soft_symcover(AbsLinear{1}(), A; maxiter = mi, starts = st, sigma = sg, rng = MersenneTwister(seed)) +@noinline _soft_symcover_call(::AbsLinear{2}, A::Matrix{Float64}, mi::Int64, st::Int64, seed::Int64) = + MatrixCovers.soft_symcover(AbsLinear{2}(), A; maxiter = mi, starts = st, rng = MersenneTwister(seed)) +@noinline _soft_symcover_call(::AbsLinear{2}, A::Matrix{Float64}, mi::Int64, st::Int64, sg::Float64, seed::Int64) = + MatrixCovers.soft_symcover(AbsLinear{2}(), A; maxiter = mi, starts = st, sigma = sg, rng = MersenneTwister(seed)) + +@noinline function _soft_cover_dispatch(A::Matrix{Float64}, penalty::Penalty, + mi::Int64, st::Int64, sg::Float64, seed::Int64) + if penalty === abslog1 + a, b = mi < 0 ? MatrixCovers.soft_cover(AbsLog{1}(), A) : + MatrixCovers.soft_cover(AbsLog{1}(), A; maxiter = mi) + return vcat(a, b) + end + if penalty === abslog2 + a, b, _ = _soft_cover_min_dispatch_kernel(A, mi < 0 ? 40 : mi) + return vcat(a, b) + end + if penalty === abslinear1 + m = mi < 0 ? 100 : mi + s = st < 0 ? 4 : st + a, b = isnan(sg) ? _soft_cover_call(AbsLinear{1}(), A, m, s, seed) : + _soft_cover_call(AbsLinear{1}(), A, m, s, sg, seed) + return vcat(a, b) + end + m = mi < 0 ? 200 : mi + s = st < 0 ? 4 : st + a, b = isnan(sg) ? _soft_cover_call(AbsLinear{2}(), A, m, s, seed) : + _soft_cover_call(AbsLinear{2}(), A, m, s, sg, seed) + return vcat(a, b) +end + +@noinline _soft_cover_min_dispatch_kernel(A::Matrix{Float64}, mi::Int64) = + MatrixCovers._cover_min_abslog2(A; maxiter = mi, κs = Float64[], boost = false) + +@noinline _soft_cover_call(::AbsLinear{1}, A::Matrix{Float64}, mi::Int64, st::Int64, seed::Int64) = + MatrixCovers.soft_cover(AbsLinear{1}(), A; maxiter = mi, starts = st, rng = MersenneTwister(seed)) +@noinline _soft_cover_call(::AbsLinear{1}, A::Matrix{Float64}, mi::Int64, st::Int64, sg::Float64, seed::Int64) = + MatrixCovers.soft_cover(AbsLinear{1}(), A; maxiter = mi, starts = st, sigma = sg, rng = MersenneTwister(seed)) +@noinline _soft_cover_call(::AbsLinear{2}, A::Matrix{Float64}, mi::Int64, st::Int64, seed::Int64) = + MatrixCovers.soft_cover(AbsLinear{2}(), A; maxiter = mi, starts = st, rng = MersenneTwister(seed)) +@noinline _soft_cover_call(::AbsLinear{2}, A::Matrix{Float64}, mi::Int64, st::Int64, sg::Float64, seed::Int64) = + MatrixCovers.soft_cover(AbsLinear{2}(), A; maxiter = mi, starts = st, sigma = sg, rng = MersenneTwister(seed)) + +@noinline function _soft_symcover_min_dispatch(A::Matrix{Float64}, penalty::Penalty, mi::Int64) + penalty === abslog2 && return _soft_symcover_min_only_kernel(A, mi < 0 ? 40 : mi) + penalty === abslog1 && throw(ArgumentError("MatrixCovers does not implement AbsLog{1} for soft_symcover_min")) + throw(ArgumentError(_EXT_IPOPT)) +end + +@noinline _soft_symcover_min_only_kernel(A::Matrix{Float64}, mi::Int64) = + MatrixCovers._symcover_min_abslog2(A; maxiter = mi, κs = Float64[], boost = false, fname = :soft_symcover_min)[1] + +@noinline function _soft_cover_min_dispatch(A::Matrix{Float64}, penalty::Penalty, mi::Int64) + if penalty === abslog2 + a, b, _ = _soft_cover_min_only_kernel(A, mi < 0 ? 40 : mi) + return vcat(a, b) + end + penalty === abslog1 && throw(ArgumentError("MatrixCovers does not implement AbsLog{1} for soft_cover_min")) + throw(ArgumentError(_EXT_IPOPT)) +end + +@noinline _soft_cover_min_only_kernel(A::Matrix{Float64}, mi::Int64) = + MatrixCovers._cover_min_abslog2(A; maxiter = mi, κs = Float64[], boost = false) diff --git a/lib/test/python/test_smoke.py b/lib/test/python/test_smoke.py index abff63e..75ff50f 100644 --- a/lib/test/python/test_smoke.py +++ b/lib/test/python/test_smoke.py @@ -26,6 +26,12 @@ def check(name, cond): a_min = mc.symcover_min(A_sym) check("symcover_min matches reference", np.allclose(a_min, [2.0, 2.0], rtol=1e-6)) +# penalty accepted as a member and as a string +a_min_member = mc.symcover_min(A_sym, penalty=mc.Penalty.abslog2) +a_min_str = mc.symcover_min(A_sym, penalty="abslog2") +check("symcover_min(penalty=Penalty.abslog2) matches reference", np.allclose(a_min_member, [2.0, 2.0], rtol=1e-6)) +check("symcover_min(penalty='abslog2') matches reference", np.allclose(a_min_str, [2.0, 2.0], rtol=1e-6)) + # cover / cover_min a_asym, b_asym = mc.cover(A_asym) check( @@ -71,29 +77,32 @@ def check(name, cond): 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: +except ValueError as e: check("unknown penalty string raises ValueError", True) + check("unknown penalty string names valid members", + all(name in str(e) for name in ("abslog1", "abslog2", "abslinear1", "abslinear2"))) try: from matrixcovers import _lowlevel - _lowlevel.mc_symcover_min( + _lowlevel.matrixcovers_symcover_min( + _lowlevel.CMatrix_borrowed_Float64.from_numpy(np.asfortranarray(A_sym)), 99, - _lowlevel.CMatrix_Float64.from_numpy(np.asfortranarray(A_sym)), - -1, 1, - _lowlevel.CVector_Float64.from_numpy(np.zeros(2)), + _lowlevel.COpt_Int64.from_optional(None), + 1, ) 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) + check("bad low-level penalty enum raises JLWError", e.code == 2) 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) + check("shape-mismatched call raises cleanly", e.code == 3) try: mc.symcover_min(A_sym, penalty="abslinear2") @@ -101,6 +110,15 @@ def check(name, cond): except mc.JLWError as e: check("extension-only penalty raises JLWError code 2", e.code == 2) +try: + mc.symcover_min(A_sym, penalty="abslinear1") + check("symcover_min(abslinear1) raises JLWError code 2 mentioning the extension", False) +except mc.JLWError as e: + check( + "symcover_min(abslinear1) raises JLWError code 2 mentioning the extension", + e.code == 2 and "extension" in e.message, + ) + if failures: print(f"\n{len(failures)} check(s) failed: {failures}") sys.exit(1)