add missing methods and fix tests

Author: RenanOD

Project.toml

@@ -8,14 +8,14 @@ FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
-NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-FastClosures = "≥ 0.2.1"
-LinearOperators = "0.5.0"
-NLPModels = "0.6.1, < 0.9.0"
-Requires = "≥ 0.3.0"
+FastClosures = "0.2.1, 0.3.0"
+LinearOperators = "0.7.0, 1"
+NLPModels = "0.10.0, 1"
+Requires = "0.3, 0.4, 0.5"
 julia = "^1.0.0"

src/qpmodel.jl

@@ -1,4 +1,3 @@
-# override some NLPModels functions
 export jac_structure!, hess_structure!, jac_coord!, hess_coord!
 
 mutable struct QPData
@@ -36,46 +35,41 @@ mutable struct QuadraticModel <: AbstractQuadraticModel
   end
 end
 
-function QuadraticModel(model :: AbstractNLPModel)
+function QuadraticModel(model :: AbstractNLPModel, x :: AbstractVector = model.meta.x0, args...; kwargs...)
     nvar = model.meta.nvar
     ncon = model.meta.ncon
-    z = zeros(nvar)
-    c = grad(model, z)
-    H = hess(model, model.meta.x0)
-    A = jac(model, z)
-    QuadraticModel(c, H, opHermitian(H), A,
-                   model.meta.lcon, model.meta.ucon,
-                   model.meta.lvar, model.meta.uvar)
+    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)
 end
 
 linobj(qp::AbstractQuadraticModel, args...) = qp.data.c
 
-function objgrad(qp :: AbstractQuadraticModel, x :: AbstractVector)
-    g = Vector{eltype(x)}(length(x))
+function NLPModels.objgrad(qp :: AbstractQuadraticModel, x :: AbstractVector)
+    g = Vector{eltype(x)}(undef, length(x))
     objgrad!(qp, x, g)
 end
 
-function objgrad!(qp :: AbstractQuadraticModel, x :: AbstractVector, g :: AbstractVector)
+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) + 0.5 * dot(v, x)
+    f = qp.data.c0 + dot(qp.data.c, x) + dot(v, x) / 2
     qp.counters.neval_hprod += 1
     (f, g)
 end
 
-function obj(qp :: AbstractQuadraticModel, x :: AbstractVector)
+function NLPModels.obj(qp :: AbstractQuadraticModel, x :: AbstractVector)
     v = qp.data.opH * x
-    f = qp.data.c0 + dot(qp.data.c, x) + 0.5 * dot(v, x)
+    f = qp.data.c0 + dot(qp.data.c, x) + dot(v, x) / 2
     qp.counters.neval_hprod += 1
     f
 end
 
-function grad(qp :: AbstractQuadraticModel, x :: AbstractVector)
-    g = Vector{eltype(x)}(undef, qp.meta.nvar)
-    grad!(qp, x, g)
-end
-
-function grad!(qp :: AbstractQuadraticModel, x :: AbstractVector, g :: AbstractVector)
+function NLPModels.grad!(qp :: AbstractQuadraticModel, x :: AbstractVector, g :: AbstractVector)
     v = qp.data.opH * x
     @. g = qp.data.c + v
     qp.counters.neval_hprod += 1
@@ -86,12 +80,10 @@ hess(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...) = qp.data.H
 
 hess_op(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...) = qp.data.opH
 
-hess_coord(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...) = findnz(qp.data.H)
-
 """
 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}; kwargs...)
+function NLPModels.hess_structure!(qp :: QuadraticModel, rows :: Vector{<: Integer}, cols :: Vector{<: Integer})
     rows .= qp.data.H.rowval
     cols .= findnz(qp.data.H)[2]
 end
@@ -104,16 +96,16 @@ function NLPModels.hess_coord!(qp :: QuadraticModel, :: AbstractVector, rows ::
     vals .= findnz(qp.data.H)[3]
 end
 
-jac(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...) = qp.data.A
+jac(qp :: AbstractQuadraticModel, ::AbstractVector) = qp.data.A
 
-jac_op(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...) = LinearOperator(qp.data.A)
+jac_op(qp :: AbstractQuadraticModel, ::AbstractVector) = LinearOperator(qp.data.A)
 
-jac_coord(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...) = findnz(qp.data.A)
+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}; kwargs...)
+function NLPModels.jac_structure!(qp :: QuadraticModel, rows :: Vector{<: Integer}, cols :: Vector{<: Integer})
     rows .= qp.data.A.rowval
     cols .= findnz(qp.data.A)[2]
 end
@@ -122,38 +114,30 @@ 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}; kwargs...)
+                              cols :: Vector{<: Integer}, vals :: Vector{<: AbstractFloat})
     vals .= findnz(qp.data.A)[3]
 end
 
-function cons(qp::AbstractQuadraticModel, x :: AbstractVector)
-    c = Vector{eltype(x)}(undef, qp.meta.ncon)
-    cons!(qp, x, c)
-end
-
-function cons!(qp :: AbstractQuadraticModel, x :: AbstractVector, c :: AbstractVector)
+function NLPModels.cons!(qp :: AbstractQuadraticModel, x :: AbstractVector, c :: AbstractVector)
     mul!(c, qp.data.A, x)
     qp.counters.neval_jprod += 1
     c
 end
 
-function hprod(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...)
-    @closure v -> begin
-        qp.counters.neval_hprod += 1
-        qp.data.opH * v
-    end
+function NLPModels.hprod!(qp :: AbstractQuadraticModel, x :: AbstractVector, v :: AbstractVector, Av :: AbstractVector)
+    qp.counters.neval_hprod += 1
+    Av .= qp.data.opH * v
+    return Av
 end
 
-function jprod(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...)
-    @closure v -> begin
-        qp.counters.neval_jprod += 1
-        qp.data.A * v
-    end
+function NLPModels.jprod!(qp :: AbstractQuadraticModel, x :: AbstractVector, v :: AbstractVector, Av :: AbstractVector)
+    qp.counters.neval_jprod += 1
+    mul!(Av, qp.data.A, v)
+    return Av
 end
 
-function jtprod(qp :: AbstractQuadraticModel, ::AbstractVector; kwargs...)
-    @closure v -> begin
-        qp.counters.neval_jtprod += 1
-        qp.data.A' * v
-    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
 end

test/runtests.jl

@@ -1,25 +1,8 @@
 # stdlib
-using SparseArrays, Test
+using Printf, SparseArrays, Test
 
 # our packages
-using LinearOperators, NLPModels, NLPModelsIpopt, QuadraticModels
+using LinearOperators, NLPModels, QuadraticModels
 
-@testset "QuadraticModelsTests" begin
-	c0        = 0.0
-	c         = [1.5; -2.0]
-	Q         = sparse([1; 2; 2], [1; 1; 2], [8.0; 2.0; 10.0])
-	A         = sparse([1; 2; 1; 2], [1; 1; 2; 2], [2.0; -1.0; 1.0; 2.0])
-	lcon      = [2.0; -Inf]
-	ucon      = [Inf; 6.0]
-	uvar      = [20; Inf]
-	lvar      = [0.0; 0.0]
-	objective = 4.3718750e+00
-
-	qp     = QuadraticModel(c, Q, opHermitian(Q), A, lcon, ucon, lvar, uvar, c0 = c0)
-	output = ipopt(qp, print_level = 0)
-
-	@test output.dual_feas                      < 1e-6
-	@test output.primal_feas                    < 1e-6
-	@test abs(output.objective - objective)     < 1e-6
-
-end # testset
+include("simpleqp.jl")
+include("test_consistency.jl")

test/simpleqp.jl

@@ -0,0 +1,14 @@
+
+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]
+  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

test/test_consistency.jl

@@ -0,0 +1,40 @@
+function check_quadratic_model(model, quadraticmodel)
+    @assert typeof(quadraticmodel) <: QuadraticModels.AbstractQuadraticModel
+    rtol  = 1e-8
+    @assert quadraticmodel.meta.nvar == model.meta.nvar
+    @assert quadraticmodel.meta.ncon == model.meta.ncon
+
+    x = [-(-1.0)^i for i = 1:quadraticmodel.meta.nvar]
+
+    @assert isapprox(obj(model, x), obj(quadraticmodel, x), rtol=rtol)
+
+    f, g = objgrad(model, x)
+    f_quad, g_quad = objgrad(quadraticmodel, x)
+
+    @assert isapprox(f, f_quad, rtol=rtol)
+    @assert isapprox(g, g_quad, rtol=rtol)
+    @assert isapprox(cons(model, x), cons(quadraticmodel, x), rtol=rtol)
+    @assert isapprox(jac(model, x), jac(quadraticmodel, x), rtol=rtol)
+
+    v = [-(-1.0)^i for i = 1:quadraticmodel.meta.nvar]
+    u = [-(-1.0)^i for i = 1:quadraticmodel.meta.ncon]
+
+    @assert isapprox(jprod(model, x, v), jprod(quadraticmodel, x, v), rtol=rtol)
+    @assert isapprox(jtprod(model, x, u), jtprod(quadraticmodel, x, u), rtol=rtol)
+
+    H = hess_op(quadraticmodel, x)
+    @assert typeof(H) <: LinearOperators.AbstractLinearOperator
+    @assert size(H) == (model.meta.nvar, model.meta.nvar)
+    @assert isapprox(H * v, hprod(model, x, v), rtol=rtol)
+
+    reset!(quadraticmodel)
+  end
+
+  for problem in ["simpleqp"]
+    problem_f = eval(Symbol(problem * "_autodiff"))
+    nlp = problem_f()
+    @printf("Checking QuadraticModel formulation of %-8s\t", problem)
+    quadratic_model = QuadraticModel(nlp, zeros(nlp.meta.nvar))
+    check_quadratic_model(nlp, quadratic_model)
+    @printf("✓\n")
+  end