Add VectorNonlinearOracle set

Author: odow

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
 

docs/src/reference/standard_form.md

@@ -104,6 +104,7 @@ ACTIVATE_ON_ONE
 Complements
 HyperRectangle
 Scaled
+VectorNonlinearOracle
 ```
 
 ## Constraint programming sets

src/Test/test_basic_constraint.jl

@@ -175,6 +175,27 @@ 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}
+    return 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,
+    )
+end
+
 function _test_function_modification(
     model::MOI.ModelLike,
     config::Config{T},

src/Test/test_nonlinear.jl

@@ -2223,3 +2223,141 @@ 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 = [0.0, 0.0],
+        u = [0.0, 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] = 2 * x[1]
+            ret[2] = 2 * x[2]
+            ret[3] = 3 * x[3]^2
+            ret[4] = -1.0
+            ret[5] = -1.0
+            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] = 2 * u[1]
+            ret[2] = 2 * u[2]
+            ret[3] = 6 * x[3] * u[2]
+            return
+        end,
+    )
+    @test MOI.dimension(set) == 5
+    x, y = MOI.add_variables(model, 3), MOI.add_variables(model, 2)
+    MOI.add_constraints.(model, x, MOI.EqualTo.(1.0:3.0))
+    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(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, 13],
+            (MOI.VectorOfVariables, MOI.VectorNonlinearOracle{T}) =>
+                zeros(T, 5),
+        ),
+    )
+    return
+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 = [0.0, 0.0],
+        u = [0.0, 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] = 2 * x[1]
+            ret[2] = 2 * x[2]
+            ret[3] = 3 * x[3]^2
+            ret[4] = -1.0
+            ret[5] = -1.0
+            return
+        end,
+    )
+    @test MOI.dimension(set) == 5
+    x, y = MOI.add_variables(model, 3), MOI.add_variables(model, 2)
+    MOI.add_constraints.(model, x, MOI.EqualTo.(1.0:3.0))
+    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(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, 13],
+            (MOI.VectorOfVariables, MOI.VectorNonlinearOracle{T}) =>
+                zeros(T, 5),
+        ),
+    )
+    return
+end
+
+version_added(::typeof(test_vector_nonlinear_oracle_no_hessian)) = v"1.46.0"

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),

src/sets.jl

@@ -2666,6 +2666,177 @@ 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(;
+    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(DimenionMismatch())
+        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.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(

test/General/sets.jl

@@ -419,6 +419,21 @@ function test_update_dimension()
     return
 end
 
+function test_VectorNonlinearOracle()
+    @test_throws(
+        DimenionMismatch,
+        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,
+        ),
+    )
+    return
+end
+
 end  # module
 
 TestSets.runtests()