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 3c4207a2..8eb47a0b 100644
--- a/benchmarks/tables/bpdn-constr-table.jl
+++ b/benchmarks/tables/bpdn-constr-table.jl
@@ -13,12 +13,13 @@ 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
+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 = ϵ,
@@ -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,41 @@ stats = benchmark_table(
   subsolvers,
   solver_options,
   subsolver_options,
-  "BPDN-cstr",
+  "BPDN-cstr, M = $Mmonotone",
   random_seed,
-  tex = false,
+  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)
+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,
+  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 dfa5768c..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
@@ -14,12 +15,13 @@ 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
+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 +30,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 +40,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 +51,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 +64,7 @@ options5 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = true,
+  Mmonotone = Mmonotone,
 )
 options5_nrTR = ROSolverOptions(
   ν = ν,
@@ -69,6 +75,7 @@ options5_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options6 = ROSolverOptions(
   ν = ν,
@@ -78,6 +85,7 @@ options6 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = false,
+  Mmonotone = Mmonotone,
 )
 options6_nrTR = ROSolverOptions(
   ν = ν,
@@ -88,6 +96,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]
@@ -139,5 +148,38 @@ stats = benchmark_table(
   subsolver_options,
   "BPDN",
   random_seed,
-  tex = false,
+  tex = true,
 );
+
+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,
+  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 4584e2e0..71b8dd7e 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
@@ -6,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)
@@ -23,8 +24,9 @@ maxIter_inner = 200 # max iter for subsolver
 ϵ = 1.0e-4
 ϵi = 1.0e-3
 ϵri = 1.0e-6
+Mmonotone = 2
 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 = ϵ,
@@ -33,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,
@@ -42,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,
@@ -52,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 = ϵ,
@@ -64,6 +69,7 @@ options5 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = true,
+  Mmonotone = Mmonotone,
 )
 options5_nrTR = ROSolverOptions(
   ν = ν,
@@ -74,6 +80,7 @@ options5_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options6 = ROSolverOptions(
   ν = ν,
@@ -83,6 +90,7 @@ options6 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = false,
+  Mmonotone = Mmonotone,
 )
 options6_nrTR = ROSolverOptions(
   ν = ν,
@@ -93,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]
@@ -113,7 +122,7 @@ solver_options = [
   options,
   options,
   options,
-  options,
+  # options,
 ]
 subsolver_options = [
   options2,
@@ -132,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,
@@ -143,9 +153,9 @@ names, stats = benchmark_table(
   subsolvers[subset],
   solver_options[subset],
   subsolver_options[subset],
-  "FH with ν = $ν, λ = $λ",
+  "FH with ν = $ν, λ = $λ, M = $Mmonotone",
   random_seed,
-  tex = false,
+  tex = true,
 );
 
 if display_sol
@@ -163,3 +173,54 @@ if display_sol
     # 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,
+  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, !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 20cabfc6..62a0604d 100644
--- a/benchmarks/tables/nnmf-table.jl
+++ b/benchmarks/tables/nnmf-table.jl
@@ -1,11 +1,12 @@
 include("regulopt-tables.jl")
+include("regulopt-plots.jl")
 
 random_seed = 1234
 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
@@ -13,9 +14,10 @@ h = NormL0(λ)
 ϵri = 1.0e-6
 maxIter = 500
 maxIter_inner = 100
+Mmonotone = 10
 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 = ϵ,
@@ -24,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,
@@ -33,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,
@@ -43,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 = ϵ,
@@ -55,6 +60,7 @@ options5 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = true,
+  Mmonotone = Mmonotone,
 )
 options5_nrTR = ROSolverOptions(
   ν = ν,
@@ -65,6 +71,7 @@ options5_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options6 = ROSolverOptions(
   ν = ν,
@@ -74,6 +81,7 @@ options6 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = false,
+  Mmonotone = Mmonotone,
 )
 options6_nrTR = ROSolverOptions(
   ν = ν,
@@ -84,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]
@@ -122,6 +131,7 @@ subsolver_options = [
   options4,
   options4_nrTR,
 ] # n'importe lequel si subsolver = :None
+
 subset = 1:length(solvers)
 
 benchmark_table(
@@ -134,7 +144,42 @@ 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 = 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)
+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]
+
+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, "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
new file mode 100644
index 00000000..34684794
--- /dev/null
+++ b/benchmarks/tables/regulopt-plots.jl
@@ -0,0 +1,118 @@
+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.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.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 = string(solver)
+    (subsolver == :None) || (out_str *= string("-", subsolver))
+  end
+  return out_str
+end
+
+function benchmark_plot(
+  f::AbstractNLPModel,
+  selected,
+  h,
+  solvers,
+  subsolvers,
+  solver_options,
+  subsolver_options,
+  random_seed::Int;
+  measured::Symbol = :grad, # set to :grad to eval grad
+  xmode::String = "log",
+  ymode::String = "log",
+)
+  solver_names = [
+    "$(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, [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) - (f + h)^*$",
+      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..34cf11a5 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"
@@ -23,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
@@ -65,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)
   ]
@@ -74,6 +76,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)
@@ -91,9 +94,11 @@ function benchmark_table(
         x0 = f.meta.x0,
         subsolver = eval(subsolver),
         subsolver_options = sub_opt,
+        # subsolver_logger= Logging.ConsoleLogger(),
         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 ed937972..77099dd4 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
@@ -17,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 = ϵ,
@@ -27,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,
@@ -36,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,
@@ -46,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 = ϵ,
@@ -58,6 +63,7 @@ options5 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = true,
+  Mmonotone = Mmonotone,
 )
 options5_nrTR = ROSolverOptions(
   ν = ν,
@@ -68,6 +74,7 @@ options5_nrTR = ROSolverOptions(
   spectral = false,
   psb = true,
   reduce_TR = false,
+  Mmonotone = Mmonotone,
 )
 options6 = ROSolverOptions(
   ν = ν,
@@ -77,6 +84,7 @@ options6 = ROSolverOptions(
   maxIter = maxIter,
   spectral = false,
   psb = false,
+  Mmonotone = Mmonotone,
 )
 options6_nrTR = ROSolverOptions(
   ν = ν,
@@ -87,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]
@@ -141,5 +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,
+  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, "svm.tikz"),
+#   p;
+#   include_preamble = true,
+# )
\ No newline at end of file
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..25fb9f61 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 = fill!(Vector{R}(undef, Mmonotone), R(-Inf))
   Complex_hist = zeros(Int, maxIter)
   if verbose > 0
     #! format: off
@@ -205,8 +207,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)
@@ -219,6 +221,7 @@ function TRDH(
     elapsed_time = time() - start_time
     Fobj_hist[k] = fk
     Hobj_hist[k] = 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
@@ -265,7 +268,9 @@ 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()
+    fhmax = Mmonotone > 0 ? maximum(FHobj_hist) : fk + hk
+    Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 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
@@ -286,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 ? "↘" : "=")
 
@@ -315,7 +320,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 +330,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
 
diff --git a/src/input_struct.jl b/src/input_struct.jl
index e9dd305c..e6126b98 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 = 5,
   ) where {R <: Real}
     @assert ϵa ≥ 0
     @assert ϵr ≥ 0
@@ -73,6 +75,7 @@ mutable struct ROSolverOptions{R}
       spectral,
       psb,
       reduce_TR,
+      Mmonotone,
     )
   end
 end