Merit functions as NLPModel

.gitignore

@@ -1,2 +1,4 @@
 docs/build
 docs/Manifest.toml
+
+Manifest.toml

docs/src/api.md

@@ -27,6 +27,17 @@ redirect!
 armijo_wolfe
 ```
 
+## Merit
+
+See also [`obj`](@ref).
+
+```@docs
+AbstractMeritModel
+derivative
+L1Merit
+AugLagMerit
+```
+
 ## Stats
 
 ```@docs

docs/src/reference.md

@@ -1,4 +1,12 @@
 # Reference
 
-```@index
+## Contents
+
+```@contents
+Pages = ["reference.md"]
 ```
+
+## Index
+
+```@index
+```

src/SolverTools.jl

@@ -16,6 +16,7 @@ include("auxiliary/logger.jl")
 include("stats/stats.jl")
 
 # Algorithmic components.
+include("merit/merit.jl")
 include("linesearch/linesearch.jl")
 include("trust-region/trust-region.jl")
 

src/linesearch/line_model.jl

@@ -1,5 +1,3 @@
-import NLPModels: obj, grad, grad!, hess
-
 export LineModel
 export obj, grad, derivative, grad!, derivative!, hess, redirect!
 
@@ -38,7 +36,7 @@ end
 
     ϕ(t) := f(x + td).
 """
-function obj(f :: LineModel, t :: AbstractFloat)
+function NLPModels.obj(f :: LineModel, t :: AbstractFloat)
   NLPModels.increment!(f, :neval_obj)
   return obj(f.nlp, f.x + t * f.d)
 end
@@ -51,7 +49,7 @@ i.e.,
 
     ϕ'(t) = ∇f(x + td)ᵀd.
 """
-function grad(f :: LineModel, t :: AbstractFloat)
+function NLPModels.grad(f :: LineModel, t :: AbstractFloat)
   NLPModels.increment!(f, :neval_grad)
   return dot(grad(f.nlp, f.x + t * f.d), f.d)
 end
@@ -67,7 +65,7 @@ i.e.,
 
 The gradient ∇f(x + td) is stored in `g`.
 """
-function grad!(f :: LineModel, t :: AbstractFloat, g :: AbstractVector)
+function NLPModels.grad!(f :: LineModel, t :: AbstractFloat, g :: AbstractVector)
   NLPModels.increment!(f, :neval_grad)
   return dot(grad!(f.nlp, f.x + t * f.d, g), f.d)
 end
@@ -81,7 +79,7 @@ i.e.,
 
     ϕ"(t) = dᵀ∇²f(x + td)d.
 """
-function hess(f :: LineModel, t :: AbstractFloat)
+function NLPModels.hess(f :: LineModel, t :: AbstractFloat)
   NLPModels.increment!(f, :neval_hess)
   return dot(f.d, hprod(f.nlp, f.x + t * f.d, f.d))
 end

src/merit/auglagmerit.jl

@@ -0,0 +1,115 @@
+export AugLagMerit
+
+@doc raw"""
+  AugLagMerit(nlp, η; kwargs...)
+
+Creates an augmented Lagrangian merit function for the equality constrained problem
+```math
+\min f(x) \quad \text{s.to} \quad c(x) = 0
+```
+defined by
+```math
+\phi(x, yₖ; η) = f(x) + yₖᵀc(x) + η ½\|c(x)\|²
+```
+
+In addition to the keyword arguments declared in [`AbstractMeritModel`](@ref), an `AugLagMerit` also
+accepts the argument `y`.
+"""
+mutable struct AugLagMerit{M <: AbstractNLPModel, T <: Real, V <: AbstractVector} <: AbstractMeritModel
+  meta :: NLPModelMeta
+  counters :: Counters
+  nlp  :: M
+  η  :: T
+  fx   :: T
+  gx   :: V
+  cx   :: V
+  Ad   :: V
+  y  :: V
+  y⁺   :: V
+  Jᵀy⁺ :: V
+  Jv   :: V
+  JᵀJv   :: V
+end
+
+function AugLagMerit(
+  nlp  :: M,
+  η  :: T;
+  fx   :: T = T(Inf),
+  gx   :: V = fill(T(Inf), nlp.meta.nvar),
+  cx   :: V = fill(T(Inf), nlp.meta.ncon),
+  Ad   :: V = fill(T(Inf), nlp.meta.ncon),
+  y  :: V = fill(T(Inf), nlp.meta.ncon),
+  y⁺   :: V = fill(T(Inf), nlp.meta.ncon), # y + η * c(x)
+  Jᵀy⁺ :: V = fill(T(Inf), nlp.meta.nvar),
+  Jv   :: V = fill(T(Inf), nlp.meta.ncon),
+  JᵀJv :: V = fill(T(Inf), nlp.meta.nvar)
+) where {M <: AbstractNLPModel, T <: Real, V <: AbstractVector{<: T}}
+  meta = NLPModelMeta(nlp.meta.nvar)
+  AugLagMerit{M,T,V}(meta, Counters(), nlp, η, fx, gx, cx, Ad, y, y⁺, Jᵀy⁺, Jv, JᵀJv)
+end
+
+function NLPModels.obj(merit :: AugLagMerit, x :: AbstractVector; update :: Bool = true)
+  @lencheck merit.meta.nvar x
+  NLPModels.increment!(merit, :neval_obj)
+  if update
+    merit.fx = obj(merit.nlp, x)
+    merit.nlp.meta.ncon > 0 && cons!(merit.nlp, x, merit.cx)
+  end
+  return merit.fx + dot(merit.y, merit.cx) + merit.η * dot(merit.cx, merit.cx) / 2
+end
+
+function derivative(merit :: AugLagMerit, x :: AbstractVector, d :: AbstractVector; update :: Bool = true)
+  @lencheck merit.meta.nvar x d
+  if update
+    grad!(merit.nlp, x, merit.gx)
+    merit.nlp.meta.ncon > 0 && jprod!(merit.nlp, x, d, merit.Ad)
+  end
+  if merit.nlp.meta.ncon == 0
+    return dot(merit.gx, d)
+  else
+    return dot(merit.gx, d) + dot(merit.y, merit.Ad) + merit.η * dot(merit.cx, merit.Ad)
+  end
+end
+
+function NLPModels.grad!(merit :: AugLagMerit, x :: AbstractVector, g :: AbstractVector; update :: Bool = true)
+  @lencheck merit.meta.nvar x g
+  NLPModels.increment!(merit, :neval_grad)
+  if update
+    grad!(nlp.model, x, merit.gx)
+    merit.nlp.meta.ncon > 0 && cons!(merit.nlp, x, merit.cx)
+    merit.y⁺ .= merit.y .+ merit.η .* merit.cx
+    merit.nlp.meta.ncon > 0 && jtprod!(merit.nlp, x, merit.y⁺, merit.Jᵀy⁺)
+  end
+  g .= merit.gx .+ merit.Jᵀy⁺
+  return g
+end
+
+function NLPModels.objgrad!(merit :: AugLagMerit, x :: AbstractVector, g :: AbstractVector)
+  @lencheck merit.meta.nvar x g
+  NLPModels.increment!(merit, :neval_obj)
+  NLPModels.increment!(merit, :neval_grad)
+  if update
+    merit.fx = obj(merit.nlp, x)
+    grad!(nlp.model, x, merit.gx)
+    merit.nlp.meta.ncon > 0 && cons!(merit.nlp, x, merit.cx)
+    merit.y⁺ .= merit.y .+ merit.η .* merit.cx
+    merit.nlp.meta.ncon > 0 && jtprod!(merit.nlp, x, merit.y⁺, merit.Jᵀy⁺)
+  end
+  f = merit.fx + dot(merit.y, merit.cx) + merit.η * dot(merit.cx, merit.cx) / 2
+  g .= merit.gx .+ merit.Jᵀy⁺
+  return f, g
+end
+
+function NLPModels.hprod!(merit :: AugLagMerit, x :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight :: Float64 = 1.0)
+  @lencheck merit.meta.nvar x v Hv
+  NLPModels.increment!(merit, :neval_hprod)
+  if update
+    merit.nlp.meta.ncon > 0 && cons!(merit.nlp, x, merit.cx)
+    merit.y⁺ .= merit.y .+ merit.η .* merit.cx
+  end
+  jprod!(merit.model, x, v, merit.Jv)
+  jtprod!(merit.model, x, merit.Jv, merit.JᵀJv)
+  hprod!(merit.model, x, merit.y⁺, v, Hv, obj_weight = obj_weight)
+  Hv .+= merit.η * merit.JᵀJv
+  return Hv
+end

src/merit/l1merit.jl

@@ -0,0 +1,68 @@
+export L1Merit
+
+@doc raw"""
+    L1Merit(nlp, η; kwargs...)
+
+Creates a ℓ₁ merit function for the equality constrained problem
+```math
+\min f(x) \quad \text{s.to} \quad c(x) = 0
+```
+defined by
+```math
+\phi_1(x; η) = f(x) + η\|c(x)\|₁
+```
+"""
+mutable struct L1Merit{M <: AbstractNLPModel, T <: Real, V <: AbstractVector} <: AbstractMeritModel
+  meta :: NLPModelMeta
+  counters :: Counters
+  nlp :: M
+  η :: T
+  fx :: T
+  gx :: V
+  cx :: V
+  Ad :: V
+end
+
+function L1Merit(
+  nlp :: M,
+  η :: T;
+  fx :: T = T(Inf),
+  gx :: V = fill(T(Inf), nlp.meta.nvar),
+  cx :: V = fill(T(Inf), nlp.meta.ncon),
+  Ad :: V = fill(T(Inf), nlp.meta.ncon)
+) where {M <: AbstractNLPModel, T <: Real, V <: AbstractVector{<: T}}
+  meta = NLPModelMeta(nlp.meta.nvar)
+  L1Merit{M,T,V}(meta, Counters(), nlp, η, fx, gx, cx, Ad)
+end
+
+function NLPModels.obj(merit :: L1Merit, x :: AbstractVector; update :: Bool = true)
+  @lencheck merit.meta.nvar x
+  NLPModels.increment!(merit, :neval_obj)
+  if update
+    merit.fx = obj(merit.nlp, x)
+    merit.nlp.meta.ncon > 0 && cons!(merit.nlp, x, merit.cx)
+  end
+  return merit.fx + merit.η * norm(merit.cx, 1)
+end
+
+function derivative(merit :: L1Merit, x :: AbstractVector, d :: AbstractVector; update :: Bool = true)
+  @lencheck merit.meta.nvar x d
+  if update
+    grad!(merit.nlp, x, merit.gx)
+    merit.nlp.meta.ncon > 0 && jprod!(merit.nlp, x, d, merit.Ad)
+  end
+  if merit.nlp.meta.ncon == 0
+    return dot(merit.gx, d)
+  else
+    return dot(merit.gx, d) + merit.η * sum(
+      if ci > 0
+        Adi
+      elseif ci < 0
+        -Adi
+      else
+        abs(Adi)
+      end
+      for (ci, Adi) in zip(merit.cx, merit.Ad)
+    )
+  end
+end

src/merit/merit.jl

@@ -0,0 +1,40 @@
+using NLPModels
+export AbstractMeritModel, obj, derivative
+
+"""
+    AbstractMeritModel
+
+Model for merit functions. All models should store
+- `nlp`: The NLP with the corresponding problem.
+- `η`: The merit parameter.
+- `fx`: The objective at some point.
+- `cx`: The constraints vector at some point.
+- `gx`: The constraints vector at some point.
+- `Ad`: The Jacobian-direction product.
+
+All models allow a constructor of form
+
+    Merit(nlp, η; fx=Inf, cx=[Inf,…,Inf], gx=[Inf,…,Inf], Ad=[Inf,…,Inf])
+
+Additional arguments and constructors may be provided.
+
+An AbstractMeritModel is an AbstractNLPModel, but the API may not be completely implemented. For
+instance, the `L1Merit` model doesn't provide any gradient function, but it provides a directional
+derivative function.
+
+Furthermore, all implemented methods accept an `update` keyword that defaults to `true`. It is used
+to determine whether the internal stored values should be updated or not.
+"""
+abstract type AbstractMeritModel <: AbstractNLPModel end
+
+"""
+    derivative(merit, x, d; update=true)
+
+Computes the directional derivative of `merit` at `x` on direction `d`.
+This will call `grad!` and `jprod` to update the internal values of `gx` and `Ad`, but will assume that `cx` is correct.
+The option exist to allow updating the `η` parameter without recomputing `fx` and `cx`.
+"""
+function derivative end
+
+include("auglagmerit.jl")
+include("l1merit.jl")