Override MOI.Nonlinear.model for ArrayDiff.Mode

.github/workflows/ci.yml

@@ -25,6 +25,13 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v3
+      - name: MOI
+        shell: julia --project=@. {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="MathOptInterface", rev="bl/linearity"),
+          ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         with:

src/ArrayDiff.jl

@@ -86,10 +86,23 @@ function from_onnx end
 
 model(::Mode{S}) where {S} = Model{eltype(S)}()
 
+# Hook so that solvers using `MOI.Nonlinear.model(backend)` (for example,
+# NLopt and NLPModelsJuMP) receive an ArrayDiff model for an ArrayDiff mode.
+# ArrayDiff handles scalar nonlinear functions. The outer layers own quadratic
+# functions, variables, bounds, parameters, and vector nonlinear oracles.
+function Nonlinear.model(mode::Mode)
+    return Nonlinear.ModelWithQuad(
+        Nonlinear.ModelWithOracles(model(mode)),
+    )
+end
+
 # Extend MOI.Nonlinear.set_objective so that solvers calling
 # MOI.Nonlinear.set_objective(arraydiff_model, snf) dispatch here.
 function Nonlinear.set_objective(model::Model, obj::MOI.ScalarNonlinearFunction)
     model.objective = parse_expression(model, obj)
+    if model.objective_sense == MOI.FEASIBILITY_SENSE
+        model.objective_sense = MOI.MIN_SENSE
+    end
     return
 end
 
@@ -98,6 +111,99 @@ function Nonlinear.set_objective(model::Model, ::Nothing)
     return
 end
 
+Nonlinear._parameter_values(model::Model) = model.parameters
+Nonlinear._has_nonlinear_data(model::Model) =
+    model.objective !== nothing || !isempty(model.constraints)
+Nonlinear._is_nonlinear_input(
+    ::Model{T},
+    ::MOI.ScalarNonlinearFunction,
+    ::Union{
+        MOI.LessThan{T},
+        MOI.GreaterThan{T},
+        MOI.EqualTo{T},
+        MOI.Interval{T},
+    },
+) where {T} = true
+Nonlinear._is_nonlinear_objective(
+    ::Model,
+    ::MOI.ScalarNonlinearFunction,
+) = true
+
+MOI.supports_incremental_interface(::Model) = true
+MOI.supports(::Model, ::MOI.ObjectiveSense) = true
+MOI.get(model::Model, ::MOI.ObjectiveSense) = model.objective_sense
+function MOI.set(model::Model, ::MOI.ObjectiveSense, sense)
+    model.objective_sense = sense
+    return
+end
+MOI.supports(
+    ::Model,
+    ::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction},
+) = true
+function MOI.set(
+    model::Model,
+    ::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction},
+    f::MOI.ScalarNonlinearFunction,
+)
+    return Nonlinear.set_objective(model, f)
+end
+
+function MOI.supports_constraint(
+    ::Model{T},
+    ::Type{MOI.ScalarNonlinearFunction},
+    ::Type{S},
+) where {
+    T,
+    S<:Union{
+        MOI.LessThan{T},
+        MOI.GreaterThan{T},
+        MOI.EqualTo{T},
+        MOI.Interval{T},
+    },
+}
+    return true
+end
+function MOI.add_constraint(
+    model::Model{T},
+    f::MOI.ScalarNonlinearFunction,
+    s::S,
+) where {
+    T,
+    S<:Union{
+        MOI.LessThan{T},
+        MOI.GreaterThan{T},
+        MOI.EqualTo{T},
+        MOI.Interval{T},
+    },
+}
+    ci = add_constraint(model, f, s)
+    return MOI.ConstraintIndex{typeof(f),S}(ci.value)
+end
+function MOI.is_valid(
+    model::Model{T},
+    ci::MOI.ConstraintIndex{MOI.ScalarNonlinearFunction,S},
+) where {
+    T,
+    S<:Union{
+        MOI.LessThan{T},
+        MOI.GreaterThan{T},
+        MOI.EqualTo{T},
+        MOI.Interval{T},
+    },
+}
+    return haskey(model.constraints, ConstraintIndex(ci.value))
+end
+
+function Nonlinear.constraint_rows(
+    model::Model,
+    ci::MOI.ConstraintIndex{MOI.ScalarNonlinearFunction,<:MOI.AbstractScalarSet},
+)
+    row = findfirst(==(ConstraintIndex(ci.value)), keys(model.constraints))
+    return [something(row)]
+end
+Nonlinear.constraint_dual_starts(model::Model{T}) where {T} =
+    fill(nothing, length(model.constraints))
+
 # Create an ArrayDiff Evaluator from an ArrayDiff Model.
 function Evaluator(
     model::ArrayDiff.Model,
@@ -118,6 +224,12 @@ function Nonlinear.Evaluator(
     return Evaluator(model, mode, ordered_variables)
 end
 
+function Nonlinear._constraint_bounds(evaluator::Evaluator)
+    return [_bound(c.set) for (_, c) in evaluator.model.constraints]
+end
+Nonlinear._has_objective(evaluator::Evaluator) =
+    evaluator.model.objective !== nothing
+
 include("JuMP/JuMP.jl")
 
 end  # module

src/evaluator.jl

@@ -30,14 +30,16 @@ end
 
 function MOI.eval_objective(evaluator::Evaluator, x)
     start = time()
-    obj = MOI.eval_objective(evaluator.backend, x)
+    obj = Nonlinear._objective_sign(evaluator.model.objective_sense) *
+          MOI.eval_objective(evaluator.backend, x)
     evaluator.eval_objective_timer += time() - start
     return obj
 end
 
 function MOI.eval_objective_gradient(evaluator::Evaluator, g, x)
     start = time()
     MOI.eval_objective_gradient(evaluator.backend, g, x)
+    g .*= Nonlinear._objective_sign(evaluator.model.objective_sense)
     evaluator.eval_objective_gradient_timer += time() - start
     return
 end
@@ -113,7 +115,8 @@ end
 
 function MOI.eval_hessian_lagrangian(evaluator::Evaluator, H, x, σ, μ)
     start = time()
-    MOI.eval_hessian_lagrangian(evaluator.backend, H, x, σ, μ)
+    sign = Nonlinear._objective_sign(evaluator.model.objective_sense)
+    MOI.eval_hessian_lagrangian(evaluator.backend, H, x, sign * σ, μ)
     evaluator.eval_hessian_lagrangian_timer += time() - start
     return
 end
@@ -146,7 +149,15 @@ function MOI.eval_hessian_lagrangian_product(
     μ,
 )
     start = time()
-    MOI.eval_hessian_lagrangian_product(evaluator.backend, H, x, v, σ, μ)
+    sign = Nonlinear._objective_sign(evaluator.model.objective_sense)
+    MOI.eval_hessian_lagrangian_product(
+        evaluator.backend,
+        H,
+        x,
+        v,
+        sign * σ,
+        μ,
+    )
     evaluator.eval_hessian_lagrangian_timer += time() - start
     return
 end

src/model.jl

@@ -3,6 +3,9 @@
 
 function set_objective(model::Model, obj)
     model.objective = parse_expression(model, obj)
+    if model.objective_sense == MOI.FEASIBILITY_SENSE
+        model.objective_sense = MOI.MIN_SENSE
+    end
     return
 end
 

src/types.jl

@@ -205,7 +205,7 @@ It has the following fields:
  * `parameters::Vector{Float64}` : holds the current values of the parameters.
  * `operators::OperatorRegistry` : stores the operators used in the model.
 """
-mutable struct Model{T}
+mutable struct Model{T} <: MOI.ModelLike
     objective::Union{Nothing,Expression{T}}
     # Vector residual for nonlinear least-squares objectives. When set, callers
     # can query its value, `J*v`, `J'*v`, etc. via the evaluator.
@@ -214,6 +214,7 @@ mutable struct Model{T}
     constraints::OrderedCollections.OrderedDict{ConstraintIndex,Constraint{T}}
     parameters::Vector{T}
     operators::OperatorRegistry
+    objective_sense::MOI.OptimizationSense
     # This is a private field, used only to increment the ConstraintIndex.
     last_constraint_index::Int64
     function Model{T}() where {T}
@@ -224,6 +225,7 @@ mutable struct Model{T}
             OrderedCollections.OrderedDict{ConstraintIndex,Constraint{T}}(),
             T[],
             OperatorRegistry(),
+            MOI.FEASIBILITY_SENSE,
             0,
         )
     end

test/IpoptBackend.jl

@@ -0,0 +1,60 @@
+module TestIpoptBackend
+
+using Test
+import ArrayDiff
+import Ipopt
+import MathOptInterface as MOI
+
+function _square_minus(x, value)
+    difference = MOI.ScalarNonlinearFunction(:-, Any[x, value])
+    return MOI.ScalarNonlinearFunction(:^, Any[difference, 2.0])
+end
+
+function test_ipopt_with_arraydiff_mode()
+    model = Ipopt.Optimizer()
+    MOI.set(model, MOI.Silent(), true)
+    mode = ArrayDiff.Mode()
+    MOI.set(model, MOI.AutomaticDifferentiationBackend(), mode)
+    @test MOI.get(model, MOI.AutomaticDifferentiationBackend()) === mode
+
+    x, y = MOI.add_variables(model, 2)
+    objective = MOI.ScalarNonlinearFunction(
+        :+,
+        Any[_square_minus(x, 1.0), _square_minus(y, 2.0)],
+    )
+    MOI.set(model, MOI.ObjectiveFunction{typeof(objective)}(), objective)
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+
+    quadratic = MOI.ScalarQuadraticFunction(
+        [
+            MOI.ScalarQuadraticTerm(2.0, x, x),
+            MOI.ScalarQuadraticTerm(2.0, y, y),
+        ],
+        MOI.ScalarAffineTerm{Float64}[],
+        0.0,
+    )
+    MOI.add_constraint(model, quadratic, MOI.LessThan(10.0))
+
+    oracle = MOI.VectorNonlinearOracle(;
+        dimension = 2,
+        l = [3.0],
+        u = [Inf],
+        eval_f = (output, input) -> (output[1] = input[1] + input[2]),
+        jacobian_structure = [(1, 1), (1, 2)],
+        eval_jacobian = (values, input) -> (values .= 1.0),
+        hessian_lagrangian_structure = Tuple{Int,Int}[],
+        eval_hessian_lagrangian = (values, input, multipliers) -> nothing,
+    )
+    MOI.add_constraint(model, MOI.VectorOfVariables([x, y]), oracle)
+
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED
+    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1.0 atol = 1e-4
+    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 2.0 atol = 1e-4
+    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol = 1e-8
+    return
+end
+
+test_ipopt_with_arraydiff_mode()
+
+end  # module

test/runtests.jl

@@ -2,6 +2,7 @@ include("ReverseAD.jl")
 include("ArrayDiff.jl")
 include("JuMP.jl")
 include("MathOptAI.jl")
+include("IpoptBackend.jl")
 include("ONNXExt.jl")
 if VERSION >= v"1.11"
     # [sources] not supported on Julia v1.10