From 433325219e0410900c9035e817699d2bdb8377b6 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Wed, 16 Aug 2023 19:56:48 +1200 Subject: [PATCH] Add support for VectorNonlinearFunction (#82) --- Project.toml | 2 +- src/MOI_wrapper.jl | 100 +++++++++++++++++++++++++++++++++++++++++--- test/MOI_wrapper.jl | 27 ++++++++++++ 3 files changed, 123 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index fdbdec3..7a9dcf9 100644 --- a/Project.toml +++ b/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" diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index ce494fd..953b096 100644 --- a/src/MOI_wrapper.jl +++ b/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; diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 4b0621d..80c58e0 100644 --- a/test/MOI_wrapper.jl +++ b/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()