diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ad997e49e..2a432487a 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -32,6 +32,7 @@ jobs: - OptimizationNOMAD - OptimizationOptimJL - OptimizationOptimisers + - OptimizationPRIMA - OptimizationQuadDIRECT - OptimizationSpeedMapping - OptimizationPolyalgorithms diff --git a/.gitignore b/.gitignore index bd2b72711..1f9682187 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,5 @@ .DS_Store -/Manifest.toml +Manifest.toml /dev/ /docs/build/ .vscode \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index e2a5b27e3..41734e54a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -62,8 +62,8 @@ packages. | Metaheuristics | ❌ | ❌ | ❌ | ✅ | ❌ | ✅ | 🟡 | | NOMAD | ❌ | ❌ | ❌ | ✅ | ❌ | ✅ | 🟡 | | NLopt | ✅ | ❌ | ✅ | ✅ | 🟡 | ✅ | 🟡 | -| Nonconvex | ✅ | ✅ | ✅ | ✅ | 🟡 | ✅ | 🟡 | | Optim | ✅ | ✅ | ✅ | ✅ | ✅ | ✅ | ✅ | +| PRIMA | ❌ | ❌ | ✅ | ✅ | ✅ | ❌ | ❌ | | QuadDIRECT | ❌ | ❌ | ❌ | ✅ | ❌ | ✅ | ❌ | ✅ = supported diff --git a/docs/src/optimization_packages/blackboxoptim.md b/docs/src/optimization_packages/blackboxoptim.md index 6d0c57f7f..ca5b2385b 100644 --- a/docs/src/optimization_packages/blackboxoptim.md +++ b/docs/src/optimization_packages/blackboxoptim.md @@ -1,6 +1,6 @@ # BlackBoxOptim.jl -[`BlackBoxOptim`](https://github.com/robertfeldt/BlackBoxOptim.jl) is a Julia package implementing **(Meta-)heuristic/stochastic algorithms** that do not require for the optimized function to be differentiable. +[`BlackBoxOptim`](https://github.com/robertfeldt/BlackBoxOptim.jl) is a Julia package implementing **(Meta-)heuristic/stochastic algorithms** that do not require differentiability. ## Installation: OptimizationBBO.jl diff --git a/docs/src/optimization_packages/nonconvex.md b/docs/src/optimization_packages/nonconvex.md deleted file mode 100644 index f68e12977..000000000 --- a/docs/src/optimization_packages/nonconvex.md +++ /dev/null @@ -1,85 +0,0 @@ -# Nonconvex.jl - -[`Nonconvex`](https://github.com/JuliaNonconvex/Nonconvex.jl) is a Julia package implementing and wrapping nonconvex constrained optimization algorithms. - -## Installation: OptimizationNonconvex.jl - -To use this package, install the OptimizationNonconvex package: - -```julia -import Pkg; -Pkg.add("OptimizationNonconvex"); -``` - -## Global Optimizer - -### Without Constraint Equations - -A `Nonconvex` algorithm is called using one of the following: - - - [Method of moving asymptotes (MMA)](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/mma/#Method-of-moving-asymptotes-(MMA)): - - + `MMA87()` - + `MMA02()` - - - [Ipopt](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/ipopt/): - - + `IpoptAlg()` - - [NLopt](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/nlopt/): - - + `NLoptAlg(solver)` where solver can be any of the `NLopt` algorithms - - [Augmented Lagrangian algorithm](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/auglag/): - - + `AugLag()` - + only works with constraints - - [Mixed integer nonlinear programming (MINLP)](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/minlp/): - - + Juniper + Ipopt: `JuniperIpoptAlg()` - + Pavito + Ipopt + Cbc: `PavitoIpoptCbcAlg()` - - [Multi-start optimization](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/hyperopt/): - - + `HyperoptAlg(subsolver)` where `subalg` can be any of the described `Nonconvex` algorithm - - [Surrogate-assisted Bayesian optimization](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/surrogate/) - - + `BayesOptAlg(subsolver)` where `subalg` can be any of the described `Nonconvex` algorithm - - [Multiple Trajectory Search](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/mts/) - - + `MTSAlg()` - -When performing optimizing a combination of integer and floating-point parameters, the `integer` keyword has to be set. It takes a boolean vector indicating which parameter is an integer. - -## Notes - -Some optimizer may require further options to be defined in order to work. - -The currently available algorithms are listed [here](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/algorithms/) - -The algorithms in [`Nonconvex`](https://julianonconvex.github.io/Nonconvex.jl/stable/algorithms/algorithms/) are performing global optimization on problems without constraint equations. However, lower and upper constraints set by `lb` and `ub` in the `OptimizationProblem` are required. - -## Examples - -The Rosenbrock function can be optimized using the Method of moving asymptotes algorithm `MMA02()` 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 = Optimization.OptimizationProblem(f, x0, p, lb = [-1.0, -1.0], ub = [1.0, 1.0]) -sol = solve(prob, MMA02(), maxiters = 100000, maxtime = 1000.0) -``` - -The options of for a sub-algorithm are passed simply as a NamedTuple and Optimization.jl infers the correct `Nonconvex` options struct: - -```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 = Optimization.OptimizationProblem(f, x0, p, lb = [-1.0, -1.0], ub = [1.0, 1.0]) -sol = solve(prob, HyperoptAlg(IpoptAlg()), sub_options = (; max_iter = 100)) -``` - -### With Constraint Equations - -While `Nonconvex.jl` supports such constraints, `Optimization.jl` currently does not relay these constraints. diff --git a/docs/src/optimization_packages/prima.md b/docs/src/optimization_packages/prima.md new file mode 100644 index 000000000..6a57c933a --- /dev/null +++ b/docs/src/optimization_packages/prima.md @@ -0,0 +1,49 @@ +# PRIMA.jl + +[PRIMA.jl](https://github.com/libprima/PRIMA.jl) is a julia wrapper for the fortran library [prima](https://github.com/libprima/prima) which implements Powell's derivative free optimization methods. + +## Installation: OptimizationPRIMA + +To use this package, install the OptimizationPRIMA package: + +```julia +import Pkg; +Pkg.add("OptimizationPRIMA"); +``` + +## Local Optimizer + +The five Powell's algorithms of the prima library are provided by the PRIMA.jl package: + +`UOBYQA`: (Unconstrained Optimization BY Quadratic Approximations) is for unconstrained optimization, that is Ω = ℝⁿ. + +`NEWUOA`: is also for unconstrained optimization. According to M.J.D. Powell, newuoa is superior to uobyqa. + +`BOBYQA`: (Bounded Optimization BY Quadratic Approximations) is for simple bound constrained problems, that is Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu }. + +`LINCOA`: (LINearly Constrained Optimization) is for constrained optimization problems with bound constraints, linear equality constraints, and linear inequality constraints. + +`COBYLA`: (Constrained Optimization BY Linear Approximations) is for general constrained problems with bound constraints, non-linear constraints, linear equality constraints, and linear inequality constraints. + +```@example PRIMA +rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 +x0 = zeros(2) +_p = [1.0, 100.0] + +prob = OptimizationProblem(rosenbrock, x0, _p) + +sol = Optimization.solve(prob, UOBYQA(), maxiters = 1000) + +sol = Optimization.solve(prob, NEWUOA(), maxiters = 1000) + +sol = Optimization.solve(prob, BOBYQA(), maxiters = 1000) + +sol = Optimization.solve(prob, LINCOA(), maxiters = 1000) + +function con2_c(res, x, p) + res .= [x[1] + x[2], x[2] * sin(x[1]) - x[1]] +end +optprob = OptimizationFunction(rosenbrock, AutoForwardDiff(), cons = con2_c) +prob = OptimizationProblem(optprob, x0, _p, lcons = [1, -100], ucons = [1, 100]) +sol = Optimization.solve(prob, COBYLA(), maxiters = 1000) +``` diff --git a/ext/OptimizationEnzymeExt.jl b/ext/OptimizationEnzymeExt.jl index e2dde3ce1..daa1d2d39 100644 --- a/ext/OptimizationEnzymeExt.jl +++ b/ext/OptimizationEnzymeExt.jl @@ -6,7 +6,7 @@ import Optimization.LinearAlgebra: I import Optimization.ADTypes: AutoEnzyme isdefined(Base, :get_extension) ? (using Enzyme) : (using ..Enzyme) -@inline function firstapply(f::F, θ, p, args...) where F +@inline function firstapply(f::F, θ, p, args...) where {F} res = f(θ, p, args...) if isa(res, AbstractFloat) res @@ -18,9 +18,8 @@ end function Optimization.instantiate_function(f::OptimizationFunction{true}, x, adtype::AutoEnzyme, p, num_cons = 0) - if f.grad === nothing - grad = let + grad = let function (res, θ, args...) res .= zero(eltype(res)) Enzyme.autodiff(Enzyme.Reverse, diff --git a/ext/OptimizationReverseDiffExt.jl b/ext/OptimizationReverseDiffExt.jl index 87e2c140e..d72640285 100644 --- a/ext/OptimizationReverseDiffExt.jl +++ b/ext/OptimizationReverseDiffExt.jl @@ -21,7 +21,7 @@ function Optimization.instantiate_function(f, x, adtype::AutoReverseDiff, p = SciMLBase.NullParameters(), num_cons = 0) _f = (θ, args...) -> first(f.f(θ, p, args...)) - + chunksize = default_chunk_size(length(x)) if f.grad === nothing @@ -33,7 +33,10 @@ function Optimization.instantiate_function(f, x, adtype::AutoReverseDiff, end else cfg = ReverseDiff.GradientConfig(x) - grad = (res, θ, args...) -> ReverseDiff.gradient!(res, x -> _f(x, args...), θ, cfg) + grad = (res, θ, args...) -> ReverseDiff.gradient!(res, + x -> _f(x, args...), + θ, + cfg) end else grad = (G, θ, args...) -> f.grad(G, θ, p, args...) @@ -41,8 +44,12 @@ function Optimization.instantiate_function(f, x, adtype::AutoReverseDiff, if f.hess === nothing if adtype.compile - T = ForwardDiff.Tag(OptimizationReverseDiffTag(),eltype(x)) - xdual = ForwardDiff.Dual{typeof(T),eltype(x),chunksize}.(x, Ref(ForwardDiff.Partials((ones(eltype(x), chunksize)...,)))) + T = ForwardDiff.Tag(OptimizationReverseDiffTag(), eltype(x)) + xdual = ForwardDiff.Dual{ + typeof(T), + eltype(x), + chunksize, + }.(x, Ref(ForwardDiff.Partials((ones(eltype(x), chunksize)...,)))) h_tape = ReverseDiff.GradientTape(_f, xdual) htape = ReverseDiff.compile(h_tape) function g(θ) @@ -110,7 +117,10 @@ function Optimization.instantiate_function(f, x, adtype::AutoReverseDiff, ReverseDiff.gradient!(res1, htape, θ) end gs = [x -> grad_cons(x, conshtapes[i]) for i in 1:num_cons] - jaccfgs = [ForwardDiff.JacobianConfig(gs[i], x, ForwardDiff.Chunk{chunksize}(), T) for i in 1:num_cons] + jaccfgs = [ForwardDiff.JacobianConfig(gs[i], + x, + ForwardDiff.Chunk{chunksize}(), + T) for i in 1:num_cons] cons_h = function (res, θ) for i in 1:num_cons ForwardDiff.jacobian!(res[i], gs[i], θ, jaccfgs[i], Val{false}()) @@ -155,7 +165,10 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, end else cfg = ReverseDiff.GradientConfig(cache.u0) - grad = (res, θ, args...) -> ReverseDiff.gradient!(res, x -> _f(x, args...), θ, cfg) + grad = (res, θ, args...) -> ReverseDiff.gradient!(res, + x -> _f(x, args...), + θ, + cfg) end else grad = (G, θ, args...) -> f.grad(G, θ, cache.p, args...) @@ -163,15 +176,22 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, if f.hess === nothing if adtype.compile - T = ForwardDiff.Tag(OptimizationReverseDiffTag(),eltype(cache.u0)) - xdual = ForwardDiff.Dual{typeof(T),eltype(cache.u0),chunksize}.(cache.u0, Ref(ForwardDiff.Partials((ones(eltype(cache.u0), chunksize)...,)))) + T = ForwardDiff.Tag(OptimizationReverseDiffTag(), eltype(cache.u0)) + xdual = ForwardDiff.Dual{ + typeof(T), + eltype(cache.u0), + chunksize, + }.(cache.u0, Ref(ForwardDiff.Partials((ones(eltype(cache.u0), chunksize)...,)))) h_tape = ReverseDiff.GradientTape(_f, xdual) htape = ReverseDiff.compile(h_tape) function g(θ) res1 = zeros(eltype(θ), length(θ)) ReverseDiff.gradient!(res1, htape, θ) end - jaccfg = ForwardDiff.JacobianConfig(g, cache.u0, ForwardDiff.Chunk{chunksize}(), T) + jaccfg = ForwardDiff.JacobianConfig(g, + cache.u0, + ForwardDiff.Chunk{chunksize}(), + T) hess = function (res, θ, args...) ForwardDiff.jacobian!(res, g, θ, jaccfg, Val{false}()) end @@ -232,7 +252,10 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, ReverseDiff.gradient!(res1, htape, θ) end gs = [x -> grad_cons(x, conshtapes[i]) for i in 1:num_cons] - jaccfgs = [ForwardDiff.JacobianConfig(gs[i], cache.u0, ForwardDiff.Chunk{chunksize}(), T) for i in 1:num_cons] + jaccfgs = [ForwardDiff.JacobianConfig(gs[i], + cache.u0, + ForwardDiff.Chunk{chunksize}(), + T) for i in 1:num_cons] cons_h = function (res, θ) for i in 1:num_cons ForwardDiff.jacobian!(res[i], gs[i], θ, jaccfgs[i], Val{false}()) diff --git a/ext/OptimizationSparseDiffExt.jl b/ext/OptimizationSparseDiffExt.jl index a115b2710..888b71d38 100644 --- a/ext/OptimizationSparseDiffExt.jl +++ b/ext/OptimizationSparseDiffExt.jl @@ -2,10 +2,14 @@ module OptimizationSparseDiffExt import Optimization, Optimization.ArrayInterface import Optimization.SciMLBase: OptimizationFunction -import Optimization.ADTypes: AutoSparseForwardDiff, AutoSparseFiniteDiff, AutoSparseReverseDiff +import Optimization.ADTypes: AutoSparseForwardDiff, + AutoSparseFiniteDiff, AutoSparseReverseDiff using Optimization.LinearAlgebra, ReverseDiff -isdefined(Base, :get_extension) ? (using SparseDiffTools, SparseDiffTools.ForwardDiff, SparseDiffTools.FiniteDiff, Symbolics) : -(using ..SparseDiffTools, ..SparseDiffTools.ForwardDiff, ..SparseDiffTools.FiniteDiff, ..Symbolics) +isdefined(Base, :get_extension) ? +(using SparseDiffTools, + SparseDiffTools.ForwardDiff, SparseDiffTools.FiniteDiff, Symbolics) : +(using ..SparseDiffTools, + ..SparseDiffTools.ForwardDiff, ..SparseDiffTools.FiniteDiff, ..Symbolics) function default_chunk_size(len) if len < ForwardDiff.DEFAULT_CHUNK_THRESHOLD @@ -247,7 +251,7 @@ function Optimization.instantiate_function(f, x, adtype::AutoSparseFiniteDiff, p else grad = (G, θ, args...) -> f.grad(G, θ, p, args...) end - + hess_sparsity = f.hess_prototype hess_colors = f.hess_colorvec if f.hess === nothing @@ -371,7 +375,7 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, else grad = (G, θ, args...) -> f.grad(G, θ, cache.p, args...) end - + hess_sparsity = f.hess_prototype hess_colors = f.hess_colorvec if f.hess === nothing @@ -406,7 +410,8 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, cons_jac_colorvec = f.cons_jac_colorvec if cons !== nothing && f.cons_j === nothing cons_jac_prototype = f.cons_jac_prototype === nothing ? - Symbolics.jacobian_sparsity(cons, zeros(eltype(cache.u0), num_cons), + Symbolics.jacobian_sparsity(cons, + zeros(eltype(cache.u0), num_cons), cache.u0) : f.cons_jac_prototype cons_jac_colorvec = f.cons_jac_colorvec === nothing ? @@ -503,7 +508,10 @@ function Optimization.instantiate_function(f, x, adtype::AutoSparseReverseDiff, end else cfg = ReverseDiff.GradientConfig(x) - grad = (res, θ, args...) -> ReverseDiff.gradient!(res, x -> _f(x, args...), θ, cfg) + grad = (res, θ, args...) -> ReverseDiff.gradient!(res, + x -> _f(x, args...), + θ, + cfg) end else grad = (G, θ, args...) -> f.grad(G, θ, p, args...) @@ -515,20 +523,32 @@ function Optimization.instantiate_function(f, x, adtype::AutoSparseReverseDiff, hess_sparsity = Symbolics.hessian_sparsity(_f, x) hess_colors = SparseDiffTools.matrix_colors(tril(hess_sparsity)) if adtype.compile - T = ForwardDiff.Tag(OptimizationSparseReverseTag(),eltype(x)) - xdual = ForwardDiff.Dual{typeof(T),eltype(x),min(chunksize, maximum(hess_colors))}.(x, Ref(ForwardDiff.Partials((ones(eltype(x), min(chunksize, maximum(hess_colors)))...,)))) + T = ForwardDiff.Tag(OptimizationSparseReverseTag(), eltype(x)) + xdual = ForwardDiff.Dual{ + typeof(T), + eltype(x), + min(chunksize, maximum(hess_colors)), + }.(x, + Ref(ForwardDiff.Partials((ones(eltype(x), + min(chunksize, maximum(hess_colors)))...,)))) h_tape = ReverseDiff.GradientTape(_f, xdual) htape = ReverseDiff.compile(h_tape) function g(res1, θ) ReverseDiff.gradient!(res1, htape, θ) end - jaccfg = ForwardColorJacCache(g, x; tag = typeof(T), colorvec = hess_colors, sparsity = hess_sparsity) + jaccfg = ForwardColorJacCache(g, + x; + tag = typeof(T), + colorvec = hess_colors, + sparsity = hess_sparsity) hess = function (res, θ, args...) SparseDiffTools.forwarddiff_color_jacobian!(res, g, θ, jaccfg) end else hess = function (res, θ, args...) - res .= SparseDiffTools.forwarddiff_color_jacobian(θ, colorvec = hess_colors, sparsity = hess_sparsity) do θ + res .= SparseDiffTools.forwarddiff_color_jacobian(θ, + colorvec = hess_colors, + sparsity = hess_sparsity) do θ ReverseDiff.gradient(x -> _f(x, args...), θ) end end @@ -561,7 +581,11 @@ function Optimization.instantiate_function(f, x, adtype::AutoSparseReverseDiff, cons_jac_prototype = f.cons_jac_prototype cons_jac_colorvec = f.cons_jac_colorvec if cons !== nothing && f.cons_j === nothing - jaccache = SparseDiffTools.sparse_jacobian_cache(AutoSparseForwardDiff(), SparseDiffTools.SymbolicsSparsityDetection(), cons_oop, x, fx = zeros(eltype(x), num_cons)) + jaccache = SparseDiffTools.sparse_jacobian_cache(AutoSparseForwardDiff(), + SparseDiffTools.SymbolicsSparsityDetection(), + cons_oop, + x, + fx = zeros(eltype(x), num_cons)) # let cons = cons, θ = cache.u0, cons_jac_colorvec = cons_jac_colorvec, cons_jac_prototype = cons_jac_prototype, num_cons = num_cons # ForwardColorJacCache(cons, θ; # colorvec = cons_jac_colorvec, @@ -570,14 +594,14 @@ function Optimization.instantiate_function(f, x, adtype::AutoSparseReverseDiff, # end cons_jac_prototype = jaccache.jac_prototype cons_jac_colorvec = jaccache.coloring - cons_j = function (J, θ, args...;cons = cons, cache = jaccache.cache) + cons_j = function (J, θ, args...; cons = cons, cache = jaccache.cache) forwarddiff_color_jacobian!(J, cons, θ, cache) return end else cons_j = (J, θ) -> f.cons_j(J, θ, p) end - + conshess_sparsity = f.cons_hess_prototype conshess_colors = f.cons_hess_colorvec if cons !== nothing && f.cons_h === nothing @@ -585,24 +609,39 @@ function Optimization.instantiate_function(f, x, adtype::AutoSparseReverseDiff, conshess_sparsity = Symbolics.hessian_sparsity.(fncs, Ref(x)) conshess_colors = SparseDiffTools.matrix_colors.(conshess_sparsity) if adtype.compile - T = ForwardDiff.Tag(OptimizationSparseReverseTag(),eltype(x)) - xduals = [ForwardDiff.Dual{typeof(T),eltype(x),min(chunksize, maximum(conshess_colors[i]))}.(x, Ref(ForwardDiff.Partials((ones(eltype(x), min(chunksize, maximum(conshess_colors[i])))...,)))) for i in 1:num_cons] - consh_tapes = [ReverseDiff.GradientTape(fncs[i], xduals[i]) for i in 1:num_cons] + T = ForwardDiff.Tag(OptimizationSparseReverseTag(), eltype(x)) + xduals = [ForwardDiff.Dual{ + typeof(T), + eltype(x), + min(chunksize, maximum(conshess_colors[i])), + }.(x, + Ref(ForwardDiff.Partials((ones(eltype(x), + min(chunksize, maximum(conshess_colors[i])))...,)))) for i in 1:num_cons] + consh_tapes = [ReverseDiff.GradientTape(fncs[i], xduals[i]) for i in 1:num_cons] conshtapes = ReverseDiff.compile.(consh_tapes) function grad_cons(res1, θ, htape) ReverseDiff.gradient!(res1, htape, θ) end gs = [(res1, x) -> grad_cons(res1, x, conshtapes[i]) for i in 1:num_cons] - jaccfgs = [ForwardColorJacCache(gs[i], x; tag = typeof(T), colorvec = conshess_colors[i], sparsity = conshess_sparsity[i]) for i in 1:num_cons] + jaccfgs = [ForwardColorJacCache(gs[i], + x; + tag = typeof(T), + colorvec = conshess_colors[i], + sparsity = conshess_sparsity[i]) for i in 1:num_cons] cons_h = function (res, θ, args...) for i in 1:num_cons - SparseDiffTools.forwarddiff_color_jacobian!(res[i], gs[i], θ, jaccfgs[i]) + SparseDiffTools.forwarddiff_color_jacobian!(res[i], + gs[i], + θ, + jaccfgs[i]) end end else cons_h = function (res, θ) for i in 1:num_cons - res[i] .= SparseDiffTools.forwarddiff_color_jacobian(θ, colorvec = conshess_colors[i], sparsity = conshess_sparsity[i]) do θ + res[i] .= SparseDiffTools.forwarddiff_color_jacobian(θ, + colorvec = conshess_colors[i], + sparsity = conshess_sparsity[i]) do θ ReverseDiff.gradient(fncs[i], θ) end end @@ -643,7 +682,10 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, end else cfg = ReverseDiff.GradientConfig(cache.u0) - grad = (res, θ, args...) -> ReverseDiff.gradient!(res, x -> _f(x, args...), θ, cfg) + grad = (res, θ, args...) -> ReverseDiff.gradient!(res, + x -> _f(x, args...), + θ, + cfg) end else grad = (G, θ, args...) -> f.grad(G, θ, cache.p, args...) @@ -655,20 +697,32 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, hess_sparsity = Symbolics.hessian_sparsity(_f, cache.u0) hess_colors = SparseDiffTools.matrix_colors(tril(hess_sparsity)) if adtype.compile - T = ForwardDiff.Tag(OptimizationSparseReverseTag(),eltype(cache.u0)) - xdual = ForwardDiff.Dual{typeof(T),eltype(cache.u0),min(chunksize, maximum(hess_colors))}.(cache.u0, Ref(ForwardDiff.Partials((ones(eltype(cache.u0), min(chunksize, maximum(hess_colors)))...,)))) + T = ForwardDiff.Tag(OptimizationSparseReverseTag(), eltype(cache.u0)) + xdual = ForwardDiff.Dual{ + typeof(T), + eltype(cache.u0), + min(chunksize, maximum(hess_colors)), + }.(cache.u0, + Ref(ForwardDiff.Partials((ones(eltype(cache.u0), + min(chunksize, maximum(hess_colors)))...,)))) h_tape = ReverseDiff.GradientTape(_f, xdual) htape = ReverseDiff.compile(h_tape) function g(res1, θ) ReverseDiff.gradient!(res1, htape, θ) end - jaccfg = ForwardColorJacCache(g, cache.u0; tag = typeof(T), colorvec = hess_colors, sparsity = hess_sparsity) + jaccfg = ForwardColorJacCache(g, + cache.u0; + tag = typeof(T), + colorvec = hess_colors, + sparsity = hess_sparsity) hess = function (res, θ, args...) SparseDiffTools.forwarddiff_color_jacobian!(res, g, θ, jaccfg) end else hess = function (res, θ, args...) - res .= SparseDiffTools.forwarddiff_color_jacobian(θ, colorvec = hess_colors, sparsity = hess_sparsity) do θ + res .= SparseDiffTools.forwarddiff_color_jacobian(θ, + colorvec = hess_colors, + sparsity = hess_sparsity) do θ ReverseDiff.gradient(x -> _f(x, args...), θ) end end @@ -694,7 +748,7 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, if f.cons === nothing cons = nothing else - cons = function (res, θ) + cons = function (res, θ) f.cons(res, θ, cache.p) return end @@ -705,10 +759,14 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, cons_jac_colorvec = f.cons_jac_colorvec if cons !== nothing && f.cons_j === nothing # cons_jac_prototype = Symbolics.jacobian_sparsity(cons, - # zeros(eltype(cache.u0), num_cons), - # cache.u0) + # zeros(eltype(cache.u0), num_cons), + # cache.u0) # cons_jac_colorvec = matrix_colors(cons_jac_prototype) - jaccache = SparseDiffTools.sparse_jacobian_cache(AutoSparseForwardDiff(), SparseDiffTools.SymbolicsSparsityDetection(), cons_oop, cache.u0, fx = zeros(eltype(cache.u0), num_cons)) + jaccache = SparseDiffTools.sparse_jacobian_cache(AutoSparseForwardDiff(), + SparseDiffTools.SymbolicsSparsityDetection(), + cons_oop, + cache.u0, + fx = zeros(eltype(cache.u0), num_cons)) # let cons = cons, θ = cache.u0, cons_jac_colorvec = cons_jac_colorvec, cons_jac_prototype = cons_jac_prototype, num_cons = num_cons # ForwardColorJacCache(cons, θ; # colorvec = cons_jac_colorvec, @@ -724,7 +782,7 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, else cons_j = (J, θ) -> f.cons_j(J, θ, cache.p) end - + conshess_sparsity = f.cons_hess_prototype conshess_colors = f.cons_hess_colorvec if cons !== nothing && f.cons_h === nothing @@ -742,30 +800,45 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, end conshess_colors = SparseDiffTools.matrix_colors.(conshess_sparsity) if adtype.compile - T = ForwardDiff.Tag(OptimizationSparseReverseTag(),eltype(cache.u0)) - xduals = [ForwardDiff.Dual{typeof(T),eltype(cache.u0),min(chunksize, maximum(conshess_colors[i]))}.(cache.u0, Ref(ForwardDiff.Partials((ones(eltype(cache.u0), min(chunksize, maximum(conshess_colors[i])))...,)))) for i in 1:num_cons] - consh_tapes = [ReverseDiff.GradientTape(fncs[i], xduals[i]) for i in 1:num_cons] + T = ForwardDiff.Tag(OptimizationSparseReverseTag(), eltype(cache.u0)) + xduals = [ForwardDiff.Dual{ + typeof(T), + eltype(cache.u0), + min(chunksize, maximum(conshess_colors[i])), + }.(cache.u0, + Ref(ForwardDiff.Partials((ones(eltype(cache.u0), + min(chunksize, maximum(conshess_colors[i])))...,)))) for i in 1:num_cons] + consh_tapes = [ReverseDiff.GradientTape(fncs[i], xduals[i]) for i in 1:num_cons] conshtapes = ReverseDiff.compile.(consh_tapes) function grad_cons(res1, θ, htape) ReverseDiff.gradient!(res1, htape, θ) end gs = let conshtapes = conshtapes - map(1:num_cons) do i - function (res1, x) + map(1:num_cons) do i + function (res1, x) grad_cons(res1, x, conshtapes[i]) end end end - jaccfgs = [ForwardColorJacCache(gs[i], cache.u0; tag = typeof(T), colorvec = conshess_colors[i], sparsity = conshess_sparsity[i]) for i in 1:num_cons] + jaccfgs = [ForwardColorJacCache(gs[i], + cache.u0; + tag = typeof(T), + colorvec = conshess_colors[i], + sparsity = conshess_sparsity[i]) for i in 1:num_cons] cons_h = function (res, θ) for i in 1:num_cons - SparseDiffTools.forwarddiff_color_jacobian!(res[i], gs[i], θ, jaccfgs[i]) + SparseDiffTools.forwarddiff_color_jacobian!(res[i], + gs[i], + θ, + jaccfgs[i]) end end else cons_h = function (res, θ) for i in 1:num_cons - res[i] .= SparseDiffTools.forwarddiff_color_jacobian(θ, colorvec = conshess_colors[i], sparsity = conshess_sparsity[i]) do θ + res[i] .= SparseDiffTools.forwarddiff_color_jacobian(θ, + colorvec = conshess_colors[i], + sparsity = conshess_sparsity[i]) do θ ReverseDiff.gradient(fncs[i], θ) end end diff --git a/ext/OptimizationZygoteExt.jl b/ext/OptimizationZygoteExt.jl index 250145f0d..3725466ce 100644 --- a/ext/OptimizationZygoteExt.jl +++ b/ext/OptimizationZygoteExt.jl @@ -85,7 +85,6 @@ end function Optimization.instantiate_function(f, cache::Optimization.ReInitCache, adtype::AutoZygote, num_cons = 0) - _f = (θ, args...) -> f(θ, cache.p, args...)[1] if f.grad === nothing grad = function (res, θ, args...) diff --git a/lib/OptimizationBBO/src/OptimizationBBO.jl b/lib/OptimizationBBO/src/OptimizationBBO.jl index 9ef8ccd61..f4d3b7eb9 100644 --- a/lib/OptimizationBBO/src/OptimizationBBO.jl +++ b/lib/OptimizationBBO/src/OptimizationBBO.jl @@ -135,7 +135,8 @@ function SciMLBase.__solve(cache::Optimization.OptimizationCache{ maxtime = Optimization._check_and_convert_maxtime(cache.solver_args.maxtime) _loss = function (θ) - if cache.callback === Optimization.DEFAULT_CALLBACK && cache.data === Optimization.DEFAULT_DATA + if cache.callback === Optimization.DEFAULT_CALLBACK && + cache.data === Optimization.DEFAULT_DATA return first(cache.f(θ, cache.p)) elseif cache.callback === Optimization.DEFAULT_CALLBACK return first(cache.f(θ, cache.p, cur...)) diff --git a/lib/OptimizationMOI/src/nlp.jl b/lib/OptimizationMOI/src/nlp.jl index d23df8a91..055bf91d2 100644 --- a/lib/OptimizationMOI/src/nlp.jl +++ b/lib/OptimizationMOI/src/nlp.jl @@ -102,7 +102,10 @@ function SciMLBase.get_paramsyms(sol::SciMLBase.OptimizationSolution{ sol.cache.evaluator.f.paramsyms end -function MOIOptimizationNLPCache(prob::OptimizationProblem, opt; callback = nothing, kwargs...) +function MOIOptimizationNLPCache(prob::OptimizationProblem, + opt; + callback = nothing, + kwargs...) reinit_cache = Optimization.ReInitCache(prob.u0, prob.p) # everything that can be changed via `reinit` num_cons = prob.ucons === nothing ? 0 : length(prob.ucons) @@ -418,7 +421,7 @@ function SciMLBase.__solve(cache::MOIOptimizationNLPCache) if cache.evaluator.callback !== nothing MOI.set(opt_setup, MOI.Silent(), true) end - + MOI.optimize!(opt_setup) if MOI.get(opt_setup, MOI.ResultCount()) >= 1 minimizer = MOI.get(opt_setup, MOI.VariablePrimal(), θ) diff --git a/lib/OptimizationMOI/test/runtests.jl b/lib/OptimizationMOI/test/runtests.jl index ba04cae51..75386710e 100644 --- a/lib/OptimizationMOI/test/runtests.jl +++ b/lib/OptimizationMOI/test/runtests.jl @@ -40,11 +40,11 @@ end callback = function (p, l) global iter iter += 1 - + display(l) return false end - + sol = solve(prob, Ipopt.Optimizer(); callback) @test 10 * sol.objective < l1 diff --git a/lib/OptimizationOptimisers/src/sophia.jl b/lib/OptimizationOptimisers/src/sophia.jl index fa9fa827b..ada54090e 100644 --- a/lib/OptimizationOptimisers/src/sophia.jl +++ b/lib/OptimizationOptimisers/src/sophia.jl @@ -71,7 +71,8 @@ function SciMLBase.__solve(cache::OptimizationCache{ maxiters = Optimization._check_and_convert_maxiters(maxiters) _loss = function (θ) - if cache.callback === Optimization.DEFAULT_CALLBACK && data === Optimization.DEFAULT_DATA + if cache.callback === Optimization.DEFAULT_CALLBACK && + data === Optimization.DEFAULT_DATA return first(cache.f(θ, cache.p)) elseif cache.callback === Optimization.DEFAULT_CALLBACK return first(cache.f(θ, cache.p, cur...)) diff --git a/lib/OptimizationPRIMA/LICENSE b/lib/OptimizationPRIMA/LICENSE new file mode 100644 index 000000000..0922eea00 --- /dev/null +++ b/lib/OptimizationPRIMA/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2023 Vaibhav Dixit and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/lib/OptimizationPRIMA/Project.toml b/lib/OptimizationPRIMA/Project.toml new file mode 100644 index 000000000..c1366d96e --- /dev/null +++ b/lib/OptimizationPRIMA/Project.toml @@ -0,0 +1,21 @@ +name = "OptimizationPRIMA" +uuid = "72f8369c-a2ea-4298-9126-56167ce9cbc2" +authors = ["Vaibhav Dixit and contributors"] +version = "0.0.1" + +[deps] +Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" + + +[compat] +julia = "1" + +[extras] +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test", "ForwardDiff", "ModelingToolkit"] diff --git a/lib/OptimizationPRIMA/src/OptimizationPRIMA.jl b/lib/OptimizationPRIMA/src/OptimizationPRIMA.jl new file mode 100644 index 000000000..9ca20ec5b --- /dev/null +++ b/lib/OptimizationPRIMA/src/OptimizationPRIMA.jl @@ -0,0 +1,188 @@ +module OptimizationPRIMA + +using Optimization, Optimization.SciMLBase, Reexport +@reexport using PRIMA + +abstract type PRIMASolvers end + +struct UOBYQA <: PRIMASolvers end +struct NEWUOA <: PRIMASolvers end +struct BOBYQA <: PRIMASolvers end +struct LINCOA <: PRIMASolvers end +struct COBYLA <: PRIMASolvers end + +SciMLBase.supports_opt_cache_interface(::PRIMASolvers) = true +SciMLBase.allowsconstraints(::Union{LINCOA, COBYLA}) = true +SciMLBase.allowsbounds(opt::Union{BOBYQA, LINCOA, COBYLA}) = true +SciMLBase.requiresconstraints(opt::COBYLA) = true + +function Optimization.OptimizationCache(prob::SciMLBase.OptimizationProblem, opt::PRIMASolvers, data; + callback = Optimization.DEFAULT_CALLBACK, + maxiters::Union{Number, Nothing} = nothing, + maxtime::Union{Number, Nothing} = nothing, + abstol::Union{Number, Nothing} = nothing, + reltol::Union{Number, Nothing} = nothing, + progress = false, + kwargs...) + reinit_cache = Optimization.ReInitCache(prob.u0, prob.p) + num_cons = prob.ucons === nothing ? 0 : length(prob.ucons) + if prob.f.adtype isa SciMLBase.NoAD && opt isa COBYLA + throw("We evaluate the jacobian and hessian of the constraints once to automatically detect + linear and nonlinear constraints, please provide a valid AD backend for using COBYLA.") + else + f = Optimization.instantiate_function(prob.f, reinit_cache, prob.f.adtype, num_cons) + end + + return Optimization.OptimizationCache(f, reinit_cache, prob.lb, prob.ub, prob.lcons, + prob.ucons, prob.sense, + opt, data, progress, callback, + merge((; maxiters, maxtime, abstol, reltol), + NamedTuple(kwargs))) +end + +function get_solve_func(opt::PRIMASolvers) + if opt isa UOBYQA + return PRIMA.uobyqa + elseif opt isa NEWUOA + return PRIMA.newuoa + elseif opt isa BOBYQA + return PRIMA.bobyqa + elseif opt isa LINCOA + return PRIMA.lincoa + elseif opt isa COBYLA + return PRIMA.cobyla + end +end + +function __map_optimizer_args!(cache::Optimization.OptimizationCache, opt::PRIMASolvers; + callback = nothing, + maxiters::Union{Number, Nothing} = nothing, + maxtime::Union{Number, Nothing} = nothing, + abstol::Union{Number, Nothing} = nothing, + reltol::Union{Number, Nothing} = nothing, + kwargs...) + kws = (; kwargs...) + + if !isnothing(maxiters) + kws = (; kws..., maxfun = maxiters) + end + + if cache.ub !== nothing + kws = (; kws..., xu = cache.ub, xl = cache.lb) + end + + if !isnothing(maxtime) || !isnothing(abstol) || !isnothing(reltol) + error("maxtime, abstol and reltol kwargs not supported in $opt") + end + + return kws +end + +function sciml_prima_retcode(rc::AbstractString) + if rc in ["SMALL_TR_RADIUS", "TRSUBP_FAILED", "NAN_INF_X", "NAN_INF_F", "NAN_INF_MODEL", + "DAMAGING_ROUNDING", "ZERO_LINEAR_CONSTRAINT", "INVALID_INPUT", "ASSERTION_FAILS", + "VALIDATION_FAILS", "MEMORY_ALLOCATION_FAILS"] + return ReturnCode.Failure + else + rc in ["FTARGET_ACHIEVED" + "MAXFUN_REACHED" + "MAXTR_REACHED" + "NO_SPACE_BETWEEN_BOUNDS"] + return ReturnCode.Success + end +end + +function SciMLBase.__solve(cache::Optimization.OptimizationCache{ + F, + RC, + LB, + UB, + LC, + UC, + S, + O, + D, + P, + C, +}) where { + F, + RC, + LB, + UB, + LC, + UC, + S, + O <: PRIMASolvers, + D, + P, + C, +} + _loss = function (θ) + x = cache.f(θ, cache.p) + if cache.callback(θ, x...) + error("Optimization halted by callback.") + end + return x[1] + end + + optfunc = get_solve_func(cache.opt) + + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) + maxtime = Optimization._check_and_convert_maxtime(cache.solver_args.maxtime) + + kws = __map_optimizer_args!(cache, cache.opt; callback = cache.callback, + maxiters = maxiters, + maxtime = maxtime, + cache.solver_args...) + + t0 = time() + if cache.opt isa COBYLA + lineqsinds = Int[] + linineqsinds = Int[] + nonlininds = Int[] + H = [zeros(length(cache.u0), length(cache.u0)) for i in 1:length(cache.lcons)] + J = zeros(length(cache.lcons), length(cache.u0)) + + cache.f.cons_h(H, ones(length(cache.u0))) + cache.f.cons_j(J, ones(length(cache.u0))) + for i in eachindex(cache.lcons) + if iszero(H[i]) && cache.lcons[i] == cache.ucons[i] + push!(lineqsinds, i) + elseif iszero(H[i]) && cache.lcons[i] != cache.ucons[i] + push!(linineqsinds, i) + else + push!(nonlininds, i) + end + end + res1 = zeros(length(cache.lcons)) + nonlincons = (res, θ) -> (cache.f.cons(res1, θ); res .= res1[nonlininds]) + A₁ = J[lineqsinds, :] + b₁ = cache.lcons[lineqsinds] + A₂ = J[linineqsinds, :] + b₂ = cache.ucons[linineqsinds] + function fwcons(θ, res) + nonlincons(res, θ) + return _loss(θ) + end + (minx, minf, nf, rc, cstrv) = optfunc(fwcons, + cache.u0; + linear_eq = (A₁, b₁), + linear_ineq = (A₂, b₂), + nonlinear_ineq = length(nonlininds), + kws...) + elseif cache.opt isa LINCOA + (minx, minf, nf, rc, cstrv) = optfunc(_loss, cache.u0; kws...) + else + (minx, minf, nf, rc) = optfunc(_loss, cache.u0; kws...) + end + t1 = time() + + retcode = sciml_prima_retcode(PRIMA.reason(rc)) + + SciMLBase.build_solution(cache, cache.opt, minx, + minf; retcode = retcode, + solve_time = t1 - t0) +end + +export UOBYQA, NEWUOA, BOBYQA, LINCOA, COBYLA +end diff --git a/lib/OptimizationPRIMA/test/runtests.jl b/lib/OptimizationPRIMA/test/runtests.jl new file mode 100644 index 000000000..0d483bf9e --- /dev/null +++ b/lib/OptimizationPRIMA/test/runtests.jl @@ -0,0 +1,50 @@ +using OptimizationPRIMA, Optimization, ForwardDiff, ModelingToolkit +using Test + +@testset "OptimizationPRIMA.jl" begin + rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 + x0 = zeros(2) + _p = [1.0, 100.0] + l1 = rosenbrock(x0, _p) + + prob = OptimizationProblem(rosenbrock, x0, _p) + sol = Optimization.solve(prob, UOBYQA(), maxiters = 1000) + @test 10 * sol.objective < l1 + sol = Optimization.solve(prob, NEWUOA(), maxiters = 1000) + @test 10 * sol.objective < l1 + sol = Optimization.solve(prob, BOBYQA(), maxiters = 1000) + @test 10 * sol.objective < l1 + sol = Optimization.solve(prob, LINCOA(), maxiters = 1000) + @test 10 * sol.objective < l1 + @test_throws SciMLBase.IncompatibleOptimizerError Optimization.solve(prob, + COBYLA(), + maxiters = 1000) + + function con2_c(res, x, p) + res .= [x[1] + x[2], x[2] * sin(x[1]) - x[1]] + end + optprob = OptimizationFunction(rosenbrock, AutoForwardDiff(), cons = con2_c) + prob = OptimizationProblem(optprob, x0, _p, lcons = [1, -100], ucons = [1, 100]) + sol = Optimization.solve(prob, COBYLA(), maxiters = 1000) + @test sol.objective < l1 + + function con2_c(res, x, p) + res .= [x[1] + x[2]] + end + optprob = OptimizationFunction(rosenbrock, AutoForwardDiff(), cons = con2_c) + prob = OptimizationProblem(optprob, x0, _p, lcons = [1], ucons = [1]) + sol = Optimization.solve(prob, COBYLA(), maxiters = 1000) + @test sol.objective < l1 + + prob = OptimizationProblem(optprob, x0, _p, lcons = [1], ucons = [5]) + sol = Optimization.solve(prob, COBYLA(), maxiters = 1000) + @test sol.objective < l1 + + function con2_c(res, x, p) + res .= [x[2] * sin(x[1]) - x[1]] + end + optprob = OptimizationFunction(rosenbrock, AutoModelingToolkit(), cons = con2_c) + prob = OptimizationProblem(optprob, x0, _p, lcons = [10], ucons = [50]) + sol = Optimization.solve(prob, COBYLA(), maxiters = 1000) + @test 10 * sol.objective < l1 +end diff --git a/test/AD_performance_regression.jl b/test/AD_performance_regression.jl index 66291f029..328003590 100644 --- a/test/AD_performance_regression.jl +++ b/test/AD_performance_regression.jl @@ -3,37 +3,152 @@ using ReverseDiff, Enzyme, BenchmarkTools, Test lookup_pg = Dict(5 => 11, 4 => 13, 2 => 15, 3 => 17, 1 => 19) ref_gen_idxs = [5, 4, 2, 3, 1] -cost_arrs = Dict(5 => [0.0, 1000.0, 0.0], 4 => [0.0, 4000.0, 0.0], 2 => [0.0, 1500.0, 0.0], 3 => [0.0, 3000.0, 0.0], 1 => [0.0, 1400.0, 0.0]) +cost_arrs = Dict(5 => [0.0, 1000.0, 0.0], + 4 => [0.0, 4000.0, 0.0], + 2 => [0.0, 1500.0, 0.0], + 3 => [0.0, 3000.0, 0.0], + 1 => [0.0, 1400.0, 0.0]) + +opf_objective = let lookup_pg = lookup_pg, ref_gen_idxs = ref_gen_idxs, + cost_arrs = cost_arrs -opf_objective = let lookup_pg=lookup_pg, ref_gen_idxs=ref_gen_idxs, cost_arrs=cost_arrs function (x, _) #start = time() cost = 0.0 for i in ref_gen_idxs pg = x[lookup_pg[i]] _cost_arr = cost_arrs[i] - cost += _cost_arr[1]*pg^2 + _cost_arr[2]*pg + _cost_arr[3] + cost += _cost_arr[1] * pg^2 + _cost_arr[2] * pg + _cost_arr[3] end #total_callback_time += time() - start return cost end end -optprob = Optimization.OptimizationFunction(opf_objective, Optimization.AutoReverseDiff(true)) +optprob = Optimization.OptimizationFunction(opf_objective, + Optimization.AutoReverseDiff(true)) -test_u0 = [0.6292298794022337, 0.30740951571225206, 0.0215258802699263, 0.38457509230779996, 0.9419186480931858, 0.34961116773074874, 0.875763562401991, 0.3203478635827923, 0.6354060958226175, 0.45537545721771266, 0.3120599359696674, 0.2421238802331842, 0.886455177641366, 0.49797378087768696, 0.652913329799645, 0.03590201299300255, 0.5618806749518928, 0.8142146688533769, 0.3973557130434364, 0.27827135011662674, 0.16456134856048643, 0.7465018431665373, 0.4898329811551083, 0.6966035226583556, 0.7419662648518377, 0.8505905798503723, 0.27102126066405097, 0.1988238097281576, 0.09684601934490256, 0.49238142828542797, 0.1366594202307445, 0.6337080281764231, 0.28814906958008235, 0.5404996094640431, 0.015153517398975858, 0.6338449294034381, 0.5165464961007717, 0.572879113636733, 0.9652420600585092, 0.26535868365228543, 0.865686920119479, 0.38426996353892773, 0.007412077949221274, 0.3889835001514599] +test_u0 = [ + 0.6292298794022337, + 0.30740951571225206, + 0.0215258802699263, + 0.38457509230779996, + 0.9419186480931858, + 0.34961116773074874, + 0.875763562401991, + 0.3203478635827923, + 0.6354060958226175, + 0.45537545721771266, + 0.3120599359696674, + 0.2421238802331842, + 0.886455177641366, + 0.49797378087768696, + 0.652913329799645, + 0.03590201299300255, + 0.5618806749518928, + 0.8142146688533769, + 0.3973557130434364, + 0.27827135011662674, + 0.16456134856048643, + 0.7465018431665373, + 0.4898329811551083, + 0.6966035226583556, + 0.7419662648518377, + 0.8505905798503723, + 0.27102126066405097, + 0.1988238097281576, + 0.09684601934490256, + 0.49238142828542797, + 0.1366594202307445, + 0.6337080281764231, + 0.28814906958008235, + 0.5404996094640431, + 0.015153517398975858, + 0.6338449294034381, + 0.5165464961007717, + 0.572879113636733, + 0.9652420600585092, + 0.26535868365228543, + 0.865686920119479, + 0.38426996353892773, + 0.007412077949221274, + 0.3889835001514599, +] test_obj = 7079.190664351089 -test_cons = [0.0215258802699263, -1.0701734802505833, -5.108902216849063, -3.49724505910433, -2.617834191007569, 0.5457423426033834, -0.7150251969424766, -2.473175092089014, -2.071687022809815, -1.5522321037165985, -1.0107399030803794, 3.0047739260369246, 0.2849522377447594, -2.8227966798520674, 3.2236954017592256, 1.0793383525116511, -1.633412293595111, -3.1618224299953224, -0.7775962590542184, 1.7252573527333024, -4.23535583005632, -1.7030832394691608, 1.5810450617647889, -0.33289810365419437, 0.19476447251065077, 1.0688558672739048, 1.563372246165339, 9.915310272572729, 1.4932615291788414, 2.0016715378998793, -1.4038702698147258, -0.8834081057449231, 0.21730536348839036, -7.40879932706212, -1.6000837514115611, 0.8542376821320647, 0.06615508569119477, -0.6077039991323074, 0.6138802155526912, 0.0061762164203837955, -0.3065125522705683, 0.5843454392910835, 0.7251928172073308, 1.2740182727083802, 0.11298343104675009, 0.2518186223833513, 0.4202616621130535, 0.3751697141306502, 0.4019890236200105, 0.5950107614751935, 1.0021074654956683, 0.897077248544158, 0.15136310228960612] +test_cons = [ + 0.0215258802699263, + -1.0701734802505833, + -5.108902216849063, + -3.49724505910433, + -2.617834191007569, + 0.5457423426033834, + -0.7150251969424766, + -2.473175092089014, + -2.071687022809815, + -1.5522321037165985, + -1.0107399030803794, + 3.0047739260369246, + 0.2849522377447594, + -2.8227966798520674, + 3.2236954017592256, + 1.0793383525116511, + -1.633412293595111, + -3.1618224299953224, + -0.7775962590542184, + 1.7252573527333024, + -4.23535583005632, + -1.7030832394691608, + 1.5810450617647889, + -0.33289810365419437, + 0.19476447251065077, + 1.0688558672739048, + 1.563372246165339, + 9.915310272572729, + 1.4932615291788414, + 2.0016715378998793, + -1.4038702698147258, + -0.8834081057449231, + 0.21730536348839036, + -7.40879932706212, + -1.6000837514115611, + 0.8542376821320647, + 0.06615508569119477, + -0.6077039991323074, + 0.6138802155526912, + 0.0061762164203837955, + -0.3065125522705683, + 0.5843454392910835, + 0.7251928172073308, + 1.2740182727083802, + 0.11298343104675009, + 0.2518186223833513, + 0.4202616621130535, + 0.3751697141306502, + 0.4019890236200105, + 0.5950107614751935, + 1.0021074654956683, + 0.897077248544158, + 0.15136310228960612, +] res = zero(test_u0) -_f = Optimization.instantiate_function(optprob, test_u0, Optimization.AutoReverseDiff(false), nothing) +_f = Optimization.instantiate_function(optprob, + test_u0, + Optimization.AutoReverseDiff(false), + nothing) _f.f(test_u0, nothing) @test @ballocated($(_f.grad)($res, $test_u0)) > 1000 -_f2 = Optimization.instantiate_function(optprob, test_u0, Optimization.AutoReverseDiff(true), nothing) +_f2 = Optimization.instantiate_function(optprob, + test_u0, + Optimization.AutoReverseDiff(true), + nothing) _f2.f(test_u0, nothing) @test @ballocated($(_f2.grad)($res, $test_u0)) == 0 -_f3 = Optimization.instantiate_function(optprob, test_u0, Optimization.AutoEnzyme(), nothing) +_f3 = Optimization.instantiate_function(optprob, + test_u0, + Optimization.AutoEnzyme(), + nothing) _f3.f(test_u0, nothing) -@test @ballocated($(_f3.grad)($res, $test_u0)) == 0 \ No newline at end of file +@test @ballocated($(_f3.grad)($res, $test_u0)) == 0 diff --git a/test/ADtests.jl b/test/ADtests.jl index c948ebf36..7cca09fc0 100644 --- a/test/ADtests.jl +++ b/test/ADtests.jl @@ -85,8 +85,8 @@ if VERSION >= v"1.9" optprob.cons_h(H3, x0) @test H3 == [[2.0 0.0; 0.0 2.0]] -G2 = Array{Float64}(undef, 2) -H2 = Array{Float64}(undef, 2, 2) + G2 = Array{Float64}(undef, 2) + H2 = Array{Float64}(undef, 2, 2) optf = OptimizationFunction(rosenbrock, Optimization.AutoEnzyme(), cons = con2_c) optprob = Optimization.instantiate_function(optf, x0, Optimization.AutoEnzyme(), @@ -126,12 +126,12 @@ H3 = [Array{Float64}(undef, 2, 2), Array{Float64}(undef, 2, 2)] optprob.cons_h(H3, x0) H3 == [[2.0 0.0; 0.0 2.0], [-0.0 1.0; 1.0 0.0]] - G2 = Array{Float64}(undef, 2) H2 = Array{Float64}(undef, 2, 2) optf = OptimizationFunction(rosenbrock, Optimization.AutoReverseDiff(), cons = con2_c) -optprob = Optimization.instantiate_function(optf, x0, Optimization.AutoReverseDiff(compile=true), +optprob = Optimization.instantiate_function(optf, x0, + Optimization.AutoReverseDiff(compile = true), nothing, 2) optprob.grad(G2, x0) @test G1 == G2 @@ -524,7 +524,8 @@ sol = solve(prob, Optimisers.ADAM(0.1), maxiters = 1000) @test 10 * sol.objective < l1 optf = OptimizationFunction(rosenbrock, Optimization.AutoSparseReverseDiff(), cons = con2_c) -optprob = Optimization.instantiate_function(optf, x0, Optimization.AutoSparseReverseDiff(true), +optprob = Optimization.instantiate_function(optf, x0, + Optimization.AutoSparseReverseDiff(true), nothing, 2) G2 = Array{Float64}(undef, 2) optprob.grad(G2, x0) @@ -544,7 +545,6 @@ H3 = [Array{Float64}(undef, 2, 2), Array{Float64}(undef, 2, 2)] optprob.cons_h(H3, x0) @test H3 ≈ [[2.0 0.0; 0.0 2.0], [-0.0 1.0; 1.0 0.0]] - optf = OptimizationFunction(rosenbrock, Optimization.AutoSparseReverseDiff(), cons = con2_c) optprob = Optimization.instantiate_function(optf, x0, Optimization.AutoSparseReverseDiff(), nothing, 2) @@ -566,7 +566,6 @@ H3 = [Array{Float64}(undef, 2, 2), Array{Float64}(undef, 2, 2)] optprob.cons_h(H3, x0) @test H3 ≈ [[2.0 0.0; 0.0 2.0], [-0.0 1.0; 1.0 0.0]] - optf = OptimizationFunction(rosenbrock, Optimization.AutoSparseReverseDiff()) optprob = Optimization.instantiate_function(optf, x0, Optimization.AutoSparseReverseDiff(), nothing)