Bump to LinearOperators 2 and NLPModels 0.17 (#21)

* compat linop and NLPModels

* add more tests

* fix issue test nls

* type issue hess multi precision

* return res

* simplify test slacknls

Author: geoffroyleconte

Project.toml

@@ -12,8 +12,8 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 FastClosures = "0.2.1, 0.3.0"
-LinearOperators = "1.1.0"
-NLPModels = "0.16"
+LinearOperators = "2.0"
+NLPModels = "0.17"
 julia = "^1.3.0"
 
 [extras]

src/feasibility-form-nls.jl

@@ -233,7 +233,7 @@ function NLPModels.hess_coord!(
   xr::AbstractVector,
   y::AbstractVector,
   vals::AbstractVector;
-  obj_weight::Float64 = 1.0,
+  obj_weight::Real = one(eltype(xr)),
 )
   @lencheck nlp.meta.nvar xr
   @lencheck nlp.meta.ncon y
@@ -261,16 +261,17 @@ function NLPModels.hprod!(
   y::AbstractVector,
   v::AbstractVector,
   hv::AbstractVector;
-  obj_weight::Float64 = 1.0,
+  obj_weight::Real = one(eltype(xr)),
 )
   @lencheck nlp.meta.nvar xr v hv
   @lencheck nlp.meta.ncon y
   n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
   x = @view xr[1:n]
+  T = eltype(xr)
   if m > 0
-    @views hprod!(nlp.internal, x, y[(ne + 1):end], v[1:n], hv[1:n], obj_weight = 0.0)
+    @views hprod!(nlp.internal, x, y[(ne + 1):end], v[1:n], hv[1:n], obj_weight = zero(T))
   else
-    fill!(hv, 0.0)
+    fill!(hv, zero(T))
   end
   for i = 1:ne
     @views hv[1:n] .+= hprod_residual(nlp.internal, x, i, v[1:n]) * y[i]
@@ -356,7 +357,7 @@ function NLPModels.jtprod_residual!(
   @lencheck nlp.nls_meta.nequ v
   increment!(nlp, :neval_jtprod_residual)
   n, ne = nlp.internal.meta.nvar, nlp.internal.nls_meta.nequ
-  Jtv[1:n] .= 0.0
+  Jtv[1:n] .= zero(eltype(x))
   Jtv[(n + 1):end] .= v
   return Jtv
 end
@@ -387,7 +388,7 @@ function NLPModels.jth_hess_residual(nlp::FeasibilityFormNLS, x::AbstractVector,
   @lencheck nlp.meta.nvar x
   increment!(nlp, :neval_jhess_residual)
   n = nlp.meta.nvar
-  return spzeros(n, n)
+  return spzeros(eltype(x), n, n)
 end
 
 function NLPModels.hprod_residual!(
@@ -399,6 +400,6 @@ function NLPModels.hprod_residual!(
 )
   @lencheck nlp.meta.nvar x v Hiv
   increment!(nlp, :neval_hprod_residual)
-  fill!(Hiv, 0.0)
+  fill!(Hiv, zero(eltype(x)))
   return Hiv
 end

src/feasibility-residual.jl

@@ -128,7 +128,7 @@ end
 
 function NLPModels.hess_residual(nls::FeasibilityResidual, x::AbstractVector, v::AbstractVector)
   increment!(nls, :neval_hess_residual)
-  return hess(nls.nlp, x, v, obj_weight = 0.0)
+  return hess(nls.nlp, x, v, obj_weight = zero(eltype(x)))
 end
 
 function NLPModels.hess_structure_residual!(
@@ -146,14 +146,15 @@ function NLPModels.hess_coord_residual!(
   vals::AbstractVector,
 )
   increment!(nls, :neval_hess_residual)
-  return hess_coord!(nls.nlp, x, v, vals, obj_weight = 0.0)
+  return hess_coord!(nls.nlp, x, v, vals, obj_weight = zero(eltype(x)))
 end
 
 function NLPModels.jth_hess_residual(nls::FeasibilityResidual, x::AbstractVector, i::Int)
   increment!(nls, :neval_jhess_residual)
-  y = zeros(nls.nls_meta.nequ)
-  y[i] = 1.0
-  return hess(nls.nlp, x, y, obj_weight = 0.0)
+  T = eltype(x)
+  y = zeros(T, nls.nls_meta.nequ)
+  y[i] = one(T)
+  return hess(nls.nlp, x, y, obj_weight = zero(T))
 end
 
 function NLPModels.hprod_residual!(
@@ -164,9 +165,10 @@ function NLPModels.hprod_residual!(
   Hiv::AbstractVector,
 )
   increment!(nls, :neval_hprod_residual)
-  y = zeros(nls.nls_meta.nequ)
-  y[i] = 1.0
-  return hprod!(nls.nlp, x, y, v, Hiv, obj_weight = 0.0)
+  T = eltype(x)
+  y = zeros(T, nls.nls_meta.nequ)
+  y[i] = one(T)
+  return hprod!(nls.nlp, x, y, v, Hiv, obj_weight = zero(T))
 end
 
 function NLPModels.hess(
@@ -178,7 +180,7 @@ function NLPModels.hess(
   cx = cons(nls.nlp, x)
   Jx = jac(nls.nlp, x)
   Hx = tril(Jx' * Jx)
-  Hx .+= hess(nls.nlp, x, cx, obj_weight = 0.0).data
+  Hx .+= hess(nls.nlp, x, cx, obj_weight = zero(eltype(x))).data
   return Symmetric(obj_weight * Hx, :L)
 end
 
@@ -193,8 +195,9 @@ function NLPModels.hprod!(
   cx = cons(nls.nlp, x)
   Jv = jprod(nls.nlp, x, v)
   jtprod!(nls.nlp, x, Jv, Hv)
-  Hiv = zeros(eltype(x), nls.meta.nvar)
-  hprod!(nls.nlp, x, cx, v, Hiv, obj_weight = 0.0)
+  T = eltype(x)
+  Hiv = zeros(T, nls.meta.nvar)
+  hprod!(nls.nlp, x, cx, v, Hiv, obj_weight = zero(T))
   Hv .+= Hiv
   Hv .*= obj_weight
   return Hv

src/quasi-newton.jl

@@ -84,7 +84,7 @@ function NLPModels.hprod!(
   Hv::AbstractVector;
   kwargs...,
 )
-  Hv[1:(nlp.meta.nvar)] .= nlp.op * v
+  @views mul!(Hv[1:(nlp.meta.nvar)], nlp.op, v)
   return Hv
 end
 

src/slack-model.jl

@@ -108,6 +108,7 @@ function slack_meta(meta::NLPModelMeta{T, S}; name = meta.name * "-slack") where
     lin = meta.lin,
     nln = meta.nln,
     name = name,
+    minimize = meta.minimize,
   )
 end
 
@@ -465,16 +466,32 @@ function NLPModels.jac_op_residual!(
 )
   @lencheck nls.meta.nvar x Jtv
   @lencheck nls.nls_meta.nequ Jv
-  prod = @closure v -> jprod_residual!(nls, x, v, Jv)
-  ctprod = @closure v -> jtprod_residual!(nls, x, v, Jtv)
+  prod! = @closure (res, v, α, β) -> begin
+    jprod_residual!(nls, x, v, Jv)
+    if β == 0
+      @. res = α * Jv
+    else
+      @. res = α * Jv + β * res
+    end
+    return res
+  end
+  ctprod! = @closure (res, v, α, β) -> begin
+    jtprod_residual!(nls, x, v, Jtv)
+    if β == 0
+      @. res = α * Jtv
+    else
+      @. res = α * Jtv + β * res
+    end
+    return res
+  end
   return LinearOperator{eltype(x)}(
     nls_meta(nls).nequ,
     nls_meta(nls).nvar,
     false,
     false,
-    prod,
-    ctprod,
-    ctprod,
+    prod!,
+    ctprod!,
+    ctprod!,
   )
 end
 
@@ -550,14 +567,22 @@ function NLPModels.hess_op_residual!(
   Hiv::AbstractVector,
 )
   @lencheck nls.meta.nvar x Hiv
-  prod = @closure v -> hprod_residual!(nls, x, i, v, Hiv)
+  prod! = @closure (res, v, α, β) -> begin
+    hprod_residual!(nls, x, i, v, Hiv)
+    if β == 0
+      @. res = α * Hiv
+    else
+      @. res = α * Hiv + β * res
+    end
+    return res
+  end
   return LinearOperator{eltype(x)}(
     nls_meta(nls).nvar,
     nls_meta(nls).nvar,
     true,
     true,
-    prod,
-    prod,
-    prod,
+    prod!,
+    prod!,
+    prod!,
   )
 end

test/nlp/quasi-newton.jl

@@ -12,21 +12,21 @@
       n = nlp.meta.nvar
       m = nlp.meta.ncon
 
-      s, y = randn(n), randn(n)
+      s, y = randn(T, n), randn(T, n)
       B = QNO(n)
       push!(B, s, y)
       push!(nlp, s, y)
       H(x) = B
       H(x, y) = B
 
-      y = randn(m)
-      x = randn(n)
-      v = randn(n)
-      w = randn(m)
-      Jv = zeros(m)
-      Jtw = zeros(n)
-      Hv = zeros(n)
-      Hvals = zeros(nlp.meta.nnzh)
+      y = randn(T, m)
+      x = randn(T, n)
+      v = randn(T, n)
+      w = randn(T, m)
+      Jv = zeros(T, m)
+      Jtw = zeros(T, n)
+      Hv = zeros(T, n)
+      Hvals = zeros(T, nlp.meta.nnzh)
 
       # Basic methods
       @test obj(nlp, x) ≈ f(x)

test/nlp/slack-model.jl

@@ -1,24 +1,24 @@
 @testset "SlackModel NLP tests" begin
   @testset "API" for T in [Float64, Float32]
     f(x) = (x[1] - 2)^2 + (x[2] - 1)^2
-    ∇f(x) = [2 * (x[1] - 2); 2 * (x[2] - 1); 0]
-    H(x) = [2.0 0 0; 0 2.0 0; 0 0 0]
-    c(x) = [x[1] - 2x[2] + 1; -x[1]^2 / 4 - x[2]^2 + 1 - x[3]]
-    J(x) = [1.0 -2.0 0; -0.5x[1] -2.0x[2] -1]
-    H(x, y) = H(x) + y[2] * [-0.5 0 0; 0 -2.0 0; 0 0 0]
+    ∇f(x) = T[2 * (x[1] - 2); 2 * (x[2] - 1); 0]
+    H(x) = T[2.0 0 0; 0 2.0 0; 0 0 0]
+    c(x) = T[x[1] - 2x[2] + 1; -x[1]^2 / 4 - x[2]^2 + 1 - x[3]]
+    J(x) = T[1.0 -2.0 0; -0.5x[1] -2.0x[2] -1]
+    H(x, y) = H(x) + y[2] * T[-0.5 0 0; 0 -2.0 0; 0 0 0]
 
     nlp = SlackModel(SimpleNLPModel(T))
     n = nlp.meta.nvar
     m = nlp.meta.ncon
 
-    x = randn(n)
-    y = randn(m)
-    v = randn(n)
-    w = randn(m)
-    Jv = zeros(m)
-    Jtw = zeros(n)
-    Hv = zeros(n)
-    Hvals = zeros(nlp.meta.nnzh)
+    x = randn(T, n)
+    y = randn(T, m)
+    v = randn(T, n)
+    w = randn(T, m)
+    Jv = zeros(T, m)
+    Jtw = zeros(T, n)
+    Hv = zeros(T, n)
+    Hvals = zeros(T, nlp.meta.nnzh)
 
     # Basic methods
     @test obj(nlp, x) ≈ f(x)
@@ -52,15 +52,23 @@
     Jop = jac_op!(nlp, x, Jv, Jtw)
     @test Jop * v ≈ J(x) * v
     @test Jop' * w ≈ J(x)' * w
+    res = J(x) * v - w
+    @test mul!(w, Jop, v, one(T), -one(T)) ≈ res
+    res = J(x)' * w - v
+    @test mul!(v, Jop', w, one(T), -one(T)) ≈ res
     Jop = jac_op!(nlp, jac_structure(nlp)..., jac_coord(nlp, x), Jv, Jtw)
     @test Jop * v ≈ J(x) * v
     @test Jop' * w ≈ J(x)' * w
+    res = J(x) * v - w
+    @test mul!(w, Jop, v, one(T), -one(T)) ≈ res
+    res = J(x)' * w - v
+    @test mul!(v, Jop', w, one(T), -one(T)) ≈ res
     Jop = jac_op!(nlp, x, jac_structure(nlp)..., Jv, Jtw)
     @test Jop * v ≈ J(x) * v
     @test Jop' * w ≈ J(x)' * w
-    ghjv = zeros(m)
+    ghjv = zeros(T, m)
     for j = 1:m
-      eⱼ = [i == j ? 1.0 : 0.0 for i = 1:m]
+      eⱼ = [i == j ? one(T) : zero(T) for i = 1:m]
       Cⱼ(x) = H(x, eⱼ) - H(x)
       ghjv[j] = dot(gx, Cⱼ(x) * v)
     end
@@ -73,14 +81,21 @@
     @test Hop * v ≈ H(x) * v
     Hop = hess_op!(nlp, x, Hv)
     @test Hop * v ≈ H(x) * v
+    z = ones(T, nlp.meta.nvar)
+    res = H(x) * v - z
+    @test mul!(z, Hop, v, one(T), -one(T)) ≈ res
     Hop = hess_op!(nlp, hess_structure(nlp)..., hess_coord(nlp, x), Hv)
     @test Hop * v ≈ H(x) * v
+    res = H(x) * v - z
+    @test mul!(z, Hop, v, one(T), -one(T)) ≈ res
     Hop = hess_op!(nlp, x, hess_structure(nlp)..., Hv)
     @test Hop * v ≈ H(x) * v
     Hop = hess_op(nlp, x, y)
     @test Hop * v ≈ H(x, y) * v
     Hop = hess_op!(nlp, x, y, Hv)
     @test Hop * v ≈ H(x, y) * v
+    res = H(x, y) * v - z
+    @test mul!(z, Hop, v, one(T), -one(T)) ≈ res
     Hop = hess_op!(nlp, hess_structure(nlp)..., hess_coord(nlp, x, y), Hv)
     @test Hop * v ≈ H(x, y) * v
     Hop = hess_op!(nlp, x, y, hess_structure(nlp)..., Hv)