Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support quadratic constraints #181

Merged
merged 21 commits into from
Jun 21, 2024
Merged
Show file tree
Hide file tree
Changes from 19 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
162 changes: 126 additions & 36 deletions src/moi_nlp_model.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,9 @@ mutable struct MathOptNLPModel <: AbstractNLPModel{Float64, Vector{Float64}}
meta::NLPModelMeta{Float64, Vector{Float64}}
eval::MOI.Nonlinear.Evaluator
lincon::LinearConstraints
quadcon::QuadraticConstraints
nlcon::NonLinearStructure
λ::Vector{Float64}
obj::Objective
counters::Counters
end
Expand All @@ -27,22 +29,23 @@ end

function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = "Generic")
index_map, nvar, lvar, uvar, x0 = parser_variables(moimodel)
nlin, lincon, lin_lcon, lin_ucon = parser_MOI(moimodel, index_map)
nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon = parser_MOI(moimodel, index_map, nvar)

nlp_data = _nlp_block(moimodel)
nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian)
λ = zeros(nnln) # Lagrange multipliers for hess_coord! and hprod! without y

if nlp_data.has_objective
obj = Objective("NONLINEAR", 0.0, spzeros(Float64, nvar), COO(), 0)
else
obj = parser_objective_MOI(moimodel, nvar, index_map)
end

ncon = nlin + nnln
lcon = vcat(lin_lcon, nl_lcon)
ucon = vcat(lin_ucon, nl_ucon)
nnzj = lincon.nnzj + nlcon.nnzj
nnzh = obj.nnzh + nlcon.nnzh
ncon = nlin + quadcon.nquad + nnln
lcon = vcat(lin_lcon, quad_lcon, nl_lcon)
ucon = vcat(lin_ucon, quad_ucon, nl_ucon)
nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj
nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh

meta = NLPModelMeta(
nvar,
Expand All @@ -57,13 +60,13 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String =
nnzh = nnzh,
lin = collect(1:nlin),
lin_nnzj = lincon.nnzj,
nln_nnzj = nlcon.nnzj,
nln_nnzj = quadcon.nnzj + nlcon.nnzj,
minimize = MOI.get(moimodel, MOI.ObjectiveSense()) == MOI.MIN_SENSE,
islp = (obj.type == "LINEAR") && (nnln == 0),
islp = (obj.type == "LINEAR") && (nnln == 0) && (quadcon.nquad == 0),
name = name,
)

return MathOptNLPModel(meta, nlp_data.evaluator, lincon, nlcon, obj, Counters()), index_map
return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, λ, obj, Counters()), index_map
end

function NLPModels.obj(nlp::MathOptNLPModel, x::AbstractVector)
Expand Down Expand Up @@ -106,7 +109,13 @@ end

function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVector)
increment!(nlp, :neval_cons_nln)
MOI.eval_constraint(nlp.eval, c, x)
for i = 1:(nlp.quadcon.nquad)
qcon = nlp.quadcon.constraints[i]
c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
end
if nlp.meta.nnln > nlp.quadcon.nquad
MOI.eval_constraint(nlp.eval, view(c, (nlp.quadcon.nquad + 1):(nlp.meta.nnln)), x)
end
return c
end

Expand All @@ -125,8 +134,20 @@ function NLPModels.jac_nln_structure!(
rows::AbstractVector{<:Integer},
cols::AbstractVector{<:Integer},
)
view(rows, 1:(nlp.nlcon.nnzj)) .= nlp.nlcon.jac_rows
view(cols, 1:(nlp.nlcon.nnzj)) .= nlp.nlcon.jac_cols
if nlp.quadcon.nquad > 0
index = 0
for i = 1:(nlp.quadcon.nquad)
qcon = nlp.quadcon.constraints[i]
view(rows, index+1:index+qcon.nnzg) .= i
view(cols, index+1:index+qcon.nnzg) .= qcon.g
index += qcon.nnzg
end
end
if nlp.meta.nnln > nlp.quadcon.nquad
ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj)
view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nlp.nlcon.jac_rows
view(cols, ind_nnln) .= nlp.nlcon.jac_cols
end
return rows, cols
end

Expand All @@ -138,7 +159,33 @@ end

function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::AbstractVector)
increment!(nlp, :neval_jac_nln)
MOI.eval_constraint_jacobian(nlp.eval, view(vals, 1:(nlp.nlcon.nnzj)), x)
if nlp.quadcon.nquad > 0
index = 0
view(vals, 1:nlp.quadcon.nnzj) .= 0.0
for i = 1:(nlp.quadcon.nquad)
qcon = nlp.quadcon.constraints[i]
for (j, ind) in enumerate(qcon.b.nzind)
k = qcon.dg[ind]
vals[index+k] += qcon.b.nzval[j]
end
for j = 1:qcon.nnzh
row = qcon.A.rows[j]
col = qcon.A.cols[j]
val = qcon.A.vals[j]
k1 = qcon.dg[row]
vals[index+k1] += val * x[col]
if row != col
k2 = qcon.dg[col]
vals[index+k2] += val * x[row]
end
end
index += qcon.nnzg
end
end
if nlp.meta.nnln > nlp.quadcon.nquad
ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj)
MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x)
end
return vals
end

Expand Down Expand Up @@ -166,7 +213,18 @@ function NLPModels.jprod_nln!(
Jv::AbstractVector,
)
increment!(nlp, :neval_jprod_nln)
MOI.eval_constraint_jacobian_product(nlp.eval, Jv, x, v)
if nlp.meta.nnln > nlp.quadcon.nquad
ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln)
MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v)
end
(nlp.meta.nnln == nlp.quadcon.nquad) && (Jv .= 0.0)
if nlp.quadcon.nquad > 0
for i = 1:(nlp.quadcon.nquad)
# Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v
qcon = nlp.quadcon.constraints[i]
Jv[i] += coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
end
end
return Jv
end

Expand Down Expand Up @@ -194,7 +252,19 @@ function NLPModels.jtprod_nln!(
Jtv::AbstractVector,
)
increment!(nlp, :neval_jtprod_nln)
MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, v)
if nlp.meta.nnln > nlp.quadcon.nquad
ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln)
MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln))
end
(nlp.meta.nnln == nlp.quadcon.nquad) && (Jtv .= 0.0)
if nlp.quadcon.nquad > 0
for i = 1:(nlp.quadcon.nquad)
# Jtv += v[i] * (Aᵢ * x + bᵢ)
qcon = nlp.quadcon.constraints[i]
coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, Jtv, v[i])
Jtv .+= v[i] .* qcon.b
end
end
return Jtv
end

Expand All @@ -207,9 +277,18 @@ function NLPModels.hess_structure!(
view(rows, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.rows
view(cols, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.cols
end
if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
view(rows, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_rows
view(cols, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_cols
if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
view(rows, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_rows
view(cols, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_cols
end
if nlp.quadcon.nquad > 0
index = nlp.obj.nnzh
for i = 1:(nlp.quadcon.nquad)
qcon = nlp.quadcon.constraints[i]
view(rows, (index + 1):(index + qcon.nnzh)) .= qcon.A.rows
view(cols, (index + 1):(index + qcon.nnzh)) .= qcon.A.cols
index += qcon.nnzh
end
end
return rows, cols
end
Expand All @@ -225,15 +304,23 @@ function NLPModels.hess_coord!(
if nlp.obj.type == "QUADRATIC"
view(vals, 1:(nlp.obj.nnzh)) .= obj_weight .* nlp.obj.hessian.vals
end
if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
MOI.eval_hessian_lagrangian(
nlp.eval,
view(vals, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)),
if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon))
MOI.eval_hessian_lagrangian(nlp.eval,
view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)),
x,
obj_weight,
view(y, nlp.meta.nln),
λ
)
end
if nlp.quadcon.nquad > 0
index = nlp.obj.nnzh
for i = 1:(nlp.quadcon.nquad)
qcon = nlp.quadcon.constraints[i]
view(vals, (index + 1):(index + qcon.nnzh)) .= y[nlp.meta.nlin + i] .* qcon.A.vals
index += qcon.nnzh
end
end
return vals
end

Expand All @@ -252,7 +339,9 @@ function NLPModels.hess_coord!(
view(vals, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= 0.0
end
if nlp.obj.type == "NONLINEAR"
MOI.eval_hessian_lagrangian(nlp.eval, vals, x, obj_weight, zeros(nlp.meta.nnln))
view(vals, 1:(nlp.obj.nnzh + nlp.quadcon.nnzh)) .= 0.0
ind_nnln = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)
MOI.eval_hessian_lagrangian(nlp.eval, view(vals, ind_nnln), x, obj_weight, nlp.λ)
end
return vals
end
Expand All @@ -269,19 +358,20 @@ function NLPModels.hprod!(
if (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == 0)
hv .= 0.0
end
if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, view(y, nlp.meta.nln))
if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon))
MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, λ)
end
if nlp.obj.type == "QUADRATIC"
nlp.meta.nnln == 0 && (hv .= 0.0)
coo_sym_add_mul!(
nlp.obj.hessian.rows,
nlp.obj.hessian.cols,
nlp.obj.hessian.vals,
v,
hv,
obj_weight,
)
(nlp.meta.nnln == nlp.quadcon.nquad) && (hv .= 0.0)
coo_sym_add_mul!(nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals, v, hv, obj_weight)
end
if nlp.quadcon.nquad > 0
(nlp.obj.type == "LINEAR") && (hv .= 0.0)
for i = 1:(nlp.quadcon.nquad)
qcon = nlp.quadcon.constraints[i]
coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, y[nlp.meta.nlin + i])
end
end
return hv
end
Expand All @@ -302,7 +392,7 @@ function NLPModels.hprod!(
coo_sym_add_mul!(nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals, v, hv, obj_weight)
end
if nlp.obj.type == "NONLINEAR"
MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, zeros(nlp.meta.nnln))
MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, nlp.λ)
end
return hv
end
Loading
Loading