From e2dfbbbbec44d71d0b2ee14422b607b293cf1a53 Mon Sep 17 00:00:00 2001 From: Ayush1299 Date: Wed, 2 Apr 2025 10:17:05 +0530 Subject: [PATCH 1/2] Fix: Resolved solver_tracker issue and converted .jmd to .jl - Fixed solver_tracker issue that was causing UndefVarError in onlinear_solver_23_tests.jl. - Successfully converted onlinear_solver_23_tests.jmd to .jl format. - Imported necessary dependencies to ensure smooth execution. Issues Faced: - The benchmark still doesn't display plots. - Investigated the final cells but couldn't find an explicit issue causing the missing plots. - Need further debugging to trace where the plotting logic is breaking. Next Steps: - Check if the plots are being suppressed or failing silently. - Verify dependencies and recent commits for changes affecting visualization. --- .gitignore | 3 +- .../nonlinear_solver_23_tests.jmd | 2 + .../output/nonlinear_solver_23_tests.jl | 561 ++++++++++++++++++ 3 files changed, 565 insertions(+), 1 deletion(-) create mode 100644 benchmarks/NonlinearProblem/output/nonlinear_solver_23_tests.jl diff --git a/.gitignore b/.gitignore index a74426f43..ca05434c9 100644 --- a/.gitignore +++ b/.gitignore @@ -22,4 +22,5 @@ pdf/ notebook/ markdown/ -.vscode \ No newline at end of file +.vscode +.qodo diff --git a/benchmarks/NonlinearProblem/nonlinear_solver_23_tests.jmd b/benchmarks/NonlinearProblem/nonlinear_solver_23_tests.jmd index ed7577df5..fef2b540e 100644 --- a/benchmarks/NonlinearProblem/nonlinear_solver_23_tests.jmd +++ b/benchmarks/NonlinearProblem/nonlinear_solver_23_tests.jmd @@ -14,6 +14,8 @@ These benchmarks compares the runtime and error for a range of nonlinear solvers - NLsolveJL's [Anderson acceleration](https://docs.sciml.ai/NonlinearSolve/stable/api/nlsolve/#Solver-API) (`NLsolveJL(method=:anderson)`). - Sundials's [Newton-Krylov](https://docs.sciml.ai/NonlinearSolve/stable/api/sundials/#Solver-API) method (`KINSOL()`). +solver_tracker = []; + Furthermore, for NonlinearSolve.jl's Newton Raphson method we try the following Line Search options (in addition to the default): - `HagerZhang` - `MoreThuente` diff --git a/benchmarks/NonlinearProblem/output/nonlinear_solver_23_tests.jl b/benchmarks/NonlinearProblem/output/nonlinear_solver_23_tests.jl new file mode 100644 index 000000000..cbc1e540d --- /dev/null +++ b/benchmarks/NonlinearProblem/output/nonlinear_solver_23_tests.jl @@ -0,0 +1,561 @@ +# Nonlinear Solver 23 Test Problems +# author Torkel Loman & Avik Pal + +# These benchmarks compares the runtime and error for a range of nonlinear solvers. The problems are a standard set of problems as described [here](https://people.sc.fsu.edu/~jburkardt/m_src/test_nonlin/test_nonlin.html). The solvers are implemented in [NonlinearProblemLibrary.jl](https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/lib/NonlinearProblemLibrary/src/NonlinearProblemLibrary.jl), where you can find the problem function declarations. For each problem we test the following solvers: +# - NonlinearSolve.jl's [Newton Raphson](https://docs.sciml.ai/NonlinearSolve/stable/api/nonlinearsolve/#NonlinearSolve.NewtonRaphson) method (`NewtonRaphson()`). +# - NonlinearSolve.jl's [Trust Region](https://docs.sciml.ai/NonlinearSolve/stable/api/nonlinearsolve/#NonlinearSolve.TrustRegion) method (`TrustRegion()`). +# - NonlinearSolve.jl's Levenberg-Marquardt method (`LevenbergMarquardt()`). +# - MINPACK's [Modified Powell](https://docs.sciml.ai/NonlinearSolve/stable/api/minpack/#NonlinearSolveMINPACK.CMINPACK) method (`CMINPACK(method=:hybr)`). +# - MINPACK's [Levenberg-Marquardt](https://docs.sciml.ai/NonlinearSolve/stable/api/minpack/#NonlinearSolveMINPACK.CMINPACK) method (`CMINPACK(method=:lm)`). +# - NLsolveJL's [Newton Raphson](https://docs.sciml.ai/NonlinearSolve/stable/api/nlsolve/#Solver-API) (`NLsolveJL(method=:newton)`). +# - NLsolveJL's [Trust Region](https://docs.sciml.ai/NonlinearSolve/stable/api/nlsolve/#Solver-API) (`NLsolveJL()`). +# - NLsolveJL's [Anderson acceleration](https://docs.sciml.ai/NonlinearSolve/stable/api/nlsolve/#Solver-API) (`NLsolveJL(method=:anderson)`). +# - Sundials's [Newton-Krylov](https://docs.sciml.ai/NonlinearSolve/stable/api/sundials/#Solver-API) method (`KINSOL()`). + +# Furthermore, for NonlinearSolve.jl's Newton Raphson method we try the following Line Search options (in addition to the default): +# - `HagerZhang` +# - `MoreThuente` +# - `BackTracking` + +# # and for NonlinearSolve.jl's Trust Region we try the following Radius Update schemes (in addition to the default): +# - `NLsolve` +# - `NocedalWright` +# - `Hei` +# - `Yuan` +# - `Bastin` +# - `Fan` +# and finally for NonlinearSolve.jl's Levenberg-Marquardt method why try using both the default `α_geodesic` value (`0.75`) and a modified value (`0.5`), and also with and without setting the `CholeskyFactorization` linear solver. + +# For each benchmarked problem, the second, third, and fourth plots compares the performance of NonlinearSolve's Newton Raphson, Trust Region, and Levenberg-Marquardt methods, respectively. The first plot compares the best methods from each of these categories to the various methods available from other packages. At the end of the benchmarks, we print a summary table of which solvers succeeded for which problems. + +# Fetch required packages. + +```julia +using NonlinearSolve, LinearSolve, StaticArrays, Sundials, Setfield, + BenchmarkTools, LinearAlgebra, DiffEqDevTools, NonlinearProblemLibrary, CairoMakie, + RecursiveFactorization +import PolyesterForwardDiff, MINPACK, NLsolve + +const RUS = RadiusUpdateSchemes; +``` + +# Declare the benchmarked solvers (and their names and plotting options). + +```julia + +solvers_all = [ + (; pkg = :nonlinearsolve, type = :general, name = "Default PolyAlg.", solver = Dict(:alg => FastShortcutNonlinearPolyalg(; u0_len = 10))), + (; pkg = :nonlinearsolve, type = :NR, name = "Newton Raphson", solver = Dict(:alg => NewtonRaphson())), + (; pkg = :nonlinearsolve, type = :NR, name = "NR (HagerZhang)", solver = Dict(:alg => NewtonRaphson(; linesearch = HagerZhang()))), + (; pkg = :nonlinearsolve, type = :NR, name = "NR (MoreThuente)", solver = Dict(:alg => NewtonRaphson(; linesearch = MoreThuente()))), + (; pkg = :nonlinearsolve, type = :NR, name = "NR (BackTracking)", solver = Dict(:alg => NewtonRaphson(; linesearch = BackTracking()))), + (; pkg = :nonlinearsolve, type = :TR, name = "Trust Region", solver = Dict(:alg => TrustRegion())), + (; pkg = :nonlinearsolve, type = :TR, name = "TR (NLsolve Update)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.NLsolve))), + (; pkg = :nonlinearsolve, type = :TR, name = "TR (Nocedal Wright)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.NocedalWright))), + (; pkg = :nonlinearsolve, type = :TR, name = "TR (Hei)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.Hei))), + (; pkg = :nonlinearsolve, type = :TR, name = "TR (Yuan)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.Yuan))), + (; pkg = :nonlinearsolve, type = :TR, name = "TR (Bastin)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.Bastin))), + (; pkg = :nonlinearsolve, type = :TR, name = "TR (Fan)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.Fan))), + (; pkg = :nonlinearsolve, type = :LM, name = "Levenberg-Marquardt", solver = Dict(:alg => LevenbergMarquardt(; linsolve = QRFactorization()))), + (; pkg = :nonlinearsolve, type = :LM, name = "LM with Cholesky", solver = Dict(:alg => LevenbergMarquardt(; linsolve = CholeskyFactorization()))), + (; pkg = :nonlinearsolve, type = :LM, name = "LM (α_geodesic=0.5)", solver = Dict(:alg => LevenbergMarquardt(; linsolve = QRFactorization(), α_geodesic=0.5))), + (; pkg = :nonlinearsolve, type = :LM, name = "LM (α_geodesic=0.5) Chol.", solver = Dict(:alg => LevenbergMarquardt(; linsolve = CholeskyFactorization(), α_geodesic=0.5))), + (; pkg = :nonlinearsolve, type = :LM, name = "LM (no Accln.)", solver = Dict(:alg => LevenbergMarquardt(; linsolve = QRFactorization(), disable_geodesic = Val(true)))), + (; pkg = :nonlinearsolve, type = :LM, name = "LM (no Accln.) Chol.", solver = Dict(:alg => LevenbergMarquardt(; linsolve = CholeskyFactorization(), disable_geodesic = Val(true)))), + (; pkg = :nonlinearsolve, type = :general, name = "Pseudo Transient", solver = Dict(:alg => PseudoTransient(; alpha_initial=10.0))), + (; pkg = :wrapper, type = :general, name = "Powell [MINPACK]", solver = Dict(:alg => CMINPACK(; method=:hybr))), + (; pkg = :wrapper, type = :general, name = "LM [MINPACK]", solver = Dict(:alg => CMINPACK(; method=:lm))), + (; pkg = :wrapper, type = :general, name = "NR [NLsolve.jl]", solver = Dict(:alg => NLsolveJL(; method=:newton))), + (; pkg = :wrapper, type = :general, name = "TR [NLsolve.jl]", solver = Dict(:alg => NLsolveJL())), + (; pkg = :wrapper, type = :general, name = "NR [Sundials]", solver = Dict(:alg => KINSOL())), + (; pkg = :wrapper, type = :general, name = "NR LineSearch [Sundials]", solver = Dict(:alg => KINSOL(; globalization_strategy=:LineSearch))) +]; + +solver_tracker = []; +wp_general_tracker = []; +``` + +# Sets tolerances. + +```julia +abstols = 1.0 ./ 10.0 .^ (4:12) +reltols = 1.0 ./ 10.0 .^ (4:12); +``` + +# Prepares various helper functions for benchmarking a specific problem. + +```julia + +function benchmark_problem!(prob_name; solver_tracker=solver_tracker) + + prob = nlprob_23_testcases[prob_name] + + successful_solvers = filter(Base.Fix1(check_solver, prob), solvers_all); + push!(solver_tracker, prob_name => successful_solvers); + + solvers_NR = filter(s -> s.type==:NR, successful_solvers) + solvers_TR = filter(s -> s.type==:TR, successful_solvers) + solvers_LM = filter(s -> s.type==:LM, successful_solvers) + wp_NR = WorkPrecisionSet(prob.prob, abstols, reltols, getfield.(solvers_NR, :solver); + names=getfield.(solvers_NR, :name), numruns=100, error_estimate=:l∞, + maxiters=1000, + termination_condition = NonlinearSolve.AbsNormTerminationMode(Base.Fix1(maximum, abs))) + wp_TR = WorkPrecisionSet(prob.prob, abstols, reltols, getfield.(solvers_TR, :solver); + names=getfield.(solvers_TR, :name), numruns=100, error_estimate=:l∞, + maxiters=1000, + termination_condition = NonlinearSolve.AbsNormTerminationMode(Base.Fix1(maximum, abs))) + wp_LM = WorkPrecisionSet(prob.prob, abstols, reltols, getfield.(solvers_LM, :solver); + names=getfield.(solvers_LM, :name), numruns=100, error_estimate=:l∞, + maxiters=1000, + termination_condition = NonlinearSolve.AbsNormTerminationMode(Base.Fix1(maximum, abs))) + + solvers_general = filter(s -> s.type==:general, successful_solvers) + add_solver!(solvers_general, nothing, solvers_TR, wp_TR) + add_solver!(solvers_general, nothing, solvers_LM, wp_LM) + add_solver!(solvers_general, nothing, solvers_NR, wp_NR) + + wp_general = WorkPrecisionSet(prob.prob, abstols, reltols, + getfield.(solvers_general, :solver); names=getfield.(solvers_general, :name), + numruns=100, error_estimate=:l∞, maxiters=1000) + + push!(wp_general_tracker, prob_name => wp_general) + + fig = plot_collective_benchmark(prob_name, wp_general, wp_NR, wp_TR, wp_LM) + + save(replace(lowercase(prob_name), " " => "_") * "_wpd.svg", fig) + + return fig +end + +function check_solver(prob, solver) + try + sol = solve(prob.prob, solver.solver[:alg]; abstol=1e-8, reltol=1e-8, + maxiters=1000000, + termination_condition=NonlinearSolve.AbsNormTerminationMode(Base.Fix1(maximum, abs))) + if norm(sol.resid, Inf) < 1e-6 + Base.printstyled("[Info] Solver $(solver.name) returned retcode $(sol.retcode) \ + with an residual norm = $(norm(sol.resid, Inf)).\n"; color=:green) + return true + else + Base.printstyled("[Warn] Solver $(solver.name) had a very large residual \ + (norm = $(norm(sol.resid, Inf))).\n"; color=:red) + return false + end + WorkPrecisionSet(prob.prob, [1e-4, 1e-12], [1e-4, 1e-12], [solver.solver]; + names=[solver.name], numruns=5, error_estimate=:l∞, maxiters=1000) + catch e + Base.printstyled("[Warn] Solver $(solver.name) threw an error: $e.\n"; color=:red) + return false + end + return true +end + +function add_solver!(solvers_general, selected_solver_name, additional_solver_set, wp) + if isnothing(selected_solver_name) + isempty(wp.wps) && return + selected_idx = argmin(median.(getfield.(wp.wps, :times))) + else + selected_idx = findfirst(s -> s.name==selected_solver_name, additional_solver_set) + isnothing(selected_solver) && error("The $(selected_solver_name) was designated to \ + be added to the general solver set, however, it seemed to fail on this \ + problem.") + end + isnothing(selected_idx) || + pushfirst!(solvers_general, additional_solver_set[selected_idx]) +end; +``` + +# Plotting related helper functions. + +```julia +__log10_zero(x) = ifelse(iszero(x), -100, log10(x)) +Makie.inverse_transform(::typeof(__log10_zero)) = exp10 +Makie.defaultlimits(::typeof(__log10_zero)) = Makie.defaultlimits(log10) +Makie.defined_interval(::typeof(__log10_zero)) = 0.0..Inf + +function Makie.get_minor_tickvalues(i::IntervalsBetween, scale::typeof(__log10_zero), + tickvalues, vmin, vmax) + return [] +end + +tickformatter(values) = map(values) do v + e = log10(v) + if isinteger(e) && e == -100 + return rich("10", superscript("-∞")) + end + sup = isinteger(e) ? Int(e) : round(e; digits=2) + return rich("10", superscript(string(sup))) +end + +function __filter_nearzero((ticks, ticklabels)) + if first(ticks) ≈ 1e-100 + idxs = findall(x -> x ≈ 1e-100 || x ≥ 10^-40, ticks) + return ticks[idxs], ticklabels[idxs] + end + return ticks, ticklabels +end + +function plot_collective_benchmark(prob_name, wp_general, wp_NR, wp_TR, wp_LM) + LINESTYLES = Dict(:nonlinearsolve => :solid, :simplenonlinearsolve => :dash, + :wrapper => :dot) + ASPECT_RATIO = 0.7 + WIDTH = 1400 + HEIGHT = round(Int, WIDTH * ASPECT_RATIO) + STROKEWIDTH = 2.5 + + colors = cgrad(:seaborn_bright, length(solvers_all); categorical = true) + cycle = Cycle([:marker], covary = true) + plot_theme = Theme(Lines = (; cycle), Scatter = (; cycle)) + + fig = with_theme(plot_theme) do + fig = Figure(; size = (WIDTH, HEIGHT)) + axs = [] + + xmin, xmax, ymin, ymax = Inf, -Inf, Inf, -Inf + + for i in 1:2, j in 1:2 + wp = (wp_general, wp_NR, wp_TR, wp_LM)[2 * (i - 1) + j] + + ax = Axis(fig[i + 1, j], ylabel = j == 1 ? L"Time $\mathbf{(s)}$" : "", + xlabelsize = 22, ylabelsize = 22, + xlabel = i == 2 ? L"Error: $\mathbf{||f(u^\ast)||_\infty}$" : "", + xscale = __log10_zero, yscale = __log10_zero, + xtickwidth = STROKEWIDTH, + ytickwidth = STROKEWIDTH, spinewidth = STROKEWIDTH, + xticklabelsize = 20, yticklabelsize = 20, + xticklabelsvisible = i == 2, yticklabelsvisible = j == 1, + xticksvisible = i == 2, yticksvisible = j == 1,) + push!(axs, ax) + + ls = [] + scs = [] + + for wpᵢ in wp.wps + idx = findfirst(s -> s.name == wpᵢ.name, solvers_all) + errs = getindex.(wpᵢ.errors, :l∞) + times = wpᵢ.times + + emin, emax = extrema(errs) + tmin, tmax = extrema(times) + emin < xmin && (xmin = emin) + emax > xmax && (xmax = emax) + tmin < ymin && (ymin = tmin) + tmax > ymax && (ymax = tmax) + + l = lines!(ax, errs, times; color = colors[idx], linewidth = 5, + linestyle = LINESTYLES[solvers_all[idx].pkg], alpha = 0.8, + label = wpᵢ.name) + sc = scatter!(ax, errs, times; color = colors[idx], markersize = 16, + strokewidth = 2, marker = Cycled(idx), alpha = 0.8, label = wpᵢ.name) + push!(ls, l) + push!(scs, sc) + end + + legend_title = ("", "Newton Raphson", "Trust Region", "Levenberg-Marquardt")[2 * (i - 1) + j] + + Legend(fig[ifelse(i == 1, 1, 4), j], [[l, sc] for (l, sc) in zip(ls, scs)], + [wpᵢ.name for wpᵢ in wp.wps], legend_title; + framevisible=true, framewidth = STROKEWIDTH, + nbanks = 3, labelsize = 16, titlesize = 16, + tellheight = true, tellwidth = false, patchsize = (40.0f0, 20.0f0)) + end + + linkaxes!(axs...) + + xmin = max(xmin, 10^-100) + + xticks = __filter_nearzero(Makie.get_ticks(LogTicks(WilkinsonTicks(10; k_min = 5)), + __log10_zero, tickformatter, xmin, xmax)) + yticks = __filter_nearzero(Makie.get_ticks(LogTicks(WilkinsonTicks(10; k_min = 5)), + __log10_zero, tickformatter, ymin, ymax)) + + foreach(axs) do ax + ax.xticks = xticks + ax.yticks = yticks + end + + fig[0, :] = Label(fig, "Work-Precision Diagram for $(prob_name)", + fontsize = 24, tellwidth = false, font = :bold) + + fig + end + + return fig +end +``` + +# We here run benchmarks for each of the 23 models. + +### Problem 1 (Generalized Rosenbrock function) + +```julia +benchmark_problem!("Generalized Rosenbrock function") +``` + +### Problem 2 (Powell singular function) + +```julia +benchmark_problem!("Powell singular function") +``` + +### Problem 3 (Powell badly scaled function) + +```julia +benchmark_problem!("Powell badly scaled function") +``` + +### Problem 4 (Wood function) + +```julia +benchmark_problem!("Wood function") +``` + +### Problem 5 (Helical valley function) + +```julia +benchmark_problem!("Helical valley function") +``` + +### Problem 6 (Watson function) + +```julia +benchmark_problem!("Watson function") +``` + +### Problem 7 (Chebyquad function) + +```julia +benchmark_problem!("Chebyquad function") +``` + +### Problem 8 (Brown almost linear function) + +```julia +benchmark_problem!("Brown almost linear function") +``` + +### Problem 9 (Discrete boundary value function) + +```julia +benchmark_problem!("Discrete boundary value function") +``` + +### Problem 10 (Discrete integral equation function) + +```julia +benchmark_problem!("Discrete integral equation function") +``` + +### Problem 11 (Trigonometric function) + +```julia +benchmark_problem!("Trigonometric function") +``` + +### Problem 12 (Variably dimensioned function) + +```julia +benchmark_problem!("Variably dimensioned function") +``` + +### Problem 13 (Broyden tridiagonal function) + +```julia +benchmark_problem!("Broyden tridiagonal function") +``` + +### Problem 14 (Broyden banded function) + +```julia +benchmark_problem!("Broyden banded function") +``` + +### Problem 15 (Hammarling 2 by 2 matrix square root problem) + +```julia +benchmark_problem!("Hammarling 2 by 2 matrix square root problem") +``` + +### Problem 16 (Hammarling 3 by 3 matrix square root problem) + +```julia +benchmark_problem!("Hammarling 3 by 3 matrix square root problem") +``` + +### Problem 17 (Dennis and Schnabel 2 by 2 example) + +```julia +benchmark_problem!("Dennis and Schnabel 2 by 2 example") +``` + +### Problem 18 (Sample problem 18) + +```julia +benchmark_problem!("Sample problem 18") +``` + +### Problem 19 (Sample problem 19) + +```julia +benchmark_problem!("Sample problem 19") +``` + +### Problem 20 (Scalar problem f(x) = x(x - 5)^2) + +```julia +benchmark_problem!("Scalar problem f(x) = x(x - 5)^2") +``` + +### Problem 21 (Freudenstein-Roth function) + +```julia +benchmark_problem!("Freudenstein-Roth function") +``` + +### Problem 22 (Boggs function) + +```julia +benchmark_problem!("Boggs function") +``` + +### Problem 23 (Chandrasekhar function) + +```julia +benchmark_problem!("Chandrasekhar function") +``` + +# Summary of successful solvers + +# Finally, we print a summary of which solvers successfully solved which problems. + +```julia +solver_successes = [(solver.name in getfield.(prob[2], :name)) ? "O" : "X" for prob in solver_tracker, solver in solvers_all] +total_successes = [sum(solver_successes[:,i] .== "O") for i in 1:length(solvers_all)] +solver_outcomes = vcat(total_successes', solver_successes); +``` + +```julia; wrap = false; results = "md2html" +using PrettyTables +io = IOBuffer() +println(io, "```@raw html") +pretty_table(io, solver_outcomes; backend = Val(:html), header = getfield.(solvers_all, :name), alignment=:c) +println(io, "```") +Base.Text(String(take!(io))) +``` + +# Summary of General Solver Performance on All Problems + +```julia +fig = begin + LINESTYLES = Dict(:nonlinearsolve => :solid, :simplenonlinearsolve => :dash, + :wrapper => :dot) + ASPECT_RATIO = 1 + WIDTH = 1800 + HEIGHT = round(Int, WIDTH * ASPECT_RATIO) + STROKEWIDTH = 2.5 + + colors = cgrad(:seaborn_bright, length(solvers_all); categorical = true) + cycle = Cycle([:marker], covary = true) + plot_theme = Theme(Lines = (; cycle), Scatter = (; cycle)) + + with_theme(plot_theme) do + fig = Figure(; size = (WIDTH, HEIGHT)) + + axs = Matrix{Any}(undef, 5, 5) + + ls = [] + scs = [] + labels = [] + solver_times = [] + + for i in 1:5, j in 1:5 + idx = 5 * (i - 1) + j + + idx > length(wp_general_tracker) && break + + prob_name, wp = wp_general_tracker[idx] + + ax = Axis(fig[i, j], + xscale = __log10_zero, yscale = __log10_zero, + xtickwidth = STROKEWIDTH, + ytickwidth = STROKEWIDTH, spinewidth = STROKEWIDTH, + title = prob_name, titlegap = 10, + xticklabelsize = 16, yticklabelsize = 16) + + xmin, xmax, ymin, ymax = Inf, -Inf, Inf, -Inf + + for wpᵢ in wp.wps + idx = findfirst(s -> s.name == wpᵢ.name, solvers_all) + errs = getindex.(wpᵢ.errors, :l∞) + times = wpᵢ.times + + emin, emax = extrema(errs) + tmin, tmax = extrema(times) + emin < xmin && (xmin = emin) + emax > xmax && (xmax = emax) + tmin < ymin && (ymin = tmin) + tmax > ymax && (ymax = tmax) + + l = lines!(ax, errs, times; color = colors[idx], linewidth = 5, + linestyle = LINESTYLES[solvers_all[idx].pkg], alpha = 0.8, + label = wpᵢ.name) + sc = scatter!(ax, errs, times; color = colors[idx], markersize = 16, + strokewidth = 2, marker = Cycled(idx), alpha = 0.8, label = wpᵢ.name) + + if wpᵢ.name ∉ labels + push!(ls, l) + push!(scs, sc) + push!(labels, wpᵢ.name) + end + + if wpᵢ.name ∈ first.(solver_times) + idxi = findfirst(x -> first(x) == wpᵢ.name, solver_times) + push!(solver_times[idxi][2], median(times) / length(wp.prob.u0)) + else + push!(solver_times, wpᵢ.name => [median(times) / length(wp.prob.u0)]) + end + end + + xmin = max(xmin, 10^-100) + + xticks = __filter_nearzero(Makie.get_ticks(LogTicks(WilkinsonTicks(5; k_min = 3)), + __log10_zero, tickformatter, xmin, xmax)) + yticks = __filter_nearzero(Makie.get_ticks(LogTicks(WilkinsonTicks(5; k_min = 3)), + __log10_zero, tickformatter, ymin, ymax)) + + ax.xticks = xticks + ax.yticks = yticks + end + + ordering = sortperm(median.(last.(solver_times))) + + fig[0, :] = Label(fig, "Work-Precision Diagram for 23 Test Problems", + fontsize = 24, tellwidth = false, font = :bold) + + fig[:, 0] = Label(fig, "Time (s)", fontsize = 20, tellheight = false, font = :bold, + rotation = π / 2) + fig[end + 1, :] = Label(fig, + L"Error: $\mathbf{||f(u^\ast)||_\infty}$", + fontsize = 20, tellwidth = false, font = :bold) + + Legend(fig[5, 4:5], [[l, sc] for (l, sc) in zip(ls[ordering], scs[ordering])], + labels[ordering], "Successful Solvers"; + framevisible=true, framewidth = STROKEWIDTH, orientation = :horizontal, + titlesize = 20, nbanks = 9, labelsize = 20, halign = :center, + tellheight = false, tellwidth = false, patchsize = (40.0f0, 20.0f0)) + + return fig + end +end +``` + +```julia +save("summary_wp_23test_problems.svg", fig) +``` + +```julia, echo = false +using SciMLBenchmarks +SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file]) +``` + +# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl From 7d1cd9866b1a50a247db9391ba535300550d7287 Mon Sep 17 00:00:00 2001 From: AyushRajput <134182912+ayush2281@users.noreply.github.com> Date: Tue, 8 Apr 2025 13:27:58 +0530 Subject: [PATCH 2/2] Delete benchmarks/NonlinearProblem/output/nonlinear_solver_23_tests.jl --- .../output/nonlinear_solver_23_tests.jl | 561 ------------------ 1 file changed, 561 deletions(-) delete mode 100644 benchmarks/NonlinearProblem/output/nonlinear_solver_23_tests.jl diff --git a/benchmarks/NonlinearProblem/output/nonlinear_solver_23_tests.jl b/benchmarks/NonlinearProblem/output/nonlinear_solver_23_tests.jl deleted file mode 100644 index cbc1e540d..000000000 --- a/benchmarks/NonlinearProblem/output/nonlinear_solver_23_tests.jl +++ /dev/null @@ -1,561 +0,0 @@ -# Nonlinear Solver 23 Test Problems -# author Torkel Loman & Avik Pal - -# These benchmarks compares the runtime and error for a range of nonlinear solvers. The problems are a standard set of problems as described [here](https://people.sc.fsu.edu/~jburkardt/m_src/test_nonlin/test_nonlin.html). The solvers are implemented in [NonlinearProblemLibrary.jl](https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/lib/NonlinearProblemLibrary/src/NonlinearProblemLibrary.jl), where you can find the problem function declarations. For each problem we test the following solvers: -# - NonlinearSolve.jl's [Newton Raphson](https://docs.sciml.ai/NonlinearSolve/stable/api/nonlinearsolve/#NonlinearSolve.NewtonRaphson) method (`NewtonRaphson()`). -# - NonlinearSolve.jl's [Trust Region](https://docs.sciml.ai/NonlinearSolve/stable/api/nonlinearsolve/#NonlinearSolve.TrustRegion) method (`TrustRegion()`). -# - NonlinearSolve.jl's Levenberg-Marquardt method (`LevenbergMarquardt()`). -# - MINPACK's [Modified Powell](https://docs.sciml.ai/NonlinearSolve/stable/api/minpack/#NonlinearSolveMINPACK.CMINPACK) method (`CMINPACK(method=:hybr)`). -# - MINPACK's [Levenberg-Marquardt](https://docs.sciml.ai/NonlinearSolve/stable/api/minpack/#NonlinearSolveMINPACK.CMINPACK) method (`CMINPACK(method=:lm)`). -# - NLsolveJL's [Newton Raphson](https://docs.sciml.ai/NonlinearSolve/stable/api/nlsolve/#Solver-API) (`NLsolveJL(method=:newton)`). -# - NLsolveJL's [Trust Region](https://docs.sciml.ai/NonlinearSolve/stable/api/nlsolve/#Solver-API) (`NLsolveJL()`). -# - NLsolveJL's [Anderson acceleration](https://docs.sciml.ai/NonlinearSolve/stable/api/nlsolve/#Solver-API) (`NLsolveJL(method=:anderson)`). -# - Sundials's [Newton-Krylov](https://docs.sciml.ai/NonlinearSolve/stable/api/sundials/#Solver-API) method (`KINSOL()`). - -# Furthermore, for NonlinearSolve.jl's Newton Raphson method we try the following Line Search options (in addition to the default): -# - `HagerZhang` -# - `MoreThuente` -# - `BackTracking` - -# # and for NonlinearSolve.jl's Trust Region we try the following Radius Update schemes (in addition to the default): -# - `NLsolve` -# - `NocedalWright` -# - `Hei` -# - `Yuan` -# - `Bastin` -# - `Fan` -# and finally for NonlinearSolve.jl's Levenberg-Marquardt method why try using both the default `α_geodesic` value (`0.75`) and a modified value (`0.5`), and also with and without setting the `CholeskyFactorization` linear solver. - -# For each benchmarked problem, the second, third, and fourth plots compares the performance of NonlinearSolve's Newton Raphson, Trust Region, and Levenberg-Marquardt methods, respectively. The first plot compares the best methods from each of these categories to the various methods available from other packages. At the end of the benchmarks, we print a summary table of which solvers succeeded for which problems. - -# Fetch required packages. - -```julia -using NonlinearSolve, LinearSolve, StaticArrays, Sundials, Setfield, - BenchmarkTools, LinearAlgebra, DiffEqDevTools, NonlinearProblemLibrary, CairoMakie, - RecursiveFactorization -import PolyesterForwardDiff, MINPACK, NLsolve - -const RUS = RadiusUpdateSchemes; -``` - -# Declare the benchmarked solvers (and their names and plotting options). - -```julia - -solvers_all = [ - (; pkg = :nonlinearsolve, type = :general, name = "Default PolyAlg.", solver = Dict(:alg => FastShortcutNonlinearPolyalg(; u0_len = 10))), - (; pkg = :nonlinearsolve, type = :NR, name = "Newton Raphson", solver = Dict(:alg => NewtonRaphson())), - (; pkg = :nonlinearsolve, type = :NR, name = "NR (HagerZhang)", solver = Dict(:alg => NewtonRaphson(; linesearch = HagerZhang()))), - (; pkg = :nonlinearsolve, type = :NR, name = "NR (MoreThuente)", solver = Dict(:alg => NewtonRaphson(; linesearch = MoreThuente()))), - (; pkg = :nonlinearsolve, type = :NR, name = "NR (BackTracking)", solver = Dict(:alg => NewtonRaphson(; linesearch = BackTracking()))), - (; pkg = :nonlinearsolve, type = :TR, name = "Trust Region", solver = Dict(:alg => TrustRegion())), - (; pkg = :nonlinearsolve, type = :TR, name = "TR (NLsolve Update)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.NLsolve))), - (; pkg = :nonlinearsolve, type = :TR, name = "TR (Nocedal Wright)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.NocedalWright))), - (; pkg = :nonlinearsolve, type = :TR, name = "TR (Hei)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.Hei))), - (; pkg = :nonlinearsolve, type = :TR, name = "TR (Yuan)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.Yuan))), - (; pkg = :nonlinearsolve, type = :TR, name = "TR (Bastin)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.Bastin))), - (; pkg = :nonlinearsolve, type = :TR, name = "TR (Fan)", solver = Dict(:alg => TrustRegion(; radius_update_scheme = RUS.Fan))), - (; pkg = :nonlinearsolve, type = :LM, name = "Levenberg-Marquardt", solver = Dict(:alg => LevenbergMarquardt(; linsolve = QRFactorization()))), - (; pkg = :nonlinearsolve, type = :LM, name = "LM with Cholesky", solver = Dict(:alg => LevenbergMarquardt(; linsolve = CholeskyFactorization()))), - (; pkg = :nonlinearsolve, type = :LM, name = "LM (α_geodesic=0.5)", solver = Dict(:alg => LevenbergMarquardt(; linsolve = QRFactorization(), α_geodesic=0.5))), - (; pkg = :nonlinearsolve, type = :LM, name = "LM (α_geodesic=0.5) Chol.", solver = Dict(:alg => LevenbergMarquardt(; linsolve = CholeskyFactorization(), α_geodesic=0.5))), - (; pkg = :nonlinearsolve, type = :LM, name = "LM (no Accln.)", solver = Dict(:alg => LevenbergMarquardt(; linsolve = QRFactorization(), disable_geodesic = Val(true)))), - (; pkg = :nonlinearsolve, type = :LM, name = "LM (no Accln.) Chol.", solver = Dict(:alg => LevenbergMarquardt(; linsolve = CholeskyFactorization(), disable_geodesic = Val(true)))), - (; pkg = :nonlinearsolve, type = :general, name = "Pseudo Transient", solver = Dict(:alg => PseudoTransient(; alpha_initial=10.0))), - (; pkg = :wrapper, type = :general, name = "Powell [MINPACK]", solver = Dict(:alg => CMINPACK(; method=:hybr))), - (; pkg = :wrapper, type = :general, name = "LM [MINPACK]", solver = Dict(:alg => CMINPACK(; method=:lm))), - (; pkg = :wrapper, type = :general, name = "NR [NLsolve.jl]", solver = Dict(:alg => NLsolveJL(; method=:newton))), - (; pkg = :wrapper, type = :general, name = "TR [NLsolve.jl]", solver = Dict(:alg => NLsolveJL())), - (; pkg = :wrapper, type = :general, name = "NR [Sundials]", solver = Dict(:alg => KINSOL())), - (; pkg = :wrapper, type = :general, name = "NR LineSearch [Sundials]", solver = Dict(:alg => KINSOL(; globalization_strategy=:LineSearch))) -]; - -solver_tracker = []; -wp_general_tracker = []; -``` - -# Sets tolerances. - -```julia -abstols = 1.0 ./ 10.0 .^ (4:12) -reltols = 1.0 ./ 10.0 .^ (4:12); -``` - -# Prepares various helper functions for benchmarking a specific problem. - -```julia - -function benchmark_problem!(prob_name; solver_tracker=solver_tracker) - - prob = nlprob_23_testcases[prob_name] - - successful_solvers = filter(Base.Fix1(check_solver, prob), solvers_all); - push!(solver_tracker, prob_name => successful_solvers); - - solvers_NR = filter(s -> s.type==:NR, successful_solvers) - solvers_TR = filter(s -> s.type==:TR, successful_solvers) - solvers_LM = filter(s -> s.type==:LM, successful_solvers) - wp_NR = WorkPrecisionSet(prob.prob, abstols, reltols, getfield.(solvers_NR, :solver); - names=getfield.(solvers_NR, :name), numruns=100, error_estimate=:l∞, - maxiters=1000, - termination_condition = NonlinearSolve.AbsNormTerminationMode(Base.Fix1(maximum, abs))) - wp_TR = WorkPrecisionSet(prob.prob, abstols, reltols, getfield.(solvers_TR, :solver); - names=getfield.(solvers_TR, :name), numruns=100, error_estimate=:l∞, - maxiters=1000, - termination_condition = NonlinearSolve.AbsNormTerminationMode(Base.Fix1(maximum, abs))) - wp_LM = WorkPrecisionSet(prob.prob, abstols, reltols, getfield.(solvers_LM, :solver); - names=getfield.(solvers_LM, :name), numruns=100, error_estimate=:l∞, - maxiters=1000, - termination_condition = NonlinearSolve.AbsNormTerminationMode(Base.Fix1(maximum, abs))) - - solvers_general = filter(s -> s.type==:general, successful_solvers) - add_solver!(solvers_general, nothing, solvers_TR, wp_TR) - add_solver!(solvers_general, nothing, solvers_LM, wp_LM) - add_solver!(solvers_general, nothing, solvers_NR, wp_NR) - - wp_general = WorkPrecisionSet(prob.prob, abstols, reltols, - getfield.(solvers_general, :solver); names=getfield.(solvers_general, :name), - numruns=100, error_estimate=:l∞, maxiters=1000) - - push!(wp_general_tracker, prob_name => wp_general) - - fig = plot_collective_benchmark(prob_name, wp_general, wp_NR, wp_TR, wp_LM) - - save(replace(lowercase(prob_name), " " => "_") * "_wpd.svg", fig) - - return fig -end - -function check_solver(prob, solver) - try - sol = solve(prob.prob, solver.solver[:alg]; abstol=1e-8, reltol=1e-8, - maxiters=1000000, - termination_condition=NonlinearSolve.AbsNormTerminationMode(Base.Fix1(maximum, abs))) - if norm(sol.resid, Inf) < 1e-6 - Base.printstyled("[Info] Solver $(solver.name) returned retcode $(sol.retcode) \ - with an residual norm = $(norm(sol.resid, Inf)).\n"; color=:green) - return true - else - Base.printstyled("[Warn] Solver $(solver.name) had a very large residual \ - (norm = $(norm(sol.resid, Inf))).\n"; color=:red) - return false - end - WorkPrecisionSet(prob.prob, [1e-4, 1e-12], [1e-4, 1e-12], [solver.solver]; - names=[solver.name], numruns=5, error_estimate=:l∞, maxiters=1000) - catch e - Base.printstyled("[Warn] Solver $(solver.name) threw an error: $e.\n"; color=:red) - return false - end - return true -end - -function add_solver!(solvers_general, selected_solver_name, additional_solver_set, wp) - if isnothing(selected_solver_name) - isempty(wp.wps) && return - selected_idx = argmin(median.(getfield.(wp.wps, :times))) - else - selected_idx = findfirst(s -> s.name==selected_solver_name, additional_solver_set) - isnothing(selected_solver) && error("The $(selected_solver_name) was designated to \ - be added to the general solver set, however, it seemed to fail on this \ - problem.") - end - isnothing(selected_idx) || - pushfirst!(solvers_general, additional_solver_set[selected_idx]) -end; -``` - -# Plotting related helper functions. - -```julia -__log10_zero(x) = ifelse(iszero(x), -100, log10(x)) -Makie.inverse_transform(::typeof(__log10_zero)) = exp10 -Makie.defaultlimits(::typeof(__log10_zero)) = Makie.defaultlimits(log10) -Makie.defined_interval(::typeof(__log10_zero)) = 0.0..Inf - -function Makie.get_minor_tickvalues(i::IntervalsBetween, scale::typeof(__log10_zero), - tickvalues, vmin, vmax) - return [] -end - -tickformatter(values) = map(values) do v - e = log10(v) - if isinteger(e) && e == -100 - return rich("10", superscript("-∞")) - end - sup = isinteger(e) ? Int(e) : round(e; digits=2) - return rich("10", superscript(string(sup))) -end - -function __filter_nearzero((ticks, ticklabels)) - if first(ticks) ≈ 1e-100 - idxs = findall(x -> x ≈ 1e-100 || x ≥ 10^-40, ticks) - return ticks[idxs], ticklabels[idxs] - end - return ticks, ticklabels -end - -function plot_collective_benchmark(prob_name, wp_general, wp_NR, wp_TR, wp_LM) - LINESTYLES = Dict(:nonlinearsolve => :solid, :simplenonlinearsolve => :dash, - :wrapper => :dot) - ASPECT_RATIO = 0.7 - WIDTH = 1400 - HEIGHT = round(Int, WIDTH * ASPECT_RATIO) - STROKEWIDTH = 2.5 - - colors = cgrad(:seaborn_bright, length(solvers_all); categorical = true) - cycle = Cycle([:marker], covary = true) - plot_theme = Theme(Lines = (; cycle), Scatter = (; cycle)) - - fig = with_theme(plot_theme) do - fig = Figure(; size = (WIDTH, HEIGHT)) - axs = [] - - xmin, xmax, ymin, ymax = Inf, -Inf, Inf, -Inf - - for i in 1:2, j in 1:2 - wp = (wp_general, wp_NR, wp_TR, wp_LM)[2 * (i - 1) + j] - - ax = Axis(fig[i + 1, j], ylabel = j == 1 ? L"Time $\mathbf{(s)}$" : "", - xlabelsize = 22, ylabelsize = 22, - xlabel = i == 2 ? L"Error: $\mathbf{||f(u^\ast)||_\infty}$" : "", - xscale = __log10_zero, yscale = __log10_zero, - xtickwidth = STROKEWIDTH, - ytickwidth = STROKEWIDTH, spinewidth = STROKEWIDTH, - xticklabelsize = 20, yticklabelsize = 20, - xticklabelsvisible = i == 2, yticklabelsvisible = j == 1, - xticksvisible = i == 2, yticksvisible = j == 1,) - push!(axs, ax) - - ls = [] - scs = [] - - for wpᵢ in wp.wps - idx = findfirst(s -> s.name == wpᵢ.name, solvers_all) - errs = getindex.(wpᵢ.errors, :l∞) - times = wpᵢ.times - - emin, emax = extrema(errs) - tmin, tmax = extrema(times) - emin < xmin && (xmin = emin) - emax > xmax && (xmax = emax) - tmin < ymin && (ymin = tmin) - tmax > ymax && (ymax = tmax) - - l = lines!(ax, errs, times; color = colors[idx], linewidth = 5, - linestyle = LINESTYLES[solvers_all[idx].pkg], alpha = 0.8, - label = wpᵢ.name) - sc = scatter!(ax, errs, times; color = colors[idx], markersize = 16, - strokewidth = 2, marker = Cycled(idx), alpha = 0.8, label = wpᵢ.name) - push!(ls, l) - push!(scs, sc) - end - - legend_title = ("", "Newton Raphson", "Trust Region", "Levenberg-Marquardt")[2 * (i - 1) + j] - - Legend(fig[ifelse(i == 1, 1, 4), j], [[l, sc] for (l, sc) in zip(ls, scs)], - [wpᵢ.name for wpᵢ in wp.wps], legend_title; - framevisible=true, framewidth = STROKEWIDTH, - nbanks = 3, labelsize = 16, titlesize = 16, - tellheight = true, tellwidth = false, patchsize = (40.0f0, 20.0f0)) - end - - linkaxes!(axs...) - - xmin = max(xmin, 10^-100) - - xticks = __filter_nearzero(Makie.get_ticks(LogTicks(WilkinsonTicks(10; k_min = 5)), - __log10_zero, tickformatter, xmin, xmax)) - yticks = __filter_nearzero(Makie.get_ticks(LogTicks(WilkinsonTicks(10; k_min = 5)), - __log10_zero, tickformatter, ymin, ymax)) - - foreach(axs) do ax - ax.xticks = xticks - ax.yticks = yticks - end - - fig[0, :] = Label(fig, "Work-Precision Diagram for $(prob_name)", - fontsize = 24, tellwidth = false, font = :bold) - - fig - end - - return fig -end -``` - -# We here run benchmarks for each of the 23 models. - -### Problem 1 (Generalized Rosenbrock function) - -```julia -benchmark_problem!("Generalized Rosenbrock function") -``` - -### Problem 2 (Powell singular function) - -```julia -benchmark_problem!("Powell singular function") -``` - -### Problem 3 (Powell badly scaled function) - -```julia -benchmark_problem!("Powell badly scaled function") -``` - -### Problem 4 (Wood function) - -```julia -benchmark_problem!("Wood function") -``` - -### Problem 5 (Helical valley function) - -```julia -benchmark_problem!("Helical valley function") -``` - -### Problem 6 (Watson function) - -```julia -benchmark_problem!("Watson function") -``` - -### Problem 7 (Chebyquad function) - -```julia -benchmark_problem!("Chebyquad function") -``` - -### Problem 8 (Brown almost linear function) - -```julia -benchmark_problem!("Brown almost linear function") -``` - -### Problem 9 (Discrete boundary value function) - -```julia -benchmark_problem!("Discrete boundary value function") -``` - -### Problem 10 (Discrete integral equation function) - -```julia -benchmark_problem!("Discrete integral equation function") -``` - -### Problem 11 (Trigonometric function) - -```julia -benchmark_problem!("Trigonometric function") -``` - -### Problem 12 (Variably dimensioned function) - -```julia -benchmark_problem!("Variably dimensioned function") -``` - -### Problem 13 (Broyden tridiagonal function) - -```julia -benchmark_problem!("Broyden tridiagonal function") -``` - -### Problem 14 (Broyden banded function) - -```julia -benchmark_problem!("Broyden banded function") -``` - -### Problem 15 (Hammarling 2 by 2 matrix square root problem) - -```julia -benchmark_problem!("Hammarling 2 by 2 matrix square root problem") -``` - -### Problem 16 (Hammarling 3 by 3 matrix square root problem) - -```julia -benchmark_problem!("Hammarling 3 by 3 matrix square root problem") -``` - -### Problem 17 (Dennis and Schnabel 2 by 2 example) - -```julia -benchmark_problem!("Dennis and Schnabel 2 by 2 example") -``` - -### Problem 18 (Sample problem 18) - -```julia -benchmark_problem!("Sample problem 18") -``` - -### Problem 19 (Sample problem 19) - -```julia -benchmark_problem!("Sample problem 19") -``` - -### Problem 20 (Scalar problem f(x) = x(x - 5)^2) - -```julia -benchmark_problem!("Scalar problem f(x) = x(x - 5)^2") -``` - -### Problem 21 (Freudenstein-Roth function) - -```julia -benchmark_problem!("Freudenstein-Roth function") -``` - -### Problem 22 (Boggs function) - -```julia -benchmark_problem!("Boggs function") -``` - -### Problem 23 (Chandrasekhar function) - -```julia -benchmark_problem!("Chandrasekhar function") -``` - -# Summary of successful solvers - -# Finally, we print a summary of which solvers successfully solved which problems. - -```julia -solver_successes = [(solver.name in getfield.(prob[2], :name)) ? "O" : "X" for prob in solver_tracker, solver in solvers_all] -total_successes = [sum(solver_successes[:,i] .== "O") for i in 1:length(solvers_all)] -solver_outcomes = vcat(total_successes', solver_successes); -``` - -```julia; wrap = false; results = "md2html" -using PrettyTables -io = IOBuffer() -println(io, "```@raw html") -pretty_table(io, solver_outcomes; backend = Val(:html), header = getfield.(solvers_all, :name), alignment=:c) -println(io, "```") -Base.Text(String(take!(io))) -``` - -# Summary of General Solver Performance on All Problems - -```julia -fig = begin - LINESTYLES = Dict(:nonlinearsolve => :solid, :simplenonlinearsolve => :dash, - :wrapper => :dot) - ASPECT_RATIO = 1 - WIDTH = 1800 - HEIGHT = round(Int, WIDTH * ASPECT_RATIO) - STROKEWIDTH = 2.5 - - colors = cgrad(:seaborn_bright, length(solvers_all); categorical = true) - cycle = Cycle([:marker], covary = true) - plot_theme = Theme(Lines = (; cycle), Scatter = (; cycle)) - - with_theme(plot_theme) do - fig = Figure(; size = (WIDTH, HEIGHT)) - - axs = Matrix{Any}(undef, 5, 5) - - ls = [] - scs = [] - labels = [] - solver_times = [] - - for i in 1:5, j in 1:5 - idx = 5 * (i - 1) + j - - idx > length(wp_general_tracker) && break - - prob_name, wp = wp_general_tracker[idx] - - ax = Axis(fig[i, j], - xscale = __log10_zero, yscale = __log10_zero, - xtickwidth = STROKEWIDTH, - ytickwidth = STROKEWIDTH, spinewidth = STROKEWIDTH, - title = prob_name, titlegap = 10, - xticklabelsize = 16, yticklabelsize = 16) - - xmin, xmax, ymin, ymax = Inf, -Inf, Inf, -Inf - - for wpᵢ in wp.wps - idx = findfirst(s -> s.name == wpᵢ.name, solvers_all) - errs = getindex.(wpᵢ.errors, :l∞) - times = wpᵢ.times - - emin, emax = extrema(errs) - tmin, tmax = extrema(times) - emin < xmin && (xmin = emin) - emax > xmax && (xmax = emax) - tmin < ymin && (ymin = tmin) - tmax > ymax && (ymax = tmax) - - l = lines!(ax, errs, times; color = colors[idx], linewidth = 5, - linestyle = LINESTYLES[solvers_all[idx].pkg], alpha = 0.8, - label = wpᵢ.name) - sc = scatter!(ax, errs, times; color = colors[idx], markersize = 16, - strokewidth = 2, marker = Cycled(idx), alpha = 0.8, label = wpᵢ.name) - - if wpᵢ.name ∉ labels - push!(ls, l) - push!(scs, sc) - push!(labels, wpᵢ.name) - end - - if wpᵢ.name ∈ first.(solver_times) - idxi = findfirst(x -> first(x) == wpᵢ.name, solver_times) - push!(solver_times[idxi][2], median(times) / length(wp.prob.u0)) - else - push!(solver_times, wpᵢ.name => [median(times) / length(wp.prob.u0)]) - end - end - - xmin = max(xmin, 10^-100) - - xticks = __filter_nearzero(Makie.get_ticks(LogTicks(WilkinsonTicks(5; k_min = 3)), - __log10_zero, tickformatter, xmin, xmax)) - yticks = __filter_nearzero(Makie.get_ticks(LogTicks(WilkinsonTicks(5; k_min = 3)), - __log10_zero, tickformatter, ymin, ymax)) - - ax.xticks = xticks - ax.yticks = yticks - end - - ordering = sortperm(median.(last.(solver_times))) - - fig[0, :] = Label(fig, "Work-Precision Diagram for 23 Test Problems", - fontsize = 24, tellwidth = false, font = :bold) - - fig[:, 0] = Label(fig, "Time (s)", fontsize = 20, tellheight = false, font = :bold, - rotation = π / 2) - fig[end + 1, :] = Label(fig, - L"Error: $\mathbf{||f(u^\ast)||_\infty}$", - fontsize = 20, tellwidth = false, font = :bold) - - Legend(fig[5, 4:5], [[l, sc] for (l, sc) in zip(ls[ordering], scs[ordering])], - labels[ordering], "Successful Solvers"; - framevisible=true, framewidth = STROKEWIDTH, orientation = :horizontal, - titlesize = 20, nbanks = 9, labelsize = 20, halign = :center, - tellheight = false, tellwidth = false, patchsize = (40.0f0, 20.0f0)) - - return fig - end -end -``` - -```julia -save("summary_wp_23test_problems.svg", fig) -``` - -```julia, echo = false -using SciMLBenchmarks -SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file]) -``` - -# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl