arc_bmark.jl

using LinearAlgebra
using Logging

using JLD2

using CUTEst
using JSOSolvers
using NLPModels
using SolverBenchmark
using SolverTools

using Printf
using DataFrames
include("run_solver_by_prob.jl")

using ARCTR
using State
using Stopping

srcpath = joinpath(@__DIR__, "..", "src")
probnames = ["ROSENBR", "WOODS", "PENALTY1"]
problems = (CUTEstModel(prob) for prob ∈ probnames)
# problems = (CUTEstModel(probname) for probname in CUTEst.select(max_var=5000, max_con=0, only_free_var=true))

const atol = 1.0e-5
const rtol = 1.0e-6
const max_time = 3600.0

for (fn, fnsimple) ∈ [(:ARCqKOp, :arcqk), (:ST_TROp, :st)]
  @eval begin

    function $fnsimple(prob)
      stop = NLPStopping(prob,
                         Stopping.unconstrained,
                         NLPAtX(prob.meta.x0),
                         meta = StoppingMeta(atol=atol,
                                             rtol=rtol,
                                             ),
                         )
      stop.meta.rtol_x = 0.0  # rtol = rtol_x is not reasonable
      stop.meta.rtol_f = 0.0  # no stalled...
      stop.meta.max_iter = typemax(Int)
      stop.meta.max_time = max_time
      stop.meta.max_eval = typemax(Int)

      Δt = @timed begin
        final_state, stop.meta.optimal = eval($fn)(prob, stop)
      end

      status = :unknown
      stop.meta.optimal && (status = :first_order)
      stop.meta.unbounded && (status = :unbounded)
      stop.meta.stalled && (status = :stalled)
      stop.meta.tired && (status = :max_eval)  # should be more specific

      return GenericExecutionStats(status, prob,
                                   solution = final_state.x,
                                   iter = stop.meta.nb_of_stop,  # not quite the number of iterations!
                                   primal_feas = 0.0,
                                   dual_feas = norm(final_state.gx),
                                   objective = final_state.fx,
                                   elapsed_time = Δt[2],
                                  )
    end

  end
end

solvers = Dict{Symbol,Function}(
               :ARCqK => arcqk,
               :SteihaugToint => st,
               # :Trunk => prob -> trunk(prob, max_time=max_time, atol=atol, rtol=rtol),
               # :TRON => prob -> tron(prob, max_time=max_time, atol=atol, rtol=rtol),
               )

stats = solve_problems(solvers, problems)

# t_bmark = @timed begin
#   stats = bmark_solvers(solvers, problems)
# end
#
# # save DataFrames to disk
# jldopen("st_stats.jld2", "w") do file
#   file["st"] = stats[:SteihaugToint]
# end
# jldopen("arcqkop_stats.jld2", "w") do file
#   file["arcqkop"] = stats[:ARCqK]
# end
# jldopen("trunk_stats.jld2", "w") do file
#   file["trunk"] = stats[:ARCqK]
# end
# jldopen("tron_stats.jld2", "w") do file
#   file["tron"] = stats[:ARCqK]
# end
#
# # to load a DataFrame from disk, use
# # arcqkop_stats = jldopen("benchmark/arcqkop_stats.jld2") do file
# #   file["arcqkop"]
# # end
#
# # extract categories of problems based on number of variables
# const nvar_S = 100
# const nvar_M = 1000
# const nvar_L = 10_000
# stats_S = Dict(key => stats[key][stats[key].nvar .≤ nvar_S, :] for key in keys(stats))
# stats_M = Dict(key => stats[key][nvar_S .< stats[key].nvar .≤ nvar_M, :] for key in keys(stats))
# stats_L = Dict(key => stats[key][nvar_M .< stats[key].nvar .≤ nvar_L, :] for key in keys(stats))
# stats_XL = Dict(key => stats[key][nvar_L .< stats[key].nvar, :] for key in keys(stats))
#
# # useful to count the number of problems with each final status
# function count_unique(X)
#   vals = Dict{eltype(X),Int}()
#   for x ∈ X
#     vals[x] = x ∈ keys(vals) ? (vals[x] + 1) : 1
#   end
#   vals
# end
#
# for solver ∈ keys(solvers)
#   display(count_unique(stats[solver].status))
# end
#
# # performance measures
solved(df) = df.status .== :first_order
costnames = ["time",
             "objective evals",
             "gradient evals",
             "hessian-vector products",
             "obj + grad + hprod"]
costs = [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,
         df -> .!solved(df) .* Inf .+ df.neval_obj .+ df.neval_grad .+ df.neval_hprod]

# plot and save performance profiles
p = profile_solvers(stats, costs, costnames)
# Plots.pdf(p, "arcqk_vs_st_all.pdf")
#
# p_S = profile_solvers(stats_S, costs, costnames)
# Plots.pdf(p_S, "arcqk_vs_st_small.pdf")
# p_M = profile_solvers(stats_M, costs, costnames)
# Plots.pdf(p_M, "arcqk_vs_st_medium.pdf")
# p_L = profile_solvers(stats_L, costs, costnames)
# Plots.pdf(p_L, "arcqk_vs_st_large.pdf")