Add scripts to reproduce paper results (#96)

Author: tmigot

paper/2022-05-16_ST_TROp_ARCqKOpShift05_cutest_277_1000000.jld2

@@ -0,0 +1,11 @@
+[deps]
+CUTEst = "1b53aba6-35b6-5f92-a507-53c67d53f819"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SolverBenchmark = "581a75fa-a23a-52d0-a590-d6201de2218a"
+
+[compat]
+DataFrames = "1.3"

paper/Project.toml

@@ -0,0 +1,16 @@
+# Scalable Adaptive Cubic Regularization Methods
+
+This folder contains a set of scripts to reproduce the results in
+
+> Dussault, J.-P., Migot, T. & Orban, D. (2023).
+> Scalable adaptive cubic regularization methods.
+> Mathematical Programming.
+> [10.1007/s10107-023-02007-6](https://doi.org/10.1007/s10107-023-02007-6)
+
+There are three scripts:
+
+- `make_problem_cutest_list.jl`: generate a data file with the list of problems solved;
+- `benchmark_arctr.jl`: run the benchmark comparing ST_TR and ARCqK, and store the result in a JLD2 file;
+- `figure_solvers.jl`: generates a set of plots.
+
+The result file `2022-05-16_ST_TROp_ARCqKOpShift05_cutest_277_1000000.jld2` can be accessed without reproducing the benchmark.

paper/README.md

@@ -0,0 +1,78 @@
+using Pkg;
+Pkg.activate("");
+using AdaptiveRegularization, CUTEst, Dates, NLPModels, SolverBenchmark
+
+nmax = 1000000
+problems = readlines("paper/list_problems_$nmax.dat")
+cutest_problems = (CUTEstModel(p) for p in problems)
+
+unsolved_problems = [
+  "CYCLOOCFLS",
+  "SBRYBND",
+  "BA-L49LS",
+  "NONMSQRT",
+  "BA-L73LS",
+  "INDEF",
+  "BA-L21LS",
+  "FLETCBV3",
+  "FLETCHBV",
+  "BA-L16LS",
+  "BA-L52LS",
+]
+
+max_time = 3600.0
+max_ev = typemax(Int)
+max_iter = typemax(Int)
+atol = 1e-5
+rtol = 1e-6
+
+eL = 10.0
+eU = 20.0
+ψ = 10 ^ (1 // 2)
+shifts = 10.0 .^ (collect(-eL:0.5:eU))
+
+ζ = 0.5
+
+γ₁ = 0.1
+γ₂ = 5
+η₁ = 0.1
+η₂ = 0.75
+TR = TrustRegion(
+    10;
+    acceptance_threshold = η₁,
+    increase_threshold = η₂,
+    increase_factor = γ₂,
+    decrease_factor = γ₁,
+)
+
+solvers = Dict(
+    :ARCqKOp =>
+        nlp -> ARCqKOp(
+            nlp,
+            TR = TR,
+            shifts = shifts,
+            verbose = false,
+            atol = atol,
+            rtol = rtol,
+            max_time = max_time,
+            max_iter = max_iter,
+            ζ = ζ,
+        ),
+    :ST_TROp =>
+        nlp -> ST_TROp(
+            nlp,
+            TR = TR,
+            verbose = false,
+            atol = atol,
+            rtol = rtol,
+            max_time = max_time,
+            max_iter = max_iter,
+            ζ = ζ,
+        ),
+)
+
+# Skipping unsolved problems to solve time :-)
+stats = bmark_solvers(solvers, cutest_problems, skipif = prob -> (prob.meta.name in unsolved_problems))
+
+using JLD2
+@save "paper/$(today())_$(prod(String.(keys(solvers)) .* "_"))cutest_$(string(length(problems)))_$(nmax).jld2" stats

paper/benchmark_arctr.jl

@@ -0,0 +1,61 @@
+using Pkg;
+Pkg.activate("paper");
+using JLD2, Plots, SolverBenchmark, DataFrames
+
+name = "2022-05-16_ST_TROp_ARCqKOpShift05_cutest_277_1000000"
+@load "paper/$name.jld2" stats
+
+solved(df) = (df.status .== :first_order) .| (df.status .== :unbounded)
+
+for solver in keys(stats)
+    open("paper/$(name)_result_$(solver).dat", "w") do io
+        print(
+            io,
+            stats[solver][
+                !,
+                [
+                    :name,
+                    :nvar,
+                    # :ncon,
+                    :status,
+                    :objective,
+                    :elapsed_time,
+                    :iter,
+                    # :primal_feas,
+                    :dual_feas,
+                    :neval_obj,
+                    :neval_grad,
+                    :neval_hprod,
+                    :neval_hess,
+                ],
+            ],
+        )
+    end
+end
+
+nmins = [0, 100, 1000, 10000]
+for nmin in nmins
+    # Same figure with minimum number of variables
+    stats2 = copy(stats)
+    for solver in keys(stats)
+        stats2[solver] = stats[solver][stats[solver].nvar.>=nmin, :]
+    end
+
+    nb_problems = length(stats2[first(keys(stats))][!, :name])
+
+    # Figures comparing two results:
+    costs_all = [
+        df -> .!solved(df) * Inf + df.elapsed_time,
+        df -> .!solved(df) * Inf + df.neval_obj,
+        df -> .!solved(df) * Inf + df.neval_grad,
+        df -> .!solved(df) * Inf + df.neval_hprod,
+    ]
+    costnames_all = [
+        "elapsed time",
+        "objective evals",
+        "gradient evals",
+        "hessian-vector products",
+    ]
+    p = profile_solvers(stats2, costs_all, costnames_all, height = 400, width = 400, margin=5Plots.mm)
+    png(p, "paper/$(name)_all($(nb_problems))_min_$(nmin)")
+end