From 07d8e1f901d57c785f66fe791efa7a0072013fe0 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 00:07:46 -0400 Subject: [PATCH 01/21] Support quadratic constraints --- src/moi_nlp_model.jl | 159 ++++++++++++++++++++++++++++--------- src/moi_nls_model.jl | 2 +- src/utils.jl | 131 ++++++++++++++++++++++++++---- test/nlp_problems/hs100.jl | 24 ++++++ test/nlp_problems/hs219.jl | 16 ++++ test/nlp_problems/hs61.jl | 12 +++ test/runtests.jl | 2 +- test/test_moi_nlp_model.jl | 30 ++++++- 8 files changed, 322 insertions(+), 54 deletions(-) create mode 100644 test/nlp_problems/hs100.jl create mode 100644 test/nlp_problems/hs219.jl create mode 100644 test/nlp_problems/hs61.jl diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index bd86ef9..d443f37 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -4,6 +4,7 @@ mutable struct MathOptNLPModel <: AbstractNLPModel{Float64, Vector{Float64}} meta::NLPModelMeta{Float64, Vector{Float64}} eval::MOI.Nonlinear.Evaluator lincon::LinearConstraints + quadcon::QuadraticConstraints nlcon::NonLinearStructure obj::Objective counters::Counters @@ -27,7 +28,7 @@ end function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = "Generic") index_map, nvar, lvar, uvar, x0 = parser_variables(moimodel) - nlin, lincon, lin_lcon, lin_ucon = parser_MOI(moimodel, index_map) + nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon = parser_MOI(moimodel, index_map, nvar) nlp_data = _nlp_block(moimodel) nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) @@ -38,11 +39,11 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = obj = parser_objective_MOI(moimodel, nvar, index_map) end - ncon = nlin + nnln - lcon = vcat(lin_lcon, nl_lcon) - ucon = vcat(lin_ucon, nl_ucon) - nnzj = lincon.nnzj + nlcon.nnzj - nnzh = obj.nnzh + nlcon.nnzh + ncon = nlin + quadcon.nquad + nnln + lcon = vcat(lin_lcon, quad_lcon, nl_lcon) + ucon = vcat(lin_ucon, quad_ucon, nl_ucon) + nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj + nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh meta = NLPModelMeta( nvar, @@ -57,13 +58,13 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = nnzh = nnzh, lin = collect(1:nlin), lin_nnzj = lincon.nnzj, - nln_nnzj = nlcon.nnzj, + nln_nnzj = quadcon.nnzj + nlcon.nnzj, minimize = MOI.get(moimodel, MOI.ObjectiveSense()) == MOI.MIN_SENSE, - islp = (obj.type == "LINEAR") && (nnln == 0), + islp = (obj.type == "LINEAR") && (nnln == 0) && (quadcon.nquad == 0), name = name, ) - return MathOptNLPModel(meta, nlp_data.evaluator, lincon, nlcon, obj, Counters()), index_map + return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, obj, Counters()), index_map end function NLPModels.obj(nlp::MathOptNLPModel, x::AbstractVector) @@ -106,7 +107,13 @@ end function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVector) increment!(nlp, :neval_cons_nln) - MOI.eval_constraint(nlp.eval, c, x) + for i = 1:(nlp.quadcon.nquad) + qcon = nlp.quadcon.constaints[i] + c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x) + end + if nlp.meta.nnln > nlp.quadcon.nquad + MOI.eval_constraint(nlp.eval, view(c, (nlp.quadcon.nquad + 1):(nlp.meta.nnln)), x) + end return c end @@ -125,8 +132,20 @@ function NLPModels.jac_nln_structure!( rows::AbstractVector{<:Integer}, cols::AbstractVector{<:Integer}, ) - view(rows, 1:(nlp.nlcon.nnzj)) .= nlp.nlcon.jac_rows - view(cols, 1:(nlp.nlcon.nnzj)) .= nlp.nlcon.jac_cols + index = 0 + if nlp.quadcon.nquad > 0 + for i = 1:(nlp.quadcon.nquad) + qcon = nlp.quadcon.constraints[i] + view(rows, index+1:index+qcon.nnzg) .= i + view(cols, index+1:index+qcon.nnzg) .= qcon.g + index += qcon.nnzg + end + end + if nlp.meta.nnln > nlp.quadcon.nquad + ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) + view(rows, ind_nnln) .= nlp.nlcon.jac_rows + view(cols, ind_nnln) .= nlp.nlcon.jac_cols + end return rows, cols end @@ -138,7 +157,33 @@ end function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::AbstractVector) increment!(nlp, :neval_jac_nln) - MOI.eval_constraint_jacobian(nlp.eval, view(vals, 1:(nlp.nlcon.nnzj)), x) + index = 0 + if nlp.quadcon.nquad > 0 + view(vals, 1:nlp.quadcon.nnzj) .= 0.0 + for i = 1:(nlp.quadcon.nquad) + qcon = nlp.quadcon.constraints[i] + for (j, ind) in enumerate(qcon.b.nzind) + k = qcon.dg[ind] + vals[index+k] += qcon.b.nzval[j] + end + for j = 1:qcon.nnzh + row = qcon.A.rows[j] + col = qcon.A.cols[j] + val = qcon.A.vals[j] + k1 = qcon.dg[row] + vals[index+k1] += val * x[col] + if row != col + k2 = qcon.dg[col] + vals[index+k2] += val * x[row] + end + end + index += qcon.nnzg + end + end + if nlp.meta.nnln > nlp.quadcon.nquad + ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) + MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) + end return vals end @@ -166,7 +211,18 @@ function NLPModels.jprod_nln!( Jv::AbstractVector, ) increment!(nlp, :neval_jprod_nln) - MOI.eval_constraint_jacobian_product(nlp.eval, Jv, x, v) + if nlp.meta.nnln > nlp.quadcon.nquad + ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) + MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v) + end + (nlp.meta.nnln == nlp.quadcon.nquad) && (Jv .= 0.0) + if nlp.quadcon.nquad > 0 + for i = 1:(nlp.quadcon.nquad) + # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v + qcon = nlp.quadcon.constraints[i] + v[i] += coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) + end + end return Jv end @@ -194,7 +250,19 @@ function NLPModels.jtprod_nln!( Jtv::AbstractVector, ) increment!(nlp, :neval_jtprod_nln) - MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, v) + if nlp.meta.nnln > nlp.quadcon.nquad + ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) + MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln)) + end + (nlp.meta.nnln == nlp.quadcon.nquad) && (Jtv .= 0.0) + if nlp.quadcon.nquad > 0 + for i = 1:(nlp.quadcon.nquad) + # Jtv += v[i] * (Aᵢ * x + bᵢ) + qcon = nlp.quadcon.constraints[i] + coo_sym_add_mul!(rows, cols, vals, x, Jtv, v[i]) + Jtv .+= v[i] .* qcon.b + end + end return Jtv end @@ -207,9 +275,18 @@ function NLPModels.hess_structure!( view(rows, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.rows view(cols, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.cols end - if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0) - view(rows, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_rows - view(cols, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_cols + if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) + view(rows, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_rows + view(cols, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_cols + end + if nlp.quadcon.nquad > 0 + index = nlp.obj.nnzh + for i = 1:(nlp.quadcon.nquad) + qcon = nlp.quadcon.constraints[i] + view(rows, (index + 1):(index + qcon.nnzh)) .= qcon.A.rows + view(cols, (index + 1):(index + qcon.nnzh)) .= qcon.A.cols + index += qcon.nnzh + end end return rows, cols end @@ -225,15 +302,22 @@ function NLPModels.hess_coord!( if nlp.obj.type == "QUADRATIC" view(vals, 1:(nlp.obj.nnzh)) .= obj_weight .* nlp.obj.hessian.vals end - if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0) - MOI.eval_hessian_lagrangian( - nlp.eval, - view(vals, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)), + if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) + MOI.eval_hessian_lagrangian(nlp.eval, + view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)), x, obj_weight, - view(y, nlp.meta.nln), + view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nnln)) ) end + if nlp.quadcon.nquad > 0 + index = nlp.obj.nnzh + for i = 1:(nlp.quadcon.nquad) + qcon = nlp.quadcon.constraints[i] + view(vals, (index + 1):(index + qcon.nnzh)) .= y[i] .* qcon.A.vals + index += qcon.nnzh + end + end return vals end @@ -252,7 +336,8 @@ function NLPModels.hess_coord!( view(vals, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= 0.0 end if nlp.obj.type == "NONLINEAR" - MOI.eval_hessian_lagrangian(nlp.eval, vals, x, obj_weight, zeros(nlp.meta.nnln)) + λ = zeros(nlp.meta.nnln - nlp.quadcon.nquad) # Should be stored in the structure MathOptNLPModel + MOI.eval_hessian_lagrangian(nlp.eval, vals, x, obj_weight, λ) end return vals end @@ -269,19 +354,20 @@ function NLPModels.hprod!( if (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == 0) hv .= 0.0 end - if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0) - MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, view(y, nlp.meta.nln)) + if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) + λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)) + MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, λ) end if nlp.obj.type == "QUADRATIC" - nlp.meta.nnln == 0 && (hv .= 0.0) - coo_sym_add_mul!( - nlp.obj.hessian.rows, - nlp.obj.hessian.cols, - nlp.obj.hessian.vals, - v, - hv, - obj_weight, - ) + (nlp.meta.nnln == nlp.quadcon.nquad) && (hv .= 0.0) + coo_sym_add_mul!(nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals, v, hv, obj_weight) + end + if nlp.quadcon.nquad > 0 + (nlp.obj.type == "LINEAR") && (hv .= 0.0) + for i = 1:(nlp.quadcon.nquad) + qcon = nlp.quadcon.constraints[i] + coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, y[nlp.meta.nlin + i]) + end end return hv end @@ -302,7 +388,8 @@ function NLPModels.hprod!( coo_sym_add_mul!(nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals, v, hv, obj_weight) end if nlp.obj.type == "NONLINEAR" - MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, zeros(nlp.meta.nnln)) + λ = zeros(nlp.meta.nnln - nlp.quadcon.nquad) # Should be stored in the structure MathOptNLPModel + MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, λ) end return hv end diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl index 0c1cd8d..b0e4dda 100644 --- a/src/moi_nls_model.jl +++ b/src/moi_nls_model.jl @@ -30,7 +30,7 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri _nlp_sync!(cmodel) moimodel = backend(cmodel) - nlin, lincon, lin_lcon, lin_ucon = parser_MOI(moimodel, index_map) + nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon = parser_MOI(moimodel, index_map, nvar) nlp_data = _nlp_block(moimodel) nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) diff --git a/src/utils.jl b/src/utils.jl index e1cf1c6..443afa5 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -6,10 +6,18 @@ import NLPModels.increment!, NLPModels.decrement! using JuMP, MathOptInterface const MOI = MathOptInterface +# VariableIndex +const VI = MOI.VariableIndex # VariableIndex(value) + # ScalarAffineFunctions and VectorAffineFunctions -const SAF = MOI.ScalarAffineFunction{Float64} -const VAF = MOI.VectorAffineFunction{Float64} -const AF = Union{SAF, VAF} +const SAF = MOI.ScalarAffineFunction{Float64} # ScalarAffineFunction{T}(terms, constant) +const VAF = MOI.VectorAffineFunction{Float64} # VectorAffineFunction{T}(terms, constants) +const AF = Union{SAF, VAF} + +# ScalarQuadraticFunctions and VectorQuadraticFunctions +const SQF = MOI.ScalarQuadraticFunction{Float64} # ScalarQuadraticFunction{T}(affine_terms, quadratic_terms, constant) +const VQF = MOI.VectorQuadraticFunction{Float64} # VectorQuadraticFunction{T}(affine_terms, quadratic_terms, constants) +const QF = Union{SQF, VQF} # AffLinSets and VecLinSets const ALS = Union{ @@ -21,10 +29,6 @@ const ALS = Union{ const VLS = Union{MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros} const LS = Union{ALS, VLS} -const VI = MOI.VariableIndex -const SQF = MOI.ScalarQuadraticFunction{Float64} -const LinQuad = Union{VI, SAF, SQF} - # Expressions const VF = VariableRef const AE = GenericAffExpr{Float64, VariableRef} @@ -32,6 +36,9 @@ const LE = Union{VF, AE} const QE = GenericQuadExpr{Float64, VariableRef} const NLE = NonlinearExpression +const LinQuad = Union{VI, SAF, SQF} + +# Sparse matrix in coordinate format mutable struct COO rows::Vector{Int} cols::Vector{Int} @@ -45,6 +52,23 @@ mutable struct LinearConstraints nnzj::Int end +# xᵀAx + bᵀx +mutable struct QuadraticConstraint + A::COO + b::SparseVector{Float64} + g::Vector{Int} + dg::Dict{Int,Int} + nnzg::Int + nnzh::Int +end + +mutable struct QuadraticConstraints + nquad::Int + constraints::Vector{QuadraticConstraint} + nnzj::Int + nnzh::Int +end + mutable struct NonLinearStructure jac_rows::Vector{Int} jac_cols::Vector{Int} @@ -200,11 +224,69 @@ function parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_uco end """ - parser_MOI(moimodel, index_map) + parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) + +Parse a `ScalarQuadraticFunction` fun with its associated set. +`qcons`, `quad_lcon`, `quad_ucon` are updated. +""" +function parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) + _index(v::MOI.VariableIndex) = index_map[v].value + + b = spzeros(Float64, nvar) + rows = Int[] + cols = Int[] + vals = Float64[] + + # Parse a ScalarAffineTerm{Float64}(coefficient, variable_index) + for term in fun.affine_terms + b[_index(term.variable)] = term.coefficient + end + + # Parse a ScalarQuadraticTerm{Float64}(coefficient, variable_index_1, variable_index_2) + for term in fun.quadratic_terms + i = _index(term.variable_1) + j = _index(term.variable_2) + if i ≥ j + push!(rows, i) + push!(cols, j) + else + push!(rows, j) + push!(cols, i) + end + push!(vals, term.coefficient) + end + + if typeof(set) in (MOI.Interval{Float64}, MOI.GreaterThan{Float64}) + push!(quad_lcon, -fun.constant + set.lower) + elseif typeof(set) == MOI.EqualTo{Float64} + push!(quad_lcon, -fun.constant + set.value) + else + push!(quad_lcon, -Inf) + end + + if typeof(set) in (MOI.Interval{Float64}, MOI.LessThan{Float64}) + push!(quad_ucon, -fun.constant + set.upper) + elseif typeof(set) == MOI.EqualTo{Float64} + push!(quad_ucon, -fun.constant + set.value) + else + push!(quad_ucon, Inf) + end + + A = COO(rows, cols, vals) + g = unique(vcat(rows, b.nzind)) # sparsity pattern of Ax + b + nnzg = length(g) + dg = Dict{Int,Int}(g[i] => i for i = 1:nnzg) + nnzh = length(vals) + qcon = QuadraticConstraint(A, b, g, dg, nnzg, nnzh) + push!(qcons, qcon) +end + +""" + parser_MOI(moimodel, index_map, nvar) Parse linear constraints of a `MOI.ModelLike`. """ -function parser_MOI(moimodel, index_map) +function parser_MOI(moimodel, index_map, nvar) # Variables associated to linear constraints nlin = 0 @@ -214,9 +296,16 @@ function parser_MOI(moimodel, index_map) lin_lcon = Float64[] lin_ucon = Float64[] + # Variables associated to quadratic constraints + nquad = 0 + qcons = QuadraticConstraint[] + quad_lcon = Float64[] + quad_ucon = Float64[] + contypes = MOI.get(moimodel, MOI.ListOfConstraintTypesPresent()) for (F, S) in contypes - F <: AF || F == MOI.ScalarNonlinearFunction || F == VI || @warn("Function $F is not supported.") + F <: AF || F <: SQF || F == MOI.ScalarNonlinearFunction || F == VI || @warn("Function $F is not supported.") + # F <: AF || F <: QF || F == MOI.ScalarNonlinearFunction || F == VI || @warn("Function $F is not supported.") S <: LS || @warn("Set $S is not supported.") conindices = MOI.get(moimodel, MOI.ListOfConstraintIndices{F, S}()) @@ -237,13 +326,25 @@ function parser_MOI(moimodel, index_map) parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon, index_map) nlin += set.dimension end + if typeof(fun) <: SQF + parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) + nquad += 1 + end + # if typeof(fun) <: VQF + # parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) + # nquad += set.dimension + # end end end coo = COO(linrows, lincols, linvals) - nnzj = length(linvals) - lincon = LinearConstraints(coo, nnzj) + lin_nnzj = length(linvals) + lincon = LinearConstraints(coo, lin_nnzj) + + quad_nnzj = sum(qcons[i].nnzg for i = 1:nquad, init=0) + quad_nnzh = sum(qcons[i].nnzh for i = 1:nquad, init=0) + quadcon = QuadraticConstraints(nquad, qcons, quad_nnzj, quad_nnzh) - return nlin, lincon, lin_lcon, lin_ucon + return nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon end # Affine or quadratic, nothing to do @@ -419,8 +520,8 @@ function parser_objective_MOI(moimodel, nvar, index_map) push!(rows, i) push!(cols, j) else - push!(cols, j) - push!(rows, i) + push!(rows, j) + push!(cols, i) end push!(vals, term.coefficient) end diff --git a/test/nlp_problems/hs100.jl b/test/nlp_problems/hs100.jl new file mode 100644 index 0000000..b4af448 --- /dev/null +++ b/test/nlp_problems/hs100.jl @@ -0,0 +1,24 @@ +function hs100(args...; kwargs...) + nlp = Model() + x0 = [1, 2, 0, 4, 0, 1, 1] + @variable(nlp, x[i = 1:7], start = x0[i]) + + @NLconstraint(nlp, 127 - 2 * x[1]^2 - 3 * x[2]^4 - x[3] - 4 * x[4]^2 - 5 * x[5] ≥ 0) + @constraint(nlp, 282 - 7 * x[1] - 3 * x[2] - 10 * x[3]^2 - x[4] + x[5] ≥ 0) + @constraint(nlp, 196 - 23 * x[1] - x[2]^2 - 6 * x[6]^2 + 8 * x[7] ≥ 0) + @constraint(nlp, -4 * x[1]^2 - x[2]^2 + 3 * x[1] * x[2] - 2 * x[3]^2 - 5 * x[6] + 11 * x[7] ≥ 0) + + @NLobjective( + nlp, + Min, + (x[1] - 10)^2 + + 5 * (x[2] - 12)^2 + + x[3]^4 + + 3 * (x[4] - 11)^2 + + 10 * x[5]^6 + + 7 * x[6]^2 + + x[7]^4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7] + ) + + return nlp +end diff --git a/test/nlp_problems/hs219.jl b/test/nlp_problems/hs219.jl new file mode 100644 index 0000000..ada2147 --- /dev/null +++ b/test/nlp_problems/hs219.jl @@ -0,0 +1,16 @@ +function hs219(args...; kwargs...) + nlp = Model() + x0 = [10, 10, 10, 10] + @variable(nlp, x[i = 1:4], start = x0[i]) + + @constraint(nlp, x[1]^2 - x[2] - x[4]^2 == 0) + @NLconstraint(nlp, x[2] - x[1]^3 - x[3]^2 == 0) + + @NLobjective( + nlp, + Min, + -x[1] + ) + + return nlp +end diff --git a/test/nlp_problems/hs61.jl b/test/nlp_problems/hs61.jl new file mode 100644 index 0000000..5e422cc --- /dev/null +++ b/test/nlp_problems/hs61.jl @@ -0,0 +1,12 @@ +"HS61 model" +function hs61(args...; kwargs...) + nlp = Model() + @variable(nlp, x[i = 1:3], start = 0) + + @constraint(nlp, 3 * x[1] - 2 * x[2]^2 - 7 == 0) + @constraint(nlp, 4 * x[1] - x[3]^2 - 11 == 0) + + @NLobjective(nlp, Min, 4 * x[1]^2 + 2 * x[2]^2 + 2 * x[3]^2 - 33 * x[1] + 16 * x[2] - 24 * x[3]) + + return nlp +end diff --git a/test/runtests.jl b/test/runtests.jl index 49d13d9..eeb3fee 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,7 @@ nlp_problems = setdiff(NLPModelsTest.nlp_problems, ["MGH01Feas"]) nls_problems = NLPModelsTest.nls_problems extra_nls_problems = ["HS30", "HS43", "MGH07", "nlsnohesspb"] -for problem in lowercase.(nlp_problems ∪ ["nohesspb"]) +for problem in lowercase.(nlp_problems ∪ ["nohesspb", "hs61", "hs100", "hs219"]) include(joinpath("nlp_problems", "$problem.jl")) end diff --git a/test/test_moi_nlp_model.jl b/test/test_moi_nlp_model.jl index 162694d..8df6224 100644 --- a/test/test_moi_nlp_model.jl +++ b/test/test_moi_nlp_model.jl @@ -11,7 +11,7 @@ println("Testing MathOptNLPModel") "‖c(x₀)‖" ) # Test that every problem can be instantiated. -for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb"])) +for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb", "hs61", "hs100"])) prob_fn = eval(prob) nlp = MathOptNLPModel(prob_fn(), hessian = (prob != :nohesspb)) n = nlp.meta.nvar @@ -23,3 +23,31 @@ for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb"])) @printf("%-15s %4d %4d %10.4e %10.4e %10s\n", prob, n, m, fx, ngx, ncx) end println() + +@testset "Testing quadratic constraints with JuMP" begin + Hessian(x) = zeros(4, 4) + + function Hessian(x, y) + H = zeros(4, 4) + H[1, 1] = 2 * y[1] - y[2] * 6 * x[1] + H[3, 3] = - 2 * y[2] + H[4, 4] = - 2 * y[1] + return H + end + + Jacobian(x) = [ + 2x[1] -1 0 -2x[4]; + -3x[1]^2 1 -2x[3] 0 + ] + + jump = hs219() + nlp = MathOptNLPModel(jump) + x1 = nlp.meta.x0 + y1 = ones(nlp.meta.ncon) + v1 = 2 * ones(nlp.meta.nvar) + @test jac(nlp, x1) ≈ Jacobian(x1) + @test hess(nlp, x1) ≈ Hessian(x1) + @test hess(nlp, x1, y1) ≈ Hessian(x1, y1) + @test hprod(nlp, x1, x1) ≈ Hessian(x1) * x1 + @test hprod(nlp, x1, y1, v1) ≈ Hessian(x1, y1) * v1 +end From 4dc88cec3eeb44e8ed7d578f1e01fc774c4dbf9a Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 00:17:43 -0400 Subject: [PATCH 02/21] Update utils.jl --- src/utils.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 443afa5..7e28aea 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -340,8 +340,12 @@ function parser_MOI(moimodel, index_map, nvar) lin_nnzj = length(linvals) lincon = LinearConstraints(coo, lin_nnzj) - quad_nnzj = sum(qcons[i].nnzg for i = 1:nquad, init=0) - quad_nnzh = sum(qcons[i].nnzh for i = 1:nquad, init=0) + quad_nnzj = 0 + quad_nnzh = 0 + for i = 1:nquad + quad_nnzj += qcons[i].nnzg + quad_nnzh += qcons[i].nnzh + end quadcon = QuadraticConstraints(nquad, qcons, quad_nnzj, quad_nnzh) return nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon From 2be947ca6917066958fda31f0e3b559084a5556d Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 00:21:35 -0400 Subject: [PATCH 03/21] Fix a typo --- src/moi_nlp_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index d443f37..92fbae6 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -108,7 +108,7 @@ end function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVector) increment!(nlp, :neval_cons_nln) for i = 1:(nlp.quadcon.nquad) - qcon = nlp.quadcon.constaints[i] + qcon = nlp.quadcon.constraints[i] c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x) end if nlp.meta.nnln > nlp.quadcon.nquad From dd738ba58e982b3ff24b14a464246b6f12ec25cb Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 01:02:59 -0400 Subject: [PATCH 04/21] Add a function parser_VQF --- src/utils.jl | 74 ++++++++++++++++++++++++++++++++++---- test/test_moi_nlp_model.jl | 28 --------------- 2 files changed, 68 insertions(+), 34 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 7e28aea..519ae7d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -281,6 +281,69 @@ function parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) push!(qcons, qcon) end +""" + parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) + +Parse a `VectorQuadraticFunction` fun with its associated set. +`qcons`, `quad_lcon`, `quad_ucon` are updated. +""" +function parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) + _index(v::MOI.VariableIndex) = index_map[v].value + + ncon = length(fun.constants) + for i = 1:ncon + b = spzeros(Float64, nvar) + rows = Int[] + cols = Int[] + vals = Float64[] + + quadratic_terms = fun.quadratic_terms[i] + affine_terms = fun.affine_terms[i] + constant = fun.constants[i] + + # Parse a VectorAffineTerm{Float64}(output_index, scalar_term) + for term in affine_terms + @assert term.output_index == i + b[_index(term.scalar_term.variable)] = term.scalar_term.coefficient + end + + # Parse a VectorQuadraticTerm{Float64}(output_index, scalar_term) + for term in quadratic_terms + @assert term.output_index == i + i = _index(term.scalar_term.variable_1) + j = _index(term.scalar_term.variable_2) + if i ≥ j + push!(rows, i) + push!(cols, j) + else + push!(rows, j) + push!(cols, i) + end + push!(vals, term.scalar_term.coefficient) + end + + if typeof(set) in (MOI.Nonnegatives, MOI.Zeros) + append!(lin_lcon, constant) + else + append!(lin_lcon, -Inf) + end + + if typeof(set) in (MOI.Nonpositives, MOI.Zeros) + append!(lin_ucon, -constant) + else + append!(lin_ucon, Inf) + end + + A = COO(rows, cols, vals) + g = unique(vcat(rows, b.nzind)) # sparsity pattern of Ax + b + nnzg = length(g) + dg = Dict{Int,Int}(g[i] => i for i = 1:nnzg) + nnzh = length(vals) + qcon = QuadraticConstraint(A, b, g, dg, nnzg, nnzh) + push!(qcons, qcon) + end +end + """ parser_MOI(moimodel, index_map, nvar) @@ -304,8 +367,7 @@ function parser_MOI(moimodel, index_map, nvar) contypes = MOI.get(moimodel, MOI.ListOfConstraintTypesPresent()) for (F, S) in contypes - F <: AF || F <: SQF || F == MOI.ScalarNonlinearFunction || F == VI || @warn("Function $F is not supported.") - # F <: AF || F <: QF || F == MOI.ScalarNonlinearFunction || F == VI || @warn("Function $F is not supported.") + F <: AF || F <: QF || F == MOI.ScalarNonlinearFunction || F == VI || @warn("Function $F is not supported.") S <: LS || @warn("Set $S is not supported.") conindices = MOI.get(moimodel, MOI.ListOfConstraintIndices{F, S}()) @@ -330,10 +392,10 @@ function parser_MOI(moimodel, index_map, nvar) parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) nquad += 1 end - # if typeof(fun) <: VQF - # parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) - # nquad += set.dimension - # end + if typeof(fun) <: VQF + parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) + nquad += set.dimension + end end end coo = COO(linrows, lincols, linvals) diff --git a/test/test_moi_nlp_model.jl b/test/test_moi_nlp_model.jl index 8df6224..4f236ba 100644 --- a/test/test_moi_nlp_model.jl +++ b/test/test_moi_nlp_model.jl @@ -23,31 +23,3 @@ for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb", "hs61", "hs100"])) @printf("%-15s %4d %4d %10.4e %10.4e %10s\n", prob, n, m, fx, ngx, ncx) end println() - -@testset "Testing quadratic constraints with JuMP" begin - Hessian(x) = zeros(4, 4) - - function Hessian(x, y) - H = zeros(4, 4) - H[1, 1] = 2 * y[1] - y[2] * 6 * x[1] - H[3, 3] = - 2 * y[2] - H[4, 4] = - 2 * y[1] - return H - end - - Jacobian(x) = [ - 2x[1] -1 0 -2x[4]; - -3x[1]^2 1 -2x[3] 0 - ] - - jump = hs219() - nlp = MathOptNLPModel(jump) - x1 = nlp.meta.x0 - y1 = ones(nlp.meta.ncon) - v1 = 2 * ones(nlp.meta.nvar) - @test jac(nlp, x1) ≈ Jacobian(x1) - @test hess(nlp, x1) ≈ Hessian(x1) - @test hess(nlp, x1, y1) ≈ Hessian(x1, y1) - @test hprod(nlp, x1, x1) ≈ Hessian(x1) * x1 - @test hprod(nlp, x1, y1, v1) ≈ Hessian(x1, y1) * v1 -end From d80355fee468864327c633bbe49a405e35072a2e Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 01:09:13 -0400 Subject: [PATCH 05/21] test_moi_nlp_model.jl --- test/test_moi_nlp_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_moi_nlp_model.jl b/test/test_moi_nlp_model.jl index 4f236ba..e7a63cb 100644 --- a/test/test_moi_nlp_model.jl +++ b/test/test_moi_nlp_model.jl @@ -11,7 +11,7 @@ println("Testing MathOptNLPModel") "‖c(x₀)‖" ) # Test that every problem can be instantiated. -for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb", "hs61", "hs100"])) +for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb", "hs61", "hs100", "hs219"])) prob_fn = eval(prob) nlp = MathOptNLPModel(prob_fn(), hessian = (prob != :nohesspb)) n = nlp.meta.nvar From 9f9ba02f105437e19fed5f564343516513586f1b Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 01:17:09 -0400 Subject: [PATCH 06/21] Fix two mistakes --- src/moi_nlp_model.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 92fbae6..f0bfdb4 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -212,7 +212,7 @@ function NLPModels.jprod_nln!( ) increment!(nlp, :neval_jprod_nln) if nlp.meta.nnln > nlp.quadcon.nquad - ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) + ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln) MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v) end (nlp.meta.nnln == nlp.quadcon.nquad) && (Jv .= 0.0) @@ -251,7 +251,7 @@ function NLPModels.jtprod_nln!( ) increment!(nlp, :neval_jtprod_nln) if nlp.meta.nnln > nlp.quadcon.nquad - ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) + ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln) MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln)) end (nlp.meta.nnln == nlp.quadcon.nquad) && (Jtv .= 0.0) From ebee588384cf228f301a99c73ea6ba84a23adaaa Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 01:25:18 -0400 Subject: [PATCH 07/21] Fix another bug --- src/moi_nlp_model.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index f0bfdb4..9f709ac 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -132,8 +132,8 @@ function NLPModels.jac_nln_structure!( rows::AbstractVector{<:Integer}, cols::AbstractVector{<:Integer}, ) - index = 0 if nlp.quadcon.nquad > 0 + index = 0 for i = 1:(nlp.quadcon.nquad) qcon = nlp.quadcon.constraints[i] view(rows, index+1:index+qcon.nnzg) .= i @@ -157,8 +157,8 @@ end function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::AbstractVector) increment!(nlp, :neval_jac_nln) - index = 0 if nlp.quadcon.nquad > 0 + index = 0 view(vals, 1:nlp.quadcon.nnzj) .= 0.0 for i = 1:(nlp.quadcon.nquad) qcon = nlp.quadcon.constraints[i] @@ -307,7 +307,7 @@ function NLPModels.hess_coord!( view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)), x, obj_weight, - view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nnln)) + view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1)::(nlp.meta.ncon)) ) end if nlp.quadcon.nquad > 0 From 7d9f5a1dd6ae462797e5dcdcfed3c79493c383b0 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 15:01:44 -0400 Subject: [PATCH 08/21] Fix MathOptNLPModel --- src/moi_nlp_model.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 9f709ac..5a7bd9a 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -303,11 +303,12 @@ function NLPModels.hess_coord!( view(vals, 1:(nlp.obj.nnzh)) .= obj_weight .* nlp.obj.hessian.vals end if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) + λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)) MOI.eval_hessian_lagrangian(nlp.eval, view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)), x, obj_weight, - view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1)::(nlp.meta.ncon)) + λ ) end if nlp.quadcon.nquad > 0 From b7154317e772047ff1d6accdea18dade6e80b01a Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 15:10:44 -0400 Subject: [PATCH 09/21] Support quadratic constraints in MathOptNLSModel --- src/moi_nls_model.jl | 147 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 122 insertions(+), 25 deletions(-) diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl index b0e4dda..6b5c991 100644 --- a/src/moi_nls_model.jl +++ b/src/moi_nls_model.jl @@ -9,6 +9,7 @@ mutable struct MathOptNLSModel <: AbstractNLSModel{Float64, Vector{Float64}} linequ::LinearEquations nlequ::NonLinearStructure lincon::LinearConstraints + quadcon::QuadraticConstraints nlcon::NonLinearStructure counters::NLSCounters end @@ -39,11 +40,11 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri Fnnzj = linequ.nnzj + nlequ.nnzj Fnnzh = nlequ.nnzh - ncon = nlin + nnln - lcon = vcat(lin_lcon, nl_lcon) - ucon = vcat(lin_ucon, nl_ucon) - cnnzj = lincon.nnzj + nlcon.nnzj - cnnzh = lls.nnzh + nlcon.nnzh + ncon = nlin + quadcon.nquad + nnln + lcon = vcat(lin_lcon, quad_lcon, nl_lcon) + ucon = vcat(lin_ucon, quad_ucon, nl_ucon) + cnnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj + cnnzh = lls.nnzh + quadcon.nnzh + nlcon.nnzh meta = NLPModelMeta( nvar, @@ -58,7 +59,7 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri nnzh = cnnzh, lin = collect(1:nlin), lin_nnzj = lincon.nnzj, - nln_nnzj = nlcon.nnzj, + nln_nnzj = quadcon.nnzj + nlcon.nnzj, minimize = objective_sense(cmodel) == MOI.MIN_SENSE, islp = false, name = name, @@ -73,6 +74,7 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri linequ, nlequ, lincon, + quadcon, nlcon, NLSCounters(), ) @@ -253,7 +255,13 @@ end function NLPModels.cons_nln!(nls::MathOptNLSModel, x::AbstractVector, c::AbstractVector) increment!(nls, :neval_cons_nln) - MOI.eval_constraint(nls.ceval, c, x) + for i = 1:(nls.quadcon.nquad) + qcon = nls.quadcon.constraints[i] + c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x) + end + if nls.meta.nnln > nls.quadcon.nquad + MOI.eval_constraint(nls.ceval, view(c, (nls.quadcon.nquad + 1):(nls.meta.nnln)), x) + end return c end @@ -272,8 +280,20 @@ function NLPModels.jac_nln_structure!( rows::AbstractVector{<:Integer}, cols::AbstractVector{<:Integer}, ) - view(rows, 1:(nls.nlcon.nnzj)) .= nls.nlcon.jac_rows - view(cols, 1:(nls.nlcon.nnzj)) .= nls.nlcon.jac_cols + if nls.quadcon.nquad > 0 + index = 0 + for i = 1:(nls.quadcon.nquad) + qcon = nls.quadcon.constraints[i] + view(rows, index+1:index+qcon.nnzg) .= i + view(cols, index+1:index+qcon.nnzg) .= qcon.g + index += qcon.nnzg + end + end + if nls.meta.nnln > nls.quadcon.nquad + ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj) + view(rows, ind_nnln) .= nls.nlcon.jac_rows + view(cols, ind_nnln) .= nls.nlcon.jac_cols + end return rows, cols end @@ -285,7 +305,33 @@ end function NLPModels.jac_nln_coord!(nls::MathOptNLSModel, x::AbstractVector, vals::AbstractVector) increment!(nls, :neval_jac_nln) - MOI.eval_constraint_jacobian(nls.ceval, view(vals, 1:(nls.nlcon.nnzj)), x) + if nls.quadcon.nquad > 0 + index = 0 + view(vals, 1:nls.quadcon.nnzj) .= 0.0 + for i = 1:(nls.quadcon.nquad) + qcon = nls.quadcon.constraints[i] + for (j, ind) in enumerate(qcon.b.nzind) + k = qcon.dg[ind] + vals[index+k] += qcon.b.nzval[j] + end + for j = 1:qcon.nnzh + row = qcon.A.rows[j] + col = qcon.A.cols[j] + val = qcon.A.vals[j] + k1 = qcon.dg[row] + vals[index+k1] += val * x[col] + if row != col + k2 = qcon.dg[col] + vals[index+k2] += val * x[row] + end + end + index += qcon.nnzg + end + end + if nls.meta.nnln > nls.quadcon.nquad + ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj) + MOI.eval_constraint_jacobian(nls.ceval, view(vals, ind_nnln), x) + end return vals end @@ -313,7 +359,18 @@ function NLPModels.jprod_nln!( Jv::AbstractVector, ) increment!(nls, :neval_jprod_nln) - MOI.eval_constraint_jacobian_product(nls.ceval, Jv, x, v) + if nls.meta.nnln > nls.quadcon.nquad + ind_nnln = (nls.quadcon.nquad + 1):(nls.meta.nnln) + MOI.eval_constraint_jacobian_product(nls.ceval, view(Jv, ind_nnln), x, v) + end + (nls.meta.nnln == nls.quadcon.nquad) && (Jv .= 0.0) + if nls.quadcon.nquad > 0 + for i = 1:(nls.quadcon.nquad) + # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v + qcon = nls.quadcon.constraints[i] + v[i] += coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) + end + end return Jv end @@ -341,7 +398,19 @@ function NLPModels.jtprod_nln!( Jtv::AbstractVector, ) increment!(nls, :neval_jtprod_nln) - MOI.eval_constraint_jacobian_transpose_product(nls.ceval, Jtv, x, v) + if nls.meta.nnln > nls.quadcon.nquad + ind_nnln = (nls.quadcon.nquad + 1):(nls.meta.nnln) + MOI.eval_constraint_jacobian_transpose_product(nls.ceval, Jtv, x, view(v, ind_nnln)) + end + (nls.meta.nnln == nls.quadcon.nquad) && (Jtv .= 0.0) + if nls.quadcon.nquad > 0 + for i = 1:(nls.quadcon.nquad) + # Jtv += v[i] * (Aᵢ * x + bᵢ) + qcon = nls.quadcon.constraints[i] + coo_sym_add_mul!(rows, cols, vals, x, Jtv, v[i]) + Jtv .+= v[i] .* qcon.b + end + end return Jtv end @@ -354,9 +423,18 @@ function NLPModels.hess_structure!( view(rows, 1:(nls.lls.nnzh)) .= nls.lls.hessian.rows view(cols, 1:(nls.lls.nnzh)) .= nls.lls.hessian.cols end - if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > 0) - view(rows, (nls.lls.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_rows - view(cols, (nls.lls.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_cols + if (nls.nls_meta.nnln > 0) || (nlp.meta.nnln > nlp.quadcon.nquad) + view(rows, (nls.lls.nnzh + nlp.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_rows + view(cols, (nls.lls.nnzh + nlp.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_cols + end + if nls.quadcon.nquad > 0 + index = nls.lls.nnzh + for i = 1:(nls.quadcon.nquad) + qcon = nls.quadcon.constraints[i] + view(rows, (index + 1):(index + qcon.nnzh)) .= qcon.A.rows + view(cols, (index + 1):(index + qcon.nnzh)) .= qcon.A.cols + index += qcon.nnzh + end end return rows, cols end @@ -372,15 +450,24 @@ function NLPModels.hess_coord!( if nls.nls_meta.nlin > 0 view(vals, 1:(nls.lls.nnzh)) .= obj_weight .* nls.lls.hessian.vals end - if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > 0) + if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad) + λ = view(y, (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon)) MOI.eval_hessian_lagrangian( nls.ceval, - view(vals, (nls.lls.nnzh + 1):(nls.meta.nnzh)), + view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)), x, obj_weight, - view(y, nls.meta.nln), + λ ) end + if nls.quadcon.nquad > 0 + index = nls.lls.nnzh + for i = 1:(nls.quadcon.nquad) + qcon = nls.quadcon.constraints[i] + view(vals, (index + 1):(index + qcon.nnzh)) .= y[i] .* qcon.A.vals + index += qcon.nnzh + end + end return vals end @@ -394,16 +481,18 @@ function NLPModels.hess_coord!( if nls.nls_meta.nlin > 0 view(vals, 1:(nls.lls.nnzh)) .= obj_weight .* nls.lls.hessian.vals end + view(vals, (nls.lls.nnzh + 1):(nls.lls.nnzh + nls.quadcon.nnzh)) .= 0.0 if nls.nls_meta.nnln > 0 + λ = zeros(nlp.meta.nnln - nlp.quadcon.nquad) # Should be stored in the structure MathOptNLSModel MOI.eval_hessian_lagrangian( nls.ceval, - view(vals, (nls.lls.nnzh + 1):(nls.meta.nnzh)), + view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)), x, obj_weight, - zeros(nls.meta.nnln), + λ, ) else - view(vals, (nls.lls.nnzh + 1):(nls.meta.nnzh)) .= 0.0 + view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= 0.0 end return vals end @@ -417,11 +506,12 @@ function NLPModels.hprod!( obj_weight::Float64 = 1.0, ) increment!(nls, :neval_hprod) - if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > 0) - MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, view(y, nls.meta.nln)) + if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad) + λ = view(y, (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon)) + MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, λ) end + (nls.nls_meta.nnln == 0) && (nls.meta.nnln == nls.quadcon.nquad) && (hv .= 0.0) if nls.nls_meta.nlin > 0 - (nls.nls_meta.nnln == 0) && (nls.meta.nnln == 0) && (hv .= 0.0) coo_sym_add_mul!( nls.lls.hessian.rows, nls.lls.hessian.cols, @@ -431,6 +521,12 @@ function NLPModels.hprod!( obj_weight, ) end + if nls.quadcon.nquad > 0 + for i = 1:(nls.quadcon.nquad) + qcon = nls.quadcon.constraints[i] + coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, y[nls.meta.nlin + i]) + end + end return hv end @@ -443,7 +539,8 @@ function NLPModels.hprod!( ) increment!(nls, :neval_hprod) if nls.nls_meta.nnln > 0 - MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, zeros(nls.meta.nnln)) + λ = zeros(nls.meta.nnln - nls.quadcon.nquad) # Should be stored in the structure MathOptNLSModel + MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, λ) end if nls.nls_meta.nlin > 0 (nls.nls_meta.nnln == 0) && (hv .= 0.0) From 0810e40ff5acfeb5099703713a456b988cce374d Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 15:51:46 -0400 Subject: [PATCH 10/21] Replace nlp by nls in MathOptNLSModel --- src/moi_nls_model.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl index 6b5c991..5c57fbb 100644 --- a/src/moi_nls_model.jl +++ b/src/moi_nls_model.jl @@ -423,9 +423,9 @@ function NLPModels.hess_structure!( view(rows, 1:(nls.lls.nnzh)) .= nls.lls.hessian.rows view(cols, 1:(nls.lls.nnzh)) .= nls.lls.hessian.cols end - if (nls.nls_meta.nnln > 0) || (nlp.meta.nnln > nlp.quadcon.nquad) - view(rows, (nls.lls.nnzh + nlp.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_rows - view(cols, (nls.lls.nnzh + nlp.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_cols + if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad) + view(rows, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_rows + view(cols, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_cols end if nls.quadcon.nquad > 0 index = nls.lls.nnzh @@ -483,7 +483,7 @@ function NLPModels.hess_coord!( end view(vals, (nls.lls.nnzh + 1):(nls.lls.nnzh + nls.quadcon.nnzh)) .= 0.0 if nls.nls_meta.nnln > 0 - λ = zeros(nlp.meta.nnln - nlp.quadcon.nquad) # Should be stored in the structure MathOptNLSModel + λ = zeros(nls.meta.nnln - nls.quadcon.nquad) # Should be stored in the structure MathOptNLSModel MOI.eval_hessian_lagrangian( nls.ceval, view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)), From 547401aa95bf625992ab831bca933f85aa811242 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 16:12:26 -0400 Subject: [PATCH 11/21] Return an error if we encounter an unsupported feature --- src/utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 519ae7d..11c6ce3 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -367,8 +367,8 @@ function parser_MOI(moimodel, index_map, nvar) contypes = MOI.get(moimodel, MOI.ListOfConstraintTypesPresent()) for (F, S) in contypes - F <: AF || F <: QF || F == MOI.ScalarNonlinearFunction || F == VI || @warn("Function $F is not supported.") - S <: LS || @warn("Set $S is not supported.") + F <: AF || F <: QF || F == MOI.ScalarNonlinearFunction || F == VI || error("Function $F is not supported.") + S <: LS || error("Set $S is not supported.") conindices = MOI.get(moimodel, MOI.ListOfConstraintIndices{F, S}()) for cidx in conindices From 9952897b707fdc1b2c10eeab8f768a612c1bfdcc Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 16:22:03 -0400 Subject: [PATCH 12/21] Add Lagrange multipliers for hess_coord! and hprod! without y in MathOptNLPModel and MathOptNLSModel --- src/moi_nlp_model.jl | 10 +++++----- src/moi_nls_model.jl | 9 +++++---- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 5a7bd9a..7fda9c8 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -6,6 +6,7 @@ mutable struct MathOptNLPModel <: AbstractNLPModel{Float64, Vector{Float64}} lincon::LinearConstraints quadcon::QuadraticConstraints nlcon::NonLinearStructure + λ::Vector{Float64} obj::Objective counters::Counters end @@ -32,6 +33,7 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = nlp_data = _nlp_block(moimodel) nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) + λ = zeros(nnln - quadcon.nquad) # Lagrange multipliers for hess_coord! and hprod! without y if nlp_data.has_objective obj = Objective("NONLINEAR", 0.0, spzeros(Float64, nvar), COO(), 0) @@ -64,7 +66,7 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = name = name, ) - return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, obj, Counters()), index_map + return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, λ, obj, Counters()), index_map end function NLPModels.obj(nlp::MathOptNLPModel, x::AbstractVector) @@ -337,8 +339,7 @@ function NLPModels.hess_coord!( view(vals, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= 0.0 end if nlp.obj.type == "NONLINEAR" - λ = zeros(nlp.meta.nnln - nlp.quadcon.nquad) # Should be stored in the structure MathOptNLPModel - MOI.eval_hessian_lagrangian(nlp.eval, vals, x, obj_weight, λ) + MOI.eval_hessian_lagrangian(nlp.eval, vals, x, obj_weight, nlp.λ) end return vals end @@ -389,8 +390,7 @@ function NLPModels.hprod!( coo_sym_add_mul!(nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals, v, hv, obj_weight) end if nlp.obj.type == "NONLINEAR" - λ = zeros(nlp.meta.nnln - nlp.quadcon.nquad) # Should be stored in the structure MathOptNLPModel - MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, λ) + MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, nlp.λ) end return hv end diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl index 5c57fbb..6012f30 100644 --- a/src/moi_nls_model.jl +++ b/src/moi_nls_model.jl @@ -11,6 +11,7 @@ mutable struct MathOptNLSModel <: AbstractNLSModel{Float64, Vector{Float64}} lincon::LinearConstraints quadcon::QuadraticConstraints nlcon::NonLinearStructure + λ::Vector{Float64} counters::NLSCounters end @@ -35,6 +36,7 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri nlp_data = _nlp_block(moimodel) nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) + λ = zeros(nnln - quadcon.nquad) # Lagrange multipliers for hess_coord! and hprod! without y nequ = nlinequ + nnlnequ Fnnzj = linequ.nnzj + nlequ.nnzj @@ -76,6 +78,7 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri lincon, quadcon, nlcon, + λ, NLSCounters(), ) end @@ -483,13 +486,12 @@ function NLPModels.hess_coord!( end view(vals, (nls.lls.nnzh + 1):(nls.lls.nnzh + nls.quadcon.nnzh)) .= 0.0 if nls.nls_meta.nnln > 0 - λ = zeros(nls.meta.nnln - nls.quadcon.nquad) # Should be stored in the structure MathOptNLSModel MOI.eval_hessian_lagrangian( nls.ceval, view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)), x, obj_weight, - λ, + nls.λ, ) else view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= 0.0 @@ -539,8 +541,7 @@ function NLPModels.hprod!( ) increment!(nls, :neval_hprod) if nls.nls_meta.nnln > 0 - λ = zeros(nls.meta.nnln - nls.quadcon.nquad) # Should be stored in the structure MathOptNLSModel - MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, λ) + MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, nls.λ) end if nls.nls_meta.nlin > 0 (nls.nls_meta.nnln == 0) && (hv .= 0.0) From 0b227b47a2ad9c99b4c7c4f9efa569f421295532 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 16:28:38 -0400 Subject: [PATCH 13/21] Fix the size of lambda in MathOptNLPModel and MathOptNLSModel --- src/moi_nlp_model.jl | 2 +- src/moi_nls_model.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 7fda9c8..674e15c 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -33,7 +33,7 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = nlp_data = _nlp_block(moimodel) nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) - λ = zeros(nnln - quadcon.nquad) # Lagrange multipliers for hess_coord! and hprod! without y + λ = zeros(nnln) # Lagrange multipliers for hess_coord! and hprod! without y if nlp_data.has_objective obj = Objective("NONLINEAR", 0.0, spzeros(Float64, nvar), COO(), 0) diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl index 6012f30..7d6a621 100644 --- a/src/moi_nls_model.jl +++ b/src/moi_nls_model.jl @@ -36,7 +36,7 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri nlp_data = _nlp_block(moimodel) nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) - λ = zeros(nnln - quadcon.nquad) # Lagrange multipliers for hess_coord! and hprod! without y + λ = zeros(nnln) # Lagrange multipliers for hess_coord! and hprod! without y nequ = nlinequ + nnlnequ Fnnzj = linequ.nnzj + nlequ.nnzj From 9b33bee67f8a78b94642a320738083598f6defe1 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 23:07:04 -0400 Subject: [PATCH 14/21] Fix many bugs... --- src/moi_nlp_model.jl | 12 +++++++----- src/moi_nls_model.jl | 21 +++++++++------------ 2 files changed, 16 insertions(+), 17 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 674e15c..4d4b623 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -145,7 +145,7 @@ function NLPModels.jac_nln_structure!( end if nlp.meta.nnln > nlp.quadcon.nquad ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) - view(rows, ind_nnln) .= nlp.nlcon.jac_rows + view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nlp.nlcon.jac_rows view(cols, ind_nnln) .= nlp.nlcon.jac_cols end return rows, cols @@ -222,7 +222,7 @@ function NLPModels.jprod_nln!( for i = 1:(nlp.quadcon.nquad) # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v qcon = nlp.quadcon.constraints[i] - v[i] += coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) + Jv[i] += coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) end end return Jv @@ -261,7 +261,7 @@ function NLPModels.jtprod_nln!( for i = 1:(nlp.quadcon.nquad) # Jtv += v[i] * (Aᵢ * x + bᵢ) qcon = nlp.quadcon.constraints[i] - coo_sym_add_mul!(rows, cols, vals, x, Jtv, v[i]) + coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, Jtv, v[i]) Jtv .+= v[i] .* qcon.b end end @@ -317,7 +317,7 @@ function NLPModels.hess_coord!( index = nlp.obj.nnzh for i = 1:(nlp.quadcon.nquad) qcon = nlp.quadcon.constraints[i] - view(vals, (index + 1):(index + qcon.nnzh)) .= y[i] .* qcon.A.vals + view(vals, (index + 1):(index + qcon.nnzh)) .= y[nlp.meta.nlin + i] .* qcon.A.vals index += qcon.nnzh end end @@ -339,7 +339,9 @@ function NLPModels.hess_coord!( view(vals, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= 0.0 end if nlp.obj.type == "NONLINEAR" - MOI.eval_hessian_lagrangian(nlp.eval, vals, x, obj_weight, nlp.λ) + view(vals, 1:(nlp.obj.nnzh + nlp.quadcon.nnzh)) .= 0.0 + ind_nnln = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh) + MOI.eval_hessian_lagrangian(nlp.eval, view(vals, ind_nnln), x, obj_weight, nlp.λ) end return vals end diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl index 7d6a621..779a97d 100644 --- a/src/moi_nls_model.jl +++ b/src/moi_nls_model.jl @@ -294,7 +294,7 @@ function NLPModels.jac_nln_structure!( end if nls.meta.nnln > nls.quadcon.nquad ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj) - view(rows, ind_nnln) .= nls.nlcon.jac_rows + view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nls.nlcon.jac_rows view(cols, ind_nnln) .= nls.nlcon.jac_cols end return rows, cols @@ -371,7 +371,7 @@ function NLPModels.jprod_nln!( for i = 1:(nls.quadcon.nquad) # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v qcon = nls.quadcon.constraints[i] - v[i] += coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) + Jv[i] += coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) end end return Jv @@ -410,7 +410,7 @@ function NLPModels.jtprod_nln!( for i = 1:(nls.quadcon.nquad) # Jtv += v[i] * (Aᵢ * x + bᵢ) qcon = nls.quadcon.constraints[i] - coo_sym_add_mul!(rows, cols, vals, x, Jtv, v[i]) + coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, Jtv, v[i]) Jtv .+= v[i] .* qcon.b end end @@ -467,7 +467,7 @@ function NLPModels.hess_coord!( index = nls.lls.nnzh for i = 1:(nls.quadcon.nquad) qcon = nls.quadcon.constraints[i] - view(vals, (index + 1):(index + qcon.nnzh)) .= y[i] .* qcon.A.vals + view(vals, (index + 1):(index + qcon.nnzh)) .= y[nls.meta.nlin + i] .* qcon.A.vals index += qcon.nnzh end end @@ -486,15 +486,12 @@ function NLPModels.hess_coord!( end view(vals, (nls.lls.nnzh + 1):(nls.lls.nnzh + nls.quadcon.nnzh)) .= 0.0 if nls.nls_meta.nnln > 0 - MOI.eval_hessian_lagrangian( - nls.ceval, - view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)), - x, - obj_weight, - nls.λ, - ) + ind_nnln = (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh) + MOI.eval_hessian_lagrangian(nls.ceval, view(vals, ind_nnln), x, obj_weight, nls.λ) else - view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= 0.0 + if nls.meta.nnln > nls.quadcon.nquad + view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= 0.0 + end end return vals end From 8d10910291428bd097f4118f9b2fc6fb9c0a4504 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 20 Jun 2024 23:26:36 -0400 Subject: [PATCH 15/21] nlp -> nls --- src/moi_nls_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl index 779a97d..c2f7650 100644 --- a/src/moi_nls_model.jl +++ b/src/moi_nls_model.jl @@ -294,7 +294,7 @@ function NLPModels.jac_nln_structure!( end if nls.meta.nnln > nls.quadcon.nquad ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj) - view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nls.nlcon.jac_rows + view(rows, ind_nnln) .= nls.quadcon.nquad .+ nls.nlcon.jac_rows view(cols, ind_nnln) .= nls.nlcon.jac_cols end return rows, cols From c8713cb98c48210deffb9f96dc2ab6f5165ce813 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Fri, 21 Jun 2024 12:39:43 -0400 Subject: [PATCH 16/21] Update problems to test quadratic constraints --- test/nlp_problems/hs10.jl | 2 +- test/nlp_problems/hs11.jl | 2 +- test/nlp_problems/hs14.jl | 2 +- test/nlp_problems/hs6.jl | 2 +- test/nls_problems/nlshs20.jl | 6 +++--- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/test/nlp_problems/hs10.jl b/test/nlp_problems/hs10.jl index e6d47bd..3e657e1 100644 --- a/test/nlp_problems/hs10.jl +++ b/test/nlp_problems/hs10.jl @@ -8,7 +8,7 @@ function hs10() @objective(nlp, Min, x[1] - x[2]) - @NLconstraint(nlp, -3 * x[1]^2 + 2 * x[1] * x[2] - x[2]^2 ≥ -1) + @constraint(nlp, -3 * x[1]^2 + 2 * x[1] * x[2] - x[2]^2 ≥ -1) return nlp end diff --git a/test/nlp_problems/hs11.jl b/test/nlp_problems/hs11.jl index 2217383..fd23caf 100644 --- a/test/nlp_problems/hs11.jl +++ b/test/nlp_problems/hs11.jl @@ -8,7 +8,7 @@ function hs11() @objective(nlp, Min, (x[1] - 5)^2 + x[2]^2 - 25) - @NLconstraint(nlp, -x[1]^2 + x[2] ≥ 0) + @constraint(nlp, -x[1]^2 + x[2] ≥ 0) return nlp end diff --git a/test/nlp_problems/hs14.jl b/test/nlp_problems/hs14.jl index 01bae89..e349d98 100644 --- a/test/nlp_problems/hs14.jl +++ b/test/nlp_problems/hs14.jl @@ -8,7 +8,7 @@ function hs14() @objective(nlp, Min, (x[1] - 2)^2 + (x[2] - 1)^2) - @NLconstraint(nlp, -x[1]^2 / 4 - x[2]^2 + 1 ≥ 0) + @constraint(nlp, -x[1]^2 / 4 - x[2]^2 + 1 ≥ 0) @constraint(nlp, x[1] - 2 * x[2] + 1 == 0) diff --git a/test/nlp_problems/hs6.jl b/test/nlp_problems/hs6.jl index dde92d9..d229738 100644 --- a/test/nlp_problems/hs6.jl +++ b/test/nlp_problems/hs6.jl @@ -8,7 +8,7 @@ function hs6() @NLobjective(nlp, Min, (1 - x[1])^2) - @NLconstraint(nlp, 10 * (x[2] - x[1]^2) == 0) + @constraint(nlp, 10 * (x[2] - x[1]^2) == 0) return nlp end diff --git a/test/nls_problems/nlshs20.jl b/test/nls_problems/nlshs20.jl index 98d500b..7b1ef17 100644 --- a/test/nls_problems/nlshs20.jl +++ b/test/nls_problems/nlshs20.jl @@ -6,9 +6,9 @@ function nlshs20() x0 = [-2.0; 1.0] @variable(model, lvar[i] ≤ x[i = 1:2] ≤ uvar[i], start = x0[i]) - @NLconstraint(model, x[1] + x[2]^2 ≥ 0) - @NLconstraint(model, x[1]^2 + x[2] ≥ 0) - @NLconstraint(model, x[1]^2 + x[2]^2 ≥ 1) + @constraint(model, x[1] + x[2]^2 ≥ 0) + @constraint(model, x[1]^2 + x[2] ≥ 0) + @constraint(model, x[1]^2 + x[2]^2 ≥ 1) @expression(model, F1, 1 - x[1]) @NLexpression(model, F2, 10 * (x[2] - x[1]^2)) From 3eb3a5fc3f29ee40f882ab64d23a07bd4e0a59f2 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Fri, 21 Jun 2024 12:44:41 -0400 Subject: [PATCH 17/21] Update a constraint in hs100.jl --- test/nlp_problems/hs100.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/nlp_problems/hs100.jl b/test/nlp_problems/hs100.jl index b4af448..92156cb 100644 --- a/test/nlp_problems/hs100.jl +++ b/test/nlp_problems/hs100.jl @@ -5,7 +5,7 @@ function hs100(args...; kwargs...) @NLconstraint(nlp, 127 - 2 * x[1]^2 - 3 * x[2]^4 - x[3] - 4 * x[4]^2 - 5 * x[5] ≥ 0) @constraint(nlp, 282 - 7 * x[1] - 3 * x[2] - 10 * x[3]^2 - x[4] + x[5] ≥ 0) - @constraint(nlp, 196 - 23 * x[1] - x[2]^2 - 6 * x[6]^2 + 8 * x[7] ≥ 0) + @constraint(nlp, -196 + 23 * x[1] + x[2]^2 + 6 * x[6]^2 - 8 * x[7] ≤ 0) @constraint(nlp, -4 * x[1]^2 - x[2]^2 + 3 * x[1] * x[2] - 2 * x[3]^2 - 5 * x[6] + 11 * x[7] ≥ 0) @NLobjective( From 6b4c1fb5dc072d5bfe3945424c20f7e9576c93d5 Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Fri, 21 Jun 2024 13:27:37 -0400 Subject: [PATCH 18/21] Update nlp_consistency.jl --- test/nlp_consistency.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/nlp_consistency.jl b/test/nlp_consistency.jl index f8108d9..b116823 100644 --- a/test/nlp_consistency.jl +++ b/test/nlp_consistency.jl @@ -4,7 +4,7 @@ for problem in nlp_problems problem_f = eval(Symbol(lowercase(problem))) nlp_moi = MathOptNLPModel(problem_f()) nlps = [nlp_manual; nlp_moi] - consistent_nlps(nlps, linear_api = true) + consistent_nlps(nlps, linear_api=true, test_slack=false) view_subarray_nlp(nlp_moi) end end From 02317b4d6fe17b6f149ecebf2b21a42a3811148f Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Fri, 21 Jun 2024 13:28:15 -0400 Subject: [PATCH 19/21] Update nls_consistency.jl --- test/nls_consistency.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/nls_consistency.jl b/test/nls_consistency.jl index a7a868f..2025c52 100644 --- a/test/nls_consistency.jl +++ b/test/nls_consistency.jl @@ -8,7 +8,7 @@ for problem in nls_problems if isdefined(Main, Symbol(spc)) push!(nlss, eval(Meta.parse(spc))()) end - consistent_nlss(nlss, linear_api = true) + consistent_nlss(nlss, linear_api=true, test_slack=false) view_subarray_nls(nls_moi) end end From 7efb9e28e9a2304fdcf513a0a2843d637086bf3c Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Fri, 21 Jun 2024 15:12:19 -0400 Subject: [PATCH 20/21] Fix errors in parser_VQF --- src/utils.jl | 48 ++++++++++++++++++------------------ test/nlp_problems/quadcon.jl | 13 ++++++++++ test/nls_problems/nlsqc.jl | 13 ++++++++++ test/runtests.jl | 6 +++-- test/test_moi_nlp_model.jl | 2 +- 5 files changed, 55 insertions(+), 27 deletions(-) create mode 100644 test/nlp_problems/quadcon.jl create mode 100644 test/nls_problems/nlsqc.jl diff --git a/src/utils.jl b/src/utils.jl index 11c6ce3..a960ced 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -291,53 +291,53 @@ function parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) _index(v::MOI.VariableIndex) = index_map[v].value ncon = length(fun.constants) - for i = 1:ncon + for k = 1:ncon b = spzeros(Float64, nvar) rows = Int[] cols = Int[] vals = Float64[] - quadratic_terms = fun.quadratic_terms[i] - affine_terms = fun.affine_terms[i] - constant = fun.constants[i] - # Parse a VectorAffineTerm{Float64}(output_index, scalar_term) - for term in affine_terms - @assert term.output_index == i - b[_index(term.scalar_term.variable)] = term.scalar_term.coefficient + for affine_term in fun.affine_terms + if affine_term.output_index == k + b[_index(affine_term.scalar_term.variable)] = affine_term.scalar_term.coefficient + end end # Parse a VectorQuadraticTerm{Float64}(output_index, scalar_term) - for term in quadratic_terms - @assert term.output_index == i - i = _index(term.scalar_term.variable_1) - j = _index(term.scalar_term.variable_2) - if i ≥ j - push!(rows, i) - push!(cols, j) - else - push!(rows, j) - push!(cols, i) + for quadratic_term in fun.quadratic_terms + if quadratic_term.output_index == k + i = _index(quadratic_term.scalar_term.variable_1) + j = _index(quadratic_term.scalar_term.variable_2) + if i ≥ j + push!(rows, i) + push!(cols, j) + else + push!(rows, j) + push!(cols, i) + end + push!(vals, quadratic_term.scalar_term.coefficient) end - push!(vals, term.scalar_term.coefficient) end + constant = fun.constants[k] + if typeof(set) in (MOI.Nonnegatives, MOI.Zeros) - append!(lin_lcon, constant) + append!(quad_lcon, constant) else - append!(lin_lcon, -Inf) + append!(quad_lcon, -Inf) end if typeof(set) in (MOI.Nonpositives, MOI.Zeros) - append!(lin_ucon, -constant) + append!(quad_ucon, -constant) else - append!(lin_ucon, Inf) + append!(quad_ucon, Inf) end A = COO(rows, cols, vals) g = unique(vcat(rows, b.nzind)) # sparsity pattern of Ax + b nnzg = length(g) - dg = Dict{Int,Int}(g[i] => i for i = 1:nnzg) + dg = Dict{Int,Int}(g[p] => p for p=1:nnzg) nnzh = length(vals) qcon = QuadraticConstraint(A, b, g, dg, nnzg, nnzh) push!(qcons, qcon) diff --git a/test/nlp_problems/quadcon.jl b/test/nlp_problems/quadcon.jl new file mode 100644 index 0000000..020ca1d --- /dev/null +++ b/test/nlp_problems/quadcon.jl @@ -0,0 +1,13 @@ +function quadcon() + nlp = Model() + + @variable(nlp, x[1:3]) + + @constraint(nlp, [x[1], x[1]^2 + x[1] * x[2] + 3.0] in MOI.Nonnegatives(2)) + @constraint(nlp, [x[2], x[2]^2 + x[2] * x[3] + 4.0] in MOI.Nonpositives(2)) + @constraint(nlp, [x[3], x[3]^2 + x[3] * x[1] + 5.0] in MOI.Zeros(2)) + + @objective(nlp, Min, x[1] * x[2] * exp(x[3])) + + return nlp +end diff --git a/test/nls_problems/nlsqc.jl b/test/nls_problems/nlsqc.jl new file mode 100644 index 0000000..3c80c81 --- /dev/null +++ b/test/nls_problems/nlsqc.jl @@ -0,0 +1,13 @@ +function nlsqc() + nls = Model() + + @variable(nls, x[1:3]) + + @constraint(nls, [x[1], x[1]^2 + x[1] * x[2] + 3.0] in MOI.Nonnegatives(2)) + @constraint(nls, [x[2], x[2]^2 + x[2] * x[3] + 4.0] in MOI.Nonpositives(2)) + @constraint(nls, [x[3], x[3]^2 + x[3] * x[1] + 5.0] in MOI.Zeros(2)) + + @NLexpression(nls, F[i = 1:3], x[i]^2 - i^2) + + return MathOptNLSModel(nls, F, name = "nlsqc") +end diff --git a/test/runtests.jl b/test/runtests.jl index eeb3fee..382de7e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,9 +4,11 @@ using Test, Printf nlp_problems = setdiff(NLPModelsTest.nlp_problems, ["MGH01Feas"]) nls_problems = NLPModelsTest.nls_problems -extra_nls_problems = ["HS30", "HS43", "MGH07", "nlsnohesspb"] -for problem in lowercase.(nlp_problems ∪ ["nohesspb", "hs61", "hs100", "hs219"]) +extra_nlp_problems = ["nohesspb", "hs61", "hs100", "hs219", "quadcon"] +extra_nls_problems = ["nlsnohesspb", "HS30", "HS43", "MGH07", "nlsqc"] + +for problem in lowercase.(nlp_problems ∪ extra_nlp_problems) include(joinpath("nlp_problems", "$problem.jl")) end diff --git a/test/test_moi_nlp_model.jl b/test/test_moi_nlp_model.jl index e7a63cb..bfafe05 100644 --- a/test/test_moi_nlp_model.jl +++ b/test/test_moi_nlp_model.jl @@ -11,7 +11,7 @@ println("Testing MathOptNLPModel") "‖c(x₀)‖" ) # Test that every problem can be instantiated. -for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb", "hs61", "hs100", "hs219"])) +for prob in Symbol.(lowercase.(nlp_problems ∪ extra_nlp_problems)) prob_fn = eval(prob) nlp = MathOptNLPModel(prob_fn(), hessian = (prob != :nohesspb)) n = nlp.meta.nvar From f39f1e47af0ede8d6a563f799d0d9d6937c00161 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Fri, 21 Jun 2024 16:14:11 -0400 Subject: [PATCH 21/21] Improve utils.jl --- src/utils.jl | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index a960ced..35d3b33 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -273,9 +273,13 @@ function parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) end A = COO(rows, cols, vals) - g = unique(vcat(rows, b.nzind)) # sparsity pattern of Ax + b + g = unique(vcat(rows, cols, b.nzind)) # sparsity pattern of Ax + b nnzg = length(g) - dg = Dict{Int,Int}(g[i] => i for i = 1:nnzg) + # dg is a dictionary where: + # - The key `r` specifies a row index in the vector Ax + b. + # - The value `dg[r]` is a position in the vector (of length nnzg) + # where the non-zero entries of the Jacobian for row `r` are stored. + dg = Dict{Int,Int}(g[p] => p for p = 1:nnzg) nnzh = length(vals) qcon = QuadraticConstraint(A, b, g, dg, nnzg, nnzh) push!(qcons, qcon) @@ -335,9 +339,13 @@ function parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map) end A = COO(rows, cols, vals) - g = unique(vcat(rows, b.nzind)) # sparsity pattern of Ax + b + g = unique(vcat(rows, cols, b.nzind)) # sparsity pattern of Ax + b nnzg = length(g) - dg = Dict{Int,Int}(g[p] => p for p=1:nnzg) + # dg is a dictionary where: + # - The key `r` specifies a row index in the vector Ax + b. + # - The value `dg[r]` is a position in the vector (of length nnzg) + # where the non-zero entries of the Jacobian for row `r` are stored. + dg = Dict{Int,Int}(g[p] => p for p = 1:nnzg) nnzh = length(vals) qcon = QuadraticConstraint(A, b, g, dg, nnzg, nnzh) push!(qcons, qcon) @@ -481,9 +489,9 @@ function _nlp_block(model::MOI.ModelLike) end """ - parser_NL(jmodel, moimodel) + parser_NL(nlp_data; hessian) -Parse nonlinear constraints of a `MOI.Nonlinear.Evaluator`. +Parse nonlinear constraints of an `nlp_data`. """ function parser_NL(nlp_data; hessian::Bool = true) nnln = length(nlp_data.constraint_bounds) @@ -540,7 +548,7 @@ function parser_variables(model::MOI.ModelLike) end """ - parser_objective_MOI(moimodel, nvar) + parser_objective_MOI(moimodel, nvar, index_map) Parse linear and quadratic objective of a `MOI.ModelLike`. """ @@ -596,7 +604,7 @@ function parser_objective_MOI(moimodel, nvar, index_map) end """ - parser_linear_expression(cmodel, nvar, F) + parser_linear_expression(cmodel, nvar, index_map, F) Parse linear expressions of type `VariableRef` and `GenericAffExpr{Float64,VariableRef}`. """ @@ -689,7 +697,7 @@ function add_constraint_model(Fmodel, Fi::AbstractArray) end """ - parser_nonlinear_expression(cmodel, nvar, F) + parser_nonlinear_expression(cmodel, nvar, F; hessian) Parse nonlinear expressions of type `NonlinearExpression`. """