From a16d9edeaff4fc2f031dbf67f040e60c24cfad6e Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 11 Dec 2025 20:23:22 -0600 Subject: [PATCH 01/19] Support MOI.VectorNonlinearOracle --- Project.toml | 4 +- src/moi_nlp_model.jl | 104 ++++++++++++++++++++++++++++++++++++++++--- src/utils.jl | 11 ++++- 3 files changed, 109 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index 446fff9..e7fb9e0 100644 --- a/Project.toml +++ b/Project.toml @@ -14,8 +14,8 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] JuMP = "1.25" LinearAlgebra = "1.10" -MathOptInterface = "1.40" -NLPModels = "0.21.5" +MathOptInterface = "1.46" +NLPModels = "0.21.6" NLPModelsTest = "0.10.3" Percival = "0.7.3" Printf = "1.10" diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index c85a7c0..3942558 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 + oracles::Vector{Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}} λ::Vector{Float64} hv::Vector{Float64} obj::Objective @@ -35,6 +36,8 @@ 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) + oracles = Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}[] + counters = Counters() λ = zeros(Float64, nnln) # Lagrange multipliers for hess_coord! and hprod! without y hv = zeros(Float64, nvar) # workspace for ghjvprod! @@ -44,6 +47,7 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = obj = parser_objective_MOI(moimodel, nvar, index_map) end + # Update ncon, lcon, ucon, nnzj and nnzh for the oracles ncon = nlin + quadcon.nquad + nnln lcon = vcat(lin_lcon, quad_lcon, nl_lcon) ucon = vcat(lin_ucon, quad_ucon, nl_ucon) @@ -67,9 +71,13 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = minimize = MOI.get(moimodel, MOI.ObjectiveSense()) == MOI.MIN_SENSE, islp = (obj.type == "LINEAR") && (nnln == 0) && (quadcon.nquad == 0), name = name, + jprod_available = isempty(oracles), + jtprod_available = isempty(oracles), + hprod_available = isempty(oracles) && hessian, + hess_available = hessian, ) - return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, λ, hv, obj, Counters()), + return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, oracles, λ, hv, obj, counters), index_map end @@ -120,7 +128,18 @@ function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::Abstrac end end if nlp.meta.nnln > nlp.quadcon.nquad - MOI.eval_constraint(nlp.eval, view(c, (nlp.quadcon.nquad + 1):(nlp.meta.nnln)), x) + offset = nlp.quadcon.nquad + for (f, s) in nlp.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + index_oracle = (offset + 1):(offset + s.set.output_dimension) + c_oracle = view(c, index_oracle) + s.set.eval_f(c_oracle, s.x) + offset += s.set.output_dimension + end + index_nnln = (offset + 1):(nlp.meta.nnln) + MOI.eval_constraint(nlp.eval, view(c, index_nnln), x) end return c end @@ -139,7 +158,17 @@ function NLPModels.cons!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVec end end if nlp.meta.nnln > nlp.quadcon.nquad - index_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon) + offset = nlp.meta.nlin + nlp.quadcon.nquad + for (f, s) in nlp.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + index_oracle = (offset + 1):(offset + s.set.output_dimension) + c_oracle = view(c, index_oracle) + s.set.eval_f(c_oracle, s.x) + offset += s.set.output_dimension + end + index_nnln = (offset + 1):(nlp.meta.ncon) MOI.eval_constraint(nlp.eval, view(c, index_nnln), x) end end @@ -174,7 +203,14 @@ function NLPModels.jac_nln_structure!( end end if nlp.meta.nnln > nlp.quadcon.nquad - ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) + offset = nlp.quadcon.nnzj + # for (f, s) in nlp.oracles + # for (i, j) in s.set.jacobian_structure + # ... + # offset += 1 + # end + # end + ind_nnln = (offset + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nlp.nlcon.jac_rows view(cols, ind_nnln) .= nlp.nlcon.jac_cols end @@ -204,7 +240,14 @@ function NLPModels.jac_structure!( end end if nlp.meta.nnln > nlp.quadcon.nquad - ind_nnln = (nlp.lincon.nnzj + nlp.quadcon.nnzj + 1):(nlp.meta.nnzj) + offset = nlp.lincon.nnzj + nlp.quadcon.nnzj + # for (f, s) in nlp.oracles + # for (i, j) in s.set.jacobian_structure + # ... + # offset += 1 + # end + # end + ind_nnln = (offset + 1):(nlp.meta.nnzj) view(rows, ind_nnln) .= nlp.meta.nlin .+ nlp.quadcon.nquad .+ nlp.nlcon.jac_rows view(cols, ind_nnln) .= nlp.nlcon.jac_cols end @@ -246,7 +289,16 @@ function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals: end end if nlp.meta.nnln > nlp.quadcon.nquad - ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) + offset = nlp.quadcon.nnzj + # for (f, s) in nlp.oracles + # for i in 1:s.set.input_dimension + # s.x[i] = x[f.variables[i].value] + # end + # nnz_oracle = length(s.set.jacobian_structure) + # s.set.eval_jacobian(view(values, (offset + 1):(oracle + nnz_oracle)), s.x) + # offset += nnz_oracle + # end + ind_nnln = (offset + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) end return vals @@ -284,7 +336,16 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs end end if nlp.meta.nnln > nlp.quadcon.nquad - ind_nnln = (nlp.lincon.nnzj + nlp.quadcon.nnzj + 1):(nlp.meta.nnzj) + offset = nlp.lincon.nnzj + nlp.quadcon.nnzj + # for (f, s) in nlp.oracles + # for i in 1:s.set.input_dimension + # s.x[i] = x[f.variables[i].value] + # end + # nnz_oracle = length(s.set.jacobian_structure) + # s.set.eval_jacobian(view(values, (offset + 1):(oracle + nnz_oracle)), s.x) + # offset += nnz_oracle + # end + ind_nnln = (offset + 1):(nlp.meta.nnzj) MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) end end @@ -297,6 +358,7 @@ function NLPModels.jprod_lin!( v::AbstractVector, Jv::AbstractVector, ) + nlp.meta.jprod_available || error("J * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jprod_lin) jprod_lin!( nlp, @@ -315,6 +377,7 @@ function NLPModels.jprod_nln!( v::AbstractVector, Jv::AbstractVector, ) + nlp.meta.jprod_available || error("J * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jprod_nln) if nlp.quadcon.nquad > 0 for i = 1:(nlp.quadcon.nquad) @@ -336,6 +399,7 @@ function NLPModels.jprod!( v::AbstractVector, Jv::AbstractVector, ) + nlp.meta.jprod_available || error("J * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jprod) if nlp.meta.nlin > 0 view(Jv, nlp.meta.lin) .= 0.0 @@ -371,6 +435,7 @@ function NLPModels.jtprod_lin!( v::AbstractVector, Jtv::AbstractVector, ) + nlp.meta.jtprod_available || error("J' * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jtprod_lin) jtprod_lin!( nlp, @@ -389,6 +454,7 @@ function NLPModels.jtprod_nln!( v::AbstractVector, Jtv::AbstractVector, ) + nlp.meta.jtprod_available || error("J' * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jtprod_nln) if nlp.meta.nnln > nlp.quadcon.nquad ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln) @@ -412,6 +478,7 @@ function NLPModels.jtprod!( v::AbstractVector, Jtv::AbstractVector, ) + nlp.meta.jtprod_available || error("J' * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jtprod) if nlp.meta.nnln > nlp.quadcon.nquad ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon) @@ -462,6 +529,13 @@ function NLPModels.hess_structure!( index += qcon.nnzh end end + # if !isempty(nlp.oracles) + # for (f, s) in nlp.oracles + # for (i, j) in s.set.hessian_lagrangian_structure + # ... + # end + # end + # end return rows, cols end @@ -494,6 +568,15 @@ function NLPModels.hess_coord!( index += qcon.nnzh end end + # if !isempty(nlp.oracles) + # for (f, s) in nlp.oracles + # for i in 1:s.set.input_dimension + # s.x[i] = x[f.variables[i].value] + # end + # H_nnz = length(s.set.hessian_lagrangian_structure) + # ... + # end + # end return vals end @@ -548,6 +631,9 @@ function NLPModels.jth_hess_coord!( ) nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 0.0 end + # if !isempty(nlp.oracles) + # ... + # end return vals end @@ -559,6 +645,7 @@ function NLPModels.hprod!( hv::AbstractVector; obj_weight::Float64 = 1.0, ) + nlp.meta.hprod_available && error("H * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_hprod) if (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == 0) hv .= 0.0 @@ -595,6 +682,7 @@ function NLPModels.hprod!( hv::AbstractVector; obj_weight::Float64 = 1.0, ) + nlp.meta.hprod_available && error("H * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_hprod) if nlp.obj.type == "LINEAR" hv .= 0.0 @@ -623,6 +711,7 @@ function NLPModels.jth_hprod!( j::Integer, hv::AbstractVector, ) + nlp.meta.hprod_available && error("H * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jhprod) @rangecheck 1 nlp.meta.ncon j hv .= 0.0 @@ -644,6 +733,7 @@ function NLPModels.ghjvprod!( v::AbstractVector, ghv::AbstractVector, ) + nlp.meta.hprod_available && error("The operation is not supported by this MathOptNLPModel.") increment!(nlp, :neval_hprod) ghv .= 0.0 for i = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad) diff --git a/src/utils.jl b/src/utils.jl index a98b54b..84fe4a7 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -24,6 +24,16 @@ const SNF = MOI.ScalarNonlinearFunction const VNF = MOI.VectorNonlinearFunction const NF = Union{SNF, VNF} +# Cache of VectorNonlinearOracle +mutable struct _VectorNonlinearOracleCache + set::MOI.VectorNonlinearOracle{Float64} + x::Vector{Float64} + + function _VectorNonlinearOracleCache(set::MOI.VectorNonlinearOracle{Float64}) + return new(set, zeros(set.input_dimension)) + end +end + # AffLinSets and VecLinSets const ALS = Union{ MOI.EqualTo{Float64}, @@ -444,7 +454,6 @@ function parser_MOI(moimodel, index_map, nvar) coo = COO(linrows, lincols, linvals) lin_nnzj = length(linvals) lincon = LinearConstraints(coo, lin_nnzj) - quad_nnzj = 0 quad_nnzh = 0 for i = 1:nquad From 1b0bd359e5ddb69ed878e5057086ed46b8523c91 Mon Sep 17 00:00:00 2001 From: Andrew Date: Fri, 12 Dec 2025 15:05:24 -0500 Subject: [PATCH 02/19] skip support check oracle parser_MOI --- src/utils.jl | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/src/utils.jl b/src/utils.jl index 84fe4a7..de7be9b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -84,6 +84,17 @@ mutable struct QuadraticConstraints nnzh::Int end +""" + NonLinearStructure + +Structure containing Jacobian and Hessian structures of nonlinear constraints: +- jac_rows: row indices of the Jacobian in Coordinate format (COO) format +- jac_cols: column indices of the Jacobian in COO format +- nnzj: number of non-zero entries in the Jacobian +- hess_rows: row indices of the Hessian in COO format +- hess_cols: column indices of the Hessian in COO format +- nnzh: number of non-zero entries in the Hessian +""" mutable struct NonLinearStructure jac_rows::Vector{Int} jac_cols::Vector{Int} @@ -417,6 +428,10 @@ function parser_MOI(moimodel, index_map, nvar) contypes = MOI.get(moimodel, MOI.ListOfConstraintTypesPresent()) for (F, S) in contypes + # Ignore VectorNonlinearOracle here, we'll parse it separately + if F == MOI.VectorOfVariables && S <: MOI.VectorNonlinearOracle{Float64} + continue + end (F == VNF) && error( "The function $F is not supported. Please use `.<=`, `.==`, and `.>=` in your constraints to ensure compatibility with ScalarNonlinearFunction.", ) @@ -538,6 +553,12 @@ end parser_NL(nlp_data; hessian) Parse nonlinear constraints of an `nlp_data`. + +Returns: +- nnln: number of nonlinear constraints +- nlcon: NonLinearStructure containing Jacobian and Hessian structures +- nl_lcon: lower bounds of nonlinear constraints +- nl_ucon: upper bounds of nonlinear constraints """ function parser_NL(nlp_data; hessian::Bool = true) nnln = length(nlp_data.constraint_bounds) From a549525f7c89eaffe698f5a01b16b1ced0957038 Mon Sep 17 00:00:00 2001 From: Andrew Date: Fri, 12 Dec 2025 15:09:52 -0500 Subject: [PATCH 03/19] add parser oracles --- src/utils.jl | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/utils.jl b/src/utils.jl index de7be9b..77d2d90 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -581,6 +581,47 @@ function parser_NL(nlp_data; hessian::Bool = true) return nnln, nlcon, nl_lcon, nl_ucon end +""" + parser_oracles(moimodel, index_map) + +Parse nonlinear oracles of a `MOI.ModelLike`. +""" +function parser_oracles(moimodel, index_map) + oracles = Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}[] + l_oracle = Float64[] + u_oracle = Float64[] + + # We know this pair exists from ListOfConstraintTypesPresent + for ci in MOI.get( + moimodel, + MOI.ListOfConstraintIndices{MOI.VectorOfVariables, MOI.VectorNonlinearOracle{Float64}}(), + ) + f = MOI.get(moimodel, MOI.ConstraintFunction(), ci) # ::MOI.VectorOfVariables + set = MOI.get(moimodel, MOI.ConstraintSet(), ci) # ::MOI.VectorNonlinearOracle{Float64} + + cache = _VectorNonlinearOracleCache(set) + push!(oracles, (f, cache)) + + # Bounds: MOI.VectorNonlinearOracle stores them internally (l, u) + append!(l_oracle, set.l) + append!(u_oracle, set.u) + end + + # Sizes: number of scalar constraints represented by all oracles + noracle = length(l_oracle) + + # Sparsity: + nnzj_oracle = 0 + nnzh_oracle = 0 + for (_, cache) in oracles + nnzj_oracle += length(cache.set.jacobian_structure) + # there may or may not be Hessian info + nnzh_oracle += length(cache.set.hessian_lagrangian_structure) + end + + return oracles, noracle, l_oracle, u_oracle, nnzj_oracle, nnzh_oracle +end + """ parser_variables(model) From 659dc56043f8ede7b4eed9bafc54fa7b54d2a6df Mon Sep 17 00:00:00 2001 From: Andrew Date: Fri, 12 Dec 2025 15:16:22 -0500 Subject: [PATCH 04/19] count number of elements considering oracles --- src/moi_nlp_model.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 3942558..762a8b0 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -36,6 +36,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) + oracles, noracle, oracle_lcon, oracle_ucon, nnzj_oracle, nnzh_oracle = parser_oracles(moimodel, index_map) oracles = Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}[] counters = Counters() λ = zeros(Float64, nnln) # Lagrange multipliers for hess_coord! and hprod! without y @@ -47,12 +48,12 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = obj = parser_objective_MOI(moimodel, nvar, index_map) end - # Update ncon, lcon, ucon, nnzj and nnzh for the oracles - 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 + # Total counts + ncon = nlin + quadcon.nquad + nnln + noracle + lcon = vcat(lin_lcon, quad_lcon, nl_lcon, oracle_lcon) + ucon = vcat(lin_ucon, quad_ucon, nl_ucon, oracle_ucon) + nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj + nnzj_oracle + nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh + nnzh_oracle meta = NLPModelMeta( nvar, @@ -67,7 +68,7 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = nnzh = nnzh, lin = collect(1:nlin), lin_nnzj = lincon.nnzj, - nln_nnzj = quadcon.nnzj + nlcon.nnzj, + nln_nnzj = quadcon.nnzj + nlcon.nnzj + nnzj_oracle, minimize = MOI.get(moimodel, MOI.ObjectiveSense()) == MOI.MIN_SENSE, islp = (obj.type == "LINEAR") && (nnln == 0) && (quadcon.nquad == 0), name = name, From a149452bf585f05709dcd29e339be6639e348592 Mon Sep 17 00:00:00 2001 From: Andrew Date: Fri, 12 Dec 2025 15:42:30 -0500 Subject: [PATCH 05/19] update NLPModels.jac_nln_structure! --- src/moi_nlp_model.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 762a8b0..1a4102e 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -205,14 +205,24 @@ function NLPModels.jac_nln_structure!( end if nlp.meta.nnln > nlp.quadcon.nquad offset = nlp.quadcon.nnzj - # for (f, s) in nlp.oracles - # for (i, j) in s.set.jacobian_structure - # ... - # offset += 1 - # end - # end - ind_nnln = (offset + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) - view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nlp.nlcon.jac_rows + # structure of oracle Jacobians + for (k, (f, s)) in enumerate(nlp.oracles) + # Shift row index by quadcon.nquad + # plus previous oracle outputs. + row_offset = nlp.quadcon.nquad + for i in 1:(k - 1) + row_offset += nlp.oracles[i][2].set.output_dimension + end + + for (r, c) in s.set.jacobian_structure + offset += 1 + rows[offset] = row_offset + r + cols[offset] = f.variables[c].value + end + end + # non-oracle nonlinear constraints + ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) + view(rows, ind_nnln) .= nlp.quadcon.nquad .+ sum(s.set.output_dimension for (_, s) in nlp.oracles) .+ nlp.nlcon.jac_rows view(cols, ind_nnln) .= nlp.nlcon.jac_cols end return rows, cols From df7269188566dac7c0b990818e3fd9ed8dfd4fcd Mon Sep 17 00:00:00 2001 From: Andrew Date: Fri, 12 Dec 2025 15:46:13 -0500 Subject: [PATCH 06/19] update NLPModels.jac_nln_coord! --- src/moi_nlp_model.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 1a4102e..1b53088 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -301,15 +301,15 @@ function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals: end if nlp.meta.nnln > nlp.quadcon.nquad offset = nlp.quadcon.nnzj - # for (f, s) in nlp.oracles - # for i in 1:s.set.input_dimension - # s.x[i] = x[f.variables[i].value] - # end - # nnz_oracle = length(s.set.jacobian_structure) - # s.set.eval_jacobian(view(values, (offset + 1):(oracle + nnz_oracle)), s.x) - # offset += nnz_oracle - # end - ind_nnln = (offset + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj) + for (f, s) in nlp.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + nnz_oracle = length(s.set.jacobian_structure) + s.set.eval_jacobian(view(values, (offset + 1):(offset + nnz_oracle)), s.x) + offset += nnz_oracle + end + ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) end return vals From 4499347154a1761c34dc6e2ad43617c99743268f Mon Sep 17 00:00:00 2001 From: Andrew Date: Fri, 12 Dec 2025 15:47:17 -0500 Subject: [PATCH 07/19] update NLPModels.jac_coord! --- src/moi_nlp_model.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 1b53088..0f58dbe 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -348,14 +348,14 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs end if nlp.meta.nnln > nlp.quadcon.nquad offset = nlp.lincon.nnzj + nlp.quadcon.nnzj - # for (f, s) in nlp.oracles - # for i in 1:s.set.input_dimension - # s.x[i] = x[f.variables[i].value] - # end - # nnz_oracle = length(s.set.jacobian_structure) - # s.set.eval_jacobian(view(values, (offset + 1):(oracle + nnz_oracle)), s.x) - # offset += nnz_oracle - # end + for (f, s) in nlp.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + nnz_oracle = length(s.set.jacobian_structure) + s.set.eval_jacobian(view(values, (offset + 1):(offset + nnz_oracle)), s.x) + offset += nnz_oracle + end ind_nnln = (offset + 1):(nlp.meta.nnzj) MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) end From df356d7f2a7cafc83426a2b5359841cc7249e657 Mon Sep 17 00:00:00 2001 From: Andrew Date: Fri, 12 Dec 2025 18:08:37 -0500 Subject: [PATCH 08/19] update NLPModels.hess_structure! --- src/moi_nlp_model.jl | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 0f58dbe..5bd6d50 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -527,12 +527,8 @@ 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 > 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 + index = nlp.obj.nnzh 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 @@ -540,13 +536,19 @@ function NLPModels.hess_structure!( index += qcon.nnzh end end - # if !isempty(nlp.oracles) - # for (f, s) in nlp.oracles - # for (i, j) in s.set.hessian_lagrangian_structure - # ... - # end - # end - # end + if !isempty(nlp.oracles) + for (f, s) in nlp.oracles + for (i_local, j_local) in s.set.hessian_lagrangian_structure + index += 1 + rows[index] = f.variables[i_local].value + cols[index] = f.variables[j_local].value + end + end + end + if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) + view(rows, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_rows + view(cols, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_cols + end return rows, cols end From ea99eb3e24ed66e45fdd2f972573fd281d3f4c0c Mon Sep 17 00:00:00 2001 From: Andrew Date: Mon, 15 Dec 2025 15:26:37 -0500 Subject: [PATCH 09/19] rearrange order nnl and oracle --- src/moi_nlp_model.jl | 95 ++++++++++++++++++++++---------------------- src/utils.jl | 33 ++++++++++++--- 2 files changed, 76 insertions(+), 52 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 5bd6d50..0146051 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -6,7 +6,7 @@ mutable struct MathOptNLPModel <: AbstractNLPModel{Float64, Vector{Float64}} lincon::LinearConstraints quadcon::QuadraticConstraints nlcon::NonLinearStructure - oracles::Vector{Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}} + oracles_data::OraclesData λ::Vector{Float64} hv::Vector{Float64} obj::Objective @@ -36,8 +36,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) - oracles, noracle, oracle_lcon, oracle_ucon, nnzj_oracle, nnzh_oracle = parser_oracles(moimodel, index_map) - oracles = Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}[] + oracles_data = parser_oracles(moimodel) counters = Counters() λ = zeros(Float64, nnln) # Lagrange multipliers for hess_coord! and hprod! without y hv = zeros(Float64, nvar) # workspace for ghjvprod! @@ -49,12 +48,11 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = end # Total counts - ncon = nlin + quadcon.nquad + nnln + noracle - lcon = vcat(lin_lcon, quad_lcon, nl_lcon, oracle_lcon) - ucon = vcat(lin_ucon, quad_ucon, nl_ucon, oracle_ucon) - nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj + nnzj_oracle - nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh + nnzh_oracle - + ncon = nlin + quadcon.nquad + nnln + oracles_data.noracle + lcon = vcat(lin_lcon, quad_lcon, nl_lcon, oracles_data.l_oracle) + ucon = vcat(lin_ucon, quad_ucon, nl_ucon, oracles_data.u_oracle) + nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj + oracles_data.nnzj_oracle + nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh + oracles_data.nnzh_oracle meta = NLPModelMeta( nvar, x0 = x0, @@ -68,17 +66,17 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = nnzh = nnzh, lin = collect(1:nlin), lin_nnzj = lincon.nnzj, - nln_nnzj = quadcon.nnzj + nlcon.nnzj + nnzj_oracle, + nln_nnzj = quadcon.nnzj + nlcon.nnzj + oracles_data.nnzj_oracle, minimize = MOI.get(moimodel, MOI.ObjectiveSense()) == MOI.MIN_SENSE, islp = (obj.type == "LINEAR") && (nnln == 0) && (quadcon.nquad == 0), name = name, - jprod_available = isempty(oracles), - jtprod_available = isempty(oracles), - hprod_available = isempty(oracles) && hessian, + jprod_available = isempty(oracles_data.oracles), + jtprod_available = isempty(oracles_data.oracles), + hprod_available = isempty(oracles_data.oracles) && hessian, hess_available = hessian, ) - return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, oracles, λ, hv, obj, counters), + return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, oracles_data, λ, hv, obj, counters), index_map end @@ -129,8 +127,9 @@ function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::Abstrac end end if nlp.meta.nnln > nlp.quadcon.nquad - offset = nlp.quadcon.nquad - for (f, s) in nlp.oracles + offset = nlp.quadcon.nquad + nlp.nlcon.nnln + MOI.eval_constraint(nlp.eval, view(c, (nlp.quadcon.nquad + 1):(offset)), x) + for (f, s) in nlp.oracles_data.oracles for i in 1:s.set.input_dimension s.x[i] = x[f.variables[i].value] end @@ -139,8 +138,6 @@ function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::Abstrac s.set.eval_f(c_oracle, s.x) offset += s.set.output_dimension end - index_nnln = (offset + 1):(nlp.meta.nnln) - MOI.eval_constraint(nlp.eval, view(c, index_nnln), x) end return c end @@ -159,8 +156,9 @@ function NLPModels.cons!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVec end end if nlp.meta.nnln > nlp.quadcon.nquad - offset = nlp.meta.nlin + nlp.quadcon.nquad - for (f, s) in nlp.oracles + offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + MOI.eval_constraint(nlp.eval, view(c, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(offset)), x) + for (f, s) in nlp.oracles_data.oracles for i in 1:s.set.input_dimension s.x[i] = x[f.variables[i].value] end @@ -169,8 +167,6 @@ function NLPModels.cons!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVec s.set.eval_f(c_oracle, s.x) offset += s.set.output_dimension end - index_nnln = (offset + 1):(nlp.meta.ncon) - MOI.eval_constraint(nlp.eval, view(c, index_nnln), x) end end return c @@ -205,13 +201,18 @@ function NLPModels.jac_nln_structure!( end if nlp.meta.nnln > nlp.quadcon.nquad offset = nlp.quadcon.nnzj + # non-oracle nonlinear constraints + ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) + view(rows, ind_nnln) .= nlp.quadcon.nquad .+ sum(s.set.output_dimension for (_, s) in nlp.oracles_data.oracles) .+ nlp.nlcon.jac_rows + view(cols, ind_nnln) .= nlp.nlcon.jac_cols + offset += nlp.nlcon.nnzj # structure of oracle Jacobians - for (k, (f, s)) in enumerate(nlp.oracles) + for (k, (f, s)) in enumerate(nlp.oracles_data.oracles) # Shift row index by quadcon.nquad # plus previous oracle outputs. row_offset = nlp.quadcon.nquad for i in 1:(k - 1) - row_offset += nlp.oracles[i][2].set.output_dimension + row_offset += nlp.oracles_data.oracles[i][2].set.output_dimension end for (r, c) in s.set.jacobian_structure @@ -220,10 +221,6 @@ function NLPModels.jac_nln_structure!( cols[offset] = f.variables[c].value end end - # non-oracle nonlinear constraints - ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) - view(rows, ind_nnln) .= nlp.quadcon.nquad .+ sum(s.set.output_dimension for (_, s) in nlp.oracles) .+ nlp.nlcon.jac_rows - view(cols, ind_nnln) .= nlp.nlcon.jac_cols end return rows, cols end @@ -252,15 +249,16 @@ function NLPModels.jac_structure!( end if nlp.meta.nnln > nlp.quadcon.nquad offset = nlp.lincon.nnzj + nlp.quadcon.nnzj - # for (f, s) in nlp.oracles - # for (i, j) in s.set.jacobian_structure - # ... - # offset += 1 - # end - # end - ind_nnln = (offset + 1):(nlp.meta.nnzj) + ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) view(rows, ind_nnln) .= nlp.meta.nlin .+ nlp.quadcon.nquad .+ nlp.nlcon.jac_rows view(cols, ind_nnln) .= nlp.nlcon.jac_cols + for (f, s) in nlp.oracles_data.oracles + for (i, j) in s.set.jacobian_structure + offset += 1 + rows[offset] = nlp.meta.nlin + nlp.quadcon.nquad + s.set.output_dimension + i + cols[offset] = f.variables[j].value + end + end end end return rows, cols @@ -301,7 +299,10 @@ function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals: end if nlp.meta.nnln > nlp.quadcon.nquad offset = nlp.quadcon.nnzj - for (f, s) in nlp.oracles + ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) + MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) + offset += nlp.nlcon.nnzj + for (f, s) in nlp.oracles_data.oracles for i in 1:s.set.input_dimension s.x[i] = x[f.variables[i].value] end @@ -309,8 +310,6 @@ function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals: s.set.eval_jacobian(view(values, (offset + 1):(offset + nnz_oracle)), s.x) offset += nnz_oracle end - ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) - MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) end return vals end @@ -348,7 +347,10 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs end if nlp.meta.nnln > nlp.quadcon.nquad offset = nlp.lincon.nnzj + nlp.quadcon.nnzj - for (f, s) in nlp.oracles + ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) + MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) + offset += nlp.nlcon.nnzj + for (f, s) in nlp.oracles_data.oracles for i in 1:s.set.input_dimension s.x[i] = x[f.variables[i].value] end @@ -356,8 +358,6 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs s.set.eval_jacobian(view(values, (offset + 1):(offset + nnz_oracle)), s.x) offset += nnz_oracle end - ind_nnln = (offset + 1):(nlp.meta.nnzj) - MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) end end return vals @@ -536,8 +536,13 @@ function NLPModels.hess_structure!( index += qcon.nnzh end end - if !isempty(nlp.oracles) - for (f, s) in nlp.oracles + if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) + view(rows, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_rows + view(cols, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_cols + index += nlp.nlcon.nnzh + end + if !isempty(nlp.oracles_data.oracles) + for (f, s) in nlp.oracles_data.oracles for (i_local, j_local) in s.set.hessian_lagrangian_structure index += 1 rows[index] = f.variables[i_local].value @@ -545,10 +550,6 @@ function NLPModels.hess_structure!( end end end - if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) - view(rows, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_rows - view(cols, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_cols - end return rows, cols end diff --git a/src/utils.jl b/src/utils.jl index 77d2d90..54f29c0 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -88,6 +88,7 @@ end NonLinearStructure Structure containing Jacobian and Hessian structures of nonlinear constraints: +- nnln: number of nonlinear constraints - jac_rows: row indices of the Jacobian in Coordinate format (COO) format - jac_cols: column indices of the Jacobian in COO format - nnzj: number of non-zero entries in the Jacobian @@ -96,6 +97,7 @@ Structure containing Jacobian and Hessian structures of nonlinear constraints: - nnzh: number of non-zero entries in the Hessian """ mutable struct NonLinearStructure + nnln::Int jac_rows::Vector{Int} jac_cols::Vector{Int} nnzj::Int @@ -118,6 +120,26 @@ mutable struct Objective nnzh::Int end +""" + OraclesData + +Structure containing nonlinear oracles data: +- oracles: vector of tuples (MOI.VectorOfVariables, _VectorNonlinearOracleCache) +- noracle: number of scalar constraints represented by all oracles +- l_oracle: lower bounds of oracle constraints +- u_oracle: upper bounds of oracle constraints +- nnzj_oracle: number of non-zero entries in the Jacobian of all oracles +- nnzh_oracle: number of non-zero entries in the Hessian of all oracles +""" +mutable struct OraclesData + oracles::Vector{Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}} + noracle::Int + l_oracle::Vector{Float64} + u_oracle::Vector{Float64} + nnzj_oracle::Int + nnzh_oracle::Int +end + """ replace!(ex, x) @@ -576,17 +598,17 @@ function parser_NL(nlp_data; hessian::Bool = true) hess_rows = hessian ? getindex.(hess, 1) : Int[] hess_cols = hessian ? getindex.(hess, 2) : Int[] nnzh = length(hess) - nlcon = NonLinearStructure(jac_rows, jac_cols, nnzj, hess_rows, hess_cols, nnzh) + nlcon = NonLinearStructure(nnln, jac_rows, jac_cols, nnzj, hess_rows, hess_cols, nnzh) return nnln, nlcon, nl_lcon, nl_ucon end """ - parser_oracles(moimodel, index_map) + parser_oracles(moimodel) Parse nonlinear oracles of a `MOI.ModelLike`. """ -function parser_oracles(moimodel, index_map) +function parser_oracles(moimodel) oracles = Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}[] l_oracle = Float64[] u_oracle = Float64[] @@ -619,7 +641,7 @@ function parser_oracles(moimodel, index_map) nnzh_oracle += length(cache.set.hessian_lagrangian_structure) end - return oracles, noracle, l_oracle, u_oracle, nnzj_oracle, nnzh_oracle + return OraclesData(oracles, noracle, l_oracle, u_oracle, nnzj_oracle, nnzh_oracle) end """ @@ -813,6 +835,7 @@ function parser_nonlinear_expression(cmodel, nvar, F; hessian::Bool = true) # Nonlinear least squares model F_is_array_of_containers = F isa Array{<:AbstractArray} + nnlnequ = 0 if F_is_array_of_containers nnlnequ = sum(sum(isa(Fi, NLE) for Fi in FF) for FF in F) if nnlnequ > 0 @@ -849,7 +872,7 @@ function parser_nonlinear_expression(cmodel, nvar, F; hessian::Bool = true) Fhess_cols = hessian && Feval ≠ nothing ? getindex.(Fhess, 2) : Int[] nl_Fnnzh = length(Fhess) - nlequ = NonLinearStructure(Fjac_rows, Fjac_cols, nl_Fnnzj, Fhess_rows, Fhess_cols, nl_Fnnzh) + nlequ = NonLinearStructure(nnlnequ, Fjac_rows, Fjac_cols, nl_Fnnzj, Fhess_rows, Fhess_cols, nl_Fnnzh) return Feval, nlequ, nnlnequ end From 532725659787b9c17227448a7bdcd07b15018085 Mon Sep 17 00:00:00 2001 From: Andrew Date: Mon, 15 Dec 2025 20:07:42 -0500 Subject: [PATCH 10/19] fix bug --- src/moi_nlp_model.jl | 108 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 93 insertions(+), 15 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 0146051..19eb9b3 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -203,7 +203,7 @@ function NLPModels.jac_nln_structure!( offset = nlp.quadcon.nnzj # non-oracle nonlinear constraints ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) - view(rows, ind_nnln) .= nlp.quadcon.nquad .+ sum(s.set.output_dimension for (_, s) in nlp.oracles_data.oracles) .+ nlp.nlcon.jac_rows + view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nlp.nlcon.jac_rows view(cols, ind_nnln) .= nlp.nlcon.jac_cols offset += nlp.nlcon.nnzj # structure of oracle Jacobians @@ -248,16 +248,21 @@ function NLPModels.jac_structure!( end end if nlp.meta.nnln > nlp.quadcon.nquad + # non-oracle nonlinear constraints offset = nlp.lincon.nnzj + nlp.quadcon.nnzj ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) view(rows, ind_nnln) .= nlp.meta.nlin .+ nlp.quadcon.nquad .+ nlp.nlcon.jac_rows view(cols, ind_nnln) .= nlp.nlcon.jac_cols + offset += nlp.nlcon.nnzj + # structure of oracle Jacobians + row_offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln for (f, s) in nlp.oracles_data.oracles for (i, j) in s.set.jacobian_structure offset += 1 - rows[offset] = nlp.meta.nlin + nlp.quadcon.nquad + s.set.output_dimension + i + rows[offset] = row_offset + i cols[offset] = f.variables[j].value end + row_offset += s.set.output_dimension end end end @@ -360,6 +365,7 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs end end end + @assert offset == nlp.meta.nnzj return vals end @@ -561,9 +567,16 @@ function NLPModels.hess_coord!( obj_weight::Float64 = 1.0, ) increment!(nlp, :neval_hess) + + # Running index over Hessian nonzeros we've filled so far. + index = 0 + + # Objective Hessian block if nlp.obj.type == "QUADRATIC" view(vals, 1:(nlp.obj.nnzh)) .= obj_weight .* nlp.obj.hessian.vals + index += nlp.obj.nnzh end + # Nonlinear objective Hessian block 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( @@ -574,23 +587,42 @@ function NLPModels.hess_coord!( λ, ) end + # Quadratic constraint Hessian blocks 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[nlp.meta.nlin + i] .* qcon.A.vals index += qcon.nnzh end end - # if !isempty(nlp.oracles) - # for (f, s) in nlp.oracles - # for i in 1:s.set.input_dimension - # s.x[i] = x[f.variables[i].value] - # end - # H_nnz = length(s.set.hessian_lagrangian_structure) - # ... - # end - # end + # Oracle Hessian blocks + if !isempty(nlp.oracles_data.oracles) + λ_oracle_all = view( + y, + (nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1): + (nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + nlp.oracles_data.noracle), + ) + + λ_offset = 0 + for (f, s) in nlp.oracles_data.oracles + # build local x for this oracle + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + + nout = s.set.output_dimension + μ = view(λ_oracle_all, (λ_offset + 1):(λ_offset + nout)) + + H_nnz = length(s.set.hessian_lagrangian_structure) + ind = (index + 1):(index + H_nnz) + s.set.eval_hessian_lagrangian(view(vals, ind), s.x, μ) + + index += H_nnz + λ_offset += nout + end + end + + @assert index == nlp.meta.nnzh return vals end @@ -645,9 +677,55 @@ function NLPModels.jth_hess_coord!( ) nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 0.0 end - # if !isempty(nlp.oracles) - # ... - # end + # check for oracle constraints + if !isempty(nlp.oracles_data.oracles) + n_nl = nlp.nlcon.nnln + noracle = nlp.oracles_data.noracle + + # check if j is an oracle constraint index + first_oracle = nlp.meta.nlin + nlp.quadcon.nquad + n_nl + 1 + last_oracle = first_oracle + noracle - 1 + + if first_oracle ≤ j ≤ last_oracle + # local oracle constraint index k ∈ 1: noracle + k = j - (nlp.meta.nlin + nlp.quadcon.nquad + n_nl) + + # starting index of this oracle's Hessian block + index = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh + + # walk through oracles to find which one owns component k + offset_outputs = 0 + for (f, s) in nlp.oracles_data.oracles + nout = s.set.output_dimension + H_nnz = length(s.set.hessian_lagrangian_structure) + + if k > offset_outputs + nout + # skip this oracle's block + index += H_nnz + offset_outputs += nout + continue + end + + # this oracle owns constraint k; local index ℓ + ℓ = k - offset_outputs + + # build local x + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + + # local multipliers μ: one-hot at ℓ + μ = zeros(eltype(x), nout) + μ[ℓ] = 1.0 + + # fill this oracle's Hessian block + ind = (index + 1):(index + H_nnz) + s.set.eval_hessian_lagrangian(view(vals, ind), s.x, μ) + + break + end + end + end return vals end From 699228fbc6eeb422bb88c38073cb5a6cfd0f2f5a Mon Sep 17 00:00:00 2001 From: Andrew Date: Tue, 16 Dec 2025 11:24:52 -0500 Subject: [PATCH 11/19] hs10 oracle tests --- test/nlp_consistency.jl | 2 +- test/nlp_problems/hs10.jl | 47 +++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 ++ test/test_moi_nlp_model.jl | 2 +- test/test_moi_nlp_oracle.jl | 21 +++++++++++++++++ 5 files changed, 72 insertions(+), 2 deletions(-) create mode 100644 test/test_moi_nlp_oracle.jl diff --git a/test/nlp_consistency.jl b/test/nlp_consistency.jl index b649357..a68ed06 100644 --- a/test/nlp_consistency.jl +++ b/test/nlp_consistency.jl @@ -1,4 +1,4 @@ -for problem in nlp_problems +for problem in Symbol.(lowercase.(nlp_problems ∪ extra_nlp_oracle_problems)) @testset "Problem $problem" begin nlp_manual = eval(Symbol(problem))() problem_f = eval(Symbol(lowercase(problem))) diff --git a/test/nlp_problems/hs10.jl b/test/nlp_problems/hs10.jl index 3e657e1..ceb69c3 100644 --- a/test/nlp_problems/hs10.jl +++ b/test/nlp_problems/hs10.jl @@ -12,3 +12,50 @@ function hs10() return nlp end + +# Same problem, but the nonlinear constraint is a VectorNonlinearOracle +function hs10_oracle() + model = Model() + + @variable(model, x[1:2]) + set_start_value(x[1], -10) + set_start_value(x[2], 10) + + @objective(model, Min, x[1] - x[2]) + + # g(x) = -3 x1^2 + 2 x1 x2 - x2^2 + 1 ≥ 0 + # Bounds: 0 ≤ g(x) ≤ +∞ + set = MOI.VectorNonlinearOracle(; + dimension = 2, # number of input variables (x1, x2) + l = [0.0], # lower bound on g(x) + u = [Inf], # upper bound on g(x) + eval_f = (ret, xv) -> begin + # ret[1] = g(x) + ret[1] = -3 * xv[1]^2 + 2 * xv[1] * xv[2] - xv[2]^2 + 1 + end, + # Jacobian of g(x): ∇g = [∂g/∂x1, ∂g/∂x2] + # ∂g/∂x1 = -6 x1 + 2 x2 + # ∂g/∂x2 = 2 x1 - 2 x2 + jacobian_structure = [(1, 1), (1, 2)], # (row, col) for each entry of ret + eval_jacobian = (ret, xv) -> begin + ret[1] = -6 * xv[1] + 2 * xv[2] # d g / d x1 + ret[2] = 2 * xv[1] - 2 * xv[2] # d g / d x2 + end, + # Hessian of g(x) (constant): + # ∂²g/∂x1² = -6 + # ∂²g/∂x1∂x2 = ∂²g/∂x2∂x1 = 2 + # ∂²g/∂x2² = -2 + hessian_lagrangian_structure = [(1, 1), (1, 2), (2, 2)], + eval_hessian_lagrangian = (ret, xv, μ) -> begin + # Hessian of μ[1] * g(x) + ret[1] = μ[1] * (-6.0) # (1,1) + ret[2] = μ[1] * ( 2.0) # (1,2) + ret[3] = μ[1] * (-2.0) # (2,2) + end, + ) + + # Same constraint, but expressed as a VectorNonlinearOracle on [x[1], x[2]] + @constraint(model, c, [x[1], x[2]] in set) + + return model +end diff --git a/test/runtests.jl b/test/runtests.jl index 7125c6d..2254815 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,8 @@ nls_problems = NLPModelsTest.nls_problems extra_nlp_problems = ["nohesspb", "hs61", "hs100", "hs219", "quadcon", "operatorspb", "nf"] extra_nls_problems = ["nlsnohesspb", "HS30", "HS43", "MGH07", "nlsqc"] +extra_nlp_oracle_problems = ["hs10_oracle"] + 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 bfafe05..536425a 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 ∪ extra_nlp_problems)) +for prob in Symbol.(lowercase.(nlp_problems ∪ extra_nlp_problems ∪ extra_nlp_oracle_problems)) prob_fn = eval(prob) nlp = MathOptNLPModel(prob_fn(), hessian = (prob != :nohesspb)) n = nlp.meta.nvar diff --git a/test/test_moi_nlp_oracle.jl b/test/test_moi_nlp_oracle.jl new file mode 100644 index 0000000..e13be42 --- /dev/null +++ b/test/test_moi_nlp_oracle.jl @@ -0,0 +1,21 @@ +@testset "Testing Oracle vs Non-Oracle NLP Models: $prob" for prob in extra_nlp_oracle_problems + prob_no_oracle = replace(prob, "_oracle" => "") + prob_fn_no_oracle = eval(Symbol(prob_no_oracle)) + prob_fn = eval(Symbol(prob)) + nlp_no_oracle = MathOptNLPModel(prob_fn_no_oracle(), hessian = true) + nlp_with_oracle = MathOptNLPModel(prob_fn(), hessian = true) + n = nlp_no_oracle.meta.nvar + m = nlp_no_oracle.meta.ncon + x = nlp_no_oracle.meta.x0 + fx_no_oracle = obj(nlp_no_oracle, x) + fx_with_oracle = obj(nlp_with_oracle, x) + @test isapprox(fx_no_oracle, fx_with_oracle; atol = 1e-8, rtol = 1e-8) + ngx_no_oracle = grad(nlp_no_oracle, x) + ngx_with_oracle = grad(nlp_with_oracle, x) + @test isapprox(ngx_no_oracle, ngx_with_oracle; atol = 1e-8, rtol = 1e-8) + if m > 0 + ncx_no_oracle = cons(nlp_no_oracle, x) + ncx_with_oracle = cons(nlp_with_oracle, x) + @test isapprox(ncx_no_oracle, ncx_with_oracle; rtol = 1e-1) + end +end \ No newline at end of file From 08c47ab2dd5730be6fd6baa146ac76355b18d470 Mon Sep 17 00:00:00 2001 From: Andrew Date: Tue, 16 Dec 2025 11:33:06 -0500 Subject: [PATCH 12/19] add hs14_oracle --- test/nlp_problems/hs10.jl | 2 +- test/nlp_problems/hs14.jl | 53 +++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 +- 3 files changed, 55 insertions(+), 2 deletions(-) diff --git a/test/nlp_problems/hs10.jl b/test/nlp_problems/hs10.jl index ceb69c3..10e098e 100644 --- a/test/nlp_problems/hs10.jl +++ b/test/nlp_problems/hs10.jl @@ -13,7 +13,7 @@ function hs10() return nlp end -# Same problem, but the nonlinear constraint is a VectorNonlinearOracle +"Problem 10 in the Hock-Schittkowski suite, but the nonlinear constraint is a VectorNonlinearOracle" function hs10_oracle() model = Model() diff --git a/test/nlp_problems/hs14.jl b/test/nlp_problems/hs14.jl index e349d98..8c68994 100644 --- a/test/nlp_problems/hs14.jl +++ b/test/nlp_problems/hs14.jl @@ -14,3 +14,56 @@ function hs14() return nlp end + +"Problem 14 in the Hock–Schittkowski suite, with inequality as a VectorNonlinearOracle" +function hs14_oracle() + model = Model() + + @variable(model, x[1:2]) + set_start_value(x[1], 2) + set_start_value(x[2], 2) + + @objective(model, Min, (x[1] - 2)^2 + (x[2] - 1)^2) + + # Inequality: -x1^2/4 - x2^2 + 1 ≥ 0 + # g(x) = -x1^2/4 - x2^2 + 1, 0 ≤ g(x) ≤ +∞ + set = MOI.VectorNonlinearOracle(; + dimension = 2, # 2 inputs: x1, x2 + l = [0.0], # lower bound on g(x) + u = [Inf], # upper bound on g(x) + eval_f = (ret, xv) -> begin + # ret[1] = g(x) + ret[1] = -0.25 * xv[1]^2 - xv[2]^2 + 1 + end, + # Jacobian of g(x): ∇g = [∂g/∂x1, ∂g/∂x2] + # ∂g/∂x1 = -0.5 x1 + # ∂g/∂x2 = -2 x2 + jacobian_structure = [(1, 1), (1, 2)], + eval_jacobian = (ret, xv) -> begin + ret[1] = -0.5 * xv[1] # d g / d x1 + ret[2] = -2.0 * xv[2] # d g / d x2 + end, + # Hessian of g(x): + # ∂²g/∂x1² = -0.5 + # ∂²g/∂x1∂x2 = 0 + # ∂²g/∂x2² = -2 + # + # We store only upper-triangular entries: (1,1), (1,2), (2,2) + hessian_lagrangian_structure = [(1, 1), (1, 2), (2, 2)], + eval_hessian_lagrangian = (ret, xv, μ) -> begin + # Hessian of μ[1] * g(x) + ret[1] = μ[1] * (-0.5) # (1,1) + ret[2] = μ[1] * 0.0 # (1,2) + ret[3] = μ[1] * (-2.0) # (2,2) + end, + ) + + # Inequality as oracle + @constraint(model, [x[1], x[2]] in set) + + # Equality stays as a standard nonlinear constraint: + # x1 - 2 x2 + 1 == 0 + @constraint(model, x[1] - 2 * x[2] + 1 == 0) + + return model +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 2254815..501d2fd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,7 +8,7 @@ nls_problems = NLPModelsTest.nls_problems extra_nlp_problems = ["nohesspb", "hs61", "hs100", "hs219", "quadcon", "operatorspb", "nf"] extra_nls_problems = ["nlsnohesspb", "HS30", "HS43", "MGH07", "nlsqc"] -extra_nlp_oracle_problems = ["hs10_oracle"] +extra_nlp_oracle_problems = ["hs10_oracle", "hs14_oracle"] for problem in lowercase.(nlp_problems ∪ extra_nlp_problems) include(joinpath("nlp_problems", "$problem.jl")) From 796df8524a4620103e04d9dd7aae20ea7778f59b Mon Sep 17 00:00:00 2001 From: Andrew Date: Tue, 16 Dec 2025 11:55:56 -0500 Subject: [PATCH 13/19] normalize constraints test --- test/nlp_problems/hs10.jl | 4 ++-- test/nlp_problems/hs14.jl | 4 ++-- test/runtests.jl | 2 ++ test/test_moi_nlp_oracle.jl | 2 +- 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/test/nlp_problems/hs10.jl b/test/nlp_problems/hs10.jl index 10e098e..5d29ad8 100644 --- a/test/nlp_problems/hs10.jl +++ b/test/nlp_problems/hs10.jl @@ -27,11 +27,11 @@ function hs10_oracle() # Bounds: 0 ≤ g(x) ≤ +∞ set = MOI.VectorNonlinearOracle(; dimension = 2, # number of input variables (x1, x2) - l = [0.0], # lower bound on g(x) + l = [-1.0], # lower bound on g(x) u = [Inf], # upper bound on g(x) eval_f = (ret, xv) -> begin # ret[1] = g(x) - ret[1] = -3 * xv[1]^2 + 2 * xv[1] * xv[2] - xv[2]^2 + 1 + ret[1] = -3 * xv[1]^2 + 2 * xv[1] * xv[2] - xv[2]^2 end, # Jacobian of g(x): ∇g = [∂g/∂x1, ∂g/∂x2] # ∂g/∂x1 = -6 x1 + 2 x2 diff --git a/test/nlp_problems/hs14.jl b/test/nlp_problems/hs14.jl index 8c68994..9eaa82d 100644 --- a/test/nlp_problems/hs14.jl +++ b/test/nlp_problems/hs14.jl @@ -29,11 +29,11 @@ function hs14_oracle() # g(x) = -x1^2/4 - x2^2 + 1, 0 ≤ g(x) ≤ +∞ set = MOI.VectorNonlinearOracle(; dimension = 2, # 2 inputs: x1, x2 - l = [0.0], # lower bound on g(x) + l = [-1.0], # lower bound on g(x) u = [Inf], # upper bound on g(x) eval_f = (ret, xv) -> begin # ret[1] = g(x) - ret[1] = -0.25 * xv[1]^2 - xv[2]^2 + 1 + ret[1] = -0.25 * xv[1]^2 - xv[2]^2 end, # Jacobian of g(x): ∇g = [∂g/∂x1, ∂g/∂x2] # ∂g/∂x1 = -0.5 x1 diff --git a/test/runtests.jl b/test/runtests.jl index 501d2fd..72d701e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,6 +21,8 @@ end include("test_moi_nlp_model.jl") include("test_moi_nls_model.jl") +include("test_moi_nlp_oracle.jl") + include("nlp_consistency.jl") include("nls_consistency.jl") diff --git a/test/test_moi_nlp_oracle.jl b/test/test_moi_nlp_oracle.jl index e13be42..0e6ede5 100644 --- a/test/test_moi_nlp_oracle.jl +++ b/test/test_moi_nlp_oracle.jl @@ -16,6 +16,6 @@ if m > 0 ncx_no_oracle = cons(nlp_no_oracle, x) ncx_with_oracle = cons(nlp_with_oracle, x) - @test isapprox(ncx_no_oracle, ncx_with_oracle; rtol = 1e-1) + @test isapprox(ncx_no_oracle, ncx_with_oracle; atol = 1e-8, rtol = 1e-8) end end \ No newline at end of file From 1c9caddbf62a55b14d5525acf721717730dca1d6 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Thu, 18 Dec 2025 11:19:20 -0500 Subject: [PATCH 14/19] More nlp tests [WIP] (#5) * hs10 oracle tests * add hs14_oracle * normalize constraints test * add more tests nlp * start checking jacobian and hessian * fix bug hessian * fix bg consistency --- src/moi_nlp_model.jl | 44 +++++++++++------- test/nlp_consistency.jl | 2 +- test/nlp_problems/hs100.jl | 90 +++++++++++++++++++++++++++++++++++++ test/nlp_problems/hs14.jl | 2 +- test/nlp_problems/hs61.jl | 73 ++++++++++++++++++++++++++++++ test/runtests.jl | 2 +- test/test_moi_nlp_oracle.jl | 41 +++++++++++++---- 7 files changed, 226 insertions(+), 28 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 19eb9b3..6a64809 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -312,7 +312,7 @@ function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals: s.x[i] = x[f.variables[i].value] end nnz_oracle = length(s.set.jacobian_structure) - s.set.eval_jacobian(view(values, (offset + 1):(offset + nnz_oracle)), s.x) + s.set.eval_jacobian(view(vals, (offset + 1):(offset + nnz_oracle)), s.x) offset += nnz_oracle end end @@ -350,8 +350,8 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs index += qcon.nnzg end end + offset = nlp.lincon.nnzj + nlp.quadcon.nnzj if nlp.meta.nnln > nlp.quadcon.nquad - offset = nlp.lincon.nnzj + nlp.quadcon.nnzj ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) offset += nlp.nlcon.nnzj @@ -360,7 +360,7 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs s.x[i] = x[f.variables[i].value] end nnz_oracle = length(s.set.jacobian_structure) - s.set.eval_jacobian(view(values, (offset + 1):(offset + nnz_oracle)), s.x) + s.set.eval_jacobian(view(vals, (offset + 1):(offset + nnz_oracle)), s.x) offset += nnz_oracle end end @@ -568,35 +568,45 @@ function NLPModels.hess_coord!( ) increment!(nlp, :neval_hess) - # Running index over Hessian nonzeros we've filled so far. - index = 0 - - # Objective Hessian block + # 1. Quadratic objective block (if any) if nlp.obj.type == "QUADRATIC" - view(vals, 1:(nlp.obj.nnzh)) .= obj_weight .* nlp.obj.hessian.vals - index += nlp.obj.nnzh + view(vals, 1:nlp.obj.nnzh) .= obj_weight .* nlp.obj.hessian.vals end - # Nonlinear objective Hessian block + + # 2. Nonlinear block (objective + JuMP @NLconstraint) if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) - λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)) + # Multipliers for the JuMP nonlinear constraints (not the oracles) + λ = view( + y, + (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon), + ) + + first_nl = nlp.obj.nnzh + nlp.quadcon.nnzh + 1 + last_nl = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh + MOI.eval_hessian_lagrangian( nlp.eval, - view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)), + view(vals, first_nl:last_nl), x, obj_weight, λ, ) end - # Quadratic constraint Hessian blocks + + # 3. Quadratic constraint Hessian blocks 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[nlp.meta.nlin + i] .* qcon.A.vals + ind = (index + 1):(index + qcon.nnzh) + view(vals, ind) .= y[nlp.meta.nlin + i] .* qcon.A.vals index += qcon.nnzh end end - # Oracle Hessian blocks + + # 4. Oracle Hessian blocks are appended at the very end if !isempty(nlp.oracles_data.oracles) + # Multipliers for oracle constraints only λ_oracle_all = view( y, (nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1): @@ -604,6 +614,9 @@ function NLPModels.hess_coord!( ) λ_offset = 0 + # Start after obj + quadcon + nonlinear block + index = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh + for (f, s) in nlp.oracles_data.oracles # build local x for this oracle for i in 1:s.set.input_dimension @@ -622,7 +635,6 @@ function NLPModels.hess_coord!( end end - @assert index == nlp.meta.nnzh return vals end diff --git a/test/nlp_consistency.jl b/test/nlp_consistency.jl index a68ed06..9d7d516 100644 --- a/test/nlp_consistency.jl +++ b/test/nlp_consistency.jl @@ -1,4 +1,4 @@ -for problem in Symbol.(lowercase.(nlp_problems ∪ extra_nlp_oracle_problems)) +for problem in [nlp_problems; extra_nlp_oracle_problems] @testset "Problem $problem" begin nlp_manual = eval(Symbol(problem))() problem_f = eval(Symbol(lowercase(problem))) diff --git a/test/nlp_problems/hs100.jl b/test/nlp_problems/hs100.jl index 92156cb..75e640b 100644 --- a/test/nlp_problems/hs100.jl +++ b/test/nlp_problems/hs100.jl @@ -22,3 +22,93 @@ function hs100(args...; kwargs...) return nlp end + +"HS100 model with the 2 middle constraints as VectorNonlinearOracle" +function hs100_oracle(args...; kwargs...) + model = Model() + x0 = [1, 2, 0, 4, 0, 1, 1] + @variable(model, x[i = 1:7], start = x0[i]) + + # 1st constraint: keep as NLconstraint + @NLconstraint(model, + 127 - 2 * x[1]^2 - 3 * x[2]^4 - x[3] - 4 * x[4]^2 - 5 * x[5] ≥ 0 + ) + + # 2nd constraint as oracle: + # Original: 282 - 7x1 - 3x2 - 10x3^2 - x4 + x5 ≥ 0 + # Canonical: f2(x) = -7x1 - 3x2 - 10x3^2 - x4 + x5, l2 = -282, u2 = +∞ + set2 = MOI.VectorNonlinearOracle(; + dimension = 5, # inputs: x1, x2, x3, x4, x5 + l = [-282.0], + u = [Inf], + eval_f = (ret, xv) -> begin + # xv = [x1, x2, x3, x4, x5] + ret[1] = -7.0 * xv[1] - 3.0 * xv[2] - 10.0 * xv[3]^2 - xv[4] + xv[5] + end, + # ∇f2 = [-7, -3, -20*x3, -1, 1] + jacobian_structure = [(1, 1), (1, 2), (1, 3), (1, 4), (1, 5)], + eval_jacobian = (ret, xv) -> begin + ret[1] = -7.0 # ∂f2/∂x1 + ret[2] = -3.0 # ∂f2/∂x2 + ret[3] = -20.0 * xv[3] # ∂f2/∂x3 + ret[4] = -1.0 # ∂f2/∂x4 + ret[5] = 1.0 # ∂f2/∂x5 + end, + # Hessian of f2: only (3,3) = -20 + hessian_lagrangian_structure = [(3, 3)], + eval_hessian_lagrangian = (ret, xv, μ) -> begin + # Hessian of μ[1] * f2(x) + ret[1] = μ[1] * (-20.0) # (3,3) + end, + ) + @constraint(model, [x[1], x[2], x[3], x[4], x[5]] in set2) + + # 3rd constraint as oracle: + # Original: -196 + 23x1 + x2^2 + 6x6^2 - 8x7 ≤ 0 + # Canonical: f3(x) = 23x1 + x2^2 + 6x6^2 - 8x7, l3 = -∞, u3 = 196 + set3 = MOI.VectorNonlinearOracle(; + dimension = 4, # inputs: x1, x2, x6, x7 + l = [-Inf], + u = [196.0], + eval_f = (ret, xv) -> begin + # xv = [x1, x2, x6, x7] + ret[1] = 23.0 * xv[1] + xv[2]^2 + 6.0 * xv[3]^2 - 8.0 * xv[4] + end, + # ∇f3 = [23, 2*x2, 12*x6, -8] + jacobian_structure = [(1, 1), (1, 2), (1, 3), (1, 4)], + eval_jacobian = (ret, xv) -> begin + ret[1] = 23.0 # ∂f3/∂x1 + ret[2] = 2.0 * xv[2] # ∂f3/∂x2 + ret[3] = 12.0 * xv[3] # ∂f3/∂x6 + ret[4] = -8.0 # ∂f3/∂x7 + end, + # Hessian of f3: (2,2) = 2, (3,3) = 12 + hessian_lagrangian_structure = [(2, 2), (3, 3)], + eval_hessian_lagrangian = (ret, xv, μ) -> begin + # Hessian of μ[1] * f3(x) + ret[1] = μ[1] * 2.0 # (2,2) + ret[2] = μ[1] * 12.0 # (3,3) + end, + ) + @constraint(model, [x[1], x[2], x[6], x[7]] in set3) + + # 4th constraint: keep as standard nonlinear + @constraint(model, + -4 * x[1]^2 - x[2]^2 + 3 * x[1] * x[2] - 2 * x[3]^2 - 5 * x[6] + 11 * x[7] ≥ 0 + ) + + # Objective: same as original + @NLobjective( + model, + 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 model +end diff --git a/test/nlp_problems/hs14.jl b/test/nlp_problems/hs14.jl index 9eaa82d..19b6588 100644 --- a/test/nlp_problems/hs14.jl +++ b/test/nlp_problems/hs14.jl @@ -66,4 +66,4 @@ function hs14_oracle() @constraint(model, x[1] - 2 * x[2] + 1 == 0) return model -end \ No newline at end of file +end diff --git a/test/nlp_problems/hs61.jl b/test/nlp_problems/hs61.jl index 5e422cc..384b223 100644 --- a/test/nlp_problems/hs61.jl +++ b/test/nlp_problems/hs61.jl @@ -10,3 +10,76 @@ function hs61(args...; kwargs...) return nlp end + +"HS61 model with both constraints as VectorNonlinearOracle" +function hs61_oracle(args...; kwargs...) + model = Model() + @variable(model, x[i = 1:3], start = 0) + + @NLobjective(model, Min, + 4 * x[1]^2 + 2 * x[2]^2 + 2 * x[3]^2 - + 33 * x[1] + 16 * x[2] - 24 * x[3] + ) + + # First equality: 3*x1 - 2*x2^2 - 7 == 0 + # Canonical form for MathOptNLPModel: + # f1(x) = 3*x1 - 2*x2^2 + # l1 = u1 = 7 + set1 = MOI.VectorNonlinearOracle(; + dimension = 2, # inputs: x1, x2 + l = [7.0], + u = [7.0], + eval_f = (ret, xv) -> begin + # xv = [x1, x2] + ret[1] = 3.0 * xv[1] - 2.0 * xv[2]^2 + end, + # ∇f1 = [3, -4*x2] + jacobian_structure = [(1, 1), (1, 2)], + eval_jacobian = (ret, xv) -> begin + ret[1] = 3.0 # d f1 / d x1 + ret[2] = -4.0 * xv[2] # d f1 / d x2 + end, + # Hessian of f1: + # ∂²f1/∂x1² = 0 + # ∂²f1/∂x2² = -4 + # (mixed derivatives are 0) + hessian_lagrangian_structure = [(2, 2)], + eval_hessian_lagrangian = (ret, xv, μ) -> begin + # Hessian of μ[1] * f1(x) + ret[1] = μ[1] * (-4.0) # (2,2) + end, + ) + + # Second equality: 4*x1 - x3^2 - 11 == 0 + # Canonical form: + # f2(x) = 4*x1 - x3^2 + # l2 = u2 = 11 + set2 = MOI.VectorNonlinearOracle(; + dimension = 2, # inputs: x1, x3 + l = [11.0], + u = [11.0], + eval_f = (ret, xv) -> begin + # xv = [x1, x3] + ret[1] = 4.0 * xv[1] - xv[2]^2 + end, + # ∇f2 = [4, -2*x3] + jacobian_structure = [(1, 1), (1, 2)], + eval_jacobian = (ret, xv) -> begin + ret[1] = 4.0 # d f2 / d x1 + ret[2] = -2.0 * xv[2] # d f2 / d x3 + end, + # Hessian of f2: + # ∂²f2/∂x1² = 0 + # ∂²f2/∂x3² = -2 + hessian_lagrangian_structure = [(2, 2)], + eval_hessian_lagrangian = (ret, xv, μ) -> begin + ret[1] = μ[1] * (-2.0) # (2,2) + end, + ) + + # Equality constraints as oracles + @constraint(model, [x[1], x[2]] in set1) + @constraint(model, [x[1], x[3]] in set2) + + return model +end diff --git a/test/runtests.jl b/test/runtests.jl index 72d701e..3f09d96 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,7 +8,7 @@ nls_problems = NLPModelsTest.nls_problems extra_nlp_problems = ["nohesspb", "hs61", "hs100", "hs219", "quadcon", "operatorspb", "nf"] extra_nls_problems = ["nlsnohesspb", "HS30", "HS43", "MGH07", "nlsqc"] -extra_nlp_oracle_problems = ["hs10_oracle", "hs14_oracle"] +extra_nlp_oracle_problems = ["hs10_oracle", "hs14_oracle", "hs61_oracle", "hs100_oracle"] for problem in lowercase.(nlp_problems ∪ extra_nlp_problems) include(joinpath("nlp_problems", "$problem.jl")) diff --git a/test/test_moi_nlp_oracle.jl b/test/test_moi_nlp_oracle.jl index 0e6ede5..d945995 100644 --- a/test/test_moi_nlp_oracle.jl +++ b/test/test_moi_nlp_oracle.jl @@ -2,20 +2,43 @@ prob_no_oracle = replace(prob, "_oracle" => "") prob_fn_no_oracle = eval(Symbol(prob_no_oracle)) prob_fn = eval(Symbol(prob)) - nlp_no_oracle = MathOptNLPModel(prob_fn_no_oracle(), hessian = true) - nlp_with_oracle = MathOptNLPModel(prob_fn(), hessian = true) + + nlp_no_oracle = MathOptNLPModel(prob_fn_no_oracle(), hessian = true) + nlp_with_oracle = MathOptNLPModel(prob_fn(), hessian = true) + n = nlp_no_oracle.meta.nvar m = nlp_no_oracle.meta.ncon x = nlp_no_oracle.meta.x0 - fx_no_oracle = obj(nlp_no_oracle, x) + + # Objective value + fx_no_oracle = obj(nlp_no_oracle, x) fx_with_oracle = obj(nlp_with_oracle, x) @test isapprox(fx_no_oracle, fx_with_oracle; atol = 1e-8, rtol = 1e-8) - ngx_no_oracle = grad(nlp_no_oracle, x) + + # Gradient of objective + ngx_no_oracle = grad(nlp_no_oracle, x) ngx_with_oracle = grad(nlp_with_oracle, x) @test isapprox(ngx_no_oracle, ngx_with_oracle; atol = 1e-8, rtol = 1e-8) - if m > 0 - ncx_no_oracle = cons(nlp_no_oracle, x) - ncx_with_oracle = cons(nlp_with_oracle, x) - @test isapprox(ncx_no_oracle, ncx_with_oracle; atol = 1e-8, rtol = 1e-8) - end + + # Constraint values (up to ordering) + ncx_no_oracle = cons(nlp_no_oracle, x) + ncx_with_oracle = cons(nlp_with_oracle, x) + @test isapprox(sort(ncx_no_oracle), sort(ncx_with_oracle); + atol = 1e-8, rtol = 1e-8) + + # Jacobian: compare J'J, which is invariant to row permutations + J_no = jac(nlp_no_oracle, x) + J_with = jac(nlp_with_oracle, x) + + G_no = Matrix(J_no)' * Matrix(J_no) + G_with = Matrix(J_with)' * Matrix(J_with) + + @test isapprox(G_no, G_with; atol = 1e-8, rtol = 1e-8) + + # Hessian of the objective: use y = 0 so constraints don't enter + λ = zeros(m) + H_no = hess(nlp_no_oracle, x, λ) + H_with = hess(nlp_with_oracle, x, λ) + + @test isapprox(Matrix(H_no), Matrix(H_with); atol = 1e-8, rtol = 1e-8) end \ No newline at end of file From 345b20e707f35a4d0c3479ab341699e7a83ac11c Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Fri, 19 Dec 2025 02:33:50 +0100 Subject: [PATCH 15/19] Full pass on the support of MOI.VectorNonlinearOracle --- Project.toml | 2 +- src/moi_nlp_model.jl | 433 ++++++++++++++++++------------------ src/moi_nls_model.jl | 34 +-- src/utils.jl | 60 ++--- test/nlp_consistency.jl | 7 +- test/nlp_problems/hs10.jl | 2 +- test/test_moi_nlp_oracle.jl | 27 ++- 7 files changed, 283 insertions(+), 282 deletions(-) diff --git a/Project.toml b/Project.toml index e7fb9e0..84a70c6 100644 --- a/Project.toml +++ b/Project.toml @@ -16,7 +16,7 @@ JuMP = "1.25" LinearAlgebra = "1.10" MathOptInterface = "1.46" NLPModels = "0.21.6" -NLPModelsTest = "0.10.3" +NLPModelsTest = "0.10.5" Percival = "0.7.3" Printf = "1.10" SolverCore = "0.3" diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 6a64809..872c574 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -6,7 +6,7 @@ mutable struct MathOptNLPModel <: AbstractNLPModel{Float64, Vector{Float64}} lincon::LinearConstraints quadcon::QuadraticConstraints nlcon::NonLinearStructure - oracles_data::OraclesData + oracles::Oracles λ::Vector{Float64} hv::Vector{Float64} obj::Objective @@ -35,10 +35,10 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = parser_MOI(moimodel, index_map, nvar) nlp_data = _nlp_block(moimodel) - nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) - oracles_data = parser_oracles(moimodel) + nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) + oracles, hessian_oracles_supported = parser_oracles(moimodel) counters = Counters() - λ = zeros(Float64, nnln) # Lagrange multipliers for hess_coord! and hprod! without y + λ = zeros(Float64, nlcon.nnln) # Lagrange multipliers for hess_coord! and hprod! without y hv = zeros(Float64, nvar) # workspace for ghjvprod! if nlp_data.has_objective @@ -48,11 +48,11 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = end # Total counts - ncon = nlin + quadcon.nquad + nnln + oracles_data.noracle - lcon = vcat(lin_lcon, quad_lcon, nl_lcon, oracles_data.l_oracle) - ucon = vcat(lin_ucon, quad_ucon, nl_ucon, oracles_data.u_oracle) - nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj + oracles_data.nnzj_oracle - nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh + oracles_data.nnzh_oracle + ncon = nlin + quadcon.nquad + nlcon.nnln + oracles.ncon + lcon = vcat(lin_lcon, quad_lcon, nl_lcon, oracles.lcon) + ucon = vcat(lin_ucon, quad_ucon, nl_ucon, oracles.ucon) + nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj + oracles.nnzj + nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh + oracles.nnzh meta = NLPModelMeta( nvar, x0 = x0, @@ -66,17 +66,17 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = nnzh = nnzh, lin = collect(1:nlin), lin_nnzj = lincon.nnzj, - nln_nnzj = quadcon.nnzj + nlcon.nnzj + oracles_data.nnzj_oracle, + nln_nnzj = quadcon.nnzj + nlcon.nnzj + oracles.nnzj, minimize = MOI.get(moimodel, MOI.ObjectiveSense()) == MOI.MIN_SENSE, - islp = (obj.type == "LINEAR") && (nnln == 0) && (quadcon.nquad == 0), + islp = (obj.type == "LINEAR") && (nlin == ncon), name = name, - jprod_available = isempty(oracles_data.oracles), - jtprod_available = isempty(oracles_data.oracles), - hprod_available = isempty(oracles_data.oracles) && hessian, - hess_available = hessian, + jprod_available = isempty(oracles.oracles), + jtprod_available = isempty(oracles.oracles), + hprod_available = isempty(oracles.oracles) && hessian, + hess_available = hessian && hessian_oracles_supported, ) - return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, oracles_data, λ, hv, obj, counters), + return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, oracles, λ, hv, obj, counters), index_map end @@ -120,16 +120,21 @@ end function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVector) increment!(nlp, :neval_cons_nln) + offset = 0 if nlp.quadcon.nquad > 0 + offset += nlp.quadcon.nquad for i = 1:(nlp.quadcon.nquad) 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 end - if nlp.meta.nnln > nlp.quadcon.nquad - offset = nlp.quadcon.nquad + nlp.nlcon.nnln - MOI.eval_constraint(nlp.eval, view(c, (nlp.quadcon.nquad + 1):(offset)), x) - for (f, s) in nlp.oracles_data.oracles + if nlp.nlcon.nnln > 0 + offset += nlp.nlcon.nnln + index_nlcon = (nlp.quadcon.nquad + 1):(offset) + MOI.eval_constraint(nlp.eval, view(c, ), x) + end + if nlp.oracles.ncon > 0 + for (f, s) in nlp.oracles.oracles for i in 1:s.set.input_dimension s.x[i] = x[f.variables[i].value] end @@ -139,36 +144,40 @@ function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::Abstrac offset += s.set.output_dimension end end + @assert offset == nlp.meta.nnln return c end function NLPModels.cons!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVector) increment!(nlp, :neval_cons) + offset = 0 if nlp.meta.nlin > 0 + offset += nlp.meta.nlin coo_prod!(nlp.lincon.jacobian.rows, nlp.lincon.jacobian.cols, nlp.lincon.jacobian.vals, x, c) end - if nlp.meta.nnln > 0 - if nlp.quadcon.nquad > 0 - for i = 1:(nlp.quadcon.nquad) - qcon = nlp.quadcon.constraints[i] - c[nlp.meta.nlin + i] = - 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x) - end + if nlp.quadcon.nquad > 0 + offset += nlp.quadcon.nquad + for i = 1:(nlp.quadcon.nquad) + qcon = nlp.quadcon.constraints[i] + c[nlp.meta.nlin + 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 - offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln - MOI.eval_constraint(nlp.eval, view(c, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(offset)), x) - for (f, s) in nlp.oracles_data.oracles - for i in 1:s.set.input_dimension - s.x[i] = x[f.variables[i].value] - end - index_oracle = (offset + 1):(offset + s.set.output_dimension) - c_oracle = view(c, index_oracle) - s.set.eval_f(c_oracle, s.x) - offset += s.set.output_dimension + end + if nlp.nlcon.nnln > 0 + offset += nlp.nlcon.nnln + MOI.eval_constraint(nlp.eval, view(c, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(offset)), x) + end + if nlp.oracles.ncon > 0 + for (f, s) in nlp.oracles.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] end + index_oracle = (offset + 1):(offset + s.set.output_dimension) + c_oracle = view(c, index_oracle) + s.set.eval_f(c_oracle, s.x) + offset += s.set.output_dimension end end + @assert offset == nlp.meta.ncon return c end @@ -188,40 +197,39 @@ function NLPModels.jac_nln_structure!( rows::AbstractVector{<:Integer}, cols::AbstractVector{<:Integer}, ) + offset = 0 if nlp.quadcon.nquad > 0 - index = 0 for i = 1:(nlp.quadcon.nquad) # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon qcon = nlp.quadcon.constraints[i] - ind_quad = (index + 1):(index + qcon.nnzg) + ind_quad = (offset + 1):(offset + qcon.nnzg) view(rows, ind_quad) .= i view(cols, ind_quad) .= qcon.g - index += qcon.nnzg + offset += qcon.nnzg end end - if nlp.meta.nnln > nlp.quadcon.nquad - offset = nlp.quadcon.nnzj - # non-oracle nonlinear constraints + @assert offset == nlp.quadcon.nnzj + if nlp.nlcon.nnln > 0 + # non-oracle nonlinear constraints ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nlp.nlcon.jac_rows view(cols, ind_nnln) .= nlp.nlcon.jac_cols offset += nlp.nlcon.nnzj - # structure of oracle Jacobians - for (k, (f, s)) in enumerate(nlp.oracles_data.oracles) - # Shift row index by quadcon.nquad - # plus previous oracle outputs. - row_offset = nlp.quadcon.nquad - for i in 1:(k - 1) - row_offset += nlp.oracles_data.oracles[i][2].set.output_dimension - end - - for (r, c) in s.set.jacobian_structure - offset += 1 - rows[offset] = row_offset + r - cols[offset] = f.variables[c].value - end + end + if nlp.oracles.ncon > 0 + # oracle nonlinear constraints + row_offset = nlp.quadcon.nquad + nlp.nlcon.nnln + for (f, s) in nlp.oracles.oracles + for (r, c) in s.set.jacobian_structure + offset += 1 + rows[offset] = row_offset + r + cols[offset] = f.variables[c].value + end + row_offset += s.set.output_dimension end + @assert row_offset == nlp.meta.nnln end + @assert offset == nlp.meta.nln_nnzj return rows, cols end @@ -230,42 +238,46 @@ function NLPModels.jac_structure!( rows::AbstractVector{<:Integer}, cols::AbstractVector{<:Integer}, ) + offset = 0 if nlp.meta.nlin > 0 index_lin = 1:(nlp.lincon.nnzj) view(rows, index_lin) .= nlp.lincon.jacobian.rows view(cols, index_lin) .= nlp.lincon.jacobian.cols + offset += nlp.lincon.nnzj end - if nlp.meta.nnln > 0 - if nlp.quadcon.nquad > 0 - index = 0 - for i = 1:(nlp.quadcon.nquad) - # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon - qcon = nlp.quadcon.constraints[i] - ind_quad = (nlp.lincon.nnzj + index + 1):(nlp.lincon.nnzj + index + qcon.nnzg) - view(rows, ind_quad) .= nlp.meta.nlin .+ i - view(cols, ind_quad) .= qcon.g - index += qcon.nnzg - end + if nlp.quadcon.nquad > 0 + for i = 1:(nlp.quadcon.nquad) + # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon + qcon = nlp.quadcon.constraints[i] + ind_quad = (offset + 1):(offset + qcon.nnzg) + view(rows, ind_quad) .= nlp.meta.nlin .+ i + view(cols, ind_quad) .= qcon.g + offset += qcon.nnzg end - if nlp.meta.nnln > nlp.quadcon.nquad - # non-oracle nonlinear constraints - offset = nlp.lincon.nnzj + nlp.quadcon.nnzj - ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) - view(rows, ind_nnln) .= nlp.meta.nlin .+ nlp.quadcon.nquad .+ nlp.nlcon.jac_rows - view(cols, ind_nnln) .= nlp.nlcon.jac_cols - offset += nlp.nlcon.nnzj - # structure of oracle Jacobians - row_offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln - for (f, s) in nlp.oracles_data.oracles - for (i, j) in s.set.jacobian_structure - offset += 1 - rows[offset] = row_offset + i - cols[offset] = f.variables[j].value - end - row_offset += s.set.output_dimension + end + @assert offset == nlp.lincon.nnzj + nlp.quadcon.nnzj + if nlp.nlcon.nnln > 0 + # non-oracle nonlinear constraints + offset = nlp.lincon.nnzj + nlp.quadcon.nnzj + ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) + view(rows, ind_nnln) .= nlp.meta.nlin .+ nlp.quadcon.nquad .+ nlp.nlcon.jac_rows + view(cols, ind_nnln) .= nlp.nlcon.jac_cols + offset += nlp.nlcon.nnzj + end + if nlp.oracles.ncon > 0 + # oracle nonlinear constraints + row_offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + for (f, s) in nlp.oracles.oracles + for (r, c) in s.set.jacobian_structure + offset += 1 + rows[offset] = row_offset + r + cols[offset] = f.variables[c].value end + row_offset += s.set.output_dimension end + @assert row_offset == nlp.meta.ncon end + @assert offset == nlp.meta.nnzj return rows, cols end @@ -278,91 +290,95 @@ end function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::AbstractVector) increment!(nlp, :neval_jac_nln) + offset = 0 if nlp.quadcon.nquad > 0 - index = 0 ind_quad = 1:(nlp.quadcon.nnzj) view(vals, ind_quad) .= 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] + vals[offset + 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] + vals[offset + k1] += val * x[col] if row != col k2 = qcon.dg[col] - vals[index + k2] += val * x[row] + vals[offset + k2] += val * x[row] end end - index += qcon.nnzg + offset += qcon.nnzg end end - if nlp.meta.nnln > nlp.quadcon.nquad - offset = nlp.quadcon.nnzj + if nlp.nlcon.nnln > 0 ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) offset += nlp.nlcon.nnzj - for (f, s) in nlp.oracles_data.oracles + end + if nlp.oracles.ncon > 0 + for (f, s) in nlp.oracles.oracles for i in 1:s.set.input_dimension s.x[i] = x[f.variables[i].value] end nnz_oracle = length(s.set.jacobian_structure) - s.set.eval_jacobian(view(vals, (offset + 1):(offset + nnz_oracle)), s.x) + ind_oracle = (offset + 1):(offset + nnz_oracle) + s.set.eval_jacobian(view(vals, ind_oracle), s.x) offset += nnz_oracle end end + @assert offset == nlp.meta.nln_nnzj return vals end function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::AbstractVector) increment!(nlp, :neval_jac) + offset = 0 if nlp.meta.nlin > 0 index_lin = 1:(nlp.lincon.nnzj) view(vals, index_lin) .= nlp.lincon.jacobian.vals + offset += nlp.lincon.nnzj end - if nlp.meta.nnln > 0 - if nlp.quadcon.nquad > 0 - ind_quad = (nlp.lincon.nnzj + 1):(nlp.lincon.nnzj + nlp.quadcon.nnzj) - view(vals, ind_quad) .= 0.0 - index = nlp.lincon.nnzj - 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 + if nlp.quadcon.nquad > 0 + ind_quad = (nlp.lincon.nnzj + 1):(nlp.lincon.nnzj + nlp.quadcon.nnzj) + view(vals, ind_quad) .= 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[offset + 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[offset + k1] += val * x[col] + if row != col + k2 = qcon.dg[col] + vals[offset + k2] += val * x[row] end - index += qcon.nnzg end + offset += qcon.nnzg end - offset = nlp.lincon.nnzj + nlp.quadcon.nnzj - if nlp.meta.nnln > nlp.quadcon.nquad - ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) - MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) - offset += nlp.nlcon.nnzj - for (f, s) in nlp.oracles_data.oracles - for i in 1:s.set.input_dimension - s.x[i] = x[f.variables[i].value] - end - nnz_oracle = length(s.set.jacobian_structure) - s.set.eval_jacobian(view(vals, (offset + 1):(offset + nnz_oracle)), s.x) - offset += nnz_oracle + end + if nlp.nlcon.nnln > 0 + ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj) + MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x) + offset += nlp.nlcon.nnzj + end + if nlp.oracles.ncon > 0 + for (f, s) in nlp.oracles.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] end + nnz_oracle = length(s.set.jacobian_structure) + ind_oracle = (offset + 1):(offset + nnz_oracle) + s.set.eval_jacobian(view(vals, ind_oracle), s.x) + offset += nnz_oracle end end @assert offset == nlp.meta.nnzj @@ -403,7 +419,7 @@ function NLPModels.jprod_nln!( Jv[i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) end end - if nlp.meta.nnln > nlp.quadcon.nquad + if nlp.nlcon.nnln > 0 ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln) MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v) end @@ -439,7 +455,7 @@ function NLPModels.jprod!( coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) end end - if nlp.meta.nnln > nlp.quadcon.nquad + if nlp.nlcon.nnln > 0 ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon) MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v) end @@ -473,7 +489,7 @@ function NLPModels.jtprod_nln!( ) nlp.meta.jtprod_available || error("J' * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jtprod_nln) - if nlp.meta.nnln > nlp.quadcon.nquad + if nlp.nlcon.nnln > 0 ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln) MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln)) end @@ -497,7 +513,7 @@ function NLPModels.jtprod!( ) nlp.meta.jtprod_available || error("J' * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jtprod) - if nlp.meta.nnln > nlp.quadcon.nquad + if nlp.nlcon.nnln > 0 ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon) MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln)) end @@ -542,17 +558,17 @@ function NLPModels.hess_structure!( index += qcon.nnzh end end - if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) + if (nlp.obj.type == "NONLINEAR") || (nlp.nlcon.nnln > 0) view(rows, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_rows view(cols, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_cols index += nlp.nlcon.nnzh end - if !isempty(nlp.oracles_data.oracles) - for (f, s) in nlp.oracles_data.oracles - for (i_local, j_local) in s.set.hessian_lagrangian_structure + if nlp.oracles.ncon > 0 + for (f, s) in nlp.oracles.oracles + for (i, j) in s.set.hessian_lagrangian_structure index += 1 - rows[index] = f.variables[i_local].value - cols[index] = f.variables[j_local].value + rows[index] = f.variables[i].value + cols[index] = f.variables[j].value end end end @@ -574,22 +590,17 @@ function NLPModels.hess_coord!( end # 2. Nonlinear block (objective + JuMP @NLconstraint) - if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) + if (nlp.obj.type == "NONLINEAR") || (nlp.nlcon.nnln > 0) # Multipliers for the JuMP nonlinear constraints (not the oracles) - λ = view( - y, - (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon), - ) - - first_nl = nlp.obj.nnzh + nlp.quadcon.nnzh + 1 - last_nl = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh - + ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln) + λ_nnln = view(y, ind_nnln) + ind_nnzh = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh) MOI.eval_hessian_lagrangian( nlp.eval, - view(vals, first_nl:last_nl), + view(vals, ind_nnzh), x, obj_weight, - λ, + λ_nnln, ) end @@ -605,34 +616,29 @@ function NLPModels.hess_coord!( end # 4. Oracle Hessian blocks are appended at the very end - if !isempty(nlp.oracles_data.oracles) - # Multipliers for oracle constraints only - λ_oracle_all = view( - y, - (nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1): - (nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + nlp.oracles_data.noracle), - ) + if nlp.oracles.ncon > 0 + offset_λ = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + offset_nnzh = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh - λ_offset = 0 - # Start after obj + quadcon + nonlinear block - index = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh - - for (f, s) in nlp.oracles_data.oracles + for (f, s) in nlp.oracles.oracles # build local x for this oracle for i in 1:s.set.input_dimension s.x[i] = x[f.variables[i].value] end - nout = s.set.output_dimension - μ = view(λ_oracle_all, (λ_offset + 1):(λ_offset + nout)) + ncon_oracle = s.set.output_dimension + index_λ = (offset_λ + 1):(offset_λ + ncon_oracle) + λ_oracles = view(y, index_λ) - H_nnz = length(s.set.hessian_lagrangian_structure) - ind = (index + 1):(index + H_nnz) - s.set.eval_hessian_lagrangian(view(vals, ind), s.x, μ) + nnzh_oracle = length(s.set.hessian_lagrangian_structure) + index_nnzh = (offset_nnzh + 1):(offset_nnzh + nnzh_oracle) + s.set.eval_hessian_lagrangian(view(vals, index_nnzh), s.x, λ_oracles) - index += H_nnz - λ_offset += nout + offset_λ += ncon_oracle + offset_nnzh += nnzh_oracle end + @assert offset_λ == nlp.meta.ncon + @assert offset_nnzh == nlp.meta.nnzh end return vals @@ -669,6 +675,8 @@ function NLPModels.jth_hess_coord!( increment!(nlp, :neval_jhess) @rangecheck 1 nlp.meta.ncon j vals .= 0.0 + + # Quadratic constraints if nlp.meta.nlin + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad index = nlp.obj.nnzh for i = 1:(nlp.quadcon.nquad) @@ -678,63 +686,50 @@ function NLPModels.jth_hess_coord!( end index += qcon.nnzh end - elseif nlp.meta.nlin + nlp.quadcon.nquad + 1 ≤ j ≤ nlp.meta.ncon + end + + # Non-oracle nonlinear constraints + if nlp.meta.nlin + nlp.quadcon.nquad + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 1.0 + index_nnzh = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh) MOI.eval_hessian_lagrangian( nlp.eval, - view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)), + view(vals, index_nnzh), x, 0.0, nlp.λ, ) nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 0.0 end - # check for oracle constraints - if !isempty(nlp.oracles_data.oracles) - n_nl = nlp.nlcon.nnln - noracle = nlp.oracles_data.noracle - - # check if j is an oracle constraint index - first_oracle = nlp.meta.nlin + nlp.quadcon.nquad + n_nl + 1 - last_oracle = first_oracle + noracle - 1 - - if first_oracle ≤ j ≤ last_oracle - # local oracle constraint index k ∈ 1: noracle - k = j - (nlp.meta.nlin + nlp.quadcon.nquad + n_nl) - - # starting index of this oracle's Hessian block - index = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh - - # walk through oracles to find which one owns component k - offset_outputs = 0 - for (f, s) in nlp.oracles_data.oracles - nout = s.set.output_dimension - H_nnz = length(s.set.hessian_lagrangian_structure) - - if k > offset_outputs + nout - # skip this oracle's block - index += H_nnz - offset_outputs += nout - continue - end - - # this oracle owns constraint k; local index ℓ - ℓ = k - offset_outputs - # build local x - for i in 1:s.set.input_dimension - s.x[i] = x[f.variables[i].value] + # Oracle nonlinear constraints + if nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1 ≤ j ≤ nlp.meta.ncon + offset_λ = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + offset_nnzh = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh + find_oracle = false + for (f, s) in nlp.oracles.oracles + if !find_oracle + ncon_oracle = s.set.output_dimension + nnzh_oracle = length(s.set.hessian_lagrangian_structure) + if offset_λ + 1 ≤ j ≤ offset_λ + find_oracle = true + # build local x + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + # We could preallocate the vector in the structure nlp + λ_oracle = zeros(Float64, s.set.output_dimension) + λ_oracle[j - offset_λ] = 1.0 + index_nnzh = (offset_nnzh + 1):(offset_nnzh + nnz_oracle) + s.set.eval_hessian_lagrangian( + view(vals, index_nnzh), + s.x, + λ_oracle, + ) + else + offset_λ += ncon_oracle + offset_nnzh += nnzh_oracle end - - # local multipliers μ: one-hot at ℓ - μ = zeros(eltype(x), nout) - μ[ℓ] = 1.0 - - # fill this oracle's Hessian block - ind = (index + 1):(index + H_nnz) - s.set.eval_hessian_lagrangian(view(vals, ind), s.x, μ) - - break end end end @@ -749,12 +744,12 @@ function NLPModels.hprod!( hv::AbstractVector; obj_weight::Float64 = 1.0, ) - nlp.meta.hprod_available && error("H * v products are not supported by this MathOptNLPModel.") + nlp.meta.hprod_available || error("H * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_hprod) if (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == 0) hv .= 0.0 end - if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad) + if (nlp.obj.type == "NONLINEAR") || (nlp.nlcon.nnln > 0) λ = 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 @@ -786,7 +781,7 @@ function NLPModels.hprod!( hv::AbstractVector; obj_weight::Float64 = 1.0, ) - nlp.meta.hprod_available && error("H * v products are not supported by this MathOptNLPModel.") + nlp.meta.hprod_available || error("H * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_hprod) if nlp.obj.type == "LINEAR" hv .= 0.0 @@ -815,7 +810,7 @@ function NLPModels.jth_hprod!( j::Integer, hv::AbstractVector, ) - nlp.meta.hprod_available && error("H * v products are not supported by this MathOptNLPModel.") + nlp.meta.hprod_available || error("H * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jhprod) @rangecheck 1 nlp.meta.ncon j hv .= 0.0 @@ -837,7 +832,7 @@ function NLPModels.ghjvprod!( v::AbstractVector, ghv::AbstractVector, ) - nlp.meta.hprod_available && error("The operation is not supported by this MathOptNLPModel.") + nlp.meta.hprod_available || error("The operation is not supported by this MathOptNLPModel.") increment!(nlp, :neval_hprod) ghv .= 0.0 for i = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad) diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl index 8eeccf9..a1846dd 100644 --- a/src/moi_nls_model.jl +++ b/src/moi_nls_model.jl @@ -36,14 +36,14 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri parser_MOI(moimodel, index_map, nvar) nlp_data = _nlp_block(moimodel) - nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) - λ = zeros(nnln) # Lagrange multipliers for hess_coord! and hprod! without y + nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) + λ = zeros(nlcon.nnln) # Lagrange multipliers for hess_coord! and hprod! without y nequ = nlinequ + nnlnequ Fnnzj = linequ.nnzj + nlequ.nnzj Fnnzh = nlequ.nnzh - ncon = nlin + quadcon.nquad + nnln + ncon = nlin + quadcon.nquad + nlcon.nnln lcon = vcat(lin_lcon, quad_lcon, nl_lcon) ucon = vcat(lin_ucon, quad_ucon, nl_ucon) cnnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj @@ -265,7 +265,7 @@ function NLPModels.cons_nln!(nls::MathOptNLSModel, x::AbstractVector, c::Abstrac c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x) end end - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 MOI.eval_constraint(nls.ceval, view(c, (nls.quadcon.nquad + 1):(nls.meta.nnln)), x) end return c @@ -284,7 +284,7 @@ function NLPModels.cons!(nls::MathOptNLSModel, x::AbstractVector, c::AbstractVec 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x) end end - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 index_nnln = (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon) MOI.eval_constraint(nls.ceval, view(c, index_nnln), x) end @@ -319,7 +319,7 @@ function NLPModels.jac_nln_structure!( index += qcon.nnzg end end - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj) view(rows, ind_nnln) .= nls.quadcon.nquad .+ nls.nlcon.jac_rows view(cols, ind_nnln) .= nls.nlcon.jac_cols @@ -349,7 +349,7 @@ function NLPModels.jac_structure!( index += qcon.nnzg end end - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 ind_nnln = (nls.lincon.nnzj + nls.quadcon.nnzj + 1):(nls.meta.nnzj) view(rows, ind_nnln) .= nls.meta.nlin .+ nls.quadcon.nquad .+ nls.nlcon.jac_rows view(cols, ind_nnln) .= nls.nlcon.jac_cols @@ -391,7 +391,7 @@ function NLPModels.jac_nln_coord!(nls::MathOptNLSModel, x::AbstractVector, vals: index += qcon.nnzg end end - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj) MOI.eval_constraint_jacobian(nls.ceval, view(vals, ind_nnln), x) end @@ -429,7 +429,7 @@ function NLPModels.jac_coord!(nls::MathOptNLSModel, x::AbstractVector, vals::Abs index += qcon.nnzg end end - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 ind_nnln = (nls.lincon.nnzj + nls.quadcon.nnzj + 1):(nls.meta.nnzj) MOI.eval_constraint_jacobian(nls.ceval, view(vals, ind_nnln), x) end @@ -469,7 +469,7 @@ function NLPModels.jprod_nln!( Jv[i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) end end - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 ind_nnln = (nls.quadcon.nquad + 1):(nls.meta.nnln) MOI.eval_constraint_jacobian_product(nls.ceval, view(Jv, ind_nnln), x, v) end @@ -504,7 +504,7 @@ function NLPModels.jprod!( coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v) end end - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 ind_nnln = (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon) MOI.eval_constraint_jacobian_product(nls.ceval, view(Jv, ind_nnln), x, v) end @@ -536,7 +536,7 @@ function NLPModels.jtprod_nln!( Jtv::AbstractVector, ) increment!(nls, :neval_jtprod_nln) - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 ind_nnln = (nls.quadcon.nquad + 1):(nls.meta.nnln) MOI.eval_constraint_jacobian_transpose_product(nls.ceval, Jtv, x, view(v, ind_nnln)) end @@ -559,7 +559,7 @@ function NLPModels.jtprod!( Jtv::AbstractVector, ) increment!(nls, :neval_jtprod) - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 ind_nnln = (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon) MOI.eval_constraint_jacobian_transpose_product(nls.ceval, Jtv, x, view(v, ind_nnln)) end @@ -595,7 +595,7 @@ 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 > nls.quadcon.nquad) + if (nls.nls_meta.nnln > 0) || (nls.nlcon.nnln > 0) 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 @@ -622,7 +622,7 @@ 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 > nls.quadcon.nquad) + if (nls.nls_meta.nnln > 0) || (nls.nlcon.nnln > 0) λ = view(y, (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon)) MOI.eval_hessian_lagrangian( nls.ceval, @@ -658,7 +658,7 @@ function NLPModels.hess_coord!( 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 - if nls.meta.nnln > nls.quadcon.nquad + if nls.nlcon.nnln > 0 view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= 0.0 end end @@ -674,7 +674,7 @@ function NLPModels.hprod!( obj_weight::Float64 = 1.0, ) increment!(nls, :neval_hprod) - if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad) + if (nls.nls_meta.nnln > 0) || (nls.nlcon.nnln > 0) λ = 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 diff --git a/src/utils.jl b/src/utils.jl index 54f29c0..acc8dd0 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -24,6 +24,9 @@ const SNF = MOI.ScalarNonlinearFunction const VNF = MOI.VectorNonlinearFunction const NF = Union{SNF, VNF} +# VectorNonlinearOracle +const ORACLE = MOI.VectorNonlinearOracle{Float64} # VectorNonlinearOracle{Float64}(input_dimension, output_dimension, l, u, eval_f, jacobian_structure, eval_jacobian, hessian_lagrangian_structure, eval_hessian_lagrangian) + # Cache of VectorNonlinearOracle mutable struct _VectorNonlinearOracleCache set::MOI.VectorNonlinearOracle{Float64} @@ -121,23 +124,23 @@ mutable struct Objective end """ - OraclesData + Oracles Structure containing nonlinear oracles data: - oracles: vector of tuples (MOI.VectorOfVariables, _VectorNonlinearOracleCache) -- noracle: number of scalar constraints represented by all oracles -- l_oracle: lower bounds of oracle constraints -- u_oracle: upper bounds of oracle constraints -- nnzj_oracle: number of non-zero entries in the Jacobian of all oracles -- nnzh_oracle: number of non-zero entries in the Hessian of all oracles +- ncon: number of scalar constraints represented by all oracles +- lcon: lower bounds of oracle constraints +- ucon: upper bounds of oracle constraints +- nnzj: number of non-zero entries in the Jacobian of all oracles +- nnzh: number of non-zero entries in the Hessian of all oracles """ -mutable struct OraclesData +mutable struct Oracles oracles::Vector{Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}} - noracle::Int - l_oracle::Vector{Float64} - u_oracle::Vector{Float64} - nnzj_oracle::Int - nnzh_oracle::Int + ncon::Int + lcon::Vector{Float64} + ucon::Vector{Float64} + nnzj::Int + nnzh::Int end """ @@ -577,7 +580,6 @@ end Parse nonlinear constraints of an `nlp_data`. Returns: -- nnln: number of nonlinear constraints - nlcon: NonLinearStructure containing Jacobian and Hessian structures - nl_lcon: lower bounds of nonlinear constraints - nl_ucon: upper bounds of nonlinear constraints @@ -600,18 +602,19 @@ function parser_NL(nlp_data; hessian::Bool = true) nnzh = length(hess) nlcon = NonLinearStructure(nnln, jac_rows, jac_cols, nnzj, hess_rows, hess_cols, nnzh) - return nnln, nlcon, nl_lcon, nl_ucon + return nlcon, nl_lcon, nl_ucon end """ - parser_oracles(moimodel) + oracles, hessian_oracles_supported = parser_oracles(moimodel) Parse nonlinear oracles of a `MOI.ModelLike`. """ function parser_oracles(moimodel) + hessian_oracles_supported = true oracles = Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}[] - l_oracle = Float64[] - u_oracle = Float64[] + lcon = Float64[] + ucon = Float64[] # We know this pair exists from ListOfConstraintTypesPresent for ci in MOI.get( @@ -625,23 +628,26 @@ function parser_oracles(moimodel) push!(oracles, (f, cache)) # Bounds: MOI.VectorNonlinearOracle stores them internally (l, u) - append!(l_oracle, set.l) - append!(u_oracle, set.u) + append!(lcon, set.l) + append!(ucon, set.u) + + # Support for the Hessian + hessian_oracles_supported = hessian_oracles_supported && !isnothing(set.eval_hessian_lagrangian) end - # Sizes: number of scalar constraints represented by all oracles - noracle = length(l_oracle) + # Number of scalar constraints represented by all oracles + ncon = length(lcon) - # Sparsity: - nnzj_oracle = 0 - nnzh_oracle = 0 + # Number of nonzeros for the Jacobian and Hessian + nnzj = 0 + nnzh = 0 for (_, cache) in oracles - nnzj_oracle += length(cache.set.jacobian_structure) + nnzj += length(cache.set.jacobian_structure) # there may or may not be Hessian info - nnzh_oracle += length(cache.set.hessian_lagrangian_structure) + nnzh += length(cache.set.hessian_lagrangian_structure) end - return OraclesData(oracles, noracle, l_oracle, u_oracle, nnzj_oracle, nnzh_oracle) + return Oracles(oracles, ncon, lcon, ucon, nnzj, nnzh), hessian_oracles_supported end """ diff --git a/test/nlp_consistency.jl b/test/nlp_consistency.jl index 9d7d516..da9a406 100644 --- a/test/nlp_consistency.jl +++ b/test/nlp_consistency.jl @@ -1,15 +1,16 @@ for problem in [nlp_problems; extra_nlp_oracle_problems] @testset "Problem $problem" begin - nlp_manual = eval(Symbol(problem))() - problem_f = eval(Symbol(lowercase(problem))) + nlp_manual = eval(replace(problem, "_oracle" => "") |> uppercase |> Symbol)() + problem_f = eval(problem |> lowercase |> Symbol) nlp_moi = MathOptNLPModel(problem_f()) nlps = [nlp_manual; nlp_moi] + exclude = (problem in extra_nlp_oracle_problems) ? [jprod, jtprod, hprod] : [] consistent_nlps( nlps, linear_api = true, test_slack = false, test_counters = false, - exclude = [], + exclude = exclude, ) coord_memory_nlp(nlp_moi, linear_api = true) view_subarray_nlp(nlp_moi) diff --git a/test/nlp_problems/hs10.jl b/test/nlp_problems/hs10.jl index 5d29ad8..70ef30e 100644 --- a/test/nlp_problems/hs10.jl +++ b/test/nlp_problems/hs10.jl @@ -27,7 +27,7 @@ function hs10_oracle() # Bounds: 0 ≤ g(x) ≤ +∞ set = MOI.VectorNonlinearOracle(; dimension = 2, # number of input variables (x1, x2) - l = [-1.0], # lower bound on g(x) + l = [-1.0], # lower bound on g(x) u = [Inf], # upper bound on g(x) eval_f = (ret, xv) -> begin # ret[1] = g(x) diff --git a/test/test_moi_nlp_oracle.jl b/test/test_moi_nlp_oracle.jl index d945995..49853a2 100644 --- a/test/test_moi_nlp_oracle.jl +++ b/test/test_moi_nlp_oracle.jl @@ -1,10 +1,10 @@ -@testset "Testing Oracle vs Non-Oracle NLP Models: $prob" for prob in extra_nlp_oracle_problems - prob_no_oracle = replace(prob, "_oracle" => "") - prob_fn_no_oracle = eval(Symbol(prob_no_oracle)) - prob_fn = eval(Symbol(prob)) +@testset "Testing Oracle vs Non-Oracle NLP Models: $prob_oracle" for prob_oracle in extra_nlp_oracle_problems + prob_no_oracle = replace(prob_oracle, "_oracle" => "") + prob_fn_no_oracle = eval(prob_no_oracle |> Symbol) + prob_fn_oracle = eval(prob_oracle |> Symbol) - nlp_no_oracle = MathOptNLPModel(prob_fn_no_oracle(), hessian = true) - nlp_with_oracle = MathOptNLPModel(prob_fn(), hessian = true) + nlp_no_oracle = MathOptNLPModel(prob_fn_no_oracle(), hessian = true) + nlp_oracle = MathOptNLPModel(prob_fn_oracle(), hessian = true) n = nlp_no_oracle.meta.nvar m = nlp_no_oracle.meta.ncon @@ -12,23 +12,22 @@ # Objective value fx_no_oracle = obj(nlp_no_oracle, x) - fx_with_oracle = obj(nlp_with_oracle, x) + fx_with_oracle = obj(nlp_oracle, x) @test isapprox(fx_no_oracle, fx_with_oracle; atol = 1e-8, rtol = 1e-8) # Gradient of objective ngx_no_oracle = grad(nlp_no_oracle, x) - ngx_with_oracle = grad(nlp_with_oracle, x) + ngx_with_oracle = grad(nlp_oracle, x) @test isapprox(ngx_no_oracle, ngx_with_oracle; atol = 1e-8, rtol = 1e-8) # Constraint values (up to ordering) ncx_no_oracle = cons(nlp_no_oracle, x) - ncx_with_oracle = cons(nlp_with_oracle, x) - @test isapprox(sort(ncx_no_oracle), sort(ncx_with_oracle); - atol = 1e-8, rtol = 1e-8) + ncx_with_oracle = cons(nlp_oracle, x) + @test isapprox(sort(ncx_no_oracle), sort(ncx_with_oracle); atol = 1e-8, rtol = 1e-8) # Jacobian: compare J'J, which is invariant to row permutations J_no = jac(nlp_no_oracle, x) - J_with = jac(nlp_with_oracle, x) + J_with = jac(nlp_oracle, x) G_no = Matrix(J_no)' * Matrix(J_no) G_with = Matrix(J_with)' * Matrix(J_with) @@ -38,7 +37,7 @@ # Hessian of the objective: use y = 0 so constraints don't enter λ = zeros(m) H_no = hess(nlp_no_oracle, x, λ) - H_with = hess(nlp_with_oracle, x, λ) + H_with = hess(nlp_oracle, x, λ) @test isapprox(Matrix(H_no), Matrix(H_with); atol = 1e-8, rtol = 1e-8) -end \ No newline at end of file +end From 18ed91abb641fdd91921fdd35328de94c66a19a7 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Sat, 20 Dec 2025 09:50:55 +0100 Subject: [PATCH 16/19] Add support for jprod, jtprod and hprod --- src/moi_nlp_model.jl | 191 +++++++++++++++++++++++++++++++++------- src/utils.jl | 14 ++- test/nlp_consistency.jl | 19 ++-- test/nls_consistency.jl | 4 +- 4 files changed, 181 insertions(+), 47 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 872c574..5354216 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -36,7 +36,7 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = nlp_data = _nlp_block(moimodel) nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) - oracles, hessian_oracles_supported = parser_oracles(moimodel) + oracles = parser_oracles(moimodel) counters = Counters() λ = zeros(Float64, nlcon.nnln) # Lagrange multipliers for hess_coord! and hprod! without y hv = zeros(Float64, nvar) # workspace for ghjvprod! @@ -70,10 +70,8 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = minimize = MOI.get(moimodel, MOI.ObjectiveSense()) == MOI.MIN_SENSE, islp = (obj.type == "LINEAR") && (nlin == ncon), name = name, - jprod_available = isempty(oracles.oracles), - jtprod_available = isempty(oracles.oracles), - hprod_available = isempty(oracles.oracles) && hessian, - hess_available = hessian && hessian_oracles_supported, + hprod_available = hessian && oracles.hessian_oracles_supported, + hess_available = hessian && oracles.hessian_oracles_supported, ) return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, oracles, λ, hv, obj, counters), @@ -324,10 +322,10 @@ function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals: for i in 1:s.set.input_dimension s.x[i] = x[f.variables[i].value] end - nnz_oracle = length(s.set.jacobian_structure) - ind_oracle = (offset + 1):(offset + nnz_oracle) + nnzj_oracle = length(s.set.jacobian_structure) + ind_oracle = (offset + 1):(offset + nnzj_oracle) s.set.eval_jacobian(view(vals, ind_oracle), s.x) - offset += nnz_oracle + offset += nnzj_oracle end end @assert offset == nlp.meta.nln_nnzj @@ -375,10 +373,10 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs for i in 1:s.set.input_dimension s.x[i] = x[f.variables[i].value] end - nnz_oracle = length(s.set.jacobian_structure) - ind_oracle = (offset + 1):(offset + nnz_oracle) + nnzj_oracle = length(s.set.jacobian_structure) + ind_oracle = (offset + 1):(offset + nnzj_oracle) s.set.eval_jacobian(view(vals, ind_oracle), s.x) - offset += nnz_oracle + offset += nnzj_oracle end end @assert offset == nlp.meta.nnzj @@ -391,7 +389,6 @@ function NLPModels.jprod_lin!( v::AbstractVector, Jv::AbstractVector, ) - nlp.meta.jprod_available || error("J * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jprod_lin) jprod_lin!( nlp, @@ -410,7 +407,6 @@ function NLPModels.jprod_nln!( v::AbstractVector, Jv::AbstractVector, ) - nlp.meta.jprod_available || error("J * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jprod_nln) if nlp.quadcon.nquad > 0 for i = 1:(nlp.quadcon.nquad) @@ -420,9 +416,29 @@ function NLPModels.jprod_nln!( end end if nlp.nlcon.nnln > 0 - ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln) + ind_nnln = (nlp.quadcon.nquad + 1):(nlp.quadcon.nquad + nlp.nlcon.nnln) MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v) end + if nlp.oracles.ncon > 0 + for i in (nlp.quadcon.nquad + nlp.nlcon.nnln + 1):(nlp.quadcon.nquad + nlp.nlcon.nnln + nlp.oracles.ncon) + Jv[i] = 0 + end + row_offset = nlp.quadcon.nquad + nlp.nlcon.nnln + for (f, s) in nlp.oracles.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + s.set.eval_jacobian(s.nzJ, s.x) + k = 0 + for (r, c) in s.set.jacobian_structure + k += 1 + row = row_offset + r + col = f.variables[c].value + Jv[row] += s.nzJ[k] * v[col] + end + row_offset += s.set.output_dimension + end + end return Jv end @@ -432,7 +448,6 @@ function NLPModels.jprod!( v::AbstractVector, Jv::AbstractVector, ) - nlp.meta.jprod_available || error("J * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jprod) if nlp.meta.nlin > 0 view(Jv, nlp.meta.lin) .= 0.0 @@ -456,9 +471,29 @@ function NLPModels.jprod!( end end if nlp.nlcon.nnln > 0 - ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon) + ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln) MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v) end + if nlp.oracles.ncon > 0 + for i in (nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1):(nlp.meta.ncon) + Jv[i] = 0 + end + row_offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + for (f, s) in nlp.oracles.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + s.set.eval_jacobian(s.nzJ, s.x) + k = 0 + for (r, c) in s.set.jacobian_structure + k += 1 + row = row_offset + r + col = f.variables[c].value + Jv[row] += s.nzJ[k] * v[col] + end + row_offset += s.set.output_dimension + end + end return Jv end @@ -468,7 +503,6 @@ function NLPModels.jtprod_lin!( v::AbstractVector, Jtv::AbstractVector, ) - nlp.meta.jtprod_available || error("J' * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jtprod_lin) jtprod_lin!( nlp, @@ -487,13 +521,12 @@ function NLPModels.jtprod_nln!( v::AbstractVector, Jtv::AbstractVector, ) - nlp.meta.jtprod_available || error("J' * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jtprod_nln) if nlp.nlcon.nnln > 0 - ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln) + ind_nnln = (nlp.quadcon.nquad + 1):(nlp.quadcon.nquad + nlp.nlcon.nnln) MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln)) end - (nlp.quadcon.nquad == nlp.meta.nnln) && (Jtv .= 0.0) + (nlp.nlcon.nnln == 0) && (Jtv .= 0.0) if nlp.quadcon.nquad > 0 for i = 1:(nlp.quadcon.nquad) # Jtv += v[i] * (Aᵢ * x + bᵢ) @@ -502,6 +535,23 @@ function NLPModels.jtprod_nln!( Jtv .+= v[i] .* qcon.b end end + if nlp.oracles.ncon > 0 + row_offset = nlp.quadcon.nquad + nlp.nlcon.nnln + for (f, s) in nlp.oracles.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + s.set.eval_jacobian(s.nzJ, s.x) + k = 0 + for (r, c) in s.set.jacobian_structure + k += 1 + row = row_offset + r + col = f.variables[c].value + Jtv[col] += s.nzJ[k] * v[row] + end + row_offset += s.set.output_dimension + end + end return Jtv end @@ -511,13 +561,12 @@ function NLPModels.jtprod!( v::AbstractVector, Jtv::AbstractVector, ) - nlp.meta.jtprod_available || error("J' * v products are not supported by this MathOptNLPModel.") increment!(nlp, :neval_jtprod) if nlp.nlcon.nnln > 0 - ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon) + ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln) MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln)) end - (nlp.quadcon.nquad == nlp.meta.nnln) && (Jtv .= 0.0) + (nlp.nlcon.nnln == 0) && (Jtv .= 0.0) if nlp.meta.nlin > 0 transpose = true coo_unsym_add_mul!( @@ -537,6 +586,23 @@ function NLPModels.jtprod!( Jtv .+= v[nlp.meta.nlin + i] .* qcon.b end end + if nlp.oracles.ncon > 0 + row_offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + for (f, s) in nlp.oracles.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + s.set.eval_jacobian(s.nzJ, s.x) + k = 0 + for (r, c) in s.set.jacobian_structure + k += 1 + row = row_offset + r + col = f.variables[c].value + Jtv[col] += s.nzJ[k] * v[row] + end + row_offset += s.set.output_dimension + end + end return Jtv end @@ -545,6 +611,7 @@ function NLPModels.hess_structure!( rows::AbstractVector{<:Integer}, cols::AbstractVector{<:Integer}, ) + nlp.meta.hess_available || error("The function hess_structure! is not supported by this MathOptNLPModel.") if nlp.obj.type == "QUADRATIC" view(rows, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.rows view(cols, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.cols @@ -582,6 +649,7 @@ function NLPModels.hess_coord!( vals::AbstractVector; obj_weight::Float64 = 1.0, ) + nlp.meta.hess_available || error("The function hess_coord! is not supported by this MathOptNLPModel.") increment!(nlp, :neval_hess) # 1. Quadratic objective block (if any) @@ -650,6 +718,7 @@ function NLPModels.hess_coord!( vals::AbstractVector; obj_weight::Float64 = 1.0, ) + nlp.meta.hess_available || error("The function hess_coord! is not supported by this MathOptNLPModel.") increment!(nlp, :neval_hess) if nlp.obj.type == "LINEAR" vals .= 0.0 @@ -660,7 +729,8 @@ function NLPModels.hess_coord!( end if nlp.obj.type == "NONLINEAR" view(vals, 1:(nlp.obj.nnzh + nlp.quadcon.nnzh)) .= 0.0 - ind_nnln = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh) + view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh + 1):(nlp.meta.nnzh)) .= 0.0 + ind_nnln = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh) MOI.eval_hessian_lagrangian(nlp.eval, view(vals, ind_nnln), x, obj_weight, nlp.λ) end return vals @@ -672,6 +742,7 @@ function NLPModels.jth_hess_coord!( j::Integer, vals::AbstractVector, ) + nlp.meta.hess_available || error("The function jth_hess_coord! is not supported by this MathOptNLPModel.") increment!(nlp, :neval_jhess) @rangecheck 1 nlp.meta.ncon j vals .= 0.0 @@ -711,7 +782,7 @@ function NLPModels.jth_hess_coord!( if !find_oracle ncon_oracle = s.set.output_dimension nnzh_oracle = length(s.set.hessian_lagrangian_structure) - if offset_λ + 1 ≤ j ≤ offset_λ + if offset_λ + 1 ≤ j ≤ offset_λ + s.set.output_dimension find_oracle = true # build local x for i in 1:s.set.input_dimension @@ -720,7 +791,7 @@ function NLPModels.jth_hess_coord!( # We could preallocate the vector in the structure nlp λ_oracle = zeros(Float64, s.set.output_dimension) λ_oracle[j - offset_λ] = 1.0 - index_nnzh = (offset_nnzh + 1):(offset_nnzh + nnz_oracle) + index_nnzh = (offset_nnzh + 1):(offset_nnzh + nnzh_oracle) s.set.eval_hessian_lagrangian( view(vals, index_nnzh), s.x, @@ -744,17 +815,18 @@ function NLPModels.hprod!( hv::AbstractVector; obj_weight::Float64 = 1.0, ) - nlp.meta.hprod_available || error("H * v products are not supported by this MathOptNLPModel.") + nlp.meta.hprod_available || error("The function hprod! is not supported by this MathOptNLPModel.") increment!(nlp, :neval_hprod) if (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == 0) hv .= 0.0 end if (nlp.obj.type == "NONLINEAR") || (nlp.nlcon.nnln > 0) - λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)) + ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln) + λ = view(y, ind_nnln) MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, λ) end if nlp.obj.type == "QUADRATIC" - (nlp.meta.nnln == nlp.quadcon.nquad) && (hv .= 0.0) + (nlp.nlcon.nnln == 0) && (hv .= 0.0) coo_sym_add_mul!( nlp.obj.hessian.rows, nlp.obj.hessian.cols, @@ -765,12 +837,32 @@ function NLPModels.hprod!( ) end if nlp.quadcon.nquad > 0 - (nlp.obj.type == "LINEAR") && (hv .= 0.0) + (nlp.obj.type == "LINEAR") && (nlp.nlcon.nnln == 0) && (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 + if nlp.oracles.ncon > 0 + (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == nlp.oracles.ncon) && (hv .= 0.0) + offset_y = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + for (f, s) in nlp.oracles.oracles + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + ind_y = (offset_y + 1):(offset_y + s.set.output_dimension) + y_oracle = view(y, ind_y) + s.set.eval_hessian_lagrangian(s.nzH, s.x, y_oracle) + k = 0 + for (r, c) in s.set.hessian_lagrangian_structure + k += 1 + row = f.variables[r].value + col = f.variables[c].value + hv[row] += s.nzH[k] * v[col] + end + offset_y += s.set.output_dimension + end + end return hv end @@ -781,7 +873,7 @@ function NLPModels.hprod!( hv::AbstractVector; obj_weight::Float64 = 1.0, ) - nlp.meta.hprod_available || error("H * v products are not supported by this MathOptNLPModel.") + nlp.meta.hprod_available || error("The function hprod! is not supported by this MathOptNLPModel.") increment!(nlp, :neval_hprod) if nlp.obj.type == "LINEAR" hv .= 0.0 @@ -810,17 +902,48 @@ function NLPModels.jth_hprod!( j::Integer, hv::AbstractVector, ) - nlp.meta.hprod_available || error("H * v products are not supported by this MathOptNLPModel.") + nlp.meta.hprod_available || error("The function jth_hprod! is not supported by this MathOptNLPModel.") increment!(nlp, :neval_jhprod) @rangecheck 1 nlp.meta.ncon j hv .= 0.0 if nlp.meta.nlin + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad qcon = nlp.quadcon.constraints[j - nlp.meta.nlin] coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, 1.0) - elseif nlp.meta.nlin + nlp.quadcon.nquad + 1 ≤ j ≤ nlp.meta.ncon + elseif nlp.meta.nlin + nlp.quadcon.nquad + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 1.0 MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, 0.0, nlp.λ) nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 0.0 + elseif nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1 ≤ j ≤ nlp.meta.ncon + find_oracle = false + offset_λ = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + for (f, s) in nlp.oracles.oracles + if !find_oracle + if offset_λ + 1 ≤ j ≤ offset_λ + s.set.output_dimension + find_oracle = true + # build local x + for i in 1:s.set.input_dimension + s.x[i] = x[f.variables[i].value] + end + # We could preallocate the vector in the structure nlp + λ_oracle = zeros(Float64, s.set.output_dimension) + λ_oracle[j - offset_λ] = 1.0 + s.set.eval_hessian_lagrangian( + s.nzH, + s.x, + λ_oracle, + ) + k = 0 + for (r, c) in s.set.hessian_lagrangian_structure + k += 1 + row = f.variables[r].value + col = f.variables[c].value + hv[row] += s.nzH[k] * v[col] + end + else + offset_λ += ncon_oracle + end + end + end end return hv end @@ -832,7 +955,7 @@ function NLPModels.ghjvprod!( v::AbstractVector, ghv::AbstractVector, ) - nlp.meta.hprod_available || error("The operation is not supported by this MathOptNLPModel.") + nlp.meta.hprod_available || error("The function ghjvprod! is not supported by this MathOptNLPModel.") increment!(nlp, :neval_hprod) ghv .= 0.0 for i = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad) diff --git a/src/utils.jl b/src/utils.jl index acc8dd0..3a8a984 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -31,9 +31,13 @@ const ORACLE = MOI.VectorNonlinearOracle{Float64} # VectorNonlinearOracle{Float mutable struct _VectorNonlinearOracleCache set::MOI.VectorNonlinearOracle{Float64} x::Vector{Float64} + nzJ::Vector{Float64} + nzH::Vector{Float64} function _VectorNonlinearOracleCache(set::MOI.VectorNonlinearOracle{Float64}) - return new(set, zeros(set.input_dimension)) + nnzj = length(set.jacobian_structure) + nnzh = length(set.hessian_lagrangian_structure) + return new(set, zeros(set.input_dimension), zeros(nnzj), zeros(nnzh)) end end @@ -133,6 +137,9 @@ Structure containing nonlinear oracles data: - ucon: upper bounds of oracle constraints - nnzj: number of non-zero entries in the Jacobian of all oracles - nnzh: number of non-zero entries in the Hessian of all oracles +- nzJ: buffer to store the nonzeros of the Jacobian for all oracles (needed for the functions jprod and jtprod) +- nzH: buffer to store the nonzeros of the Hessian for all oracles (needed for the function hprod) +- hessian_oracles_supported: support of the Hessian for all oracles """ mutable struct Oracles oracles::Vector{Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}} @@ -141,6 +148,7 @@ mutable struct Oracles ucon::Vector{Float64} nnzj::Int nnzh::Int + hessian_oracles_supported::Bool end """ @@ -606,7 +614,7 @@ function parser_NL(nlp_data; hessian::Bool = true) end """ - oracles, hessian_oracles_supported = parser_oracles(moimodel) + oracles = parser_oracles(moimodel) Parse nonlinear oracles of a `MOI.ModelLike`. """ @@ -647,7 +655,7 @@ function parser_oracles(moimodel) nnzh += length(cache.set.hessian_lagrangian_structure) end - return Oracles(oracles, ncon, lcon, ucon, nnzj, nnzh), hessian_oracles_supported + return Oracles(oracles, ncon, lcon, ucon, nnzj, nnzh, hessian_oracles_supported) end """ diff --git a/test/nlp_consistency.jl b/test/nlp_consistency.jl index da9a406..8c66b1e 100644 --- a/test/nlp_consistency.jl +++ b/test/nlp_consistency.jl @@ -4,15 +4,16 @@ for problem in [nlp_problems; extra_nlp_oracle_problems] problem_f = eval(problem |> lowercase |> Symbol) nlp_moi = MathOptNLPModel(problem_f()) nlps = [nlp_manual; nlp_moi] - exclude = (problem in extra_nlp_oracle_problems) ? [jprod, jtprod, hprod] : [] - consistent_nlps( - nlps, - linear_api = true, - test_slack = false, - test_counters = false, - exclude = exclude, - ) - coord_memory_nlp(nlp_moi, linear_api = true) + @testset "linear_api = $linear_api" for linear_api in (false, true) + consistent_nlps( + nlps, + linear_api = linear_api, + test_slack = false, + test_counters = false, + exclude = [], + ) + coord_memory_nlp(nlp_moi, linear_api = linear_api) + end view_subarray_nlp(nlp_moi) end end diff --git a/test/nls_consistency.jl b/test/nls_consistency.jl index 38a2159..353d856 100644 --- a/test/nls_consistency.jl +++ b/test/nls_consistency.jl @@ -8,7 +8,9 @@ for problem in nls_problems if isdefined(Main, Symbol(spc)) push!(nlss, eval(Meta.parse(spc))()) end - consistent_nlss(nlss, linear_api = true, test_slack = false, test_counters = false) + @testset "linear_api = $linear_api" for linear_api in (false, true) + consistent_nlss(nlss, linear_api = linear_api, test_slack = false, test_counters = false) + end view_subarray_nls(nls_moi) end end From 8614e026755a52b83cbd4a98c5bd18457e606fcc Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Mon, 22 Dec 2025 01:47:26 +0100 Subject: [PATCH 17/19] Fix the last bugs! --- src/moi_nlp_model.jl | 6 ++++++ test/nlp_consistency.jl | 2 +- test/nlp_problems/hs10.jl | 4 ++-- test/nlp_problems/hs14.jl | 6 +++--- test/runtests.jl | 1 + 5 files changed, 13 insertions(+), 6 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 5354216..d1e7e7f 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -859,6 +859,9 @@ function NLPModels.hprod!( row = f.variables[r].value col = f.variables[c].value hv[row] += s.nzH[k] * v[col] + if row != col + hv[col] += s.nzH[k] * v[row] + end end offset_y += s.set.output_dimension end @@ -938,6 +941,9 @@ function NLPModels.jth_hprod!( row = f.variables[r].value col = f.variables[c].value hv[row] += s.nzH[k] * v[col] + if row != col + hv[col] += s.nzH[k] * v[row] + end end else offset_λ += ncon_oracle diff --git a/test/nlp_consistency.jl b/test/nlp_consistency.jl index 8c66b1e..e3739d0 100644 --- a/test/nlp_consistency.jl +++ b/test/nlp_consistency.jl @@ -1,4 +1,4 @@ -for problem in [nlp_problems; extra_nlp_oracle_problems] +for problem in [nlp_problems; nlp_oracle_problems] @testset "Problem $problem" begin nlp_manual = eval(replace(problem, "_oracle" => "") |> uppercase |> Symbol)() problem_f = eval(problem |> lowercase |> Symbol) diff --git a/test/nlp_problems/hs10.jl b/test/nlp_problems/hs10.jl index 70ef30e..200e8e4 100644 --- a/test/nlp_problems/hs10.jl +++ b/test/nlp_problems/hs10.jl @@ -45,11 +45,11 @@ function hs10_oracle() # ∂²g/∂x1² = -6 # ∂²g/∂x1∂x2 = ∂²g/∂x2∂x1 = 2 # ∂²g/∂x2² = -2 - hessian_lagrangian_structure = [(1, 1), (1, 2), (2, 2)], + hessian_lagrangian_structure = [(1, 1), (2, 1), (2, 2)], eval_hessian_lagrangian = (ret, xv, μ) -> begin # Hessian of μ[1] * g(x) ret[1] = μ[1] * (-6.0) # (1,1) - ret[2] = μ[1] * ( 2.0) # (1,2) + ret[2] = μ[1] * ( 2.0) # (2,1) ret[3] = μ[1] * (-2.0) # (2,2) end, ) diff --git a/test/nlp_problems/hs14.jl b/test/nlp_problems/hs14.jl index 19b6588..68f78d2 100644 --- a/test/nlp_problems/hs14.jl +++ b/test/nlp_problems/hs14.jl @@ -48,12 +48,12 @@ function hs14_oracle() # ∂²g/∂x1∂x2 = 0 # ∂²g/∂x2² = -2 # - # We store only upper-triangular entries: (1,1), (1,2), (2,2) - hessian_lagrangian_structure = [(1, 1), (1, 2), (2, 2)], + # We store only lower-triangular entries: (1,1), (2,1), (2,2) + hessian_lagrangian_structure = [(1, 1), (2, 1), (2, 2)], eval_hessian_lagrangian = (ret, xv, μ) -> begin # Hessian of μ[1] * g(x) ret[1] = μ[1] * (-0.5) # (1,1) - ret[2] = μ[1] * 0.0 # (1,2) + ret[2] = μ[1] * 0.0 # (2,1) ret[3] = μ[1] * (-2.0) # (2,2) end, ) diff --git a/test/runtests.jl b/test/runtests.jl index 3f09d96..7588866 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ nls_problems = NLPModelsTest.nls_problems extra_nlp_problems = ["nohesspb", "hs61", "hs100", "hs219", "quadcon", "operatorspb", "nf"] extra_nls_problems = ["nlsnohesspb", "HS30", "HS43", "MGH07", "nlsqc"] +nlp_oracle_problems = ["hs10_oracle", "hs14_oracle"] extra_nlp_oracle_problems = ["hs10_oracle", "hs14_oracle", "hs61_oracle", "hs100_oracle"] for problem in lowercase.(nlp_problems ∪ extra_nlp_problems) From da200656110c5289abe0bf458719531b078638fc Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Mon, 22 Dec 2025 02:44:20 +0100 Subject: [PATCH 18/19] Use the last release of NLPModelsTest.jl --- Project.toml | 4 ++-- src/moi_nlp_model.jl | 6 +++--- src/moi_nls_model.jl | 6 +++--- src/utils.jl | 12 +++++++----- 4 files changed, 15 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index 84a70c6..5566c2d 100644 --- a/Project.toml +++ b/Project.toml @@ -16,10 +16,10 @@ JuMP = "1.25" LinearAlgebra = "1.10" MathOptInterface = "1.46" NLPModels = "0.21.6" -NLPModelsTest = "0.10.5" +NLPModelsTest = "0.10.6" Percival = "0.7.3" Printf = "1.10" -SolverCore = "0.3" +SolverCore = "0.3.9" SparseArrays = "1.10" Test = "1.10" julia = "1.10" diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index d1e7e7f..0a557fb 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -35,7 +35,7 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = parser_MOI(moimodel, index_map, nvar) nlp_data = _nlp_block(moimodel) - nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) + nlcon = parser_NL(nlp_data, hessian = hessian) oracles = parser_oracles(moimodel) counters = Counters() λ = zeros(Float64, nlcon.nnln) # Lagrange multipliers for hess_coord! and hprod! without y @@ -49,8 +49,8 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = # Total counts ncon = nlin + quadcon.nquad + nlcon.nnln + oracles.ncon - lcon = vcat(lin_lcon, quad_lcon, nl_lcon, oracles.lcon) - ucon = vcat(lin_ucon, quad_ucon, nl_ucon, oracles.ucon) + lcon = vcat(lin_lcon, quad_lcon, nlcon.nl_lcon, oracles.lcon) + ucon = vcat(lin_ucon, quad_ucon, nlcon.nl_ucon, oracles.ucon) nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj + oracles.nnzj nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh + oracles.nnzh meta = NLPModelMeta( diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl index a1846dd..e70d701 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 parser_MOI(moimodel, index_map, nvar) nlp_data = _nlp_block(moimodel) - nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian) + nlcon = parser_NL(nlp_data, hessian = hessian) λ = zeros(nlcon.nnln) # Lagrange multipliers for hess_coord! and hprod! without y nequ = nlinequ + nnlnequ @@ -44,8 +44,8 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri Fnnzh = nlequ.nnzh ncon = nlin + quadcon.nquad + nlcon.nnln - lcon = vcat(lin_lcon, quad_lcon, nl_lcon) - ucon = vcat(lin_ucon, quad_ucon, nl_ucon) + lcon = vcat(lin_lcon, quad_lcon, nlcon.nl_lcon) + ucon = vcat(lin_ucon, quad_ucon, nlcon.nl_ucon) cnnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj cnnzh = lls.nnzh + quadcon.nnzh + nlcon.nnzh diff --git a/src/utils.jl b/src/utils.jl index 3a8a984..5ff1ef8 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -96,6 +96,8 @@ end Structure containing Jacobian and Hessian structures of nonlinear constraints: - nnln: number of nonlinear constraints +- nl_lcon: lower bounds of nonlinear constraints +- nl_ucon: upper bounds of nonlinear constraints - jac_rows: row indices of the Jacobian in Coordinate format (COO) format - jac_cols: column indices of the Jacobian in COO format - nnzj: number of non-zero entries in the Jacobian @@ -105,6 +107,8 @@ Structure containing Jacobian and Hessian structures of nonlinear constraints: """ mutable struct NonLinearStructure nnln::Int + nl_lcon::Vector{Float64} + nl_ucon::Vector{Float64} jac_rows::Vector{Int} jac_cols::Vector{Int} nnzj::Int @@ -589,8 +593,6 @@ Parse nonlinear constraints of an `nlp_data`. Returns: - nlcon: NonLinearStructure containing Jacobian and Hessian structures -- nl_lcon: lower bounds of nonlinear constraints -- nl_ucon: upper bounds of nonlinear constraints """ function parser_NL(nlp_data; hessian::Bool = true) nnln = length(nlp_data.constraint_bounds) @@ -608,9 +610,9 @@ function parser_NL(nlp_data; hessian::Bool = true) hess_rows = hessian ? getindex.(hess, 1) : Int[] hess_cols = hessian ? getindex.(hess, 2) : Int[] nnzh = length(hess) - nlcon = NonLinearStructure(nnln, jac_rows, jac_cols, nnzj, hess_rows, hess_cols, nnzh) + nlcon = NonLinearStructure(nnln, nl_lcon, nl_ucon, jac_rows, jac_cols, nnzj, hess_rows, hess_cols, nnzh) - return nlcon, nl_lcon, nl_ucon + return nlcon end """ @@ -886,7 +888,7 @@ function parser_nonlinear_expression(cmodel, nvar, F; hessian::Bool = true) Fhess_cols = hessian && Feval ≠ nothing ? getindex.(Fhess, 2) : Int[] nl_Fnnzh = length(Fhess) - nlequ = NonLinearStructure(nnlnequ, Fjac_rows, Fjac_cols, nl_Fnnzj, Fhess_rows, Fhess_cols, nl_Fnnzh) + nlequ = NonLinearStructure(nnlnequ, Float64[], Float64[], Fjac_rows, Fjac_cols, nl_Fnnzj, Fhess_rows, Fhess_cols, nl_Fnnzh) return Feval, nlequ, nnlnequ end From 9b8c8c548bb566c27f260d6de8e8f23041a0253b Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Mon, 22 Dec 2025 03:08:15 +0100 Subject: [PATCH 19/19] Fix an error in cons_nln! --- src/moi_nlp_model.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl index 0a557fb..ed53b66 100644 --- a/src/moi_nlp_model.jl +++ b/src/moi_nlp_model.jl @@ -129,7 +129,7 @@ function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::Abstrac if nlp.nlcon.nnln > 0 offset += nlp.nlcon.nnln index_nlcon = (nlp.quadcon.nquad + 1):(offset) - MOI.eval_constraint(nlp.eval, view(c, ), x) + MOI.eval_constraint(nlp.eval, view(c, index_nlcon), x) end if nlp.oracles.ncon > 0 for (f, s) in nlp.oracles.oracles @@ -162,7 +162,8 @@ function NLPModels.cons!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVec end if nlp.nlcon.nnln > 0 offset += nlp.nlcon.nnln - MOI.eval_constraint(nlp.eval, view(c, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(offset)), x) + index_nlcon = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(offset) + MOI.eval_constraint(nlp.eval, view(c, index_nlcon), x) end if nlp.oracles.ncon > 0 for (f, s) in nlp.oracles.oracles