From 94b6edf3f7f0e6f34abd23963d6fec6491fb7202 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 18 Sep 2026 13:28:03 +0200 Subject: [PATCH 1/5] Add Pennon --- README.md | 16 +- src/Pennon.jl | 508 +++++++++++++++++++++++++++++++++++++++++++++++ src/Penopt.jl | 20 ++ test/pennon.jl | 86 ++++++++ test/runtests.jl | 1 + 5 files changed, 630 insertions(+), 1 deletion(-) create mode 100644 src/Pennon.jl create mode 100644 test/pennon.jl diff --git a/README.md b/README.md index 7652e1d..0ba0f84 100644 --- a/README.md +++ b/README.md @@ -60,6 +60,19 @@ location of the library, update `PENOPT_LIBPENBMI`, re-run Whether PENBMI is available is given by `Penopt.has_penbmi()`. +### PENNON + +The nonlinear PENNON backend is enabled by setting `PENOPT_LIBPENNON` to a +PENNON shared library or static archive and rebuilding the package. On Linux +x86-64, `deps/LX/libpennon_LX64.a` is detected automatically when present. +Availability is reported by `Penopt.has_pennon()`. + +`Penopt.Pennon.Optimizer` accepts nonlinear scalar objectives and constraints, +and accepts nonlinear semidefinite constraints as an +`MOI.VectorNonlinearFunction` in an +`MOI.PositiveSemidefiniteConeTriangle`. Derivatives are evaluated with MOI's +sparse reverse-mode automatic differentiation. + You can test the installation with `using Pkg; Pkg.test("Penopt")` in a Julia session. @@ -67,7 +80,8 @@ session. Pick the solver explicitly: `Penopt.SDP.Optimizer` solves semidefinite programs with PENSDP, `Penopt.BMI.Optimizer` solves bilinear matrix inequalities and -quadratic objectives with PENBMI. +quadratic objectives with PENBMI, and `Penopt.Pennon.Optimizer` solves general +nonlinear semidefinite programs with PENNON. ```julia using JuMP, Penopt diff --git a/src/Pennon.jl b/src/Pennon.jl new file mode 100644 index 0000000..83518a6 --- /dev/null +++ b/src/Pennon.jl @@ -0,0 +1,508 @@ +# Copyright (c) 2019: Benoît Legat and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +# PENNON uses callbacks with Fortran-style output arguments. The C interface is +# not re-entrant, so keeping the callback state in the opaque user-data pointer +# also documents that one PENNON solve is active at a time. +mutable struct _PennonCallbackState + evaluator::Any + n::Int + constraint_values::Vector{Cdouble} + callback_error::Any +end + +function _pennon_state(data::Ptr{Cvoid}) + return unsafe_pointer_to_objref(data)::_PennonCallbackState +end + +function _pennon_x(state, x) + return unsafe_wrap(Vector{Cdouble}, x, state.n; own = false) +end + +function _pennon_callback(f, state, default) + try + return f() + catch err + state.callback_error = (err, catch_backtrace()) + return default + end +end + +function _pennon_f(x, value, data)::Cvoid + state = _pennon_state(data) + v = _pennon_callback(state, NaN) do + MOI.eval_objective(state.evaluator, _pennon_x(state, x)) + end + unsafe_store!(value, v) + return +end + +function _pennon_df(x, nnz, index, value, data)::Cvoid + state = _pennon_state(data) + _pennon_callback(state, nothing) do + structure = 1:state.n + gradient = unsafe_wrap(Vector{Cdouble}, value, length(structure); own = false) + MOI.eval_objective_gradient(state.evaluator, gradient, _pennon_x(state, x)) + unsafe_store!(nnz, Cint(length(structure))) + for (i, j) in enumerate(structure) + unsafe_store!(index, Cint(j), i) + end + end + return +end + +function _pennon_hessian!(evaluate, state, row, col, value, structure) + values = unsafe_wrap(Vector{Cdouble}, value, length(structure); own = false) + evaluate(values) + for (k, (i, j)) in enumerate(structure) + unsafe_store!(row, Cint(i), k) + unsafe_store!(col, Cint(j), k) + end + return +end + +function _pennon_hf(x, nnz, row, col, value, data)::Cvoid + state = _pennon_state(data) + _pennon_callback(state, nothing) do + structure = MOI.hessian_objective_structure(state.evaluator) + unsafe_store!(nnz, Cint(length(structure))) + _pennon_hessian!(state, row, col, value, structure) do values + MOI.eval_hessian_objective( + state.evaluator, + values, + _pennon_x(state, x), + ) + end + end + return +end + +function _pennon_g(index, x, value, data)::Cvoid + state = _pennon_state(data) + _pennon_callback(state, nothing) do + MOI.eval_constraint( + state.evaluator, + state.constraint_values, + _pennon_x(state, x), + ) + unsafe_store!(value, state.constraint_values[unsafe_load(index)+1]) + end + return +end + +function _pennon_dg(index, x, nnz, variable, value, data)::Cvoid + state = _pennon_state(data) + _pennon_callback(state, nothing) do + i = unsafe_load(index) + 1 + structure = MOI.constraint_gradient_structure(state.evaluator, i) + gradient = unsafe_wrap(Vector{Cdouble}, value, length(structure); own = false) + MOI.eval_constraint_gradient( + state.evaluator, + gradient, + _pennon_x(state, x), + i, + ) + unsafe_store!(nnz, Cint(length(structure))) + for (k, j) in enumerate(structure) + unsafe_store!(variable, Cint(j), k) + end + end + return +end + +function _pennon_hg(index, x, nnz, row, col, value, data)::Cvoid + state = _pennon_state(data) + _pennon_callback(state, nothing) do + i = unsafe_load(index) + 1 + structure = MOI.hessian_constraint_structure(state.evaluator, i) + unsafe_store!(nnz, Cint(length(structure))) + _pennon_hessian!(state, row, col, value, structure) do values + MOI.eval_hessian_constraint( + state.evaluator, + values, + _pennon_x(state, x), + i, + ) + end + end + return +end + +const _PENNON_F = Ref{Ptr{Cvoid}}(C_NULL) +const _PENNON_DF = Ref{Ptr{Cvoid}}(C_NULL) +const _PENNON_HF = Ref{Ptr{Cvoid}}(C_NULL) +const _PENNON_G = Ref{Ptr{Cvoid}}(C_NULL) +const _PENNON_DG = Ref{Ptr{Cvoid}}(C_NULL) +const _PENNON_HG = Ref{Ptr{Cvoid}}(C_NULL) +const _PENNON_LOCK = ReentrantLock() + +function _init_pennon_callbacks!() + _PENNON_F[] = @cfunction( + _pennon_f, + Cvoid, + (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}), + ) + _PENNON_DF[] = @cfunction( + _pennon_df, + Cvoid, + (Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cvoid}), + ) + _PENNON_HF[] = @cfunction( + _pennon_hf, + Cvoid, + (Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cvoid}), + ) + _PENNON_G[] = @cfunction( + _pennon_g, + Cvoid, + (Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}), + ) + _PENNON_DG[] = @cfunction( + _pennon_dg, + Cvoid, + (Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cvoid}), + ) + _PENNON_HG[] = @cfunction( + _pennon_hg, + Cvoid, + ( + Ptr{Cint}, + Ptr{Cdouble}, + Ptr{Cint}, + Ptr{Cint}, + Ptr{Cint}, + Ptr{Cdouble}, + Ptr{Cvoid}, + ), + ) + return +end + +const PENNON_INFINITY = 1.0e38 +const DEFAULT_PENNON_IOPTIONS = Cint[ + 100, 100, 2, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, +] +const DEFAULT_PENNON_DOPTIONS = Cdouble[ + 1e-2, 1.0, 1.0, 1e-2, 0.5, 0.5, 1e-6, 1e-12, 1e-7, 0.05, 1.0, 1.0, + 1.0, 1.0, +] + +mutable struct PennonOptimizer <: MOI.AbstractOptimizer + num_variables::Int + x0::Vector{Cdouble} + objective::MOI.ScalarNonlinearFunction + objective_sign::Cdouble + constraints::Vector{Tuple{MOI.ScalarNonlinearFunction,Any}} + matrix_functions::Vector{MOI.VectorNonlinearFunction} + matrix_sizes::Vector{Cint} + x::Vector{Cdouble} + objective_value::Cdouble + info::Cint + ioptions::Vector{Cint} + doptions::Vector{Cdouble} + silent::Bool + function PennonOptimizer() + return new( + 0, + Cdouble[], + convert(MOI.ScalarNonlinearFunction, 0.0), + 1.0, + Tuple{MOI.ScalarNonlinearFunction,Any}[], + MOI.VectorNonlinearFunction[], + Cint[], + Cdouble[], + NaN, + -1, + copy(DEFAULT_PENNON_IOPTIONS), + copy(DEFAULT_PENNON_DOPTIONS), + false, + ) + end +end + +MOI.get(::PennonOptimizer, ::MOI.SolverName) = "Pennon" +MOI.supports_incremental_interface(::PennonOptimizer) = true +MOI.copy_to(dest::PennonOptimizer, src::MOI.ModelLike) = + MOI.Utilities.default_copy_to(dest, src) + +function MOI.is_empty(model::PennonOptimizer) + return model.num_variables == 0 && isempty(model.constraints) && + isempty(model.matrix_functions) +end + +function MOI.empty!(model::PennonOptimizer) + model.num_variables = 0 + empty!(model.x0) + model.objective = convert(MOI.ScalarNonlinearFunction, 0.0) + model.objective_sign = 1.0 + empty!(model.constraints) + empty!(model.matrix_functions) + empty!(model.matrix_sizes) + empty!(model.x) + model.objective_value = NaN + model.info = -1 + return +end + +function MOI.add_variable(model::PennonOptimizer) + model.num_variables += 1 + push!(model.x0, 0.0) + return MOI.VariableIndex(model.num_variables) +end +MOI.get(model::PennonOptimizer, ::MOI.NumberOfVariables) = model.num_variables + +MOI.supports(::PennonOptimizer, ::MOI.Silent) = true +MOI.get(model::PennonOptimizer, ::MOI.Silent) = model.silent +MOI.set(model::PennonOptimizer, ::MOI.Silent, value::Bool) = model.silent = value + +MOI.supports(::PennonOptimizer, ::MOI.ObjectiveSense) = true +function MOI.set(model::PennonOptimizer, ::MOI.ObjectiveSense, sense) + model.objective_sign = sense == MOI.MAX_SENSE ? -1.0 : 1.0 + if sense == MOI.FEASIBILITY_SENSE + model.objective = convert(MOI.ScalarNonlinearFunction, 0.0) + end + return +end + +for F in ( + MOI.VariableIndex, + MOI.ScalarAffineFunction{Cdouble}, + MOI.ScalarQuadraticFunction{Cdouble}, + MOI.ScalarNonlinearFunction, +) + @eval MOI.supports(::PennonOptimizer, ::MOI.ObjectiveFunction{$F}) = true + @eval function MOI.set( + model::PennonOptimizer, + ::MOI.ObjectiveFunction{$F}, + f::$F, + ) + model.objective = convert(MOI.ScalarNonlinearFunction, f) + return + end +end + +function MOI.supports( + ::PennonOptimizer, + ::MOI.VariablePrimalStart, + ::Type{MOI.VariableIndex}, +) + return true +end +function MOI.set( + model::PennonOptimizer, + ::MOI.VariablePrimalStart, + x::MOI.VariableIndex, + value, +) + model.x0[x.value] = value === nothing ? 0.0 : value + return +end + +const _PENNON_SCALAR_FUNCTIONS = Union{ + MOI.VariableIndex, + MOI.ScalarAffineFunction{Cdouble}, + MOI.ScalarQuadraticFunction{Cdouble}, + MOI.ScalarNonlinearFunction, +} +const _PENNON_SCALAR_SETS = Union{ + MOI.LessThan{Cdouble}, + MOI.GreaterThan{Cdouble}, + MOI.EqualTo{Cdouble}, + MOI.Interval{Cdouble}, +} + +function MOI.supports_constraint( + ::PennonOptimizer, + ::Type{F}, + ::Type{S}, +) where {F<:_PENNON_SCALAR_FUNCTIONS,S<:_PENNON_SCALAR_SETS} + return true +end + +function MOI.add_constraint( + model::PennonOptimizer, + f::F, + set::S, +) where {F<:_PENNON_SCALAR_FUNCTIONS,S<:_PENNON_SCALAR_SETS} + push!(model.constraints, (convert(MOI.ScalarNonlinearFunction, f), set)) + return MOI.ConstraintIndex{F,S}(length(model.constraints)) +end + +function MOI.supports_constraint( + ::PennonOptimizer, + ::Type{MOI.VectorNonlinearFunction}, + ::Type{MOI.PositiveSemidefiniteConeTriangle}, +) + return true +end + +function MOI.add_constraint( + model::PennonOptimizer, + f::MOI.VectorNonlinearFunction, + set::MOI.PositiveSemidefiniteConeTriangle, +) + MOI.output_dimension(f) == MOI.dimension(set) || + throw(MOI.DimensionMismatch("Function dimension does not match cone")) + push!(model.matrix_functions, f) + push!(model.matrix_sizes, set.side_dimension) + return MOI.ConstraintIndex{typeof(f),typeof(set)}(length(model.matrix_sizes)) +end + +function _pennon_bounds(set::MOI.LessThan) + return -PENNON_INFINITY, set.upper +end +function _pennon_bounds(set::MOI.GreaterThan) + return set.lower, PENNON_INFINITY +end +_pennon_bounds(set::MOI.EqualTo) = (set.value, set.value) +_pennon_bounds(set::MOI.Interval) = (set.lower, set.upper) + +function _pennon_evaluator(model::PennonOptimizer) + data = MOI.Nonlinear.Model() + objective = MOI.ScalarNonlinearFunction( + :*, + Any[model.objective_sign, model.objective], + ) + MOI.Nonlinear.set_objective(data, objective) + lower = Cdouble[] + upper = Cdouble[] + for (f, set) in model.constraints + MOI.Nonlinear.add_constraint(data, f, set) + l, u = _pennon_bounds(set) + push!(lower, l) + push!(upper, u) + end + slack = model.num_variables + for f in model.matrix_functions + for row in f.rows + slack += 1 + equality = MOI.ScalarNonlinearFunction( + :-, + Any[row, MOI.VariableIndex(slack)], + ) + MOI.Nonlinear.add_constraint(data, equality, MOI.EqualTo(0.0)) + push!(lower, 0.0) + push!(upper, 0.0) + end + end + variables = MOI.VariableIndex.(1:slack) + evaluator = MOI.Nonlinear.Evaluator( + data, + MOI.Nonlinear.SparseReverseMode(), + variables, + ) + MOI.initialize(evaluator, [:Grad, :Jac, :Hess]) + return evaluator, lower, upper +end + +function _pennon_call!(model::PennonOptimizer, evaluator, lower, upper) + has_pennon() || error( + "PENNON is not available. Set `PENOPT_LIBPENNON` to the path of " * + "`libpennon` and re-run `Pkg.build(\"Penopt\")`.", + ) + n = model.num_variables + sum(_tridim, model.matrix_sizes; init = 0) + nconstr = length(lower) + nlin = 0 + nsdp = length(model.matrix_sizes) + lbv = fill(Cdouble(-PENNON_INFINITY), n) + ubv = fill(Cdouble(PENNON_INFINITY), n) + lbmv = zeros(Cdouble, nsdp) + ubmv = fill(Cdouble(PENNON_INFINITY), nsdp) + # Keep the matrix variables explicit. PENNON's automatic slack-removal + # mode assumes a library-specific linear callback representation, whereas + # these equalities are supplied through the general nonlinear callbacks. + mtype = zeros(Cint, nsdp) + mnzs = Cint[_tridim(s) for s in model.matrix_sizes] + mrow = Cint[] + mcol = Cint[] + for side in model.matrix_sizes, col in 0:(side-1), row in 0:col + push!(mrow, row) + push!(mcol, col) + end + model.x = vcat(model.x0, zeros(Cdouble, n - model.num_variables)) + multipliers = zeros(Cdouble, max(1, 2n + 2nconstr)) + state = _PennonCallbackState(evaluator, n, zeros(nconstr), nothing) + ioptions = copy(model.ioptions) + if model.silent + ioptions[3] = 0 + end + data = pointer_from_objref(state) + lock(_PENNON_LOCK) do + model.info = GC.@preserve state ccall( + (:pennlp, libpennon), + Cint, + ( + Cint, Cint, Cint, Cint, Ptr{Cint}, Cint, Cint, + Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, + Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, + Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, + Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, + Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, + Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cvoid}, Ptr{Cvoid}, + Ptr{Cvoid}, Ptr{Cvoid}, + ), + n, nlin, nconstr, nsdp, model.matrix_sizes, n, n * (n + 1) ÷ 2, + lbv, ubv, lower, upper, lbmv, ubmv, mtype, mnzs, mrow, mcol, + model.x, multipliers, + _PENNON_F[], _PENNON_DF[], _PENNON_HF[], C_NULL, + _PENNON_G[], _PENNON_DG[], _PENNON_HG[], C_NULL, + C_NULL, C_NULL, C_NULL, ioptions, model.doptions, + C_NULL, C_NULL, data, C_NULL, + ) + end + if state.callback_error !== nothing + err, bt = state.callback_error + throw(CapturedException(err, bt)) + end + model.objective_value = MOI.eval_objective( + evaluator, + view(model.x, 1:n), + ) * model.objective_sign + return +end + +function MOI.optimize!(model::PennonOptimizer) + evaluator, lower, upper = _pennon_evaluator(model) + _pennon_call!(model, evaluator, lower, upper) + return +end + +MOI.get(model::PennonOptimizer, ::MOI.RawStatusString) = + model.info == -1 ? "Optimize not called" : "PENNON return code $(model.info)" + +function MOI.get(model::PennonOptimizer, ::MOI.TerminationStatus) + model.info == -1 && return MOI.OPTIMIZE_NOT_CALLED + model.info == 0 && return MOI.LOCALLY_SOLVED + return MOI.OTHER_ERROR +end + +MOI.get(model::PennonOptimizer, ::MOI.ResultCount) = model.info == -1 ? 0 : 1 +function MOI.get(model::PennonOptimizer, attr::MOI.ObjectiveValue) + MOI.check_result_index_bounds(model, attr) + return model.objective_value +end +function MOI.get(model::PennonOptimizer, attr::MOI.VariablePrimal, x::MOI.VariableIndex) + MOI.check_result_index_bounds(model, attr) + return model.x[x.value] +end +function MOI.get(model::PennonOptimizer, attr::MOI.PrimalStatus) + attr.result_index > MOI.get(model, MOI.ResultCount()) && return MOI.NO_SOLUTION + return model.info == 0 ? MOI.FEASIBLE_POINT : MOI.UNKNOWN_RESULT_STATUS +end + +module Pennon + +import ..Penopt + +""" + Penopt.Pennon.Optimizer() + +Optimizer for nonlinear semidefinite programs. Matrix constraints are accepted +as `MOI.VectorNonlinearFunction`-in-`MOI.PositiveSemidefiniteConeTriangle` and +their derivatives are supplied to PENNON by MOI's sparse reverse-mode AD. +""" +const Optimizer = Penopt.PennonOptimizer + +end # module Pennon diff --git a/src/Penopt.jl b/src/Penopt.jl index db32a8e..44311ad 100644 --- a/src/Penopt.jl +++ b/src/Penopt.jl @@ -15,6 +15,12 @@ else ) end +# Keep source updates loadable with a `deps.jl` generated by an older Penopt +# build. Rebuilding will replace this fallback with the detected library path. +if !@isdefined(libpennon) + const libpennon = "" +end + """ has_penbmi() @@ -24,6 +30,14 @@ Return whether the commercial PENBMI library, required by [`penbmi`](@ref) and """ has_penbmi() = !isempty(libpenbmi) +""" + has_pennon() + +Return whether the PENNON library used by [`Penopt.Pennon.Optimizer`](@ref) was +found by `Pkg.build("Penopt")`. +""" +has_pennon() = !isempty(libpennon) + const DEFAULT_IOPTIONS = Cint[1, 50, 100, 2, 0, 0, 0, 0, 0, 0, 1, 0] const DEFAULT_FOPTIONS = Cdouble[ 1.0, @@ -46,5 +60,11 @@ _tridim(n) = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n)) include("MOI_wrapper.jl") include("SDP.jl") include("BMI.jl") +include("Pennon.jl") + +function __init__() + _init_pennon_callbacks!() + return +end end # module diff --git a/test/pennon.jl b/test/pennon.jl new file mode 100644 index 0000000..b3a8de7 --- /dev/null +++ b/test/pennon.jl @@ -0,0 +1,86 @@ +# Copyright (c) 2019: Benoît Legat and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +module TestPennon + +using Test +import MathOptInterface as MOI +import Penopt + +function test_interface() + model = Penopt.Pennon.Optimizer() + @test MOI.get(model, MOI.SolverName()) == "Pennon" + @test MOI.supports_incremental_interface(model) + @test MOI.supports( + model, + MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), + ) + @test MOI.supports_constraint( + model, + MOI.VectorNonlinearFunction, + MOI.PositiveSemidefiniteConeTriangle, + ) + MOI.set(model, MOI.Silent(), true) + MOI.add_variable(model) + MOI.empty!(model) + @test MOI.is_empty(model) + @test MOI.get(model, MOI.Silent()) + return +end + +function test_scalar_constraint_dispatch() + model = Penopt.Pennon.Optimizer() + x = MOI.add_variable(model) + @test MOI.get(model, MOI.NumberOfVariables()) == 1 + variable_ci = MOI.add_constraint(model, x, MOI.LessThan(2.0)) + @test variable_ci == + MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(1) + f = MOI.ScalarNonlinearFunction(:sin, Any[x]) + ci = MOI.add_constraint(model, f, MOI.GreaterThan(0.0)) + @test ci == MOI.ConstraintIndex{ + MOI.ScalarNonlinearFunction, + MOI.GreaterThan{Float64}, + }(2) + return +end + +function test_vector_nonlinear_function() + Penopt.has_pennon() || return + model = Penopt.Pennon.Optimizer() + x = MOI.add_variable(model) + MOI.set(model, MOI.VariablePrimalStart(), x, 0.5) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x) + square = MOI.ScalarNonlinearFunction(:^, Any[x, 2]) + entry = MOI.ScalarNonlinearFunction(:-, Any[1.0, square]) + matrix = MOI.VectorNonlinearFunction([entry]) + MOI.add_constraint( + model, + matrix, + MOI.PositiveSemidefiniteConeTriangle(1), + ) + MOI.set(model, MOI.Silent(), true) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED + @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT + @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.0 atol = 1e-5 + @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1.0 atol = 1e-5 + return +end + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$(name)", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end + return +end + +end # module TestPennon + +TestPennon.runtests() diff --git a/test/runtests.jl b/test/runtests.jl index 12961bc..b1d2b08 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,4 +10,5 @@ include("sdp.jl") if Penopt.has_penbmi() include("bmi.jl") end +include("pennon.jl") include("MOI_wrapper.jl") From fd3ae3a90b10f561efc39a18f1b18acafe9185e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 18 Sep 2026 13:29:56 +0200 Subject: [PATCH 2/5] Simplify --- README.md | 6 ++--- deps/build.jl | 67 +++++++++++++++++++++++++++++++++++++-------------- 2 files changed, 52 insertions(+), 21 deletions(-) diff --git a/README.md b/README.md index 0ba0f84..a26fde0 100644 --- a/README.md +++ b/README.md @@ -63,9 +63,9 @@ Whether PENBMI is available is given by `Penopt.has_penbmi()`. ### PENNON The nonlinear PENNON backend is enabled by setting `PENOPT_LIBPENNON` to a -PENNON shared library or static archive and rebuilding the package. On Linux -x86-64, `deps/LX/libpennon_LX64.a` is detected automatically when present. -Availability is reported by `Penopt.has_pennon()`. +PENNON shared library or static archive, re-running `Pkg.build("Penopt")`, and +restarting Julia, as for PENBMI. Availability is reported by +`Penopt.has_pennon()`. `Penopt.Pennon.Optimizer` accepts nonlinear scalar objectives and constraints, and accepts nonlinear semidefinite constraints as an diff --git a/deps/build.jl b/deps/build.jl index 7a1e222..168a94a 100644 --- a/deps/build.jl +++ b/deps/build.jl @@ -7,12 +7,10 @@ # https://github.com/kocvara/pensdp # # PENSDP is the free SDP-only variant of the PENOPT family of solvers. The -# BMI-capable PENBMI is a commercial product; if the user has a PENBMI license, -# the path to `libpenbmi` can be provided via the `PENOPT_LIBPENBMI` -# environment variable and it will be used in addition to the auto-installed -# PENSDP library. The path may point either to a shared library, which is used -# as is, or to the static archive shipped by PENOPT, in which case a shared -# library is linked from it here. +# BMI-capable PENBMI and nonlinear PENNON libraries may be provided with the +# `PENOPT_LIBPENBMI` and `PENOPT_LIBPENNON` environment variables. A path may +# point either to a shared library, which is used as is, or to the static +# archive shipped by PENOPT, in which case a shared library is linked here. using Libdl @@ -120,27 +118,31 @@ end # shared library with an unusual name is still recognized as such. _is_static_library(path) = open(io -> read(io, 8), path) == b"!\n" -# PENBMI is distributed as a static library. Julia can only `ccall` into a -# shared library so the archive is relinked into one, with all its objects -# force-loaded so that the `penbmi` symbol is exported. -function _build_libpenbmi(archive) +# Julia can only `ccall` into a shared library, so PENOPT static archives are +# relinked here. `entrypoint` avoids force-loading optional objects that depend +# on interfaces, such as AMPL, that are not used by this wrapper. +function _build_penopt_library(archive, name; entrypoint = nothing) mkpath(LIB_DIR) - output = joinpath(LIB_DIR, "libpenbmi.$(Libdl.dlext)") + output = joinpath(LIB_DIR, "lib$(name).$(Libdl.dlext)") julia_libdir = joinpath(Sys.BINDIR, Base.LIBDIR, "julia") os, _ = _platform() if os == :linux # PENBMI calls the LP64 BLAS and LAPACK, so the 32-bit-integer OpenBLAS # is linked in rather than the ILP64 one that Julia itself uses. libopenblas = OpenBLAS32_jll.libopenblas_path + archive_args = entrypoint === nothing ? + `-Wl,--whole-archive $(archive) -Wl,--no-whole-archive` : + `-Wl,-u,$(entrypoint) $(archive)` cmd = `gcc -shared -fPIC -o $(output) - -Wl,--no-undefined - -Wl,--whole-archive $(archive) -Wl,--no-whole-archive + -Wl,--no-undefined $(archive_args) -L$(julia_libdir) -Wl,-rpath,$(julia_libdir) $(libopenblas) -Wl,-rpath,$(dirname(libopenblas)) -lgfortran -lpthread -ldl -lm` elseif os == :apple + archive_args = entrypoint === nothing ? + `-Wl,-force_load,$(archive)` : `-Wl,-u,_$(entrypoint) $(archive)` cmd = `gcc -dynamiclib -o $(output) - -Wl,-force_load,$(archive) + $(archive_args) -framework Accelerate -L$(julia_libdir) -Wl,-rpath,$(julia_libdir) -lgfortran -lpthread -lm` @@ -148,8 +150,8 @@ function _build_libpenbmi(archive) error( "Linking `$(archive)` into a shared library is currently " * "implemented for Linux and macOS only. On Windows, build " * - "`penbmi.dll` manually, point `PENOPT_LIBPENBMI` at it and " * - "re-run `Pkg.build(\"Penopt\")`.", + "`lib$(name).dll` manually, point the corresponding PENOPT " * + "environment variable at it, and re-run `Pkg.build(\"Penopt\")`.", ) end @info "Linking shared library" cmd @@ -157,6 +159,8 @@ function _build_libpenbmi(archive) return output end +_build_libpenbmi(archive) = _build_penopt_library(archive, "penbmi") + function _libpenbmi_path() path = get(ENV, "PENOPT_LIBPENBMI", "") if isempty(path) @@ -170,7 +174,21 @@ function _libpenbmi_path() return _is_static_library(path) ? _build_libpenbmi(path) : path end -function _write_deps(libpensdp_path, libpenbmi_path) +function _libpennon_path() + path = get(ENV, "PENOPT_LIBPENNON", "") + if isempty(path) + return path + elseif !isfile(path) + error( + "The `PENOPT_LIBPENNON` environment variable is set to " * + "`$(path)` but no such file exists.", + ) + end + return _is_static_library(path) ? + _build_penopt_library(path, "pennon"; entrypoint = "pennlp") : path +end + +function _write_deps(libpensdp_path, libpenbmi_path, libpennon_path) open(joinpath(DEPS_DIR, "deps.jl"), "w") do io println(io, "# Auto-generated by deps/build.jl - do not edit.") println(io, "import Libdl") @@ -180,6 +198,7 @@ function _write_deps(libpensdp_path, libpenbmi_path) # string. `has_penbmi()` guards every ccall in `Penopt.jl`, so the # `ccall((:penbmi, libpenbmi), ...)` site is never reached. println(io, "const libpenbmi = ", repr(libpenbmi_path)) + println(io, "const libpennon = ", repr(libpennon_path)) println(io) println(io, """ function check_deps() @@ -203,6 +222,14 @@ function _write_deps(libpensdp_path, libpenbmi_path) error(\"\$(libpenbmi) cannot be opened.\") end end + if !isempty(libpennon) + if !isfile(libpennon) + error("\$(libpennon) does not exist.") + end + if Libdl.dlopen_e(libpennon) == C_NULL + error("\$(libpennon) cannot be opened.") + end + end return end""") end @@ -211,13 +238,17 @@ end function main() libpenbmi_path = _libpenbmi_path() + libpennon_path = _libpennon_path() extract_dir = _download_and_extract() libpensdp_path = _build_shared_library(extract_dir) - _write_deps(libpensdp_path, libpenbmi_path) + _write_deps(libpensdp_path, libpenbmi_path, libpennon_path) @info "Penopt: PENSDP installed at $(libpensdp_path)" if !isempty(libpenbmi_path) @info "Penopt: PENBMI installed at $(libpenbmi_path)" end + if !isempty(libpennon_path) + @info "Penopt: PENNON installed at $(libpennon_path)" + end return end From 44433e03416887e7d381c2fd114e0de38f1a7d53 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 18 Sep 2026 14:07:23 +0200 Subject: [PATCH 3/5] Add tests --- Project.toml | 2 +- test/Project.toml | 2 + test/pennon.jl | 171 +++++++++++++++++++++++++++++++++++++++++----- 3 files changed, 157 insertions(+), 18 deletions(-) diff --git a/Project.toml b/Project.toml index c6640ab..adcf845 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Penopt" uuid = "fb8f6f4e-b8f1-11e9-2b95-bdba3fec731e" authors = ["Michal Kocvara ", "Benoît Legat ", "Michael Stingl "] -version = "0.0.0" +version = "0.1.0" [deps] Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" diff --git a/test/Project.toml b/test/Project.toml index adf82b4..45191e5 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,6 @@ [deps] +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Penopt = "fb8f6f4e-b8f1-11e9-2b95-bdba3fec731e" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/pennon.jl b/test/pennon.jl index b3a8de7..5e293cf 100644 --- a/test/pennon.jl +++ b/test/pennon.jl @@ -6,6 +6,8 @@ module TestPennon using Test +using LinearAlgebra +using JuMP import MathOptInterface as MOI import Penopt @@ -48,25 +50,160 @@ end function test_vector_nonlinear_function() Penopt.has_pennon() || return - model = Penopt.Pennon.Optimizer() - x = MOI.add_variable(model) - MOI.set(model, MOI.VariablePrimalStart(), x, 0.5) - MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) - MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x) - square = MOI.ScalarNonlinearFunction(:^, Any[x, 2]) - entry = MOI.ScalarNonlinearFunction(:-, Any[1.0, square]) - matrix = MOI.VectorNonlinearFunction([entry]) - MOI.add_constraint( + model = Model(Penopt.Pennon.Optimizer) + set_silent(model) + @variable(model, x, start = 0.5) + @objective(model, Max, x) + @expression(model, entry, @force_nonlinear(1 - x^2)) + @constraint(model, [entry] in MOI.PositiveSemidefiniteConeTriangle(1)) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test primal_status(model) == MOI.FEASIBLE_POINT + @test objective_value(model) ≈ 1.0 atol = 1e-5 + @test value(x) ≈ 1.0 atol = 1e-5 + return +end + +function _solve_and_check(model) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test primal_status(model) == MOI.FEASIBLE_POINT + @test result_count(model) == 1 + return +end + +# The next two problems use (39) and (40) in Kočvara and Stingl, +# "PENNON: Software for linear and nonlinear matrix inequalities", §5.4: +# https://arxiv.org/pdf/1504.07212 +# Use a small, symmetric instance with an analytic solution instead of the +# paper's rounded 6-by-6 numerical results. Conjugating by D=diag(1,-1,1) preserves +# the eigenvalues and makes the packed off-diagonal entries distinguishable. +function test_nearest_correlation_matrix() + Penopt.has_pennon() || return + model = Model(Penopt.Pennon.Optimizer) + set_silent(model) + @variable(model, X[i = 1:3, j = 1:3], Symmetric, start = Float64(i == j)) + H = [1.0 1.0 -1.0; 1.0 1.0 1.0; -1.0 1.0 1.0] + @objective( model, - matrix, - MOI.PositiveSemidefiniteConeTriangle(1), + Min, + @force_nonlinear( + sum((X[i, j] - H[i, j])^2 for i in 1:3, j in 1:3), + ), ) - MOI.set(model, MOI.Silent(), true) - MOI.optimize!(model) - @test MOI.get(model, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED - @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT - @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.0 atol = 1e-5 - @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1.0 atol = 1e-5 + @constraint(model, [i in 1:3], X[i, i] == 1) + # Explicit nonlinear entries exercise the VectorNonlinearFunction backend + # even when a particular PSD constraint happens to be affine. + @constraint(model, Symmetric(convert.(NonlinearExpr, X)) in PSDCone()) + _solve_and_check(model) + solution = value.(X) + # In the sign-transformed coordinates all off-diagonals equal -1/2: + # the PSD bound 1 + 2r >= 0 is active. The squared Frobenius distance is 3/2. + expected = [1.0 0.5 -0.5; 0.5 1.0 0.5; -0.5 0.5 1.0] + @test solution ≈ expected atol = 1e-5 + @test diag(solution) ≈ ones(3) atol = 1e-6 + @test eigmin(solution) >= -1e-6 + @test objective_value(model) ≈ 1.5 atol = 1e-5 + @test objective_value(model) ≈ sum(abs2, solution - H) atol = 1e-8 + return +end + +function test_conditioned_correlation_matrix() + Penopt.has_pennon() || return + model = Model(Penopt.Pennon.Optimizer) + set_silent(model) + @variable(model, Y[i = 1:3, j = 1:3], Symmetric, start = 2.0 * (i == j)) + @variable(model, 1 <= zeta <= 4, start = 2.0) + H = [1.0 1.0 -1.0; 1.0 1.0 1.0; -1.0 1.0 1.0] + @objective( + model, + Min, + sum((Y[i, j] / zeta - H[i, j])^2 for i in 1:3, j in 1:3), + ) + @constraint(model, [i in 1:3], Y[i, i] == zeta) + # Both spectral bounds are separate PSD blocks, as in (40). + @constraint(model, Symmetric(convert.(NonlinearExpr, Y - I)) in PSDCone()) + @constraint(model, Symmetric(convert.(NonlinearExpr, 4I - Y)) in PSDCone()) + _solve_and_check(model) + solution = value.(Y) + X = solution / value(zeta) + # Eigenvalues of the sign-transformed correlation matrix are 1+2r and + # 1-r. Their ratio is 4 at r=-1/3, giving zeta=3 and squared distance 8/3. + expected = [1.0 1/3 -1/3; 1/3 1.0 1/3; -1/3 1/3 1.0] + @test X ≈ expected atol = 1e-5 + @test value(zeta) ≈ 3.0 atol = 1e-5 + @test diag(X) ≈ ones(3) atol = 1e-6 + @test eigmin(solution) >= 1.0 - 1e-5 + @test eigmax(solution) <= 4.0 + 1e-5 + @test cond(X) ≈ 4.0 atol = 1e-4 + @test objective_value(model) ≈ 8 / 3 atol = 1e-5 + @test objective_value(model) ≈ sum(abs2, X - H) atol = 1e-8 + return +end + +# A one-interval instance of the nonnegative spline formulation (45)-(50), +# §5.5 of https://arxiv.org/pdf/1504.07212. The data are exact samples of +# (t-1/2)^2 + 1/4, so the independent optimum is zero residual. +function test_nonnegative_spline() + Penopt.has_pennon() || return + model = Model(Penopt.Pennon.Optimizer) + set_silent(model) + @variable(model, X[i = 1:2, j = 1:2], Symmetric, start = Float64(i == j)) + @variable(model, S[i = 1:2, j = 1:2], Symmetric, start = Float64(i == j)) + @constraint(model, Symmetric(convert.(NonlinearExpr, X)) in PSDCone()) + @constraint(model, Symmetric(convert.(NonlinearExpr, S)) in PSDCone()) + @expression( + model, + coefficients, + [ + S[1, 1], + X[1, 1] - S[1, 1] + 2S[1, 2], + 2X[1, 2] - 2S[1, 2] + S[2, 2], + X[2, 2] - S[2, 2], + ], + ) + samples = collect(0.0:0.25:1.0) + @expression( + model, + prediction[t in samples], + sum(coefficients[k] * t^(k - 1) for k in 1:4), + ) + @objective( + model, + Min, + @force_nonlinear( + sum((prediction[t] - ((t - 0.5)^2 + 0.25))^2 for t in samples), + ), + ) + _solve_and_check(model) + c = value.(coefficients) + @test c ≈ [0.5, -1.0, 1.0, 0.0] atol = 1e-4 + @test eigmin(value.(X)) >= -1e-6 + @test eigmin(value.(S)) >= -1e-6 + residual = sum((evalpoly(t, c) - ((t - 0.5)^2 + 0.25))^2 for t in samples) + @test residual <= 1e-8 + @test objective_value(model) ≈ residual atol = 1e-10 + @test minimum(evalpoly(t, c) for t in 0.0:0.01:1.0) >= 0.0 + return +end + +function test_nonlinear_matrix_off_diagonal() + Penopt.has_pennon() || return + model = Model(Penopt.Pennon.Optimizer) + set_silent(model) + @variable(model, x, start = 0.2) + @variable(model, y, start = 0.4) + @objective(model, Max, x + 2y) + @expression(model, matrix, [1-x^2 x*y; x*y 1-y^2]) + @constraint(model, Symmetric(convert.(NonlinearExpr, matrix)) in PSDCone()) + # This matrix is PSD exactly when x^2 + y^2 <= 1. Cauchy-Schwarz + # gives the unique maximizer (1,2)/sqrt(5) and objective sqrt(5). + _solve_and_check(model) + xv, yv = value(x), value(y) + @test [xv, yv] ≈ [1.0, 2.0] / sqrt(5.0) atol = 1e-5 + @test xv^2 + yv^2 <= 1.0 + 1e-6 + @test eigmin(Symmetric([1-xv^2 xv*yv; xv*yv 1-yv^2])) >= -1e-6 + @test objective_value(model) ≈ sqrt(5.0) atol = 1e-5 return end From cf77a697542e35ca93bf98fa08c66730ddbd539c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 18 Sep 2026 14:58:19 +0200 Subject: [PATCH 4/5] Use a bridge --- src/Pennon.jl | 105 +++++++++++++++++++++++++++++++++++++------------ test/pennon.jl | 90 ++++++++++++++++++++++++++++-------------- 2 files changed, 139 insertions(+), 56 deletions(-) diff --git a/src/Pennon.jl b/src/Pennon.jl index 83518a6..e139ff3 100644 --- a/src/Pennon.jl +++ b/src/Pennon.jl @@ -223,6 +223,61 @@ mutable struct PennonOptimizer <: MOI.AbstractOptimizer end MOI.get(::PennonOptimizer, ::MOI.SolverName) = "Pennon" + +# MOI's FunctionConversionBridge does not yet implement conversion to +# VectorNonlinearFunction. Keep this conversion in a bridge, not the optimizer. +struct _PennonNonlinearPSDBridge <: + MOI.Bridges.Constraint.AbstractFunctionConversionBridge{ + MOI.VectorNonlinearFunction, + MOI.PositiveSemidefiniteConeTriangle, +} + constraint::MOI.ConstraintIndex{ + MOI.VectorNonlinearFunction, + MOI.PositiveSemidefiniteConeTriangle, + } +end + +function MOI.get( + ::PennonOptimizer, + ::MOI.Bridges.ListOfNonstandardBridges{Cdouble}, +) + return Type[_PennonNonlinearPSDBridge] +end + +function MOI.supports_constraint( + ::Type{_PennonNonlinearPSDBridge}, + ::Type{<:Union{ + MOI.VectorOfVariables, + MOI.VectorAffineFunction{Cdouble}, + MOI.VectorQuadraticFunction{Cdouble}, + }}, + ::Type{MOI.PositiveSemidefiniteConeTriangle}, +) + return true +end + +function MOI.Bridges.Constraint.concrete_bridge_type( + ::Type{_PennonNonlinearPSDBridge}, + ::Type{<:MOI.AbstractVectorFunction}, + ::Type{MOI.PositiveSemidefiniteConeTriangle}, +) + return _PennonNonlinearPSDBridge +end + +function MOI.Bridges.Constraint.bridge_constraint( + ::Type{_PennonNonlinearPSDBridge}, + model::MOI.ModelLike, + f::MOI.AbstractVectorFunction, + set::MOI.PositiveSemidefiniteConeTriangle, +) + g = MOI.VectorNonlinearFunction( + MOI.ScalarNonlinearFunction[row for row in MOI.Utilities.eachscalar(f)], + ) + return _PennonNonlinearPSDBridge(MOI.add_constraint(model, g, set)) +end + +MOI.Bridges.bridging_cost(::Type{_PennonNonlinearPSDBridge}) = 100.0 + MOI.supports_incremental_interface(::PennonOptimizer) = true MOI.copy_to(dest::PennonOptimizer, src::MOI.ModelLike) = MOI.Utilities.default_copy_to(dest, src) @@ -266,21 +321,23 @@ function MOI.set(model::PennonOptimizer, ::MOI.ObjectiveSense, sense) return end -for F in ( - MOI.VariableIndex, - MOI.ScalarAffineFunction{Cdouble}, - MOI.ScalarQuadraticFunction{Cdouble}, - MOI.ScalarNonlinearFunction, +# Only accept nonlinear functions and let MOI bridges convert affine and +# quadratic inputs. TODO: add native quadratic support once +# https://github.com/jump-dev/MathOptInterface.jl/pull/3048 lands, so that we can +# use specialized quadratic AD instead of the slightly slower nonlinear AD. +function MOI.supports( + ::PennonOptimizer, + ::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}, ) - @eval MOI.supports(::PennonOptimizer, ::MOI.ObjectiveFunction{$F}) = true - @eval function MOI.set( - model::PennonOptimizer, - ::MOI.ObjectiveFunction{$F}, - f::$F, - ) - model.objective = convert(MOI.ScalarNonlinearFunction, f) - return - end + return true +end +function MOI.set( + model::PennonOptimizer, + ::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}, + f::MOI.ScalarNonlinearFunction, +) + model.objective = f + return end function MOI.supports( @@ -300,12 +357,6 @@ function MOI.set( return end -const _PENNON_SCALAR_FUNCTIONS = Union{ - MOI.VariableIndex, - MOI.ScalarAffineFunction{Cdouble}, - MOI.ScalarQuadraticFunction{Cdouble}, - MOI.ScalarNonlinearFunction, -} const _PENNON_SCALAR_SETS = Union{ MOI.LessThan{Cdouble}, MOI.GreaterThan{Cdouble}, @@ -315,19 +366,21 @@ const _PENNON_SCALAR_SETS = Union{ function MOI.supports_constraint( ::PennonOptimizer, - ::Type{F}, + ::Type{MOI.ScalarNonlinearFunction}, ::Type{S}, -) where {F<:_PENNON_SCALAR_FUNCTIONS,S<:_PENNON_SCALAR_SETS} +) where {S<:_PENNON_SCALAR_SETS} return true end function MOI.add_constraint( model::PennonOptimizer, - f::F, + f::MOI.ScalarNonlinearFunction, set::S, -) where {F<:_PENNON_SCALAR_FUNCTIONS,S<:_PENNON_SCALAR_SETS} - push!(model.constraints, (convert(MOI.ScalarNonlinearFunction, f), set)) - return MOI.ConstraintIndex{F,S}(length(model.constraints)) +) where {S<:_PENNON_SCALAR_SETS} + push!(model.constraints, (f, set)) + return MOI.ConstraintIndex{MOI.ScalarNonlinearFunction,S}( + length(model.constraints), + ) end function MOI.supports_constraint( diff --git a/test/pennon.jl b/test/pennon.jl index 5e293cf..62c1313 100644 --- a/test/pennon.jl +++ b/test/pennon.jl @@ -11,10 +11,18 @@ using JuMP import MathOptInterface as MOI import Penopt +function _model() + model = Model(Penopt.Pennon.Optimizer) + set_silent(model) + return model +end + function test_interface() model = Penopt.Pennon.Optimizer() @test MOI.get(model, MOI.SolverName()) == "Pennon" @test MOI.supports_incremental_interface(model) + @test MOI.get(model, MOI.Bridges.ListOfNonstandardBridges{Float64}()) == + [Penopt._PennonNonlinearPSDBridge] @test MOI.supports( model, MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), @@ -24,6 +32,30 @@ function test_interface() MOI.VectorNonlinearFunction, MOI.PositiveSemidefiniteConeTriangle, ) + for F in ( + MOI.VariableIndex, + MOI.ScalarAffineFunction{Float64}, + MOI.ScalarQuadraticFunction{Float64}, + ) + @test !MOI.supports(model, MOI.ObjectiveFunction{F}()) + @test !MOI.supports_constraint(model, F, MOI.LessThan{Float64}) + end + for F in ( + MOI.VectorOfVariables, + MOI.VectorAffineFunction{Float64}, + MOI.VectorQuadraticFunction{Float64}, + ) + @test !MOI.supports_constraint( + model, + F, + MOI.PositiveSemidefiniteConeTriangle, + ) + @test MOI.supports_constraint( + MOI.Bridges.full_bridge_optimizer(model, Float64), + F, + MOI.PositiveSemidefiniteConeTriangle, + ) + end MOI.set(model, MOI.Silent(), true) MOI.add_variable(model) MOI.empty!(model) @@ -36,25 +68,21 @@ function test_scalar_constraint_dispatch() model = Penopt.Pennon.Optimizer() x = MOI.add_variable(model) @test MOI.get(model, MOI.NumberOfVariables()) == 1 - variable_ci = MOI.add_constraint(model, x, MOI.LessThan(2.0)) - @test variable_ci == - MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(1) f = MOI.ScalarNonlinearFunction(:sin, Any[x]) ci = MOI.add_constraint(model, f, MOI.GreaterThan(0.0)) @test ci == MOI.ConstraintIndex{ MOI.ScalarNonlinearFunction, MOI.GreaterThan{Float64}, - }(2) + }(1) return end -function test_vector_nonlinear_function() +function test_quadratic_matrix_constraint() Penopt.has_pennon() || return - model = Model(Penopt.Pennon.Optimizer) - set_silent(model) + model = _model() @variable(model, x, start = 0.5) @objective(model, Max, x) - @expression(model, entry, @force_nonlinear(1 - x^2)) + @expression(model, entry, 1 - x^2) @constraint(model, [entry] in MOI.PositiveSemidefiniteConeTriangle(1)) optimize!(model) @test termination_status(model) == MOI.LOCALLY_SOLVED @@ -64,6 +92,18 @@ function test_vector_nonlinear_function() return end +function test_quadratic_scalar_constraint() + Penopt.has_pennon() || return + model = _model() + @variable(model, x >= 0, start = 1.5) + @objective(model, Min, x^2) + @constraint(model, (x - 2)^2 <= 1) + _solve_and_check(model) + @test value(x) ≈ 1.0 atol = 1e-5 + @test objective_value(model) ≈ 1.0 atol = 1e-5 + return +end + function _solve_and_check(model) optimize!(model) @test termination_status(model) == MOI.LOCALLY_SOLVED @@ -80,21 +120,16 @@ end # the eigenvalues and makes the packed off-diagonal entries distinguishable. function test_nearest_correlation_matrix() Penopt.has_pennon() || return - model = Model(Penopt.Pennon.Optimizer) - set_silent(model) + model = _model() @variable(model, X[i = 1:3, j = 1:3], Symmetric, start = Float64(i == j)) H = [1.0 1.0 -1.0; 1.0 1.0 1.0; -1.0 1.0 1.0] @objective( model, Min, - @force_nonlinear( - sum((X[i, j] - H[i, j])^2 for i in 1:3, j in 1:3), - ), + sum((X[i, j] - H[i, j])^2 for i in 1:3, j in 1:3), ) @constraint(model, [i in 1:3], X[i, i] == 1) - # Explicit nonlinear entries exercise the VectorNonlinearFunction backend - # even when a particular PSD constraint happens to be affine. - @constraint(model, Symmetric(convert.(NonlinearExpr, X)) in PSDCone()) + @constraint(model, X in PSDCone()) _solve_and_check(model) solution = value.(X) # In the sign-transformed coordinates all off-diagonals equal -1/2: @@ -110,8 +145,7 @@ end function test_conditioned_correlation_matrix() Penopt.has_pennon() || return - model = Model(Penopt.Pennon.Optimizer) - set_silent(model) + model = _model() @variable(model, Y[i = 1:3, j = 1:3], Symmetric, start = 2.0 * (i == j)) @variable(model, 1 <= zeta <= 4, start = 2.0) H = [1.0 1.0 -1.0; 1.0 1.0 1.0; -1.0 1.0 1.0] @@ -122,8 +156,8 @@ function test_conditioned_correlation_matrix() ) @constraint(model, [i in 1:3], Y[i, i] == zeta) # Both spectral bounds are separate PSD blocks, as in (40). - @constraint(model, Symmetric(convert.(NonlinearExpr, Y - I)) in PSDCone()) - @constraint(model, Symmetric(convert.(NonlinearExpr, 4I - Y)) in PSDCone()) + @constraint(model, Y - I in PSDCone()) + @constraint(model, 4I - Y in PSDCone()) _solve_and_check(model) solution = value.(Y) X = solution / value(zeta) @@ -146,12 +180,11 @@ end # (t-1/2)^2 + 1/4, so the independent optimum is zero residual. function test_nonnegative_spline() Penopt.has_pennon() || return - model = Model(Penopt.Pennon.Optimizer) - set_silent(model) + model = _model() @variable(model, X[i = 1:2, j = 1:2], Symmetric, start = Float64(i == j)) @variable(model, S[i = 1:2, j = 1:2], Symmetric, start = Float64(i == j)) - @constraint(model, Symmetric(convert.(NonlinearExpr, X)) in PSDCone()) - @constraint(model, Symmetric(convert.(NonlinearExpr, S)) in PSDCone()) + @constraint(model, X in PSDCone()) + @constraint(model, S in PSDCone()) @expression( model, coefficients, @@ -171,9 +204,7 @@ function test_nonnegative_spline() @objective( model, Min, - @force_nonlinear( - sum((prediction[t] - ((t - 0.5)^2 + 0.25))^2 for t in samples), - ), + sum((prediction[t] - ((t - 0.5)^2 + 0.25))^2 for t in samples), ) _solve_and_check(model) c = value.(coefficients) @@ -189,13 +220,12 @@ end function test_nonlinear_matrix_off_diagonal() Penopt.has_pennon() || return - model = Model(Penopt.Pennon.Optimizer) - set_silent(model) + model = _model() @variable(model, x, start = 0.2) @variable(model, y, start = 0.4) @objective(model, Max, x + 2y) @expression(model, matrix, [1-x^2 x*y; x*y 1-y^2]) - @constraint(model, Symmetric(convert.(NonlinearExpr, matrix)) in PSDCone()) + @constraint(model, Symmetric(matrix) in PSDCone()) # This matrix is PSD exactly when x^2 + y^2 <= 1. Cauchy-Schwarz # gives the unique maximizer (1,2)/sqrt(5) and objective sqrt(5). _solve_and_check(model) From 0645b0411927c5abe0bef52c7e3bd2cdbc4c5866 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 22 Sep 2026 08:06:00 +0200 Subject: [PATCH 5/5] Just use function conversion bridge --- Project.toml | 2 +- src/Pennon.jl | 54 ++++++-------------------------------------------- test/pennon.jl | 5 ++++- 3 files changed, 11 insertions(+), 50 deletions(-) diff --git a/Project.toml b/Project.toml index adcf845..98eb675 100644 --- a/Project.toml +++ b/Project.toml @@ -12,7 +12,7 @@ OpenBLAS32_jll = "656ef2d0-ae68-5445-9ca0-591084a874a2" [compat] Libdl = "1.10" LinearAlgebra = "1.10" -MathOptInterface = "1" +MathOptInterface = "1.54" OpenBLAS32_jll = "0.3" julia = "1.10" diff --git a/src/Pennon.jl b/src/Pennon.jl index e139ff3..8af8e70 100644 --- a/src/Pennon.jl +++ b/src/Pennon.jl @@ -224,60 +224,18 @@ end MOI.get(::PennonOptimizer, ::MOI.SolverName) = "Pennon" -# MOI's FunctionConversionBridge does not yet implement conversion to -# VectorNonlinearFunction. Keep this conversion in a bridge, not the optimizer. -struct _PennonNonlinearPSDBridge <: - MOI.Bridges.Constraint.AbstractFunctionConversionBridge{ - MOI.VectorNonlinearFunction, - MOI.PositiveSemidefiniteConeTriangle, -} - constraint::MOI.ConstraintIndex{ - MOI.VectorNonlinearFunction, - MOI.PositiveSemidefiniteConeTriangle, - } -end - function MOI.get( ::PennonOptimizer, ::MOI.Bridges.ListOfNonstandardBridges{Cdouble}, ) - return Type[_PennonNonlinearPSDBridge] -end - -function MOI.supports_constraint( - ::Type{_PennonNonlinearPSDBridge}, - ::Type{<:Union{ - MOI.VectorOfVariables, - MOI.VectorAffineFunction{Cdouble}, - MOI.VectorQuadraticFunction{Cdouble}, - }}, - ::Type{MOI.PositiveSemidefiniteConeTriangle}, -) - return true + return Type[ + MOI.Bridges.Constraint.FunctionConversionBridge{ + Cdouble, + MOI.VectorNonlinearFunction, + }, + ] end -function MOI.Bridges.Constraint.concrete_bridge_type( - ::Type{_PennonNonlinearPSDBridge}, - ::Type{<:MOI.AbstractVectorFunction}, - ::Type{MOI.PositiveSemidefiniteConeTriangle}, -) - return _PennonNonlinearPSDBridge -end - -function MOI.Bridges.Constraint.bridge_constraint( - ::Type{_PennonNonlinearPSDBridge}, - model::MOI.ModelLike, - f::MOI.AbstractVectorFunction, - set::MOI.PositiveSemidefiniteConeTriangle, -) - g = MOI.VectorNonlinearFunction( - MOI.ScalarNonlinearFunction[row for row in MOI.Utilities.eachscalar(f)], - ) - return _PennonNonlinearPSDBridge(MOI.add_constraint(model, g, set)) -end - -MOI.Bridges.bridging_cost(::Type{_PennonNonlinearPSDBridge}) = 100.0 - MOI.supports_incremental_interface(::PennonOptimizer) = true MOI.copy_to(dest::PennonOptimizer, src::MOI.ModelLike) = MOI.Utilities.default_copy_to(dest, src) diff --git a/test/pennon.jl b/test/pennon.jl index 62c1313..64065ab 100644 --- a/test/pennon.jl +++ b/test/pennon.jl @@ -22,7 +22,10 @@ function test_interface() @test MOI.get(model, MOI.SolverName()) == "Pennon" @test MOI.supports_incremental_interface(model) @test MOI.get(model, MOI.Bridges.ListOfNonstandardBridges{Float64}()) == - [Penopt._PennonNonlinearPSDBridge] + [MOI.Bridges.Constraint.FunctionConversionBridge{ + Float64, + MOI.VectorNonlinearFunction, + }] @test MOI.supports( model, MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(),