Modify QuadraticModels to NLPModels 0.12

Author: abelsiqueira

src/qpmodel.jl

@@ -1,143 +1,182 @@
 export jac_structure!, hess_structure!, jac_coord!, hess_coord!
 
 mutable struct QPData
-    c0  :: Float64                  # constant term in objective
-    c   :: AbstractVector{Float64}  # linear term
-    H   :: AbstractMatrix{Float64}  # quadratic term
-    opH :: AbstractLinearOperator   # assumed with preallocation!
-    A   :: AbstractMatrix{Float64}  # constraint matrix
+  c0    :: Float64          # constant term in objective
+  c     :: Vector           # linear term
+  Hrows :: Vector{Int}      # quadratic term
+  Hcols :: Vector{Int}
+  Hvals :: Vector
+  Arows :: Vector{Int}      # constraints matrix
+  Acols :: Vector{Int}
+  Avals :: Vector
 end
 
 abstract type AbstractQuadraticModel <: AbstractNLPModel end
 
 mutable struct QuadraticModel <: AbstractQuadraticModel
-    meta     :: NLPModelMeta
-    counters :: Counters
-    data     :: QPData
-
-  function QuadraticModel(c :: AbstractVector{Float64}, H :: AbstractMatrix{Float64},
-                          opH :: AbstractLinearOperator,
-                          A :: AbstractMatrix{Float64},
-                          lcon :: AbstractVector{Float64}, ucon :: AbstractVector{Float64},
-                          lvar :: AbstractVector{Float64}, uvar :: AbstractVector{Float64};
-                          c0 :: Float64=0.0, kwargs...)
-    ncon, nvar = size(A)
-    nnzh = issparse(H) ? nnz(H) : (nvar * (nvar + 1) / 2)
-    nnzj = issparse(A) ? nnz(A) : (nvar * ncon)
-    new(NLPModelMeta(nvar,
-                     lvar=lvar, uvar=uvar,
-                     ncon=size(A,1), lcon=lcon, ucon=ucon,
-                     nnzj=nnzj,
-                     nnzh=nnzh,
-                     lin=1:ncon, nln=Int[], islp=(ncon == 0); kwargs...),
-        Counters(),
-        QPData(c0, c, H, opH, A))
+  meta     :: NLPModelMeta
+  counters :: Counters
+  data     :: QPData
+end
+
+function QuadraticModel(c :: AbstractVector,
+                        Hrows :: AbstractVector{<: Integer}, Hcols :: AbstractVector{<: Integer}, Hvals :: AbstractVector;
+                        Arows :: AbstractVector{<: Integer} = Int[], Acols :: AbstractVector{<: Integer} = Int[], Avals :: AbstractVector = Float64[],
+                        lcon :: AbstractVector = Float64[], ucon :: AbstractVector = Float64[],
+                        lvar :: AbstractVector = fill(-Inf, length(c)), uvar :: AbstractVector = fill(Int, length(c)),
+                        c0 :: Float64=0.0, kwargs...)
+  nnzh = length(Hvals)
+  if !(nnzh == length(Hrows) == length(Hcols))
+    error("The length of Hrows, Hcols and Hvals must be the same")
+  end
+  nnzj = length(Avals)
+  if !(nnzj == length(Arows) == length(Acols))
+    error("The length of Arows, Acols and Avals must be the same")
+  end
+  ncon = length(lcon)
+  if ncon != length(ucon)
+    error("The length of lcon and ucon must be the same")
   end
+  nvar = length(c)
+  if !(nvar == length(lvar) == length(uvar))
+    error("The length of c, lvar and uvar must be the same")
+  end
+  QuadraticModel(NLPModelMeta(length(c), lvar=lvar, uvar=uvar,
+                              ncon=ncon, lcon=lcon, ucon=ucon,
+                              nnzj=nnzj, nnzh=nnzh,
+                              lin=1:ncon, nln=Int[], islp=(ncon == 0); kwargs...),
+                 Counters(),
+                 QPData(c0, c, Hrows, Hcols, Hvals, Arows, Acols, Avals))
 end
 
-function QuadraticModel(model :: AbstractNLPModel, x :: AbstractVector = model.meta.x0, args...; kwargs...)
-    nvar = model.meta.nvar
-    ncon = model.meta.ncon
-    g = grad(model, x)
-    H = hess(model, x, args...; kwargs...)
-    c = cons(model, x)
-    A = jac(model, x)
-    QuadraticModel(g, H, opHermitian(H), A,
-                   model.meta.lcon .- c, model.meta.ucon .- c,
-                   model.meta.lvar .- x, model.meta.uvar .- x)
+function QuadraticModel(c :: AbstractVector, H :: AbstractMatrix;
+                        A :: AbstractMatrix=zeros(0,length(c)), lcon :: AbstractVector=zeros(0), ucon :: AbstractVector=zeros(0),
+                        lvar :: AbstractVector = fill(-Inf, length(c)), uvar :: AbstractVector = fill(Inf, length(c)),
+                        c0 :: Float64=0.0, kwargs...)
+  ncon, nvar = size(A)
+  nnzh, Hrows, Hcols, Hvals = if issparse(H)
+    nnz(tril(H)), findnz(tril(H))...
+  else
+    I = ((i,j,H[i,j]) for i = 1:nvar, j = 1:nvar if i ≥ j)
+    div(nvar * (nvar + 1), 2), getindex.(I, 1), getindex.(I, 2), getindex.(I, 3)
+  end
+  nnzj, Arows, Acols, Avals = if issparse(A)
+    nnz(A), findnz(A)...
+  else
+    I = ((i,j,A[i,j]) for i = 1:ncon, j = 1:nvar)
+    nvar * ncon, getindex.(I, 1)[:], getindex.(I, 2)[:], getindex.(I, 3)[:]
+  end
+  QuadraticModel(NLPModelMeta(nvar, lvar=lvar, uvar=uvar,
+                              ncon=size(A,1), lcon=lcon, ucon=ucon,
+                              nnzj=nnzj, nnzh=nnzh,
+                              lin=1:ncon, nln=Int[], islp=(ncon == 0); kwargs...),
+                 Counters(),
+                 QPData(c0, c, Hrows, Hcols, Hvals, Arows, Acols, Avals))
 end
 
-linobj(qp::AbstractQuadraticModel, args...) = qp.data.c
+"""
+    QuadraticModel(nlp, x)
 
-function NLPModels.objgrad(qp :: AbstractQuadraticModel, x :: AbstractVector)
-    g = Vector{eltype(x)}(undef, length(x))
-    objgrad!(qp, x, g)
+Creates a quadratic Taylor model of `nlp` around `x`.
+"""
+function QuadraticModel(model :: AbstractNLPModel, x :: AbstractVector; kwargs...)
+  nvar = model.meta.nvar
+  ncon = model.meta.ncon
+  c0 = obj(model, x)
+  g = grad(model, x)
+  Hrows, Hcols = hess_structure(model)
+  Hvals = hess_coord(model, x)
+  if model.meta.ncon > 0
+    c = cons(model, x)
+    Arows, Acols = jac_structure(model)
+    Avals = jac_coord(model, x)
+    QuadraticModel(g, Hrows, Hcols, Hvals, c0=c0,
+                   Arows=Arows, Acols=Acols, Avals=Avals, lcon=model.meta.lcon .- c, ucon=model.meta.ucon .- c,
+                   lvar=model.meta.lvar .- x, uvar=model.meta.uvar .- x, x0=zeros(model.meta.nvar))
+  else
+    QuadraticModel(g, Hrows, Hcols, Hvals, c0=c0,
+                   lvar=model.meta.lvar .- x, uvar=model.meta.uvar .- x, x0=zeros(model.meta.nvar))
+  end
 end
 
+linobj(qp::AbstractQuadraticModel, args...) = qp.data.c
+
 function NLPModels.objgrad!(qp :: AbstractQuadraticModel, x :: AbstractVector, g :: AbstractVector)
-    v = qp.data.opH * x
-    @. g = qp.data.c + v
-    f = qp.data.c0 + dot(qp.data.c, x) + dot(v, x) / 2
-    qp.counters.neval_hprod += 1
-    (f, g)
+  NLPModels.increment!(qp, :neval_obj)
+  NLPModels.increment!(qp, :neval_grad)
+  coo_sym_prod!(qp.data.Hrows, qp.data.Hcols, qp.data.Hvals, x, g)
+  f = qp.data.c0 + dot(qp.data.c, x) + dot(g, x) / 2
+  @. g .+= qp.data.c
+  return f, g
 end
 
 function NLPModels.obj(qp :: AbstractQuadraticModel, x :: AbstractVector)
-    v = qp.data.opH * x
-    f = qp.data.c0 + dot(qp.data.c, x) + dot(v, x) / 2
-    qp.counters.neval_hprod += 1
-    f
+  NLPModels.increment!(qp, :neval_obj)
+  Hx = zeros(qp.meta.nvar)
+  coo_sym_prod!(qp.data.Hrows, qp.data.Hcols, qp.data.Hvals, x, Hx)
+  return qp.data.c0 + dot(qp.data.c, x) + dot(Hx, x) / 2
 end
 
 function NLPModels.grad!(qp :: AbstractQuadraticModel, x :: AbstractVector, g :: AbstractVector)
-    v = qp.data.opH * x
-    @. g = qp.data.c + v
-    qp.counters.neval_hprod += 1
-    g
+  NLPModels.increment!(qp, :neval_grad)
+  coo_sym_prod!(qp.data.Hrows, qp.data.Hcols, qp.data.Hvals, x, g)
+  g .+= qp.data.c
+  return g
 end
 
-hess(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...) = qp.data.H
-
-hess_op(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...) = qp.data.opH
+# TODO: Better hess_op
 
-"""
-Return the structure of the Lagrangian Hessian in sparse coordinate format in place.
-"""
-function NLPModels.hess_structure!(qp :: QuadraticModel, rows :: Vector{<: Integer}, cols :: Vector{<: Integer})
-    rows .= qp.data.H.rowval
-    cols .= findnz(qp.data.H)[2]
+function NLPModels.hess_structure!(qp :: QuadraticModel, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
+  rows .= qp.data.Hrows
+  cols .= qp.data.Hcols
+  return rows, cols
 end
 
-"""
-Evaluate the Lagrangian Hessian at `x` in sparse coordinate format. Only the lower triangle is returned.
-"""
-function NLPModels.hess_coord!(qp :: QuadraticModel, :: AbstractVector, rows :: AbstractVector{<: Integer},
-                               cols :: AbstractVector{<: Integer}, vals :: Vector{<: AbstractFloat}; kwargs...)
-    vals .= findnz(qp.data.H)[3]
+function NLPModels.hess_coord!(qp :: QuadraticModel, x :: AbstractVector, vals :: AbstractVector{<: AbstractFloat}; obj_weight :: Real=one(eltype(x)))
+  NLPModels.increment!(qp, :neval_hess)
+  vals .= obj_weight * qp.data.Hvals
+  return vals
 end
 
-jac(qp :: AbstractQuadraticModel, ::AbstractVector) = qp.data.A
-
-jac_op(qp :: AbstractQuadraticModel, ::AbstractVector) = LinearOperator(qp.data.A)
+NLPModels.hess_coord!(qp :: QuadraticModel, x :: AbstractVector, y :: AbstractVector, vals :: AbstractVector{<: AbstractFloat}; obj_weight :: Real=one(eltype(x))) = hess_coord!(qp, x, vals, obj_weight=obj_weight)
 
-jac_coord(qp :: AbstractQuadraticModel, ::AbstractVector) = findnz(qp.data.A)
-
-"""
-Return the structure of the constraints Jacobian in sparse coordinate format in place.
-"""
-function NLPModels.jac_structure!(qp :: QuadraticModel, rows :: Vector{<: Integer}, cols :: Vector{<: Integer})
-    rows .= qp.data.A.rowval
-    cols .= findnz(qp.data.A)[2]
+function NLPModels.jac_structure!(qp :: QuadraticModel, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
+  rows .= qp.data.Arows
+  cols .= qp.data.Acols
+  return rows, cols
 end
 
-"""
-Return the structure of the constraints Jacobian in sparse coordinate format in place.
-"""
-function NLPModels.jac_coord!(qp :: QuadraticModel, x :: AbstractVector, rows :: Vector{<: Integer},
-                              cols :: Vector{<: Integer}, vals :: Vector{<: AbstractFloat})
-    vals .= findnz(qp.data.A)[3]
+function NLPModels.jac_coord!(qp :: QuadraticModel, x :: AbstractVector, vals :: AbstractVector{<: AbstractFloat})
+  NLPModels.increment!(qp, :neval_jac)
+  vals .= qp.data.Avals
+  return vals
 end
 
 function NLPModels.cons!(qp :: AbstractQuadraticModel, x :: AbstractVector, c :: AbstractVector)
-    mul!(c, qp.data.A, x)
-    qp.counters.neval_jprod += 1
-    c
+  NLPModels.increment!(qp, :neval_cons)
+  coo_prod!(qp.data.Arows, qp.data.Acols, qp.data.Avals, x, c)
+  return c
 end
 
-function NLPModels.hprod!(qp :: AbstractQuadraticModel, x :: AbstractVector, v :: AbstractVector, Av :: AbstractVector)
-    qp.counters.neval_hprod += 1
-    Av .= qp.data.opH * v
-    return Av
+function NLPModels.hprod!(qp :: AbstractQuadraticModel, x :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight :: Real = one(eltype(x)))
+  NLPModels.increment!(qp, :neval_hprod)
+  coo_sym_prod!(qp.data.Hrows, qp.data.Hcols, qp.data.Hvals, v, Hv)
+  if obj_weight != 1
+    Hv .*= obj_weight
+  end
+  return Hv
 end
 
+NLPModels.hprod!(qp :: AbstractQuadraticModel, x :: AbstractVector, y :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight :: Real = one(eltype(x))) = hprod!(qp, x, v, Hv, obj_weight=obj_weight)
+
 function NLPModels.jprod!(qp :: AbstractQuadraticModel, x :: AbstractVector, v :: AbstractVector, Av :: AbstractVector)
-    qp.counters.neval_jprod += 1
-    mul!(Av, qp.data.A, v)
-    return Av
+  NLPModels.increment!(qp, :neval_jprod)
+  coo_prod!(qp.data.Arows, qp.data.Acols, qp.data.Avals, v, Av)
+  return Av
 end
 
 function NLPModels.jtprod!(qp :: AbstractQuadraticModel, x :: AbstractVector, v :: AbstractVector, Atv :: AbstractVector)
-    qp.counters.neval_jtprod += 1
-    mul!(Atv, qp.data.A', v)
-    return Atv
+  NLPModels.increment!(qp, :neval_jtprod)
+  coo_prod!(qp.data.Acols, qp.data.Arows, qp.data.Avals, v, Atv)
+  return Atv
 end

test/runtests.jl

@@ -2,7 +2,13 @@
 using Printf, SparseArrays, Test
 
 # our packages
-using LinearOperators, NLPModels, QuadraticModels
+using LinearAlgebra, LinearOperators, NLPModels, QuadraticModels
 
+nlpmodels_path = joinpath(dirname(pathof(NLPModels)), "..", "test")
+nlpmodels_problems_path = joinpath(nlpmodels_path, "problems")
+
+# Definition of the quadratic problem in ADNLPModel
 include("simpleqp.jl")
+
+include(joinpath(nlpmodels_path, "consistency.jl"))
 include("test_consistency.jl")

test/simpleqp.jl

@@ -2,13 +2,26 @@
 function simpleqp_autodiff()
 
   x0   = [-1.; 1.]
-  f(x) = 4*x[1]^2 + 5*x[2]^2 + x[1]*x[2] + 1.5*x[1] - 2*x[2]
+  f(x) = 4*x[1]^2 + 5*x[2]^2 + x[1]*x[2] + 1.5*x[1] - 2*x[2] + 1.0
   c(x) = [2*x[1] + x[2]; -x[1] + 2*x[2]]
   lcon = [2.0; -Inf]
   ucon = [Inf; 6.0]
   uvar = [20; Inf]
   lvar = [0.0; 0.0]
 
   return ADNLPModel(f, x0, c=c, lcon=lcon, ucon=ucon, lvar=lvar, uvar=uvar)
+end
+
+function simpleqp_QP()
+  c     = [1.5; -2.0]
+  H     = [8.0 0.0; 1.0 10.0]
+  A     = [2.0 1.0; -1.0 2.0]
+  c0    = 1.0
+  lcon  = [2.0; -Inf]
+  ucon  = [Inf; 6.0]
+  uvar  = [20; Inf]
+  lvar  = [0.0; 0.0]
+  x0    = [-1.0; 1.0]
 
+  return QuadraticModel(c, H, A=A, lcon=lcon, ucon=ucon, lvar=lvar, uvar=uvar, c0=c0, x0=x0)
 end