Add support for VectorNonlinearFunction (#82)

Project.toml

@@ -12,5 +12,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 DataDeps = "0.7"
-MathOptInterface = "0.10, 1"
+MathOptInterface = "1.19"
 julia = "1.6"

src/MOI_wrapper.jl

@@ -11,8 +11,8 @@ MOI.Utilities.@model(
  (), # Typed vector sets
  (), # Scalar functions
  (), # Typed scalar functions
- (), # Vector functions
- (MOI.VectorAffineFunction,), # Typed vector functions
+ (MOI.VectorOfVariables, MOI.VectorNonlinearFunction), # Vector functions
+ (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction), # Typed vector functions
  true, # is_optimizer
 )
 
@@ -220,6 +220,88 @@ function _F_linear_operator(model::Optimizer)
  return M, q, SparseArrays.nnz(M)
 end
 
+_to_f(f) = convert(MOI.ScalarNonlinearFunction, f)
+
+_to_x(f::MOI.VariableIndex) = f
+
+_to_x(f::MOI.ScalarAffineFunction) = convert(MOI.VariableIndex, f)
+
+_to_x(f::MOI.ScalarQuadraticFunction) = convert(MOI.VariableIndex, f)
+
+function _to_x(f::MOI.ScalarNonlinearFunction)
+ # Hacky way to ensure that f is a standalone variable
+ @assert f isa MOI.ScalarNonlinearFunction
+ @assert f.head == :+ && length(f.args) == 1
+ @assert f.args[1] isa MOI.VariableIndex
+ return return f.args[1]
+end
+
+function _F_nonlinear_operator(model::Optimizer)
+ x = MOI.get(model, MOI.ListOfVariableIndices())
+ f_map = Vector{MOI.ScalarNonlinearFunction}(undef, length(x))
+ for (FType, SType) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
+ if SType != MOI.Complements
+ continue
+ end
+ for ci in MOI.get(model, MOI.ListOfConstraintIndices{FType,SType}())
+ f = MOI.get(model, MOI.ConstraintFunction(), ci)
+ s = MOI.get(model, MOI.ConstraintSet(), ci)
+ N = div(MOI.dimension(s), 2)
+ scalars = MOI.Utilities.scalarize(f)
+ for i in 1:N
+ fi, xi = _to_f(scalars[i]), _to_x(scalars[i+N])
+ f_map[xi.value] = fi
+ end
+ end
+ end
+ nlp = MOI.Nonlinear.Model()
+ for fi in f_map
+ MOI.Nonlinear.add_constraint(nlp, fi, MOI.EqualTo(0.0))
+ end
+ evaluator =
+ MOI.Nonlinear.Evaluator(nlp, MOI.Nonlinear.SparseReverseMode(), x)
+ MOI.initialize(evaluator, [:Jac])
+ J_structure = MOI.jacobian_structure(evaluator)
+ forward_perm = sortperm(J_structure; by = reverse)
+ inverse_perm = invperm(forward_perm)
+ jacobian_called = false
+ function F(::Cint, x::Vector{Cdouble}, f::Vector{Cdouble})
+ MOI.eval_constraint(evaluator, f, x)
+ return Cint(0)
+ end
+ function J(
+ ::Cint,
+ nnz::Cint,
+ x::Vector{Cdouble},
+ col::Vector{Cint},
+ len::Vector{Cint},
+ row::Vector{Cint},
+ data::Vector{Cdouble},
+ )
+ if !jacobian_called
+ k = 1
+ last_col = 0
+ # We need to zero all entries up front in case some rows do not
+ # appear in the Jacobian.
+ len .= Cint(0)
+ for p in forward_perm
+ r, c = J_structure[p]
+ if c != last_col
+ col[c], last_col = k, c
+ end
+ len[c] += 1
+ row[k] = r
+ k += 1
+ end
+ jacobian_called = true
+ @assert nnz == k - 1
+ end
+ MOI.eval_constraint_jacobian(evaluator, view(data, inverse_perm), x)
+ return Cint(0)
+ end
+ return F, J, length(J_structure)
+end
+
 _finite(x, y) = isfinite(x) ? x : y
 
 function _bounds_and_starting(model::Optimizer)
@@ -238,12 +320,20 @@ function _bounds_and_starting(model::Optimizer)
 end
 
 function MOI.optimize!(model::Optimizer)
+ con_types = MOI.get(model, MOI.ListOfConstraintTypesPresent())
+ is_nlp =
+ (MOI.VectorNonlinearFunction, MOI.Complements) in con_types ||
+ (MOI.VectorQuadraticFunction{Float64}, MOI.Complements) in con_types
+ F, J, nnz = if is_nlp
+ _F_nonlinear_operator(model)
+ else
+ _F_linear_operator(model)
+ end
  model.ext[:solution] = nothing
  lower, upper, initial = _bounds_and_starting(model)
- M, q, nnz = _F_linear_operator(model)
  status, x, info = solve_mcp(
- M,
- q,
+ F,
+ J,
  lower,
  upper,
  initial;

test/MOI_wrapper.jl

@@ -364,6 +364,33 @@ function test_supports()
  return
 end
 
+function test_VectorQuadraticFunction_VectorNonlinearFunction()
+ model = PATHSolver.Optimizer()
+ MOI.Utilities.loadfromstring!(
+ model,
+ """
+ variables: v, w, x, y, z
+ VectorNonlinearFunction([0, v]) in Complements(2)
+ [-1.0 * y * y + -1.0 * z + 2, w] in Complements(2)
+ VectorNonlinearFunction([ScalarNonlinearFunction(y^3 - 2z^2 + 2), ScalarNonlinearFunction(w^5 - x + 2y - 2z - 2), x, y]) in Complements(4)
+ VectorNonlinearFunction([ScalarNonlinearFunction(w + 2x^3 - 2y + 4z - 6), z]) in Complements(2)
+ v in EqualTo(1.0)
+ w in Interval(0.0, 10.0)
+ x in GreaterThan(0.0)
+ z in Interval(1.0, 2.0)
+ """,
+ )
+ x = MOI.get(model, MOI.ListOfVariableIndices())
+ MOI.set.(model, MOI.VariablePrimalStart(), x, 1.0)
+ MOI.optimize!(model)
+ @test ≈(
+ MOI.get.(model, MOI.VariablePrimal(), x),
+ [1.0, 1.2847523, 0.9729164, 0.9093761, 1.1730350];
+ atol = 1e-6,
+ )
+ return
+end
+
 end
 
 TestMOIWrapper.runtests()