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 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/Project.toml b/lib/Project.toml new file mode 100644 index 0000000..86b3df0 --- /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.2" +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..992ddc2 --- /dev/null +++ b/lib/build.jl @@ -0,0 +1,48 @@ +# 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) + # 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 + 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..635adee --- /dev/null +++ b/lib/python/_facade.py @@ -0,0 +1,300 @@ +"""Python bindings for MatrixCovers.jl. + +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`. +""" +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 _as_matrix(A): + return np.asfortranarray(A, dtype=np.float64) + + +def _as_vector(v): + return np.ascontiguousarray(v, dtype=np.float64) + + +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: + _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 + + +# 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 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 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 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 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])`. + + `b=None` (the default) tests the symmetric cover `a*a'`, and requires `A` + to be square. + """ + if b is None: + 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: + 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): + """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") + if w is not None: + 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", "Penalty", "Linsolve", + "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..e24b793 --- /dev/null +++ b/lib/src/matrixcovers.jl @@ -0,0 +1,247 @@ +# 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: MatrixCovers, AbsLog, AbsLinear +using Random: MersenneTwister + +@export_release_entrypoints + +""" + Penalty + +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 + +""" + Linsolve + +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 + +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)" + +_linsolve_symbol(ls::Linsolve) = ls === auto ? :auto : ls === dense ? :dense : :lsqr + +include("trimmability.jl") + +# --- Hard covers ------------------------------------------------------------ + +"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) + +@api symcover(A::Matrix{Float64}; maxiter::Union{Int64, Nothing} = nothing)::Vector{Float64} + +"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) + +@api cover_ab(A::Matrix{Float64}; maxiter::Union{Int64, Nothing} = nothing)::Vector{Float64} + +""" +ϕ-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 + +@api symcover_min(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing, linsolve::Linsolve = auto)::Vector{Float64} + +""" +ϕ-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 + +@api cover_min_ab(A::Matrix{Float64}; penalty::Penalty = abslog2, + maxiter::Union{Int64, Nothing} = nothing, linsolve::Linsolve = auto)::Vector{Float64} + +# --- 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 + +@api cover_objective_ab(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}; + penalty::Penalty = abslog2)::Float64 + +# --- Gram covers -------------------------------------------------------------- + +"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 + +@api gramcover(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64})::Vector{Float64} + +"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 + +@api gramcover_weighted(a::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}, w::Vector{Float64})::Vector{Float64} + +"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 + +@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 new file mode 100644 index 0000000..75ff50f --- /dev/null +++ b/lib/test/python/test_smoke.py @@ -0,0 +1,125 @@ +"""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)) + +# 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( + "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 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.matrixcovers_symcover_min( + _lowlevel.CMatrix_borrowed_Float64.from_numpy(np.asfortranarray(A_sym)), + 99, + _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 == 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 == 3) + +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) + +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) +print("\nALL CHECKS PASSED")