diff --git a/docs/src/manual/standard_form.md b/docs/src/manual/standard_form.md index 57b8a8e2d2..1700236547 100644 --- a/docs/src/manual/standard_form.md +++ b/docs/src/manual/standard_form.md @@ -82,6 +82,7 @@ The vector-valued set types implemented in MathOptInterface.jl are: | [`RelativeEntropyCone(d)`](@ref MathOptInterface.RelativeEntropyCone) | ``\{ (u, v, w) \in \mathbb{R}^{d} : u \ge \sum_i w_i \log (\frac{w_i}{v_i}), v_i \ge 0, w_i \ge 0 \}`` | | [`HyperRectangle(l, u)`](@ref MathOptInterface.HyperRectangle) | ``\{x \in \bar{\mathbb{R}}^d: x_i \in [l_i, u_i] \forall i=1,\ldots,d\}`` | | [`NormCone(p, d)`](@ref MathOptInterface.NormCone) | ``\{ (t,x) \in \mathbb{R}^{d} : t \ge \left(\sum\limits_i \lvert x_i \rvert^p\right)^{\frac{1}{p}} \}`` | +| [`VectorNonlinearOracle`](@ref MathOptInterface.VectorNonlinearOracle)| ``\{x \in \mathbb{R}^{dimension}: l \le f(x) \le u \}`` | ## Matrix cones diff --git a/docs/src/reference/standard_form.md b/docs/src/reference/standard_form.md index ffe855b531..d78681869e 100644 --- a/docs/src/reference/standard_form.md +++ b/docs/src/reference/standard_form.md @@ -104,6 +104,7 @@ ACTIVATE_ON_ONE Complements HyperRectangle Scaled +VectorNonlinearOracle ``` ## Constraint programming sets diff --git a/src/Test/test_basic_constraint.jl b/src/Test/test_basic_constraint.jl index bc73cf2632..6e36bf7867 100644 --- a/src/Test/test_basic_constraint.jl +++ b/src/Test/test_basic_constraint.jl @@ -175,6 +175,31 @@ function _set(::Type{T}, ::Type{MOI.HyperRectangle}) where {T} return MOI.HyperRectangle(zeros(T, 3), ones(T, 3)) end +function _set(::Type{T}, ::Type{MOI.VectorNonlinearOracle}) where {T} + set = MOI.VectorNonlinearOracle(; + dimension = 3, + l = T[0, 0], + u = T[1, 0], + eval_f = (ret, x) -> begin + ret[1] = x[2]^2 + ret[2] = x[3]^2 + x[4]^3 - x[1] + return + end, + jacobian_structure = [(1, 2), (2, 1), (2, 3), (2, 4)], + eval_jacobian = (ret, x) -> begin + ret[1] = T(2) * x[2] + ret[2] = -T(1) + ret[3] = T(2) * x[3] + ret[4] = T(3) * x[4]^2 + return + end, + ) + x, ret_f, ret_J = T[1, 2, 3, 4, 5], T[0, 0], T[0, 0, 0, 0] + set.eval_f(ret_f, x) + set.eval_jacobian(ret_J, x) + return set +end + function _test_function_modification( model::MOI.ModelLike, config::Config{T}, @@ -392,6 +417,7 @@ for s in [ :Table, :Path, :HyperRectangle, + :VectorNonlinearOracle, ] S = getfield(MOI, s) functions = if S <: MOI.AbstractScalarSet diff --git a/src/Test/test_nonlinear.jl b/src/Test/test_nonlinear.jl index 79829fbaf7..609dccb5bc 100644 --- a/src/Test/test_nonlinear.jl +++ b/src/Test/test_nonlinear.jl @@ -2223,3 +2223,162 @@ function setup_test( end version_added(::typeof(test_nonlinear_quadratic_4)) = v"1.35.0" + +function test_vector_nonlinear_oracle( + model::MOI.ModelLike, + config::Config{T}, +) where {T} + @requires _supports(config, MOI.optimize!) + @requires MOI.supports_constraint( + model, + MOI.VectorOfVariables, + MOI.VectorNonlinearOracle{T}, + ) + set = MOI.VectorNonlinearOracle(; + dimension = 5, + l = T[0, 0], + u = T[0, 0], + eval_f = (ret, x) -> begin + @test length(ret) == 2 + @test length(x) == 5 + ret[1] = x[1]^2 - x[4] + ret[2] = x[2]^2 + x[3]^3 - x[5] + return + end, + jacobian_structure = [(1, 1), (2, 2), (2, 3), (1, 4), (2, 5)], + eval_jacobian = (ret, x) -> begin + @test length(ret) == 5 + @test length(x) == 5 + ret[1] = T(2) * x[1] + ret[2] = T(2) * x[2] + ret[3] = T(3) * x[3]^2 + ret[4] = -T(1) + ret[5] = -T(1) + return + end, + hessian_lagrangian_structure = [(1, 1), (2, 2), (3, 3)], + eval_hessian_lagrangian = (ret, x, u) -> begin + @test length(ret) == 3 + @test length(x) == 5 + @test length(u) == 2 + ret[1] = T(2) * u[1] + ret[2] = T(2) * u[2] + ret[3] = T(6) * x[3] * u[2] + return + end, + ) + @test MOI.dimension(set) == 5 + x = T[1, 2, 3, 4, 5] + ret = T[0, 0] + set.eval_f(ret, x) + @test ret == T[-3, 26] + ret = T[0, 0, 0, 0, 0] + set.eval_jacobian(ret, x) + @test ret == T[2, 4, 27, -1, -1] + ret = T[0, 0, 0] + set.eval_hessian_lagrangian(ret, x, T[2, 3]) + @test ret == T[4, 6, 54] + x, y = MOI.add_variables(model, 3), MOI.add_variables(model, 2) + MOI.add_constraints.(model, x, MOI.EqualTo.(T(1):T(3))) + c = MOI.add_constraint(model, MOI.VectorOfVariables([x; y]), set) + MOI.optimize!(model) + x_v = MOI.get.(model, MOI.VariablePrimal(), x) + y_v = MOI.get.(model, MOI.VariablePrimal(), y) + @test ≈(y_v, [x_v[1]^2, x_v[2]^2 + x_v[3]^3], config) + @test ≈(MOI.get(model, MOI.ConstraintPrimal(), c), [x_v; y_v], config) + @test ≈(MOI.get(model, MOI.ConstraintDual(), c), zeros(T, 5), config) + return +end + +function setup_test( + ::typeof(test_vector_nonlinear_oracle), + model::MOIU.MockOptimizer, + config::Config{T}, +) where {T} + MOI.Utilities.set_mock_optimize!( + model, + mock -> MOI.Utilities.mock_optimize!( + mock, + config.optimal_status, + T[1, 2, 3, 1, 31], + (MOI.VectorOfVariables, MOI.VectorNonlinearOracle{T}) => + [zeros(T, 5)], + ), + ) + model.eval_variable_constraint_dual = false + return () -> model.eval_variable_constraint_dual = true +end + +version_added(::typeof(test_vector_nonlinear_oracle)) = v"1.46.0" + +function test_vector_nonlinear_oracle_no_hessian( + model::MOI.ModelLike, + config::Config{T}, +) where {T} + @requires _supports(config, MOI.optimize!) + @requires MOI.supports_constraint( + model, + MOI.VectorOfVariables, + MOI.VectorNonlinearOracle{T}, + ) + set = MOI.VectorNonlinearOracle(; + dimension = 5, + l = T[0, 0], + u = T[0, 0], + eval_f = (ret, x) -> begin + ret[1] = x[1]^2 - x[4] + ret[2] = x[2]^2 + x[3]^3 - x[5] + return + end, + jacobian_structure = [(1, 1), (2, 2), (2, 3), (1, 4), (2, 5)], + eval_jacobian = (ret, x) -> begin + ret[1] = T(2) * x[1] + ret[2] = T(2) * x[2] + ret[3] = T(3) * x[3]^2 + ret[4] = -T(1) + ret[5] = -T(1) + return + end, + ) + @test MOI.dimension(set) == 5 + x = T[1, 2, 3, 4, 5] + ret = T[0, 0] + set.eval_f(ret, x) + @test ret == T[-3, 26] + ret = T[0, 0, 0, 0, 0] + set.eval_jacobian(ret, x) + @test ret == T[2, 4, 27, -1, -1] + @test isempty(set.hessian_lagrangian_structure) + @test set.eval_hessian_lagrangian === nothing + x, y = MOI.add_variables(model, 3), MOI.add_variables(model, 2) + MOI.add_constraints.(model, x, MOI.EqualTo.(T(1):T(3))) + c = MOI.add_constraint(model, MOI.VectorOfVariables([x; y]), set) + MOI.optimize!(model) + x_v = MOI.get.(model, MOI.VariablePrimal(), x) + y_v = MOI.get.(model, MOI.VariablePrimal(), y) + @test ≈(y_v, [x_v[1]^2, x_v[2]^2 + x_v[3]^3], config) + @test ≈(MOI.get(model, MOI.ConstraintPrimal(), c), [x_v; y_v], config) + @test ≈(MOI.get(model, MOI.ConstraintDual(), c), zeros(T, 5), config) + return +end + +function setup_test( + ::typeof(test_vector_nonlinear_oracle_no_hessian), + model::MOIU.MockOptimizer, + config::Config{T}, +) where {T} + MOI.Utilities.set_mock_optimize!( + model, + mock -> MOI.Utilities.mock_optimize!( + mock, + config.optimal_status, + T[1, 2, 3, 1, 31], + (MOI.VectorOfVariables, MOI.VectorNonlinearOracle{T}) => + [zeros(T, 5)], + ), + ) + model.eval_variable_constraint_dual = false + return () -> model.eval_variable_constraint_dual = true +end + +version_added(::typeof(test_vector_nonlinear_oracle_no_hessian)) = v"1.46.0" diff --git a/src/Utilities/model.jl b/src/Utilities/model.jl index 184488c3b8..600e97158a 100644 --- a/src/Utilities/model.jl +++ b/src/Utilities/model.jl @@ -843,6 +843,7 @@ const EqualToIndicatorZero{T} = MOI.Table, MOI.BinPacking, MOI.HyperRectangle, + MOI.VectorNonlinearOracle, ), (MOI.ScalarNonlinearFunction,), (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction), diff --git a/src/sets.jl b/src/sets.jl index b564a338b9..857b85774c 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -2666,6 +2666,192 @@ function Base.:(==)(x::Reified{S}, y::Reified{S}) where {S} return x.set == y.set end +""" + VectorNonlinearOracle(; + dimension::Int, + l::Vector{Float64}, + u::Vector{Float64}, + eval_f::Function, + jacobian_structure::Vector{Tuple{Int,Int}}, + eval_jacobian::Function, + hessian_lagrangian_structure::Vector{Tuple{Int,Int}} = Tuple{Int,Int}[], + eval_hessian_lagrangian::Union{Nothing,Function} = nothing, + ) <: AbstractVectorSet + +The set: +```math +S = \\{x \\in \\mathbb{R}^{dimension}: l \\le f(x) \\le u\\} +``` +where ``f`` is defined by the vectors `l` and `u`, and the callback oracles +`eval_f`, `eval_jacobian`, and `eval_hessian_lagrangian`. + +## f + +The `eval_f` function must have the signature +```julia +eval_f(ret::AbstractVector, x::AbstractVector)::Nothing +``` +which fills ``f(x)`` into the dense vector `ret`. + +## Jacobian + +The `eval_jacobian` function must have the signature +```julia +eval_jacobian(ret::AbstractVector, x::AbstractVector)::Nothing +``` +which fills the sparse Jacobian ``\\nabla f(x)`` into `ret`. + +The one-indexed sparsity structure must be provided in the `jacobian_structure` +argument. + +## Hessian + +The `eval_hessian_lagrangian` function is optional. + +If `eval_hessian_lagrangian === nothing`, Ipopt will use a Hessian approximation +instead of the exact Hessian. + +If `eval_hessian_lagrangian` is a function, it must have the signature +```julia +eval_hessian_lagrangian( + ret::AbstractVector, + x::AbstractVector, + μ::AbstractVector, +)::Nothing +``` +which fills the sparse Hessian of the Lagrangian ``\\sum \\mu_i \\nabla^2 f_i(x)`` +into `ret`. + +The one-indexed sparsity structure must be provided in the +`hessian_lagrangian_structure` argument. + +## Example + +To model the set: +```math +\\begin{align} +0 \\le & x^2 \\le 1 +0 \\le & y^2 + z^3 - w \\le 0 +\\end{align} +``` +do +```jldoctest +julia> import MathOptInterface as MOI + +julia> set = MOI.VectorNonlinearOracle(; + dimension = 3, + l = [0.0, 0.0], + u = [1.0, 0.0], + eval_f = (ret, x) -> begin + ret[1] = x[2]^2 + ret[2] = x[3]^2 + x[4]^3 - x[1] + return + end, + jacobian_structure = [(1, 2), (2, 1), (2, 3), (2, 4)], + eval_jacobian = (ret, x) -> begin + ret[1] = 2.0 * x[2] + ret[2] = -1.0 + ret[3] = 2.0 * x[3] + ret[4] = 3.0 * x[4]^2 + return + end, + hessian_lagrangian_structure = [(2, 2), (3, 3), (4, 4)], + eval_hessian_lagrangian = (ret, x, u) -> begin + ret[1] = 2.0 * u[1] + ret[2] = 2.0 * u[2] + ret[3] = 6.0 * x[4] * u[2] + return + end, + ); + +julia> set +VectorNonlinearOracle{Float64}(; + dimension = 3, + l = [0.0, 0.0], + u = [1.0, 0.0], + ..., +) +``` +""" +struct VectorNonlinearOracle{T} <: AbstractVectorSet + input_dimension::Int + output_dimension::Int + l::Vector{T} + u::Vector{T} + eval_f::Function + jacobian_structure::Vector{Tuple{Int,Int}} + eval_jacobian::Function + hessian_lagrangian_structure::Vector{Tuple{Int,Int}} + eval_hessian_lagrangian::Union{Nothing,Function} + + function VectorNonlinearOracle(; + dimension::Int, + l::Vector{T}, + u::Vector{T}, + eval_f::Function, + jacobian_structure::Vector{Tuple{Int,Int}}, + eval_jacobian::Function, + # The hessian_lagrangian is optional. + hessian_lagrangian_structure::Vector{Tuple{Int,Int}} = Tuple{Int,Int}[], + eval_hessian_lagrangian::Union{Nothing,Function} = nothing, + ) where {T} + if length(l) != length(u) + throw(DimensionMismatch()) + end + return new{T}( + dimension, + length(l), + l, + u, + eval_f, + jacobian_structure, + eval_jacobian, + hessian_lagrangian_structure, + eval_hessian_lagrangian, + ) + end +end + +dimension(s::VectorNonlinearOracle) = s.input_dimension + +function Base.copy(s::VectorNonlinearOracle) + return VectorNonlinearOracle(; + dimension = s.input_dimension, + l = copy(s.l), + u = copy(s.u), + eval_f = s.eval_f, + jacobian_structure = copy(s.jacobian_structure), + eval_jacobian = s.eval_jacobian, + hessian_lagrangian_structure = copy(s.hessian_lagrangian_structure), + eval_hessian_lagrangian = s.eval_hessian_lagrangian, + ) +end + +function Base.:(==)( + x::VectorNonlinearOracle{T}, + y::VectorNonlinearOracle{T}, +) where {T} + return x.input_dimension == y.input_dimension && + x.output_dimension == y.output_dimension && + x.l == y.l && + x.u == y.u && + x.eval_f == y.eval_f && + x.jacobian_structure == y.jacobian_structure && + x.eval_jacobian == y.eval_jacobian && + x.hessian_lagrangian_structure == y.hessian_lagrangian_structure && + x.eval_hessian_lagrangian == y.eval_hessian_lagrangian +end + +function Base.show(io::IO, s::VectorNonlinearOracle{T}) where {T} + println(io, "VectorNonlinearOracle{$T}(;") + println(io, " dimension = ", s.input_dimension, ",") + println(io, " l = ", s.l, ",") + println(io, " u = ", s.u, ",") + println(io, " ...,") + print(io, ")") + return +end + # TODO(odow): these are not necessarily isbits. They may not be safe to return # without copying if the number is BigFloat, for example. function Base.copy( diff --git a/test/General/sets.jl b/test/General/sets.jl index fdc32648a4..cffd2d985d 100644 --- a/test/General/sets.jl +++ b/test/General/sets.jl @@ -419,6 +419,54 @@ function test_update_dimension() return end +function test_VectorNonlinearOracle() + @test_throws( + DimensionMismatch, + MOI.VectorNonlinearOracle(; + dimension = 3, + l = [0.0, 0.0, 1.0], + u = [1.0, 0.0], + eval_f = (ret, x) -> nothing, + jacobian_structure = Tuple{Int,Int}[], + eval_jacobian = (ret, x) -> nothing, + ), + ) + set = MOI.VectorNonlinearOracle(; + dimension = 3, + l = [0.0, 0.0], + u = [1.0, 0.0], + eval_f = (ret, x) -> begin + ret[1] = x[2]^2 + ret[2] = x[3]^2 + x[4]^3 - x[1] + return + end, + jacobian_structure = [(1, 2), (2, 1), (2, 3), (2, 4)], + eval_jacobian = (ret, x) -> begin + ret[1] = 2.0 * x[2] + ret[2] = -1.0 + ret[3] = 2.0 * x[3] + ret[4] = 3.0 * x[4]^2 + return + end, + hessian_lagrangian_structure = [(2, 2), (3, 3), (4, 4)], + eval_hessian_lagrangian = (ret, x, u) -> begin + ret[1] = 2.0 * u[1] + ret[2] = 2.0 * u[2] + ret[3] = 6.0 * x[4] * u[2] + return + end, + ) + contents = """ + VectorNonlinearOracle{Float64}(; + dimension = 3, + l = [0.0, 0.0], + u = [1.0, 0.0], + ..., + )""" + @test sprint(show, set) == contents + return +end + end # module TestSets.runtests()