Add Parallel Particle Swarm Optimizer Benchmarks

benchmarks/GlobalOptimization/CondaPkg.toml

@@ -0,0 +1,2 @@
+[deps]
+python = ">=3.10,<3.13"

benchmarks/GlobalOptimization/pso_global_optimizers.jmd

@@ -0,0 +1,377 @@
+---
+title: PSO Global Optimizer Benchmarks
+author: Utkarsh, Chris Rackauckas
+---
+
+This benchmark evaluates Particle Swarm Optimization (PSO) variants from
+[ParallelParticleSwarms.jl](https://github.com/SciML/ParallelParticleSwarms.jl)
+against established global optimizers on the
+[BlackBoxOptimizationBenchmarking.jl](https://github.com/jonathanBieler/BlackBoxOptimizationBenchmarking.jl)
+suite, using the [Optimization.jl](https://github.com/SciML/Optimization.jl) interface.
+This tests both iterations and wall-clock time vs accuracy, i.e. for a given budget
+(in iterations or time), what percentage of problems from the set is a solver able to solve.
+
+## Setup
+```julia
+using Random
+Random.seed!(42)
+
+using BlackBoxOptimizationBenchmarking, Plots, Optimization, Memoize, Statistics
+import BlackBoxOptimizationBenchmarking: Chain, BenchmarkSetup, BenchmarkResults,
+    BBOBFunction, FunctionCallsCounter, solve_problem, pinit, compute_CI
+const BBOB = BlackBoxOptimizationBenchmarking
+
+ENV["GKSwstype"] = "nul"
+gr()
+
+using OptimizationBBO, OptimizationOptimJL, OptimizationNLopt, OptimizationSciPy
+
+using ParallelParticleSwarms
+using KernelAbstractions: CPU
+using StaticArrays, LinearAlgebra
+
+const PSOKernel     = ParallelParticleSwarms.ParallelPSOKernel
+const SyncPSOKernel = ParallelParticleSwarms.ParallelSyncPSOKernel
+const SerPSO        = ParallelParticleSwarms.SerialPSO
+const HPso          = ParallelParticleSwarms.HybridPSO
+
+const BACKEND = CPU()
+```
+```julia
+function make_success_tracker(f_raw, f_opt, Δf)
+    t0 = Ref(time())
+    time_to_success = Ref(Inf)
+    function tracked_f(u)
+        val = f_raw(u)
+        if val < Δf + f_opt && time_to_success[] == Inf
+            time_to_success[] = time() - t0[]
+        end
+        return val
+    end
+    return tracked_f, t0, time_to_success
+end
+
+function solve_problem_timed(optimizer::BenchmarkSetup, tracked_f, D::Int, run_length::Int;
+        u0 = pinit(D))
+    method = optimizer.method
+    optf = OptimizationFunction((u, _) -> tracked_f(u), AutoForwardDiff())
+    if optimizer.isboxed
+        prob = OptimizationProblem(optf, u0, lb = fill(-5.5, D), ub = fill(5.5, D))
+    else
+        prob = OptimizationProblem(optf, u0)
+    end
+    solve(prob, method; maxiters = run_length)
+end
+
+function solve_problem_timed(m::Chain, tracked_f, D::Int, run_length::Int)
+    rl1 = round(Int, m.p * run_length)
+    rl2 = run_length - rl1
+    sol = solve_problem_timed(m.first, tracked_f, D, rl1)
+    solve_problem_timed(m.second, tracked_f, D, rl2; u0 = sol.u)
+end
+
+function benchmark_time_to_success(
+    optimizer::Union{Chain, BenchmarkSetup}, funcs::Vector{BBOBFunction};
+    Ntrials::Int = 20, dimension::Int = 3, Δf::Real = 1e-6, max_run_length::Int = 100_000
+)
+    all_times = Float64[]
+    for f in funcs
+        for _ in 1:Ntrials
+            tracked_f, t0_ref, tts_ref = make_success_tracker(f.f, f.f_opt, Δf)
+            try
+                t0_ref[] = time()
+                solve_problem_timed(optimizer, tracked_f, dimension, max_run_length)
+                push!(all_times, tts_ref[])
+            catch err
+                push!(all_times, Inf)
+                @warn(string(optimizer, " failed: ", err))
+            end
+        end
+    end
+    return all_times
+end
+
+benchmark_time_to_success(optimizer, funcs::Vector{BBOBFunction}; kwargs...) =
+    benchmark_time_to_success(BenchmarkSetup(optimizer), funcs; kwargs...)
+
+function success_rate_cdf(all_times::Vector{Float64}, time_thresholds::AbstractVector{Float64})
+    N = length(all_times)
+    return [count(x -> x <= t, all_times) / N for t in time_thresholds]
+end
+```
+
+PSO variants require an `SVector` interface and a custom solve loop.
+```julia
+function pso_solve(opt, f::BBOBFunction, D::Int, maxiters::Int; local_maxiters::Int = 20)
+    obj  = (x, p) -> f.f(Vector(x))
+    ad   = opt isa HPso ? AutoForwardDiff() : Optimization.SciMLBase.NoAD()
+    optf = OptimizationFunction{false}(obj, ad)
+    lb   = SVector{D}(ntuple(_ -> -5.0, Val(D)))
+    ub   = SVector{D}(ntuple(_ ->  5.0, Val(D)))
+    x0   = SVector{D}(ntuple(_ -> -5.0 + rand() * 10.0, Val(D)))
+    prob = OptimizationProblem{false}(optf, x0, nothing; lb, ub)
+    if opt isa HPso
+        solve(prob, opt; maxiters, local_maxiters, abstol = 1e-8, reltol = 1e-8)
+    else
+        solve(prob, opt; maxiters)
+    end
+end
+
+function _extract_u(sol, D)
+    u = sol.u
+    u isa SVector && return u
+    u[]
+end
+
+function pso_benchmark(opt, funcs, run_length;
+        Ntrials = 20, dimension = 3, local_maxiters = 20, Δf = 1e-6, CI_quantile = 0.25)
+    Nf = length(funcs)
+    Nr = length(run_length)
+    success = zeros(Float64, Nf, Nr)
+    dist    = zeros(Float64, Nf, Nr)
+    fmin    = zeros(Float64, Nf, Nr)
+    # warmup
+    for f in funcs
+        try pso_solve(opt, f, dimension, 3; local_maxiters) catch; end
+    end
+    t0 = time()
+    for (fi, f) in enumerate(funcs)
+        xopt = SVector{dimension}(f.x_opt[1:dimension])
+        for (ri, rl) in enumerate(run_length)
+            hits = 0; dsum = 0.0; fsum = 0.0
+            for _ in 1:Ntrials
+                sol = try pso_solve(opt, f, dimension, rl; local_maxiters) catch; nothing end
+                if sol !== nothing
+                    u    = _extract_u(sol, dimension)
+                    fval = sol.objective isa Real ? Float64(sol.objective) : Float64(sol.objective[])
+                    hits += abs(fval - f.f_opt) < Δf ? 1 : 0
+                    dsum += norm(u .- xopt)
+                    fsum += fval - f.f_opt
+                end
+            end
+            success[fi, ri] = hits / Ntrials
+            dist[fi, ri]    = dsum / Ntrials
+            fmin[fi, ri]    = fsum / Ntrials
+        end
+    end
+    elapsed = time() - t0
+    Neff = Ntrials * Nf
+    sr   = vec(mean(success, dims = 1))
+    sc   = vec(sum(success .* Ntrials, dims = 1)) .|> round .|> Int
+    ci   = BBOB.compute_CI(sr, Neff, CI_quantile)
+    BenchmarkResults(
+        run_length                = collect(run_length),
+        success_count             = sc,
+        success_rate              = sr,
+        success_rate_qlow         = ci.success_rate_qlow,
+        success_rate_qhigh        = ci.success_rate_qhigh,
+        distance_to_minimizer     = vec(mean(dist, dims = 1)),
+        minimum                   = vec(mean(fmin, dims = 1)),
+        runtime                   = elapsed,
+        Neffective                = Neff,
+        callcount                 = Float64.(run_length),
+        success_rate_per_function = [success[fi, end] for fi in 1:Nf],
+    )
+end
+
+function pso_tts(opt, funcs; Ntrials = 20, dimension = 3, Δf = 1e-6,
+        local_maxiters = 20, max_run_length = 100_000)
+    all_times = Float64[]
+    D  = dimension
+    lb = SVector{D}(ntuple(_ -> -5.0, Val(D)))
+    ub = SVector{D}(ntuple(_ ->  5.0, Val(D)))
+    for f in funcs
+        for _ in 1:Ntrials
+            tracked_f, t0_ref, tts_ref = make_success_tracker(f.f, f.f_opt, Δf)
+            obj  = (x, p) -> tracked_f(Vector(x))
+            ad   = opt isa HPso ? AutoForwardDiff() : Optimization.SciMLBase.NoAD()
+            optf = OptimizationFunction{false}(obj, ad)
+            x0   = SVector{D}(ntuple(_ -> -5.0 + rand() * 10.0, Val(D)))
+            prob = OptimizationProblem{false}(optf, x0, nothing; lb, ub)
+            try
+                t0_ref[] = time()
+                if opt isa HPso
+                    solve(prob, opt; maxiters = max_run_length, local_maxiters,
+                        abstol = 1e-8, reltol = 1e-8)
+                else
+                    solve(prob, opt; maxiters = max_run_length)
+                end
+                push!(all_times, tts_ref[])
+            catch err
+                push!(all_times, Inf)
+                @warn(string(opt, " failed: ", err))
+            end
+        end
+    end
+    all_times
+end
+```
+```julia
+chain = (t; isboxed = false) -> Chain(
+    BenchmarkSetup(t, isboxed = isboxed),
+    BenchmarkSetup(NelderMead(), isboxed = false),
+    0.9)
+
+test_functions = BBOB.list_functions()
+dimension      = 3
+run_length     = round.(Int, 10 .^ LinRange(1, 5, 30))
+num_particles  = 1024
+
+PSO_KEYS = Set(["SerialPSO", "PSOKernel", "SyncPSOKernel", "HybridPSO_LBFGS"])
+
+setup = Dict(
+    # Baseline global optimizers 
+    "NelderMead"                               => NelderMead(),
+    "NLopt_GN_CRS2_LM"                        => chain(NLopt.GN_CRS2_LM(), isboxed = true),
+    "BBO_adaptive_de_rand_1_bin"               => chain(BBO_adaptive_de_rand_1_bin(), isboxed = true),
+    "BBO_adaptive_de_rand_1_bin_radiuslimited" => chain(BBO_adaptive_de_rand_1_bin_radiuslimited(), isboxed = true),
+    "BBO_de_rand_2_bin"                        => chain(BBO_de_rand_2_bin(), isboxed = true),
+    "ScipyDifferentialEvolution"               => chain(ScipyDifferentialEvolution(), isboxed = true),
+    # PSO variants 
+    "SerialPSO"        => SerPSO(num_particles),
+    "PSOKernel"        => PSOKernel(num_particles; backend = BACKEND, global_update = true),
+    "SyncPSOKernel"    => SyncPSOKernel(num_particles; backend = BACKEND),
+    "HybridPSO_LBFGS"  => HPso(pso = PSOKernel(num_particles; backend = BACKEND); backend = BACKEND),
+)
+
+@memoize run_bench(algo) = algo in PSO_KEYS ?
+    pso_benchmark(setup[algo], test_functions, run_length; Ntrials = 20, dimension) :
+    BBOB.benchmark(setup[algo], test_functions, run_length; Ntrials = 40, dimension)
+
+@memoize run_tts(algo) = algo in PSO_KEYS ?
+    pso_tts(setup[algo], test_functions; Ntrials = 20, dimension) :
+    benchmark_time_to_success(setup[algo], test_functions; Ntrials = 40, dimension)
+```
+
+## Test all (iterations)
+```julia
+const MARKERS = [:circle, :diamond, :utriangle, :square, :star5, :dtriangle, :pentagon,
+    :hexagon, :cross, :xcross, :rtriangle, :ltriangle, :star4, :star8, :heptagon, :octagon,
+    :vline, :hline, :+, :x]
+const LINESTYLES = [:solid, :dash, :dot, :dashdot, :dashdotdot]
+
+labels  = collect(keys(setup))
+results = Array{BBOB.BenchmarkResults}(undef, length(setup))
+
+# PSO variants first
+for (i, algo) in enumerate(labels)
+    algo in PSO_KEYS || continue
+    @info "PSO: $algo ..."
+    try
+        results[i] = run_bench(algo)
+        @info "  done" success = round(results[i].success_rate[end], digits = 3)
+    catch err
+        @warn "$algo failed" exception = (err, catch_backtrace())
+    end
+end
+
+# Baseline optimizers (threaded)
+Threads.@threads for (i, algo) in collect(enumerate(labels))
+    algo in PSO_KEYS && continue
+    results[i] = run_bench(algo)
+end
+
+results
+```
+
+## Success Rate vs. Iterations
+```julia
+idx = sortperm([b.success_rate[end] for b in results], rev = true)
+
+p = plot(xscale = :log10, legend = :outerright,
+    size = (700, 350), margin = 10Plots.px, dpi = 200)
+for (j, i) in enumerate(idx)
+    plot!(results[i], label = labels[i], showribbon = false,
+        lw = 2.5, xlim = (1, 1e5), x = :run_length,
+        markershape = MARKERS[mod1(j, length(MARKERS))],
+        linestyle = LINESTYLES[mod1(j, length(LINESTYLES))],
+        markersize = 4, markerstrokewidth = 0)
+end
+p
+```
+
+## Test all (wall-clock time to success)
+
+For the time-based benchmark, each optimizer is run once with a large iteration budget
+(100,000 iterations) per (function, trial) pair. The objective function is wrapped to
+detect the first evaluation that achieves the success criterion (objective < Δf + f_opt)
+and record the wall-clock time at that moment. This gives a true "time to success" for
+each trial, from which we build a CDF.
+```julia
+tts_results = Dict{String, Vector{Float64}}()
+
+# PSO variants first
+for algo in labels
+    algo in PSO_KEYS || continue
+    tts_results[algo] = run_tts(algo)
+end
+
+# Baseline optimizers
+for algo in labels
+    algo in PSO_KEYS && continue
+    tts_results[algo] = run_tts(algo)
+end
+```
+
+## Success Rate vs. Wall-Clock Time
+```julia
+all_finite      = filter(isfinite, vcat(values(tts_results)...))
+t_lo            = minimum(all_finite) / 2
+t_hi            = maximum(all_finite) * 2
+time_thresholds = 10 .^ range(log10(t_lo), log10(t_hi), length = 50)
+
+cdfs = Dict(algo => success_rate_cdf(tts_results[algo], time_thresholds) for algo in labels)
+idx  = sortperm([cdfs[l][end] for l in labels], rev = true)
+
+p = plot(xscale = :log10, legend = :outerright,
+    size = (700, 350), margin = 10Plots.px, dpi = 200,
+    xlabel = "Wall time (s)", ylabel = "Success rate", ylim = (0, 1))
+for (j, i) in enumerate(idx)
+    plot!(time_thresholds, cdfs[labels[i]], label = labels[i], lw = 2.5,
+        markershape = MARKERS[mod1(j, length(MARKERS))],
+        linestyle = LINESTYLES[mod1(j, length(LINESTYLES))],
+        markersize = 4, markerstrokewidth = 0)
+end
+p
+```
+
+## Success Rate per Function Heatmap
+```julia
+success_rate_per_function = reduce(hcat, b.success_rate_per_function for b in results)
+
+idx     = sortperm(mean(success_rate_per_function, dims = 1)[:], rev = false)
+idxfunc = 1:length(test_functions)
+
+p = heatmap(
+    string.(test_functions)[idxfunc], labels[idx], success_rate_per_function[idxfunc, idx]',
+    cmap = :RdYlGn, xticks = :all, yticks = :all, xrotation = 45, dpi = 200)
+```
+
+## Distance to Minimizer vs. Iterations
+```julia
+idx = sortperm([b.distance_to_minimizer[end] for b in results], rev = false)
+
+p = plot(xscale = :log10, legend = :outerright,
+    size = (900, 500), margin = 10Plots.px, ylim = (0, 5))
+for (j, i) in enumerate(idx)
+    plot!(results[i].run_length, results[i].distance_to_minimizer, label = labels[i],
+        showribbon = false, lw = 2, xlim = (1, 1e5),
+        xlabel = "Iterations", ylabel = "Mean distance to minimum",
+        markershape = MARKERS[mod1(j, length(MARKERS))],
+        linestyle = LINESTYLES[mod1(j, length(LINESTYLES))],
+        markersize = 4, markerstrokewidth = 0)
+end
+p
+```
+
+## Relative Runtime
+```julia
+ref      = findfirst("NelderMead" .== labels)
+runtimes = getfield.(results, :runtime)
+runtimes = runtimes ./ runtimes[ref]
+
+bar(
+    labels, runtimes, xrotation = 45, xticks = :all, ylabel = "Run time relative to NM",
+    yscale = :log10, yticks = [0.1, 1, 10, 100],
+    legend = false, margin = 25Plots.px)
+```