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/.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/README.md b/README.md index 77f4fb3..b8d13a6 100644 --- a/README.md +++ b/README.md @@ -74,3 +74,28 @@ true See the [documentation](https://HolyLab.github.io/MatrixCovers.jl/dev/) for motivation, examples, and a full API reference. + +## Python + +A subset of MatrixCovers is available as a compiled Python package that does +not require Julia. Install a wheel from a +[GitHub release](https://github.com/HolyLab/MatrixCovers.jl/releases): + +``` +pip install https://github.com/HolyLab/MatrixCovers.jl/releases/download/vX.Y.Z/matrixcovers-X.Y.Z-py3-none-manylinux_2_35_x86_64.whl +``` + +Replace `X.Y.Z` with a released version. Wheels support Linux x86_64 with +glibc >= 2.35 (Ubuntu 22.04+, +Debian 12+, Fedora 36+; not RHEL/Rocky 9). + +```python +import numpy as np +import matrixcovers as mc + +A = np.array([[4.0, 2.0], [2.0, 16.0]]) +a = mc.symcover(A) # a[i] * a[j] >= abs(A[i, j]) +print(mc.iscover(a, A)) # True +``` + +See `lib/python/_facade.py` for the Python API. diff --git a/lib/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..26b0f62 --- /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.1.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")