Merge branch 'master' into multistartoptimization

Author: ChrisRackauckas

docs/make.jl

@@ -31,6 +31,7 @@ makedocs(
             "CMAEvolutionStrategy.jl" => "optimization_packages/cmaevolutionstrategy.md",
             "Evolutionary.jl" => "optimization_packages/evolutionary.md",
             "Flux.jl" => "optimization_packages/flux.md",
+            "GCMAES.jl" => "optimization_packages/gcmaes.md",
             "MathOptInterface.jl" => "optimization_packages/mathoptinterface.md",
             "MultistartOptimization.jl" => "optimization_packages/multistartoptimization.md",
             "Metaheuristics.jl" => "optimization_packages/metaheuristics.md",

docs/src/optimization_packages/gcmaes.md

@@ -0,0 +1,32 @@
+# GCMAES.jl
+[`GCMAES`](https://github.com/AStupidBear/GCMAES.jl) is a Julia package implementing the **Gradient-based Covariance Matrix Adaptation Evolutionary Strategy** which can utilize the gradient information to speed up the optimization process.
+
+The GCMAES algorithm is called by `GCMAESOpt()` and the initial search variance is set as a keyword argument `σ0` (default: `σ0 = 0.2`)
+
+## Global Optimizer
+### Without Constraint Equations
+
+The method in [`GCMAES`](https://github.com/AStupidBear/GCMAES.jl) is performing global optimization on problems without
+constraint equations. However, lower and upper constraints set by `lb` and `ub` in the `OptimizationProblem` are required.
+
+The Rosenbrock function can optimized using the `GCMAESOpt()` without utilizing the gradient information as follows:
+
+```julia
+rosenbrock(x, p) =  (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
+x0 = zeros(2)
+p  = [1.0, 100.0]
+f = OptimizationFunction(rosenbrock)
+prob = GalacticOptim.OptimizationProblem(f, x0, p, lb = [-1.0,-1.0], ub = [1.0,1.0])
+sol = solve(prob, GCMAESOpt())
+```
+
+We can also utilise the gradient information of the optimization problem to aid the optimization as follows:
+
+```julia
+rosenbrock(x, p) =  (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
+x0 = zeros(2)
+p  = [1.0, 100.0]
+f = OptimizationFunction(rosenbrock, GalacticOptim.ForwardDiff)
+prob = GalacticOptim.OptimizationProblem(f, x0, p, lb = [-1.0,-1.0], ub = [1.0,1.0])
+sol = solve(prob, GCMAESOpt())
+```

docs/src/optimization_packages/speedmapping.md

@@ -0,0 +1,26 @@
+# SpeedMapping.jl
+[`SpeedMapping`](https://github.com/NicolasL-S/SpeedMapping.jl) accelerates the convergence of a mapping to a fixed point by the Alternating cyclic extrapolation algorithm which can also perform multivariate optimization based on the gradient function.
+
+The SpeedMapping algorithm is called by `SpeedMappingOpt()`
+
+## Global Optimizer
+### Without Constraint Equations
+
+The method in [`SpeedMapping`](https://github.com/NicolasL-S/SpeedMapping.jl) is performing optimization on problems without
+constraint equations. Lower and upper constraints set by `lb` and `ub` in the `OptimizationProblem` are optional.
+
+If no AD backend is defined via `OptimizationFunction` the gradient is calculated via `SpeedMapping`'s ForwardDiff AD backend.
+
+The Rosenbrock function can be optimized using the `SpeedMappingOpt()` with and without bound as follows:
+
+```julia
+rosenbrock(x, p) =  (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
+x0 = zeros(2)
+p  = [1.0, 100.0]
+f = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
+prob = OptimizationProblem(f, x0, _p)
+sol = solve(prob,SpeedMappingOpt())
+
+prob = OptimizationProblem(f, x0, _p;lb=[0.0,0.0], ub=[1.0,1.0])
+sol = solve(prob,SpeedMappingOpt())
+```

src/GalacticOptim.jl

@@ -24,6 +24,7 @@ function __init__()
     @require CMAEvolutionStrategy="8d3b24bd-414e-49e0-94fb-163cc3a3e411" include("solve/cmaevolutionstrategy.jl")
     @require Evolutionary="86b6b26d-c046-49b6-aa0b-5f0f74682bd6" include("solve/evolutionary.jl")
     @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" include("solve/flux.jl")
+    @require GCMAES="4aa9d100-eb0f-11e8-15f1-25748831eb3b" include("solve/gcmaes.jl")
     @require MathOptInterface="b8f27783-ece8-5eb3-8dc8-9495eed66fee" include("solve/moi.jl")
     @require MultistartOptimization="3933049c-43be-478e-a8bb-6e0f7fd53575" include("solve/multistartoptimization.jl")
     @require NLopt="76087f3c-5699-56af-9a33-bf431cd00edd" begin
@@ -45,6 +46,7 @@ function __init__()
         end
     end
     @require Metaheuristics="bcdb8e00-2c21-11e9-3065-2b553b22f898" include("solve/metaheuristics.jl")
+    @require SpeedMapping="f1835b91-879b-4a3f-a438-e4baacf14412" include("solve/speedmapping.jl")
 
     # AD backends
     @require FiniteDiff="6a86dc24-6348-571c-b903-95158fe2bd41" include("function/finitediff.jl")

src/solve/gcmaes.jl

@@ -0,0 +1,75 @@
+export GCMAESOpt
+
+struct GCMAESOpt end
+
+function __map_optimizer_args(prob::OptimizationProblem, opt::GCMAESOpt; 
+    cb=nothing, 
+    maxiters::Union{Number, Nothing}=nothing, 
+    maxtime::Union{Number, Nothing}=nothing, 
+    abstol::Union{Number, Nothing}=nothing, 
+    reltol::Union{Number, Nothing}=nothing,
+    kwargs...)
+
+    # add optimiser options from kwargs
+    mapped_args = (; kwargs...)
+
+    if !(isnothing(maxiters))
+        mapped_args = (; mapped_args..., maxiter = maxiters)
+    end
+
+    if !(isnothing(maxtime))
+        @warn "common maxtime is currently not used by $(opt)"
+    end
+
+    if !isnothing(abstol)
+        @warn "common abstol is currently not used by $(opt)"
+    end
+
+    if !isnothing(reltol)
+        @warn "common reltol is currently not used by $(opt)"
+    end
+
+    return mapped_args
+end
+
+
+function __solve(prob::OptimizationProblem, opt::GCMAESOpt;
+                 maxiters::Union{Number, Nothing} = nothing,
+                 maxtime::Union{Number, Nothing} = nothing,
+                 abstol::Union{Number, Nothing}=nothing,
+                 reltol::Union{Number, Nothing}=nothing,
+                 progress = false,
+                 σ0=0.2,
+                 kwargs...)
+    local x
+    local G = similar(prob.u0)
+
+    maxiters = _check_and_convert_maxiters(maxiters)
+    maxtime = _check_and_convert_maxtime(maxtime)
+
+    f = instantiate_function(prob.f,prob.u0,prob.f.adtype,prob.p)
+
+    _loss = function(θ)
+        x = f.f(θ, prob.p)
+        return x[1]
+    end
+
+    if !isnothing(f.grad)
+        g = function (θ)
+            f.grad(G, θ)
+            return G
+        end
+    end
+    
+    opt_args = _map_optimizer_args(prob,opt, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)
+    
+    t0 = time()
+    if prob.sense === MaxSense 
+        opt_xmin, opt_fmin, opt_ret = GCMAES.maximize(isnothing(f.grad) ? _loss : (_loss,g), prob.u0, σ0, prob.lb, prob.ub; opt_args...)
+    else
+        opt_xmin, opt_fmin, opt_ret = GCMAES.minimize(isnothing(f.grad) ? _loss : (_loss,g), prob.u0, σ0, prob.lb, prob.ub; opt_args...)
+    end
+    t1 = time()
+
+    SciMLBase.build_solution(prob, opt, opt_xmin, opt_fmin; retcode=Symbol(Bool(opt_ret)))
+end

src/solve/speedmapping.jl

@@ -0,0 +1,68 @@
+export SpeedMappingOpt
+
+struct SpeedMappingOpt end
+
+function __map_optimizer_args(prob::OptimizationProblem, opt::SpeedMappingOpt; 
+    cb=nothing, 
+    maxiters::Union{Number, Nothing}=nothing, 
+    maxtime::Union{Number, Nothing}=nothing, 
+    abstol::Union{Number, Nothing}=nothing, 
+    reltol::Union{Number, Nothing}=nothing,
+    kwargs...)
+
+    # add optimiser options from kwargs
+    mapped_args = (; kwargs...)
+
+    if !(isnothing(maxiters))
+        @info "maxiters defines maximum gradient calls for $(opt)"
+        mapped_args = (; mapped_args..., maps_limit = maxiters)
+    end
+
+    if !(isnothing(maxtime))
+        mapped_args = (; mapped_args..., time_limit = maxtime)
+    end
+
+    if !isnothing(abstol)
+        @warn "common abstol is currently not used by $(opt)"
+    end
+
+    if !isnothing(reltol)
+        @warn "common reltol is currently not used by $(opt)"
+    end
+
+    return mapped_args
+end
+
+
+function __solve(prob::OptimizationProblem, opt::SpeedMappingOpt;
+                 maxiters::Union{Number, Nothing} = nothing,
+                 maxtime::Union{Number, Nothing} = nothing,
+                 abstol::Union{Number, Nothing}=nothing,
+                 reltol::Union{Number, Nothing}=nothing,
+                 progress = false,
+                 kwargs...)
+    local x
+
+    maxiters = _check_and_convert_maxiters(maxiters)
+    maxtime = _check_and_convert_maxtime(maxtime)
+
+    f = instantiate_function(prob.f,prob.u0,prob.f.adtype,prob.p)
+
+    _loss = function(θ)
+        x = f.f(θ, prob.p)
+        return first(x)
+    end
+    
+    if isnothing(f.grad)
+        @info "SpeedMapping's ForwardDiff AD backend is used to calculate the gradient information."
+    end 
+
+    opt_args = _map_optimizer_args(prob,opt, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)
+    
+    t0 = time()
+    opt_res = SpeedMapping.speedmapping(prob.u0; f = _loss, g! =f.grad, lower=prob.lb, upper=prob.ub, opt_args...)
+    t1 = time()
+    opt_ret = Symbol(opt_res.converged)
+
+    SciMLBase.build_solution(prob, opt, opt_res.minimizer, _loss(opt_res.minimizer); original=opt_res, retcode=opt_ret)
+end

test/Project.toml

@@ -7,6 +7,7 @@ Evolutionary = "86b6b26d-c046-49b6-aa0b-5f0f74682bd6"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GCMAES = "4aa9d100-eb0f-11e8-15f1-25748831eb3b"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
@@ -30,6 +31,7 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
@@ -42,6 +44,7 @@ Evolutionary = ">= 0.10.0"
 FiniteDiff = ">= 2.8.1"
 Flux = ">= 0.12.6"
 ForwardDiff = ">= 0.10.19"
+GCMAES = ">= 0.1.25"
 Ipopt = ">= 0.7.0"
 IterTools = ">= 1.3.0"
 MathOptInterface = ">= 0.9.22"
@@ -63,6 +66,7 @@ Optim = ">= 1.4.1"
 OrdinaryDiffEq = ">= 5"
 ReverseDiff = ">= 1.9.0"
 SafeTestsets = ">= 0.0.1"
+SpeedMapping = ">= 0.3.0"
 Tracker = ">= 0.2"
 Zygote = ">= 0.5"
-julia = "1.5"
+julia = "1.5"

test/rosenbrock.jl

@@ -169,7 +169,7 @@ end
 # sol = solve(prob, QuadDirect(); splits = ([-0.5, 0.0, 0.5],[-0.5, 0.0, 0.5]))
 # @test 10*sol.minimum < l1
 
-@testset "Evolutionary, BlackBoxOptim, Metaheuristics, Nonconvex, MultistartOptimization" begin
+@testset "Evolutionary, BlackBoxOptim, Metaheuristics, Nonconvex, GCMAES, SpeedMapping, MultistartOptimization" begin
     optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())
     using Evolutionary
     prob = GalacticOptim.OptimizationProblem(optprob, x0, _p)
@@ -423,6 +423,33 @@ end
     sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(;maxeval=100))
     @test 10*sol.minimum < l1
 
+    using GCMAES
+    f_ad = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
+    f_noad = OptimizationFunction(rosenbrock)
+
+    prob = GalacticOptim.OptimizationProblem(f_ad, x0, _p, lb=[-1.0, -1.0], ub=[1.0, 1.0])
+    sol = solve(prob, GCMAESOpt(), maxiters=1000)
+    @test 10*sol.minimum < l1
+
+    prob = GalacticOptim.OptimizationProblem(f_noad, x0, _p, lb=[-1.0, -1.0], ub=[1.0, 1.0])
+    sol = solve(prob, GCMAESOpt(), maxiters=1000)
+    @test 10*sol.minimum < l1
+
+    using SpeedMapping
+    f = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
+    prob = OptimizationProblem(f, x0, _p)
+    sol = solve(prob,SpeedMappingOpt())
+
+    prob = OptimizationProblem(f, x0, _p;lb=[0.0,0.0], ub=[1.0,1.0])
+    sol = solve(prob,SpeedMappingOpt())
+
+    f = OptimizationFunction(rosenbrock)
+    prob = OptimizationProblem(f, x0, _p)
+    sol = solve(prob,SpeedMappingOpt())
+
+    prob = OptimizationProblem(f, x0, _p;lb=[0.0,0.0], ub=[1.0,1.0])
+    sol = solve(prob,SpeedMappingOpt())
+  
     using MultistartOptimization
     f = OptimizationFunction(rosenbrock,GalacticOptim.AutoForwardDiff())
     prob = GalacticOptim.OptimizationProblem(f, x0, _p, lb = [-1.0,-1.0], ub = [1.5,1.5])