Support quadratic constraints (#181)

Author: amontoison

src/moi_nlp_model.jl

@@ -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
@@ -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,
@@ -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)
@@ -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
 
@@ -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
 
@@ -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
 
@@ -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
 
@@ -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
 
@@ -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
@@ -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
 
@@ -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
@@ -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
@@ -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