From c386150c23091a0d3cb0253a1805795fefc8ef5a Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Sat, 1 Jun 2024 13:42:10 +0200
Subject: [PATCH 01/24] Add RegularizedNLPModels

---
 src/R2_alg.jl | 18 ++++++++++++++----
 1 file changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 88481fc9..0d543814 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -100,16 +100,26 @@ In the second form, instead of `nlp`, the user may pass in
 * `Hobj_hist`: an array with the history of values of the nonsmooth objective
 * `Complex_hist`: an array with the history of number of inner iterations.
 """
-function R2(nlp::AbstractNLPModel, args...; kwargs...)
+
+function R2(nlp::AbstractNLPModel, args...; kwargs)
+  kwargs_dict = Dict(kwargs...)
+  selected = pop!(kwargs_dict, :selected, 1:nlp.meta.nvar)
+  reg_nlp = RegularizedNLPModel(nlp, args.h, selected)
+  return R2(reg_nlp, args...; kwargs...)
+end
+
+function R2(reg_nlp::AbstractRegularizedNLPModel, args...; kwargs...)
   kwargs_dict = Dict(kwargs...)
-  x0 = pop!(kwargs_dict, :x0, nlp.meta.x0)
+  x0 = pop!(kwargs_dict, :x0, reg_nlp.model.meta.x0)
   xk, k, outdict = R2(
-    x -> obj(nlp, x),
-    (g, x) -> grad!(nlp, x, g),
+    x -> obj(reg_nlp.model, x),
+    (g, x) -> grad!(reg_nlp.model, x, g),
+    reg_nlp.h,
     args...,
     x0,
     nlp.meta.lvar,
     nlp.meta.uvar;
+    selected = reg_nlp.selected,
     kwargs_dict...,
   )
   ξ = outdict[:ξ]

From 6ed13e601c83b7ed3463814c9c988d0261311463 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Wed, 5 Jun 2024 11:17:59 +0200
Subject: [PATCH 02/24] Add RegularizedNLPModels to R2

---
 Project.toml                   |  3 ++-
 src/R2_alg.jl                  | 12 ++++++------
 src/RegularizedOptimization.jl |  2 +-
 3 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/Project.toml b/Project.toml
index 06e2ee59..d5ba1185 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,6 +11,7 @@ NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 NLPModelsModifiers = "e01155f1-5c6f-4375-a9d8-616dd036575f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537"
+RegularizedProblems = "ea076b23-609f-44d2-bb12-a4ae45328278"
 ShiftedProximalOperators = "d4fd37fa-580c-4e43-9b30-361c21aae263"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
 TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
@@ -20,7 +21,7 @@ LinearOperators = "2.7"
 NLPModels = "0.19, 0.20"
 NLPModelsModifiers = "0.7"
 ProximalOperators = "0.15"
-RegularizedProblems = "0.1"
+RegularizedProblems = "0.1.1"
 ShiftedProximalOperators = "0.2"
 SolverCore = "0.3.0"
 TSVD = "0.4"
diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 0d543814..23ab4c41 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -101,10 +101,10 @@ In the second form, instead of `nlp`, the user may pass in
 * `Complex_hist`: an array with the history of number of inner iterations.
 """
 
-function R2(nlp::AbstractNLPModel, args...; kwargs)
+function R2(nlp::AbstractNLPModel, h, args...; kwargs...)
   kwargs_dict = Dict(kwargs...)
-  selected = pop!(kwargs_dict, :selected, 1:nlp.meta.nvar)
-  reg_nlp = RegularizedNLPModel(nlp, args.h, selected)
+  selected = pop!(kwargs_dict, :selected, 1:nlp.meta.nvar) 
+  reg_nlp = RegularizedNLPModel(nlp, h, selected)
   return R2(reg_nlp, args...; kwargs...)
 end
 
@@ -117,13 +117,13 @@ function R2(reg_nlp::AbstractRegularizedNLPModel, args...; kwargs...)
     reg_nlp.h,
     args...,
     x0,
-    nlp.meta.lvar,
-    nlp.meta.uvar;
+    reg_nlp.model.meta.lvar,
+    reg_nlp.model.meta.uvar;
     selected = reg_nlp.selected,
     kwargs_dict...,
   )
   ξ = outdict[:ξ]
-  stats = GenericExecutionStats(nlp)
+  stats = GenericExecutionStats(reg_nlp.model)
   set_status!(stats, outdict[:status])
   set_solution!(stats, xk)
   set_objective!(stats, outdict[:fk] + outdict[:hk])
diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl
index 458a8cce..e92f21da 100644
--- a/src/RegularizedOptimization.jl
+++ b/src/RegularizedOptimization.jl
@@ -7,7 +7,7 @@ using LinearAlgebra, Logging, Printf
 using ProximalOperators, TSVD
 
 # dependencies from us
-using LinearOperators, NLPModels, NLPModelsModifiers, ShiftedProximalOperators, SolverCore
+using LinearOperators, NLPModels, NLPModelsModifiers, RegularizedProblems, ShiftedProximalOperators, SolverCore
 
 include("utils.jl")
 include("input_struct.jl")

From 435838a7e832ed014b2e4ccb8139eb047d8ff7b0 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Wed, 5 Jun 2024 11:45:11 +0200
Subject: [PATCH 03/24] Update documentation

---
 src/R2_alg.jl | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 23ab4c41..a9951374 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -57,6 +57,7 @@ function R2Solver(
 end
 
 """
+    R2(reg_nlp, options)
     R2(nlp, h, options)
     R2(f, ∇f!, h, options, x0)
 
@@ -75,20 +76,20 @@ where φ(s ; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs is the Taylor linear approximation
 ψ(s; xₖ) = h(xₖ + s), ‖⋅‖ is a user-defined norm and σₖ > 0 is the regularization parameter.
 
 ### Arguments
-
-* `nlp::AbstractNLPModel`: a smooth optimization problem
-* `h`: a regularizer such as those defined in ProximalOperators
+* `reg_nlp::AbstractRegularizedNLPModel`: a smooth, regularized optimization problem: min f(x) + h(x)
 * `options::ROSolverOptions`: a structure containing algorithmic parameters
-* `x0::AbstractVector`: an initial guess (in the second calling form)
+* `nlp::AbstractNLPModel`: a smooth optimization problem (in the second calling form)
+* `h`: a regularizer such as those defined in ProximalOperators (in the second calling form)
+* `x0::AbstractVector`: an initial guess (in the third calling form)
 
 ### Keyword Arguments
 
-* `x0::AbstractVector`: an initial guess (in the first calling form: default = `nlp.meta.x0`)
-* `selected::AbstractVector{<:Integer}`: (default `1:length(x0)`).
+* `x0::AbstractVector`: an initial guess (in the first and second calling form resp.: default = `reg_nlp.model.meta.x0` and `nlp.meta.x0`)
+* `selected::AbstractVector{<:Integer}`: (default `1:length(x0)`) (in the first calling form, this should be stored in reg_nlp)
 
-The objective and gradient of `nlp` will be accessed.
+The objective and gradient of the smooth part will be accessed.
 
-In the second form, instead of `nlp`, the user may pass in
+In the third form, instead of `nlp`, the user may pass in
 
 * `f` a function such that `f(x)` returns the value of f at x
 * `∇f!` a function to evaluate the gradient in place, i.e., such that `∇f!(g, x)` store ∇f(x) in `g`.

From 9e65f8279b4283f2c5e84bbe174dcbcb810181cb Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Sat, 8 Jun 2024 14:42:35 +0200
Subject: [PATCH 04/24] Change R2_alg with SolverCore.solve!

---
 src/R2_alg.jl | 487 ++++++++++++++++++++++++++++++++------------------
 1 file changed, 311 insertions(+), 176 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index a9951374..90fa9ee4 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -1,9 +1,10 @@
-export R2
+export R2,R2Solver,solve!
 
-mutable struct R2Solver{R, S <: AbstractVector{R}} <: AbstractOptimizationSolver
+mutable struct R2Solver{R <: Real,G <: Union{ShiftedProximableFunction, Nothing}, S <: AbstractVector{R}} <: AbstractOptimizationSolver
   xk::S
   ∇fk::S
   mν∇fk::S
+  ψ::G
   xkn::S
   s::S
   has_bnds::Bool
@@ -20,7 +21,8 @@ function R2Solver(
   x0::S,
   options::ROSolverOptions,
   l_bound::S,
-  u_bound::S,
+  u_bound::S;
+  ψ = nothing
 ) where {R <: Real, S <: AbstractVector{R}}
   maxIter = options.maxIter
   xk = similar(x0)
@@ -43,6 +45,7 @@ function R2Solver(
     xk,
     ∇fk,
     mν∇fk,
+    ψ,
     xkn,
     s,
     has_bnds,
@@ -56,6 +59,53 @@ function R2Solver(
   )
 end
 
+function R2Solver(
+  reg_nlp::AbstractRegularizedNLPModel{T,V};
+  max_iter::Int = 10000,
+  kwargs...
+) where {T,V}
+  kwargs_dict = Dict(kwargs...)
+  x0 = pop!(kwargs_dict, :x0, reg_nlp.model.meta.x0)
+  l_bound = reg_nlp.model.meta.lvar
+  u_bound = reg_nlp.model.meta.uvar
+
+  xk = similar(x0)
+  ∇fk = similar(x0)
+  mν∇fk = similar(x0)
+  xkn = similar(x0)
+  s = zero(x0)
+  has_bnds = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf))
+  if has_bnds
+    l_bound_m_x = similar(xk)
+    u_bound_m_x = similar(xk)
+  else
+    l_bound_m_x = similar(xk, 0)
+    u_bound_m_x = similar(xk, 0)
+  end
+  Fobj_hist = zeros(T, max_iter)
+  Hobj_hist = zeros(T, max_iter)
+  Complex_hist = zeros(Int, max_iter)
+  
+  ψ = has_bnds ? shifted(reg_nlp.h, x0, l_bound_m_x, u_bound_m_x, reg_nlp.selected) : shifted(reg_nlp.h, x0)
+  return R2Solver(
+    xk,
+    ∇fk,
+    mν∇fk,
+    ψ,
+    xkn,
+    s,
+    has_bnds,
+    l_bound,
+    u_bound,
+    l_bound_m_x,
+    u_bound_m_x,
+    Fobj_hist,
+    Hobj_hist,
+    Complex_hist,
+  )
+end
+
+
 """
     R2(reg_nlp, options)
     R2(nlp, h, options)
@@ -102,43 +152,33 @@ In the third form, instead of `nlp`, the user may pass in
 * `Complex_hist`: an array with the history of number of inner iterations.
 """
 
-function R2(nlp::AbstractNLPModel, h, args...; kwargs...)
+function R2(
+  nlp::AbstractNLPModel{R, V}, 
+  h, 
+  options::ROSolverOptions{R};
+  kwargs...) where{ R <: Real, V }
   kwargs_dict = Dict(kwargs...)
   selected = pop!(kwargs_dict, :selected, 1:nlp.meta.nvar) 
+  x0 = pop!(kwargs_dict, :x0, nlp.meta.x0)
   reg_nlp = RegularizedNLPModel(nlp, h, selected)
-  return R2(reg_nlp, args...; kwargs...)
-end
-
-function R2(reg_nlp::AbstractRegularizedNLPModel, args...; kwargs...)
-  kwargs_dict = Dict(kwargs...)
-  x0 = pop!(kwargs_dict, :x0, reg_nlp.model.meta.x0)
-  xk, k, outdict = R2(
-    x -> obj(reg_nlp.model, x),
-    (g, x) -> grad!(reg_nlp.model, x, g),
-    reg_nlp.h,
-    args...,
-    x0,
-    reg_nlp.model.meta.lvar,
-    reg_nlp.model.meta.uvar;
-    selected = reg_nlp.selected,
-    kwargs_dict...,
+  return R2(
+    reg_nlp,
+    x = x0,
+    a_tol = options.ϵa,
+    r_tol = options.ϵr,
+    neg_tol = options.neg_tol,
+    verbose = options.verbose,
+    max_iter = options.maxIter,
+    max_time = options.maxTime,
+    σmin = options.σmin,
+    η1 = options.η1,
+    η2 = options.η2,
+    ν = options.ν,
+    γ = options.γ,
   )
-  ξ = outdict[:ξ]
-  stats = GenericExecutionStats(reg_nlp.model)
-  set_status!(stats, outdict[:status])
-  set_solution!(stats, xk)
-  set_objective!(stats, outdict[:fk] + outdict[:hk])
-  set_residuals!(stats, zero(eltype(xk)), ξ)
-  set_iter!(stats, k)
-  set_time!(stats, outdict[:elapsed_time])
-  set_solver_specific!(stats, :Fhist, outdict[:Fhist])
-  set_solver_specific!(stats, :Hhist, outdict[:Hhist])
-  set_solver_specific!(stats, :NonSmooth, outdict[:NonSmooth])
-  set_solver_specific!(stats, :SubsolverCounter, outdict[:Chist])
   return stats
 end
 
-# method without bounds
 function R2(
   f::F,
   ∇f!::G,
@@ -148,23 +188,35 @@ function R2(
   selected::AbstractVector{<:Integer} = 1:length(x0),
   kwargs...,
 ) where {F <: Function, G <: Function, H, R <: Real}
-  start_time = time()
-  elapsed_time = 0.0
-  solver = R2Solver(x0, options, similar(x0, 0), similar(x0, 0))
-  k, status, fk, hk, ξ = R2!(solver, f, ∇f!, h, options, x0; selected = selected)
-  elapsed_time = time() - start_time
-  outdict = Dict(
-    :Fhist => solver.Fobj_hist[1:k],
-    :Hhist => solver.Hobj_hist[1:k],
-    :Chist => solver.Complex_hist[1:k],
+  nlp = FirstOrderModel(f,∇f!,x0)
+  reg_nlp = RegularizedNLPModel(nlp,h,selected) 
+  stats = R2(
+  reg_nlp,
+  x=x0,
+  a_tol = options.ϵa,
+  r_tol = options.ϵr,
+  neg_tol = options.neg_tol,
+  verbose = options.verbose,
+  max_iter = options.maxIter,
+  max_time = options.maxTime,
+  σmin = options.σmin,
+  η1 = options.η1,
+  η2 = options.η2,
+  ν = options.ν,
+  γ = options.γ,
+)
+outdict = Dict(
+    :Fhist => stats.solver_specific[:Fhist],
+    :Hhist => stats.solver_specific[:Hhist],
+    :Chist => stats.solver_specific[:SubsolverCounter],
     :NonSmooth => h,
-    :status => status,
-    :fk => fk,
-    :hk => hk,
-    :ξ => ξ,
-    :elapsed_time => elapsed_time,
+    :status => stats.status,
+    :fk => stats.solver_specific[:smooth_obj],
+    :hk => stats.solver_specific[:nonsmooth_obj],
+    :ξ => stats.solver_specific[:xi],
+    :elapsed_time => stats.elapsed_time,
   )
-  return solver.xk, k, outdict
+return stats.solution,stats.iter,outdict
 end
 
 function R2(
@@ -178,53 +230,93 @@ function R2(
   selected::AbstractVector{<:Integer} = 1:length(x0),
   kwargs...,
 ) where {F <: Function, G <: Function, H, R <: Real}
-  start_time = time()
-  elapsed_time = 0.0
-  solver = R2Solver(x0, options, l_bound, u_bound)
-  k, status, fk, hk, ξ = R2!(solver, f, ∇f!, h, options, x0; selected = selected)
-  elapsed_time = time() - start_time
+  nlp = FirstOrderModel(f,∇f!,x0,lcon = l_bound, ucon = u_bound)
+  reg_nlp = RegularizedNLPModel(nlp,h,selected) 
+  stats = R2(
+    reg_nlp,
+    x=x0,
+    a_tol = options.ϵa,
+    r_tol = options.ϵr,
+    neg_tol = options.neg_tol,
+    verbose = options.verbose,
+    max_iter = options.maxIter,
+    max_time = options.maxTime,
+    σmin = options.σmin,
+    η1 = options.η1,
+    η2 = options.η2,
+    ν = options.ν,
+    γ = options.γ,
+  )
   outdict = Dict(
-    :Fhist => solver.Fobj_hist[1:k],
-    :Hhist => solver.Hobj_hist[1:k],
-    :Chist => solver.Complex_hist[1:k],
+    :Fhist => stats.solver_specific[:Fhist],
+    :Hhist => stats.solver_specific[:Hhist],
+    :Chist => stats.solver_specific[:SubsolverCounter],
     :NonSmooth => h,
-    :status => status,
-    :fk => fk,
-    :hk => hk,
-    :ξ => ξ,
-    :elapsed_time => elapsed_time,
+    :status => stats.status,
+    :fk => stats.solver_specific[:smooth_obj],
+    :hk => stats.solver_specific[:nonsmooth_obj],
+    :ξ => stats.solver_specific[:xi],
+    :elapsed_time => stats.elapsed_time,
   )
-  return solver.xk, k, outdict
+  return stats.solution,stats.iter,outdict
 end
 
-function R2!(
-  solver::R2Solver{R},
-  f::F,
-  ∇f!::G,
-  h::H,
-  options::ROSolverOptions{R},
-  x0::AbstractVector{R};
-  selected::AbstractVector{<:Integer} = 1:length(x0),
-) where {F <: Function, G <: Function, H, R <: Real}
-  start_time = time()
-  elapsed_time = 0.0
-  ϵ = options.ϵa
-  ϵr = options.ϵr
-  neg_tol = options.neg_tol
-  verbose = options.verbose
-  maxIter = options.maxIter
-  maxTime = options.maxTime
-  σmin = options.σmin
-  η1 = options.η1
-  η2 = options.η2
-  ν = options.ν
-  γ = options.γ
-
-  # retrieve workspace
-  xk = solver.xk
-  xk .= x0
+
+function R2(reg_nlp::AbstractRegularizedNLPModel; kwargs...)
+  kwargs_dict = Dict(kwargs...)
+  x0 = pop!(kwargs_dict, :x, reg_nlp.model.meta.x0)
+  solver = R2Solver(reg_nlp,x0=x0)
+  stats = GenericExecutionStats(reg_nlp.model)
+  cb = (nlp, solver, stats) -> begin
+    solver.Fobj_hist[stats.iter+1] = stats.solver_specific[:smooth_obj]
+    solver.Hobj_hist[stats.iter+1] = stats.solver_specific[:nonsmooth_obj]
+    solver.Complex_hist[stats.iter+1] += 1
+  end
+  solve!(
+    solver,
+    reg_nlp,
+    stats;
+    x = x0,
+    callback = cb,
+    kwargs...
+  )
+  set_solver_specific!(stats, :Fhist, solver.Fobj_hist[1:stats.iter+1])
+  set_solver_specific!(stats, :Hhist, solver.Hobj_hist[1:stats.iter+1])
+  set_solver_specific!(stats, :SubsolverCounter, solver.Complex_hist[1:stats.iter+1])
+  return stats
+end
+
+function SolverCore.solve!(
+  solver::R2Solver{T}, 
+  reg_nlp::AbstractRegularizedNLPModel{T, V}, 
+  stats::GenericExecutionStats{T, V};
+  callback = (args...) -> nothing,
+  x::V = reg_nlp.model.meta.x0,
+  a_tol::T = √eps(T),
+  r_tol::T = √eps(T),
+  neg_tol::T = eps(T)^(1 / 4),
+  verbose::Int = 0,
+  max_iter::Int = 10000,
+  max_time::Float64 = 30.0,
+  max_eval::Int = -1,
+  σmin::T = eps(T),
+  η1::T = √√eps(T),
+  η2::T = T(0.9),
+  ν::T = eps(T)^(1 / 5),
+  γ::T = T(3),
+  ) where {T, V}
+
+  reset!(stats)
+
+  # Retrieve workspace
+  selected = reg_nlp.selected
+  h = reg_nlp.h
+  nlp = reg_nlp.model
+
+  xk = solver.xk .= x
   ∇fk = solver.∇fk
   mν∇fk = solver.mν∇fk
+  ψ = solver.ψ
   xkn = solver.xkn
   s = solver.s
   has_bnds = solver.has_bnds
@@ -234,104 +326,102 @@ function R2!(
     l_bound_m_x = solver.l_bound_m_x
     u_bound_m_x = solver.u_bound_m_x
   end
-  Fobj_hist = solver.Fobj_hist
-  Hobj_hist = solver.Hobj_hist
-  Complex_hist = solver.Complex_hist
-
-  if verbose == 0
-    ptf = Inf
-  elseif verbose == 1
-    ptf = round(maxIter / 10)
-  elseif verbose == 2
-    ptf = round(maxIter / 100)
-  else
-    ptf = 1
-  end
 
   # initialize parameters
+  improper = false
   hk = @views h(xk[selected])
   if hk == Inf
     verbose > 0 && @info "R2: finding initial guess where nonsmooth term is finite"
-    prox!(xk, h, x0, one(eltype(x0)))
+    prox!(xk, h, xk, one(eltype(x0)))
     hk = @views h(xk[selected])
     hk < Inf || error("prox computation must be erroneous")
     verbose > 0 && @debug "R2: found point where h has value" hk
   end
-  hk == -Inf && error("nonsmooth term is not proper")
+   improper = (hk == -Inf)
 
   if has_bnds
     @. l_bound_m_x = l_bound - xk
     @. u_bound_m_x = u_bound - xk
-    ψ = shifted(h, xk, l_bound_m_x, u_bound_m_x, selected)
-  else
-    ψ = shifted(h, xk)
   end
 
   if verbose > 0
-    #! format: off
-    @info @sprintf "%6s %8s %8s %7s %8s %7s %7s %7s %1s" "iter" "f(x)" "h(x)" "√(ξ/ν)" "ρ" "σ" "‖x‖" "‖s‖" ""
-    #! format: off
+    @info log_header(
+      [:iter, :fx, :hx, :xi, :ρ, :σ, :normx, :norms, :arrow],
+      [Int, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Char],
+      hdr_override = Dict{Symbol,String}(   # TODO: Add this as constant dict elsewhere
+        :iter => "iter",
+        :fx => "f(x)",
+        :hx => "h(x)",
+        :xi => "√(ξ/ν)",
+        :ρ => "ρ",
+        :σ => "σ",
+        :normx => "‖x‖",
+        :norms => "‖s‖",
+        :arrow => " "
+      ),
+      colsep = 1,
+    )
   end
 
-  local ξ::R
-  k = 0
+  local ξ::T
   σk = max(1 / ν, σmin)
   ν = 1 / σk
-  sqrt_ξ_νInv = one(R)
+  sqrt_ξ_νInv = one(T)
 
-  fk = f(xk)
-  ∇f!(∇fk, xk)
+  fk = obj(nlp, xk)
+  grad!(nlp, xk, ∇fk)
   @. mν∇fk = -ν * ∇fk
 
-  optimal = false
-  tired = maxIter > 0 && k ≥ maxIter || elapsed_time > maxTime
-
-  while !(optimal || tired)
-    k = k + 1
-    elapsed_time = time() - start_time
-    Fobj_hist[k] = fk
-    Hobj_hist[k] = hk
+  set_iter!(stats, 0)
+  start_time = time()
+  set_time!(stats, 0.0)
+  set_objective!(stats,fk + hk)
+  set_solver_specific!(stats,:smooth_obj,fk)
+  set_solver_specific!(stats,:nonsmooth_obj, hk)
+
+  φk(d) = dot(∇fk, d)   
+  mk(d)::T = φk(d) + ψ(d)::T
+
+  prox!(s, ψ, mν∇fk, ν)
+  mks = mk(s)
+
+  ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
+  ξ > 0 || error("R2: prox-gradient step should produce a decrease but ξ = $(ξ)")
+  sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
+  a_tol += r_tol * sqrt_ξ_νInv # make stopping test absolute and relative
+
+  solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ a_tol)
+
+  set_solver_specific!(stats,:xi,sqrt_ξ_νInv)
+  set_status!(
+    stats,
+    get_status(
+      reg_nlp,
+      elapsed_time = stats.elapsed_time,
+      iter = stats.iter,
+      optimal = solved,
+      improper = improper,
+      max_eval = max_eval,
+      max_time = max_time,
+      max_iter = max_iter
+    ),
+  )
 
-    # define model
-    φk(d) = dot(∇fk, d)
-    mk(d)::R = φk(d) + ψ(d)::R
+  callback(nlp, solver, stats)
 
-    prox!(s, ψ, mν∇fk, ν)
-    Complex_hist[k] += 1
-    mks = mk(s)
-    ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
-    sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
+  done = stats.status != :unknown
 
-    if ξ ≥ 0 && k == 1
-      ϵ += ϵr * sqrt_ξ_νInv # make stopping test absolute and relative
-    end
+  while !done
     
-    if (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ ϵ)
-      optimal = true
-      continue
-    end
-
-    ξ > 0 || error("R2: prox-gradient step should produce a decrease but ξ = $(ξ)")
+    # Update xk, sigma_k
     xkn .= xk .+ s
-    fkn = f(xkn)
+    fkn = obj(nlp, xkn)
     hkn = @views h(xkn[selected])
-    hkn == -Inf && error("nonsmooth term is not proper")
+    improper = (hkn == -Inf)
 
     Δobj = (fk + hk) - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
     ρk = Δobj / ξ
 
-    if (verbose > 0) && (k % ptf == 0)
-      #! format: off
-      σ_stat = (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")
-      @info @sprintf "%6d %8.1e %8.1e %7.1e %8.1e %7.1e %7.1e %7.1e %1s" k fk hk sqrt_ξ_νInv ρk σk norm(xk) norm(s) σ_stat
-
-      #! format: on
-    end
-
-    if η2 ≤ ρk < Inf
-      σk = max(σk / γ, σmin)
-    end
-
     if η1 ≤ ρk < Inf
       xk .= xkn
       if has_bnds
@@ -341,41 +431,86 @@ function R2!(
       end
       fk = fkn
       hk = hkn
-      ∇f!(∇fk, xk)
+      grad!(nlp, xk, ∇fk)
       shift!(ψ, xk)
     end
 
+    if η2 ≤ ρk < Inf
+      σk = max(σk / γ, σmin)
+    end
     if ρk < η1 || ρk == Inf
       σk = σk * γ
     end
 
     ν = 1 / σk
-    tired = maxIter > 0 && k ≥ maxIter
-    if !tired
-      @. mν∇fk = -ν * ∇fk
-    end
-  end
+    @. mν∇fk = -ν * ∇fk
+  
+    set_objective!(stats, fk + hk)
+    set_solver_specific!(stats,:smooth_obj,fk)
+    set_solver_specific!(stats,:nonsmooth_obj, hk)
+    set_iter!(stats, stats.iter + 1)
+    set_time!(stats, time() - start_time)
 
-  if verbose > 0
-    if k == 1
-      @info @sprintf "%6d %8.1e %8.1e" k fk hk
-    elseif optimal
-      #! format: off
-      @info @sprintf "%6d %8.1e %8.1e %7.1e %8s %7.1e %7.1e %7.1e" k fk hk sqrt(ξ/ν) "" σk norm(xk) norm(s)
-      #! format: on
-      @info "R2: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)"
-    end
+    prox!(s, ψ, mν∇fk, ν)
+    mks = mk(s)
+
+    ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
+    sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
+    solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ a_tol)
+  
+    verbose > 0 && 
+      stats.iter % verbose == 0 &&
+        @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
+
+    set_solver_specific!(stats,:xi,sqrt_ξ_νInv)
+    set_status!(
+      stats,
+      get_status(
+        reg_nlp,
+        elapsed_time = stats.elapsed_time,
+        iter = stats.iter,
+        optimal = solved,
+        improper = improper,
+        max_eval = max_eval,
+        max_time = max_time,
+        max_iter = max_iter
+      ),
+    )
+
+    callback(nlp, solver, stats)
+
+    done = stats.status != :unknown
   end
 
-  status = if optimal
+  verbose > 0 &&
+    stats.status == :first_order &&
+        @info "R2: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)"
+
+  set_solution!(stats,xk)
+  return stats
+end
+
+function get_status(
+  reg_nlp;
+  elapsed_time = 0.0,
+  iter = 0,
+  optimal = false,
+  improper = false,
+  max_eval = Inf,
+  max_time = Inf,
+  max_iter = Inf,
+)
+  if optimal
     :first_order
-  elseif elapsed_time > maxTime
-    :max_time
-  elseif tired
+  elseif improper
+    :improper
+  elseif iter > max_iter
     :max_iter
+  elseif elapsed_time > max_time
+    :max_time
+  elseif neval_obj(reg_nlp.model) > max_eval && max_eval != -1
+    :max_eval
   else
-    :exception
+    :unknown
   end
-
-  return k, status, fk, hk, sqrt_ξ_νInv
-end
+end
\ No newline at end of file

From 5895a858160bebab36eac30f0083456adc53bbe0 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Sat, 8 Jun 2024 14:51:39 +0200
Subject: [PATCH 05/24] Fix bound issue, all tests pass

---
 src/R2_alg.jl | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 90fa9ee4..f70e1fbb 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -78,6 +78,8 @@ function R2Solver(
   if has_bnds
     l_bound_m_x = similar(xk)
     u_bound_m_x = similar(xk)
+    @. l_bound_m_x = l_bound - x0
+    @. u_bound_m_x = u_bound - x0
   else
     l_bound_m_x = similar(xk, 0)
     u_bound_m_x = similar(xk, 0)
@@ -337,12 +339,7 @@ function SolverCore.solve!(
     hk < Inf || error("prox computation must be erroneous")
     verbose > 0 && @debug "R2: found point where h has value" hk
   end
-   improper = (hk == -Inf)
-
-  if has_bnds
-    @. l_bound_m_x = l_bound - xk
-    @. u_bound_m_x = u_bound - xk
-  end
+  improper = (hk == -Inf)
 
   if verbose > 0
     @info log_header(

From 9c5217b30c352f110291ce0c60e794632e86ad04 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Sat, 8 Jun 2024 16:06:16 +0200
Subject: [PATCH 06/24] Update documentation

---
 src/R2_alg.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index f70e1fbb..af8621e7 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -109,9 +109,7 @@ end
 
 
 """
-    R2(reg_nlp, options)
-    R2(nlp, h, options)
-    R2(f, ∇f!, h, options, x0)
+    R2(reg_nlp)
 
 A first-order quadratic regularization method for the problem
 
@@ -127,6 +125,8 @@ About each iterate xₖ, a step sₖ is computed as a solution of
 where φ(s ; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs is the Taylor linear approximation of f about xₖ,
 ψ(s; xₖ) = h(xₖ + s), ‖⋅‖ is a user-defined norm and σₖ > 0 is the regularization parameter.
 
+For advanced usage, first define a solver "R2Solver" to preallocate the memory used in the algorithm, and then call `solve!`:
+
 ### Arguments
 * `reg_nlp::AbstractRegularizedNLPModel`: a smooth, regularized optimization problem: min f(x) + h(x)
 * `options::ROSolverOptions`: a structure containing algorithmic parameters

From 3ee28f6c7b39e969a76442be949cab2b7d4d2029 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Sat, 8 Jun 2024 16:31:22 +0200
Subject: [PATCH 07/24] Minor bug fixes

---
 src/R2_alg.jl | 110 ++++++++++++++++++++++++++++++--------------------
 1 file changed, 66 insertions(+), 44 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index af8621e7..4199c6b7 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -1,5 +1,7 @@
 export R2,R2Solver,solve!
 
+import SolverCore.solve!
+
 mutable struct R2Solver{R <: Real,G <: Union{ShiftedProximableFunction, Nothing}, S <: AbstractVector{R}} <: AbstractOptimizationSolver
   xk::S
   ∇fk::S
@@ -61,11 +63,9 @@ end
 
 function R2Solver(
   reg_nlp::AbstractRegularizedNLPModel{T,V};
-  max_iter::Int = 10000,
-  kwargs...
-) where {T,V}
-  kwargs_dict = Dict(kwargs...)
-  x0 = pop!(kwargs_dict, :x0, reg_nlp.model.meta.x0)
+  max_iter::Int = 10000
+  ) where {T,V}
+  x0 = reg_nlp.model.meta.x0
   l_bound = reg_nlp.model.meta.lvar
   u_bound = reg_nlp.model.meta.uvar
 
@@ -127,31 +127,52 @@ where φ(s ; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs is the Taylor linear approximation
 
 For advanced usage, first define a solver "R2Solver" to preallocate the memory used in the algorithm, and then call `solve!`:
 
-### Arguments
-* `reg_nlp::AbstractRegularizedNLPModel`: a smooth, regularized optimization problem: min f(x) + h(x)
-* `options::ROSolverOptions`: a structure containing algorithmic parameters
-* `nlp::AbstractNLPModel`: a smooth optimization problem (in the second calling form)
-* `h`: a regularizer such as those defined in ProximalOperators (in the second calling form)
-* `x0::AbstractVector`: an initial guess (in the third calling form)
-
-### Keyword Arguments
-
-* `x0::AbstractVector`: an initial guess (in the first and second calling form resp.: default = `reg_nlp.model.meta.x0` and `nlp.meta.x0`)
-* `selected::AbstractVector{<:Integer}`: (default `1:length(x0)`) (in the first calling form, this should be stored in reg_nlp)
-
-The objective and gradient of the smooth part will be accessed.
-
-In the third form, instead of `nlp`, the user may pass in
+    solver = R2Solver(reg_nlp)
+    solve!(solver, reg_nlp)
 
-* `f` a function such that `f(x)` returns the value of f at x
-* `∇f!` a function to evaluate the gradient in place, i.e., such that `∇f!(g, x)` store ∇f(x) in `g`.
+    stats = GenericExecutionStats(reg_nlp.model)
+    solver = R2Solver(reg_nlp)
+    solve!(solver, reg_nlp, stats)
 
-### Return values
-
-* `xk`: the final iterate
-* `Fobj_hist`: an array with the history of values of the smooth objective
-* `Hobj_hist`: an array with the history of values of the nonsmooth objective
-* `Complex_hist`: an array with the history of number of inner iterations.
+### Arguments
+* `reg_nlp::AbstractRegularizedNLPModel{T, V}`: the problem to solve, see `RegularizedProblems.jl`, `NLPModels.jl`.
+
+# Keyword arguments 
+- `x::V = nlp.meta.x0`: the initial guess;
+- `atol::T = √eps(T)`: absolute tolerance;
+- `rtol::T = √eps(T)`: relative tolerance;
+- `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance
+- `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited);
+- `max_time::Float64 = 30.0`: maximum time limit in seconds;
+- `max_iter::Int = 10000`: maximum number of iterations;
+- `verbose::Int = 0`: if > 0, display iteration details every `verbose` iteration;
+- `σmin::T = eps(T)`: minimum value of the regularization parameter;
+- `η1::T = √√eps(T)`: very successful iteration threshold;
+- `η2::T = T(0.9)`: successful iteration threshold;
+- `ν::T = eps(T)^(1 / 5)`: multiplicative inverse of the regularization parameter: ν = 1/σ    ;
+- `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful.
+
+The algorithm stops when ``√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) `` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ)
+
+# Output
+The value returned is a `GenericExecutionStats`, see `SolverCore.jl`.
+
+# Callback
+The callback is called at each iteration.
+The expected signature of the callback is `callback(nlp, solver, stats)`, and its output is ignored.
+Changing any of the input arguments will affect the subsequent iterations.
+In particular, setting `stats.status = :user` will stop the algorithm.
+All relevant information should be available in `nlp` and `solver`.
+Notably, you can access, and modify, the following:
+- `solver.xk`: current iterate;
+- `solver.∇fk`: current gradient;
+- `stats`: structure holding the output of the algorithm (`GenericExecutionStats`), which contains, among other things:
+  - `stats.iter`: current iteration counter;
+  - `stats.objective`: current objective function value;
+  - `stats.solver_specific[:smooth_obj]`: current smooth part of the objective function
+  - `stats.solver_specific[:nonsmooth_objoth]`: current nonsmooth part of the objective function
+  - `stats.status`: current status of the algorithm. Should be `:unknown` unless the algorithm has attained a stopping criterion. Changing this to anything will stop the algorithm, but you should use `:user` to properly indicate the intention.
+  - `stats.elapsed_time`: elapsed time in seconds.
 """
 
 function R2(
@@ -166,8 +187,8 @@ function R2(
   return R2(
     reg_nlp,
     x = x0,
-    a_tol = options.ϵa,
-    r_tol = options.ϵr,
+    atol = options.ϵa,
+    rtol = options.ϵr,
     neg_tol = options.neg_tol,
     verbose = options.verbose,
     max_iter = options.maxIter,
@@ -195,8 +216,8 @@ function R2(
   stats = R2(
   reg_nlp,
   x=x0,
-  a_tol = options.ϵa,
-  r_tol = options.ϵr,
+  atol = options.ϵa,
+  rtol = options.ϵr,
   neg_tol = options.neg_tol,
   verbose = options.verbose,
   max_iter = options.maxIter,
@@ -237,8 +258,8 @@ function R2(
   stats = R2(
     reg_nlp,
     x=x0,
-    a_tol = options.ϵa,
-    r_tol = options.ϵr,
+    atol = options.ϵa,
+    rtol = options.ϵr,
     neg_tol = options.neg_tol,
     verbose = options.verbose,
     max_iter = options.maxIter,
@@ -265,9 +286,7 @@ end
 
 
 function R2(reg_nlp::AbstractRegularizedNLPModel; kwargs...)
-  kwargs_dict = Dict(kwargs...)
-  x0 = pop!(kwargs_dict, :x, reg_nlp.model.meta.x0)
-  solver = R2Solver(reg_nlp,x0=x0)
+  solver = R2Solver(reg_nlp)
   stats = GenericExecutionStats(reg_nlp.model)
   cb = (nlp, solver, stats) -> begin
     solver.Fobj_hist[stats.iter+1] = stats.solver_specific[:smooth_obj]
@@ -278,7 +297,6 @@ function R2(reg_nlp::AbstractRegularizedNLPModel; kwargs...)
     solver,
     reg_nlp,
     stats;
-    x = x0,
     callback = cb,
     kwargs...
   )
@@ -294,8 +312,8 @@ function SolverCore.solve!(
   stats::GenericExecutionStats{T, V};
   callback = (args...) -> nothing,
   x::V = reg_nlp.model.meta.x0,
-  a_tol::T = √eps(T),
-  r_tol::T = √eps(T),
+  atol::T = √eps(T),
+  rtol::T = √eps(T),
   neg_tol::T = eps(T)^(1 / 4),
   verbose::Int = 0,
   max_iter::Int = 10000,
@@ -314,7 +332,11 @@ function SolverCore.solve!(
   selected = reg_nlp.selected
   h = reg_nlp.h
   nlp = reg_nlp.model
-
+  
+  # Make sure ψ is well shifted (on nlp.meta.x0 in R2Solver Constructor but might be different)
+  if solver.xk != x
+    shift!(solver.ψ,x)
+  end
   xk = solver.xk .= x
   ∇fk = solver.∇fk
   mν∇fk = solver.mν∇fk
@@ -385,9 +407,9 @@ function SolverCore.solve!(
   ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
   ξ > 0 || error("R2: prox-gradient step should produce a decrease but ξ = $(ξ)")
   sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
-  a_tol += r_tol * sqrt_ξ_νInv # make stopping test absolute and relative
+  atol += rtol * sqrt_ξ_νInv # make stopping test absolute and relative
 
-  solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ a_tol)
+  solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol)
 
   set_solver_specific!(stats,:xi,sqrt_ξ_νInv)
   set_status!(
@@ -453,7 +475,7 @@ function SolverCore.solve!(
 
     ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
     sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
-    solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ a_tol)
+    solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol)
   
     verbose > 0 && 
       stats.iter % verbose == 0 &&

From c0d1905189fe9d161e960fd9748da2cabd629b33 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Sat, 8 Jun 2024 22:06:31 +0200
Subject: [PATCH 08/24] Spacing and documentation update

---
 src/R2_alg.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 4199c6b7..7a9b48ce 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -1,8 +1,8 @@
-export R2,R2Solver,solve!
+export R2, R2Solver, solve!
 
 import SolverCore.solve!
 
-mutable struct R2Solver{R <: Real,G <: Union{ShiftedProximableFunction, Nothing}, S <: AbstractVector{R}} <: AbstractOptimizationSolver
+mutable struct R2Solver{R <: Real, G <: Union{ShiftedProximableFunction, Nothing}, S <: AbstractVector{R}} <: AbstractOptimizationSolver
   xk::S
   ∇fk::S
   mν∇fk::S
@@ -149,10 +149,10 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u
 - `σmin::T = eps(T)`: minimum value of the regularization parameter;
 - `η1::T = √√eps(T)`: very successful iteration threshold;
 - `η2::T = T(0.9)`: successful iteration threshold;
-- `ν::T = eps(T)^(1 / 5)`: multiplicative inverse of the regularization parameter: ν = 1/σ    ;
+- `ν::T = eps(T)^(1 / 5)`: multiplicative inverse of the regularization parameter: ν = 1/σ;
 - `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful.
 
-The algorithm stops when ``√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) `` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ)
+The algorithm stops either when `√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ).
 
 # Output
 The value returned is a `GenericExecutionStats`, see `SolverCore.jl`.

From 86d2271cc599ec51e3f5db584cf8585164ed3115 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 10 Jun 2024 15:13:31 +0200
Subject: [PATCH 09/24] Update documentation

---
 src/R2_alg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 7a9b48ce..b8065014 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -170,7 +170,7 @@ Notably, you can access, and modify, the following:
   - `stats.iter`: current iteration counter;
   - `stats.objective`: current objective function value;
   - `stats.solver_specific[:smooth_obj]`: current smooth part of the objective function
-  - `stats.solver_specific[:nonsmooth_objoth]`: current nonsmooth part of the objective function
+  - `stats.solver_specific[:nonsmooth_obj]`: current value of the nonsmooth part of the objective function
   - `stats.status`: current status of the algorithm. Should be `:unknown` unless the algorithm has attained a stopping criterion. Changing this to anything will stop the algorithm, but you should use `:user` to properly indicate the intention.
   - `stats.elapsed_time`: elapsed time in seconds.
 """

From 61d7d2064ef9a66aaf36912eecf7b20751edcb45 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Mon, 10 Jun 2024 15:12:46 +0200
Subject: [PATCH 10/24] Update src/R2_alg.jl

Co-authored-by: Dominique <dominique.orban@gmail.com>
---
 src/R2_alg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index b8065014..c6fcb28c 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -169,7 +169,7 @@ Notably, you can access, and modify, the following:
 - `stats`: structure holding the output of the algorithm (`GenericExecutionStats`), which contains, among other things:
   - `stats.iter`: current iteration counter;
   - `stats.objective`: current objective function value;
-  - `stats.solver_specific[:smooth_obj]`: current smooth part of the objective function
+  - `stats.solver_specific[:smooth_obj]`: current value of the smooth part of the objective function
   - `stats.solver_specific[:nonsmooth_obj]`: current value of the nonsmooth part of the objective function
   - `stats.status`: current status of the algorithm. Should be `:unknown` unless the algorithm has attained a stopping criterion. Changing this to anything will stop the algorithm, but you should use `:user` to properly indicate the intention.
   - `stats.elapsed_time`: elapsed time in seconds.

From 47247e7ba3347e9da077d98ebcb8acf8a61bfe0c Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Mon, 10 Jun 2024 15:13:05 +0200
Subject: [PATCH 11/24] Update src/R2_alg.jl

Co-authored-by: Dominique <dominique.orban@gmail.com>
---
 src/R2_alg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index c6fcb28c..6cf92aae 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -109,7 +109,7 @@ end
 
 
 """
-    R2(reg_nlp)
+    R2(reg_nlp; kwargs…)
 
 A first-order quadratic regularization method for the problem
 

From 4faf4edc640a6e3af2ef94a564143ae289fa8c36 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 10 Jun 2024 15:16:09 +0200
Subject: [PATCH 12/24] Update documentation

---
 src/R2_alg.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 6cf92aae..da097087 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -174,7 +174,6 @@ Notably, you can access, and modify, the following:
   - `stats.status`: current status of the algorithm. Should be `:unknown` unless the algorithm has attained a stopping criterion. Changing this to anything will stop the algorithm, but you should use `:user` to properly indicate the intention.
   - `stats.elapsed_time`: elapsed time in seconds.
 """
-
 function R2(
   nlp::AbstractNLPModel{R, V}, 
   h, 
@@ -334,7 +333,7 @@ function SolverCore.solve!(
   nlp = reg_nlp.model
   
   # Make sure ψ is well shifted (on nlp.meta.x0 in R2Solver Constructor but might be different)
-  if solver.xk != x
+  if !all(solver.xk .== x)
     shift!(solver.ψ,x)
   end
   xk = solver.xk .= x

From dbc1189d0681c7e9e4e24d36f60eb5f07bc1ff9d Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 10 Jun 2024 16:28:38 +0200
Subject: [PATCH 13/24] Update @info log_row position

---
 src/R2_alg.jl | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index da097087..a6dcdcee 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -425,6 +425,10 @@ function SolverCore.solve!(
     ),
   )
 
+  verbose > 0 && 
+      stats.iter % verbose == 0 &&
+        @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
+  
   callback(nlp, solver, stats)
 
   done = stats.status != :unknown
@@ -475,10 +479,6 @@ function SolverCore.solve!(
     ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
     sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
     solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol)
-  
-    verbose > 0 && 
-      stats.iter % verbose == 0 &&
-        @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
 
     set_solver_specific!(stats,:xi,sqrt_ξ_νInv)
     set_status!(
@@ -495,6 +495,11 @@ function SolverCore.solve!(
       ),
     )
 
+    verbose > 0 && 
+      stats.iter % verbose == 0 &&
+        @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
+    
+
     callback(nlp, solver, stats)
 
     done = stats.status != :unknown

From ebf1457106dfd4237d4d79aa2bf78f71397c78b2 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 10 Jun 2024 16:34:40 +0200
Subject: [PATCH 14/24] Update @info log_row position

---
 src/R2_alg.jl | 13 ++++---------
 1 file changed, 4 insertions(+), 9 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index a6dcdcee..03de1dfe 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -424,10 +424,6 @@ function SolverCore.solve!(
       max_iter = max_iter
     ),
   )
-
-  verbose > 0 && 
-      stats.iter % verbose == 0 &&
-        @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
   
   callback(nlp, solver, stats)
 
@@ -444,6 +440,10 @@ function SolverCore.solve!(
     Δobj = (fk + hk) - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
     ρk = Δobj / ξ
 
+    verbose > 0 && 
+    stats.iter % verbose == 0 &&
+      @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
+
     if η1 ≤ ρk < Inf
       xk .= xkn
       if has_bnds
@@ -495,11 +495,6 @@ function SolverCore.solve!(
       ),
     )
 
-    verbose > 0 && 
-      stats.iter % verbose == 0 &&
-        @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
-    
-
     callback(nlp, solver, stats)
 
     done = stats.status != :unknown

From ef3de43cb39fd9ade1d1ed8c8612b5fab033ea92 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 10 Jun 2024 16:28:38 +0200
Subject: [PATCH 15/24] Update @info log_row position

---
 src/R2_alg.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 03de1dfe..a75c1878 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -495,6 +495,11 @@ function SolverCore.solve!(
       ),
     )
 
+    verbose > 0 && 
+      stats.iter % verbose == 0 &&
+        @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
+    
+
     callback(nlp, solver, stats)
 
     done = stats.status != :unknown

From 3dcc66b71d585230281b80e2d1f94755cf2e0c3d Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 10 Jun 2024 16:34:40 +0200
Subject: [PATCH 16/24] Update @info log_row position

---
 src/R2_alg.jl | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index a75c1878..03de1dfe 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -495,11 +495,6 @@ function SolverCore.solve!(
       ),
     )
 
-    verbose > 0 && 
-      stats.iter % verbose == 0 &&
-        @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
-    
-
     callback(nlp, solver, stats)
 
     done = stats.status != :unknown

From 2c4b99290a272840ae754d778fd20846c5ba43dd Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Thu, 13 Jun 2024 13:50:50 +0200
Subject: [PATCH 17/24] Fix R2Solver maxIter constructor

---
 src/R2_alg.jl | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 03de1dfe..452f5f63 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -40,9 +40,9 @@ function R2Solver(
     l_bound_m_x = similar(xk, 0)
     u_bound_m_x = similar(xk, 0)
   end
-  Fobj_hist = zeros(R, maxIter)
-  Hobj_hist = zeros(R, maxIter)
-  Complex_hist = zeros(Int, maxIter)
+  Fobj_hist = zeros(R, maxIter+2)
+  Hobj_hist = zeros(R, maxIter+2)
+  Complex_hist = zeros(Int, maxIter+2)
   return R2Solver(
     xk,
     ∇fk,
@@ -84,9 +84,9 @@ function R2Solver(
     l_bound_m_x = similar(xk, 0)
     u_bound_m_x = similar(xk, 0)
   end
-  Fobj_hist = zeros(T, max_iter)
-  Hobj_hist = zeros(T, max_iter)
-  Complex_hist = zeros(Int, max_iter)
+  Fobj_hist = zeros(T, max_iter+2)
+  Hobj_hist = zeros(T, max_iter+2)
+  Complex_hist = zeros(Int, max_iter+2)
   
   ψ = has_bnds ? shifted(reg_nlp.h, x0, l_bound_m_x, u_bound_m_x, reg_nlp.selected) : shifted(reg_nlp.h, x0)
   return R2Solver(
@@ -285,7 +285,9 @@ end
 
 
 function R2(reg_nlp::AbstractRegularizedNLPModel; kwargs...)
-  solver = R2Solver(reg_nlp)
+  kwargs_dict = Dict(kwargs...)
+  max_iter = pop!(kwargs_dict, :max_iter, 10000) 
+  solver = R2Solver(reg_nlp,max_iter = max_iter)
   stats = GenericExecutionStats(reg_nlp.model)
   cb = (nlp, solver, stats) -> begin
     solver.Fobj_hist[stats.iter+1] = stats.solver_specific[:smooth_obj]
@@ -297,6 +299,7 @@ function R2(reg_nlp::AbstractRegularizedNLPModel; kwargs...)
     reg_nlp,
     stats;
     callback = cb,
+    max_iter = max_iter,
     kwargs...
   )
   set_solver_specific!(stats, :Fhist, solver.Fobj_hist[1:stats.iter+1])

From cff9f50d1a0158647a9fe03f95d72bdc5d4122fa Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Wed, 19 Jun 2024 10:37:55 +0200
Subject: [PATCH 18/24] R2 now prints last iteration pass with final
 stationnarity measure

---
 src/R2_alg.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 452f5f63..c3f5a97f 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -441,7 +441,7 @@ function SolverCore.solve!(
     improper = (hkn == -Inf)
 
     Δobj = (fk + hk) - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
-    ρk = Δobj / ξ
+    global ρk = Δobj / ξ
 
     verbose > 0 && 
     stats.iter % verbose == 0 &&
@@ -505,7 +505,8 @@ function SolverCore.solve!(
 
   verbose > 0 &&
     stats.status == :first_order &&
-        @info "R2: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)"
+      @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
+      @info "R2: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)"
 
   set_solution!(stats,xk)
   return stats

From b1ffc31c44d908f9ed618e6fc65d43c3923ee683 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Wed, 19 Jun 2024 21:59:19 +0200
Subject: [PATCH 19/24] Update src/R2_alg.jl

Co-authored-by: Dominique <dominique.orban@gmail.com>
---
 src/R2_alg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index c3f5a97f..72950c59 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -335,7 +335,7 @@ function SolverCore.solve!(
   h = reg_nlp.h
   nlp = reg_nlp.model
   
-  # Make sure ψ is well shifted (on nlp.meta.x0 in R2Solver Constructor but might be different)
+  # Make sure ψ has the correct shift (on nlp.meta.x0 in R2Solver Constructor but might be different)
   if !all(solver.xk .== x)
     shift!(solver.ψ,x)
   end

From 6dfcb5355b9591f0092994afe549f13effefe188 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Wed, 19 Jun 2024 21:59:30 +0200
Subject: [PATCH 20/24] Update src/R2_alg.jl

Co-authored-by: Dominique <dominique.orban@gmail.com>
---
 src/R2_alg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 72950c59..351bec63 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -152,7 +152,7 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u
 - `ν::T = eps(T)^(1 / 5)`: multiplicative inverse of the regularization parameter: ν = 1/σ;
 - `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful.
 
-The algorithm stops either when `√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ).
+The algorithm stops either when `√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is stationarity measure.
 
 # Output
 The value returned is a `GenericExecutionStats`, see `SolverCore.jl`.

From a4a77e9edd6abd724e31c7f6e0252a58b3670358 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Sat, 22 Jun 2024 17:19:41 +0200
Subject: [PATCH 21/24] Update src/R2_alg.jl

Co-authored-by: Dominique <dominique.orban@gmail.com>
---
 src/R2_alg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 351bec63..8110e574 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -134,7 +134,7 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u
     solver = R2Solver(reg_nlp)
     solve!(solver, reg_nlp, stats)
 
-### Arguments
+# Arguments
 * `reg_nlp::AbstractRegularizedNLPModel{T, V}`: the problem to solve, see `RegularizedProblems.jl`, `NLPModels.jl`.
 
 # Keyword arguments 

From f0a273cb1784e6832ba67dba14eda966922a9301 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Wed, 26 Jun 2024 14:01:08 +0200
Subject: [PATCH 22/24] Solve nlp.meta.x0 bug

---
 src/R2_alg.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 8110e574..afbfe787 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -88,7 +88,7 @@ function R2Solver(
   Hobj_hist = zeros(T, max_iter+2)
   Complex_hist = zeros(Int, max_iter+2)
   
-  ψ = has_bnds ? shifted(reg_nlp.h, x0, l_bound_m_x, u_bound_m_x, reg_nlp.selected) : shifted(reg_nlp.h, x0)
+  ψ = has_bnds ? shifted(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) : shifted(reg_nlp.h, xk)
   return R2Solver(
     xk,
     ∇fk,
@@ -335,11 +335,11 @@ function SolverCore.solve!(
   h = reg_nlp.h
   nlp = reg_nlp.model
   
-  # Make sure ψ has the correct shift (on nlp.meta.x0 in R2Solver Constructor but might be different)
-  if !all(solver.xk .== x)
-    shift!(solver.ψ,x)
-  end
   xk = solver.xk .= x
+
+  # Make sure ψ has the correct shift 
+  shift!(solver.ψ,xk)
+
   ∇fk = solver.∇fk
   mν∇fk = solver.mν∇fk
   ψ = solver.ψ

From 95b7106b0e578214cb35f4740d7ebd5992c9b105 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Wed, 3 Jul 2024 13:47:47 +0200
Subject: [PATCH 23/24] Minor changes in R2

---
 src/R2_alg.jl | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index afbfe787..3f494ad6 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -407,11 +407,12 @@ function SolverCore.solve!(
   mks = mk(s)
 
   ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
-  ξ > 0 || error("R2: prox-gradient step should produce a decrease but ξ = $(ξ)")
+
   sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
   atol += rtol * sqrt_ξ_νInv # make stopping test absolute and relative
 
   solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol)
+  (ξ < 0 && sqrt_ξ_νInv > neg_tol) && error("R2: prox-gradient step should produce a decrease but ξ = $(ξ)")
 
   set_solver_specific!(stats,:xi,sqrt_ξ_νInv)
   set_status!(
@@ -482,6 +483,7 @@ function SolverCore.solve!(
     ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
     sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
     solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol)
+    (ξ < 0 && sqrt_ξ_νInv > neg_tol) && error("R2: prox-gradient step should produce a decrease but ξ = $(ξ)")
 
     set_solver_specific!(stats,:xi,sqrt_ξ_νInv)
     set_status!(
@@ -503,17 +505,17 @@ function SolverCore.solve!(
     done = stats.status != :unknown
   end
 
-  verbose > 0 &&
-    stats.status == :first_order &&
-      @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
-      @info "R2: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)"
+  if verbose > 0 && stats.status == :first_order
+    @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")], colsep = 1)
+    @info "R2: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)"
+  end
 
   set_solution!(stats,xk)
   return stats
 end
 
 function get_status(
-  reg_nlp;
+  reg_nlp::M;
   elapsed_time = 0.0,
   iter = 0,
   optimal = false,
@@ -521,7 +523,7 @@ function get_status(
   max_eval = Inf,
   max_time = Inf,
   max_iter = Inf,
-)
+) where{ M <: AbstractRegularizedNLPModel }
   if optimal
     :first_order
   elseif improper
@@ -530,7 +532,7 @@ function get_status(
     :max_iter
   elseif elapsed_time > max_time
     :max_time
-  elseif neval_obj(reg_nlp.model) > max_eval && max_eval != -1
+  elseif neval_obj(reg_nlp.model) > max_eval && max_eval > -1
     :max_eval
   else
     :unknown

From d088f8c21187d07a8eb228bb0a94163910e25eaf Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Wed, 3 Jul 2024 15:43:11 +0200
Subject: [PATCH 24/24] Update documentation

---
 src/R2_alg.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/R2_alg.jl b/src/R2_alg.jl
index 3f494ad6..ede95e9c 100644
--- a/src/R2_alg.jl
+++ b/src/R2_alg.jl
@@ -123,7 +123,7 @@ About each iterate xₖ, a step sₖ is computed as a solution of
     min  φ(s; xₖ) + ½ σₖ ‖s‖² + ψ(s; xₖ)
 
 where φ(s ; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs is the Taylor linear approximation of f about xₖ,
-ψ(s; xₖ) = h(xₖ + s), ‖⋅‖ is a user-defined norm and σₖ > 0 is the regularization parameter.
+ψ(s; xₖ) is either h(xₖ + s) or an approximation of h(xₖ + s), ‖⋅‖ is a user-defined norm and σₖ > 0 is the regularization parameter.
 
 For advanced usage, first define a solver "R2Solver" to preallocate the memory used in the algorithm, and then call `solve!`:
 
@@ -152,7 +152,7 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u
 - `ν::T = eps(T)^(1 / 5)`: multiplicative inverse of the regularization parameter: ν = 1/σ;
 - `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful.
 
-The algorithm stops either when `√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is stationarity measure.
+The algorithm stops either when `√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure.
 
 # Output
 The value returned is a `GenericExecutionStats`, see `SolverCore.jl`.