From 2f79960c54a47628bd0d02e5a6890b0a5cd4266b Mon Sep 17 00:00:00 2001
From: Geoffroy Leconte <gleconte50@gmail.com>
Date: Fri, 15 Mar 2024 13:39:19 -0400
Subject: [PATCH 1/9] plots

---
 benchmarks/tables/bpdn-constr-table.jl |  37 +++++---
 benchmarks/tables/bpdn-table.jl        |  15 ++++
 benchmarks/tables/fh-table.jl          |  63 +++++++++-----
 benchmarks/tables/nnmf-table.jl        |  40 +++++++--
 benchmarks/tables/regulopt-plots.jl    | 115 +++++++++++++++++++++++++
 benchmarks/tables/regulopt-tables.jl   |   1 +
 benchmarks/tables/svm-table.jl         |  35 ++++++--
 7 files changed, 255 insertions(+), 51 deletions(-)
 create mode 100644 benchmarks/tables/regulopt-plots.jl

diff --git a/benchmarks/tables/bpdn-constr-table.jl b/benchmarks/tables/bpdn-constr-table.jl
index 3c4207a2..3f9cb181 100644
--- a/benchmarks/tables/bpdn-constr-table.jl
+++ b/benchmarks/tables/bpdn-constr-table.jl
@@ -125,17 +125,32 @@ subsolver_options = [
   options4,
   options4_nrTR,
 ] # n'importe lequel si subsolver = :None
-stats = benchmark_table(
+# stats = benchmark_table(
+#   f,
+#   1:(f.meta.nvar),
+#   sol,
+#   h,
+#   λ,
+#   solvers,
+#   subsolvers,
+#   solver_options,
+#   subsolver_options,
+#   "BPDN-cstr",
+#   random_seed,
+#   tex = false,
+# );
+
+subset = [1, 2, 3, 4, 5, 6, 7]
+
+p = benchmark_plot(
   f,
   1:(f.meta.nvar),
-  sol,
   h,
-  λ,
-  solvers,
-  subsolvers,
-  solver_options,
-  subsolver_options,
-  "BPDN-cstr",
-  random_seed,
-  tex = false,
-);
+  solvers[subset],
+  subsolvers[subset],
+  solver_options[subset],
+  subsolver_options[subset],
+  random_seed;
+  xmode = "linear",
+  ymode = "log", 
+)
\ No newline at end of file
diff --git a/benchmarks/tables/bpdn-table.jl b/benchmarks/tables/bpdn-table.jl
index dfa5768c..ec2877ae 100644
--- a/benchmarks/tables/bpdn-table.jl
+++ b/benchmarks/tables/bpdn-table.jl
@@ -141,3 +141,18 @@ stats = benchmark_table(
   random_seed,
   tex = false,
 );
+
+subset = [1,2,3,4,5,6, 7]
+
+p = benchmark_plot(
+  f,
+  1:(f.meta.nvar),
+  h,
+  solvers[subset],
+  subsolvers[subset],
+  solver_options[subset],
+  subsolver_options[subset],
+  random_seed;
+  xmode = "linear",
+  ymode = "log", 
+)
\ No newline at end of file
diff --git a/benchmarks/tables/fh-table.jl b/benchmarks/tables/fh-table.jl
index 4584e2e0..635614c7 100644
--- a/benchmarks/tables/fh-table.jl
+++ b/benchmarks/tables/fh-table.jl
@@ -1,4 +1,5 @@
 include("regulopt-tables.jl")
+include("regulopt-plots.jl")
 using ADNLPModels, DifferentialEquations
 
 display_sol = true
@@ -133,33 +134,49 @@ subsolver_options = [
 ] # n'importe lequel si subsolver = :None
 subset = 2:length(solvers) # issues with R2 alone
 
-names, stats = benchmark_table(
+# names, stats = benchmark_table(
+#   f,
+#   1:(f.meta.nvar),
+#   x0,
+#   h,
+#   λ,
+#   solvers[subset],
+#   subsolvers[subset],
+#   solver_options[subset],
+#   subsolver_options[subset],
+#   "FH with ν = $ν, λ = $λ",
+#   random_seed,
+#   tex = false,
+# );
+
+# if display_sol
+#   data = zeros(length(subset) + 1, 5)
+#   data[1, :] .= x0
+#   for i = 1:length(subset)
+#     data[i + 1, :] .= stats[i].solution
+#   end
+#   pretty_table(
+#     data;
+#     header = [L"$x_1$", L"$x_2$", L"$x_3$", L"$x_4$", L"$x_5$"],
+#     row_names = vcat(["True"], names),
+#     title = "Solution FH",
+#     formatters = ft_printf("%1.2f"),
+#     # backend = Val(:latex),
+#   )
+# end
+
+subset = [8, 9, 10, 11, 12, 13, 14]
+
+p = benchmark_plot(
   f,
   1:(f.meta.nvar),
-  x0,
   h,
-  λ,
   solvers[subset],
   subsolvers[subset],
   solver_options[subset],
   subsolver_options[subset],
-  "FH with ν = $ν, λ = $λ",
-  random_seed,
-  tex = false,
-);
-
-if display_sol
-  data = zeros(length(subset) + 1, 5)
-  data[1, :] .= x0
-  for i = 1:length(subset)
-    data[i + 1, :] .= stats[i].solution
-  end
-  pretty_table(
-    data;
-    header = [L"$x_1$", L"$x_2$", L"$x_3$", L"$x_4$", L"$x_5$"],
-    row_names = vcat(["True"], names),
-    title = "Solution FH",
-    formatters = ft_printf("%1.2f"),
-    # backend = Val(:latex),
-  )
-end
+  random_seed;
+  measured = :grad,
+  xmode = "linear",
+  ymode = "log", 
+)
\ No newline at end of file
diff --git a/benchmarks/tables/nnmf-table.jl b/benchmarks/tables/nnmf-table.jl
index 20cabfc6..4a434cba 100644
--- a/benchmarks/tables/nnmf-table.jl
+++ b/benchmarks/tables/nnmf-table.jl
@@ -1,4 +1,5 @@
 include("regulopt-tables.jl")
+include("regulopt-plots.jl")
 
 random_seed = 1234
 Random.seed!(random_seed)
@@ -122,19 +123,42 @@ subsolver_options = [
   options4,
   options4_nrTR,
 ] # n'importe lequel si subsolver = :None
-subset = 1:length(solvers)
 
-benchmark_table(
+# benchmark_table(
+#   f,
+#   selected,
+#   [],
+#   h,
+#   λ,
+#   solvers[subset],
+#   subsolvers[subset],
+#   solver_options[subset],
+#   subsolver_options[subset],
+#   "NNMF with m = $m, n = $n, k = $k, ν = $ν, λ = $λ",
+#   random_seed,
+#   tex = false,
+# );
+
+# subset = 1:length(solvers)
+subset = [8, 9, 10, 11, 12, 13, 14]
+
+p = benchmark_plot(
   f,
   selected,
-  [],
   h,
-  λ,
   solvers[subset],
   subsolvers[subset],
   solver_options[subset],
   subsolver_options[subset],
-  "NNMF with m = $m, n = $n, k = $k, ν = $ν, λ = $λ",
-  random_seed,
-  tex = false,
-);
+  random_seed;
+  measured = :grad,
+  xmode = "linear",
+  ymode = "log", 
+)
+
+# path = raw"C:\Users\Geoffroy Leconte\.julia\dev\RegularizedOptimization\benchmarks"
+# pgfsave(
+#   joinpath(path, "test.tikz"),
+#   p;
+#   include_preamble = true,
+# )
\ No newline at end of file
diff --git a/benchmarks/tables/regulopt-plots.jl b/benchmarks/tables/regulopt-plots.jl
new file mode 100644
index 00000000..49919929
--- /dev/null
+++ b/benchmarks/tables/regulopt-plots.jl
@@ -0,0 +1,115 @@
+using PGFPlotsX
+using Colors
+using LaTeXStrings
+using PrettyTables, LaTeXStrings
+using Random
+using LinearAlgebra
+using ProximalOperators
+using NLPModels,
+  NLPModelsModifiers,
+  RegularizedProblems,
+  RegularizedOptimization,
+  ShiftedProximalOperators,
+  SolverBenchmark
+using Printf
+
+# utils for extracting stats / display table
+modelname(nlp::LSR1Model) = "LSR1"
+modelname(nlp::LBFGSModel) = "LBFGS"
+modelname(nlp::SpectralGradientModel) = "SpectralGradient"
+modelname(nlp::DiagonalQNModel) = "DiagonalQN"
+subsolvername(subsolver::Symbol) = subsolver == :None ? "" : string("-", subsolver)
+function options_str(
+  options::ROSolverOptions,
+  solver::Symbol,
+  subsolver_options::ROSolverOptions,
+  subsolver::Symbol,
+)
+  if solver == :TRDH
+    out_str = !options.spectral ? (options.psb ? "-PSB" : "-Andrei") : "-Spec"
+    out_str = (options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+  elseif solver == :TR && subsolver == :TRDH
+    out_str = !subsolver_options.spectral ? (subsolver_options.psb ? "-PSB" : "-Andrei") : "-Spec"
+    out_str = (subsolver_options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+  else
+    out_str = ""
+  end
+  return out_str
+end
+
+function benchmark_plot(
+  f::AbstractNLPModel,
+  selected,
+  h,
+  solvers,
+  subsolvers,
+  solver_options,
+  subsolver_options,
+  random_seed::Int;
+  measured::Symbol = :obj, # set to :grad to eval grad
+  xmode::String = "log",
+  ymode::String = "log",
+)
+  solver_names = [
+    "$(solver)$(subsolvername(subsolver))$(options_str(opt, solver, subsolver_opt, subsolver))"
+    for (solver, opt, subsolver, subsolver_opt) in
+    zip(solvers, solver_options, subsolvers, subsolver_options)
+  ]
+  n_solvers = length(solver_names)
+  objdecs = Vector{Float64}[] 
+  coords = Coordinates{2}[]
+  obj_min = Float64(Inf)
+
+  for (solver, subsolver, opt, sub_opt) in
+      zip(solvers, subsolvers, solver_options, subsolver_options)
+    @info " using $solver with subsolver = $subsolver"
+    args = solver == :R2 ? () : (NormLinf(1.0),)
+    Random.seed!(random_seed)
+    if subsolver == :None
+      solver_out = eval(solver)(f, h, args..., opt, x0 = f.meta.x0, selected = selected)
+    else
+      solver_out = eval(solver)(
+        f,
+        h,
+        args...,
+        opt,
+        x0 = f.meta.x0,
+        subsolver = eval(subsolver),
+        subsolver_options = sub_opt,
+        selected = selected,
+      )
+    end
+    objdec = solver_out.solver_specific[:Fhist] + solver_out.solver_specific[:Hhist]
+    measured == :grad && unique!(objdec)
+    obj_min = min(minimum(objdec), obj_min)
+    push!(objdecs, objdec)
+    reset!(f)
+  end
+  for i in 1:length(objdecs)
+    objdec = objdecs[i]
+    push!(
+      coords,
+      Coordinates([(k-1, objdec[k] - obj_min) for k in 1:length(objdec)]),
+    )
+  end
+
+  colors = distinguishable_colors(n_solvers)
+  l_plots = [@pgf Plot({color = colors[i]}, coords[i]) for i in 1:n_solvers]
+  
+  @pgf Axis(
+    {
+      xlabel = "iterations",
+      ylabel = L"$(f + h)(x_k)$",
+      ymode = ymode,
+      xmode = xmode,
+      no_markers,
+      legend_style = {
+        nodes={scale=0.8},
+        font = "\\tiny",
+      },
+      # legend_pos="south west",
+    },
+    Tuple(l_plots)...,
+    Legend(solver_names),
+  )
+end
diff --git a/benchmarks/tables/regulopt-tables.jl b/benchmarks/tables/regulopt-tables.jl
index ffaddff1..f5b7d9d2 100644
--- a/benchmarks/tables/regulopt-tables.jl
+++ b/benchmarks/tables/regulopt-tables.jl
@@ -74,6 +74,7 @@ function benchmark_table(
   n∇f_evals = []
   nprox_evals = []
   solver_stats = []
+  reset!(f)
 
   for (solver, subsolver, opt, sub_opt) in
       zip(solvers, subsolvers, solver_options, subsolver_options)
diff --git a/benchmarks/tables/svm-table.jl b/benchmarks/tables/svm-table.jl
index ed937972..0f22a106 100644
--- a/benchmarks/tables/svm-table.jl
+++ b/benchmarks/tables/svm-table.jl
@@ -1,4 +1,5 @@
 include("regulopt-tables.jl")
+include("regulopt-plots.jl")
 using MLDatasets
 
 random_seed = 1234
@@ -127,19 +128,35 @@ subsolver_options = [
 ] # n'importe lequel si subsolver = :None
 subset = 1:length(solvers)
 
-benchmark_table(
+# benchmark_table(
+#   f,
+#   1:(f.meta.nvar),
+#   (sol_train, sol_test),
+#   h,
+#   λ,
+#   solvers[subset],
+#   subsolvers[subset],
+#   solver_options[subset],
+#   subsolver_options[subset],
+#   "SVM with ν = $ν, λ = $λ",
+#   random_seed,
+#   nls_train = nls_train,
+#   nls_test = nls_test,
+#   tex = false,
+# );
+
+subset = [8, 9, 10, 11, 12, 13, 14]
+
+p = benchmark_plot(
   f,
   1:(f.meta.nvar),
-  (sol_train, sol_test),
   h,
-  λ,
   solvers[subset],
   subsolvers[subset],
   solver_options[subset],
   subsolver_options[subset],
-  "SVM with ν = $ν, λ = $λ",
-  random_seed,
-  nls_train = nls_train,
-  nls_test = nls_test,
-  tex = false,
-);
+  random_seed;
+  measured = :grad,
+  xmode = "linear",
+  ymode = "log", 
+)
\ No newline at end of file

From dc0972ae6de44ba1cb8f4bfd6ef25cdeefb9291b Mon Sep 17 00:00:00 2001
From: Geoffroy Leconte <gleconte50@gmail.com>
Date: Wed, 27 Mar 2024 12:23:51 -0400
Subject: [PATCH 2/9] update benchmarks

---
 benchmarks/tables/benchmark-profile.jl | 149 +++++++++++++++++++++++++
 benchmarks/tables/bpdn-constr-table.jl |  30 ++---
 benchmarks/tables/bpdn-table.jl        |   4 +-
 benchmarks/tables/fh-table.jl          |  60 +++++-----
 benchmarks/tables/nnmf-table.jl        |  49 ++++----
 benchmarks/tables/regulopt-tables.jl   |   1 +
 benchmarks/tables/svm-table.jl         |  32 +++---
 examples/plot-utils-nnmf.jl            |  47 ++++----
 src/TRDH_alg.jl                        |   8 +-
 9 files changed, 269 insertions(+), 111 deletions(-)
 create mode 100644 benchmarks/tables/benchmark-profile.jl

diff --git a/benchmarks/tables/benchmark-profile.jl b/benchmarks/tables/benchmark-profile.jl
new file mode 100644
index 00000000..f8813660
--- /dev/null
+++ b/benchmarks/tables/benchmark-profile.jl
@@ -0,0 +1,149 @@
+using LaTeXStrings
+using PrettyTables, LaTeXStrings
+using Random
+using LinearAlgebra
+using ProximalOperators
+using Plots
+using BenchmarkProfiles
+using NLPModels,
+  NLPModelsModifiers,
+  RegularizedProblems,
+  RegularizedOptimization,
+  ShiftedProximalOperators,
+  SolverBenchmark
+using Printf
+
+# utils for extracting stats / display table
+modelname(nlp::LSR1Model) = "LSR1"
+modelname(nlp::LBFGSModel) = "LBFGS"
+modelname(nlp::SpectralGradientModel) = "SpectralGradient"
+modelname(nlp::DiagonalQNModel) = "DiagonalQN"
+subsolvername(subsolver::Symbol) = subsolver == :None ? "" : string("-", subsolver)
+function options_str(
+  options::ROSolverOptions,
+  solver::Symbol,
+  subsolver_options::ROSolverOptions,
+  subsolver::Symbol,
+)
+  if solver == :TRDH
+    out_str = !options.spectral ? (options.psb ? "-PSB" : "-Andrei") : "-Spec"
+    out_str = (options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+  elseif solver == :TR && subsolver == :TRDH
+    out_str = !subsolver_options.spectral ? (subsolver_options.psb ? "-PSB" : "-Andrei") : "-Spec"
+    out_str = (subsolver_options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+  else
+    out_str = ""
+  end
+  return out_str
+end
+
+function benchmark_prof(
+  pb::Symbol, #:nnmf, :bpdn
+  solvers,
+  solver_names,
+  nb_prob::Int,
+  random_seed::Int;
+  measured::Symbol = :obj, # set to :grad to eval grad
+)
+
+  if pb == :nnmf
+    m, n, k = 100, 50, 5
+    λ = 1.0e-1
+  elseif pb == :bpdn
+    compound = 1
+  else
+    error("Problem not supported")
+  end
+  n_solvers = length(solvers)
+  data = zeros(nb_prob, n_solvers)
+
+  current_seed = random_seed
+  nb_data_min = Int(10^16) # min number of data for each solver on every pb
+  for i=1:nb_prob
+    current_seed = random_seed + i - 1
+    Random.seed!(current_seed)
+    if pb == :nnmf
+      model, nls_model, A, selected = nnmf_model(m, n, k)
+      h = NormL0(λ)
+    elseif pb == :bpdn
+      model, nls_model, sol = bpdn_model(compound, bounds = false)
+      selected = 1:length(sol)
+      λ = norm(grad(model, zeros(model.meta.nvar)), Inf) / 10
+      reset!(model)
+      h = NormL0(λ)
+    end
+    f = LSR1Model(model)
+    @info "pb $i"
+    for (j, solver, name) in zip(1:n_solvers, solvers, solver_names)
+      solver_out = solver(f, h, selected)
+      @info "pb: $i  solver: $name  status = $(solver_out.status)  obj = $(solver_out.objective)"
+      if solver_out.status == :first_order
+        data[i, j] = neval_grad(f)
+      else
+        data[i, j] = +Inf
+      end
+      reset!(f)
+    end
+  end
+
+  performance_profile(PlotsBackend(), data, solver_names, legend=:bottomright, title = String(measured))
+end
+
+ν = 1.0
+ϵ = 1.0e-5
+ϵi = 1.0e-3
+ϵri = 1.0e-6
+maxIter = 2000
+maxIter_inner = 100
+function TR_R2(f, h, selected; ϵ = ϵ, ϵi = ϵi, ϵri = ϵri, maxIter = maxIter, maxIter_inner = maxIter_inner)
+  opt = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = 0, maxIter = maxIter)
+  sub_opt = ROSolverOptions(ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+  TR(f, h, NormLinf(1.0), opt, x0 = f.meta.x0, subsolver_options = sub_opt, selected = selected)
+end
+
+function TR_TRDH_PSB(f, h, selected; ϵ = ϵ, ϵi = ϵi, ϵri = ϵri, maxIter = maxIter, maxIter_inner = maxIter_inner)
+  opt = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = 0, maxIter = maxIter)
+  sub_opt = ROSolverOptions(ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, psb = true, reduce_TR = false)
+  TR(f, h, NormLinf(1.0), opt, x0 = f.meta.x0, subsolver_options = sub_opt, selected = selected, subsolver = TRDH)
+end
+
+function TR_TRDH_Andrei(f, h, selected; ϵ = ϵ, ϵi = ϵi, ϵri = ϵri, maxIter = maxIter, maxIter_inner = maxIter_inner)
+  opt = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = 0, maxIter = maxIter)
+  sub_opt = ROSolverOptions(ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, spectral = false, psb = false, reduce_TR = false)
+  TR(f, h, NormLinf(1.0), opt, x0 = f.meta.x0, subsolver_options = sub_opt, selected = selected, subsolver = TRDH)
+end
+
+function TR_TRDH_Spec(f, h, selected; ϵ = ϵ, ϵi = ϵi, ϵri = ϵri, maxIter = maxIter, maxIter_inner = maxIter_inner)
+  opt = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = 0, maxIter = maxIter)
+  sub_opt = ROSolverOptions(ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, spectral = true, reduce_TR = false)
+  TR(f, h, NormLinf(1.0), opt, x0 = f.meta.x0, subsolver_options = sub_opt, selected = selected, subsolver = TRDH)
+end
+
+function TRDH_Spec(f, h, selected; ϵ = ϵ, ϵi = ϵi, ϵri = ϵri, maxIter = maxIter, maxIter_inner = maxIter_inner)
+  opt = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, maxIter = maxIter, spectral = true, reduce_TR = false)
+  TRDH(f, h, NormLinf(1.0), opt, x0 = f.meta.x0, selected = selected)
+end
+
+function R2_None(f, h, selected; ϵ = ϵ, ϵi = ϵi, ϵri = ϵri, maxIter = maxIter, maxIter_inner = maxIter_inner)
+  opt = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, maxIter = maxIter, spectral = true, reduce_TR = false)
+  R2(f, h, opt; x0 = f.meta.x0, selected = selected)
+end
+
+
+benchmark_prof(
+  :nnmf,
+  [TRDH_Spec, TR_R2, TR_TRDH_PSB, TR_TRDH_Andrei, TR_TRDH_Spec],
+  ["TRDH_Spec", "TR-R2", "TR-TRDH-PSB", "TR_TRDH_Andrei", "TR_TRDH_Spec"],
+  50,
+  1234;
+  measured = :grad, # set to :grad to eval grad
+)
+
+benchmark_prof(
+  :bpdn,
+  [R2_None, TRDH_Spec, TR_R2, TR_TRDH_PSB],
+  ["R2", "TRDH-Spec", "TR-R2", "TR-TRDH-PSB"],
+  50,
+  1234;
+  measured = :grad, # set to :grad to eval grad
+)
\ No newline at end of file
diff --git a/benchmarks/tables/bpdn-constr-table.jl b/benchmarks/tables/bpdn-constr-table.jl
index 3f9cb181..dd1a0aa8 100644
--- a/benchmarks/tables/bpdn-constr-table.jl
+++ b/benchmarks/tables/bpdn-constr-table.jl
@@ -13,7 +13,7 @@ h = NormL0(λ)
 verbose = 0 # 10
 ν = 1.0
 ϵ = 1.0e-5
-ϵi = 1.0e-5
+ϵi = 1.0e-3
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
@@ -125,20 +125,20 @@ subsolver_options = [
   options4,
   options4_nrTR,
 ] # n'importe lequel si subsolver = :None
-# stats = benchmark_table(
-#   f,
-#   1:(f.meta.nvar),
-#   sol,
-#   h,
-#   λ,
-#   solvers,
-#   subsolvers,
-#   solver_options,
-#   subsolver_options,
-#   "BPDN-cstr",
-#   random_seed,
-#   tex = false,
-# );
+stats = benchmark_table(
+  f,
+  1:(f.meta.nvar),
+  sol,
+  h,
+  λ,
+  solvers,
+  subsolvers,
+  solver_options,
+  subsolver_options,
+  "BPDN-cstr",
+  random_seed,
+  tex = true,
+);
 
 subset = [1, 2, 3, 4, 5, 6, 7]
 
diff --git a/benchmarks/tables/bpdn-table.jl b/benchmarks/tables/bpdn-table.jl
index ec2877ae..3a22e874 100644
--- a/benchmarks/tables/bpdn-table.jl
+++ b/benchmarks/tables/bpdn-table.jl
@@ -14,7 +14,7 @@ h = NormL0(λ)
 verbose = 0 # 10
 ν = 1.0
 ϵ = 1.0e-5
-ϵi = 1.0e-5
+ϵi = 1.0e-3
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
@@ -139,7 +139,7 @@ stats = benchmark_table(
   subsolver_options,
   "BPDN",
   random_seed,
-  tex = false,
+  tex = true,
 );
 
 subset = [1,2,3,4,5,6, 7]
diff --git a/benchmarks/tables/fh-table.jl b/benchmarks/tables/fh-table.jl
index 635614c7..96ba529a 100644
--- a/benchmarks/tables/fh-table.jl
+++ b/benchmarks/tables/fh-table.jl
@@ -7,7 +7,7 @@ display_sol = true
 random_seed = 12345
 Random.seed!(random_seed)
 
-cstr = false
+cstr = true
 ctr_val = cstr ? 0.5 : -Inf
 lvar = [-Inf, ctr_val, -Inf, -Inf, -Inf]
 uvar = fill(Inf, 5)
@@ -134,36 +134,36 @@ subsolver_options = [
 ] # n'importe lequel si subsolver = :None
 subset = 2:length(solvers) # issues with R2 alone
 
-# names, stats = benchmark_table(
-#   f,
-#   1:(f.meta.nvar),
-#   x0,
-#   h,
-#   λ,
-#   solvers[subset],
-#   subsolvers[subset],
-#   solver_options[subset],
-#   subsolver_options[subset],
-#   "FH with ν = $ν, λ = $λ",
-#   random_seed,
-#   tex = false,
-# );
+names, stats = benchmark_table(
+  f,
+  1:(f.meta.nvar),
+  x0,
+  h,
+  λ,
+  solvers[subset],
+  subsolvers[subset],
+  solver_options[subset],
+  subsolver_options[subset],
+  "FH with ν = $ν, λ = $λ",
+  random_seed,
+  tex = true,
+);
 
-# if display_sol
-#   data = zeros(length(subset) + 1, 5)
-#   data[1, :] .= x0
-#   for i = 1:length(subset)
-#     data[i + 1, :] .= stats[i].solution
-#   end
-#   pretty_table(
-#     data;
-#     header = [L"$x_1$", L"$x_2$", L"$x_3$", L"$x_4$", L"$x_5$"],
-#     row_names = vcat(["True"], names),
-#     title = "Solution FH",
-#     formatters = ft_printf("%1.2f"),
-#     # backend = Val(:latex),
-#   )
-# end
+if display_sol
+  data = zeros(length(subset) + 1, 5)
+  data[1, :] .= x0
+  for i = 1:length(subset)
+    data[i + 1, :] .= stats[i].solution
+  end
+  pretty_table(
+    data;
+    header = [L"$x_1$", L"$x_2$", L"$x_3$", L"$x_4$", L"$x_5$"],
+    row_names = vcat(["True"], names),
+    title = "Solution FH",
+    formatters = ft_printf("%1.2f"),
+    # backend = Val(:latex),
+  )
+end
 
 subset = [8, 9, 10, 11, 12, 13, 14]
 
diff --git a/benchmarks/tables/nnmf-table.jl b/benchmarks/tables/nnmf-table.jl
index 4a434cba..b61312b8 100644
--- a/benchmarks/tables/nnmf-table.jl
+++ b/benchmarks/tables/nnmf-table.jl
@@ -6,7 +6,7 @@ Random.seed!(random_seed)
 m, n, k = 100, 50, 5
 model, nls_model, A, selected = nnmf_model(m, n, k)
 f = LSR1Model(model)
-λ = 1.0e-1 # norm(grad(model, rand(model.meta.nvar)), Inf) / 100
+λ = 1.0e-3 # norm(grad(model, rand(model.meta.nvar)), Inf) / 100
 h = NormL0(λ)
 ν = 1.0
 ϵ = 1.0e-5
@@ -124,37 +124,38 @@ subsolver_options = [
   options4_nrTR,
 ] # n'importe lequel si subsolver = :None
 
-# benchmark_table(
-#   f,
-#   selected,
-#   [],
-#   h,
-#   λ,
-#   solvers[subset],
-#   subsolvers[subset],
-#   solver_options[subset],
-#   subsolver_options[subset],
-#   "NNMF with m = $m, n = $n, k = $k, ν = $ν, λ = $λ",
-#   random_seed,
-#   tex = false,
-# );
+subset = 1:length(solvers)
 
-# subset = 1:length(solvers)
-subset = [8, 9, 10, 11, 12, 13, 14]
-
-p = benchmark_plot(
+benchmark_table(
   f,
   selected,
+  [],
   h,
+  λ,
   solvers[subset],
   subsolvers[subset],
   solver_options[subset],
   subsolver_options[subset],
-  random_seed;
-  measured = :grad,
-  xmode = "linear",
-  ymode = "log", 
-)
+  "NNMF with m = $m, n = $n, k = $k, ν = $ν, λ = $λ",
+  random_seed,
+  tex = true,
+);
+
+subset = [8, 9, 10, 11, 12, 13, 14]
+
+# p = benchmark_plot(
+#   f,
+#   selected,
+#   h,
+#   solvers[subset],
+#   subsolvers[subset],
+#   solver_options[subset],
+#   subsolver_options[subset],
+#   random_seed;
+#   measured = :grad,
+#   xmode = "linear",
+#   ymode = "log", 
+# )
 
 # path = raw"C:\Users\Geoffroy Leconte\.julia\dev\RegularizedOptimization\benchmarks"
 # pgfsave(
diff --git a/benchmarks/tables/regulopt-tables.jl b/benchmarks/tables/regulopt-tables.jl
index f5b7d9d2..4c5e0f47 100644
--- a/benchmarks/tables/regulopt-tables.jl
+++ b/benchmarks/tables/regulopt-tables.jl
@@ -95,6 +95,7 @@ function benchmark_table(
         selected = selected,
       )
     end
+    println(length(findall(isequal(0), solver_out.solution[selected])))
     push!(nf_evals, obj_evals(f))
     push!(n∇f_evals, grad_evals(f))
     push!(nprox_evals, nb_prox_evals(solver_out, solver))
diff --git a/benchmarks/tables/svm-table.jl b/benchmarks/tables/svm-table.jl
index 0f22a106..d864d9f3 100644
--- a/benchmarks/tables/svm-table.jl
+++ b/benchmarks/tables/svm-table.jl
@@ -128,22 +128,22 @@ subsolver_options = [
 ] # n'importe lequel si subsolver = :None
 subset = 1:length(solvers)
 
-# benchmark_table(
-#   f,
-#   1:(f.meta.nvar),
-#   (sol_train, sol_test),
-#   h,
-#   λ,
-#   solvers[subset],
-#   subsolvers[subset],
-#   solver_options[subset],
-#   subsolver_options[subset],
-#   "SVM with ν = $ν, λ = $λ",
-#   random_seed,
-#   nls_train = nls_train,
-#   nls_test = nls_test,
-#   tex = false,
-# );
+benchmark_table(
+  f,
+  1:(f.meta.nvar),
+  (sol_train, sol_test),
+  h,
+  λ,
+  solvers[subset],
+  subsolvers[subset],
+  solver_options[subset],
+  subsolver_options[subset],
+  "SVM with ν = $ν, λ = $λ",
+  random_seed,
+  nls_train = nls_train,
+  nls_test = nls_test,
+  tex = true,
+);
 
 subset = [8, 9, 10, 11, 12, 13, 14]
 
diff --git a/examples/plot-utils-nnmf.jl b/examples/plot-utils-nnmf.jl
index ed36381e..ada66794 100644
--- a/examples/plot-utils-nnmf.jl
+++ b/examples/plot-utils-nnmf.jl
@@ -9,44 +9,51 @@ function plot_nnmf(outstruct, Avec, m, n, k, name = "tr-qr")
   H = reshape(x[(m * k + 1):end], k, n)
   WH = W * H
 
-  a = GroupPlot(2, 2, groupStyle = "horizontal sep = 2.5cm")
+  a = PGFPlots.GroupPlot(2, 2, groupStyle = "horizontal sep = 2.5cm")
   push!(
     a,
-    Axis(
+    PGFPlots.Axis(
       Plots.Image(A, (1, m), (1, n), colormap = ColorMaps.Named("Jet")),
       xlabel = "A matrix (reference)",
     ),
   )
   push!(
     a,
-    Axis(Plots.Image(WH, (1, m), (1, n), colormap = ColorMaps.Named("Jet")), xlabel = "WH matrix"),
+    PGFPlots.Axis(Plots.Image(WH, (1, m), (1, n), colormap = ColorMaps.Named("Jet")), xlabel = "WH matrix"),
   )
   push!(
     a,
-    Axis(Plots.Image(H, (1, k), (1, n), colormap = ColorMaps.Named("Jet")), xlabel = "H matrix"),
+    PGFPlots.Axis(Plots.Image(H, (1, k), (1, n), colormap = ColorMaps.Named("Jet")), xlabel = "H matrix"),
   )
   push!(
     a,
-    Axis(
+    PGFPlots.Axis(
       Plots.Image(abs.(A - WH), (1, m), (1, n), colormap = ColorMaps.Named("Jet")),
       xlabel = "|A-WH| matrix",
     ),
   )
-  save("nnmf-$(name).pdf", a)
+  a = Plots.MatrixPlot(A, (1, m), (1, n), colormap = ColorMaps.Named("Jet"))
+  save("nnmf-A-$(name).tikz", a)
+  b = Plots.MatrixPlot(WH, (1, m), (1, n), colormap = ColorMaps.Named("Jet"))
+  save("nnmf-WH-$(name).tikz", b)
+  c = Plots.Image(H, (1, k), (1, n), colormap = ColorMaps.Named("Jet"))
+  save("nnmf-H-$(name).tikz", c)
+  d = Plots.MatrixPlot(abs.(A - WH), (1, m), (1, n), colormap = ColorMaps.Named("Jet"))
+  save("nnmf-A-WH-$(name).tikz", d)
 
-  b = Axis(
-    Plots.Linear(1:length(Comp_pg), Comp_pg, mark = "none"),
-    xlabel = "outer iterations",
-    ylabel = "inner iterations",
-    ymode = "log",
-  )
-  save("nnmf-inner-outer-$(name).pdf", b)
+  # b = Axis(
+  #   Plots.Linear(1:length(Comp_pg), Comp_pg, mark = "none"),
+  #   xlabel = "outer iterations",
+  #   ylabel = "inner iterations",
+  #   ymode = "log",
+  # )
+  # save("nnmf-inner-outer-$(name).pdf", b)
 
-  c = Axis(
-    Plots.Linear(1:length(objdec), objdec, mark = "none"),
-    xlabel = "\$ k^{th}\$  \$ \\nabla f \$ Call",
-    ylabel = "Objective Value",
-    ymode = "log",
-  )
-  save("nnmf-objdec-$(name).pdf", c)
+  # c = Axis(
+  #   Plots.Linear(1:length(objdec), objdec, mark = "none"),
+  #   xlabel = "\$ k^{th}\$  \$ \\nabla f \$ Call",
+  #   ylabel = "Objective Value",
+  #   ymode = "log",
+  # )
+  # save("nnmf-objdec-$(name).pdf", c)
 end
diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl
index 8fef5479..d72453e0 100644
--- a/src/TRDH_alg.jl
+++ b/src/TRDH_alg.jl
@@ -205,8 +205,8 @@ function TRDH(
   ∇fk⁻ = copy(∇fk)
   Dk = spectral ? SpectralGradient(hess_init_val, length(xk)) :
     ((Bk isa UniformScaling) ? DiagonalQN(fill!(similar(xk), hess_init_val), psb) : DiagonalQN(diag(Bk), psb))
-  DkNorm = norm(Dk.d, Inf)
-  νInv = (DkNorm + one(R) / (α * Δk))
+  DkNorm⁺ = norm(max.(Dk.d, 0), Inf)
+  νInv = (DkNorm⁺ + one(R) / (α * Δk))
   ν = one(R) / νInv
   mν∇fk = -ν .* ∇fk
   sqrt_ξ_νInv = one(R)
@@ -315,7 +315,7 @@ function TRDH(
       shift!(ψ, xk)
       ∇f!(∇fk, xk)
       push!(Dk, s, ∇fk - ∇fk⁻) # update QN operator
-      DkNorm = norm(Dk.d, Inf)
+      DkNorm⁺ = norm(max.(Dk.d, 0), Inf)
       ∇fk⁻ .= ∇fk
     end
 
@@ -325,7 +325,7 @@ function TRDH(
       has_bnds ? set_bounds!(ψ, l_bound_k, u_bound_k) : set_radius!(ψ, Δk)
     end
 
-    νInv = reduce_TR ? (DkNorm + one(R) / (α * Δk)) : (DkNorm + one(R) / α)
+    νInv = reduce_TR ? (DkNorm⁺ + one(R) / (α * Δk)) : (DkNorm⁺ + one(R) / α)
     ν = one(R) / νInv
     mν∇fk .= -ν .* ∇fk
 

From c96b216748524ebf95478810d81b64a544cf505d Mon Sep 17 00:00:00 2001
From: Geoffroy Leconte <gleconte50@gmail.com>
Date: Tue, 6 Jun 2023 16:43:57 -0400
Subject: [PATCH 3/9] WIP nonmonotone

---
 benchmarks/tables/bpdn-constr-table.jl | 20 ++++++++++++++------
 benchmarks/tables/bpdn-table.jl        | 18 +++++++++++++-----
 benchmarks/tables/fh-table.jl          | 20 ++++++++++++++------
 benchmarks/tables/nnmf-table.jl        | 20 ++++++++++++++------
 benchmarks/tables/regulopt-tables.jl   |  2 ++
 benchmarks/tables/svm-table.jl         | 20 ++++++++++++++------
 src/TRDH_alg.jl                        | 12 +++++++++++-
 src/input_struct.jl                    |  3 +++
 8 files changed, 85 insertions(+), 30 deletions(-)

diff --git a/benchmarks/tables/bpdn-constr-table.jl b/benchmarks/tables/bpdn-constr-table.jl
index dd1a0aa8..d3787584 100644
--- a/benchmarks/tables/bpdn-constr-table.jl
+++ b/benchmarks/tables/bpdn-constr-table.jl
@@ -17,8 +17,9 @@ verbose = 0 # 10
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
+Mmonotone = 0
 options =
-  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -27,8 +28,9 @@ options_nrTR = ROSolverOptions(
   maxIter = maxIter,
   spectral = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options2_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
@@ -36,9 +38,10 @@ options2_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options3 =
-  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options3_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
@@ -46,10 +49,11 @@ options3_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options4_nrTR =
-  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false)
+  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false, Mmonotone = Mmonotone)
 options5 = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -58,6 +62,7 @@ options5 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = true,
+  Mmonotone = Mmonotone,
 )
 options5_nrTR = ROSolverOptions(
   ν = ν,
@@ -68,6 +73,7 @@ options5_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options6 = ROSolverOptions(
   ν = ν,
@@ -77,6 +83,7 @@ options6 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = false,
+  Mmonotone = Mmonotone,
 )
 options6_nrTR = ROSolverOptions(
   ν = ν,
@@ -87,6 +94,7 @@ options6_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 
 solvers = [:R2, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR, :TR, :TR]
@@ -135,7 +143,7 @@ stats = benchmark_table(
   subsolvers,
   solver_options,
   subsolver_options,
-  "BPDN-cstr",
+  "BPDN-cstr, M = $Mmonotone",
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/bpdn-table.jl b/benchmarks/tables/bpdn-table.jl
index 3a22e874..4b07e844 100644
--- a/benchmarks/tables/bpdn-table.jl
+++ b/benchmarks/tables/bpdn-table.jl
@@ -18,8 +18,9 @@ verbose = 0 # 10
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
+Mmonotone = 0
 options =
-  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -28,8 +29,9 @@ options_nrTR = ROSolverOptions(
   maxIter = maxIter,
   spectral = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options2_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
@@ -37,9 +39,10 @@ options2_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone
 )
 options3 =
-  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options3_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
@@ -47,10 +50,11 @@ options3_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options4_nrTR =
-  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false)
+  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false, Mmonotone = Mmonotone)
 options5 = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -59,6 +63,7 @@ options5 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = true,
+  Mmonotone = Mmonotone,
 )
 options5_nrTR = ROSolverOptions(
   ν = ν,
@@ -69,6 +74,7 @@ options5_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options6 = ROSolverOptions(
   ν = ν,
@@ -78,6 +84,7 @@ options6 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = false,
+  Mmonotone = Mmonotone,
 )
 options6_nrTR = ROSolverOptions(
   ν = ν,
@@ -88,6 +95,7 @@ options6_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 
 solvers = [:R2, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR, :TR, :TR]
diff --git a/benchmarks/tables/fh-table.jl b/benchmarks/tables/fh-table.jl
index 96ba529a..ebdfd04c 100644
--- a/benchmarks/tables/fh-table.jl
+++ b/benchmarks/tables/fh-table.jl
@@ -24,8 +24,9 @@ maxIter_inner = 200 # max iter for subsolver
 ϵ = 1.0e-4
 ϵi = 1.0e-3
 ϵri = 1.0e-6
+Mmonotone = 5
 options =
-  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -34,8 +35,9 @@ options_nrTR = ROSolverOptions(
   maxIter = maxIter,
   spectral = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options2_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
@@ -43,9 +45,10 @@ options2_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options3 =
-  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options3_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
@@ -53,10 +56,11 @@ options3_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options4_nrTR =
-  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false)
+  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false, Mmonotone = Mmonotone)
 options5 = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -65,6 +69,7 @@ options5 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = true,
+  Mmonotone = Mmonotone,
 )
 options5_nrTR = ROSolverOptions(
   ν = ν,
@@ -75,6 +80,7 @@ options5_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options6 = ROSolverOptions(
   ν = ν,
@@ -84,6 +90,7 @@ options6 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = false,
+  Mmonotone = Mmonotone,
 )
 options6_nrTR = ROSolverOptions(
   ν = ν,
@@ -94,6 +101,7 @@ options6_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 
 solvers = [:R2, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR, :TR, :TR]
@@ -144,7 +152,7 @@ names, stats = benchmark_table(
   subsolvers[subset],
   solver_options[subset],
   subsolver_options[subset],
-  "FH with ν = $ν, λ = $λ",
+  "FH with ν = $ν, λ = $λ, M = $Mmonotone",
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/nnmf-table.jl b/benchmarks/tables/nnmf-table.jl
index b61312b8..d69ba6c4 100644
--- a/benchmarks/tables/nnmf-table.jl
+++ b/benchmarks/tables/nnmf-table.jl
@@ -14,9 +14,10 @@ h = NormL0(λ)
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
+Mmonotone = 0
 verbose = 0 #10
 options =
-  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -25,8 +26,9 @@ options_nrTR = ROSolverOptions(
   maxIter = maxIter,
   spectral = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options2_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
@@ -34,9 +36,10 @@ options2_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options3 =
-  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options3_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
@@ -44,10 +47,11 @@ options3_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options4_nrTR =
-  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false)
+  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false, Mmonotone = Mmonotone)
 options5 = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -56,6 +60,7 @@ options5 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = true,
+  Mmonotone = Mmonotone,
 )
 options5_nrTR = ROSolverOptions(
   ν = ν,
@@ -66,6 +71,7 @@ options5_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options6 = ROSolverOptions(
   ν = ν,
@@ -75,6 +81,7 @@ options6 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = false,
+  Mmonotone = Mmonotone,
 )
 options6_nrTR = ROSolverOptions(
   ν = ν,
@@ -85,6 +92,7 @@ options6_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 
 solvers = [:R2, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR, :TR, :TR]
@@ -136,7 +144,7 @@ benchmark_table(
   subsolvers[subset],
   solver_options[subset],
   subsolver_options[subset],
-  "NNMF with m = $m, n = $n, k = $k, ν = $ν, λ = $λ",
+  "NNMF with m = $m, n = $n, k = $k, ν = $ν, λ = $λ, M = $Mmonotone",
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/regulopt-tables.jl b/benchmarks/tables/regulopt-tables.jl
index 4c5e0f47..7cff31e4 100644
--- a/benchmarks/tables/regulopt-tables.jl
+++ b/benchmarks/tables/regulopt-tables.jl
@@ -9,6 +9,7 @@ using NLPModels,
   ShiftedProximalOperators,
   SolverBenchmark
 using Printf
+using Logging
 
 # utils for extracting stats / display table
 modelname(nlp::LSR1Model) = "LSR1"
@@ -92,6 +93,7 @@ function benchmark_table(
         x0 = f.meta.x0,
         subsolver = eval(subsolver),
         subsolver_options = sub_opt,
+        # subsolver_logger= Logging.ConsoleLogger(),
         selected = selected,
       )
     end
diff --git a/benchmarks/tables/svm-table.jl b/benchmarks/tables/svm-table.jl
index d864d9f3..4112803c 100644
--- a/benchmarks/tables/svm-table.jl
+++ b/benchmarks/tables/svm-table.jl
@@ -18,8 +18,9 @@ verbose = 0 #10
 ϵri = 1.0e-6
 maxIter = 1000
 maxIter_inner = 100
+Mmonotone = 10
 options =
-  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -28,8 +29,9 @@ options_nrTR = ROSolverOptions(
   maxIter = maxIter,
   spectral = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options2_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
@@ -37,9 +39,10 @@ options2_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options3 =
-  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+  ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options3_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
@@ -47,10 +50,11 @@ options3_nrTR = ROSolverOptions(
   ϵr = ϵri,
   maxIter = maxIter_inner,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
-options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = Mmonotone)
 options4_nrTR =
-  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false)
+  ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, reduce_TR = false, Mmonotone = Mmonotone)
 options5 = ROSolverOptions(
   ν = ν,
   ϵa = ϵ,
@@ -59,6 +63,7 @@ options5 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = true,
+  Mmonotone = Mmonotone,
 )
 options5_nrTR = ROSolverOptions(
   ν = ν,
@@ -69,6 +74,7 @@ options5_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options6 = ROSolverOptions(
   ν = ν,
@@ -78,6 +84,7 @@ options6 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = false,
+  Mmonotone = Mmonotone,
 )
 options6_nrTR = ROSolverOptions(
   ν = ν,
@@ -88,6 +95,7 @@ options6_nrTR = ROSolverOptions(
   spectral = false,
   psb = false,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 
 solvers = [:R2, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR, :TR, :TR]
@@ -126,7 +134,7 @@ subsolver_options = [
   options4,
   options4_nrTR,
 ] # n'importe lequel si subsolver = :None
-subset = 1:length(solvers)
+subset = 2:length(solvers)
 
 benchmark_table(
   f,
diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl
index d72453e0..adcbab72 100644
--- a/src/TRDH_alg.jl
+++ b/src/TRDH_alg.jl
@@ -122,6 +122,7 @@ function TRDH(
   psb = options.psb
   hess_init_val = (Bk isa UniformScaling) ? Bk.λ : (one(R) / options.ν)
   reduce_TR = options.reduce_TR
+  Mmonotone = options.Mmonotone
 
   local l_bound, u_bound
   has_bnds = false
@@ -184,6 +185,7 @@ function TRDH(
 
   Fobj_hist = zeros(maxIter)
   Hobj_hist = zeros(maxIter)
+  FHobj_hist = zeros(maxIter)
   Complex_hist = zeros(Int, maxIter)
   if verbose > 0
     #! format: off
@@ -219,6 +221,7 @@ function TRDH(
     elapsed_time = time() - start_time
     Fobj_hist[k] = fk
     Hobj_hist[k] = hk
+    FHobj_hist[k] = fk + hk
 
     # model for prox-gradient step to update Δk if ||s|| is too small and ξ1
     φ1(d) = ∇fk' * d
@@ -265,7 +268,14 @@ function TRDH(
     hkn = h(xkn[selected])
     hkn == -Inf && error("nonsmooth term is not proper")
 
-    Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
+    # if Mmonotone == 0
+    #   Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
+    # else
+      Mlast = max(1, k - Mmonotone)
+      fhmax = @views maximum(FHobj_hist[Mlast:k])
+      Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 10 * eps()
+      # Δmod = FHobj_hist[imax] - mk(s) - fkn + max(1, abs(hk)) * 10 * eps()
+    # end
     ξ = hk - mk(s) + max(1, abs(hk)) * 10 * eps()
 
     if !reduce_TR
diff --git a/src/input_struct.jl b/src/input_struct.jl
index e9dd305c..0859d74b 100644
--- a/src/input_struct.jl
+++ b/src/input_struct.jl
@@ -19,6 +19,7 @@ mutable struct ROSolverOptions{R}
   spectral::Bool # for TRDH: use spectral gradient update if true, otherwise DiagonalQN
   psb::Bool # for TRDH with DiagonalQN (spectral = false): use PSB update if true, otherwise Andrei update
   reduce_TR::Bool
+  Mmonotone::Int # last coefs to consider for model decrease computation
 
   function ROSolverOptions{R}(;
     ϵa::R = √eps(R),
@@ -39,6 +40,7 @@ mutable struct ROSolverOptions{R}
     spectral::Bool = false,
     psb::Bool = false,
     reduce_TR::Bool = true,
+    Mmonotone::Int = 0,
   ) where {R <: Real}
     @assert ϵa ≥ 0
     @assert ϵr ≥ 0
@@ -73,6 +75,7 @@ mutable struct ROSolverOptions{R}
       spectral,
       psb,
       reduce_TR,
+      Mmonotone,
     )
   end
 end

From ee77b2dd9e6f28503ea7403f76aa71b37e4a6eb8 Mon Sep 17 00:00:00 2001
From: Geoffroy Leconte <gleconte50@gmail.com>
Date: Tue, 27 Jun 2023 14:04:45 -0400
Subject: [PATCH 4/9] Mmonotone benchmarks

---
 benchmarks/tables/bpdn-constr-table.jl |  1 +
 benchmarks/tables/bpdn-table.jl        |  3 ++-
 benchmarks/tables/fh-table.jl          |  3 ++-
 benchmarks/tables/nnmf-table.jl        |  3 ++-
 benchmarks/tables/regulopt-tables.jl   |  2 +-
 benchmarks/tables/svm-table.jl         |  2 +-
 src/TRDH_alg.jl                        | 14 ++++----------
 7 files changed, 13 insertions(+), 15 deletions(-)

diff --git a/benchmarks/tables/bpdn-constr-table.jl b/benchmarks/tables/bpdn-constr-table.jl
index d3787584..5a2c3e8c 100644
--- a/benchmarks/tables/bpdn-constr-table.jl
+++ b/benchmarks/tables/bpdn-constr-table.jl
@@ -144,6 +144,7 @@ stats = benchmark_table(
   solver_options,
   subsolver_options,
   "BPDN-cstr, M = $Mmonotone",
+<<<<<<< master
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/bpdn-table.jl b/benchmarks/tables/bpdn-table.jl
index 4b07e844..4f5100a2 100644
--- a/benchmarks/tables/bpdn-table.jl
+++ b/benchmarks/tables/bpdn-table.jl
@@ -18,7 +18,7 @@ verbose = 0 # 10
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
-Mmonotone = 0
+Mmonotone = 5
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
@@ -146,6 +146,7 @@ stats = benchmark_table(
   solver_options,
   subsolver_options,
   "BPDN",
+<<<<<<< master
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/fh-table.jl b/benchmarks/tables/fh-table.jl
index ebdfd04c..a20c4b62 100644
--- a/benchmarks/tables/fh-table.jl
+++ b/benchmarks/tables/fh-table.jl
@@ -153,6 +153,7 @@ names, stats = benchmark_table(
   solver_options[subset],
   subsolver_options[subset],
   "FH with ν = $ν, λ = $λ, M = $Mmonotone",
+<<<<<<< master
   random_seed,
   tex = true,
 );
@@ -169,7 +170,7 @@ if display_sol
     row_names = vcat(["True"], names),
     title = "Solution FH",
     formatters = ft_printf("%1.2f"),
-    # backend = Val(:latex),
+    backend = Val(:latex),
   )
 end
 
diff --git a/benchmarks/tables/nnmf-table.jl b/benchmarks/tables/nnmf-table.jl
index d69ba6c4..93307a01 100644
--- a/benchmarks/tables/nnmf-table.jl
+++ b/benchmarks/tables/nnmf-table.jl
@@ -14,7 +14,7 @@ h = NormL0(λ)
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
-Mmonotone = 0
+Mmonotone = 5
 verbose = 0 #10
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
@@ -145,6 +145,7 @@ benchmark_table(
   solver_options[subset],
   subsolver_options[subset],
   "NNMF with m = $m, n = $n, k = $k, ν = $ν, λ = $λ, M = $Mmonotone",
+<<<<<<< master
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/regulopt-tables.jl b/benchmarks/tables/regulopt-tables.jl
index 7cff31e4..530b89d4 100644
--- a/benchmarks/tables/regulopt-tables.jl
+++ b/benchmarks/tables/regulopt-tables.jl
@@ -9,7 +9,7 @@ using NLPModels,
   ShiftedProximalOperators,
   SolverBenchmark
 using Printf
-using Logging
+# using Logging
 
 # utils for extracting stats / display table
 modelname(nlp::LSR1Model) = "LSR1"
diff --git a/benchmarks/tables/svm-table.jl b/benchmarks/tables/svm-table.jl
index 4112803c..4f3973f1 100644
--- a/benchmarks/tables/svm-table.jl
+++ b/benchmarks/tables/svm-table.jl
@@ -18,7 +18,7 @@ verbose = 0 #10
 ϵri = 1.0e-6
 maxIter = 1000
 maxIter_inner = 100
-Mmonotone = 10
+Mmonotone = 5
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl
index adcbab72..4b97fe7f 100644
--- a/src/TRDH_alg.jl
+++ b/src/TRDH_alg.jl
@@ -185,7 +185,7 @@ function TRDH(
 
   Fobj_hist = zeros(maxIter)
   Hobj_hist = zeros(maxIter)
-  FHobj_hist = zeros(maxIter)
+  FHobj_hist = fill!(Vector{R}(undef, Mmonotone), R(-Inf))
   Complex_hist = zeros(Int, maxIter)
   if verbose > 0
     #! format: off
@@ -221,7 +221,7 @@ function TRDH(
     elapsed_time = time() - start_time
     Fobj_hist[k] = fk
     Hobj_hist[k] = hk
-    FHobj_hist[k] = fk + hk
+    Mmonotone > 0 && (FHobj_hist[mod(k-1, Mmonotone) + 1] = fk + hk)
 
     # model for prox-gradient step to update Δk if ||s|| is too small and ξ1
     φ1(d) = ∇fk' * d
@@ -268,14 +268,8 @@ function TRDH(
     hkn = h(xkn[selected])
     hkn == -Inf && error("nonsmooth term is not proper")
 
-    # if Mmonotone == 0
-    #   Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
-    # else
-      Mlast = max(1, k - Mmonotone)
-      fhmax = @views maximum(FHobj_hist[Mlast:k])
-      Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 10 * eps()
-      # Δmod = FHobj_hist[imax] - mk(s) - fkn + max(1, abs(hk)) * 10 * eps()
-    # end
+    fhmax = Mmonotone > 0 ? maximum(FHobj_hist) : fk + hk
+    Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 10 * eps()
     ξ = hk - mk(s) + max(1, abs(hk)) * 10 * eps()
 
     if !reduce_TR

From 97c7b5b1e03ba82b25d086187e25aa2491d9a86e Mon Sep 17 00:00:00 2001
From: Geoffroy Leconte <gleconte50@gmail.com>
Date: Tue, 27 Jun 2023 14:10:58 -0400
Subject: [PATCH 5/9] remove latex outputs

---
 benchmarks/tables/bpdn-constr-table.jl | 1 -
 benchmarks/tables/bpdn-table.jl        | 1 -
 benchmarks/tables/fh-table.jl          | 1 -
 benchmarks/tables/nnmf-table.jl        | 1 -
 benchmarks/tables/svm-table.jl         | 4 ++--
 5 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/benchmarks/tables/bpdn-constr-table.jl b/benchmarks/tables/bpdn-constr-table.jl
index 5a2c3e8c..d3787584 100644
--- a/benchmarks/tables/bpdn-constr-table.jl
+++ b/benchmarks/tables/bpdn-constr-table.jl
@@ -144,7 +144,6 @@ stats = benchmark_table(
   solver_options,
   subsolver_options,
   "BPDN-cstr, M = $Mmonotone",
-<<<<<<< master
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/bpdn-table.jl b/benchmarks/tables/bpdn-table.jl
index 4f5100a2..b629d1dd 100644
--- a/benchmarks/tables/bpdn-table.jl
+++ b/benchmarks/tables/bpdn-table.jl
@@ -146,7 +146,6 @@ stats = benchmark_table(
   solver_options,
   subsolver_options,
   "BPDN",
-<<<<<<< master
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/fh-table.jl b/benchmarks/tables/fh-table.jl
index a20c4b62..39e6b3da 100644
--- a/benchmarks/tables/fh-table.jl
+++ b/benchmarks/tables/fh-table.jl
@@ -153,7 +153,6 @@ names, stats = benchmark_table(
   solver_options[subset],
   subsolver_options[subset],
   "FH with ν = $ν, λ = $λ, M = $Mmonotone",
-<<<<<<< master
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/nnmf-table.jl b/benchmarks/tables/nnmf-table.jl
index 93307a01..b86b199c 100644
--- a/benchmarks/tables/nnmf-table.jl
+++ b/benchmarks/tables/nnmf-table.jl
@@ -145,7 +145,6 @@ benchmark_table(
   solver_options[subset],
   subsolver_options[subset],
   "NNMF with m = $m, n = $n, k = $k, ν = $ν, λ = $λ, M = $Mmonotone",
-<<<<<<< master
   random_seed,
   tex = true,
 );
diff --git a/benchmarks/tables/svm-table.jl b/benchmarks/tables/svm-table.jl
index 4f3973f1..5e1b17d5 100644
--- a/benchmarks/tables/svm-table.jl
+++ b/benchmarks/tables/svm-table.jl
@@ -134,7 +134,7 @@ subsolver_options = [
   options4,
   options4_nrTR,
 ] # n'importe lequel si subsolver = :None
-subset = 2:length(solvers)
+subset = 1:length(solvers)
 
 benchmark_table(
   f,
@@ -150,7 +150,7 @@ benchmark_table(
   random_seed,
   nls_train = nls_train,
   nls_test = nls_test,
-  tex = true,
+  tex = false,
 );
 
 subset = [8, 9, 10, 11, 12, 13, 14]

From b26513bb15f704c5d260d3a7e052cfe27c47b3ce Mon Sep 17 00:00:00 2001
From: Geoffroy Leconte <gleconte50@gmail.com>
Date: Tue, 27 Jun 2023 14:13:12 -0400
Subject: [PATCH 6/9] change Mmonotone bpdn-cstr

---
 benchmarks/tables/bpdn-constr-table.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmarks/tables/bpdn-constr-table.jl b/benchmarks/tables/bpdn-constr-table.jl
index d3787584..c4936e05 100644
--- a/benchmarks/tables/bpdn-constr-table.jl
+++ b/benchmarks/tables/bpdn-constr-table.jl
@@ -17,7 +17,7 @@ verbose = 0 # 10
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
-Mmonotone = 0
+Mmonotone = 5
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(

From 043a0fc10ffecbd1b5c1ae634f29ae121e76d00e Mon Sep 17 00:00:00 2001
From: Geoffroy Leconte <gleconte50@gmail.com>
Date: Mon, 22 Apr 2024 15:17:35 +0200
Subject: [PATCH 7/9] fix denominator

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

diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl
index 4b97fe7f..3d00bb14 100644
--- a/src/TRDH_alg.jl
+++ b/src/TRDH_alg.jl
@@ -270,7 +270,7 @@ function TRDH(
 
     fhmax = Mmonotone > 0 ? maximum(FHobj_hist) : fk + hk
     Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 10 * eps()
-    ξ = hk - mk(s) + max(1, abs(hk)) * 10 * eps()
+    ξ = fhmax - fk - mk(s) + max(1, abs(hk)) * 10 * eps()
 
     if !reduce_TR
 

From 6f417a2fad17500003522846b2e04f22401057c8 Mon Sep 17 00:00:00 2001
From: Geoffroy Leconte <gleconte50@gmail.com>
Date: Tue, 23 Apr 2024 00:27:12 +0200
Subject: [PATCH 8/9] fix iTRDH stop crit

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

diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl
index 3d00bb14..25fb9f61 100644
--- a/src/TRDH_alg.jl
+++ b/src/TRDH_alg.jl
@@ -270,7 +270,8 @@ function TRDH(
 
     fhmax = Mmonotone > 0 ? maximum(FHobj_hist) : fk + hk
     Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 10 * eps()
-    ξ = fhmax - fk - mk(s) + max(1, abs(hk)) * 10 * eps()
+    Δmod = fhmax - (fk + mk(s)) + max(1, abs(hk)) * 10 * eps()
+    ξ = hk - mk(s) + max(1, abs(hk)) * 10 * eps()
 
     if !reduce_TR
 
@@ -290,7 +291,7 @@ function TRDH(
       error("TRDH: failed to compute a step: ξ = $ξ")
     end
 
-    ρk = Δobj / ξ
+    ρk = Δobj / Δmod
 
     TR_stat = (η2 ≤ ρk < Inf) ? "↗" : (ρk < η1 ? "↘" : "=")
 

From 11a5636b9a6efe206ca7eae0f86bea68f97d54be Mon Sep 17 00:00:00 2001
From: Geoffroy Leconte <gleconte50@gmail.com>
Date: Wed, 29 May 2024 15:46:10 +0200
Subject: [PATCH 9/9] nonmonotone tests

---
 benchmarks/tables/bpdn-constr-table.jl | 29 +++++++++++---
 benchmarks/tables/bpdn-table.jl        | 31 ++++++++++++---
 benchmarks/tables/fh-table.jl          | 54 +++++++++++++++++++++-----
 benchmarks/tables/nnmf-table.jl        | 44 +++++++++++++--------
 benchmarks/tables/regulopt-plots.jl    | 21 +++++-----
 benchmarks/tables/regulopt-tables.jl   | 13 ++++---
 benchmarks/tables/svm-table.jl         | 35 +++++++++++++----
 src/input_struct.jl                    |  2 +-
 8 files changed, 169 insertions(+), 60 deletions(-)

diff --git a/benchmarks/tables/bpdn-constr-table.jl b/benchmarks/tables/bpdn-constr-table.jl
index c4936e05..8eb47a0b 100644
--- a/benchmarks/tables/bpdn-constr-table.jl
+++ b/benchmarks/tables/bpdn-constr-table.jl
@@ -17,7 +17,7 @@ verbose = 0 # 10
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
-Mmonotone = 5
+Mmonotone = 10
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
@@ -150,15 +150,34 @@ stats = benchmark_table(
 
 subset = [1, 2, 3, 4, 5, 6, 7]
 
+opt =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 0)
+optPSB =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true, Mmonotone = 0)
+opt5 =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 5)
+sopt = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+solvers2 = [:R2, :TRDH, :TRDH, :TR, :TR]
+subsolvers2 = [:None, :None, :None, :R2, :TRDH]
+solver_options2 = [opt, opt, optPSB, opt, opt]
+subsolver_options2 = [sopt, sopt, sopt, sopt, sopt]
+
 p = benchmark_plot(
   f,
   1:(f.meta.nvar),
   h,
-  solvers[subset],
-  subsolvers[subset],
-  solver_options[subset],
-  subsolver_options[subset],
+  solvers2,
+  subsolvers2,
+  solver_options2,
+  subsolver_options2,
   random_seed;
   xmode = "linear",
   ymode = "log", 
+)
+
+path = raw"C:\Users\Geoffroy Leconte\Documents\doctorat\biblio\papiers\indef-pg\figs\objplots"
+pgfsave(
+  joinpath(path, "bpdn-cstr.tikz"),
+  p;
+  include_preamble = true,
 )
\ No newline at end of file
diff --git a/benchmarks/tables/bpdn-table.jl b/benchmarks/tables/bpdn-table.jl
index b629d1dd..61a4efb5 100644
--- a/benchmarks/tables/bpdn-table.jl
+++ b/benchmarks/tables/bpdn-table.jl
@@ -1,4 +1,5 @@
 include("regulopt-tables.jl")
+include("regulopt-plots.jl")
 
 # model
 random_seed = 1234
@@ -18,7 +19,7 @@ verbose = 0 # 10
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
-Mmonotone = 5
+Mmonotone = 10
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
@@ -150,17 +151,35 @@ stats = benchmark_table(
   tex = true,
 );
 
-subset = [1,2,3,4,5,6, 7]
+opt =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 0)
+optPSB =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true, Mmonotone = 0)
+opt5 =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 5)
+sopt = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+# subset = [1,2,3,4,5,6,7]
+solvers2 = [:R2, :TRDH, :TRDH, :TR, :TR]
+subsolvers2 = [:None, :None, :None, :R2, :TRDH]
+solver_options2 = [opt, opt, optPSB, opt, opt]
+subsolver_options2 = [sopt, sopt, sopt, sopt, sopt]
 
 p = benchmark_plot(
   f,
   1:(f.meta.nvar),
   h,
-  solvers[subset],
-  subsolvers[subset],
-  solver_options[subset],
-  subsolver_options[subset],
+  solvers2,
+  subsolvers2,
+  solver_options2,
+  subsolver_options2,
   random_seed;
   xmode = "linear",
   ymode = "log", 
+)
+
+path = raw"C:\Users\Geoffroy Leconte\Documents\doctorat\biblio\papiers\indef-pg\figs\objplots"
+pgfsave(
+  joinpath(path, "bpdn.tikz"),
+  p;
+  include_preamble = true,
 )
\ No newline at end of file
diff --git a/benchmarks/tables/fh-table.jl b/benchmarks/tables/fh-table.jl
index 39e6b3da..71b8dd7e 100644
--- a/benchmarks/tables/fh-table.jl
+++ b/benchmarks/tables/fh-table.jl
@@ -24,7 +24,7 @@ maxIter_inner = 200 # max iter for subsolver
 ϵ = 1.0e-4
 ϵi = 1.0e-3
 ϵri = 1.0e-6
-Mmonotone = 5
+Mmonotone = 2
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
@@ -122,7 +122,7 @@ solver_options = [
   options,
   options,
   options,
-  options,
+  # options,
 ]
 subsolver_options = [
   options2,
@@ -141,6 +141,7 @@ subsolver_options = [
   options4_nrTR,
 ] # n'importe lequel si subsolver = :None
 subset = 2:length(solvers) # issues with R2 alone
+# subset = [2, 3, 4, 5, 6, 7, 8, 9]
 
 names, stats = benchmark_table(
   f,
@@ -169,22 +170,57 @@ if display_sol
     row_names = vcat(["True"], names),
     title = "Solution FH",
     formatters = ft_printf("%1.2f"),
-    backend = Val(:latex),
+    # backend = Val(:latex),
   )
 end
 
 subset = [8, 9, 10, 11, 12, 13, 14]
 
+if !cstr
+  opt = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 0)
+  optPSB =
+    ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true, Mmonotone = 0)
+  opt5 =
+    ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 5)
+  sopt = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+  soptPSB = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+  soptAndrei = ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+  solvers2 = [:TR, :TR, :TR]
+  subsolvers2 = [:R2, :TRDH, :TRDH]
+  solver_options2 = [opt, opt, opt]
+  subsolver_options2 = [sopt, soptPSB, soptAndrei]
+
+else
+  opt = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 0)
+  optPSB10 =
+    ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true, Mmonotone = 10)
+  opt5 =
+    ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 5)
+  sopt = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+  soptPSB = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+  soptAndrei = ROSolverOptions(spectral = false, psb = false, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+  solvers2 = [:TRDH, :TR, :TR, :TR]
+  subsolvers2 = [:None, :R2, :TRDH, :TRDH]
+  solver_options2 = [optPSB10, opt, opt, opt]
+  subsolver_options2 = [sopt, sopt, soptPSB, soptAndrei]
+end
+
 p = benchmark_plot(
   f,
   1:(f.meta.nvar),
   h,
-  solvers[subset],
-  subsolvers[subset],
-  solver_options[subset],
-  subsolver_options[subset],
+  solvers2,
+  subsolvers2,
+  solver_options2,
+  subsolver_options2,
   random_seed;
-  measured = :grad,
   xmode = "linear",
   ymode = "log", 
-)
\ No newline at end of file
+)
+
+# path = raw"C:\Users\Geoffroy Leconte\Documents\doctorat\biblio\papiers\indef-pg\figs\objplots"
+# pgfsave(
+#   joinpath(path, !cstr ? "fh.tikz" : "fh-cstr.tikz"),
+#   p;
+#   include_preamble = true,
+# )
\ No newline at end of file
diff --git a/benchmarks/tables/nnmf-table.jl b/benchmarks/tables/nnmf-table.jl
index b86b199c..62a0604d 100644
--- a/benchmarks/tables/nnmf-table.jl
+++ b/benchmarks/tables/nnmf-table.jl
@@ -14,7 +14,7 @@ h = NormL0(λ)
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
-Mmonotone = 5
+Mmonotone = 10
 verbose = 0 #10
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
@@ -151,23 +151,35 @@ benchmark_table(
 
 subset = [8, 9, 10, 11, 12, 13, 14]
 
-# p = benchmark_plot(
-#   f,
-#   selected,
-#   h,
-#   solvers[subset],
-#   subsolvers[subset],
-#   solver_options[subset],
-#   subsolver_options[subset],
-#   random_seed;
-#   measured = :grad,
-#   xmode = "linear",
-#   ymode = "log", 
-# )
+opt =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 0)
+optPSB =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true, Mmonotone = 0)
+opt2 =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 2)
+sopt = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+# subset = [1,2,3,4,5,6,7]
+solvers2 = [:TRDH, :TRDH, :TR, :TR]
+subsolvers2 = [:None, :None, :R2, :TRDH]
+solver_options2 = [opt, opt2, opt, opt]
+subsolver_options2 = [sopt, sopt, sopt, sopt]
 
-# path = raw"C:\Users\Geoffroy Leconte\.julia\dev\RegularizedOptimization\benchmarks"
+p = benchmark_plot(
+  f,
+  1:(f.meta.nvar),
+  h,
+  solvers2,
+  subsolvers2,
+  solver_options2,
+  subsolver_options2,
+  random_seed;
+  xmode = "linear",
+  ymode = "log", 
+)
+
+# path = raw"C:\Users\Geoffroy Leconte\Documents\doctorat\biblio\papiers\indef-pg\figs\objplots"
 # pgfsave(
-#   joinpath(path, "test.tikz"),
+#   joinpath(path, "nnmf.tikz"),
 #   p;
 #   include_preamble = true,
 # )
\ No newline at end of file
diff --git a/benchmarks/tables/regulopt-plots.jl b/benchmarks/tables/regulopt-plots.jl
index 49919929..34684794 100644
--- a/benchmarks/tables/regulopt-plots.jl
+++ b/benchmarks/tables/regulopt-plots.jl
@@ -26,13 +26,16 @@ function options_str(
   subsolver::Symbol,
 )
   if solver == :TRDH
-    out_str = !options.spectral ? (options.psb ? "-PSB" : "-Andrei") : "-Spec"
-    out_str = (options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+    out_str = (options.reduce_TR) ? string(solver) : string("i", solver)
+    (options.Mmonotone != 0) && (out_str *= string("-", options.Mmonotone)) 
+    out_str *= !options.spectral ? (options.psb ? "-PSB" : "-Andrei") : "-Spec"
   elseif solver == :TR && subsolver == :TRDH
-    out_str = !subsolver_options.spectral ? (subsolver_options.psb ? "-PSB" : "-Andrei") : "-Spec"
-    out_str = (subsolver_options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+    out_str = (subsolver_options.reduce_TR) ? string(solver, "-", subsolver) : string(solver, "-i", subsolver)
+    (subsolver_options.Mmonotone != 0) && (out_str *= string("-", subsolver_options.Mmonotone))
+    out_str *= !subsolver_options.spectral ? (subsolver_options.psb ? "-PSB" : "-Andrei") : "-Spec"
   else
-    out_str = ""
+    out_str = string(solver)
+    (subsolver == :None) || (out_str *= string("-", subsolver))
   end
   return out_str
 end
@@ -46,12 +49,12 @@ function benchmark_plot(
   solver_options,
   subsolver_options,
   random_seed::Int;
-  measured::Symbol = :obj, # set to :grad to eval grad
+  measured::Symbol = :grad, # set to :grad to eval grad
   xmode::String = "log",
   ymode::String = "log",
 )
   solver_names = [
-    "$(solver)$(subsolvername(subsolver))$(options_str(opt, solver, subsolver_opt, subsolver))"
+    "$(options_str(opt, solver, subsolver_opt, subsolver))"
     for (solver, opt, subsolver, subsolver_opt) in
     zip(solvers, solver_options, subsolvers, subsolver_options)
   ]
@@ -93,13 +96,13 @@ function benchmark_plot(
     )
   end
 
-  colors = distinguishable_colors(n_solvers)
+  colors = distinguishable_colors(n_solvers, [RGB(1,1,1)], dropseed=true)
   l_plots = [@pgf Plot({color = colors[i]}, coords[i]) for i in 1:n_solvers]
   
   @pgf Axis(
     {
       xlabel = "iterations",
-      ylabel = L"$(f + h)(x_k)$",
+      ylabel = L"$(f + h)(x_k) - (f + h)^*$",
       ymode = ymode,
       xmode = xmode,
       no_markers,
diff --git a/benchmarks/tables/regulopt-tables.jl b/benchmarks/tables/regulopt-tables.jl
index 530b89d4..34cf11a5 100644
--- a/benchmarks/tables/regulopt-tables.jl
+++ b/benchmarks/tables/regulopt-tables.jl
@@ -24,13 +24,14 @@ function options_str(
   subsolver::Symbol,
 )
   if solver == :TRDH
-    out_str = !options.spectral ? (options.psb ? "-PSB" : "-Andrei") : "-Spec"
-    out_str = (options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+    out_str = (options.reduce_TR) ? string(solver) : string("i", solver)
+    out_str *= !options.spectral ? (options.psb ? "-PSB" : "-Andrei") : "-Spec"
   elseif solver == :TR && subsolver == :TRDH
-    out_str = !subsolver_options.spectral ? (subsolver_options.psb ? "-PSB" : "-Andrei") : "-Spec"
-    out_str = (subsolver_options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+    out_str = (subsolver_options.reduce_TR) ? string(solver, "-", subsolver) : string(solver, "-i", subsolver)
+    out_str *= !subsolver_options.spectral ? (subsolver_options.psb ? "-PSB" : "-Andrei") : "-Spec"
   else
-    out_str = ""
+    out_str = string(solver)
+    (subsolver == :None) || (out_str *= string("-", subsolver))
   end
   return out_str
 end
@@ -66,7 +67,7 @@ function benchmark_table(
   nls_test::Union{Nothing, AbstractNLSModel} = nothing, # for SVM
 )
   solver_names = [
-    "$(solver)$(subsolvername(subsolver))$(options_str(opt, solver, subsolver_opt, subsolver))"
+    "$(options_str(opt, solver, subsolver_opt, subsolver))"
     for (solver, opt, subsolver, subsolver_opt) in
     zip(solvers, solver_options, subsolvers, subsolver_options)
   ]
diff --git a/benchmarks/tables/svm-table.jl b/benchmarks/tables/svm-table.jl
index 5e1b17d5..77099dd4 100644
--- a/benchmarks/tables/svm-table.jl
+++ b/benchmarks/tables/svm-table.jl
@@ -18,7 +18,7 @@ verbose = 0 #10
 ϵri = 1.0e-6
 maxIter = 1000
 maxIter_inner = 100
-Mmonotone = 5
+Mmonotone = 10
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = Mmonotone)
 options_nrTR = ROSolverOptions(
@@ -150,21 +150,40 @@ benchmark_table(
   random_seed,
   nls_train = nls_train,
   nls_test = nls_test,
-  tex = false,
+  tex = true,
 );
 
 subset = [8, 9, 10, 11, 12, 13, 14]
 
+opt =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 0)
+optPSB =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true, Mmonotone = 0)
+opt5 =
+  ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true, Mmonotone = 5)
+sopt = ROSolverOptions(spectral = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+soptPSB = ROSolverOptions(spectral = false, psb = true, ϵa = ϵi, ϵr = ϵri, maxIter = maxIter_inner, Mmonotone = 0)
+solvers2 = [:R2, :TRDH, :TRDH, :TR, :TR]
+subsolvers2 = [:None, :None, :None, :R2, :TRDH]
+solver_options2 = [opt, opt, opt5, opt, opt]
+subsolver_options2 = [sopt, sopt, sopt, sopt, soptPSB]
+
 p = benchmark_plot(
   f,
   1:(f.meta.nvar),
   h,
-  solvers[subset],
-  subsolvers[subset],
-  solver_options[subset],
-  subsolver_options[subset],
+  solvers2,
+  subsolvers2,
+  solver_options2,
+  subsolver_options2,
   random_seed;
-  measured = :grad,
   xmode = "linear",
   ymode = "log", 
-)
\ No newline at end of file
+)
+
+# path = raw"C:\Users\Geoffroy Leconte\Documents\doctorat\biblio\papiers\indef-pg\figs\objplots"
+# pgfsave(
+#   joinpath(path, "svm.tikz"),
+#   p;
+#   include_preamble = true,
+# )
\ No newline at end of file
diff --git a/src/input_struct.jl b/src/input_struct.jl
index 0859d74b..e6126b98 100644
--- a/src/input_struct.jl
+++ b/src/input_struct.jl
@@ -40,7 +40,7 @@ mutable struct ROSolverOptions{R}
     spectral::Bool = false,
     psb::Bool = false,
     reduce_TR::Bool = true,
-    Mmonotone::Int = 0,
+    Mmonotone::Int = 5,
   ) where {R <: Real}
     @assert ϵa ≥ 0
     @assert ϵr ≥ 0